Amesos2 - Direct Sparse Solver Interfaces Version of the Day
basker_scalartraits.hpp
1// @HEADER
2// *****************************************************************************
3// Basker: A Direct Linear Solver package
4//
5// Copyright 2011 NTESS and the Basker contributors.
6// SPDX-License-Identifier: LGPL-2.1-or-later
7// *****************************************************************************
8// @HEADER
9
10#ifndef BASKER_SCALARTRAITS_HPP
11#define BASKER_SCALARTRAITS_HPP
12
13#define MY_SCALAR_ABS(x) (((x) < 0) ? -(x) : (x))
14
15
16template <class T>
17struct BASKER_ScalarTraits
18{
19 typedef T magnitudeType;
20 static inline double reciprocal(double c ){return 0;}
21 static inline double divide(double a, double b){return 0;}
22 static inline magnitudeType approxABS(double a){return 0;}
23 static inline magnitudeType abs(double a){return 0;}
24 static inline bool gt (double a, double b){return 0;}
25};
26
27//double
28template <>
29struct BASKER_ScalarTraits<double>
30{
31 typedef double magnitudeType;
32 static inline double reciprocal(double c){return 1.0/c;}
33 static inline double divide(double a, double b){return a/b;}
34 static inline magnitudeType approxABS(double a)
35 { return (MY_SCALAR_ABS(a));}
36 static inline magnitudeType abs(double a)
37 { return (MY_SCALAR_ABS(a));}
38 static inline bool gt (double a, double b){return (a>b);}
39
40};
41
42//float
43template <>
44struct BASKER_ScalarTraits<float>
45{
46 typedef float magnitudeType;
47 static inline float reciprocal(float c){return 1.0/c;}
48 static inline float divide(float a, float b){return a/b;}
49 static inline magnitudeType approxABS(float a)
50 { return (MY_SCALAR_ABS(a));}
51 static inline magnitudeType abs(float a)
52 { return (MY_SCALAR_ABS(a));}
53 static inline bool gt(float a, float b){return (a>b);}
54
55
56};
57
58
59//complex
60#ifdef HAVE_TEUCHOS_COMPLEX
61
62template <class T>
63struct BASKER_ScalarTraits< std::complex<T> >
64{
65 typedef std::complex<T> ComplexT ;
66 typedef typename BASKER_ScalarTraits<T>::magnitudeType magnitudeType ;
67
68 static inline ComplexT reciprocal (ComplexT c)
69 {
70 T r, den, cr, ci ;
71 ComplexT ret ;
72 cr = (Teuchos::ScalarTraits<ComplexT>::real(c)) ;
73 ci = (Teuchos::ScalarTraits<ComplexT>::imag(c)) ;
74 if (MY_SCALAR_ABS (cr) >= MY_SCALAR_ABS (ci))
75 {
76 r = ci / cr ;
77 den = cr + r * ci ;
78 ret = std::complex<T>(1.0 / den, -r / den) ;
79 }
80 else
81 {
82 r = cr / ci ;
83 den = r * cr + ci ;
84 ret = std::complex<T>(r / den, -1.0 / den) ;
85 }
86 return ret;
87 }
88
89 static inline ComplexT divide (ComplexT a, ComplexT b)
90 {
91 T r, den, ar, ai, br, bi ;
92 ComplexT ret;
93
94 br = (Teuchos::ScalarTraits<ComplexT>::real(b)) ;
95 bi = (Teuchos::ScalarTraits<ComplexT>::imag(b)) ;
96 ar = (Teuchos::ScalarTraits<ComplexT>::real(a)) ;
97 ai = (Teuchos::ScalarTraits<ComplexT>::imag(a)) ;
98 if (MY_SCALAR_ABS (br) >= MY_SCALAR_ABS (bi))
99 {
100 r = bi / br ;
101 den = br + r * bi ;
102 ret = std::complex<T>((ar + ai * r) / den, (ai - ar * r) / den) ;
103 }
104 else
105 {
106 r = br / bi ;
107 den = r * br + bi ;
108 ret = std::complex<T>((ar * r + ai) / den, (ai * r - ar) / den) ;
109 }
110 return ret;
111 }
112
113 static inline magnitudeType approxABS (ComplexT a)
114 {
115 return ( MY_SCALAR_ABS (Teuchos::ScalarTraits<ComplexT>::real(a)) +
116 MY_SCALAR_ABS (Teuchos::ScalarTraits<ComplexT>::imag(a)) ) ;
117 }
118
119 static inline magnitudeType abs (ComplexT a)
120 {
121 T r, ar, ai ;
122 magnitudeType s;
123
124 ar = MY_SCALAR_ABS (Teuchos::ScalarTraits<ComplexT>::real(a)) ;
125 ai = MY_SCALAR_ABS (Teuchos::ScalarTraits<ComplexT>::imag(a)) ;
126 if (ar >= ai)
127 {
128 if (ar + ai == ar)
129 {
130 (s) = ar ;
131 }
132 else
133 {
134 r = ai / ar ;
135 (s) = ar * sqrt (1.0 + r*r) ;
136 }
137 }
138 else
139 {
140 if (ai + ar == ai)
141 {
142 (s) = ai ;
143 }
144 else
145 {
146 r = ar / ai ;
147 (s) = ai * sqrt (1.0 + r*r) ;
148 }
149 }
150 return s;
151 }
152 static inline bool gt(ComplexT a, ComplexT b)
153 {
154 return( (Teuchos::ScalarTraits<ComplexT>::real(a)+Teuchos::ScalarTraits<ComplexT>::imag(a)) > (Teuchos::ScalarTraits<ComplexT>::real(b) + Teuchos::ScalarTraits<ComplexT>::imag(b)));
155 }
156
157
158};
159
160#else //C++ complexx
161#include <complex>
162
163template <class T>
164struct BASKER_ScalarTraits< std::complex<T> >
165{
166 typedef std::complex<T> ComplexT ;
167 typedef typename BASKER_ScalarTraits<T>::magnitudeType magnitudeType ;
168
169 static inline ComplexT reciprocal (ComplexT c)
170 {
171 T r, den, cr, ci ;
172 ComplexT ret ;
173 cr = (std::real(c)) ;
174 ci = (std::imag(c)) ;
175 if (MY_SCALAR_ABS (cr) >= MY_SCALAR_ABS (ci))
176 {
177 r = ci / cr ;
178 den = cr + r * ci ;
179 ret = std::complex<T>(1.0 / den, -r / den) ;
180 }
181 else
182 {
183 r = cr / ci ;
184 den = r * cr + ci ;
185 ret = std::complex<T>(r / den, -1.0 / den) ;
186 }
187 return ret;
188 }
189
190 static inline ComplexT divide (ComplexT a, ComplexT b)
191 {
192 T r, den, ar, ai, br, bi ;
193 ComplexT ret;
194
195 br = (std::real(b)) ;
196 bi = (std::imag(b)) ;
197 ar = (std::real(a)) ;
198 ai = (std::imag(a)) ;
199 if (MY_SCALAR_ABS (br) >= MY_SCALAR_ABS (bi))
200 {
201 r = bi / br ;
202 den = br + r * bi ;
203 ret = std::complex<T>((ar + ai * r) / den, (ai - ar * r) / den) ;
204 }
205 else
206 {
207 r = br / bi ;
208 den = r * br + bi ;
209 ret = std::complex<T>((ar * r + ai) / den, (ai * r - ar) / den) ;
210 }
211 return ret;
212 }
213
214 static inline magnitudeType approxABS (ComplexT a)
215 {
216 return ( MY_SCALAR_ABS (std::real(a)) +
217 MY_SCALAR_ABS (std::imag(a)) ) ;
218 }
219
220 static inline magnitudeType abs (ComplexT a)
221 {
222 T r, ar, ai ;
223 magnitudeType s;
224
225 ar = MY_SCALAR_ABS (std::real(a)) ;
226 ai = MY_SCALAR_ABS (std::imag(a)) ;
227 if (ar >= ai)
228 {
229 if (ar + ai == ar)
230 {
231 (s) = ar ;
232 }
233 else
234 {
235 r = ai / ar ;
236 (s) = ar * sqrt (1.0 + r*r) ;
237 }
238 }
239 else
240 {
241 if (ai + ar == ai)
242 {
243 (s) = ai ;
244 }
245 else
246 {
247 r = ar / ai ;
248 (s) = ai * sqrt (1.0 + r*r) ;
249 }
250 }
251 return s;
252 }
253 static inline bool gt(ComplexT a, ComplexT b)
254 {
255 return ((std::real(a)+std::imag(a)) > (std::real(b) + std::imag(b)));
256 }
257
258};
259
260
261#endif // HAVE _TEUCHOS_COMPLEX
262
263#endif
Generated by doxygen 1.9.8