Changeset 1494 in Sophya for trunk/SophyaExt/LinAlg

Timestamp:

May 15, 2001, 10:16:59 AM (24 years ago)

Author:

ansari

Message:

Ajout LeastSquareSolve - Reza 15/5/2001

Location:

trunk/SophyaExt/LinAlg

Files:

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

trunk/SophyaExt/LinAlg/intflapack.cc

@@ -62 +62 @@
               int_4* ipiv, complex<r_8>* b, int_4* ldb, 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);
+  void dgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, r_8* a, int_4* lda, 
+              r_8* b, int_4* ldb, r_8* work, int_4* lwork, int_4* info);
+  void cgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_4>* a, int_4* lda, 
+              complex<r_4>* b, int_4* ldb, complex<r_4>* work, int_4* lwork, int_4* info);
+  void zgels_(char * trans, int_4* m, int_4* n, int_4* nrhs, complex<r_8>* a, int_4* lda, 
+              complex<r_8>* b, int_4* ldb, complex<r_8>* work, int_4* lwork, 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,
@@ -143 +153 @@
   return(info);
 }
+
+template <class T>
+int LapackServer<T>::LeastSquareSolve(TArray<T>& a, TArray<T> & b)
+{
+  if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
+    throw(SzMismatchError("LapackServer::LinSolve(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) !=  b.Size(rowb)) 
+    throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) RowSize(a <> b) "));
+
+  if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
+     throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b Not Column Packed"));
+
+  if ( a.Size(rowa) <  a.Size(cola)) 
+    throw(SzMismatchError("LapackServer::LeastSquareSolve(a,b) NRows<NCols "));
+
+  int_4 m = a.Size(rowa);
+  int_4 n = a.Size(cola);
+  int_4 nrhs = b.Size(colb);
+
+  int_4 lda = a.Step(cola);
+  int_4 ldb = b.Step(colb);
+  int_4 info;
+
+  int_4 minmn = (m < n) ? m : n;
+  int_4 maxmn = (m > n) ? m : n;
+  int_4 maxmnrhs = (nrhs > maxmn) ? nrhs : maxmn;
+  if (maxmnrhs < 1) maxmnrhs = 1;
+  
+  int_4 lwork = minmn+maxmnrhs*5;
+  T * work = new T[lwork];
+
+  char trans = 'N';
+  
+  if (typeid(T) == typeid(r_4) )
+    sgels_(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda, 
+           (r_4 *)b.Data(), &ldb, (r_4 *)work, &lwork, &info);
+  else if (typeid(T) == typeid(r_8) ) 
+    dgels_(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda, 
+           (r_8 *)b.Data(), &ldb, (r_8 *)work, &lwork, &info);
+  else if (typeid(T) == typeid(complex<r_4>) ) 
+    cgels_(&trans, &m, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, 
+           (complex<r_4> *)b.Data(), &ldb, (complex<r_4> *)work, &lwork, &info);
+  else if (typeid(T) == typeid(complex<r_8>) ) 
+    zgels_(&trans, &m, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, 
+           (complex<r_8> *)b.Data(), &ldb, (complex<r_8> *)work, &lwork, &info);
+  else {
+    delete[] work;
+    string tn = typeid(T).name();
+    cerr << " LapackServer::LeastSquareSolve(a,b) - Unsupported DataType T = " << tn << endl;
+    throw TypeMismatchExc("LapackServer::LeastSquareSolve(a,b) - Unsupported DataType (T)");
+  }
+  delete[] work;
+  return(info);
+}
+
 
 //! Interface to Lapack SVD driver s/d/c/zgesv().

trunk/SophyaExt/LinAlg/intflapack.h

@@ -14 +14 @@
 
   virtual int LinSolve(TArray<T>& a, TArray<T> & b); 
+  virtual int LeastSquareSolve(TArray<T>& a, TArray<T> & b); 
+
   virtual int SVD(TArray<T>& a, TArray<T> & s); 
   virtual int SVD(TArray<T>& a, TArray<T> & s, TArray<T> & u, TArray<T> & vt); 
@@ -40 +42 @@
 { LapackServer<T> lps; return( lps.LinSolve(a, b) );  }
 
+template <class T>
+inline int LapackLeastSquareSolve(TArray<T>& a, TArray<T> & b)
+{ LapackServer<T> lps; return( lps.LeastSquareSolve(a, b) );  }
+
 /*! \ingroup LinAlg
     \fn LapackSVD(TArray<T>&, TArray<T> &)