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

Timestamp:

Nov 24, 2000, 10:47:37 AM (25 years ago)

Author:

ansari

Message:

Ajout de SVD ds LapackServer - Reza 24/11/2000

File:

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

trunk/SophyaExt/LinAlg/intflapack.cc

@@ -1 +1 @@
 #include <iostream.h>
 #include "intflapack.h"
+#include "tvector.h"
+#include "tmatrix.h"
 #include <typeinfo>
 
 extern "C" {
-void sgesv_(int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
-            int_4* ipiv, r_4* b, int_4* ldb, int_4* info);
-void dgesv_(int_4* n, int_4* nrhs, r_8* a, int_4* lda, 
-            int_4* ipiv, r_8* b, int_4* ldb, int_4* info);
-void cgesv_(int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 
-            int_4* ipiv, complex<r_4>* b, int_4* ldb, int_4* info);
-void zgesv_(int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 
-            int_4* ipiv, complex<r_8>* b, int_4* ldb, int_4* info);
+// Drivers pour resolution de systemes lineaires
+  void sgesv_(int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
+              int_4* ipiv, r_4* b, int_4* ldb, int_4* info);
+  void dgesv_(int_4* n, int_4* nrhs, r_8* a, int_4* lda, 
+              int_4* ipiv, r_8* b, int_4* ldb, int_4* info);
+  void cgesv_(int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 
+              int_4* ipiv, complex<r_4>* b, int_4* ldb, int_4* info);
+  void zgesv_(int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 
+              int_4* ipiv, complex<r_8>* b, int_4* ldb, int_4* info);
+
+// Driver pour decomposition SVD 
+  void sgesvd_(char* jobu, char* jobvt, 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* info);
+  void dgesvd_(char* jobu, char* jobvt, 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* info);
+  void cgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_4>* a, int_4* lda,
+               complex<r_4>* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt, 
+               complex<r_4>* work, int_4* lwork, int_4* info);
+  void zgesvd_(char* jobu, char* jobvt, int_4* m, int_4* n, complex<r_8>* a, int_4* lda,
+               complex<r_8>* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt, 
+               complex<r_8>* work, int_4* lwork, int_4* info);
+               
+}
+
+
+//   -------------- Classe LapackServer<T> --------------
+
+template <class T>
+LapackServer<T>::LapackServer<T>()
+{
+  SetWorkSpaceSizeFactor();
+}
+
+template <class T>
+LapackServer<T>::~LapackServer<T>()
+{
 }
 
@@ -20 +52 @@
     throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b NbDimensions() != 2"));
   
-  uint_4 rowa = a.RowsKA();
-  uint_4 cola = a.ColsKA();
-  uint_4 rowb = b.RowsKA();
-  uint_4 colb = b.ColsKA();
+  int_4 rowa = a.RowsKA();
+  int_4 cola = a.ColsKA();
+  int_4 rowb = b.RowsKA();
+  int_4 colb = b.ColsKA();
   if ( a.Size(rowa) !=  a.Size(cola)) 
     throw(SzMismatchError("LapackServer::LinSolve(a,b) a Not a square Array"));
@@ -30 +62 @@
 
   if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
-     throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b Not Packed Columns "));
+     throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b Not Column Packed"));
 
   int_4 n = a.Size(rowa);
@@ -56 +88 @@
   }
   delete[] ipiv;
