Index: trunk/SophyaExt/LinAlg/intflapack.cc
===================================================================
--- trunk/SophyaExt/LinAlg/intflapack.cc	(revision 2563)
+++ trunk/SophyaExt/LinAlg/intflapack.cc	(revision 2567)
@@ -1,7 +1,10 @@
 #include <iostream>
+#include <math.h>
 #include "intflapack.h"
 #include "tvector.h"
 #include "tmatrix.h"
 #include <typeinfo>
+
+#define GARDMEM 5
 
 /*************** Pour memoire  (Christophe) ***************
@@ -114,4 +117,18 @@
 	      complex<r_8>* b, int_4* ldb, complex<r_8>* work, int_4* lwork, int_4* info);
 
+// Driver pour resolution de systemes au sens de Xi2 par SVD Divide & Conquer
+  void sgelsd_(int_4* m,int_4* n,int_4* nrhs,r_4* a,int_4* lda, 
+	      r_4* b,int_4* ldb,r_4* s,r_4* rcond,int_4* rank,
+              r_4* work,int_4* lwork,int_4* iwork,int_4* info);
+  void dgelsd_(int_4* m,int_4* n,int_4* nrhs,r_8* a,int_4* lda, 
+	      r_8* b,int_4* ldb,r_8* s,r_8* rcond,int_4* rank,
+              r_8* work,int_4* lwork,int_4* iwork,int_4* info);
+  void cgelsd_(int_4* m,int_4* n,int_4* nrhs,complex<r_4>* a,int_4* lda, 
+	      complex<r_4>* b,int_4* ldb,r_4* s,r_4* rcond,int_4* rank,
+              complex<r_4>* work,int_4* lwork,r_4* rwork,int_4* iwork,int_4* info);
+  void zgelsd_(int_4* m,int_4* n,int_4* nrhs,complex<r_8>* a,int_4* lda, 
+	      complex<r_8>* b,int_4* ldb,r_8* s,r_8* rcond,int_4* rank,
+              complex<r_8>* work,int_4* lwork,r_8* rwork,int_4* iwork,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,
@@ -305,22 +322,22 @@
 
   if (typeid(T) == typeid(r_4) ) {
-    lwork = ilaenv_en_C(1,"SSYTRF",struplo,n,-1,-1,-1) * n  +5;
-    work = new T[lwork];
+    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 *)work, &lwork, &info);
   } else if (typeid(T) == typeid(r_8) )  {
-    lwork = ilaenv_en_C(1,"DSYTRF",struplo,n,-1,-1,-1) * n  +5;
-    work = new T[lwork];
+    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 *)work, &lwork, &info);
   } else if (typeid(T) == typeid(complex<r_4>) )  {
-    lwork = ilaenv_en_C(1,"CSYTRF",struplo,n,-1,-1,-1) * n  +5;
-    work = new T[lwork];
+    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> *)work, &lwork, &info);
   } else if (typeid(T) == typeid(complex<r_8>) )  {
-    lwork = ilaenv_en_C(1,"ZSYTRF",struplo,n,-1,-1,-1) * n  +5;
-    work = new T[lwork];
+    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,
@@ -420,4 +437,113 @@
 }
 
+////////////////////////////////////////////////////////////////////////////////////
+//! Interface to Lapack least squares solver driver s/d/c/zgelsd().
+/*! Solves the linear least squares problem defined by an m-by-n matrix 
+  \b a and an m element vector \b b , using SVD factorization Divide and Conquer.
+  Inout arrays should have FortranMemory mapping (column packed).
+  \param rcond : definition of zero value (S(i) <= RCOND*S(0) are treated as zero).
+                 If RCOND < 0, machine precision is used instead.
+  \param a : input matrix, overwritten on output
+  \param b : input vector b overwritten by solution on output (beware of size changing)
+  \param x : output matrix of solutions.
+  \return : return code from lapack driver _gelsd()
+  \warning : b is not resized.
+ */
+template <class T>
+int LapackServer<T>::LeastSquareSolveSVD_DC(TMatrix<T>& a,TMatrix<T>& b,TVector<r_8>& s,int_4& rank,r_8 rcond)
+{
+  if ( ( a.NbDimensions() != 2 ) )
+    throw(SzMismatchError("LapackServer::LeastSquareSolveSVD_DC(a,b) a != 2"));
+  
+  if (!a.IsPacked() || !b.IsPacked())
+     throw(SzMismatchError("LapackServer::LeastSquareSolveSVD_DC(a,b) a Or b Not Packed"));
+
+  int_4 m = a.NRows();
+  int_4 n = a.NCols();
+
+  if(b.NRows() != n)
+     throw(SzMismatchError("LapackServer::LeastSquareSolveSVD_DC(a,b) bad matching dim between a and b"));
+
+  int_4 nrhs = b.NCols();
+  int_4 minmn = (m < n) ? m : n;
+  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);
+  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);
+  }
+  s.ReSize(minmn);
+  
+  int_4 smlsiz = 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);
+  else if(typeid(T) == typeid(complex<r_4>) ) smlsiz = ilaenv_en_C(9,"CGELSD"," ",0,0,0,0);
+  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;
+
+  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;
+    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;
+  } 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];
+    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);
+    r_4* rwork = new r_4[lrwork +GARDMEM];
+    int_4* 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>*)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;
+  } 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);
+    r_8* rwork = new r_8[lrwork +GARDMEM];
+    int_4* iwork = new int_4[3*minmn*nlvl+11*minmn  +GARDMEM];
+    r_8 srcond = rcond;
+    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>*)work,&lwork,(r_8*)rwork,(int_4*)iwork,&info);
+    delete [] work; delete [] rwork; delete [] iwork;
+  } else {
+    string tn = typeid(T).name();
+    cerr << " LapackServer::LeastSquareSolveSVD_DC(a,b) - Unsupported DataType T = " << tn << endl;
+    throw TypeMismatchExc("LapackServer::LeastSquareSolveSVD_DC(a,b) - Unsupported DataType (T)");
+  }
+
+  return(info);
+}
+
 
 ////////////////////////////////////////////////////////////////////////////////////
