Changeset 2572 in Sophya for trunk/SophyaExt/LinAlg

Timestamp:

Jul 28, 2004, 4:21:21 PM (21 years ago)

Author:

cmv

Message:

implementation calcul lwork avec lwork=-1 cmv 28/07/04

Location:

trunk/SophyaExt/LinAlg

Files:

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

trunk/SophyaExt/LinAlg/intflapack.cc

@@ -199 +199 @@
 }
 
-// --- ATTENTION BUG POSSIBLE dans l'avenir (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;
+ int_4 nc1 = strlen(name), nc2 = strlen(opts), 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;
+}
+
+template <class T>
+int_4 LapackServer<T>::type2i4(void *val,int nbytes)
+// Retourne un entier contenant la valeur contenue dans val
+// - nbytes = nombre de bytes dans le contenu de val
+// ex: r_4 x = 3.4;              type2i4(&x,4) -> 3
+// ex: r_8 x = 3.4;              type2i4(&x,8) -> 3
+// ex: complex<r_4> x(3.4,7.8);  type2i4(&x,4) -> 3
+// ex: complex<r_8> x(3.4,7.8);  type2i4(&x,8) -> 3
+{
+  r_4* x4; r_8* x8; int_4 lw=0;
+  if(nbytes==4) {x4 = (r_4*)val; lw = (int_4)(*x4);}
+    else        {x8 = (r_8*)val; lw = (int_4)(*x8);}
+  return lw;
 }
 
@@ -317 +317 @@
   int_4 lwork = -1;
   T * work = NULL;
+  T wkget[2];
 
   char uplo = 'U'; // char uplo = 'L';
@@ -322 +323 @@
 
   if (typeid(T) == typeid(r_4) ) {
-    lwork = ilaenv_en_C(1,"SSYTRF",struplo,n,-1,-1,-1) * n;
-    work = new T[lwork  +GARDMEM];
+    ssysv_(&uplo, &n, &nrhs, (r_4 *)a.Data(), &lda, ipiv, (r_4 *)b.Data(), &ldb,
+          (r_4 *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork  +GARDMEM];
     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  +GARDMEM];
+    dsysv_(&uplo, &n, &nrhs, (r_8 *)a.Data(), &lda, ipiv, (r_8 *)b.Data(), &ldb,
+          (r_8 *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork  +GARDMEM];
     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  +GARDMEM];
+    csysv_(&uplo, &n, &nrhs, (complex<r_4> *)a.Data(), &lda, ipiv,
+          (complex<r_4> *)b.Data(), &ldb,
+          (complex<r_4> *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork  +GARDMEM];
     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  +GARDMEM];
+    zsysv_(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
+          (complex<r_8> *)b.Data(), &ldb,
+          (complex<r_8> *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork  +GARDMEM];
     zsysv_(&uplo, &n, &nrhs, (complex<r_8> *)a.Data(), &lda, ipiv,
           (complex<r_8> *)b.Data(), &ldb,
@@ -410 +417 @@
   if (maxmnrhs < 1) maxmnrhs = 1;
   
-  int_4 lwork = minmn+maxmnrhs*5;
-  T * work = new T[lwork];
+  int_4 lwork = -1; //minmn+maxmnrhs*5;
+  T * work = NULL;
+  T wkget[2];
 
   char trans = 'N';
   
-  if (typeid(T) == typeid(r_4) )
+  if (typeid(T) == typeid(r_4) ) {
+    sgels_(&trans, &m, &n, &nrhs, (r_4 *)a.Data(), &lda, 
+           (r_4 *)b.Data(), &ldb, (r_4 *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork  +GARDMEM];
     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) ) 
+  } else if (typeid(T) == typeid(r_8) )  {
+    dgels_(&trans, &m, &n, &nrhs, (r_8 *)a.Data(), &lda, 
+           (r_8 *)b.Data(), &ldb, (r_8 *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork  +GARDMEM];
     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>) ) 
+  } 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> *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork  +GARDMEM];
     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>) ) 
+  } 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> *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork  +GARDMEM];
     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;
