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

Timestamp:

Jul 19, 2004, 6:09:30 PM (21 years ago)

Author:

cmv

Message:

intro inversion matrices symetriques + ilaenv (cmv 19/07/04)

File:

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

trunk/SophyaExt/LinAlg/intflapack.cc

@@ -52 +52 @@
 
 extern "C" {
+// Le calculateur de workingspace
+  int_4 ilaenv_(int_4 *ispec,char *name,char *opts,int_4 *n1,int_4 *n2,int_4 *n3,int_4 *n4,
+                int_4 nc1,int_4 nc2);
+
 // Drivers pour resolution de systemes lineaires
   void sgesv_(int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
@@ -62 +66 @@
               int_4* ipiv, complex<r_8>* b, int_4* ldb, int_4* info);
 
-  // Driver pour resolution de systemes au sens de Xi2
+// Drivers pour resolution de systemes lineaires symetriques
+  void ssysv_(char* uplo, int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
+              int_4* ipiv, r_4* b, int_4* ldb,
+              r_4* work, int_4* lwork, int_4* info);
+  void dsysv_(char* uplo, int_4* n, int_4* nrhs, r_8* a, int_4* lda, 
+              int_4* ipiv, r_8* b, int_4* ldb,
+              r_8* work, int_4* lwork, int_4* info);
+  void csysv_(char* uplo, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 
+              int_4* ipiv, complex<r_4>* b, int_4* ldb,
+              complex<r_4>* work, int_4* lwork, int_4* info);
+  void zsysv_(char* uplo, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 
+              int_4* ipiv, complex<r_8>* b, int_4* ldb,
+              complex<r_8>* work, int_4* lwork, int_4* info);
+
+// Driver pour resolution de systemes au sens de Xi2
   void sgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, r_4* a, int_4* lda, 
               r_4* b, int_4* ldb, r_4* work, int_4* lwork, int_4* info);
@@ -100 +118 @@
 LapackServer<T>::~LapackServer()
 {
+}
+
+// --- ATTENTION BUG PREVISIBLE (CMV) --- REZA A LIRE S.T.P.
+// -> Cette connerie de Fortran/C interface
+// Dans les routines fortran de lapack:
+// Appel depuis le C avec:
+//   int_4 lwork = -1;
+// SUBROUTINE SSYSV( UPLO,N,NRHS,A,LDA,IPIV,B,LDB,WORK,LWORK,INFO)
+//   INTEGER            INFO, LDA, LDB, LWORK, N, NRHS
+//   LOGICAL            LQUERY
+//   LQUERY = ( LWORK.EQ.-1 )
+//   ELSE IF( LWORK.LT.1 .AND. .NOT.LQUERY ) THEN
+//    ==> le test est bien interprete sous Linux mais pas sous OSF
+//    ==> Sous OSF "LWORK.EQ.-1" est FALSE quand on passe lwork=-1 par argument
+//    ==> POUR REZA: confusion entier 4 / 8 bits ??? (bizarre on l'aurait vu avant?)
+template <class T>
+int_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)
+{
+ int_4 nc1 = strlen(name);
+ int_4 nc2 = strlen(opts);
+ int_4 rc=0;
+ rc = ilaenv_(&ispec,name,opts,&n1,&n2,&n3,&n4,nc1,nc2);
+ //cout<<"ilaenv_en_C("<<ispec<<","<<name<<"("<<nc1<<"),"<<opts<<"("<<nc2<<"),"
+ //    <<n1<<","<<n2<<","<<n3<<","<<n4<<") = "<<rc<<endl;
+ return rc;
 }
 
@@ -150 +193 @@
     throw TypeMismatchExc("LapackServer::LinSolve(a,b) - Unsupported DataType (T)");
   }
+  delete[] ipiv;
+  return(info);
+}
+
+//! 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).
+  \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()
+ */
+template <class T>
+int LapackServer<T>::LinSolveSym(TArray<T>& a, TArray<T> & b)
+// --- REMARQUES DE CMV ---
+// 1./ contrairement a ce qui est dit dans la doc, il s'agit
+//     de matrices SYMETRIQUES complexes et non HERMITIENNES !!!
+// 2./ pourquoi les routines de LinSolve pour des matrices symetriques
+//     sont plus de deux fois plus lentes que les LinSolve generales ???
+{
+  if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
+    throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Or b NbDimensions() != 2"));
+  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::LinSolveSym(a,b) a Not a square Array"));
+  if ( a.Size(rowa) !=  b.Size(rowb)) 
+    throw(SzMismatchError("LapackServer::LinSolveSym(a,b) RowSize(a <> b) "));
+
+  if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
+     throw(SzMismatchError("LapackServer::LinSolveSym(a,b) a Or b Not Column Packed"));
+
+  int_4 n = a.Size(rowa);
+  int_4 nrhs = b.Size(colb);
+  int_4 lda = a.Step(cola);
+  int_4 ldb = b.Step(colb);
+  int_4 info = 0;
+  int_4* ipiv = new int_4[n];
+  int_4 lwork = -1;
+  T * work = NULL;
+
+  char uplo = 'U'; // char uplo = 'L';
+  char struplo[5]; struplo[0] = uplo; struplo[1] = '\0';
+
+  if (typeid(T) == typeid(r_4) ) {
+    lwork = ilaenv_en_C(1,"SSYTRF",struplo,n,-1,-1,-1) * n;
+    work = new T[lwork+5];
+    ssysv_(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,
+          (r_4 *)work, &lwork, &info);
+  } else if (typeid(T) == typeid(r_8) )  {
+    lwork = ilaenv_en_C(1,"DSYTRF",struplo,n,-1,-1,-1) * n;
+    work = new T[lwork+5];
+    dsysv_(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,
+          (r_8 *)work, &lwork, &info);
+  } else if (typeid(T) == typeid(complex<r_4>) )  {
+    lwork = ilaenv_en_C(1,"CSYTRF",struplo,n,-1,-1,-1) * n;
+    work = new T[lwork+5];
+    csysv_(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
+          (complex<r_4> *)b.Data(), &ldb,
+          (complex<r_4> *)work, &lwork, &info);
+  } else if (typeid(T) == typeid(complex<r_8>) )  {
+    lwork = ilaenv_en_C(1,"ZSYTRF",struplo,n,-1,-1,-1) * n;
+    work = new T[lwork+5];
+    zsysv_(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
+          (complex<r_8> *)b.Data(), &ldb,
+          (complex<r_8> *)work, &lwork, &info);
+  } else {
+    delete[] work;
+    delete[] ipiv;
+    string tn = typeid(T).name();
+    cerr << " LapackServer::LinSolveSym(a,b) - Unsupported DataType T = " << tn << endl;
+    throw TypeMismatchExc("LapackServer::LinSolveSym(a,b) - Unsupported DataType (T)");
+  }
+  delete[] work;
   delete[] ipiv;
   return(info);