DMRITool  v0.1.1-139-g860d86b4
Diffusion MRI Tool
utlLapack.h
Go to the documentation of this file.
1 
18 #ifndef __utlLapack_h
19 #define __utlLapack_h
20 
24 #include "utlCoreMacro.h"
25 
26 #ifdef __cplusplus
27 extern "C"
28 {
29 #endif //__cplusplus
30 
31 #ifdef UTL_USE_MKL
32  #include "utlMKL.h"
33  #include <mkl_lapacke.h>
34  #include <mkl_lapack.h>
35  #include <mkl.h>
36  #include <complex>
37 #else
38 
39  #define lapack_complex_double std::complex<double>
40  #define lapack_complex_float std::complex<float>
41  #include <lapacke.h>
42 #endif
43 
44 
45 // #if defined(UTL_USE_MKL)
46 // #include <mkl_lapack.h>
47 // #else
48 // #if !defined(_MKL_LAPACK_H_) && !defined(_LAPACKE_H_)
49 
50 // [>* http:www.netlib.org/lapack/explore-html/dd/d4c/dsyev_8f.html <]
51 // extern void dsyev_(char* JOBZ, char* UPLO, int* N, double *A, int* LDA, double* W, double* WORK, int* LWORK, int* INFO);
52 // extern void ssyev_(char* JOBZ, char* UPLO, int* N, float *A, int* LDA, float* W, float* WORK, int* LWORK, int* INFO);
53 // [>* http:www.netlib.org/lapack/explore-html/d1/da2/dsyevd_8f.html <]
54 // extern void dsyevd_(char* JOBZ, char* UPLO, int* N, double *A, int* LDA, double* W, double* WORK, int* LWORK, int* IWORK, int* LIWORK, int* INFO);
55 // extern void ssyevd_(char* JOBZ, char* UPLO, int* N, float *A, int* LDA, float* W, float* WORK, int* LWORK, int* IWORK, int* LIWORK, int* INFO);
56 
57 // [>* http:www.netlib.org/lapack/explore-html/d8/d2d/dgesvd_8f.html <]
58 // extern void dgesvd_(char* JOBU, char* JOBVT, int* M, int* N, double *A, int* LDA, double* S, double* U, int* LDU, double* VT, int* LDVT, double* WORK, int* LWORK, int* INFO);
59 // extern void sgesvd_(char* JOBU, char* JOBVT, int* M, int* N, float *A, int* LDA, float* S, float* U, int* LDU, float* VT, int* LDVT, float* WORK, int* LWORK, int* INFO);
60 // [>* http:www.netlib.org/lapack/explore-html/db/db4/dgesdd_8f.html <]
61 // extern void dgesdd_(char* JOBZ, int* M, int* N, double *A, int* LDA, double* S, double* U, int* LDU, double* VT, int* LDVT, double* WORK, int* LWORK, int* IWORK, int* INFO);
62 // extern void sgesdd_(char* JOBZ, int* M, int* N, float *A, int* LDA, float* S, float* U, int* LDU, float* VT, int* LDVT, float* WORK, int* LWORK, int* IWORK, int* INFO);
63 
64 // [>* http:www.netlib.org/lapack/explore-html/d5/d8f/dsytri_8f.html <]
65 // extern void dsytri_(char* UPLO, int* N, double *A, int* LDA, int* IPIV, double* WORK, int* INFO);
66 // extern void ssytri_(char* UPLO, int* N, float *A, int* LDA, int* IPIV, float* WORK, int* INFO);
67 
68 // [>* http:www.netlib.org/lapack/explore-html/dd/df4/dsytrf_8f.html <]
69 // extern void dsytrf_(char* UPLO, int* N, double *A, int* LDA, int* IPIV, double* WORK, int* LWORK, int* INFO);
70 // extern void ssytrf_(char* UPLO, int* N, float *A, int* LDA, int* IPIV, float* WORK, int* LWORK, int* INFO);
71 // #endif
72 
73 // #endif
74 
75 float LAPACKE_slange( int matrix_order, char norm, lapack_int m, lapack_int n, const float* a, lapack_int lda );
76 double LAPACKE_dlange( int matrix_order, char norm, lapack_int m, lapack_int n, const double* a, lapack_int lda );
77 double LAPACKE_zlange( int matrix_order, char norm, lapack_int m, lapack_int n, const std::complex<double>* a, lapack_int lda );
78 float LAPACKE_clange( int matrix_order, char norm, lapack_int m, lapack_int n, const std::complex<float>* a, lapack_int lda );
79 
80 
81 #ifdef __cplusplus
82  }
83 #endif /* __cplusplus */
84 
85 namespace utl
86 {
87 
88 // [>************************************************************************************<]
89 // [>* template function definitions <]
90 
91 // template <class T> inline void
92 // syev(char* JOBZ, char* UPLO, int* N, T *A, int* LDA, T* W, T* WORK, int* LWORK, int* INFO);
93 // template <class T> inline void
94 // syevd(char* JOBZ, char* UPLO, int* N, T *A, int* LDA, T* W, T* WORK, int* LWORK, int* IWORK, int* LIWORK, int* INFO);
95 
96 // template <class T> inline void
97 // gesvd(char* JOBU, char* JOBVT, int* M, int* N, T *A, int* LDA, T* S, T* U, int* LDU, T* VT, int* LDVT, T* WORK, int* LWORK, int* INFO);
98 // template <class T> inline void
99 // gesdd(char* JOBZ, int* M, int* N, T *A, int* LDA, T* S, T* U, int* LDU, T* VT, int* LDVT, T* WORK, int* LWORK, int* IWORK, int* INFO);
100 
101 // template <class T> inline void
102 // sytri(char* UPLO, int* N, T *A, int* LDA, int* IPIV, T* WORK, int* INFO);
103 
104 // template <class T> inline void
105 // sytrf(char* UPLO, int* N, T* A, int* LDA, int* IPIV, T* WORK, int* LWORK, int* INFO);
106 
107 // template <class T> inline T
108 // lange(char* NORM, int* M, int* N, const T* A, int* LDA, T* WORK);
109 
110 
111 
112 // // ************************************************************************************
113 // [>* function implementations using double and float <]
114 
115 // template <> inline void
116 // syev<double>(char* JOBZ, char* UPLO, int* N, double *A, int* LDA, double* W, double* WORK, int* LWORK, int* INFO)
117 // { dsyev_(JOBZ,UPLO,N,A,LDA,W,WORK,LWORK,INFO);}
118 // template <> inline void
119 // syev<float>(char* JOBZ, char* UPLO, int* N, float *A, int* LDA, float* W, float* WORK, int* LWORK, int* INFO)
120 // { ssyev_(JOBZ,UPLO,N,A,LDA,W,WORK,LWORK,INFO);}
121 
122 // template <> inline void
123 // syevd<double>(char* JOBZ, char* UPLO, int* N, double *A, int* LDA, double* W, double* WORK, int* LWORK, int* IWORK, int* LIWORK, int* INFO)
124 // { dsyevd_(JOBZ,UPLO,N,A,LDA,W,WORK,LWORK,IWORK,LIWORK,INFO);}
125 // template <> inline void
126 // syevd<float>(char* JOBZ, char* UPLO, int* N, float *A, int* LDA, float* W, float* WORK, int* LWORK, int* IWORK, int* LIWORK, int* INFO)
127 // { ssyevd_(JOBZ,UPLO,N,A,LDA,W,WORK,LWORK,IWORK,LIWORK,INFO);}
128 
129 // template <> inline void
130 // gesvd<double>(char* JOBU, char* JOBVT, int* M, int* N, double *A, int* LDA, double* S, double* U, int* LDU, double* VT, int* LDVT, double* WORK, int* LWORK, int* INFO)
131 // { dgesvd_(JOBU,JOBVT,M,N,A,LDA,S,U,LDU,VT,LDVT,WORK,LWORK,INFO);}
132 // template <> inline void
133 // gesvd<float>(char* JOBU, char* JOBVT, int* M, int* N, float *A, int* LDA, float* S, float* U, int* LDU, float* VT, int* LDVT, float* WORK, int* LWORK, int* INFO)
134 // { sgesvd_(JOBU,JOBVT,M,N,A,LDA,S,U,LDU,VT,LDVT,WORK,LWORK,INFO);}
135 
136 // template <> inline void
137 // gesdd<double>(char* JOBZ, int* M, int* N, double *A, int* LDA, double* S, double* U, int* LDU, double* VT, int* LDVT, double* WORK, int* LWORK, int* IWORK, int* INFO)
138 // { dgesdd_(JOBZ,M,N,A,LDA,S,U,LDU,VT,LDVT,WORK,LWORK,IWORK,INFO);}
139 // template <> inline void
140 // gesdd<float>(char* JOBZ, int* M, int* N, float *A, int* LDA, float* S, float* U, int* LDU, float* VT, int* LDVT, float* WORK, int* LWORK, int* IWORK, int* INFO)
141 // { sgesdd_(JOBZ,M,N,A,LDA,S,U,LDU,VT,LDVT,WORK,LWORK,IWORK,INFO);}
142 
143 // template <> inline void
144 // sytri<double>(char* UPLO, int* N, double *A, int* LDA, int* IPIV, double* WORK, int* INFO)
145 // { dsytri_(UPLO,N,A,LDA,IPIV,WORK,INFO);}
146 // template <> inline void
147 // sytri<float>(char* UPLO, int* N, float *A, int* LDA, int* IPIV, float* WORK, int* INFO)
148 // { ssytri_(UPLO,N,A,LDA,IPIV,WORK,INFO);}
149 
150 // template <class T> inline void
151 // sytrf(char* UPLO, int* N, double* A, int* LDA, int* IPIV, double* WORK, int* LWORK, int* INFO)
152 // { dsytrf_(UPLO,N,A,LDA,IPIV,WORK,LWORK,INFO);}
153 
154 // template <class T> inline void
155 // sytrf(char* UPLO, int* N, float* A, int* LDA, int* IPIV, float* WORK, int* LWORK, int* INFO)
156 // { ssytrf_(UPLO,N,A,LDA,IPIV,WORK,LWORK,INFO);}
157 
158 // template <> inline double
159 // lange<double>( char* norm, int* m, int* n, const double* a, int* lda, double* work )
160 // { return dlange_(norm,m,n,a,lda,work); }
161 // template <> inline float
162 // lange<float>( char* norm, int* m, int* n, const float* a, int* lda, float* work )
163 // { return slange_(norm,m,n,a,lda,work); }
164 
165 
166 /**************************************************************************************/
169 template <class T> inline int
170 syev(int matrix_order, char JOBZ, char UPLO, int N, T *A, int LDA, T* W);
171 template <class T> inline int
172 syevd(int matrix_order, char JOBZ, char UPLO, int N, T *A, int LDA, T* W);
173 template <class T> inline int
174 geev( int matrix_layout, char jobvl, char jobvr, int n, T* a, int lda, T* wr, T* wi, T* vl, int ldvl, T* vr, int ldvr );
175 
176 template <class T> inline int
177 gesvd(int matrix_order, char JOBU, char JOBVT, int M, int N, T *A, int LDA, T* S, T* U, int LDU, T* VT, int LDVT, T* superb);
178 template <class T> inline int
179 gesdd(int matrix_order, char JOBZ, int M, int N, T *A, int LDA, T* S, T* U, int LDU, T* VT, int LDVT);
180 
181 template <class T> inline int
182 sytri(int matrix_order, char UPLO, int N, T *A, int LDA, const int* IPIV);
183 template <class T> inline int
184 sytrf(int matrix_order, char UPLO, int N, T* A, int LDA, int* IPIV);
185 
186 template <class T> inline T
187 lange(int matrix_order, char norm, int m, int n, const T* A, int LDA);
188 
189 
190 template <class T> inline int
191 getrf (int matrix_layout , int m , int n , T * a , int lda, int * ipiv );
192 
193 template <class T> inline int
194 getri(int matrix_layout , int n , T * a , int lda , const int * ipiv);
195 
196 /**************************************************************************************/
199 template <> inline int
200 syev<double>(int matrix_order, char JOBZ, char UPLO, int N, double *A, int LDA, double* W)
201 { return LAPACKE_dsyev(matrix_order,JOBZ,UPLO,N,A,LDA,W);}
202 template <> inline int
203 syev<float>(int matrix_order, char JOBZ, char UPLO, int N, float *A, int LDA, float* W)
204 { return LAPACKE_ssyev(matrix_order,JOBZ,UPLO,N,A,LDA,W);}
205 
206 template <> inline int
207 syevd<double>(int matrix_order, char JOBZ, char UPLO, int N, double *A, int LDA, double* W)
208 { return LAPACKE_dsyevd(matrix_order,JOBZ,UPLO,N,A,LDA,W);}
209 template <> inline int
210 syevd<float>(int matrix_order, char JOBZ, char UPLO, int N, float *A, int LDA, float* W)
211 { return LAPACKE_ssyevd(matrix_order,JOBZ,UPLO,N,A,LDA,W);}
212 template <> inline int
213 syevd<std::complex<double> >(int matrix_order, char JOBZ, char UPLO, int N, std::complex<double> *A, int LDA, std::complex<double>* W)
214 { utlGlobalException(true, "connot use complex values"); return 0;}
215 template <> inline int
216 syevd<std::complex<float> >(int matrix_order, char JOBZ, char UPLO, int N, std::complex<float> *A, int LDA, std::complex<float>* W)
217 { utlGlobalException(true, "connot use complex values"); return 0;}
218 
219 template <> inline int
220 geev<double>( int matrix_layout, char jobvl, char jobvr, int n, double* a, int lda, double* wr, double* wi, double* vl, int ldvl, double* vr, int ldvr )
221 { return LAPACKE_dgeev( matrix_layout, jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr ); }
222 template <> inline int
223 geev<float>( int matrix_layout, char jobvl, char jobvr, int n, float* a, int lda, float* wr, float* wi, float* vl, int ldvl, float* vr, int ldvr )
224 { return LAPACKE_sgeev( matrix_layout, jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr ); }
225 template <> inline int
226 geev<std::complex<double> >( int matrix_layout, char jobvl, char jobvr, int n, std::complex<double>* a, int lda, std::complex<double>* wr, std::complex<double>* wi, std::complex<double>* vl, int ldvl, std::complex<double>* vr, int ldvr )
227 { utlGlobalException(true, "connot use complex values"); return 0; }
228 template <> inline int
229 geev<std::complex<float> >( int matrix_layout, char jobvl, char jobvr, int n, std::complex<float>* a, int lda, std::complex<float>* wr, std::complex<float>* wi, std::complex<float>* vl, int ldvl, std::complex<float>* vr, int ldvr )
230 { utlGlobalException(true, "connot use complex values"); return 0; }
231 inline int
232 geev( int matrix_layout, char jobvl, char jobvr, int n, std::complex<double>* a, int lda, std::complex<double>* w, std::complex<double>* vl, int ldvl, std::complex<double>* vr, int ldvr )
233 { return LAPACKE_zgeev( matrix_layout, jobvl, jobvr, n, a, lda, w, vl, ldvl, vr, ldvr ); }
234 
235 template <> inline int
236 gesvd<double>(int matrix_order, char JOBU, char JOBVT, int M, int N, double *A, int LDA, double* S, double* U, int LDU, double* VT, int LDVT, double* superb)
237 { return LAPACKE_dgesvd(matrix_order,JOBU,JOBVT,M,N,A,LDA,S,U,LDU,VT,LDVT,superb);}
238 template <> inline int
239 gesvd<float>(int matrix_order, char JOBU, char JOBVT, int M, int N, float *A, int LDA, float* S, float* U, int LDU, float* VT, int LDVT, float* superb)
240 { return LAPACKE_sgesvd(matrix_order,JOBU,JOBVT,M,N,A,LDA,S,U,LDU,VT,LDVT,superb);}
241 
242 template <> inline int
243 gesdd<double>(int matrix_order, char JOBZ, int M, int N, double *A, int LDA, double* S, double* U, int LDU, double* VT, int LDVT)
244 { return LAPACKE_dgesdd(matrix_order,JOBZ,M,N,A,LDA,S,U,LDU,VT,LDVT);}
245 template <> inline int
246 gesdd<float>(int matrix_order, char JOBZ, int M, int N, float *A, int LDA, float* S, float* U, int LDU, float* VT, int LDVT)
247 { return LAPACKE_sgesdd(matrix_order,JOBZ,M,N,A,LDA,S,U,LDU,VT,LDVT);}
248 inline int
249 gesdd(int matrix_order, char JOBZ, int M, int N, std::complex<double> *A, int LDA, double* S, std::complex<double>* U, int LDU, std::complex<double>* VT, int LDVT)
250 { return LAPACKE_zgesdd(matrix_order,JOBZ,M,N,A,LDA,S,U,LDU,VT,LDVT);}
251 inline int
252 gesdd(int matrix_order, char JOBZ, int M, int N, std::complex<float> *A, int LDA, float* S, std::complex<float>* U, int LDU, std::complex<float>* VT, int LDVT)
253 { return LAPACKE_cgesdd(matrix_order,JOBZ,M,N,A,LDA,S,U,LDU,VT,LDVT);}
254 
255 template <> inline int
256 sytri<double>(int matrix_order, char UPLO, int N, double *A, int LDA, const int* IPIV)
257 { return LAPACKE_dsytri(matrix_order,UPLO,N,A,LDA,IPIV);}
258 template <> inline int
259 sytri<float>(int matrix_order, char UPLO, int N, float *A, int LDA, const int* IPIV)
260 { return LAPACKE_ssytri(matrix_order,UPLO,N,A,LDA,IPIV); }
261 template <> inline int
262 sytri<std::complex<double> >(int matrix_order, char UPLO, int N, std::complex<double> *A, int LDA, const int* IPIV)
263 { return LAPACKE_zsytri(matrix_order,UPLO,N,A,LDA,IPIV);}
264 template <> inline int
265 sytri<std::complex<float> >(int matrix_order, char UPLO, int N, std::complex<float> *A, int LDA, const int* IPIV)
266 { return LAPACKE_csytri(matrix_order,UPLO,N,A,LDA,IPIV);}
267 
268 template <> inline int
269 sytrf<double>(int matrix_order, char UPLO, int N, double* A, int LDA, int* IPIV)
270 { return LAPACKE_dsytrf(matrix_order,UPLO,N,A,LDA,IPIV); }
271 template <> inline int
272 sytrf<float>(int matrix_order, char UPLO, int N, float* A, int LDA, int* IPIV)
273 { return LAPACKE_ssytrf(matrix_order,UPLO,N,A,LDA,IPIV); }
274 template <> inline int
275 sytrf<std::complex<double> >(int matrix_order, char UPLO, int N, std::complex<double>* A, int LDA, int* IPIV)
276 { return LAPACKE_zsytrf(matrix_order,UPLO,N,A,LDA,IPIV); }
277 template <> inline int
278 sytrf<std::complex<float> >(int matrix_order, char UPLO, int N, std::complex<float>* A, int LDA, int* IPIV)
279 { return LAPACKE_csytrf(matrix_order,UPLO,N,A,LDA,IPIV); }
280 
281 template <> inline double
282 lange<double>(int matrix_order, char norm, int m, int n, const double* A, int LDA)
283 { return LAPACKE_dlange(matrix_order,norm,m,n,A,LDA); }
284 template <> inline float
285 lange<float>(int matrix_order, char norm, int m, int n, const float* A, int LDA)
286 { return LAPACKE_slange(matrix_order,norm,m,n,A,LDA); }
287 inline double
288 lange(int matrix_order, char norm, int m, int n, const std::complex<double>* A, int LDA)
289 { return LAPACKE_zlange(matrix_order,norm,m,n,A,LDA); }
290 inline float
291 lange(int matrix_order, char norm, int m, int n, const std::complex<float>* A, int LDA)
292 { return LAPACKE_clange(matrix_order,norm,m,n,A,LDA); }
293 
294 template <> inline int
295 getrf<double>(int matrix_layout , int m , int n , double * a , int lda, int * ipiv )
296 { return LAPACKE_dgetrf(matrix_layout, m, n, a, lda, ipiv); }
297 template <> inline int
298 getrf<float>(int matrix_layout , int m , int n , float * a , int lda, int * ipiv )
299 { return LAPACKE_sgetrf(matrix_layout, m, n, a, lda, ipiv); }
300 template <> inline int
301 getrf<std::complex<double> >(int matrix_layout , int m , int n , std::complex<double> * a , int lda, int * ipiv )
302 { return LAPACKE_zgetrf(matrix_layout, m, n, a, lda, ipiv); }
303 template <> inline int
304 getrf<std::complex<float> >(int matrix_layout , int m , int n , std::complex<float> * a , int lda, int * ipiv )
305 { return LAPACKE_cgetrf(matrix_layout, m, n, a, lda, ipiv); }
306 
307 template <> inline int
308 getri<double>(int matrix_layout , int n , double * a , int lda , const int * ipiv)
309 { return LAPACKE_dgetri(matrix_layout, n, a, lda, ipiv); }
310 template <> inline int
311 getri<float>(int matrix_layout , int n , float * a , int lda , const int * ipiv)
312 { return LAPACKE_sgetri(matrix_layout, n, a, lda, ipiv); }
313 template <> inline int
314 getri<std::complex<double> >(int matrix_layout , int n , std::complex<double> * a , int lda , const int * ipiv)
315 { return LAPACKE_zgetri(matrix_layout, n, a, lda, ipiv); }
316 template <> inline int
317 getri<std::complex<float> >(int matrix_layout , int n , std::complex<float> * a , int lda , const int * ipiv)
318 { return LAPACKE_cgetri(matrix_layout, n, a, lda, ipiv); }
319 
320 
321 #define __utl_getri_Matrix(T, FuncName, RowMajorMatrixName, GetRowsFuncName, GetColsFuncName, MatrixGetDataFuncName, ReSizeFuncName) \
322 template <class T> inline void \
323 FuncName(const RowMajorMatrixName& mat, RowMajorMatrixName& result) \
324 { \
325  int N = mat.GetRowsFuncName(),INFO=0; \
326  int ipiv[N]; \
327  result = mat; \
328  INFO = getrf<T>(LAPACK_ROW_MAJOR, N,N,result.MatrixGetDataFuncName(),N,ipiv); \
329  INFO = getri<T>(LAPACK_ROW_MAJOR, N,result.MatrixGetDataFuncName(),N,ipiv); \
330 }
331 
332 
355 #define __utl_geev_Matrix(T, FuncName, RowMajorMatrixName, GetRowsFuncName, GetColsFuncName, MatrixGetDataFuncName, VectorName, VectorGetDataFuncName, ReSizeFuncName) \
356 template <class T> inline void \
357 FuncName ( const RowMajorMatrixName& mat, VectorName& valReal, VectorName& valImg, RowMajorMatrixName& vecRealR, RowMajorMatrixName& vecImgR, RowMajorMatrixName& vecRealL, RowMajorMatrixName& vecImgL) \
358 { \
359  utlException(mat.GetRowsFuncName()!=mat.GetColsFuncName(), "The matrix should be square"); \
360  char JOBVL='V', JOBVR='V'; \
361  int N = mat.GetRowsFuncName(),INFO=0; \
362  int LDA=N, LDVL=N, LDVR=N; \
363  \
364  RowMajorMatrixName matCopy(mat), vecR(N,N), vecL(N,N); \
365  valReal.ReSizeFuncName(N); \
366  valImg.ReSizeFuncName(N); \
367  vecRealR.ReSizeFuncName(N,N); \
368  vecImgR.ReSizeFuncName(N,N); \
369  vecRealL.ReSizeFuncName(N,N); \
370  vecImgL.ReSizeFuncName(N,N); \
371  \
372  INFO = geev<T>( LAPACK_ROW_MAJOR, JOBVL, JOBVR, N, matCopy.MatrixGetDataFuncName(), LDA, valReal.VectorGetDataFuncName(), valImg.VectorGetDataFuncName(), vecL.MatrixGetDataFuncName(), LDVL, vecR.MatrixGetDataFuncName(), LDVR ); \
373  \
374  T const * vecR_data=vecR.MatrixGetDataFuncName(); \
375  T* vecRealR_data=vecRealR.MatrixGetDataFuncName(); \
376  T* vecImgR_data=vecImgR.MatrixGetDataFuncName(); \
377  T const * vecL_data=vecL.MatrixGetDataFuncName(); \
378  T* vecRealL_data=vecRealL.MatrixGetDataFuncName(); \
379  T* vecImgL_data=vecImgL.MatrixGetDataFuncName(); \
380  for ( int i = 0; i < N; ++i ) \
381  { \
382  int j=0; \
383  int k0 = i*LDVR; \
384  while( j < N ) \
385  { \
386  int kk=k0+j; \
387  if( valImg[j] == (double)0.0 ) \
388  { \
389  vecRealR_data[kk] = vecR_data[kk]; \
390  vecImgR_data[kk] = 0.0; \
391  vecRealL_data[kk] = vecL_data[kk]; \
392  vecImgL_data[kk] = 0.0; \
393  j++; \
394  } \
395  else \
396  { \
397  vecRealR_data[kk] = vecR_data[kk]; \
398  vecImgR_data[kk] = vecR_data[kk+1]; \
399  \
400  vecRealR_data[kk+1] = vecR_data[kk]; \
401  vecImgR_data[kk+1] = -vecR_data[kk+1]; \
402  \
403  vecRealL_data[kk] = vecL_data[kk]; \
404  vecImgL_data[kk] = vecL_data[kk+1]; \
405  \
406  vecRealL_data[kk+1] = vecL_data[kk]; \
407  vecImgL_data[kk+1] = -vecL_data[kk+1]; \
408  \
409  j += 2; \
410  } \
411  } \
412  } \
413 } \
414  \
415 template <class T> inline void \
416 FuncName ( const RowMajorMatrixName& mat, VectorName& valReal, VectorName& valImg) \
417 { \
418  utlException(mat.GetRowsFuncName()!=mat.GetColsFuncName(), "The matrix should be square"); \
419  char JOBVL='N', JOBVR='N'; \
420  int N = mat.GetRowsFuncName(),INFO=0; \
421  int LDA=N, LDVL=N, LDVR=N; \
422  \
423  RowMajorMatrixName matCopy(mat); \
424  valReal.ReSizeFuncName(N); \
425  valImg.ReSizeFuncName(N); \
426  \
427  INFO = geev<T>( LAPACK_ROW_MAJOR, JOBVL, JOBVR, N, matCopy.MatrixGetDataFuncName(), LDA, valReal.VectorGetDataFuncName(), valImg.VectorGetDataFuncName(), NULL, LDVL, NULL, LDVR ); \
428 } \
429  \
430 template <class T> inline void \
431 FuncName ( const RowMajorMatrixName& mat, VectorName& valReal, VectorName& valImg, RowMajorMatrixName& vecRealR, RowMajorMatrixName& vecImgR) \
432 { \
433  utlException(mat.GetRowsFuncName()!=mat.GetColsFuncName(), "The matrix should be square"); \
434  char JOBVL='N', JOBVR='V'; \
435  int N = mat.GetRowsFuncName(),INFO=0; \
436  int LDA=N, LDVL=N, LDVR=N; \
437  \
438  RowMajorMatrixName matCopy(mat), eigenVectorCopy(N,N); \
439  valReal.ReSizeFuncName(N); \
440  valImg.ReSizeFuncName(N); \
441  vecRealR.ReSizeFuncName(N,N); \
442  vecImgR.ReSizeFuncName(N,N); \
443  \
444  INFO = geev<T>( LAPACK_ROW_MAJOR, JOBVL, JOBVR, N, matCopy.MatrixGetDataFuncName(), LDA, valReal.VectorGetDataFuncName(), valImg.VectorGetDataFuncName(), NULL, LDVL, eigenVectorCopy.MatrixGetDataFuncName(), LDVR ); \
445  \
446  T const * vecR_data=eigenVectorCopy.MatrixGetDataFuncName(); \
447  T* vecRealR_data=vecRealR.MatrixGetDataFuncName(); \
448  T* vecImgR_data=vecImgR.MatrixGetDataFuncName(); \
449  for ( int i = 0; i < N; ++i ) \
450  { \
451  int j=0; \
452  int k0 = i*LDVR; \
453  while( j < N ) \
454  { \
455  int kk=k0+j; \
456  if( valImg[j] == (double)0.0 ) \
457  { \
458  vecRealR_data[kk] = vecR_data[kk]; \
459  vecImgR_data[kk] = 0.0; \
460  j++; \
461  } \
462  else \
463  { \
464  vecRealR_data[kk] = vecR_data[kk]; \
465  vecImgR_data[kk] = vecR_data[kk+1]; \
466  \
467  vecRealR_data[kk+1] = vecR_data[kk]; \
468  vecImgR_data[kk+1] = -vecR_data[kk+1]; \
469  j += 2; \
470  } \
471  } \
472  } \
473 } \
474 
475 
477 }
478 
479 #endif
utl::gesvd< float >
int gesvd< float >(int matrix_order, char JOBU, char JOBVT, int M, int N, float *A, int LDA, float *S, float *U, int LDU, float *VT, int LDVT, float *superb)
Definition: utlLapack.h:239
utl::sytrf
int sytrf(int matrix_order, char UPLO, int N, T *A, int LDA, int *IPIV)
utl::getri< double >
int getri< double >(int matrix_layout, int n, double *a, int lda, const int *ipiv)
Definition: utlLapack.h:308
utl::sytrf< double >
int sytrf< double >(int matrix_order, char UPLO, int N, double *A, int LDA, int *IPIV)
Definition: utlLapack.h:269
utl::syevd< float >
int syevd< float >(int matrix_order, char JOBZ, char UPLO, int N, float *A, int LDA, float *W)
Definition: utlLapack.h:210
LAPACKE_slange
float LAPACKE_slange(int matrix_order, char norm, lapack_int m, lapack_int n, const float *a, lapack_int lda)
utl::sytri
int sytri(int matrix_order, char UPLO, int N, T *A, int LDA, const int *IPIV)
utl::getri
int getri(int matrix_layout, int n, T *a, int lda, const int *ipiv)
utl::sytri< float >
int sytri< float >(int matrix_order, char UPLO, int N, float *A, int LDA, const int *IPIV)
Definition: utlLapack.h:259
utl::lange< float >
float lange< float >(int matrix_order, char norm, int m, int n, const float *A, int LDA)
Definition: utlLapack.h:285
utl::geev< double >
int geev< double >(int matrix_layout, char jobvl, char jobvr, int n, double *a, int lda, double *wr, double *wi, double *vl, int ldvl, double *vr, int ldvr)
Definition: utlLapack.h:220
utl::sytri< double >
int sytri< double >(int matrix_order, char UPLO, int N, double *A, int LDA, const int *IPIV)
Definition: utlLapack.h:256
utl::sytrf< float >
int sytrf< float >(int matrix_order, char UPLO, int N, float *A, int LDA, int *IPIV)
Definition: utlLapack.h:272
utl::syevd
int syevd(int matrix_order, char JOBZ, char UPLO, int N, T *A, int LDA, T *W)
utl::gesdd
int gesdd(int matrix_order, char JOBZ, int M, int N, T *A, int LDA, T *S, T *U, int LDU, T *VT, int LDVT)
utlGlobalException
#define utlGlobalException(cond, expout)
Definition: utlCoreMacro.h:372
utl::getrf
int getrf(int matrix_layout, int m, int n, T *a, int lda, int *ipiv)
utl::syev< float >
int syev< float >(int matrix_order, char JOBZ, char UPLO, int N, float *A, int LDA, float *W)
Definition: utlLapack.h:203
utl
Definition: utl.h:90
LAPACKE_dlange
double LAPACKE_dlange(int matrix_order, char norm, lapack_int m, lapack_int n, const double *a, lapack_int lda)
utl::getrf< double >
int getrf< double >(int matrix_layout, int m, int n, double *a, int lda, int *ipiv)
Definition: utlLapack.h:295
utl::getrf< float >
int getrf< float >(int matrix_layout, int m, int n, float *a, int lda, int *ipiv)
Definition: utlLapack.h:298
utl::lange< double >
double lange< double >(int matrix_order, char norm, int m, int n, const double *A, int LDA)
Definition: utlLapack.h:282
utl::gesvd< double >
int gesvd< double >(int matrix_order, char JOBU, char JOBVT, int M, int N, double *A, int LDA, double *S, double *U, int LDU, double *VT, int LDVT, double *superb)
Definition: utlLapack.h:236
utl::gesvd
int gesvd(int matrix_order, char JOBU, char JOBVT, int M, int N, T *A, int LDA, T *S, T *U, int LDU, T *VT, int LDVT, T *superb)
utl::lange
T lange(int matrix_order, char norm, int m, int n, const T *A, int LDA)
LAPACKE_zlange
double LAPACKE_zlange(int matrix_order, char norm, lapack_int m, lapack_int n, const std::complex< double > *a, lapack_int lda)
utl::syevd< double >
int syevd< double >(int matrix_order, char JOBZ, char UPLO, int N, double *A, int LDA, double *W)
Definition: utlLapack.h:207
utl::gesdd< float >
int gesdd< float >(int matrix_order, char JOBZ, int M, int N, float *A, int LDA, float *S, float *U, int LDU, float *VT, int LDVT)
Definition: utlLapack.h:246
utl::getri< float >
int getri< float >(int matrix_layout, int n, float *a, int lda, const int *ipiv)
Definition: utlLapack.h:311
utl::gesdd< double >
int gesdd< double >(int matrix_order, char JOBZ, int M, int N, double *A, int LDA, double *S, double *U, int LDU, double *VT, int LDVT)
Definition: utlLapack.h:243
utl::syev
int syev(int matrix_order, char JOBZ, char UPLO, int N, T *A, int LDA, T *W)
utl::geev
int geev(int matrix_layout, char jobvl, char jobvr, int n, T *a, int lda, T *wr, T *wi, T *vl, int ldvl, T *vr, int ldvr)
utlCoreMacro.h
macros for utlCore
utl::syev< double >
int syev< double >(int matrix_order, char JOBZ, char UPLO, int N, double *A, int LDA, double *W)
Definition: utlLapack.h:200
utl::geev< float >
int geev< float >(int matrix_layout, char jobvl, char jobvr, int n, float *a, int lda, float *wr, float *wi, float *vl, int ldvl, float *vr, int ldvr)
Definition: utlLapack.h:223
LAPACKE_clange
float LAPACKE_clange(int matrix_order, char norm, lapack_int m, lapack_int n, const std::complex< float > *a, lapack_int lda)
Generated by   doxygen 1.8.11