@@ -521,8 +647,8 @@
   int_4 ldu = up->Step(up->ColsKA());
   int_4 ldvt = vtp->Step(vtp->ColsKA());
+  int_4 info;
 
   int_4 lwork = maxmn*5*wspace_size_factor;
   T * work = new T[lwork];
-  int_4 info;
 
   if (typeid(T) == typeid(r_4) ) {
@@ -535,5 +661,5 @@
 	    (r_8 *)work, &lwork, &info);
   } else if (typeid(T) == typeid(complex<r_4>) ) {
-    r_4 * rwork = new r_4[5*minmn +5];
+    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,
@@ -544,5 +670,5 @@
     delete [] rwork; delete [] sloc;
   } else if (typeid(T) == typeid(complex<r_8>) )  {
-    r_8 * rwork = new r_8[5*minmn +5];
+    r_8 * rwork = new r_8[5*minmn +GARDMEM];
     r_8 * sloc  = new r_8[minmn];
     zgesvd_(&jobu, &jobvt, &m, &n, (complex<r_8> *)a.Data(), &lda,
@@ -595,6 +721,6 @@
     r_4* sloc = new r_4[minmn];
     int_4 lwork = 3*minmn*minmn + supermax;
-    r_4* work = new r_4[lwork +5];
-    int_4* iwork = new int_4[8*minmn +5];
+    r_4* work = new r_4[lwork +GARDMEM];
+    int_4* 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, 
@@ -604,6 +730,6 @@
   } else if(typeid(T) == typeid(r_8) ) {
     int_4 lwork = 3*minmn*minmn + supermax;
-    r_8* work = new r_8[lwork +5];
-    int_4* iwork = new int_4[8*minmn +5];
+    r_8* work = new r_8[lwork +GARDMEM];
+    int_4* 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, 
@@ -613,7 +739,7 @@
     r_4* sloc = new r_4[minmn];
     int_4 lwork = minmn*minmn+2*minmn+maxmn;
-    complex<r_4>* work = new complex<r_4>[lwork +5];
-    r_4* rwork = new r_4[5*minmn*minmn+5*minmn +5];
-    int_4* iwork = new int_4[8*minmn +5];
+    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];
     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, 
@@ -623,7 +749,7 @@
   } 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 +5];
