Changeset 2563 in Sophya

Timestamp:

Jul 26, 2004, 12:55:16 PM (21 years ago)

Author:

cmv

Message:

change arg SVD_DC -> TVector<r_8> cmv 26/07/04

Location:

trunk/SophyaExt/LinAlg

Files:

• intflapack.cc (modified) (7 diffs)
• intflapack.h (modified) (2 diffs)

trunk/SophyaExt/LinAlg/intflapack.cc

@@ -209 +209 @@
 
 ////////////////////////////////////////////////////////////////////////////////////
-//! Interface to Lapack linear system solver driver s/d/c/zgesvd().
-/*! Solve the linear system a * x = b. Input arrays
-  should have FortranMemory mapping (column packed).
+//! Interface to Lapack linear system solver driver s/d/c/zgesv().
+/*! Solve the linear system a * x = b using LU factorization.
+  Input arrays should have FortranMemory mapping (column packed).
   \param a : input matrix, overwritten on output
   \param b : input-output, input vector b, contains x on exit
@@ -262 +262 @@
 
 ////////////////////////////////////////////////////////////////////////////////////
-//! Interface to Lapack linear system solver driver s/d/c/zsysvd().
-/*! Solve the linear system a * x = b with a symetric. Input arrays
-  should have FortranMemory mapping (column packed).
+//! Interface to Lapack linear system solver driver s/d/c/zsysv().
+/*! Solve the linear system a * x = b with a symetric matrix using LU factorization.
+  Input arrays should have FortranMemory mapping (column packed).
   \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
+  \return : return code from lapack driver _sysv()
  */
 template <class T>
@@ -341 +341 @@
 //! Interface to Lapack least squares solver driver s/d/c/zgels().
 /*! Solves the linear least squares problem defined by an m-by-n matrix 
-  \b a and an m element vector \b b .
+  \b a and an m element vector \b b , using QR or LQ factorization .
   A solution \b x to the overdetermined system of linear equations
   b = a * x is computed, minimizing the norm of b-a*x. 
@@ -569 +569 @@
 /*! 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)
+int LapackServer<T>::SVD_DC(TMatrix<T>& a, TVector<r_8>& s, TMatrix<T>& u, TMatrix<T>& vt)
 {
 
@@ -593 +593 @@
 
   if(typeid(T) == typeid(r_4) ) {
+    r_4* sloc = new r_4[minmn];
     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*)sloc,(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt, 
            (r_4*)work,&lwork,(int_4*)iwork,&info);
-    delete [] work; delete [] iwork;
+    for(int_4 i=0;i<minmn;i++) s(i) = (r_8) sloc[i];
+    delete [] sloc; delete [] work; delete [] iwork;
   } else if(typeid(T) == typeid(r_8) ) {
     int_4 lwork = 3*minmn*minmn + supermax;
@@ -617 +619 @@
            (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];
+    for(int_4 i=0;i<minmn;i++) s(i) = (r_8) 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];
@@ -626 +627 @@
     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, 
+           (r_8*)s.Data(),(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;
+    delete [] work; delete [] rwork; delete [] iwork;
   } else {
     string tn = typeid(T).name();

trunk/SophyaExt/LinAlg/intflapack.h

@@ -20 +20 @@
   virtual int SVD(TArray<T>& a, TArray<T> & s); 
   virtual int SVD(TArray<T>& a, TArray<T> & s, TArray<T> & u, TArray<T> & vt);
-  virtual int SVD_DC(TMatrix<T>& a, TVector<T>& s, TMatrix<T>& u, TMatrix<T>& vt);
+  virtual int SVD_DC(TMatrix<T>& a, TVector<r_8>& s, TMatrix<T>& u, TMatrix<T>& vt);
  
   virtual int LapackEigenSym(TArray<T>& a, TVector<r_8>& b, bool eigenvector=true);
@@ -84 +84 @@
 
 /*! \ingroup LinAlg
-    \fn LapackSVD_DC(TMatrix<T>&, TVector<T>&, TMatrix<T>&, TMatrix<T>&)
+    \fn LapackSVD_DC(TMatrix<T>&, TVector<r_8>&, TMatrix<T>&, TMatrix<T>&)
     \brief SVD decomposition DC using LapackServer. 
 */
 template <class T>
-inline int LapackSVD_DC(TMatrix<T>& a, TVector<T>& s, TMatrix<T>& u, TMatrix<T>& vt)
+inline int LapackSVD_DC(TMatrix<T>& a, TVector<r_8>& s, TMatrix<T>& u, TMatrix<T>& vt)
 { LapackServer<T> lps; return( lps.SVD_DC(a, s, u, vt) ); }