Amesos2 - Direct Sparse Solver Interfaces  Version of the Day
basker_scalartraits.hpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Basker: A Direct Linear Solver package
5 // Copyright 2011 Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, with Sandia Corporation, the
8 // U.S. Government retains certain rights in this software.
9 //
10 // This library is free software; you can redistribute it and/or modify
11 // it under the terms of the GNU Lesser General Public License as
12 // published by the Free Software Foundation; either version 2.1 of the
13 // License, or (at your option) any later version.
14 //
15 // This library is distributed in the hope that it will be useful, but
16 // WITHOUT ANY WARRANTY; without even the implied warranty of
17 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 // Lesser General Public License for more details.
19 //
20 // You should have received a copy of the GNU Lesser General Public
21 // License along with this library; if not, write to the Free Software
22 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
23 // USA
24 // Questions? Contact Mike A. Heroux (maherou@sandia.gov)
25 //
26 // ***********************************************************************
27 // @HEADER
28 
29 #ifndef BASKER_SCALARTRAITS_HPP
30 #define BASKER_SCALARTRAITS_HPP
31 
32 #define MY_SCALAR_ABS(x) (((x) < 0) ? -(x) : (x))
33 
34 
35 template <class T>
36 struct BASKER_ScalarTraits
37 {
38  typedef T magnitudeType;
39  static inline double reciprocal(double c ){return 0;}
40  static inline double divide(double a, double b){return 0;}
41  static inline magnitudeType approxABS(double a){return 0;}
42  static inline magnitudeType abs(double a){return 0;}
43  static inline bool gt (double a, double b){return 0;}
44 };
45 
46 //double
47 template <>
48 struct BASKER_ScalarTraits<double>
49 {
50  typedef double magnitudeType;
51  static inline double reciprocal(double c){return 1.0/c;}
52  static inline double divide(double a, double b){return a/b;}
53  static inline magnitudeType approxABS(double a)
54  { return (MY_SCALAR_ABS(a));}
55  static inline magnitudeType abs(double a)
56  { return (MY_SCALAR_ABS(a));}
57  static inline bool gt (double a, double b){return (a>b);}
58 
59 };
60 
61 //float
62 template <>
63 struct BASKER_ScalarTraits<float>
64 {
65  typedef float magnitudeType;
66  static inline float reciprocal(float c){return 1.0/c;}
67  static inline float divide(float a, float b){return a/b;}
68  static inline magnitudeType approxABS(float a)
69  { return (MY_SCALAR_ABS(a));}
70  static inline magnitudeType abs(float a)
71  { return (MY_SCALAR_ABS(a));}
72  static inline bool gt(float a, float b){return (a>b);}
73 
74 
75 };
76 
77 
78 //complex
79 #ifdef HAVE_TEUCHOS_COMPLEX
80 
81 template <class T>
82 struct BASKER_ScalarTraits< std::complex<T> >
83 {
84  typedef std::complex<T> ComplexT ;
85  typedef typename BASKER_ScalarTraits<T>::magnitudeType magnitudeType ;
86 
87  static inline ComplexT reciprocal (ComplexT c)
88  {
89  T r, den, cr, ci ;
90  ComplexT ret ;
91  cr = (Teuchos::ScalarTraits<ComplexT>::real(c)) ;
92  ci = (Teuchos::ScalarTraits<ComplexT>::imag(c)) ;
93  if (MY_SCALAR_ABS (cr) >= MY_SCALAR_ABS (ci))
94  {
95  r = ci / cr ;
96  den = cr + r * ci ;
97  ret = std::complex<T>(1.0 / den, -r / den) ;
98  }
99  else
100  {
101  r = cr / ci ;
102  den = r * cr + ci ;
103  ret = std::complex<T>(r / den, -1.0 / den) ;
104  }
105  return ret;
106  }
107 
108  static inline ComplexT divide (ComplexT a, ComplexT b)
109  {
110  T r, den, ar, ai, br, bi ;
111  ComplexT ret;
112 
113  br = (Teuchos::ScalarTraits<ComplexT>::real(b)) ;
114  bi = (Teuchos::ScalarTraits<ComplexT>::imag(b)) ;
115  ar = (Teuchos::ScalarTraits<ComplexT>::real(a)) ;
116  ai = (Teuchos::ScalarTraits<ComplexT>::imag(a)) ;
117  if (MY_SCALAR_ABS (br) >= MY_SCALAR_ABS (bi))
118  {
119  r = bi / br ;
120  den = br + r * bi ;
121  ret = std::complex<T>((ar + ai * r) / den, (ai - ar * r) / den) ;
122  }
123  else
124  {
125  r = br / bi ;
126  den = r * br + bi ;
127  ret = std::complex<T>((ar * r + ai) / den, (ai * r - ar) / den) ;
128  }
129  return ret;
130  }
131 
132  static inline magnitudeType approxABS (ComplexT a)
133  {
134  return ( MY_SCALAR_ABS (Teuchos::ScalarTraits<ComplexT>::real(a)) +
135  MY_SCALAR_ABS (Teuchos::ScalarTraits<ComplexT>::imag(a)) ) ;
136  }
137 
138  static inline magnitudeType abs (ComplexT a)
139  {
140  T r, ar, ai ;
141  magnitudeType s;
142 
143  ar = MY_SCALAR_ABS (Teuchos::ScalarTraits<ComplexT>::real(a)) ;
144  ai = MY_SCALAR_ABS (Teuchos::ScalarTraits<ComplexT>::imag(a)) ;
145  if (ar >= ai)
146  {
147  if (ar + ai == ar)
148  {
149  (s) = ar ;
150  }
151  else
152  {
153  r = ai / ar ;
154  (s) = ar * sqrt (1.0 + r*r) ;
155  }
156  }
157  else
158  {
159  if (ai + ar == ai)
160  {
161  (s) = ai ;
162  }
163  else
164  {
165  r = ar / ai ;
166  (s) = ai * sqrt (1.0 + r*r) ;
167  }
168  }
169  return s;
170  }
171  static inline bool gt(ComplexT a, ComplexT b)
172  {
173  return( (Teuchos::ScalarTraits<ComplexT>::real(a)+Teuchos::ScalarTraits<ComplexT>::imag(a)) > (Teuchos::ScalarTraits<ComplexT>::real(b) + Teuchos::ScalarTraits<ComplexT>::imag(b)));
174  }
175 
176 
177 };
178 
179 #else //C++ complexx
180 #include <complex>
181 
182 template <class T>
183 struct BASKER_ScalarTraits< std::complex<T> >
184 {
185  typedef std::complex<T> ComplexT ;
186  typedef typename BASKER_ScalarTraits<T>::magnitudeType magnitudeType ;
187 
188  static inline ComplexT reciprocal (ComplexT c)
189  {
190  T r, den, cr, ci ;
191  ComplexT ret ;
192  cr = (std::real(c)) ;
193  ci = (std::imag(c)) ;
194  if (MY_SCALAR_ABS (cr) >= MY_SCALAR_ABS (ci))
195  {
196  r = ci / cr ;
197  den = cr + r * ci ;
198  ret = std::complex<T>(1.0 / den, -r / den) ;
199  }
200  else
201  {
202  r = cr / ci ;
203  den = r * cr + ci ;
204  ret = std::complex<T>(r / den, -1.0 / den) ;
205  }
206  return ret;
207  }
208 
209  static inline ComplexT divide (ComplexT a, ComplexT b)
210  {
211  T r, den, ar, ai, br, bi ;
212  ComplexT ret;
213 
214  br = (std::real(b)) ;
215  bi = (std::imag(b)) ;
216  ar = (std::real(a)) ;
217  ai = (std::imag(a)) ;
218  if (MY_SCALAR_ABS (br) >= MY_SCALAR_ABS (bi))
219  {
220  r = bi / br ;
221  den = br + r * bi ;
222  ret = std::complex<T>((ar + ai * r) / den, (ai - ar * r) / den) ;
223  }
224  else
225  {
226  r = br / bi ;
227  den = r * br + bi ;
228  ret = std::complex<T>((ar * r + ai) / den, (ai * r - ar) / den) ;
229  }
230  return ret;
231  }
232 
233  static inline magnitudeType approxABS (ComplexT a)
234  {
235  return ( MY_SCALAR_ABS (std::real(a)) +
236  MY_SCALAR_ABS (std::imag(a)) ) ;
237  }
238 
239  static inline magnitudeType abs (ComplexT a)
240  {
241  T r, ar, ai ;
242  magnitudeType s;
243 
244  ar = MY_SCALAR_ABS (std::real(a)) ;
245  ai = MY_SCALAR_ABS (std::imag(a)) ;
246  if (ar >= ai)
247  {
248  if (ar + ai == ar)
249  {
250  (s) = ar ;
251  }
252  else
253  {
254  r = ai / ar ;
255  (s) = ar * sqrt (1.0 + r*r) ;
256  }
257  }
258  else
259  {
260  if (ai + ar == ai)
261  {
262  (s) = ai ;
263  }
264  else
265  {
266  r = ar / ai ;
267  (s) = ai * sqrt (1.0 + r*r) ;
268  }
269  }
270  return s;
271  }
272  static inline bool gt(ComplexT a, ComplexT b)
273  {
274  return ((std::real(a)+std::imag(a)) > (std::real(b) + std::imag(b)));
275  }
276 
277 };
278 
279 
280 #endif // HAVE _TEUCHOS_COMPLEX
281 
282 #endif
Generated on Thu Apr 18 2024 09:18:29 for Amesos2 - Direct Sparse Solver Interfaces by   doxygen 1.8.5