-    r_8* rwork = new r_8[5*minmn*minmn+5*minmn +5];
-    int_4* iwork = new int_4[8*minmn +5];
+    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];
     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, 
@@ -672,5 +798,5 @@
 
   if (typeid(T) == typeid(r_4) ) {
-    int_4 lwork = 3*n-1 +5; r_4* work = new r_4[lwork];
+    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 *)work,&lwork,&info);
@@ -678,5 +804,5 @@
     delete [] work; delete [] w;
   } else if (typeid(T) == typeid(r_8) )  {
-    int_4 lwork = 3*n-1 +5; r_8* work = new r_8[lwork];
+    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 *)work,&lwork,&info);
@@ -684,6 +810,6 @@
     delete [] work; delete [] w;
   } else if (typeid(T) == typeid(complex<r_4>) )  {
-    int_4 lwork = 2*n-1 +5; complex<r_4>* work = new complex<r_4>[lwork];
-    r_4* rwork = new r_4[3*n-2  +5]; r_4* w = new r_4[n];
+    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> *)work,&lwork,(r_4 *)rwork,&info);
@@ -691,6 +817,6 @@
     delete [] work; delete [] rwork; delete [] w;
   } else if (typeid(T) == typeid(complex<r_8>) )  {
-    int_4 lwork = 2*n-1 +5; complex<r_8>* work = new complex<r_8>[lwork];
-    r_8* rwork = new r_8[3*n-2  +5]; r_8* w = new r_8[n];
+    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> *)work,&lwork,(r_8 *)rwork,&info);
@@ -755,5 +881,5 @@
 
   if (typeid(T) == typeid(r_4) ) {
-    int_4 lwork = 4*n +5; r_4* work = new r_4[lwork];
+    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,
@@ -763,5 +889,5 @@
     delete [] work; delete [] wr; delete [] wi;
   } else if (typeid(T) == typeid(r_8) )  {
-    int_4 lwork = 4*n +5; r_8* work = new r_8[lwork];
+    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,
@@ -771,6 +897,6 @@
     delete [] work; delete [] wr; delete [] wi;
   } else if (typeid(T) == typeid(complex<r_4>) )  {
-    int_4 lwork = 2*n +5; complex<r_4>* work = new complex<r_4>[lwork];
-    r_4* rwork = new r_4[2*n+5]; r_4* vl = NULL; TVector< complex<r_4> > w(n);
+    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,
@@ -779,6 +905,6 @@
     delete [] work; delete [] rwork;
   } else if (typeid(T) == typeid(complex<r_8>) )  {
-    int_4 lwork = 2*n +5; complex<r_8>* work = new complex<r_8>[lwork];
-    r_8* rwork = new r_8[2*n+5]; r_8* vl = NULL;
+    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,
Index: trunk/SophyaExt/LinAlg/intflapack.h
===================================================================
--- trunk/SophyaExt/LinAlg/intflapack.h	(revision 2563)
+++ trunk/SophyaExt/LinAlg/intflapack.h	(revision 2567)
@@ -17,4 +17,5 @@
   virtual int LinSolveSym(TArray<T>& a, TArray<T> & b); 
   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.); 
 
   virtual int SVD(TArray<T>& a, TArray<T> & s); 
@@ -66,4 +67,12 @@
 
 /*! \ingroup LinAlg
+    \fn LapackLeastSquareSolveSVD_DC
+    \brief Solves the linear least squares problem A*X = B by SVD
+*/
+template <class T>
+inline int LapackLeastSquareSolveSVD_DC(TMatrix<T>& a,TMatrix<T>& b,TVector<r_8>& s,int_4& rank,r_8 rcond=-1.)
+{ LapackServer<T> lps; return( lps.LeastSquareSolveSVD_DC(a,b,s,rank,rcond) );  }
+
+/*! \ingroup LinAlg
     \fn LapackSVD(TArray<T>&, TArray<T> &)
     \brief SVD decomposition using LapackServer. 