+  return(info);
+}
+
+template <class T>
+int LapackServer<T>::SVD(TArray<T>& a, TArray<T> & s)
+{
+  return (SVDDriver(a, s, NULL, NULL) );
+}
+
+template <class T>
+int LapackServer<T>::SVD(TArray<T>& a, TArray<T> & s, TArray<T> & u, TArray<T> & vt)
+{
+  return (SVDDriver(a, s, &u, &vt) );
+}
+
+template <class T>
+int LapackServer<T>::SVDDriver(TArray<T>& a, TArray<T> & s, TArray<T>* up, TArray<T>* vtp)
+{
+  if ( ( a.NbDimensions() != 2 )  )
+    throw(SzMismatchError("LapackServer::SVD(a, ...) a.NbDimensions() != 2"));
+
+  int_4 rowa = a.RowsKA();
+  int_4 cola = a.ColsKA();
+
+  if ( !a.IsPacked(rowa) )
+     throw(SzMismatchError("LapackServer::SVD(a, ...) a Not Column Packed "));
+
+  int_4 m = a.Size(rowa);
+  int_4 n = a.Size(cola);
+  int_4 maxmn = (m > n) ? m : n;
+  int_4 minmn = (m < n) ? m : n;
+
+  char jobu, jobvt;
+  jobu = 'N';
+  jobvt = 'N';
+
+  sa_size_t sz[2];
+  if ( up != NULL) {
+    if ( dynamic_cast< TVector<T> * > (vtp) ) 
+      throw( TypeMismatchExc("LapackServer::SVD() Wrong type (=TVector<T>) for u !") );
+    up->SetMemoryMapping(BaseArray::FortranMemoryMapping);
+    sz[0] = sz[1] = m;
+    up->ReSize(2, sz );
+    jobu = 'A';
+  }
+  else { 
+    up = new TMatrix<T>(1,1);
+    jobu = 'N';
+  }
+  if ( vtp != NULL) {
+    if ( dynamic_cast< TVector<T> * > (vtp) ) 
+      throw( TypeMismatchExc("LapackServer::SVD() Wrong type (=TVector<T>) for vt !") );
+    vtp->SetMemoryMapping(BaseArray::FortranMemoryMapping);
+    sz[0] = sz[1] = n;
+    vtp->ReSize(2, sz );
+    jobvt = 'A';
+  }
+  else { 
+    vtp = new TMatrix<T>(1,1);
+    jobvt = 'N';
+  }
+
+  TVector<T> *vs = dynamic_cast< TVector<T> * > (&s);
+  if (vs) vs->ReSize(minmn);
+  else {
+    TMatrix<T> *ms = dynamic_cast< TMatrix<T> * > (&s);
+    if (ms) ms->ReSize(minmn,1);
+    else  { 
+      sz[0] = minmn; sz[1] = 1;
+      s.ReSize(1, sz);
+    }
+  }
+  
+  int_4 lda = a.Step(a.ColsKA());
+  int_4 ldu = up->Step(up->ColsKA());
+  int_4 ldvt = vtp->Step(vtp->ColsKA());
+
+  int_4 lwork = maxmn*5*wspace_size_factor;
+  T * work = new T[lwork];
+  int_4 info;  
+
+  if (typeid(T) == typeid(r_4) )
+    sgesvd_(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda,
+            (r_4 *)s.Data(), (r_4 *) up->Data(), &ldu, (r_4 *)vtp->Data(), &ldvt, 
+            (r_4 *)work, &lwork, &info);
+  else if (typeid(T) == typeid(r_8) ) 
+    dgesvd_(&jobu, &jobvt, &m, &n, (r_8 *)a.Data(), &lda,
+            (r_8 *)s.Data(), (r_8 *) up->Data(), &ldu, (r_8 *)vtp->Data(), &ldvt, 
+            (r_8 *)work, &lwork, &info);
+  else if (typeid(T) == typeid(complex<r_4>) ) 
+    cgesvd_(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,
+            (complex<r_4> *)s.Data(), (complex<r_4> *) up->Data(), &ldu, 
+            (complex<r_4> *)vtp->Data(), &ldvt, 
+            (complex<r_4> *)work, &lwork, &info);
+  else if (typeid(T) == typeid(complex<r_8>) ) 
+    zgesvd_(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,
+            (complex<r_8> *)s.Data(), (complex<r_8> *) up->Data(), &ldu, 
+            (complex<r_8> *)vtp->Data(), &ldvt, 
+            (complex<r_8> *)work, &lwork, &info);
+  else {
+    if (jobu == 'N') delete up;
+    if (jobvt == 'N') delete vtp;
+    string tn = typeid(T).name();
+    cerr << " LapackServer::SVDDriver(...) - Unsupported DataType T = " << tn << endl;
+    throw TypeMismatchExc("LapackServer::LinSolve(a,b) - Unsupported DataType (T)");
+  }
+
+  if (jobu == 'N') delete up;
+  if (jobvt == 'N') delete vtp;
   return(info);
 }