Changeset 2646 in Sophya for trunk/SophyaExt/LinAlg

Timestamp:

Feb 7, 2005, 5:42:27 PM (21 years ago)

Author:

ansari

Message:

Ajout methode calcul de la matrice inverse par Lapack en utilisant la resolution de systeme + MAJ doc - Reza 7 Fev 2005

Location:

trunk/SophyaExt/LinAlg

Files:

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

trunk/SophyaExt/LinAlg/intflapack.cc

@@ -48 +48 @@
   its is available from http://www.netlib.org.
   
-  The present version of our LapackServer (Feb 2001) provides only 
-  interfaces for the linear system solver and singular value 
-  decomposition (SVD). Only arrays with BaseArray::FortranMemoryMapping 
+  The present version of LapackServer (Feb 2005) provides  
+  interfaces for the linear system solver, singular value 
+  decomposition (SVD), Least square solver and  
+  eigen value / eigen vector decomposition. 
+  Only arrays with BaseArray::FortranMemoryMapping 
   can be handled by LapackServer. LapackServer can be instanciated 
   for simple and double precision real or complex array types.
-
+  \warning The input array is overwritten in most cases.
   The example below shows solving a linear system A*X = B
 
@@ -456 +458 @@
   if(work) delete [] work;
   return(info);
+}
+
+////////////////////////////////////////////////////////////////////////////////////
+//! Square matrix inversion using Lapack linear system solver
+/*! Compute the inverse of a square matrix using linear system solver routine
+  Input arrays should have FortranMemory mapping (column packed).
+  \param a : input matrix, overwritten on output
+  \param ainv : output matrix, contains inverse(a) on exit. 
+  ainv is allocated if it has size 0 
+  If not allocated, ainv is automatically 
+  \return : return code from LapackServer::LinSolve()
+  \sa LapackServer::LinSolve()
+ */
+template <class T>
+int LapackServer<T>::ComputeInverse(TMatrix<T>& a, TMatrix<T> & ainv)
+{  
+  if ( a.NbDimensions() != 2 )
+    throw(SzMismatchError("LapackServer::Inverse() NDim(a) != 2"));
+  if ( a.GetMemoryMapping() != BaseArray::FortranMemoryMapping )
+    throw(SzMismatchError("LapackServer::Inverse() a NOT in FortranMemoryMapping"));
+  if ( a.NRows() != a.NCols() )
+    throw(SzMismatchError("LapackServer::Inverse() a NOT square matrix (a.NRows!=a.NCols)"));
+  if (ainv.IsAllocated()) {
+    bool smo, ssz;
+    ssz = a.CompareSizes(ainv, smo);
+    if ( (ssz == false) || (smo == false) ) 
+      throw(SzMismatchError("LapackServer::Inverse() ainv<>a Size/MemOrg mismatch "));   
+  }
+  else  ainv.SetSize(a.NRows(), a.NCols(), BaseArray::FortranMemoryMapping, false);
+  ainv = IdentityMatrix(); 
+  return LinSolve(a, ainv);
 }
 
@@ -1026 +1059 @@
 
 
-
-////////////////////////////////////////////////////////////////////////////////////
-void rztest_lapack(TArray<r_4>& aa, TArray<r_4>& bb)
-{
-  if ( aa.NbDimensions() != 2 ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a Matrix"));
-  if ( aa.SizeX() !=  aa.SizeY()) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a square Matrix"));
-  if ( bb.NbDimensions() != 2 ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A Not a Matrix"));
-  if ( bb.SizeX() !=  aa.SizeX() ) throw(SzMismatchError("rztest_lapack(TMatrix<r_4> A <> B "));
-  if ( !bb.IsPacked() || !bb.IsPacked() ) 
-    throw(SzMismatchError("rztest_lapack(TMatrix<r_4> Not packed A or B "));
-  
-  int_4 n = aa.SizeX();
-  int_4 nrhs = bb.SizeY();
-  int_4 lda = n;
-  int_4 ldb = bb.SizeX();
-  int_4 info;
-  int_4* ipiv = new int_4[n];
-  sgesv_(&n, &nrhs, aa.Data(), &lda, ipiv, bb.Data(), &ldb, &info);
-  delete[] ipiv;
-  cout << "rztest_lapack/Info= " << info << endl;
-  cout << aa << "\n" << bb << endl;
-  return;
-}
-
 ///////////////////////////////////////////////////////////////
 #ifdef __CXX_PRAGMA_TEMPLATES__

trunk/SophyaExt/LinAlg/intflapack.h

@@ -18 +18 @@
   virtual int LeastSquareSolve(TArray<T>& a, TArray<T> & b); 
   virtual int LeastSquareSolveSVD_DC(TMatrix<T>& a,TMatrix<T>& b,TVector<r_8>& s,int_4& rank,r_8 rcond=-1.); 
+
+  // Calcul de la matrice inverse en utilisant la resolution de syst. lineaire
+  virtual int ComputeInverse(TMatrix<T>& a, TMatrix<T>& ainv);
 
   virtual int SVD(TArray<T>& a, TArray<T> & s); 
@@ -66 +69 @@
 inline int LapackLeastSquareSolve(TArray<T>& a, TArray<T> & b)
 { LapackServer<T> lps; return( lps.LeastSquareSolve(a, b) );  }
+
+/*! \ingroup LinAlg
+    \fn LapackInverse(TMatrix<T>&)
+    \brief Computes the inverse matrix using linear system solver LapackServer::LinSolve. 
+*/
+template <class T>
+inline TMatrix<T> LapackInverse(TMatrix<T>& a)
+{ LapackServer<T> lps; TMatrix<T> ainv; lps.ComputeInverse(a, ainv);  return ainv; }
 
 /*! \ingroup LinAlg