NGSolve  4.9
basiclinalg/ng_lapack.hpp
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