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
35template <class T>
36struct 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
47template <>
48struct 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
62template <>
63struct 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
81template <class T>
82struct 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
182template <class T>
183struct 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