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 
16 template <class T>
17 struct 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
28 template <>
29 struct 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
43 template <>
44 struct 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 
62 template <class T>
63 struct 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 
163 template <class T>
164 struct 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 on Fri Nov 22 2024 09:16:49 for Amesos2 - Direct Sparse Solver Interfaces by   doxygen 1.8.5