+  } else {
+    if(work) delete [] work; work=NULL;
     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;
+  if(work) delete [] work;
   return(info);
 }
@@ -447 +467 @@
   \param b : input vector b overwritten by solution on output (beware of size changing)
   \param x : output matrix of solutions.
+  \param rank : output the rank of the matrix.
   \return : return code from lapack driver _gelsd()
   \warning : b is not resized.
@@ -462 +483 @@
   int_4 n = a.NCols();
 
-  if(b.NRows() != n)
+  if(b.NRows() != m)
      throw(SzMismatchError("LapackServer::LeastSquareSolveSVD_DC(a,b) bad matching dim between a and b"));
 
@@ -469 +490 @@
   int_4 maxmn = (m > n) ? m : n;
 
-  if(b.NRows() != n)
-     throw(SzMismatchError("LapackServer::LeastSquareSolveSVD_DC(a,b) bad matching dim between a and b"));
-
   int_4 lda = m;
   int_4 ldb = maxmn;
   int_4 info;
 
-  { // Use {} for automatic des-alloc bsave
-  TMatrix<T> bsave(n,nrhs); bsave.SetMemoryMapping(BaseArray::FortranMemoryMapping);
+  { // Use {} for automatic des-allocation of "bsave"
+  TMatrix<T> bsave(m,nrhs); bsave.SetMemoryMapping(BaseArray::FortranMemoryMapping);
   bsave = b;
   b.ReSize(maxmn,nrhs); b = (T) 0;
-  for(int i=0;i<n;i++) for(int j=0;j<nrhs;j++) b(i,j) = bsave(i,j);
-  }
+  for(int i=0;i<m;i++) for(int j=0;j<nrhs;j++) b(i,j) = bsave(i,j);
+  } // Use {} for automatic des-allocation of "bsave"
   s.ReSize(minmn);
   
-  int_4 smlsiz = 25;
+  int_4 smlsiz = 25;  // Normallement ilaenv_en_C(9,...) renvoie toujours 25
   if(typeid(T) == typeid(r_4) )               smlsiz = ilaenv_en_C(9,"SGELSD"," ",0,0,0,0);
   else if(typeid(T) == typeid(r_8) )          smlsiz = ilaenv_en_C(9,"DGELSD"," ",0,0,0,0);
@@ -490 +508 @@
   else if(typeid(T) == typeid(complex<r_8>) ) smlsiz = ilaenv_en_C(9,"ZGELSD"," ",0,0,0,0);
   if(smlsiz<0) smlsiz = 0;
-
   r_8 dum = log((r_8)minmn/(r_8)(smlsiz+1.)) / log(2.);
   int_4 nlvl = int_4(dum) + 1; if(nlvl<0) nlvl = 0;
 
