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

Timestamp:

Jul 23, 2004, 12:53:35 PM (21 years ago)

Author:

cmv

Message:

add SVD decomp by Divide and Conquer (cmv 23/07/04)

File:

• trunk/SophyaExt/LinAlg/intflapack.cc (modified) (21 diffs)

trunk/SophyaExt/LinAlg/intflapack.cc

@@ -128 +128 @@
                complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* info);
 
+// Driver pour decomposition SVD Divide and Conquer
+  void sgesdd_(char* jobz, int_4* m, int_4* n, r_4* a, int_4* lda,
+               r_4* s, r_4* u, int_4* ldu, r_4* vt, int_4* ldvt, 
+               r_4* work, int_4* lwork, int_4* iwork, int_4* info);
+  void dgesdd_(char* jobz, int_4* m, int_4* n, r_8* a, int_4* lda,
+               r_8* s, r_8* u, int_4* ldu, r_8* vt, int_4* ldvt, 
+               r_8* work, int_4* lwork, int_4* iwork, int_4* info);
+  void cgesdd_(char* jobz, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,
+               r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt, 
+               complex<r_4>* work, int_4* lwork, r_4* rwork, int_4* iwork, int_4* info);
+  void zgesdd_(char* jobz, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,
+               r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt, 
+               complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* iwork, int_4* info);
+
 // Driver pour eigen decomposition for symetric/hermitian matrices
   void ssyev_(char* jobz, char* uplo, int_4* n, r_4* a, int_4* lda, r_4* w,
@@ -253 +267 @@
   \param a : input matrix symetric , overwritten on output
   \param b : input-output, input vector b, contains x on exit
-  \return : return code from lapack driver _gesv()
+  \return : return code from lapack driver
  */
 template <class T>
@@ -347 +361 @@
 {
   if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
-    throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b NbDimensions() != 2"));
+    throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) a Or b NbDimensions() != 2"));
   
   int_4 rowa = a.RowsKA();
@@ -450 +464 @@
 {
   if ( ( a.NbDimensions() != 2 )  )
-    throw(SzMismatchError("LapackServer::SVD(a, ...) a.NbDimensions() != 2"));
+    throw(SzMismatchError("LapackServer::SVDDriver(a, ...) a.NbDimensions() != 2"));
 
   int_4 rowa = a.RowsKA();
@@ -456 +470 @@
 
   if ( !a.IsPacked(rowa) )
-     throw(SzMismatchError("LapackServer::SVD(a, ...) a Not Column Packed "));
+     throw(SzMismatchError("LapackServer::SVDDriver(a, ...) a Not Column Packed "));
 
   int_4 m = a.Size(rowa);
@@ -470 +484 @@
   if ( up != NULL) {
     if ( dynamic_cast< TVector<T> * > (vtp) ) 
-      throw( TypeMismatchExc("LapackServer::SVD() Wrong type (=TVector<T>) for u !") );
+      throw( TypeMismatchExc("LapackServer::SVDDriver() Wrong type (=TVector<T>) for u !") );
     up->SetMemoryMapping(BaseArray::FortranMemoryMapping);
     sz[0] = sz[1] = m;
@@ -482 +496 @@
   if ( vtp != NULL) {
     if ( dynamic_cast< TVector<T> * > (vtp) ) 
-      throw( TypeMismatchExc("LapackServer::SVD() Wrong type (=TVector<T>) for vt !") );
+      throw( TypeMismatchExc("LapackServer::SVDDriver() Wrong type (=TVector<T>) for vt !") );
     vtp->SetMemoryMapping(BaseArray::FortranMemoryMapping);
     sz[0] = sz[1] = n;
@@ -510 +524 @@
   int_4 lwork = maxmn*5*wspace_size_factor;
   T * work = new T[lwork];
-  int_4 info;  
+  int_4 info;
 
   if (typeid(T) == typeid(r_4) ) {
@@ -543 +557 @@
     string tn = typeid(T).name();
     cerr << " LapackServer::SVDDriver(...) - Unsupported DataType T = " << tn << endl;
-    throw TypeMismatchExc("LapackServer::LinSolve(a,b) - Unsupported DataType (T)");
+    throw TypeMismatchExc("LapackServer::SVDDriver(a,b) - Unsupported DataType (T)");
   }
 
   if (jobu == 'N') delete up;
   if (jobvt == 'N') delete vtp;
+  return(info);
+}
+
+
+//! Interface to Lapack SVD driver s/d/c/zgesdd().
+/*! Same as SVD but with Divide and Conquer method */
+template <class T>
+int LapackServer<T>::SVD_DC(TMatrix<T>& a, TVector<T>& s, TMatrix<T>& u, TMatrix<T>& vt)
+{
+
+  if ( !a.IsPacked() )
+     throw(SzMismatchError("LapackServer::SVD_DC(a, ...) a Not Packed "));
+
+  int_4 m = a.NRows();
+  int_4 n = a.NCols();
+  int_4 maxmn = (m > n) ? m : n;
+  int_4 minmn = (m < n) ? m : n;
+  int_4 supermax = 4*minmn*minmn+4*minmn; if(maxmn>supermax) supermax=maxmn;
+
+  char jobz = 'A';
+
+  s.ReSize(minmn);
+  u.ReSize(m,m);
+  vt.ReSize(n,n);
+  
+  int_4 lda = n;
+  int_4 ldu = u.NCols();
+  int_4 ldvt = vt.NCols();
+  int_4 info;  
+
+  if(typeid(T) == typeid(r_4) ) {
+    int_4 lwork = 3*minmn*minmn + supermax;
+    r_4* work = new r_4[lwork +5];
+    int_4* iwork = new int_4[8*minmn +5];
+    sgesdd_(&jobz,&m,&n,(r_4*)a.Data(),&lda,
+           (r_4*)s.Data(),(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt, 
+           (r_4*)work,&lwork,(int_4*)iwork,&info);
+    delete [] work; delete [] iwork;
+  } else if(typeid(T) == typeid(r_8) ) {
+    int_4 lwork = 3*minmn*minmn + supermax;
+    r_8* work = new r_8[lwork +5];
+    int_4* iwork = new int_4[8*minmn +5];
+    dgesdd_(&jobz,&m,&n,(r_8*)a.Data(),&lda,
+           (r_8*)s.Data(),(r_8*)u.Data(),&ldu,(r_8*)vt.Data(),&ldvt, 
+           (r_8*)work,&lwork,(int_4*)iwork,&info);
+    delete [] work; delete [] iwork;
+  } else if(typeid(T) == typeid(complex<r_4>) ) {
+    r_4* sloc = new r_4[minmn];
+    int_4 lwork = minmn*minmn+2*minmn+maxmn;
+    complex<r_4>* work = new complex<r_4>[lwork +5];
+    r_4* rwork = new r_4[5*minmn*minmn+5*minmn +5];
+    int_4* iwork = new int_4[8*minmn +5];
+    cgesdd_(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda,
+           (r_4*)sloc,(complex<r_4>*)u.Data(),&ldu,(complex<r_4>*)vt.Data(),&ldvt, 
+           (complex<r_4>*)work,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
+    for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
+    delete [] sloc; delete [] work; delete [] rwork; delete [] iwork;
+  } else if(typeid(T) == typeid(complex<r_8>) )  {
+    r_8* sloc = new r_8[minmn];
+    int_4 lwork = minmn*minmn+2*minmn+maxmn;
+    complex<r_8>* work = new complex<r_8>[lwork +5];
+    r_8* rwork = new r_8[5*minmn*minmn+5*minmn +5];
+    int_4* iwork = new int_4[8*minmn +5];
+    zgesdd_(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda,
+           (r_8*)sloc,(complex<r_8>*)u.Data(),&ldu,(complex<r_8>*)vt.Data(),&ldvt, 
+           (complex<r_8>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
+    for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
+    delete [] sloc; delete [] work; delete [] rwork; delete [] iwork;
+  } else {
+    string tn = typeid(T).name();
+    cerr << " LapackServer::SVD_DC(...) - Unsupported DataType T = " << tn << endl;
+    throw TypeMismatchExc("LapackServer::SVD_DC - Unsupported DataType (T)");
+  }
+
   return(info);
 }
@@ -559 +647 @@
   \param eigenvector : if true compute eigenvectors, if not only eigenvalues
   \param a : on return array of eigenvectors (same order than eval, one vector = one column)
-  \return : return code from lapack driver _gesvd()
+  \return : return code from lapack driver
  */
 
@@ -574 +662 @@
     throw(SzMismatchError("LapackServer::LapackEigenSym(a,b) a Not Column Packed"));
 
-  char uplo='U'; char struplo[5]; struplo[0]=uplo; struplo[1]='\0';
+  char uplo='U';
   char jobz='N'; if(eigenvector) jobz='V';
-  char strjobz[5]; strjobz[0]=jobz; strjobz[1]='\0';
 
   int_4 n = a.Size(rowa);
@@ -587 +674 @@
     int_4 lwork = 3*n-1 +5; r_4* work = new r_4[lwork];
     r_4* w = new r_4[n];
-    ssyev_(strjobz,struplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info);
+    ssyev_(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info);
     if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
     delete [] work; delete [] w;
@@ -593 +680 @@
     int_4 lwork = 3*n-1 +5; r_8* work = new r_8[lwork];
     r_8* w = new r_8[n];
-    dsyev_(strjobz,struplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info);
+    dsyev_(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info);
     if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
     delete [] work; delete [] w;
@@ -599 +686 @@
     int_4 lwork = 2*n-1 +5; complex<r_4>* work = new complex<r_4>[lwork];
     r_4* rwork = new r_4[3*n-2  +5]; r_4* w = new r_4[n];
-    cheev_(strjobz,struplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
+    cheev_(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
           ,(complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
     if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
@@ -606 +693 @@
     int_4 lwork = 2*n-1 +5; complex<r_8>* work = new complex<r_8>[lwork];
     r_8* rwork = new r_8[3*n-2  +5]; r_8* w = new r_8[n];
-    zheev_(strjobz,struplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
+    zheev_(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
           ,(complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
     if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
@@ -627 +714 @@
   \param eigenvector : if true compute (right) eigenvectors, if not only eigenvalues
   \param a : on return array of eigenvectors
-  \return : return code from lapack driver _gesvd()
+  \return : return code from lapack driver
   \verbatim
   eval : contains the computed eigenvalues.
@@ -656 +743 @@
     throw(SzMismatchError("LapackServer::LapackEigen(a,b) a Not Column Packed"));
 
-  char jobvl = 'N'; char strjobvl[5]; strjobvl[0] = jobvl; strjobvl[1] = '\0';
+  char jobvl = 'N';
   char jobvr = 'N'; if(eigenvector) jobvr='V';
-  char strjobvr[5]; strjobvr[0] = jobvr; strjobvr[1] = '\0';
 
   int_4 n = a.Size(rowa);
@@ -671 +757 @@
     int_4 lwork = 4*n +5; r_4* work = new r_4[lwork];
     r_4* wr = new r_4[n]; r_4* wi = new r_4[n]; r_4* vl = NULL;
-    sgeev_(strjobvl,strjobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
+    sgeev_(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
            (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr,
            (r_4 *)work,&lwork,&info);
@@ -679 +765 @@
     int_4 lwork = 4*n +5; r_8* work = new r_8[lwork];
     r_8* wr = new r_8[n]; r_8* wi = new r_8[n]; r_8* vl = NULL;
-    dgeev_(strjobvl,strjobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
+    dgeev_(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
            (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr,
            (r_8 *)work,&lwork,&info);
@@ -687 +773 @@
     int_4 lwork = 2*n +5; complex<r_4>* work = new complex<r_4>[lwork];
     r_4* rwork = new r_4[2*n+5]; r_4* vl = NULL; TVector< complex<r_4> > w(n);
-    cgeev_(strjobvl,strjobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
+    cgeev_(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
            (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr,
            (complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
@@ -695 +781 @@
     int_4 lwork = 2*n +5; complex<r_8>* work = new complex<r_8>[lwork];
     r_8* rwork = new r_8[2*n+5]; r_8* vl = NULL;
-    zgeev_(strjobvl,strjobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
+    zgeev_(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
            (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr,
            (complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);