Changeset 814 in Sophya for trunk/SophyaExt/LinAlg

Timestamp:

Apr 5, 2000, 5:45:21 PM (25 years ago)

Author:

ansari

Message:

Debut LapackServer avec LinSolve - Reza 5/4/2000

Location:

trunk/SophyaExt/LinAlg

Files:

• intflapack.cc (modified) (3 diffs)
• intflapack.h (modified) (1 diff)

trunk/SophyaExt/LinAlg/intflapack.cc

@@ +1 @@
+#include <iostream.h>
 #include "intflapack.h"
+#include <typeinfo>
 
 extern "C" {
@@ -10 +12 @@
 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);
+}
+
+template <class T>
+TArray<T>& LapackServer<T>::LinSolve(TArray<T>& a, TArray<T> & b)
+{
+  if ( ( a.NbDimensions() != 2 ) || ( b.NbDimensions() != 2 ) )
+    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();
+  if ( a.Size(rowa) !=  a.Size(cola)) 
+    throw(SzMismatchError("LapackServer::LinSolve(a,b) a Not a square Array"));
+  if ( a.Size(rowa) !=  b.Size(rowa)) 
+    throw(SzMismatchError("LapackServer::LinSolve(a,b) RowSize(a <> b) "));
+
+  if (!a.IsPacked(rowa) || !b.IsPacked(rowb))
+     throw(SzMismatchError("LapackServer::LinSolve(a,b) a Or b Not Packed Columns "));
+
+  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;
+  int_4* ipiv = new int_4[n];
+
+  if (typeid(T) == typeid(r_4) )
+    sgesv_(&n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb, &info);
+  else if (typeid(T) == typeid(r_8) ) 
+    dgesv_(&n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb, &info);
+  else if (typeid(T) == typeid(complex<r_4>) ) 
+    cgesv_(&n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv, 
+           (complex<r_4> *)b.Data(), &ldb, &info);
+  else if (typeid(T) == typeid(complex<r_8>) ) 
+    zgesv_(&n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv, 
+           (complex<r_8> *)b.Data(), &ldb, &info);
+  else {
+    delete[] ipiv;
+    string tn = typeid(T).name();
+    cerr << " LapackServer::LinSolve(a,b) - Unsupported DataType T = " << tn << endl;
+    throw TypeMismatchExc("LapackServer::LinSolve(a,b) - Unsupported DataType (T)");
+  }
+  delete[] ipiv;
+  return(b);
 }
 
@@ -28 +75 @@
   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__
+#pragma define_template LapackServer<r_4>
+#pragma define_template LapackServer<r_8>
+#pragma define_template LapackServer< complex<r_4> >
+#pragma define_template LapackServer< complex<r_8> >
+#endif
+
+#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
+template class LapackServer<r_4>;
+template class LapackServer<r_8>;
+template class LapackServer< complex<r_4> >;
+template class LapackServer< complex<r_8> >;
+#endif
+
+#if defined(OS_LINUX)
+//  Pour le link avec f2c sous Linux
+extern "C" {
+    void MAIN__();
+}
+
+void MAIN__()
+{
+  cerr << "MAIN__() function for linking with libf2c.a " << endl;
+  cerr << "  This function should never be called !!! " << endl;
+  throw PError("MAIN__() should not be called - see intflapack.cc");
+}
+#endif

trunk/SophyaExt/LinAlg/intflapack.h

@@ -5 +5 @@
 #include "tarray.h"
 
+namespace SOPHYA {
+
+template <class T>
+class LapackServer { 
+public:
+  TArray<T>& LinSolve(TArray<T>& a, TArray<T> & b); 
+};
+
+template <class T>
+inline TArray<T>& LapackLinSolve(TArray<T>& a, TArray<T> & b)
+  { LapackServer<T> lps; return( lps.LinSolve(a, b) );  }
+
+
+} // Fin du namespace
 
 void rztest_lapack(TArray<r_4>& a, TArray<r_4>& b);
-// void rztest_lapack(TMatrix<r_8>& a, TMatrix<r_8>& b);
 
 #endif