Changeset 2559 in Sophya for trunk/SophyaExt/LinAlg

Timestamp:

Jul 22, 2004, 3:38:49 PM (21 years ago)

Author:

cmv

Message:

correction bug svddriver pour complex<r_4> et complex<r_8> (cmv 22/07/04)

File:

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

trunk/SophyaExt/LinAlg/intflapack.cc

@@ -122 +122 @@
                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);
+               r_4* s, complex<r_4>* u, int_4* ldu, complex<r_4>* vt, int_4* ldvt, 
+               complex<r_4>* work, int_4* lwork, r_4* rwork, 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);
+               r_8* s, complex<r_8>* u, int_4* ldu, complex<r_8>* vt, int_4* ldvt, 
+               complex<r_8>* work, int_4* lwork, r_8* rwork, int_4* info);
 
 // Driver pour eigen decomposition for symetric/hermitian matrices
@@ -512 +512 @@
   int_4 info;  
 
-  if (typeid(T) == typeid(r_4) )
+  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) ) 
+  } 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>) ) 
+  } else if (typeid(T) == typeid(complex<r_4>) ) {
+    r_4 * rwork = new r_4[5*minmn +5];
+    r_4 * sloc  = new r_4[minmn];
     cgesvd_(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,
-            (complex<r_4> *)s.Data(), (complex<r_4> *) up->Data(), &ldu, 
+            (r_4 *)sloc, (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>) ) 
+            (complex<r_4> *)work, &lwork, (r_4 *)rwork, &info);
+    for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
+    delete [] rwork; delete [] sloc;
+  } else if (typeid(T) == typeid(complex<r_8>) )  {
+    r_8 * rwork = new r_8[5*minmn +5];
+    r_8 * sloc  = new r_8[minmn];
     zgesvd_(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,
-            (complex<r_8> *)s.Data(), (complex<r_8> *) up->Data(), &ldu, 
+            (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu, 
             (complex<r_8> *)vtp->Data(), &ldvt, 
-            (complex<r_8> *)work, &lwork, &info);
-  else {
+            (complex<r_8> *)work, &lwork, (r_8 *)rwork, &info);
+    for(int_4 i=0;i<minmn;i++) s[i] = sloc[i];
+    delete [] rwork; delete [] sloc;
+  } else {
     if (jobu == 'N') delete up;
     if (jobvt == 'N') delete vtp;