NGSolve
4.9
|
00001 #ifndef FILE_NG_LAPACK 00002 #define FILE_NG_LAPACK 00003 00004 // #include <mkl_cblas.h> 00005 00006 /****************************************************************************/ 00007 /* File: ng_lapack.hpp */ 00008 /* Author: Joachim Schoeberl */ 00009 /* Date: 21. Nov. 2004 */ 00010 /****************************************************************************/ 00011 00012 00013 00014 namespace ngbla 00015 { 00016 // Interface to lapack functions 00017 00018 00019 #ifdef LAPACK 00020 00021 extern "C" { 00022 #ifdef MKL_ILP64 00023 typedef long int integer; 00024 #else 00025 typedef int integer; 00026 #endif 00027 // typedef char logical; 00028 typedef integer logical; 00029 typedef float real; 00030 typedef double doublereal; 00031 typedef Complex doublecomplex; 00032 typedef complex<float> singlecomplex; 00033 00034 00035 // Windows SDK defines VOID in the file WinNT.h 00036 #ifndef VOID 00037 typedef void VOID; 00038 #endif 00039 00040 typedef int ftnlen; 00041 typedef int L_fp; // ? 00042 00043 00044 #include "clapack.h" 00045 } 00046 00047 00048 00049 // BLAS 1 00050 00051 inline double LapackDot (FlatVector<double> x, FlatVector<double> y) 00052 { 00053 integer n = x.Size(); 00054 integer incx = 1; 00055 integer incy = 1; 00056 return ddot_ (&n, &x(0), &incx, &y(0), &incy); 00057 } 00058 00059 00060 00061 // BLAS 2 00062 00063 /* 00064 extern "C" 00065 void dgemv_ (char & trans, int & m, int & n, double & alpha, double & a, int & lda, double & x, int & incx, 00066 double & beta, double & y, int & incy); 00067 */ 00068 00069 inline void LapackMultAx (ngbla::FlatMatrix<double> a, 00070 ngbla::FlatVector<double> x, 00071 ngbla::FlatVector<double> y) 00072 { 00073 char trans = 'T'; 00074 integer m = a.Width(); 00075 integer n = a.Height(); 00076 double alpha = 1; 00077 integer lda = a.Width(); 00078 integer incx = 1; 00079 double beta = 0; 00080 integer incy = 1; 00081 dgemv_ (&trans, &m, &n, &alpha, &a(0,0), &lda, &x(0), &incx, &beta, &y(0), &incy); 00082 } 00083 00084 inline void LapackMultAx (ngbla::FlatMatrix<Complex> a, 00085 ngbla::FlatVector<Complex> x, 00086 ngbla::FlatVector<Complex> y) 00087 { 00088 char trans = 'T'; 00089 integer m = a.Width(); 00090 integer n = a.Height(); 00091 Complex alpha(1,0); 00092 integer lda = a.Width(); 00093 integer incx = 1; 00094 Complex beta(0, 0); 00095 integer incy = 1; 00096 zgemv_ (&trans, &m, &n, &alpha, &a(0,0), &lda, &x(0), &incx, &beta, &y(0), &incy); 00097 } 00098 00099 00100 00101 inline void LapackMultAtx (ngbla::FlatMatrix<double> a, 00102 ngbla::FlatVector<double> x, 00103 ngbla::FlatVector<double> y) 00104 { 00105 char trans = 'N'; 00106 integer m = a.Width(); 00107 integer n = a.Height(); 00108 double alpha = 1; 00109 integer lda = a.Width(); 00110 integer incx = 1; 00111 double beta = 0; 00112 integer incy = 1; 00113 dgemv_ (&trans, &m, &n, &alpha, &a(0,0), &lda, &x(0), &incx, &beta, &y(0), &incy); 00114 } 00115 00116 00117 /* 00118 extern "C" 00119 void dger_ (int & m, int & n, double & alpha, double & x, int & incx, double & y, int & incy, double & a, int & lda); 00120 */ 00121 00122 inline void LapackAddxyt (ngbla::FlatMatrix<double> a, 00123 double fac, 00124 ngbla::FlatVector<double> x, 00125 ngbla::FlatVector<double> y) 00126 { 00127 integer m = a.Width(); 00128 integer n = a.Height(); 00129 double alpha = fac; 00130 integer lda = a.Width(); 00131 integer incx = 1; 00132 integer incy = 1; 00133 00134 dger_ (&m, &n, &alpha, &y(0), &incx, &x(0), &incy, &a(0,0), &lda); 00135 } 00136 00137 00138 00139 inline int gemm(char *transa, char *transb, integer *m, integer * 00140 n, integer *k, doublereal *alpha, doublereal *a, integer *lda, 00141 doublereal *b, integer *ldb, doublereal *beta, doublereal *c__, 00142 integer *ldc) 00143 { 00144 return dgemm_ (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c__, ldc); 00145 } 00146 00147 inline int gemm(char *transa, char *transb, integer *m, integer * 00148 n, integer *k, doublecomplex *alpha, doublecomplex *a, integer *lda, 00149 doublecomplex *b, integer *ldb, doublecomplex *beta, doublecomplex * 00150 c__, integer *ldc) 00151 { 00152 return zgemm_ (transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c__, ldc); 00153 } 00154 00155 00156 template <typename SCAL> 00157 inline void BASE_LapackMult (SliceMatrix<SCAL> a, bool transa, 00158 SliceMatrix<SCAL> b, bool transb, 00159 SliceMatrix<SCAL> c) 00160 { 00161 char transa_ = transa ? 'T' : 'N'; 00162 char transb_ = transb ? 'T' : 'N'; 00163 00164 integer n = c.Width(); 00165 integer m = c.Height(); 00166 if (n == 0 || m == 0) return; 00167 integer k = transa ? a.Height() : a.Width(); 00168 SCAL alpha = 1.0; 00169 SCAL beta = 0; 00170 integer lda = a.Dist(); 00171 integer ldb = b.Dist(); 00172 integer ldc = c.Dist(); 00173 00174 gemm (&transb_, &transa_, &n, &m, &k, &alpha, 00175 &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00176 } 00177 00178 template <typename SCAL> 00179 inline void BASE_LapackMultAdd (SliceMatrix<SCAL> a, bool transa, 00180 SliceMatrix<SCAL> b, bool transb, 00181 SCAL aalpha, 00182 SliceMatrix<SCAL> c, 00183 SCAL abeta) 00184 { 00185 char transa_ = transa ? 'T' : 'N'; 00186 char transb_ = transb ? 'T' : 'N'; 00187 00188 integer n = c.Width(); 00189 integer m = c.Height(); 00190 if (n == 0 || m == 0) return; 00191 integer k = transa ? a.Height() : a.Width(); 00192 SCAL alpha = aalpha; 00193 SCAL beta = abeta; 00194 integer lda = a.Dist(); 00195 integer ldb = b.Dist(); 00196 integer ldc = c.Dist(); 00197 00198 gemm (&transb_, &transa_, &n, &m, &k, &alpha, 00199 &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00200 } 00201 00202 00203 00204 template <typename TA, typename TB, typename TC> 00205 inline void LapackMult (const TA & a, 00206 const TB & b, 00207 const TC & c) 00208 { BASE_LapackMult<typename TC::TSCAL> (a, false, b, false, c); } 00209 00210 template <typename TA, typename TB, typename TC> 00211 inline void LapackMult (const TA & a, 00212 TransExpr<TB> b, 00213 const TC & c) 00214 { BASE_LapackMult<typename TC::TSCAL> (a, false, b.A(), true, c); } 00215 00216 template <typename TA, typename TB, typename TC> 00217 inline void LapackMult (TransExpr<TA> a, 00218 const TB & b, 00219 const TC & c) 00220 { BASE_LapackMult<typename TC::TSCAL> (a.A(), true, b, false, c); } 00221 00222 template <typename TA, typename TB, typename TC> 00223 inline void LapackMult (TransExpr<TA> a, 00224 TransExpr<TB> b, 00225 const TC & c) 00226 { BASE_LapackMult<typename TC::TSCAL> (a.A(), true, b.A(), true, c); } 00227 00228 00229 00230 /* 00231 inline void LapackMult (SliceMatrix<double> a, 00232 SliceMatrix<double> b, 00233 SliceMatrix<double> c) 00234 { BASE_LapackMult<double> (a, false, b, false, c); } 00235 00236 template <typename TA> 00237 inline void LapackMult (SliceMatrix<double> a, 00238 TransExpr<TA> b, 00239 SliceMatrix<double> c) 00240 { BASE_LapackMult<double> (a, false, b.A(), true, c); } 00241 00242 template <typename TA> 00243 inline void LapackMult (TransExpr<TA> a, 00244 SliceMatrix<double> b, 00245 SliceMatrix<double> c) 00246 { BASE_LapackMult<double> (a.A(), true, b, false, c); } 00247 00248 template <typename TA, typename TB> 00249 inline void LapackMult (TransExpr<TA> a, 00250 TransExpr<TB> b, 00251 SliceMatrix<double> c) 00252 { BASE_LapackMult<double> (a.A(), true, b.A(), true, c); } 00253 */ 00254 00255 00256 /* 00257 inline void LapackMult (SliceMatrix<Complex> a, 00258 SliceMatrix<Complex> b, 00259 SliceMatrix<Complex> c) 00260 { BASE_LapackMult<Complex> (a, false, b, false, c); } 00261 00262 template <typename TA> 00263 inline void LapackMult (SliceMatrix<Complex> a, 00264 TransExpr<TA> b, 00265 SliceMatrix<Complex> c) 00266 { BASE_LapackMult<Complex> (a, false, b.A(), true, c); } 00267 00268 template <typename TA> 00269 inline void LapackMult (TransExpr<TA> a, 00270 SliceMatrix<Complex> b, 00271 SliceMatrix<Complex> c) 00272 { BASE_LapackMult<Complex> (a.A(), true, b, false, c); } 00273 00274 template <typename TA, typename TB> 00275 inline void LapackMult (TransExpr<TA> a, 00276 TransExpr<TB> b, 00277 SliceMatrix<Complex> c) 00278 { BASE_LapackMult<Complex> (a.A(), true, b.A(), true, c); } 00279 */ 00280 00281 00282 00283 00284 00285 inline void LapackMultAdd (SliceMatrix<double> a, 00286 SliceMatrix<double> b, 00287 double alpha, 00288 SliceMatrix<double> c, 00289 double beta) 00290 { BASE_LapackMultAdd<double> (a, false, b, false, alpha, c, beta); } 00291 00292 template <typename TA> 00293 inline void LapackMultAdd (SliceMatrix<double> a, 00294 TransExpr<TA> b, 00295 double alpha, 00296 SliceMatrix<double> c, 00297 double beta = 1.0) 00298 { BASE_LapackMultAdd<double> (a, false, b.A(), true, alpha, c, beta); } 00299 00300 template <typename TA> 00301 inline void LapackMultAdd (TransExpr<TA> a, 00302 SliceMatrix<double> b, 00303 double alpha, 00304 SliceMatrix<double> c, 00305 double beta = 1.0) 00306 { BASE_LapackMultAdd<double> (a.A(), true, b, false, alpha, c, beta); } 00307 00308 template <typename TA, typename TB> 00309 inline void LapackMultAdd (TransExpr<TA> a, 00310 TransExpr<TB> b, 00311 double alpha, 00312 SliceMatrix<double> c, 00313 double beta = 1.0) 00314 { BASE_LapackMultAdd<double> (a.A(), true, b.A(), true, alpha, c, beta); } 00315 00316 00317 00318 inline void LapackMultAdd (SliceMatrix<Complex> a, 00319 SliceMatrix<Complex> b, 00320 Complex alpha, 00321 SliceMatrix<Complex> c, 00322 Complex beta = 1.0) 00323 { BASE_LapackMultAdd<Complex> (a, false, b, false, alpha, c, beta); } 00324 00325 template <typename TA> 00326 inline void LapackMultAdd (SliceMatrix<Complex> a, 00327 TransExpr<TA> b, 00328 Complex alpha, 00329 SliceMatrix<Complex> c, 00330 Complex beta = 1.0) 00331 { BASE_LapackMultAdd<Complex> (a, false, b.A(), true, alpha, c, beta); } 00332 00333 template <typename TA> 00334 inline void LapackMultAdd (TransExpr<TA> a, 00335 SliceMatrix<Complex> b, 00336 Complex alpha, 00337 SliceMatrix<Complex> c, 00338 Complex beta) 00339 { BASE_LapackMultAdd<Complex> (a.A(), true, b, false, alpha, c, beta); } 00340 00341 template <typename TA, typename TB> 00342 inline void LapackMultAdd (TransExpr<TA> a, 00343 TransExpr<TB> b, 00344 Complex alpha, 00345 SliceMatrix<Complex> c, 00346 Complex beta) 00347 { BASE_LapackMultAdd<Complex> (a.A(), true, b.A(), true, alpha, c, beta); } 00348 00349 00350 00351 00352 00353 00354 00355 /* 00356 00357 // old-style function names, new driver 00358 inline void LapackMultABt (SliceMatrix<double> a, 00359 SliceMatrix<double> b, 00360 SliceMatrix<double> c) 00361 { BASE_LapackMult (a, false, b, true, c); } 00362 00363 inline void LapackMultABt (SliceMatrix<Complex> a, 00364 SliceMatrix<Complex> b, 00365 SliceMatrix<Complex> c) 00366 { BASE_LapackMult (a, false, b, true, c); } 00367 00368 inline void LapackMultAB (SliceMatrix<double> a, 00369 SliceMatrix<double> b, 00370 SliceMatrix<double> c) 00371 { BASE_LapackMult (a, false, b, false, c); } 00372 00373 inline void LapackMultAB (SliceMatrix<Complex> a, 00374 SliceMatrix<Complex> b, 00375 SliceMatrix<Complex> c) 00376 { BASE_LapackMult (a, false, b, false, c); } 00377 00378 00379 inline void LapackMultAtB (SliceMatrix<double> a, 00380 SliceMatrix<double> b, 00381 SliceMatrix<double> c) 00382 { BASE_LapackMult (a, true, b, false, c); } 00383 00384 inline void LapackMultAtB (SliceMatrix<Complex> a, 00385 SliceMatrix<Complex> b, 00386 SliceMatrix<Complex> c) 00387 { BASE_LapackMult (a, true, b, false, c); } 00388 00389 00390 inline void LapackMultAtBt (SliceMatrix<double> a, 00391 SliceMatrix<double> b, 00392 SliceMatrix<double> c) 00393 { BASE_LapackMult (a, true, b, true, c); } 00394 inline void LapackMultAtBt (SliceMatrix<Complex> a, 00395 SliceMatrix<Complex> b, 00396 SliceMatrix<Complex> c) 00397 { BASE_LapackMult (a, true, b, true, c); } 00398 00399 */ 00400 00401 00402 00403 00404 00405 00406 // old functions for compatibilty 00407 inline void LapackMultABt (SliceMatrix<double> a, 00408 SliceMatrix<double> b, 00409 SliceMatrix<double> c) 00410 { 00411 char transa = 'T'; 00412 char transb = 'N'; 00413 integer m = c.Height(); 00414 integer n = c.Width(); 00415 integer k = a.Width(); 00416 double alpha = 1.0; 00417 double beta = 0; 00418 integer lda = a.Dist(); 00419 integer ldb = b.Dist(); 00420 integer ldc = c.Dist(); 00421 00422 dgemm_ (&transa, &transb, &n, &m, &k, &alpha, &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00423 } 00424 00425 00426 inline void LapackMultAB (SliceMatrix<double> a, 00427 SliceMatrix<double> b, 00428 SliceMatrix<double> c) 00429 { 00430 char transa = 'N'; 00431 char transb = 'N'; 00432 integer m = c.Height(); 00433 integer n = c.Width(); 00434 integer k = a.Width(); 00435 double alpha = 1.0; 00436 double beta = 0; 00437 integer lda = a.Dist(); 00438 integer ldb = b.Dist(); 00439 integer ldc = c.Dist(); 00440 dgemm_ (&transa, &transb, &n, &m, &k, &alpha, &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00441 } 00442 00443 inline void LapackMultAtB (ngbla::FlatMatrix<double> a, 00444 ngbla::FlatMatrix<double> b, 00445 ngbla::SliceMatrix<double> c) 00446 { 00447 char transa = 'N'; 00448 char transb = 'T'; 00449 integer m = c.Height(); // changed n,m 00450 integer n = c.Width(); 00451 integer k = a.Height(); 00452 double alpha = 1.0; 00453 double beta = 0; 00454 integer lda = a.Width(); 00455 integer ldb = b.Width(); 00456 integer ldc = c.Dist(); // c.Width(); 00457 00458 dgemm_ (&transa, &transb, &n, &m, &k, &alpha, &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00459 } 00460 00461 inline void LapackMultAtBt (ngbla::FlatMatrix<double> a, 00462 ngbla::FlatMatrix<double> b, 00463 ngbla::FlatMatrix<double> c) 00464 { 00465 char transa = 'T'; 00466 char transb = 'T'; 00467 integer m = c.Height(); 00468 integer n = c.Width(); 00469 integer k = a.Height(); 00470 double alpha = 1.0; 00471 double beta = 0; 00472 integer lda = a.Width(); 00473 integer ldb = b.Width(); 00474 integer ldc = c.Width(); 00475 00476 dgemm_ (&transa, &transb, &n, &m, &k, &alpha, &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00477 } 00478 00479 inline void LapackMultABt (ngbla::FlatMatrix<ngbla::Complex> a, 00480 ngbla::FlatMatrix<ngbla::Complex> b, 00481 ngbla::FlatMatrix<ngbla::Complex> c) 00482 { 00483 // c = a * Trans (b); 00484 00485 char transa = 'T'; 00486 char transb = 'N'; 00487 integer m = c.Height(); 00488 integer n = c.Width(); 00489 integer k = a.Width(); 00490 Complex alpha(1,0); // double alpha[2] = { 1.0, 0.0 }; 00491 Complex beta(0,0); // double beta[2] = { 0.0, 0.0 }; 00492 integer lda = a.Width(); 00493 integer ldb = b.Width(); 00494 integer ldc = c.Width(); 00495 00496 zgemm_ (&transa, &transb, &n, &m, &k, &alpha, 00497 &b(0,0), &ldb, 00498 &a(0,0), &lda, &beta, 00499 &c(0,0), &ldc); 00500 00501 } 00502 00503 00504 00505 inline void LapackMultAddABt (ngbla::FlatMatrix<double> a, 00506 ngbla::FlatMatrix<double> b, 00507 double fac, 00508 ngbla::FlatMatrix<double> c) 00509 { 00510 char transa = 'T'; 00511 char transb = 'N'; 00512 integer m = c.Height(); 00513 integer n = c.Width(); 00514 integer k = a.Width(); 00515 double alpha = fac; 00516 double beta = 1.0; 00517 integer lda = a.Width(); 00518 integer ldb = b.Width(); 00519 integer ldc = c.Width(); 00520 00521 dgemm_ (&transa, &transb, &n, &m, &k, &alpha, &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00522 } 00523 00524 00525 inline void LapackMultAddAB (ngbla::FlatMatrix<double> a, 00526 ngbla::FlatMatrix<double> b, 00527 double fac, 00528 ngbla::FlatMatrix<double> c) 00529 { 00530 char transa = 'N'; 00531 char transb = 'N'; 00532 integer m = c.Height(); 00533 integer n = c.Width(); 00534 integer k = a.Width(); 00535 double alpha = fac; 00536 double beta = 1.0; 00537 integer lda = a.Width(); 00538 integer ldb = b.Width(); 00539 integer ldc = c.Width(); 00540 00541 dgemm_ (&transa, &transb, &n, &m, &k, &alpha, &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00542 } 00543 00544 inline void LapackMultAddAtB (ngbla::FlatMatrix<double> a, 00545 ngbla::FlatMatrix<double> b, 00546 double fac, 00547 ngbla::FlatMatrix<double> c) 00548 { 00549 char transa = 'N'; 00550 char transb = 'T'; 00551 integer m = c.Height(); // changed n,m 00552 integer n = c.Width(); 00553 integer k = a.Height(); 00554 double alpha = fac; 00555 double beta = 1.0; 00556 integer lda = a.Width(); 00557 integer ldb = b.Width(); 00558 integer ldc = c.Width(); 00559 00560 dgemm_ (&transa, &transb, &n, &m, &k, &alpha, &b(0,0), &ldb, &a(0,0), &lda, &beta, &c(0,0), &ldc); 00561 } 00562 00563 00564 inline void LapackMultAddABt (ngbla::FlatMatrix<ngbla::Complex> a, 00565 ngbla::FlatMatrix<ngbla::Complex> b, 00566 double fac, 00567 ngbla::FlatMatrix<ngbla::Complex> c) 00568 { 00569 // c += fac * a * Trans (b); 00570 char transa = 'T'; 00571 char transb = 'N'; 00572 integer m = c.Height(); 00573 integer n = c.Width(); 00574 integer k = a.Width(); 00575 Complex alpha(fac, 0); 00576 Complex beta(1,0); 00577 integer lda = a.Width(); 00578 integer ldb = b.Width(); 00579 integer ldc = c.Width(); 00580 00581 zgemm_ (&transa, &transb, &n, &m, &k, &alpha, 00582 &b(0,0), &ldb, 00583 &a(0,0), &lda, &beta, 00584 &c(0,0), &ldc); 00585 } 00586 00587 00588 00589 inline void LapackMultAddAtB (ngbla::FlatMatrix<ngbla::Complex> a, 00590 ngbla::FlatMatrix<ngbla::Complex> b, 00591 double fac, 00592 ngbla::FlatMatrix<ngbla::Complex> c) 00593 { 00594 // c += fac * a * Trans (b); 00595 char transa = 'N'; 00596 char transb = 'T'; 00597 integer m = c.Height(); 00598 integer n = c.Width(); 00599 integer k = a.Height(); 00600 Complex alpha(fac, 0); // double alpha[2] = { fac, 0 }; 00601 Complex beta(1,0); // double beta[2] = { 1.0, 0 }; 00602 integer lda = a.Width(); 00603 integer ldb = b.Width(); 00604 integer ldc = c.Width(); 00605 00606 zgemm_ (&transa, &transb, &n, &m, &k, &alpha, 00607 &b(0,0), &ldb, 00608 &a(0,0), &lda, &beta, 00609 &c(0,0), &ldc); 00610 } 00611 00612 00613 00614 00615 inline void LapackMultAddAB (ngbla::FlatMatrix<ngbla::Complex> a, 00616 ngbla::FlatMatrix<ngbla::Complex> b, 00617 double fac, 00618 ngbla::FlatMatrix<ngbla::Complex> c) 00619 { 00620 char transa = 'N'; 00621 char transb = 'N'; 00622 integer m = c.Height(); 00623 integer n = c.Width(); 00624 integer k = a.Width(); 00625 Complex alpha(fac, 0); // double alpha[2] = { fac, 0 }; 00626 Complex beta(1,0); // double beta[2] = { 1.0, 0 }; 00627 integer lda = a.Width(); 00628 integer ldb = b.Width(); 00629 integer ldc = c.Width(); 00630 00631 zgemm_ (&transa, &transb, &n, &m, &k, &alpha, 00632 &b(0,0), &ldb, 00633 &a(0,0), &lda, &beta, 00634 &c(0,0), &ldc); 00635 } 00636 00637 00638 00639 00640 /* 00641 extern "C" 00642 void dgetri_ (int & n, double & a, int & lda, 00643 int & ipiv, double & work, int & lwork, int & info); 00644 00645 extern "C" 00646 void dgetrf_ (int & n, int & m, double & a, int & lda, 00647 int & ipiv, int & info); 00648 00649 extern "C" 00650 void dgetrs_ (char & trans, int & n, int & nrhs, 00651 double & a, int & lda, int & ipiv, 00652 double & b, int & ldb, int & info); 00653 00654 00655 00656 extern "C" 00657 void dpotrf_ (char & uplo, int & n, double & a, int & lda, 00658 int & info); 00659 00660 extern "C" 00661 void dpotrs_ (char & uplo, int & n, int & nrhs, 00662 double & a, int & lda, 00663 double & b, int & ldb, int & info); 00664 */ 00665 00666 00667 00668 inline void LapackInverse (ngbla::FlatMatrix<double> a) 00669 { 00670 integer m = a.Height(); 00671 if (m == 0) return; 00672 integer n = a.Width(); 00673 integer lda = a.Width(); 00674 00675 Array<integer> ipiv(n); 00676 integer info; 00677 00678 dgetrf_ (&n, &m, &a(0,0), &lda, &ipiv[0], &info); 00679 00680 double hwork; 00681 integer lwork = -1; 00682 dgetri_ (&n, &a(0,0), &lda, &ipiv[0], &hwork, &lwork, &info); 00683 lwork = integer(hwork); 00684 00685 Array<double> work(lwork); 00686 dgetri_ (&n, &a(0,0), &lda, &ipiv[0], &work[0], &lwork, &info); 00687 } 00688 00689 00690 00691 inline void LapackInverseSPD (ngbla::FlatMatrix<double> a) 00692 { 00693 integer n = a.Width(); 00694 if (n == 0) return; 00695 integer lda = n; 00696 00697 integer info; 00698 char uplo = 'L'; 00699 00700 dpotrf_ (&uplo, &n, &a(0,0), &lda, &info); 00701 dpotri_ (&uplo, &n, &a(0,0), &lda, &info); 00702 for (int i = 0; i < n; i++) 00703 for (int j = 0; j < i; j++) 00704 a(i,j) = a(j,i); 00705 } 00706 00707 00708 00709 00710 00711 00712 00713 00714 00715 /* 00716 Compoutes B <--- B A^{-1} (trans = 'N') 00717 Compoutes B <--- B A^{-T} (trans = 'T') 00718 Compoutes B <--- B A^{-H} (trans = 'H') 00719 00720 // trans = 'N': solve A x = b^T 00721 // trans = 'T': solve A^T x = b^T 00722 // trans = 'C': solve A^H x = b^T 00723 */ 00724 inline void LapackAInvBt (ngbla::FlatMatrix<double> a, ngbla::FlatMatrix<double> b, char trans = 'N') 00725 { 00726 integer m = a.Height(); 00727 integer n = a.Width(); 00728 integer lda = a.Width(); 00729 integer ldb = b.Width(); 00730 integer nrhs = b.Height(); 00731 00732 ArrayMem<integer,100> ipiv(n); 00733 integer info; 00734 // char uplo = 'L'; 00735 00736 dgetrf_ (&n, &m, &a(0,0), &lda, &ipiv[0], &info); 00737 dgetrs_ (&trans, &n, &nrhs, &a(0,0), &lda, &ipiv[0], &b(0,0), &ldb, &info); 00738 00739 /* 00740 // symmetric, non-spd 00741 int lwork = -1; 00742 double hwork[1] = { 1.0 }; 00743 dsytrf_ (&uplo, &n, &a(0,0), &lda, &ipiv[0], &hwork[0], &lwork, &info); 00744 lwork = int(hwork[0]); 00745 ArrayMem<double, 1000> work(lwork); 00746 dsytrf_ (&uplo, &n, &a(0,0), &lda, &ipiv[0], &work[0], &lwork, &info); 00747 dsytrs_ (&uplo, &n, &nrhs, &a(0,0), &lda, &ipiv[0], &b(0,0), &ldb, &info); 00748 */ 00749 00750 /* 00751 // spd 00752 dpotrf_ (&uplo, &n, &a(0,0), &lda, &info); 00753 dpotrs_ (&uplo, &n, &nrhs, &a(0,0), &lda, &b(0,0), &ldb, &info); 00754 */ 00755 } 00756 00757 00758 00759 00760 /* 00761 extern "C" 00762 void zgetri_ (int & n, double & a, int & lda, 00763 int & ipiv, double & work, int & lwork, int & info); 00764 00765 extern "C" 00766 void zgetrf_ (int & n, int & m, double & a, int & lda, 00767 int & ipiv, int & info); 00768 00769 extern "C" 00770 void zgetrs_ (char & trans, int & m, int & nrhs, double & a, int & lda, 00771 int & ipiv, 00772 double & b, int & ldb, 00773 int & info); 00774 */ 00775 00776 00777 inline void LapackInverse (ngbla::FlatMatrix<ngbla::Complex> a) 00778 { 00779 integer m = a.Height(); 00780 if (m == 0) return; 00781 00782 integer n = a.Width(); 00783 integer lda = a.Width(); 00784 integer * ipiv = new integer[n]; 00785 integer lwork = 100*n; 00786 Complex * work = new Complex[lwork]; 00787 integer info; 00788 00789 // std::cout << "a = " << std::endl << a << std::endl; 00790 zgetrf_ (&n, &m, &a(0,0), &lda, ipiv, &info); 00791 // std::cout << "factors = " << std::endl << a << std::endl; 00792 00793 if (info != 0) 00794 { 00795 std::cout << "ZGETRF::info = " << info << std::endl; 00796 // *testout << "ZGETRF::info = " << info << std::endl; 00797 // *testout << "a = " << endl << a << endl; 00798 } 00799 zgetri_ (&n, &a(0,0), &lda, ipiv, work, &lwork, &info); 00800 if (info != 0) 00801 std::cout << "ZGETRI::info = " << info << std::endl; 00802 00803 delete [] work; 00804 delete [] ipiv; 00805 00806 } 00807 00808 00809 /* 00810 trans = 'N': solve A x = b^T 00811 trans = 'T': solve A^T x = b^T 00812 trans = 'C': solve A^H x = b^T 00813 */ 00814 00815 inline void LapackAInvBt (ngbla::FlatMatrix<ngbla::Complex> a, ngbla::FlatMatrix<ngbla::Complex> b, char trans = 'N') 00816 { 00817 integer m = a.Height(); 00818 integer n = a.Width(); 00819 integer lda = a.Width(); 00820 integer ldb = b.Width(); 00821 integer nrhs = b.Height(); 00822 integer * ipiv = new integer[n]; 00823 integer lwork = 100*n; 00824 double * work = new double[2*lwork]; 00825 integer info; 00826 00827 00828 zgetrf_ (&n, &m,&a(0,0), &lda, ipiv, &info); 00829 zgetrs_ (&trans, &n, &nrhs, &a(0,0), &lda, ipiv, &b(0,0), &ldb, &info); 00830 00831 delete [] work; 00832 delete [] ipiv; 00833 // std::cerr << "complex LapackAInvBt not implemented" << std::endl; 00834 } 00835 00836 00837 /* 00838 extern "C" 00839 void dsyev_(char & jobz, char & uplo, int & n , double & A , int & lda, double & w, 00840 double & work, int & lwork, int & info); 00841 00842 extern "C" 00843 void zgeev_( char *jobvl, char *jobvr, int *n, std::complex<double> *A, int * lda, std::complex<double>* lami, 00844 std::complex<double> * vl, int * nvl, std::complex<double> * vr, int * nvr, 00845 std::complex<double> * work, int * lwork, double * rwork, int * info); 00846 */ 00847 00848 // extern "C" 00849 // void dgeev_( char *jobvl, char *jobvr, int *n, double *A, int * lda, double* lami_re, double *lami_im, 00850 // double * vl, int * nvl, double * vr, int * nvr, 00851 // double * work, int * lwork, /* double * rwork, */ int * info); 00852 00853 00854 00855 inline void LapackEigenValuesSymmetric (ngbla::FlatMatrix<double> a, 00856 ngbla::FlatVector<double> lami, 00857 ngbla::FlatMatrix<double> evecs = ngbla::FlatMatrix<double>(0,0)){ 00858 char jobz, uplo = 'U'; 00859 integer n = a.Height(); 00860 integer lwork=(n+2)*n+1; 00861 00862 double* work = new double[lwork]; 00863 integer info; 00864 00865 double * matA; 00866 00867 if ( evecs.Height() ) 00868 { 00869 // eigenvectors are calculated 00870 evecs = a; 00871 jobz = 'V'; 00872 matA = &evecs(0,0); 00873 } 00874 else 00875 { 00876 // only eigenvalues are calculated, matrix a is destroyed!! 00877 jobz = 'N'; 00878 matA = &a(0,0); 00879 } 00880 dsyev_(&jobz, &uplo , &n , matA, &n, &lami(0), work, &lwork, &info); 00881 00882 if (info) 00883 std::cerr << "LapackEigenValuesSymmetric, info = " << info << std::endl; 00884 00885 delete [] work; 00886 } 00887 00888 00889 00890 inline void LapackEigenValues (ngbla::FlatMatrix<double> a, 00891 ngbla::FlatVector<ngbla::Complex> lami, 00892 ngbla::FlatMatrix<double> eveci ) 00893 { 00894 char jobvr = 'V' , jobvl= 'N'; 00895 00896 integer n = a.Height(); 00897 integer nvl = 1; 00898 integer nvr = eveci.Width() ; 00899 00900 double * vl = 0; 00901 double * vr;// = new std::complex<double> [nvr*n]; 00902 double * lami_re = new double[n], * lami_im = new double[n]; 00903 00904 integer lwork = 8*n; 00905 double * work = new double[lwork]; 00906 double *rwork = new double[8*n]; 00907 integer info = 0; 00908 00909 if ( eveci.Width() ) 00910 { 00911 vr = &eveci(0,0); 00912 } 00913 else 00914 { 00915 nvr = n; 00916 vr = new double [nvr*n]; 00917 } 00918 00919 dgeev_(&jobvl, &jobvr, &n, &a(0,0), &n, lami_re, lami_im, vl, &nvl, vr, &nvr, work, &lwork, /* rwork, */ &info); 00920 00921 if(info != 0) 00922 { 00923 std::cout << "**** Error in zggev_, info = " << info << " *****" << std::endl; 00924 return; 00925 } 00926 00927 for ( int i = 0; i < lami.Size(); i++ ) 00928 lami(i) = ngbla::Complex (lami_re[i], lami_im[i]); 00929 00930 delete[] work; 00931 delete[] rwork; 00932 if ( !eveci.Width() ) 00933 delete[] vr; 00934 delete [] lami_re; 00935 delete [] lami_im; 00936 } 00937 00938 inline void LapackEigenValues (ngbla::FlatMatrix<ngbla::Complex> a, 00939 ngbla::FlatVector<ngbla::Complex> lami, 00940 ngbla::FlatMatrix<ngbla::Complex> eveci ) 00941 { 00942 char jobvr = 'V' , jobvl= 'N'; 00943 00944 integer n = a.Height(); 00945 integer nvl = 1; 00946 integer nvr = eveci.Width() ; 00947 00948 std::complex<double> * vl = 0; 00949 std::complex<double> * vr;// = new std::complex<double> [nvr*n]; 00950 00951 integer lwork = 8*n; 00952 std::complex<double> * work = new std::complex<double>[lwork]; 00953 double *rwork = new double[8*n]; 00954 integer info = 0; 00955 00956 if ( eveci.Width() ) 00957 { 00958 vr = &eveci(0,0); 00959 } 00960 else 00961 { 00962 nvr = n; 00963 vr = new std::complex<double> [nvr*n]; 00964 } 00965 00966 zgeev_(&jobvl, &jobvr, &n, (std::complex<double>*)(void*)&a(0,0), &n, (std::complex<double>*)(void*)&lami(0), vl, &nvl, vr, &nvr, work, &lwork, rwork, &info); 00967 // alpha, beta, &vl, &nvl, vr, &nvr, 00968 // work , &lwork, rwork, &info); 00969 00970 if(info != 0) 00971 { 00972 std::cout << "**** Error in zggev_, info = " << info << " *****" << std::endl; 00973 return; 00974 } 00975 00976 delete[] work; 00977 delete[] rwork; 00978 if ( !eveci.Width() ) 00979 delete[] vr; 00980 } 00981 00982 00983 00984 00985 00986 inline void LapackEigenValuesSymmetric (ngbla::FlatMatrix<ngbla::Complex> a, 00987 ngbla::FlatVector<ngbla::Complex> lami, 00988 ngbla::FlatMatrix<ngbla::Complex> eveci = ngbla::FlatMatrix<ngbla::Complex> (0,0) ) 00989 { 00990 // std::cerr << "complex evp not implemented" << std::endl; 00991 LapackEigenValues ( a, lami, eveci ); 00992 } 00993 00994 00995 00996 // Solve complex generalized eigenvalue problem (QZ) 00997 /* 00998 extern "C" 00999 void zggev_(char *jobvl,char* jobvr,int* N,std::complex<double>* A, int* lda,std::complex<double>* B, int* ldb, std::complex<double>* alpha, std::complex<double>* beta, std::complex<double>* vl, int* ldvl, std::complex<double>* vr, int* ldvr, std::complex<double>* work, int* lwork, double* rwork, int* info); 01000 */ 01001 01002 01003 // Attention A,B are overwritten !!! 01004 inline void LapackEigenValues (ngbla::FlatMatrix<ngbla::Complex> a, 01005 ngbla::FlatMatrix<ngbla::Complex> b, 01006 ngbla::FlatVector<ngbla::Complex> lami) 01007 01008 { 01009 integer n = a.Height(); 01010 // integer evecs_bool = 0; 01011 // std::complex<double> * evecs, * dummy; 01012 01013 char jobvr = 'N', jobvl= 'N'; 01014 // bool balancing = 0; 01015 01016 std::complex<double> * alpha= new std::complex<double>[n]; 01017 std::complex<double> * beta = new std::complex<double>[n]; 01018 std::complex<double> vl=0.; 01019 01020 integer nvl = 1; 01021 std::complex<double> * vr ; 01022 01023 std::complex<double> * work = new std::complex<double>[8*n]; 01024 integer lwork = 8*n; 01025 double *rwork = new double[8*n]; 01026 01027 integer nvr = n ; 01028 01029 //std::complex<double> * A1,*B1; 01030 integer i; 01031 01032 // char job=balance_type; // Permute and Scale in Balancing 01033 // integer ihi,ilo; 01034 double * lscale, *rscale; 01035 lscale = new double[n]; 01036 rscale = new double[n]; 01037 double * work2; 01038 work2 = new double[6*n]; 01039 01040 // char side = 'R'; 01041 01042 integer info = 0; 01043 01044 // integer ii; 01045 // ii=0; 01046 01047 // if(balancing) zggbal_(&job,&n, A, &n , B, &n, &ilo, &ihi, lscale, rscale, work2, &info) ; 01048 01049 // if(info == 0 ) 01050 if (1) 01051 { 01052 zggev_(&jobvl, &jobvr, &n, &a(0,0), &n, &b(0,0), &n, alpha, beta, &vl, &nvl, vr, &nvr, 01053 work , &lwork, rwork, &info); 01054 01055 if(info==0) 01056 { 01057 /* 01058 if(jobvr == 'V' && balancing) 01059 { 01060 zggbak_(&job, &side, &n, &ilo, &ihi, lscale, rscale, &n, vr, &n,&info) ; 01061 01062 if(info!=0) 01063 { 01064 std::cout << "***** Error in zggbak_ **** " << endl; 01065 return; 01066 } 01067 } 01068 */ 01069 } 01070 else 01071 { 01072 std::cout << "**** Error in zggev_, info = " << info << " *****" << std::endl; 01073 // return; 01074 } 01075 } 01076 else 01077 { 01078 std::cout << "**** Error in zggbal_ **** " << std::endl; 01079 return; 01080 } 01081 01082 delete [] work; 01083 delete [] rwork; 01084 01085 delete [] lscale; 01086 delete [] rscale; 01087 delete [] work2; 01088 01089 for(i=0;i<n;i++) 01090 { 01091 if(abs(beta[i]) >= 1.e-30) 01092 lami[i]=std::complex<double>(alpha[i]/beta[i]); 01093 else 01094 { 01095 lami[i] = std::complex<double>(100.,100.); 01096 } 01097 } 01098 01099 /* 01100 std:: complex<double> resid[n]; 01101 double error; 01102 for(int k=0; k<n;k++) 01103 { 01104 for(i=0; i<n;i++) 01105 { 01106 resid[i] = 0.; 01107 for(int j=0;j<n;j++) 01108 { 01109 resid[i] += A[j*n + i]*evecs[k*n+j]; 01110 resid[i] -=B[j*n+i]*evecs[k*n+j]*lami[k]; 01111 } 01112 01113 error = abs(resid[i]) * abs(resid[i]) ; 01114 } 01115 error = sqrt(error); 01116 cout << " lami (i) " << lami[k] << "\t" << alpha[k] << "\t" << beta[k] << endl; 01117 cout << " error lapack " << k << " \t " << error << endl; 01118 } 01119 */ 01120 delete [] alpha; 01121 delete [] beta; 01122 01123 } 01124 01125 01126 01127 01128 01129 01130 01131 01132 01133 01134 /* 01135 extern "C" 01136 void dsygv_(int & itype, char & jobzm, char & uplo , int & n1, double & A , int & n, 01137 double & B, int & nb, double & lami, double & work, int & lwork, int & info); 01138 */ 01139 01140 inline void LapackEigenValuesSymmetric (ngbla::FlatMatrix<double> a, 01141 ngbla::FlatMatrix<double> b, 01142 ngbla::FlatVector<double> lami, 01143 ngbla::FlatMatrix<double> evecs = ngbla::FlatMatrix<double>(0,0)) 01144 { 01145 char jobz = 'N' , uplo = 'U'; 01146 integer n = a.Height(); 01147 01148 integer lwork=(n+2)*n+1; 01149 double* work = new double[lwork]; 01150 01151 integer info; 01152 integer itype =1; 01153 01154 if ( evecs.Height() ) 01155 jobz = 'V'; 01156 else 01157 jobz = 'N'; 01158 01159 dsygv_(&itype, &jobz, &uplo , &n , &a(0,0), &n, &b(0,0), &n, 01160 &lami(0), work, &lwork, &info); 01161 01162 if ( evecs.Height() ) 01163 evecs = a; 01164 01165 if (info) 01166 std::cerr << "LapackEigenValuesSymmetric, info = " << info << std::endl; 01167 01168 delete [] work; 01169 } 01170 01171 01172 #else 01173 01174 typedef int integer; 01175 01176 inline void LapackMultAtx (ngbla::FlatMatrix<double> a, 01177 ngbla::FlatVector<double> x, 01178 ngbla::FlatVector<double> y) 01179 { y = Trans (a) * x; } 01180 01181 01182 inline void LapackAddxyt (ngbla::FlatMatrix<double> a, 01183 ngbla::FlatVector<double> x, 01184 ngbla::FlatVector<double> y) 01185 { y += Trans (a) * x; } 01186 01187 01188 template <typename TA, typename TB> 01189 inline void LapackMult (const TA & a, const TB & b, 01190 ngbla::SliceMatrix<double> c) 01191 { c = a * b; } 01192 01193 template <typename TA, typename TB> 01194 inline void LapackMult (const TA & a, const TB & b, 01195 ngbla::SliceMatrix<Complex> c) 01196 { c = a * b; } 01197 01198 template <typename TA, typename TB> 01199 inline void LapackMultAdd (const TA & a, 01200 const TB & b, 01201 double alpha, 01202 SliceMatrix<double> c, 01203 double beta) 01204 { c *= beta; c += alpha * a * b; } 01205 01206 template <typename TA, typename TB> 01207 inline void LapackMultAdd (const TA & a, 01208 const TB & b, 01209 Complex alpha, 01210 SliceMatrix<Complex> c, 01211 Complex beta) 01212 { c *= beta; c += alpha * a * b; } 01213 01214 01215 01216 inline void LapackMultABt (ngbla::FlatMatrix<double> a, 01217 ngbla::FlatMatrix<double> b, 01218 ngbla::FlatMatrix<double> c) 01219 { c = a * Trans (b); } 01220 01221 inline void LapackMultAtB (ngbla::FlatMatrix<double> a, 01222 ngbla::FlatMatrix<double> b, 01223 ngbla::FlatMatrix<double> c) 01224 { c = Trans(a) * b; } 01225 01226 01227 inline void LapackMultAB (ngbla::FlatMatrix<double> a, 01228 ngbla::FlatMatrix<double> b, 01229 ngbla::FlatMatrix<double> c) 01230 { c = a * b; } 01231 01232 01233 inline void LapackMultABt (ngbla::FlatMatrix<ngbla::Complex> a, 01234 ngbla::FlatMatrix<ngbla::Complex> b, 01235 ngbla::FlatMatrix<ngbla::Complex> c) 01236 { c = a * Trans (b); } 01237 01238 inline void LapackMultAtB (ngbla::FlatMatrix<Complex> a, 01239 ngbla::FlatMatrix<Complex> b, 01240 ngbla::FlatMatrix<Complex> c) 01241 { c = Trans(a) * b; } 01242 01243 01244 01245 inline void LapackMultAddAB (ngbla::FlatMatrix<double> a, 01246 ngbla::FlatMatrix<double> b, 01247 double fac, 01248 ngbla::FlatMatrix<double> c) 01249 { c += fac * a * b; } 01250 01251 inline void LapackMultAddABt (ngbla::FlatMatrix<double> a, 01252 ngbla::FlatMatrix<double> b, 01253 double fac, 01254 ngbla::FlatMatrix<double> c) 01255 { c += fac * a * Trans (b); } 01256 01257 inline void LapackMultAddAtB (ngbla::FlatMatrix<double> a, 01258 ngbla::FlatMatrix<double> b, 01259 double fac, 01260 ngbla::FlatMatrix<double> c) 01261 { c += fac * Trans(a) * b; } 01262 01263 01264 01265 inline void LapackMultAddAB (ngbla::FlatMatrix<ngbla::Complex> a, 01266 ngbla::FlatMatrix<ngbla::Complex> b, 01267 double fac, 01268 ngbla::FlatMatrix<ngbla::Complex> c) 01269 01270 { c += fac * a * b; } 01271 01272 inline void LapackMultAddABt (ngbla::FlatMatrix<ngbla::Complex> a, 01273 ngbla::FlatMatrix<ngbla::Complex> b, 01274 double fac, 01275 ngbla::FlatMatrix<ngbla::Complex> c) 01276 01277 { c += fac * a * Trans (b); } 01278 01279 01280 01281 01282 inline void LapackInverse (ngbla::FlatMatrix<double> a) 01283 { 01284 CalcInverse (a); 01285 /* 01286 ngbla::Matrix<> hm(a.Height()); 01287 CalcInverse (a, hm); 01288 a = hm; 01289 */ 01290 } 01291 01292 inline void LapackInverse (ngbla::FlatMatrix<ngbla::Complex> a) 01293 { 01294 CalcInverse (a); 01295 /* 01296 // std::cerr << "sorry, Inverse not available without LAPACK" << std::endl; 01297 ngbla::Matrix<Complex> hm(a.Height()); 01298 CalcInverse (a, hm); 01299 a = hm; 01300 */ 01301 } 01302 01303 inline void LapackAInvBt (ngbla::FlatMatrix<double> a, ngbla::FlatMatrix<double> b, char trans = 'N') 01304 { 01305 LapackInverse (a); 01306 ngbla::Matrix<> hb (b.Height(), b.Width()); 01307 if (trans == 'T') 01308 hb = b * Trans(a); 01309 else 01310 hb = b * a; 01311 b = hb; 01312 } 01313 01314 inline void LapackAInvBt (ngbla::FlatMatrix<Complex> a, ngbla::FlatMatrix<Complex> b, char trans = 'N') 01315 { 01316 LapackInverse (a); 01317 ngbla::Matrix<Complex> hb (b.Height(), b.Width()); 01318 if (trans == 'T') 01319 hb = b * Trans(a); 01320 else 01321 hb = b * a; 01322 b = hb; 01323 } 01324 01325 01326 inline void LapackEigenValuesSymmetric (ngbla::FlatMatrix<double> a, 01327 ngbla::FlatVector<double> lami) 01328 { 01329 std::cerr << "sorry, EVP not available without LAPACK" << std::endl; 01330 } 01331 01332 inline void LapackEigenValuesSymmetric (ngbla::FlatMatrix<double> a, 01333 ngbla::FlatMatrix<double> b, 01334 ngbla::FlatVector<double> lami) 01335 { 01336 std::cerr << "sorry, EVP not available without LAPACK" << std::endl; 01337 } 01338 01339 01340 inline void LapackEigenValuesSymmetric (ngbla::FlatMatrix<ngbla::Complex> a, 01341 ngbla::FlatVector<ngbla::Complex> lami) 01342 { 01343 std::cerr << "sorry, EVP not available without LAPACK" << std::endl; 01344 } 01345 01346 01347 #endif 01348 01349 01350 01351 01352 // several LAPACK eigenvalue solvers 01353 01354 #ifdef LAPACK 01355 01356 void LaEigNSSolve(int n, double * A, double * B, std::complex<double> * lami, int evecs_bool, double *evecs_re, double *evecs_im, char balance_type); 01357 void LaEigNSSolve(int n, std::complex<double> * A, std::complex<double> * B, std::complex<double> * lami, int evecs_bool, std::complex<double> *evecs, std::complex<double> *dummy, char balance_type); 01358 01359 void LapackSSEP(int n, double* A, double* lami, double* evecs); 01360 01361 void LapackHessenbergEP (int n, std::complex<double> * H, std::complex<double> * lami, std::complex<double> * evecs); 01362 void LapackGHEP(int n, double* A, double* B, double* lami) ; 01363 int LapackGHEPEPairs(int n, double* A, double* B, double* lami); 01364 int LapackGHEPEPairs(int n, std::complex<double>* A, std::complex<double>* B, double* lami); 01365 // A,B overwritten in A eigenvectors z^H B z = 1 01366 01367 //void LaEigNSSolve(const LaGenMatDouble &A, LaVectorDouble &eigvals); 01368 //void LaEigNSSolveIP(LaGenMatDouble &A, LaVectorDouble &eigvals); 01369 01370 void LaEigNSSolveTest(); 01371 void LaLinearSolveComplex(int n, std::complex<double> * A, std::complex<double> * F); 01372 void LaLinearSolve(int n, double * A, double * F); 01373 void LaLinearSolveRHS(int n, double * A, double * F); 01374 01375 01376 void LaEigNSSolveX(int n, std::complex<double> * A, std::complex<double> * B, std::complex<double> * lami, int evecs_bool, std::complex<double> * evecs, std::complex<double> * dummy, char balance_type); 01377 void LaEigNSSolveX(int n, double * A, double * B, std::complex<double> * lami, int evecs_bool, double * evecs, double * dummy, char balance_type); 01378 01379 #else 01380 01381 inline void LapackHessenbergEP (int n, std::complex<double> * H, std::complex<double> * lami, std::complex<double> * evecs) 01382 { 01383 cerr << "Sorry, HessebergEP not available without Lapack" << endl; 01384 } 01385 01386 01387 #endif 01388 01389 } 01390 01391 01392 #endif