TBCI Numerical high perf. C++ Library: lapack

00001 
00005 /* $Id: lapack_stdcplx.cpp,v 1.1.2.13 2008/07/19 09:38:44 garloff Exp $ */
00006 
00007 //#define NO_NS
00008 
00010 #define LAPACK_INLINE inline
00011 
00012 #include "std_cplx.h"
00013 #include "f_matrix.h"
00014 #include "f_bandmatrix.h"
00015 #include "lapack/lapack_lib.h"
00016 
00017 
00018 NAMESPACE_TBCI
00019 
00020 //********************************************************************************
00021 // Solution of linear eqution systems, partial pivoting
00022 //********************************************************************************
00023 TVector<double> lu_solve_expert(const F_BandMatrix<double>& A, const Vector<double>& B, int equi);
00024 TVector<double> lu_solve(const F_BandMatrix<double>& A, const Vector<double>& B);
00025 TVector<double> lu_solve(F_BandMatrix<double>& A, const Vector<double>& B, integer kl, integer ku);
00026 TVector<double> lu_solve(const F_Matrix<double>& A, const Vector<double>& B, int overwriteA);
00027 F_TMatrix<double> lu_solve(const F_Matrix<double>& A, const F_Matrix<double>& B, int overwriteA);
00028 
00029 F_TMatrix<double> inv(const F_Matrix<double>& A, int overwriteA);
00030 
00031 
00032 F_TMatrix<CPLX__ complex<double> > inv(const F_Matrix<CPLX__ complex<double> >& A, int overwriteA)
00033 {
00034         integer info = 0;
00035         integer n=A.columns(); // Dimension des Problems
00036         BVector<integer> ipiv(n);
00037         F_Matrix<CPLX__ complex<double> > temp;
00038         if (!overwriteA) 
00039         {
00040                 temp.resize(n,n); 
00041                 temp=A;
00042         }
00043         // Rechte Seite ist die Einheitsmatrix!
00044         F_TMatrix<CPLX__ complex<double> > invA(0.0, n, n); 
00045         for (int i=0; i<n; i++) 
00046                 invA.setval(CPLX__ complex<double>(1.0), i, i);
00047 
00048         if (overwriteA)
00049                 zgesv_(&n, &n, (doublecomplex*) A.get_fortran_matrix(), &n, ipiv.get_fortran_vector(), 
00050                 (doublecomplex*) invA.get_fortran_matrix(), &n, &info);
00051         else
00052                 zgesv_(&n, &n, (doublecomplex*) temp.get_fortran_matrix(), &n, ipiv.get_fortran_vector(), 
00053                 (doublecomplex*) invA.get_fortran_matrix(), &n, &info);
00054         if (info) 
00055                 STD__ cerr << "Error inverting Matrix A!" << STD__ endl;
00056         return invA;
00057 }
00058 
00059 F_TMatrix<CPLX__ complex<double> > lu_solve(const F_Matrix<CPLX__ complex<double> >& A, const F_Matrix<CPLX__ complex<double> >& B, int overwriteA)
00060 {
00061         integer info = 0;
00062         integer n=A.columns(); // Dimension des Problems
00063         BVector<integer> ipiv(n);
00064         F_Matrix<CPLX__ complex<double> > temp;
00065         if (!overwriteA) 
00066         {
00067                 temp.resize(n,n); 
00068                 temp=A;
00069         }
00070         F_TMatrix<CPLX__ complex<double> > invA(B);
00071 
00072         if (overwriteA)
00073                 zgesv_(&n, &n, (doublecomplex*) A.get_fortran_matrix(), &n, ipiv.get_fortran_vector(), 
00074                 (doublecomplex*) invA.get_fortran_matrix(), &n, &info);
00075         else
00076                 zgesv_(&n, &n, (doublecomplex*) temp.get_fortran_matrix(), &n, ipiv.get_fortran_vector(), 
00077                 (doublecomplex*) invA.get_fortran_matrix(), &n, &info);
00078         if (info) 
00079                 STD__ cerr << "Error inverting Matrix A!" << STD__ endl;
00080         return invA;
00081 }
00082 
00083 
00084 //********************************************************************************
00085 // eigenvalues (and eigenvectors) \a A*EV=EW*EV of symmetric double Matrix \a A
00086 //********************************************************************************
00087 int eig(const F_Matrix<double>& A, Vector<double>& EW);
00088 
00089 int eig(const F_Matrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV);
00090 
00091 //symmtric Band Matrix
00092 int eig(const F_BandMatrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV);
00093 
00094 int eig(const F_BandMatrix<double>& A, Vector<double>& EW, F_Matrix<double>& EV, 
00095         double ew_min, double ew_max);
00096 
00097 //********************************************************************************
00099 //********************************************************************************
00100 int eig(const F_Matrix<CPLX__ complex<double> >& A, Vector<CPLX__ complex<double> >& EW)
00101 {
00102         char jobvl='N';         //linke Eigenwerte und Eigenvektoren
00103         char jobvr='N';         //rechte Eigenwerte und Eigenvektoren
00104         integer N=A.columns();  // Dimension des Problems
00105         doublecomplex* vl=0;      // Keine linken Eigenvektoren
00106         doublecomplex* vr=0;      // Keine linken Eigenvektoren
00107         integer lwork=2*N;
00108         doublecomplex* work =NEW(doublecomplex,lwork);
00109         doublereal*    rwork=NEW(doublereal,lwork);
00110         integer info = 0;
00111         zgeev_(&jobvl, &jobvr, &N,
00112                (doublecomplex*) A.get_fortran_matrix(), &N,  
00113                (doublecomplex*) EW.get_fortran_vector(), vl, 
00114                &N, vr, &N, work, &lwork, rwork, &info);
00115         TBCIDELETE(doublecomplex, work, lwork);
00116         TBCIDELETE(doublereal, rwork, lwork);
00117         return (int) info;
00118 }
00119 
00120 
00121 int eig(const F_Matrix<CPLX__ complex<double> >& A, Vector<CPLX__ complex<double> >& EW, F_Matrix<CPLX__ complex<double> >& EV)
00122 {
00123         char jobvl='N';         //linke Eigenwerte und Eigenvektoren
00124         char jobvr='V';         //rechte Eigenwerte und Eigenvektoren
00125         integer N=A.columns();  // Dimension des Problems
00126         doublecomplex* vl=0;      // Keine linken Eigenvektoren
00127         integer lwork=2*N;
00128         doublecomplex* work =NEW(doublecomplex,lwork);
00129         doublereal*    rwork=NEW(doublereal,lwork);
00130         integer info = 0;
00131         zgeev_(&jobvl, &jobvr, &N, (doublecomplex*) A.get_fortran_matrix(), &N,  
00132                (doublecomplex*) EW.get_fortran_vector(), vl, 
00133                &N, (doublecomplex*) EV.get_fortran_matrix(), &N, work, &lwork, rwork, &info);
00134         TBCIDELETE(doublecomplex, work, lwork);
00135         TBCIDELETE(doublereal, rwork, lwork);
00136         return (int) info;
00137 }
00138 
00139 
00140 
00141 
00142 //********************************************************************************
00144 //********************************************************************************
00145 int eig(const F_BandMatrix<CPLX__ complex<double> >& A, 
00146         const F_BandMatrix<CPLX__ complex<double> >& B,
00147         double EW_min, double EW_max, Vector<double>& Eigenwerte, 
00148         F_Matrix<CPLX__ complex<double> >& Eigenvectoren)
00149 {
00150         // Benutzt wird nur die obere Haelfte der Bandmatrix A!!!!
00151         long uplo=1;
00152         // Ordnung des EWP bestimmen
00153         long n=A.columns();
00154         // Anzahl der Oberen Diagonalen der Matrizen A und B bestimmen
00155         long ldab=A.numSuper()+1;
00156         long ldbb=B.numSuper()+1;
00157         // Allg. Variablen
00158         long Anzahl_Eigenwerte;
00159         long BerechneEV=1; // Eigenvektoren sollen berechnet werden
00160 
00161         // Laufvariablen
00162         int i,j, x,y;
00163 
00164 
00165         // *****************************************************************************
00166         // Arrays im Fortran-Stil deklarieren,
00167         // diese Arrays werden mit den Werten aus den Matrizen besetzt
00168         //cout << "ab.bb..." << flush;
00169         // *****************************************************************************
00170         doublecomplex *ab = new doublecomplex[ldab*n];
00171         doublecomplex *bb = new doublecomplex[ldbb*n];
00172 
00173         // Arrays kopieren
00174         for (i=0; i<n; i++)
00175         {
00176                 // Array A nach ab kopieren
00177                 for (j=i; j<i+ldab && j<n; j++)
00178                 {
00179                         // x und y sind die Koordinaten in Lapack-Bandstruktur
00180                         x=ldab+i-j-1;
00181                         y=j;
00182                         ab[x+ldab*y].r = CPLX__ real(A(i,j));
00183                         ab[x+ldab*y].i = CPLX__ imag(A(i,j));
00184                 }
00185 
00186                 // Array B nach bb kopieren
00187                 for (j=i; j<i+ldbb && j<n; j++)
00188                 {
00189                         // x und y sind die Koordinaten in Lapack-Bandstruktur
00190                         x=ldbb+i-j-1;
00191                         y=j;
00192                         bb[x+ldbb*y].r = CPLX__ real(B(i,j));
00193                         bb[x+ldbb*y].i = CPLX__ imag(B(i,j));
00194                 }
00195         }
00196 
00197         // *******************************************************************************
00198         // Hilfsvariablen fuer die Eigenwertroutine deklarieren
00199 
00200         double* ew = new double[n];  // Temporaerer Eigenwertspeicher
00201         long *iwork = new long[5*n];
00202         double *rwork = new double[7*n];
00203         doublecomplex *work = new doublecomplex[n];
00204         doublecomplex *q = new doublecomplex[n*n];
00205  
00206         // *******************************************************************************
00207         // Eigenwertroutine aufrufen:
00208         long info=own_ew_(&BerechneEV, &uplo,  &n, &Anzahl_Eigenwerte, ab, &ldab, bb, &ldbb,
00209                           &EW_min, &EW_max, ew, q, &n, work, rwork, iwork);
00210         if (!info==0)
00211         {
00212                 STD__ cerr << " Error calculating eigenvalues: (status " << info << ")" << STD__ endl;
00213                 return info;
00214         }
00215         delete [] ab ;
00216         delete [] bb ;
00217         // Wenn Eigenvektoren berechnen
00218         if (BerechneEV==1 && Anzahl_Eigenwerte>0)
00219         {
00220           // neue Hilfsvariablen initialisieren
00221           doublecomplex *EV = new doublecomplex[n*Anzahl_Eigenwerte];
00222           long *ifail = new long[Anzahl_Eigenwerte];
00223           // Eigenvektoren zu 0 setzen
00224           for (i=0;i<n*Anzahl_Eigenwerte;i++)
00225           {
00226             EV[i].r=0;
00227             EV[i].i=0;
00228           }
00229           info=own_ev_(&n, &Anzahl_Eigenwerte, ew, EV, &n, q, &n, work, rwork, iwork, ifail);
00230           if (!info==0)
00231           {
00232                    STD__ cerr << " Error calculating eigenvectors" << STD__ endl;
00233                    return info;
00234           }
00235           // *******************************************************************************
00236           // Eigenvektoren in die Matrix Eigenvektoren eintragen
00237           Eigenvectoren.resize(n,Anzahl_Eigenwerte);
00238           for (j=0;j<Anzahl_Eigenwerte;j++)
00239           {
00240                 for (int i=0; i<n; i++)
00241                 {
00242                 Eigenvectoren(i,j) = CPLX__ complex<double> (EV[i+n*j].r, EV[i+n*j].i);
00243                 }
00244           }
00245           delete[] ifail;
00246           delete[] EV;
00247         } //endif Eigenvektorberechnung
00248 
00249         // Eigenwerte in Eigenwertvektor eintragen
00250         Eigenwerte.resize(Anzahl_Eigenwerte);
00251         for(i=0; i<Anzahl_Eigenwerte; i++)      Eigenwerte(i) = ew[i];
00252 
00253         delete[] ew;
00254         delete[] q;
00255         delete[] iwork;
00256         delete[] rwork;
00257         delete[] work;
00258         return 0;
00259 }
00260 
00261 
00262 
00263 // ***************************************************************************
00265 // ***************************************************************************
00266 int eig(const F_Matrix<CPLX__ complex<double> >& A, Vector<double>& EW, F_Matrix<CPLX__ complex<double> >& EVec)
00267 {
00268         char jobz='V';                  // Eigenvektoren berechnen
00269         char uplo='U';                  // obere Haelfte der Matrix wird gespeichert
00270         integer n=A.rows();             // Ordnung des EWP bestimmen
00271 
00272         integer lwork =2*n-1;
00273         doublecomplex *work = new doublecomplex[lwork];
00274         double *rwork       = new double[3*n-2];
00275         integer info=0;
00276 
00277         EVec=A;
00278         zheev_(&jobz, &uplo, &n, (doublecomplex*) EVec.get_fortran_matrix(), &n,
00279                EW.get_fortran_vector(), work, &lwork, rwork, &info);
00280 
00281         if (!info==0)
00282         {
00283                 STD__ cerr << " Error calculating eigenvalues: (status " << info << ")" << STD__ endl;
00284                 return info;
00285         }
00286 
00287         // Alle alten Arrays loeschen
00288         delete[] work;
00289         delete[] rwork;
00290         return 0;
00291 }
00292 
00293 
00294 int eig(F_Matrix<CPLX__ complex<double> > A, Vector<double>& EW)
00295 {
00296         char jobz='N';                  // Eigenvektoren berechnen
00297         char uplo='U';                  // obere Haelfte der Matrix wird gespeichert
00298         integer n=A.rows();             // Ordnung des EWP bestimmen
00299 
00300         integer lwork =2*n-1;
00301         doublecomplex *work = new doublecomplex[lwork];
00302         double *rwork       = new double[3*n-2];
00303         integer info=0;
00304 
00305         zheev_(&jobz, &uplo, &n, (doublecomplex*) A.get_fortran_matrix(), &n,
00306                EW.get_fortran_vector(), work, &lwork, rwork, &info);
00307 
00308         if (!info==0)
00309         {
00310                 STD__ cerr << " Error calculating eigenvalues: (status " << info << ")" << STD__ endl;
00311                 return info;
00312         }
00313 
00314         // Alle alten Arrays loeschen
00315         delete[] work;
00316         delete[] rwork;
00317         return 0;
00318 }
00319 
00320 template class tbci_memalloc<doublecomplex>;
00321 template class tbci_memalloc_cache<doublecomplex>;
00322 
00323 NAMESPACE_END