+  T * work = NULL;
+  int_4 * iwork = NULL;
+  int_4 lwork=-1, lrwork;
+  T wkget[2];
+
   if(typeid(T) == typeid(r_4) ) {
-    int_4 lwork = 12*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + minmn*nrhs + (smlsiz+1)*(smlsiz+1);
-    r_4* work = new r_4[lwork +GARDMEM];
-    int_4* iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
     r_4* sloc = new r_4[minmn];
     r_4 srcond = rcond;
+    iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
+    sgelsd_(&m,&n,&nrhs,(r_4*)a.Data(),&lda, 
+           (r_4*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank,
+           (r_4*)wkget,&lwork,(int_4*)iwork,&info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
     sgelsd_(&m,&n,&nrhs,(r_4*)a.Data(),&lda, 
            (r_4*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank,
            (r_4*)work,&lwork,(int_4*)iwork,&info);
     for(int_4 i=0;i<minmn;i++) s(i) = sloc[i];
-    delete [] sloc; delete [] work; delete [] iwork;
+    delete [] sloc;
   } else if(typeid(T) == typeid(r_8) )  {
-    int_4 lwork = 12*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + minmn*nrhs + (smlsiz+1)*(smlsiz+1);
-    r_8* work = new r_8[lwork +GARDMEM];
-    int_4* iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
+    iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
+    dgelsd_(&m,&n,&nrhs,(r_8*)a.Data(),&lda, 
+           (r_8*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank,
+           (r_8*)wkget,&lwork,(int_4*)iwork,&info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
     dgelsd_(&m,&n,&nrhs,(r_8*)a.Data(),&lda, 
            (r_8*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank,
            (r_8*)work,&lwork,(int_4*)iwork,&info);
-    delete [] work; delete [] iwork;
   } else if(typeid(T) == typeid(complex<r_4>) )  {
-    int_4 lwork = 2*minmn + minmn*nrhs;
-    complex<r_4>* work = new complex<r_4>[lwork +GARDMEM];
-    int_4 lrwork = 10*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
+    // Cf meme remarque que ci-dessous (complex<r_8)
+    lrwork = 10*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
+    int_4 lrwork_d = 12*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + minmn*nrhs + (smlsiz+1)*(smlsiz+1);
+    if(lrwork_d > lrwork) lrwork = lrwork_d;
     r_4* rwork = new r_4[lrwork +GARDMEM];
-    int_4* iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
+    iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
     r_4* sloc = new r_4[minmn];
     r_4 srcond = rcond;
     cgelsd_(&m,&n,&nrhs,(complex<r_4>*)a.Data(),&lda, 
            (complex<r_4>*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank,
+           (complex<r_4>*)wkget,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
+    cgelsd_(&m,&n,&nrhs,(complex<r_4>*)a.Data(),&lda, 
+           (complex<r_4>*)b.Data(),&ldb,(r_4*)sloc,&srcond,&rank,
            (complex<r_4>*)work,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
     for(int_4 i=0;i<minmn;i++) s(i) = sloc[i];
-    delete [] sloc; delete [] work; delete [] rwork; delete [] iwork;
+    delete [] sloc; delete [] rwork;
   } else if(typeid(T) == typeid(complex<r_8>) )  {
-    int_4 lwork = 2*minmn + minmn*nrhs;
-    complex<r_8>* work = new complex<r_8>[lwork +GARDMEM];
-    int_4 lrwork = 10*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
+    // CMV: Bizarrement, la formule donnee dans zgelsd() plante pour des N grands (500)
+    // On prend (par analogie) la formule pour "lwork" de dgelsd()
+    lrwork = 10*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + 3*smlsiz*nrhs + (smlsiz+1)*(smlsiz+1);
+    int_4 lrwork_d = 12*minmn + 2*minmn*smlsiz + 8*minmn*nlvl + minmn*nrhs + (smlsiz+1)*(smlsiz+1);
+    if(lrwork_d > lrwork) lrwork = lrwork_d;
     r_8* rwork = new r_8[lrwork +GARDMEM];
-    int_4* iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
-    r_8 srcond = rcond;
+    iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
     zgelsd_(&m,&n,&nrhs,(complex<r_8>*)a.Data(),&lda, 
-           (complex<r_8>*)b.Data(),&ldb,(r_8*)s.Data(),&srcond,&rank,
+           (complex<r_8>*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank,
+           (complex<r_8>*)wkget,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
+    zgelsd_(&m,&n,&nrhs,(complex<r_8>*)a.Data(),&lda, 
+           (complex<r_8>*)b.Data(),&ldb,(r_8*)s.Data(),&rcond,&rank,
            (complex<r_8>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
-    delete [] work; delete [] rwork; delete [] iwork;
+    delete [] rwork;
   } else {
+    if(work) delete [] work; work=NULL;
+    if(iwork) delete [] iwork; iwork=NULL;
     string tn = typeid(T).name();
     cerr << " LapackServer::LeastSquareSolveSVD_DC(a,b) - Unsupported DataType T = " << tn << endl;
@@ -543 +580 @@
   }
 
+  if(work) delete [] work; if(iwork) delete [] iwork;
   return(info);
 }
@@ -572 +610 @@
   \f]
   U and V are orthogonal (unitary) matrices.
-  \param a : input m-by-n matrix (in FotranMemoryMapping)
+  \param a : input m-by-n matrix (in FortranMemoryMapping)
   \param s : Vector of min(m,n) singular values (descending order)
-  \param u : Matrix of left singular vectors
-  \param vt : Transpose of right singular vectors.
+  \param u : m-by-m Matrix of left singular vectors 
+  \param vt : Transpose of right singular vectors (n-by-n matrix).
   \return : return code from lapack driver _gesvd()
  */
@@ -649 +687 @@
   int_4 info;
 
-  int_4 lwork = maxmn*5*wspace_size_factor;
-  T * work = new T[lwork];
+  int_4 lwork = -1;   // maxmn*5  *wspace_size_factor;
+  T * work = NULL;    // = new T[lwork];
+  T wkget[2];
 
   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 *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
     sgesvd_(&jobu, &jobvt, &m, &n, (r_4 *)a.Data(), &lda,
             (r_4 *)s.Data(), (r_4 *) up->Data(), &ldu, (r_4 *)vtp->Data(), &ldvt, 
@@ -659 +702 @@
     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 *)wkget, &lwork, &info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
+    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>) ) {
     r_4 * rwork = new r_4[5*minmn +GARDMEM];
     r_4 * sloc  = new r_4[minmn];
+    cgesvd_(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,
+            (r_4 *)sloc, (complex<r_4> *) up->Data(), &ldu, 
+            (complex<r_4> *)vtp->Data(), &ldvt, 
+            (complex<r_4> *)wkget, &lwork, (r_4 *)rwork, &info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
     cgesvd_(&jobu, &jobvt, &m, &n, (complex<r_4> *)a.Data(), &lda,
             (r_4 *)sloc, (complex<r_4> *) up->Data(), &ldu, 
@@ -675 +727 @@
             (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu, 
             (complex<r_8> *)vtp->Data(), &ldvt, 
+            (complex<r_8> *)wkget, &lwork, (r_8 *)rwork, &info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
+    zgesvd_(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,
+            (r_8 *)sloc, (complex<r_8> *) up->Data(), &ldu, 
+            (complex<r_8> *)vtp->Data(), &ldvt, 
             (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(work) delete [] work; work=NULL;
     if (jobu == 'N') delete up;
     if (jobvt == 'N') delete vtp;
@@ -686 +744 @@
   }
 
+  if(work) delete [] work;
   if (jobu == 'N') delete up;
   if (jobvt == 'N') delete vtp;
@@ -712 +771 @@
   u.ReSize(m,m);
   vt.ReSize(n,n);
-  
-  int_4 lda = n;
-  int_4 ldu = u.NCols();
-  int_4 ldvt = vt.NCols();
-  int_4 info;  
+
+  int_4 lda = m;
+  int_4 ldu = m;
+  int_4 ldvt = n;
+  int_4 info;
+  int_4 lwork=-1;
+  T * work = NULL;
+  int_4 * iwork = NULL;
+  T wkget[2];
 
   if(typeid(T) == typeid(r_4) ) {
     r_4* sloc = new r_4[minmn];
-    int_4 lwork = 3*minmn*minmn + supermax;
-    r_4* work = new r_4[lwork +GARDMEM];
-    int_4* iwork = new int_4[8*minmn +GARDMEM];
+    iwork = new int_4[8*minmn +GARDMEM];
+    sgesdd_(&jobz,&m,&n,(r_4*)a.Data(),&lda,
+           (r_4*)sloc,(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt, 
+           (r_4*)wkget,&lwork,(int_4*)iwork,&info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
     sgesdd_(&jobz,&m,&n,(r_4*)a.Data(),&lda,
            (r_4*)sloc,(r_4*)u.Data(),&ldu,(r_4*)vt.Data(),&ldvt, 
            (r_4*)work,&lwork,(int_4*)iwork,&info);
     for(int_4 i=0;i<minmn;i++) s(i) = (r_8) sloc[i];
-    delete [] sloc; delete [] work; delete [] iwork;
+    delete [] sloc;
   } else if(typeid(T) == typeid(r_8) ) {
-    int_4 lwork = 3*minmn*minmn + supermax;
-    r_8* work = new r_8[lwork +GARDMEM];
-    int_4* iwork = new int_4[8*minmn +GARDMEM];
+    iwork = new int_4[8*minmn +GARDMEM];
+    dgesdd_(&jobz,&m,&n,(r_8*)a.Data(),&lda,
+           (r_8*)s.Data(),(r_8*)u.Data(),&ldu,(r_8*)vt.Data(),&ldvt, 
+           (r_8*)wkget,&lwork,(int_4*)iwork,&info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
     dgesdd_(&jobz,&m,&n,(r_8*)a.Data(),&lda,
            (r_8*)s.Data(),(r_8*)u.Data(),&ldu,(r_8*)vt.Data(),&ldvt, 
            (r_8*)work,&lwork,(int_4*)iwork,&info);
-    delete [] work; delete [] iwork;
   } else if(typeid(T) == typeid(complex<r_4>) ) {
     r_4* sloc = new r_4[minmn];
-    int_4 lwork = minmn*minmn+2*minmn+maxmn;
-    complex<r_4>* work = new complex<r_4>[lwork +GARDMEM];
     r_4* rwork = new r_4[5*minmn*minmn+5*minmn +GARDMEM];
-    int_4* iwork = new int_4[8*minmn +GARDMEM];
+    iwork = new int_4[8*minmn +GARDMEM];
+    cgesdd_(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda,
+           (r_4*)sloc,(complex<r_4>*)u.Data(),&ldu,(complex<r_4>*)vt.Data(),&ldvt, 
+           (complex<r_4>*)wkget,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
+    lwork = type2i4(&wkget[0],4); work = new T[lwork +GARDMEM];
     cgesdd_(&jobz,&m,&n,(complex<r_4>*)a.Data(),&lda,
            (r_4*)sloc,(complex<r_4>*)u.Data(),&ldu,(complex<r_4>*)vt.Data(),&ldvt, 
            (complex<r_4>*)work,&lwork,(r_4*)rwork,(int_4*)iwork,&info);
     for(int_4 i=0;i<minmn;i++) s(i) = (r_8) sloc[i];
-    delete [] sloc; delete [] work; delete [] rwork; delete [] iwork;
+    delete [] sloc; delete [] rwork;
   } else if(typeid(T) == typeid(complex<r_8>) )  {
-    int_4 lwork = minmn*minmn+2*minmn+maxmn;
-    complex<r_8>* work = new complex<r_8>[lwork +GARDMEM];
     r_8* rwork = new r_8[5*minmn*minmn+5*minmn +GARDMEM];
-    int_4* iwork = new int_4[8*minmn +GARDMEM];
+    iwork = new int_4[8*minmn +GARDMEM];
+    zgesdd_(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda,
+           (r_8*)s.Data(),(complex<r_8>*)u.Data(),&ldu,(complex<r_8>*)vt.Data(),&ldvt, 
+           (complex<r_8>*)wkget,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
+    lwork = type2i4(&wkget[0],8); work = new T[lwork +GARDMEM];
     zgesdd_(&jobz,&m,&n,(complex<r_8>*)a.Data(),&lda,
            (r_8*)s.Data(),(complex<r_8>*)u.Data(),&ldu,(complex<r_8>*)vt.Data(),&ldvt, 
            (complex<r_8>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
-    delete [] work; delete [] rwork; delete [] iwork;
+    delete [] rwork;
   } else {
+    if(work) delete [] work; work=NULL;
+    if(iwork) delete [] iwork; iwork=NULL;
     string tn = typeid(T).name();
     cerr << " LapackServer::SVD_DC(...) - Unsupported DataType T = " << tn << endl;
@@ -762 +834 @@
   }
 
+  if(work) delete [] work; if(iwork) delete [] iwork;
   return(info);
 }
@@ -794 +867 @@
   int_4 lda = a.Step(cola);
   int_4 info = 0;
+  int_4 lwork = -1;
+  T * work = NULL;
+  T wkget[2];
 
   b.ReSize(n); b = 0.;
 
   if (typeid(T) == typeid(r_4) ) {
-    int_4 lwork = 3*n-1; r_4* work = new r_4[lwork  +GARDMEM];
     r_4* w = new r_4[n];
+    ssyev_(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)wkget,&lwork,&info);
+    lwork = type2i4(&wkget[0],4); /* 3*n-1;*/ work = new T[lwork  +GARDMEM];
     ssyev_(&jobz,&uplo,&n,(r_4 *)a.Data(),&lda,(r_4 *)w,(r_4 *)work,&lwork,&info);
     if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
-    delete [] work; delete [] w;
+    delete [] w;
   } else if (typeid(T) == typeid(r_8) )  {
-    int_4 lwork = 3*n-1; r_8* work = new r_8[lwork  +GARDMEM];
     r_8* w = new r_8[n];
+    dsyev_(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)wkget,&lwork,&info);
+    lwork = type2i4(&wkget[0],8); /* 3*n-1;*/ work = new T[lwork  +GARDMEM];
     dsyev_(&jobz,&uplo,&n,(r_8 *)a.Data(),&lda,(r_8 *)w,(r_8 *)work,&lwork,&info);
     if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
-    delete [] work; delete [] w;
+    delete [] w;
   } else if (typeid(T) == typeid(complex<r_4>) )  {
-    int_4 lwork = 2*n-1; complex<r_4>* work = new complex<r_4>[lwork  +GARDMEM];
     r_4* rwork = new r_4[3*n-2  +GARDMEM]; r_4* w = new r_4[n];
+    cheev_(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
+          ,(complex<r_4> *)wkget,&lwork,(r_4 *)rwork,&info);
+    lwork = type2i4(&wkget[0],4); /* 2*n-1;*/ work = new T[lwork  +GARDMEM];
     cheev_(&jobz,&uplo,&n,(complex<r_4> *)a.Data(),&lda,(r_4 *)w
           ,(complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
     if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
-    delete [] work; delete [] rwork; delete [] w;
+    delete [] rwork; delete [] w;
   } else if (typeid(T) == typeid(complex<r_8>) )  {
-    int_4 lwork = 2*n-1; complex<r_8>* work = new complex<r_8>[lwork  +GARDMEM];
     r_8* rwork = new r_8[3*n-2  +GARDMEM]; r_8* w = new r_8[n];
+    zheev_(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
+          ,(complex<r_8> *)wkget,&lwork,(r_8 *)rwork,&info);
+    lwork = type2i4(&wkget[0],8); /* 2*n-1;*/ work = new T[lwork  +GARDMEM];
     zheev_(&jobz,&uplo,&n,(complex<r_8> *)a.Data(),&lda,(r_8 *)w
           ,(complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
     if(info==0) for(int i=0;i<n;i++) b(i) = w[i];
-    delete [] work; delete [] rwork; delete [] w;
+    delete [] rwork; delete [] w;
   } else {
+    if(work) delete [] work; work=NULL;
     string tn = typeid(T).name();
     cerr << " LapackServer::LapackEigenSym(a,b) - Unsupported DataType T = " << tn << endl;
@@ -829 +912 @@
   }
 
+  if(work) delete [] work;
   return(info);
 }
@@ -880 +964 @@
   int_4 ldvr = n, ldvl = 1;
 
+  int_4 lwork = -1;
+  T * work = NULL;
+  T wkget[2];
+
   if (typeid(T) == typeid(r_4) ) {
-    int_4 lwork = 4*n; r_4* work = new r_4[lwork  +GARDMEM];
     r_4* wr = new r_4[n]; r_4* wi = new r_4[n]; r_4* vl = NULL;
+    sgeev_(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
+           (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr,
+           (r_4 *)wkget,&lwork,&info);
+    lwork = type2i4(&wkget[0],4); /* 4*n;*/ work = new T[lwork  +GARDMEM];
     sgeev_(&jobvl,&jobvr,&n,(r_4 *)a.Data(),&lda,(r_4 *)wr,(r_4 *)wi,
            (r_4 *)vl,&ldvl,(r_4 *)evec.Data(),&ldvr,
            (r_4 *)work,&lwork,&info);
     if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
-    delete [] work; delete [] wr; delete [] wi;
+    delete [] wr; delete [] wi;
   } else if (typeid(T) == typeid(r_8) )  {
-    int_4 lwork = 4*n; r_8* work = new r_8[lwork  +GARDMEM];
     r_8* wr = new r_8[n]; r_8* wi = new r_8[n]; r_8* vl = NULL;
+    dgeev_(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
+           (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr,
+           (r_8 *)wkget,&lwork,&info);
+    lwork = type2i4(&wkget[0],8); /* 4*n;*/ work = new T[lwork  +GARDMEM];
     dgeev_(&jobvl,&jobvr,&n,(r_8 *)a.Data(),&lda,(r_8 *)wr,(r_8 *)wi,
            (r_8 *)vl,&ldvl,(r_8 *)evec.Data(),&ldvr,
            (r_8 *)work,&lwork,&info);
     if(info==0) for(int i=0;i<n;i++) eval(i) = complex<r_8>(wr[i],wi[i]);
-    delete [] work; delete [] wr; delete [] wi;
+    delete [] wr; delete [] wi;
   } else if (typeid(T) == typeid(complex<r_4>) )  {
-    int_4 lwork = 2*n; complex<r_4>* work = new complex<r_4>[lwork  +GARDMEM];
     r_4* rwork = new r_4[2*n  +GARDMEM]; r_4* vl = NULL; TVector< complex<r_4> > w(n);
+    cgeev_(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
+           (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr,
+           (complex<r_4> *)wkget,&lwork,(r_4 *)rwork,&info);
+    lwork = type2i4(&wkget[0],4); /* 2*n;*/ work = new T[lwork  +GARDMEM];
     cgeev_(&jobvl,&jobvr,&n,(complex<r_4> *)a.Data(),&lda,(complex<r_4> *)w.Data(),
            (complex<r_4> *)vl,&ldvl,(complex<r_4> *)evec.Data(),&ldvr,
            (complex<r_4> *)work,&lwork,(r_4 *)rwork,&info);
     if(info==0) for(int i=0;i<n;i++) eval(i) = w(i);
-    delete [] work; delete [] rwork;
+    delete [] rwork;
   } else if (typeid(T) == typeid(complex<r_8>) )  {
-    int_4 lwork = 2*n; complex<r_8>* work = new complex<r_8>[lwork  +GARDMEM];
     r_8* rwork = new r_8[2*n  +GARDMEM]; r_8* vl = NULL;
     zgeev_(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
            (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr,
+           (complex<r_8> *)wkget,&lwork,(r_8 *)rwork,&info);
+    lwork = type2i4(&wkget[0],8); /* 2*n;*/ work = new T[lwork  +GARDMEM];
+    zgeev_(&jobvl,&jobvr,&n,(complex<r_8> *)a.Data(),&lda,(complex<r_8> *)eval.Data(),
+           (complex<r_8> *)vl,&ldvl,(complex<r_8> *)evec.Data(),&ldvr,
            (complex<r_8> *)work,&lwork,(r_8 *)rwork,&info);
-    delete [] work; delete [] rwork;
+    delete [] rwork;
   } else {
+    if(work) delete [] work; work=NULL;
     string tn = typeid(T).name();
     cerr << " LapackServer::LapackEigen(a,b) - Unsupported DataType T = " << tn << endl;
@@ -917 +1018 @@
   }
 
+  if(work) delete [] work;
   return(info);
 }

trunk/SophyaExt/LinAlg/intflapack.h

@@ -38 +38 @@
                 TArray<T>* up=NULL, TArray<T> * vtp=NULL);
   int_4 ilaenv_en_C(int_4 ispec,char *name,char *opts,int_4 n1,int_4 n2,int_4 n3,int_4 n4);
+  int_4 type2i4(void *val,int nbytes);
 
   int wspace_size_factor;