Changeset 2556 in Sophya for trunk/SophyaExt/LinAlg/intflapack.cc


Ignore:
Timestamp:
Jul 21, 2004, 5:52:36 PM (21 years ago)
Author:
cmv
Message:
  • Introduction de l'interface Lapack d'inversion des matrices symetriques
  • Introduction de l'interface Lapack de recherche de valeurs et vecteurs propres (cas general, symetrique et hermitique)
  • Introduction d'un fonction d'interface pour le calculateur de workspace (ilaenv_)
  • Commentaires sur les diverses methodes et sur les matrices FORTRAN
  • Pour tester cf Tests/tsttminv.cc

(cmv, 21/07/04)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaExt/LinAlg/intflapack.cc

    r2554 r2556  
    44#include "tmatrix.h"
    55#include <typeinfo>
     6
     7/*************** Pour memoire  (Christophe) ***************
     8Les dispositions memoires (FORTRAN) pour les vecteurs et matrices LAPACK:
     9
     101./ --- REAL X(N):
     11    if an array X of dimension (N) holds a vector x,
     12    then X(i) holds "x_i" for i=1,...,N
     13
     142./ --- REAL A(LDA,N):
     15    if a two-dimensional array A of dimension (LDA,N) holds an m-by-n matrix A,
     16    then A(i,j) holds "a_ij" for i=1,...,m et j=1,...,n  (LDA must be at least m).
     17    Note that array arguments are usually declared in the software as assumed-size
     18    arrays (last dimension *), for example: REAL A(LDA,*)
     19    --- Rangement en memoire:
     20                            | 11 12 13 14 |
     21    Ex: Real A(4,4):    A = | 21 22 23 24 |
     22                            | 31 32 33 34 |
     23                            | 41 42 43 44 |
     24        memoire: {11 21 31 41} {12 22 32 42} {13 23 33 43} {14 24 34 44}
     25                 First indice (line) "i" varies then the second (column):
     26                 (put all the first column, then put all the second column,
     27                  ..., then put all the last column)
     28***********************************************************/
    629
    730/*!
    … …  
    5174*/
    5275
     76////////////////////////////////////////////////////////////////////////////////////
    5377extern "C" {
    5478// Le calculateur de workingspace
    … …  
    103127               complex<r_8>* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt,
    104128               complex<r_8>* work, int_4* lwork, int_4* info);
    105                
    106 }
    107 
     129
     130// Driver pour eigen decomposition for symetric/hermitian matrices
     131  void ssyev_(char* jobz, char* uplo, int_4* n, r_4* a, int_4* lda, r_4* w,
     132              r_4* work, int_4 *lwork, int_4* info);
     133  void dsyev_(char* jobz, char* uplo, int_4* n, r_8* a, int_4* lda, r_8* w,
     134              r_8* work, int_4 *lwork, int_4* info);
     135  void cheev_(char* jobz, char* uplo, int_4* n, complex<r_4>* a, int_4* lda, r_4* w,
     136              complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info);
     137  void zheev_(char* jobz, char* uplo, int_4* n, complex<r_8>* a, int_4* lda, r_8* w,
     138              complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info);
     139
     140// Driver pour eigen decomposition for general squared matrices
     141  void sgeev_(char* jobl, char* jobvr, int_4* n, r_4* a, int_4* lda, r_4* wr, r_4* wi,
     142              r_4* vl, int_4* ldvl, r_4* vr, int_4* ldvr,
     143              r_4* work, int_4 *lwork, int_4* info);
     144  void dgeev_(char* jobl, char* jobvr, int_4* n, r_8* a, int_4* lda, r_8* wr, r_8* wi,
     145              r_8* vl, int_4* ldvl, r_8* vr, int_4* ldvr,
     146              r_8* work, int_4 *lwork, int_4* info);
     147  void cgeev_(char* jobl, char* jobvr, int_4* n,  complex<r_4>* a, int_4* lda, complex<r_4>* w,
     148              complex<r_4>* vl, int_4* ldvl, complex<r_4>* vr, int_4* ldvr,
     149              complex<r_4>* work, int_4 *lwork, r_4* rwork, int_4* info);
     150  void zgeev_(char* jobl, char* jobvr, int_4* n,  complex<r_8>* a, int_4* lda, complex<r_8>* w,
     151              complex<r_8>* vl, int_4* ldvl, complex<r_8>* vr, int_4* ldvr,
     152              complex<r_8>* work, int_4 *lwork, r_8* rwork, int_4* info);
     153
     154}
    108155
    109156//   -------------- Classe LapackServer<T> --------------
    110157
     158////////////////////////////////////////////////////////////////////////////////////
    111159template <class T>
    112160LapackServer<T>::LapackServer()
    … …  
    120168}
    121169
    122 // --- ATTENTION BUG PREVISIBLE (CMV) --- REZA A LIRE S.T.P.
     170// --- ATTENTION BUG POSSIBLE dans l'avenir (CMV) --- REZA A LIRE S.T.P.
    123171// -> Cette connerie de Fortran/C interface
    124172// Dans les routines fortran de lapack:
    … …  
    133181//    ==> Sous OSF "LWORK.EQ.-1" est FALSE quand on passe lwork=-1 par argument
    134182//    ==> POUR REZA: confusion entier 4 / 8 bits ??? (bizarre on l'aurait vu avant?)
     183////////////////////////////////////////////////////////////////////////////////////
    135184template <class T>
    136185int_4 LapackServer<T>::ilaenv_en_C(int_4 ispec,char *name,char *opts,int_4 n1,int_4 n2,int_4 n3,int_4 n4)
    … …  
    145194}
    146195
     196////////////////////////////////////////////////////////////////////////////////////
    147197//! Interface to Lapack linear system solver driver s/d/c/zgesvd().
    148198/*! Solve the linear system a * x = b. Input arrays
    … …  
    197247}
    198248
     249////////////////////////////////////////////////////////////////////////////////////
    199250//! Interface to Lapack linear system solver driver s/d/c/zsysvd().
    200251/*! Solve the linear system a * x = b with a symetric. Input arrays
    … …  
    210261//     de matrices SYMETRIQUES complexes et non HERMITIENNES !!!
    211262// 2./ pourquoi les routines de LinSolve pour des matrices symetriques
    212 //     sont plus de deux fois plus lentes que les LinSolve generales ???
     263//     sont plus de deux fois plus lentes que les LinSolve generales sur OSF
     264//     et sensiblement plus lentes sous Linux ???
    213265{
    214266  if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
    … …  
    239291
    240292  if (typeid(T) == typeid(r_4) ) {
    241     lwork = ilaenv_en_C(1,"SSYTRF",struplo,n,-1,-1,-1) * n;
    242     work = new T[lwork+5];
     293    lwork = ilaenv_en_C(1,"SSYTRF",struplo,n,-1,-1,-1) * n  +5;
     294    work = new T[lwork];
    243295    ssysv_(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,
    244296          (r_4 *)work, &lwork, &info);
    245297  } else if (typeid(T) == typeid(r_8) )  {
    246     lwork = ilaenv_en_C(1,"DSYTRF",struplo,n,-1,-1,-1) * n;
    247     work = new T[lwork+5];
     298    lwork = ilaenv_en_C(1,"DSYTRF",struplo,n,-1,-1,-1) * n  +5;
     299    work = new T[lwork];
    248300    dsysv_(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,
    249301          (r_8 *)work, &lwork, &info);
    250302  } else if (typeid(T) == typeid(complex<r_4>) )  {
    251     lwork = ilaenv_en_C(1,"CSYTRF",struplo,n,-1,-1,-1) * n;
    252     work = new T[lwork+5];
     303    lwork = ilaenv_en_C(1,"CSYTRF",struplo,n,-1,-1,-1) * n  +5;
     304    work = new T[lwork];
    253305    csysv_(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
    254306          (complex<r_4> *)b.Data(), &ldb,
    255307          (complex<r_4> *)work, &lwork, &info);
    256308  } else if (typeid(T) == typeid(complex<r_8>) )  {
    257     lwork = ilaenv_en_C(1,"ZSYTRF",struplo,n,-1,-1,-1) * n;
    258     work = new T[lwork+5];
     309    lwork = ilaenv_en_C(1,"ZSYTRF",struplo,n,-1,-1,-1) * n  +5;
     310    work = new T[lwork];
    259311    zsysv_(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
    260312          (complex<r_8> *)b.Data(), &ldb,
    261313          (complex<r_8> *)work, &lwork, &info);
    262314  } else {
    263     delete[] work;
     315    if(work) delete[] work;
    264316    delete[] ipiv;
    265317    string tn = typeid(T).name();
    … …  
    267319    throw TypeMismatchExc("LapackServer::LinSolveSym(a,b) - Unsupported DataType (T)");
    268320  }
    269   delete[] work;
     321  if(work) delete[] work;
    270322  delete[] ipiv;
    271323  return(info);
    272324}
    273325
     326////////////////////////////////////////////////////////////////////////////////////
    274327//! Interface to Lapack least squares solver driver s/d/c/zgels().
    275328/*! Solves the linear least squares problem defined by an m-by-n matrix
    … …  
    354407
    355408
     409////////////////////////////////////////////////////////////////////////////////////
    356410//! Interface to Lapack SVD driver s/d/c/zgesv().
    357411/*! Computes the vector of singular values of \b a. Input arrays
    … …  
    489543}
    490544
     545
     546////////////////////////////////////////////////////////////////////////////////////
     547/*! Computes the eigen values and eigen vectors of a symetric (or hermitian) matrix \b a.
     548  Input arrays should have FortranMemoryMapping (column packed).
     549  \param a : input symetric (or hermitian) n-by-n matrix
     550  \param b : Vector of eigenvalues (descending order)
     551  \param eigenvector : if true compute eigenvectors, if not only eigenvalues
     552  \param a : on return array of eigenvectors (same order than eval, one vector = one column)
     553  \return : return code from lapack driver _gesvd()
     554 */
     555
     556template <class T>
     557int LapackServer<T>::LapackEigenSym(TArray<T>& a, TVector<r_8>& b, bool eigenvector)
     558{
     559  if ( a.NbDimensions() != 2 )
     560    throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a NbDimensions() != 2"));
     561  int_4 rowa = a.RowsKA();
     562  int_4 cola = a.ColsKA();
     563  if ( a.Size(rowa) !=  a.Size(cola))
     564    throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not a square Array"));
     565  if (!a.IsPacked(rowa))
     566    throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not Column Packed"));
     567
     568  char uplo='U'; char struplo[5]; struplo[0]=uplo; struplo[1]='\0';
     569  char jobz='N'; if(eigenvector) jobz='V';
     570  char strjobz[5]; strjobz[0]=jobz; strjobz[1]='\0';
     571
     572  int_4 n = a.Size(rowa);
     573  int_4 lda = a.Step(cola);
     574  int_4 info = 0;
     575
     576  b.ReSize(n); b = 0.;
     577
     578  if (typeid(T) == typeid(r_4) ) {
     579    int_4 lwork = 3*n-1 +5; r_4* work = new r_4[lwork];
     580    r_4* w = new r_4[n];
     581    ssyev_(strjobz,struplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info);
     582    if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
     583    delete [] work; delete [] w;
     584  } else if (typeid(T) == typeid(r_8) )  {
     585    int_4 lwork = 3*n-1 +5; r_8* work = new r_8[lwork];
     586    r_8* w = new r_8[n];
     587    dsyev_(strjobz,struplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info);
     588    if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
     589    delete [] work; delete [] w;
     590  } else if (typeid(T) == typeid(complex<r_4>) )  {
     591    int_4 lwork = 2*n-1 +5; complex<r_4>* work = new complex<r_4>[lwork];
     592    r_4* rwork = new r_4[3*n-2  +5]; r_4* w = new r_4[n];
     593    cheev_(strjobz,struplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
     594          ,(complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
     595    if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
     596    delete [] work; delete [] rwork; delete [] w;
     597  } else if (typeid(T) == typeid(complex<r_8>) )  {
     598    int_4 lwork = 2*n-1 +5; complex<r_8>* work = new complex<r_8>[lwork];
     599    r_8* rwork = new r_8[3*n-2  +5]; r_8* w = new r_8[n];
     600    zheev_(strjobz,struplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
     601          ,(complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
     602    if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
     603    delete [] work; delete [] rwork; delete [] w;
     604  } else {
     605    string tn = typeid(T).name();
     606    cerr << " LapackServer::LapackEigenSym(a,b) - Unsupported DataType T = " << tn << endl;
     607    throw TypeMismatchExc("LapackServer::LapackEigenSym(a,b) - Unsupported DataType (T)");
     608  }
     609
     610  return(info);
     611}
     612
     613////////////////////////////////////////////////////////////////////////////////////
     614/*! Computes the eigen values and eigen vectors of a general squared matrix \b a.
     615  Input arrays should have FortranMemoryMapping (column packed).
     616  \param a : input general n-by-n matrix
     617  \param eval : Vector of eigenvalues (complex double precision)
     618  \param evec : Matrix of eigenvector (same order than eval, one vector = one column)
     619  \param eigenvector : if true compute (right) eigenvectors, if not only eigenvalues
     620  \param a : on return array of eigenvectors
     621  \return : return code from lapack driver _gesvd()
     622  \verbatim
     623  eval : contains the computed eigenvalues.
     624         --- For real matrices "a" :
     625         Complex conjugate pairs of eigenvalues appear consecutively
     626         with the eigenvalue having the positive imaginary part first.
     627  evec : the right eigenvectors v(j) are stored one after another
     628         in the columns of evec, in the same order as their eigenvalues.
     629         --- For real matrices "a" :
     630         If the j-th eigenvalue is real, then v(j) = evec(:,j),
     631           the j-th column of evec.
     632         If the j-th and (j+1)-st eigenvalues form a complex
     633           conjugate pair, then v(j) = evec(:,j) + i*evec(:,j+1) and
     634           v(j+1) = evec(:,j) - i*evec(:,j+1).
     635  \endverbatim
     636*/
     637
     638template <class T>
     639int LapackServer<T>::LapackEigen(TArray<T>& a, TVector< complex<r_8> >& eval, TMatrix<T>& evec, bool eigenvector)
     640{
     641  if ( a.NbDimensions() != 2 )
     642    throw(SzMismatchError("LapackServer::LapackEigen(a,b) a NbDimensions() != 2"));
     643  int_4 rowa = a.RowsKA();
     644  int_4 cola = a.ColsKA();
     645  if ( a.Size(rowa) !=  a.Size(cola))
     646    throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not a square Array"));
     647  if (!a.IsPacked(rowa))
     648    throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not Column Packed"));
     649
     650  char jobvl = 'N'; char strjobvl[5]; strjobvl[0] = jobvl; strjobvl[1] = '\0';
     651  char jobvr = 'N'; if(eigenvector) jobvr='V';
     652  char strjobvr[5]; strjobvr[0] = jobvr; strjobvr[1] = '\0';
     653
     654  int_4 n = a.Size(rowa);
     655  int_4 lda = a.Step(cola);
     656  int_4 info = 0;
     657
     658  eval.ReSize(n); eval = complex<r_8>(0.,0.);
     659  if(eigenvector) {evec.ReSize(n,n); evec = (T) 0.;}
     660  int_4 ldvr = n, ldvl = 1;
     661
     662  if (typeid(T) == typeid(r_4) ) {
     663    int_4 lwork = 4*n +5; r_4* work = new r_4[lwork];
     664    r_4* wr = new r_4[n]; r_4* wi = new r_4[n]; r_4* vl = NULL;
     665    sgeev_(strjobvl,strjobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
     666           (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr,
     667           (r_4 *)work,&lwork,&info);
     668    if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
     669    delete [] work; delete [] wr; delete [] wi;
     670  } else if (typeid(T) == typeid(r_8) )  {
     671    int_4 lwork = 4*n +5; r_8* work = new r_8[lwork];
     672    r_8* wr = new r_8[n]; r_8* wi = new r_8[n]; r_8* vl = NULL;
     673    dgeev_(strjobvl,strjobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
     674           (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr,
     675           (r_8 *)work,&lwork,&info);
     676    if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
     677    delete [] work; delete [] wr; delete [] wi;
     678  } else if (typeid(T) == typeid(complex<r_4>) )  {
     679    int_4 lwork = 2*n +5; complex<r_4>* work = new complex<r_4>[lwork];
     680    r_4* rwork = new r_4[2*n+5]; r_4* vl = NULL; TVector< complex<r_4> > w(n);
     681    cgeev_(strjobvl,strjobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
     682           (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr,
     683           (complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
     684    if(info==0) for(int i=0;i<n;i++) eval(i) = w(i);
     685    delete [] work; delete [] rwork;
     686  } else if (typeid(T) == typeid(complex<r_8>) )  {
     687    int_4 lwork = 2*n +5; complex<r_8>* work = new complex<r_8>[lwork];
     688    r_8* rwork = new r_8[2*n+5]; r_8* vl = NULL;
     689    zgeev_(strjobvl,strjobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
     690           (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr,
     691           (complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
     692    delete [] work; delete [] rwork;
     693  } else {
     694    string tn = typeid(T).name();
     695    cerr << " LapackServer::LapackEigen(a,b) - Unsupported DataType T = " << tn << endl;
     696    throw TypeMismatchExc("LapackServer::LapackEigen(a,b) - Unsupported DataType (T)");
     697  }
     698
     699  return(info);
     700}
     701
     702
     703
     704
     705
     706////////////////////////////////////////////////////////////////////////////////////
    491707void rztest_lapack(TArray<r_4>& aa, TArray<r_4>& bb)
    492708{
Note: See TracChangeset for help on using the changeset viewer.