- Timestamp:
- Jan 3, 2006, 3:34:15 PM (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/SophyaExt/LinAlg/intflapack.cc
r2646 r2875 80 80 */ 81 81 82 /* 83 Decembre 2005 : Suite portage AIX xlC 84 On declare des noms en majuscule pour les routines fortran - 85 avec ou sans underscore _ , suivant les systemes 86 */ 87 #ifdef AIX 88 89 #define sgesv sgesv 90 #define dgesv dgesv 91 #define cgesv cgesv 92 #define zgesv zgesv 93 94 #define ssysv ssysv 95 #define dsysv dsysv 96 #define csysv csysv 97 #define zsysv zsysv 98 99 #define sgels sgels 100 #define dgels dgels 101 #define cgels cgels 102 #define zgels zgels 103 104 #define sgelsd sgelsd 105 #define dgelsd dgelsd 106 #define cgelsd cgelsd 107 #define zgelsd zgelsd 108 109 #define sgesvd sgesvd 110 #define dgesvd dgesvd 111 #define cgesvd cgesvd 112 #define zgesvd zgesvd 113 114 #define sgesdd sgesdd 115 #define dgesdd dgesdd 116 #define cgesdd cgesdd 117 #define zgesdd zgesdd 118 119 #define ssyev ssyev 120 #define dsyev dsyev 121 #define cheev cheev 122 #define zheev zheev 123 124 #define sgeev sgeev 125 #define dgeev dgeev 126 #define cgeev cgeev 127 #define zgeev zgeev 128 129 #else 130 131 #define sgesv sgesv_ 132 #define dgesv dgesv_ 133 #define cgesv cgesv_ 134 #define zgesv zgesv_ 135 136 #define ssysv ssysv_ 137 #define dsysv dsysv_ 138 #define csysv csysv_ 139 #define zsysv zsysv_ 140 141 #define sgels sgels_ 142 #define dgels dgels_ 143 #define cgels cgels_ 144 #define zgels zgels_ 145 146 #define sgelsd sgelsd_ 147 #define dgelsd dgelsd_ 148 #define cgelsd cgelsd_ 149 #define zgelsd zgelsd_ 150 151 #define sgesvd sgesvd_ 152 #define dgesvd dgesvd_ 153 #define cgesvd cgesvd_ 154 #define zgesvd zgesvd_ 155 156 #define sgesdd sgesdd_ 157 #define dgesdd dgesdd_ 158 #define cgesdd cgesdd_ 159 #define zgesdd zgesdd_ 160 161 #define ssyev ssyev_ 162 #define dsyev dsyev_ 163 #define cheev cheev_ 164 #define zheev zheev_ 165 166 #define sgeev sgeev_ 167 #define dgeev dgeev_ 168 #define cgeev cgeev_ 169 #define zgeev zgeev_ 170 171 #endif 82 172 //////////////////////////////////////////////////////////////////////////////////// 83 173 extern "C" { 84 174 // Le calculateur de workingspace 85 int_4 ilaenv _(int_4 *ispec,char *name,char *opts,int_4 *n1,int_4 *n2,int_4 *n3,int_4 *n4,175 int_4 ilaenv(int_4 *ispec,char *name,char *opts,int_4 *n1,int_4 *n2,int_4 *n3,int_4 *n4, 86 176 int_4 nc1,int_4 nc2); 87 177 88 178 // Drivers pour resolution de systemes lineaires 89 void sgesv _(int_4* n, int_4* nrhs, r_4* a, int_4* lda,179 void sgesv(int_4* n, int_4* nrhs, r_4* a, int_4* lda, 90 180 int_4* ipiv, r_4* b, int_4* ldb, int_4* info); 91 void dgesv _(int_4* n, int_4* nrhs, r_8* a, int_4* lda,181 void dgesv(int_4* n, int_4* nrhs, r_8* a, int_4* lda, 92 182 int_4* ipiv, r_8* b, int_4* ldb, int_4* info); 93 void cgesv _(int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,183 void cgesv(int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 94 184 int_4* ipiv, complex<r_4>* b, int_4* ldb, int_4* info); 95 void zgesv _(int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,185 void zgesv(int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 96 186 int_4* ipiv, complex<r_8>* b, int_4* ldb, int_4* info); 97 187 98 188 // Drivers pour resolution de systemes lineaires symetriques 99 void ssysv _(char* uplo, int_4* n, int_4* nrhs, r_4* a, int_4* lda,189 void ssysv(char* uplo, int_4* n, int_4* nrhs, r_4* a, int_4* lda, 100 190 int_4* ipiv, r_4* b, int_4* ldb, 101 191 r_4* work, int_4* lwork, int_4* info); 102 void dsysv _(char* uplo, int_4* n, int_4* nrhs, r_8* a, int_4* lda,192 void dsysv(char* uplo, int_4* n, int_4* nrhs, r_8* a, int_4* lda, 103 193 int_4* ipiv, r_8* b, int_4* ldb, 104 194 r_8* work, int_4* lwork, int_4* info); 105 void csysv _(char* uplo, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,195 void csysv(char* uplo, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 106 196 int_4* ipiv, complex<r_4>* b, int_4* ldb, 107 197 complex<r_4>* work, int_4* lwork, int_4* info); 108 void zsysv _(char* uplo, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,198 void zsysv(char* uplo, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 109 199 int_4* ipiv, complex<r_8>* b, int_4* ldb, 110 200 complex<r_8>* work, int_4* lwork, int_4* info); 111 201 112 202 // Driver pour resolution de systemes au sens de Xi2 113 void sgels _(char * trans, int_4* m, int_4* n, int_4* nrhs, r_4* a, int_4* lda,203 void sgels(char * trans, int_4* m, int_4* n, int_4* nrhs, r_4* a, int_4* lda, 114 204 r_4* b, int_4* ldb, r_4* work, int_4* lwork, int_4* info); 115 void dgels _(char * trans, int_4* m, int_4* n, int_4* nrhs, r_8* a, int_4* lda,205 void dgels(char * trans, int_4* m, int_4* n, int_4* nrhs, r_8* a, int_4* lda, 116 206 r_8* b, int_4* ldb, r_8* work, int_4* lwork, int_4* info); 117 void cgels _(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda,207 void cgels(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 118 208 complex<r_4>* b, int_4* ldb, complex<r_4>* work, int_4* lwork, int_4* info); 119 void zgels _(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda,209 void zgels(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 120 210 complex<r_8>* b, int_4* ldb, complex<r_8>* work, int_4* lwork, int_4* info); 121 211 122 212 // Driver pour resolution de systemes au sens de Xi2 par SVD Divide & Conquer 123 void sgelsd _(int_4* m,int_4* n,int_4* nrhs,r_4* a,int_4* lda,213 void sgelsd(int_4* m,int_4* n,int_4* nrhs,r_4* a,int_4* lda, 124 214 r_4* b,int_4* ldb,r_4* s,r_4* rcond,int_4* rank, 125 215 r_4* work,int_4* lwork,int_4* iwork,int_4* info); 126 void dgelsd _(int_4* m,int_4* n,int_4* nrhs,r_8* a,int_4* lda,216 void dgelsd(int_4* m,int_4* n,int_4* nrhs,r_8* a,int_4* lda, 127 217 r_8* b,int_4* ldb,r_8* s,r_8* rcond,int_4* rank, 128 218 r_8* work,int_4* lwork,int_4* iwork,int_4* info); 129 void cgelsd _(int_4* m,int_4* n,int_4* nrhs,complex<r_4>* a,int_4* lda,219 void cgelsd(int_4* m,int_4* n,int_4* nrhs,complex<r_4>* a,int_4* lda, 130 220 complex<r_4>* b,int_4* ldb,r_4* s,r_4* rcond,int_4* rank, 131 221 complex<r_4>* work,int_4* lwork,r_4* rwork,int_4* iwork,int_4* info); 132 void zgelsd _(int_4* m,int_4* n,int_4* nrhs,complex<r_8>* a,int_4* lda,222 void zgelsd(int_4* m,int_4* n,int_4* nrhs,complex<r_8>* a,int_4* lda, 133 223 complex<r_8>* b,int_4* ldb,r_8* s,r_8* rcond,int_4* rank, 134 224 complex<r_8>* work,int_4* lwork,r_8* rwork,int_4* iwork,int_4* info); 135 225 136 226 // Driver pour decomposition SVD 137 void sgesvd _(char* jobu, char* jobvt, int_4* m, int_4* n, r_4* a, int_4* lda,227 void sgesvd(char* jobu, char* jobvt, int_4* m, int_4* n, r_4* a, int_4* lda, 138 228 r_4* s, r_4* u, int_4* ldu, r_4* vt, int_4* ldvt, 139 229 r_4* work, int_4* lwork, int_4* info); 140 void dgesvd _(char* jobu, char* jobvt, int_4* m, int_4* n, r_8* a, int_4* lda,230 void dgesvd(char* jobu, char* jobvt, int_4* m, int_4* n, r_8* a, int_4* lda, 141 231 r_8* s, r_8* u, int_4* ldu, r_8* vt, int_4* ldvt, 142 232 r_8* work, int_4* lwork, int_4* info); 143 void cgesvd _(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,233 void cgesvd(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_4>* a, int_4* lda, 144 234 r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt, 145 235 complex<r_4>* work, int_4* lwork, r_4* rwork, int_4* info); 146 void zgesvd _(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,236 void zgesvd(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_8>* a, int_4* lda, 147 237 r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt, 148 238 complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* info); 149 239 150 240 // Driver pour decomposition SVD Divide and Conquer 151 void sgesdd _(char* jobz, int_4* m, int_4* n, r_4* a, int_4* lda,241 void sgesdd(char* jobz, int_4* m, int_4* n, r_4* a, int_4* lda, 152 242 r_4* s, r_4* u, int_4* ldu, r_4* vt, int_4* ldvt, 153 243 r_4* work, int_4* lwork, int_4* iwork, int_4* info); 154 void dgesdd _(char* jobz, int_4* m, int_4* n, r_8* a, int_4* lda,244 void dgesdd(char* jobz, int_4* m, int_4* n, r_8* a, int_4* lda, 155 245 r_8* s, r_8* u, int_4* ldu, r_8* vt, int_4* ldvt, 156 246 r_8* work, int_4* lwork, int_4* iwork, int_4* info); 157 void cgesdd _(char* jobz, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,247 void cgesdd(char* jobz, int_4* m, int_4* n, complex<r_4>* a, int_4* lda, 158 248 r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt, 159 249 complex<r_4>* work, int_4* lwork, r_4* rwork, int_4* iwork, int_4* info); 160 void zgesdd _(char* jobz, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,250 void zgesdd(char* jobz, int_4* m, int_4* n, complex<r_8>* a, int_4* lda, 161 251 r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt, 162 252 complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* iwork, int_4* info); 163 253 164 254 // Driver pour eigen decomposition for symetric/hermitian matrices 165 void ssyev _(char* jobz, char* uplo, int_4* n, r_4* a, int_4* lda, r_4* w,255 void ssyev(char* jobz, char* uplo, int_4* n, r_4* a, int_4* lda, r_4* w, 166 256 r_4* work, int_4 *lwork, int_4* info); 167 void dsyev _(char* jobz, char* uplo, int_4* n, r_8* a, int_4* lda, r_8* w,257 void dsyev(char* jobz, char* uplo, int_4* n, r_8* a, int_4* lda, r_8* w, 168 258 r_8* work, int_4 *lwork, int_4* info); 169 void cheev _(char* jobz, char* uplo, int_4* n, complex<r_4>* a, int_4* lda, r_4* w,259 void cheev(char* jobz, char* uplo, int_4* n, complex<r_4>* a, int_4* lda, r_4* w, 170 260 complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info); 171 void zheev _(char* jobz, char* uplo, int_4* n, complex<r_8>* a, int_4* lda, r_8* w,261 void zheev(char* jobz, char* uplo, int_4* n, complex<r_8>* a, int_4* lda, r_8* w, 172 262 complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info); 173 263 174 264 // Driver pour eigen decomposition for general squared matrices 175 void sgeev _(char* jobl, char* jobvr, int_4* n, r_4* a, int_4* lda, r_4* wr, r_4* wi,265 void sgeev(char* jobl, char* jobvr, int_4* n, r_4* a, int_4* lda, r_4* wr, r_4* wi, 176 266 r_4* vl, int_4* ldvl, r_4* vr, int_4* ldvr, 177 267 r_4* work, int_4 *lwork, int_4* info); 178 void dgeev _(char* jobl, char* jobvr, int_4* n, r_8* a, int_4* lda, r_8* wr, r_8* wi,268 void dgeev(char* jobl, char* jobvr, int_4* n, r_8* a, int_4* lda, r_8* wr, r_8* wi, 179 269 r_8* vl, int_4* ldvl, r_8* vr, int_4* ldvr, 180 270 r_8* work, int_4 *lwork, int_4* info); 181 void cgeev _(char* jobl, char* jobvr, int_4* n, complex<r_4>* a, int_4* lda, complex<r_4>* w,271 void cgeev(char* jobl, char* jobvr, int_4* n, complex<r_4>* a, int_4* lda, complex<r_4>* w, 182 272 complex<r_4>* vl, int_4* ldvl, complex<r_4>* vr, int_4* ldvr, 183 273 complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info); 184 void zgeev _(char* jobl, char* jobvr, int_4* n, complex<r_8>* a, int_4* lda, complex<r_8>* w,274 void zgeev(char* jobl, char* jobvr, int_4* n, complex<r_8>* a, int_4* lda, complex<r_8>* w, 185 275 complex<r_8>* vl, int_4* ldvl, complex<r_8>* vr, int_4* ldvr, 186 276 complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info); … … 207 297 { 208 298 int_4 nc1 = strlen(name), nc2 = strlen(opts), rc=0; 209 rc = ilaenv _(&ispec,name,opts,&n1,&n2,&n3,&n4,nc1,nc2);299 rc = ilaenv(&ispec,name,opts,&n1,&n2,&n3,&n4,nc1,nc2); 210 300 //cout<<"ilaenv_en_C("<<ispec<<","<<name<<"("<<nc1<<"),"<<opts<<"("<<nc2<<")," 211 301 // <<n1<<","<<n2<<","<<n3<<","<<n4<<") = "<<rc<<endl; … … 262 352 263 353 if (typeid(T) == typeid(r_4) ) 264 sgesv _(&n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb, &info);354 sgesv(&n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb, &info); 265 355 else if (typeid(T) == typeid(r_8) ) 266 dgesv _(&n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb, &info);356 dgesv(&n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb, &info); 267 357 else if (typeid(T) == typeid(complex<r_4>) ) 268 cgesv _(&n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,358 cgesv(&n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv, 269 359 (complex<r_4> *)b.Data(), &ldb, &info); 270 360 else if (typeid(T) == typeid(complex<r_8>) ) 271 zgesv _(&n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,361 zgesv(&n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv, 272 362 (complex<r_8> *)b.Data(), &ldb, &info); 273 363 else { … … 326 416 327 417 if (typeid(T) == typeid(r_4) ) { 328 ssysv _(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,418 ssysv(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb, 329 419 (r_4 *)wkget, &lwork, &info); 330 420 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 331 ssysv _(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,421 ssysv(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb, 332 422 (r_4 *)work, &lwork, &info); 333 423 } else if (typeid(T) == typeid(r_8) ) { 334 dsysv _(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,424 dsysv(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb, 335 425 (r_8 *)wkget, &lwork, &info); 336 426 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 337 dsysv _(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,427 dsysv(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb, 338 428 (r_8 *)work, &lwork, &info); 339 429 } else if (typeid(T) == typeid(complex<r_4>) ) { 340 csysv _(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,430 csysv(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv, 341 431 (complex<r_4> *)b.Data(), &ldb, 342 432 (complex<r_4> *)wkget, &lwork, &info); 343 433 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 344 csysv _(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,434 csysv(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv, 345 435 (complex<r_4> *)b.Data(), &ldb, 346 436 (complex<r_4> *)work, &lwork, &info); 347 437 } else if (typeid(T) == typeid(complex<r_8>) ) { 348 zsysv _(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,438 zsysv(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv, 349 439 (complex<r_8> *)b.Data(), &ldb, 350 440 (complex<r_8> *)wkget, &lwork, &info); 351 441 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 352 zsysv _(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,442 zsysv(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv, 353 443 (complex<r_8> *)b.Data(), &ldb, 354 444 (complex<r_8> *)work, &lwork, &info); … … 427 517 428 518 if (typeid(T) == typeid(r_4) ) { 429 sgels _(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda,519 sgels(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda, 430 520 (r_4 *)b.Data(), &ldb, (r_4 *)wkget, &lwork, &info); 431 521 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 432 sgels _(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda,522 sgels(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda, 433 523 (r_4 *)b.Data(), &ldb, (r_4 *)work, &lwork, &info); 434 524 } else if (typeid(T) == typeid(r_8) ) { 435 dgels _(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda,525 dgels(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda, 436 526 (r_8 *)b.Data(), &ldb, (r_8 *)wkget, &lwork, &info); 437 527 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 438 dgels _(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda,528 dgels(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda, 439 529 (r_8 *)b.Data(), &ldb, (r_8 *)work, &lwork, &info); 440 530 } else if (typeid(T) == typeid(complex<r_4>) ) { 441 cgels _(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda,531 cgels(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, 442 532 (complex<r_4> *)b.Data(), &ldb, (complex<r_4> *)wkget, &lwork, &info); 443 533 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 444 cgels _(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda,534 cgels(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, 445 535 (complex<r_4> *)b.Data(), &ldb, (complex<r_4> *)work, &lwork, &info); 446 536 } else if (typeid(T) == typeid(complex<r_8>) ) { 447 zgels _(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda,537 zgels(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, 448 538 (complex<r_8> *)b.Data(), &ldb, (complex<r_8> *)wkget, &lwork, &info); 449 539 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 450 zgels _(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda,540 zgels(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, 451 541 (complex<r_8> *)b.Data(), &ldb, (complex<r_8> *)work, &lwork, &info); 452 542 } else { … … 554 644 r_4 srcond = rcond; 555 645 iwork = new int_4[3*minmn*nlvl+11*minmn +GARDMEM]; 556 sgelsd _(&m,&n,&nrhs,(r_4*)a.Data(),&lda,646 sgelsd(&m,&n,&nrhs,(r_4*)a.Data(),&lda, 557 647 (r_4*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank, 558 648 (r_4*)wkget,&lwork,(int_4*)iwork,&info); 559 649 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 560 sgelsd _(&m,&n,&nrhs,(r_4*)a.Data(),&lda,650 sgelsd(&m,&n,&nrhs,(r_4*)a.Data(),&lda, 561 651 (r_4*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank, 562 652 (r_4*)work,&lwork,(int_4*)iwork,&info); … … 565 655 } else if(typeid(T) == typeid(r_8) ) { 566 656 iwork = new int_4[3*minmn*nlvl+11*minmn +GARDMEM]; 567 dgelsd _(&m,&n,&nrhs,(r_8*)a.Data(),&lda,657 dgelsd(&m,&n,&nrhs,(r_8*)a.Data(),&lda, 568 658 (r_8*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank, 569 659 (r_8*)wkget,&lwork,(int_4*)iwork,&info); 570 660 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 571 dgelsd _(&m,&n,&nrhs,(r_8*)a.Data(),&lda,661 dgelsd(&m,&n,&nrhs,(r_8*)a.Data(),&lda, 572 662 (r_8*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank, 573 663 (r_8*)work,&lwork,(int_4*)iwork,&info); … … 581 671 r_4* sloc = new r_4[minmn]; 582 672 r_4 srcond = rcond; 583 cgelsd _(&m,&n,&nrhs,(complex<r_4>*)a.Data(),&lda,673 cgelsd(&m,&n,&nrhs,(complex<r_4>*)a.Data(),&lda, 584 674 (complex<r_4>*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank, 585 675 (complex<r_4>*)wkget,&lwork,(r_4*)rwork,(int_4*)iwork,&info); 586 676 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 587 cgelsd _(&m,&n,&nrhs,(complex<r_4>*)a.Data(),&lda,677 cgelsd(&m,&n,&nrhs,(complex<r_4>*)a.Data(),&lda, 588 678 (complex<r_4>*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank, 589 679 (complex<r_4>*)work,&lwork,(r_4*)rwork,(int_4*)iwork,&info); … … 598 688 r_8* rwork = new r_8[lrwork +GARDMEM]; 599 689 iwork = new int_4[3*minmn*nlvl+11*minmn +GARDMEM]; 600 zgelsd _(&m,&n,&nrhs,(complex<r_8>*)a.Data(),&lda,690 zgelsd(&m,&n,&nrhs,(complex<r_8>*)a.Data(),&lda, 601 691 (complex<r_8>*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank, 602 692 (complex<r_8>*)wkget,&lwork,(r_8*)rwork,(int_4*)iwork,&info); 603 693 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 604 zgelsd _(&m,&n,&nrhs,(complex<r_8>*)a.Data(),&lda,694 zgelsd(&m,&n,&nrhs,(complex<r_8>*)a.Data(),&lda, 605 695 (complex<r_8>*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank, 606 696 (complex<r_8>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info); … … 726 816 727 817 if (typeid(T) == typeid(r_4) ) { 728 sgesvd _(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda,818 sgesvd(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda, 729 819 (r_4 *)s.Data(), (r_4 *) up->Data(), &ldu, (r_4 *)vtp->Data(), &ldvt, 730 820 (r_4 *)wkget, &lwork, &info); 731 821 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 732 sgesvd _(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda,822 sgesvd(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda, 733 823 (r_4 *)s.Data(), (r_4 *) up->Data(), &ldu, (r_4 *)vtp->Data(), &ldvt, 734 824 (r_4 *)work, &lwork, &info); 735 825 } else if (typeid(T) == typeid(r_8) ) { 736 dgesvd _(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda,826 dgesvd(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda, 737 827 (r_8 *)s.Data(), (r_8 *) up->Data(), &ldu, (r_8 *)vtp->Data(), &ldvt, 738 828 (r_8 *)wkget, &lwork, &info); 739 829 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 740 dgesvd _(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda,830 dgesvd(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda, 741 831 (r_8 *)s.Data(), (r_8 *) up->Data(), &ldu, (r_8 *)vtp->Data(), &ldvt, 742 832 (r_8 *)work, &lwork, &info); … … 744 834 r_4 * rwork = new r_4[5*minmn +GARDMEM]; 745 835 r_4 * sloc = new r_4[minmn]; 746 cgesvd _(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,836 cgesvd(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda, 747 837 (r_4 *)sloc, (complex<r_4> *) up->Data(), &ldu, 748 838 (complex<r_4> *)vtp->Data(), &ldvt, 749 839 (complex<r_4> *)wkget, &lwork, (r_4 *)rwork, &info); 750 840 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 751 cgesvd _(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,841 cgesvd(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda, 752 842 (r_4 *)sloc, (complex<r_4> *) up->Data(), &ldu, 753 843 (complex<r_4> *)vtp->Data(), &ldvt, … … 758 848 r_8 * rwork = new r_8[5*minmn +GARDMEM]; 759 849 r_8 * sloc = new r_8[minmn]; 760 zgesvd _(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,850 zgesvd(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda, 761 851 (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu, 762 852 (complex<r_8> *)vtp->Data(), &ldvt, 763 853 (complex<r_8> *)wkget, &lwork, (r_8 *)rwork, &info); 764 854 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 765 zgesvd _(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,855 zgesvd(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda, 766 856 (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu, 767 857 (complex<r_8> *)vtp->Data(), &ldvt, … … 818 908 r_4* sloc = new r_4[minmn]; 819 909 iwork = new int_4[8*minmn +GARDMEM]; 820 sgesdd _(&jobz,&m,&n,(r_4*)a.Data(),&lda,910 sgesdd(&jobz,&m,&n,(r_4*)a.Data(),&lda, 821 911 (r_4*)sloc,(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt, 822 912 (r_4*)wkget,&lwork,(int_4*)iwork,&info); 823 913 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 824 sgesdd _(&jobz,&m,&n,(r_4*)a.Data(),&lda,914 sgesdd(&jobz,&m,&n,(r_4*)a.Data(),&lda, 825 915 (r_4*)sloc,(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt, 826 916 (r_4*)work,&lwork,(int_4*)iwork,&info); … … 829 919 } else if(typeid(T) == typeid(r_8) ) { 830 920 iwork = new int_4[8*minmn +GARDMEM]; 831 dgesdd _(&jobz,&m,&n,(r_8*)a.Data(),&lda,921 dgesdd(&jobz,&m,&n,(r_8*)a.Data(),&lda, 832 922 (r_8*)s.Data(),(r_8*)u.Data(),&ldu,(r_8*)vt.Data(),&ldvt, 833 923 (r_8*)wkget,&lwork,(int_4*)iwork,&info); 834 924 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 835 dgesdd _(&jobz,&m,&n,(r_8*)a.Data(),&lda,925 dgesdd(&jobz,&m,&n,(r_8*)a.Data(),&lda, 836 926 (r_8*)s.Data(),(r_8*)u.Data(),&ldu,(r_8*)vt.Data(),&ldvt, 837 927 (r_8*)work,&lwork,(int_4*)iwork,&info); … … 840 930 r_4* rwork = new r_4[5*minmn*minmn+5*minmn +GARDMEM]; 841 931 iwork = new int_4[8*minmn +GARDMEM]; 842 cgesdd _(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda,932 cgesdd(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda, 843 933 (r_4*)sloc,(complex<r_4>*)u.Data(),&ldu,(complex<r_4>*)vt.Data(),&ldvt, 844 934 (complex<r_4>*)wkget,&lwork,(r_4*)rwork,(int_4*)iwork,&info); 845 935 lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM]; 846 cgesdd _(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda,936 cgesdd(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda, 847 937 (r_4*)sloc,(complex<r_4>*)u.Data(),&ldu,(complex<r_4>*)vt.Data(),&ldvt, 848 938 (complex<r_4>*)work,&lwork,(r_4*)rwork,(int_4*)iwork,&info); … … 852 942 r_8* rwork = new r_8[5*minmn*minmn+5*minmn +GARDMEM]; 853 943 iwork = new int_4[8*minmn +GARDMEM]; 854 zgesdd _(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda,944 zgesdd(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda, 855 945 (r_8*)s.Data(),(complex<r_8>*)u.Data(),&ldu,(complex<r_8>*)vt.Data(),&ldvt, 856 946 (complex<r_8>*)wkget,&lwork,(r_8*)rwork,(int_4*)iwork,&info); 857 947 lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM]; 858 zgesdd _(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda,948 zgesdd(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda, 859 949 (r_8*)s.Data(),(complex<r_8>*)u.Data(),&ldu,(complex<r_8>*)vt.Data(),&ldvt, 860 950 (complex<r_8>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info); … … 909 999 if (typeid(T) == typeid(r_4) ) { 910 1000 r_4* w = new r_4[n]; 911 ssyev _(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)wkget,&lwork,&info);1001 ssyev(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)wkget,&lwork,&info); 912 1002 lwork = type2i4(&wkget[0],4); /* 3*n-1;*/ work = new T[lwork +GARDMEM]; 913 ssyev _(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info);1003 ssyev(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info); 914 1004 if(info==0) for(int i=0;i<n;i++) b(i) = w[i]; 915 1005 delete [] w; 916 1006 } else if (typeid(T) == typeid(r_8) ) { 917 1007 r_8* w = new r_8[n]; 918 dsyev _(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)wkget,&lwork,&info);1008 dsyev(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)wkget,&lwork,&info); 919 1009 lwork = type2i4(&wkget[0],8); /* 3*n-1;*/ work = new T[lwork +GARDMEM]; 920 dsyev _(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info);1010 dsyev(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info); 921 1011 if(info==0) for(int i=0;i<n;i++) b(i) = w[i]; 922 1012 delete [] w; 923 1013 } else if (typeid(T) == typeid(complex<r_4>) ) { 924 1014 r_4* rwork = new r_4[3*n-2 +GARDMEM]; r_4* w = new r_4[n]; 925 cheev _(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w1015 cheev(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w 926 1016 ,(complex<r_4> *)wkget,&lwork,(r_4 *)rwork,&info); 927 1017 lwork = type2i4(&wkget[0],4); /* 2*n-1;*/ work = new T[lwork +GARDMEM]; 928 cheev _(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w1018 cheev(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w 929 1019 ,(complex<r_4> *)work,&lwork,(r_4 *)rwork,&info); 930 1020 if(info==0) for(int i=0;i<n;i++) b(i) = w[i]; … … 932 1022 } else if (typeid(T) == typeid(complex<r_8>) ) { 933 1023 r_8* rwork = new r_8[3*n-2 +GARDMEM]; r_8* w = new r_8[n]; 934 zheev _(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w1024 zheev(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w 935 1025 ,(complex<r_8> *)wkget,&lwork,(r_8 *)rwork,&info); 936 1026 lwork = type2i4(&wkget[0],8); /* 2*n-1;*/ work = new T[lwork +GARDMEM]; 937 zheev _(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w1027 zheev(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w 938 1028 ,(complex<r_8> *)work,&lwork,(r_8 *)rwork,&info); 939 1029 if(info==0) for(int i=0;i<n;i++) b(i) = w[i]; … … 1004 1094 if (typeid(T) == typeid(r_4) ) { 1005 1095 r_4* wr = new r_4[n]; r_4* wi = new r_4[n]; r_4* vl = NULL; 1006 sgeev _(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,1096 sgeev(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi, 1007 1097 (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr, 1008 1098 (r_4 *)wkget,&lwork,&info); 1009 1099 lwork = type2i4(&wkget[0],4); /* 4*n;*/ work = new T[lwork +GARDMEM]; 1010 sgeev _(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,1100 sgeev(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi, 1011 1101 (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr, 1012 1102 (r_4 *)work,&lwork,&info); … … 1015 1105 } else if (typeid(T) == typeid(r_8) ) { 1016 1106 r_8* wr = new r_8[n]; r_8* wi = new r_8[n]; r_8* vl = NULL; 1017 dgeev _(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,1107 dgeev(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi, 1018 1108 (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr, 1019 1109 (r_8 *)wkget,&lwork,&info); 1020 1110 lwork = type2i4(&wkget[0],8); /* 4*n;*/ work = new T[lwork +GARDMEM]; 1021 dgeev _(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,1111 dgeev(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi, 1022 1112 (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr, 1023 1113 (r_8 *)work,&lwork,&info); … … 1026 1116 } else if (typeid(T) == typeid(complex<r_4>) ) { 1027 1117 r_4* rwork = new r_4[2*n +GARDMEM]; r_4* vl = NULL; TVector< complex<r_4> > w(n); 1028 cgeev _(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),1118 cgeev(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(), 1029 1119 (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr, 1030 1120 (complex<r_4> *)wkget,&lwork,(r_4 *)rwork,&info); 1031 1121 lwork = type2i4(&wkget[0],4); /* 2*n;*/ work = new T[lwork +GARDMEM]; 1032 cgeev _(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),1122 cgeev(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(), 1033 1123 (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr, 1034 1124 (complex<r_4> *)work,&lwork,(r_4 *)rwork,&info); … … 1037 1127 } else if (typeid(T) == typeid(complex<r_8>) ) { 1038 1128 r_8* rwork = new r_8[2*n +GARDMEM]; r_8* vl = NULL; 1039 zgeev _(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),1129 zgeev(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(), 1040 1130 (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr, 1041 1131 (complex<r_8> *)wkget,&lwork,(r_8 *)rwork,&info); 1042 1132 lwork = type2i4(&wkget[0],8); /* 2*n;*/ work = new T[lwork +GARDMEM]; 1043 zgeev _(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),1133 zgeev(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(), 1044 1134 (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr, 1045 1135 (complex<r_8> *)work,&lwork,(r_8 *)rwork,&info); … … 1068 1158 1069 1159 #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES) 1160 namespace SOPHYA { 1070 1161 template class LapackServer<r_4>; 1071 1162 template class LapackServer<r_8>; 1072 1163 template class LapackServer< complex<r_4> >; 1073 1164 template class LapackServer< complex<r_8> >; 1165 } 1074 1166 #endif 1075 1167
Note:
See TracChangeset
for help on using the changeset viewer.