2 files changed, 41 insertions, 55 deletions

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 80e5720..1018126 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -114,12 +114,12 @@ for(k=0; k<(M##r * M##c); k++) {             \
 ////////////////////////////////////////////////////////////////////////////////
 //////////////////// real svd ///////////////////////////////////////////////////
 
-/* Subroutine */ int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
+int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
         doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
         ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
         integer *info);
 
-int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
+int svd_l_R(ODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
@@ -145,15 +145,12 @@ int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
         }
     }
     DEBUGMSG("svd_l_R");
-    double *B = (double*)malloc(m*n*sizeof(double));
-    CHECK(!B,MEM);
-    memcpy(B,ap,m*n*sizeof(double));
     integer lwork = -1;
     integer res;
     // ask for optimal lwork
     double ans;
     dgesvd_ (jobu,jobvt,
-             &m,&n,B,&m,
+             &m,&n,ap,&m,
              sp,
              up,&m,
              vp,&ldvt,
@@ -163,7 +160,7 @@ int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
     double * work = (double*)malloc(lwork*sizeof(double));
     CHECK(!work,MEM);
     dgesvd_ (jobu,jobvt,
-             &m,&n,B,&m,
+             &m,&n,ap,&m,
              sp,
              up,&m,
              vp,&ldvt,
@@ -171,18 +168,17 @@ int svd_l_R(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
              &res);
     CHECK(res,res);
     free(work);
-    free(B);
     OK
 }
 
 // (alternative version)
 
-/* Subroutine */ int dgesdd_(char *jobz, integer *m, integer *n, doublereal *
+int dgesdd_(char *jobz, integer *m, integer *n, doublereal *
         a, integer *lda, doublereal *s, doublereal *u, integer *ldu,
         doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
         integer *iwork, integer *info);
 
-int svd_l_Rdd(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
+int svd_l_Rdd(ODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
@@ -202,37 +198,31 @@ int svd_l_Rdd(KODMAT(a),ODMAT(u), DVEC(s),ODMAT(v)) {
         }
     }
     DEBUGMSG("svd_l_Rdd");
-    double *B = (double*)malloc(m*n*sizeof(double));
-    CHECK(!B,MEM);
-    memcpy(B,ap,m*n*sizeof(double));
     integer* iwk = (integer*) malloc(8*q*sizeof(integer));
     CHECK(!iwk,MEM);
     integer lwk = -1;
     integer res;
     // ask for optimal lwk
     double ans;
-    dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
+    dgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,iwk,&res);
     lwk = ans;
     double * workv = (double*)malloc(lwk*sizeof(double));
     CHECK(!workv,MEM);
-    dgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
+    dgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,iwk,&res);
     CHECK(res,res);
     free(iwk);
     free(workv);
-    free(B);
     OK
 }
 
 //////////////////// complex svd ////////////////////////////////////
 
-// not in clapack.h
-
 int zgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
     doublecomplex *a, integer *lda, doublereal *s, doublecomplex *u,
     integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
     integer *lwork, doublereal *rwork, integer *info);
 
-int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
+int svd_l_C(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
@@ -257,10 +247,7 @@ int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
             ldvt = q;
         }
     }DEBUGMSG("svd_l_C");
-    doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
-    CHECK(!B,MEM);
-    memcpy(B,ap,m*n*sizeof(doublecomplex));
-
+    
     double *rwork = (double*) malloc(5*q*sizeof(double));
     CHECK(!rwork,MEM);
     integer lwork = -1;
@@ -268,7 +255,7 @@ int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
     // ask for optimal lwork
     doublecomplex ans;
     zgesvd_ (jobu,jobvt,
-             &m,&n,B,&m,
+             &m,&n,ap,&m,
              sp,
              up,&m,
              vp,&ldvt,
@@ -279,7 +266,7 @@ int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
     CHECK(!work,MEM);
     zgesvd_ (jobu,jobvt,
-             &m,&n,B,&m,
+             &m,&n,ap,&m,
              sp,
              up,&m,
              vp,&ldvt,
@@ -289,7 +276,6 @@ int svd_l_C(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
     CHECK(res,res);
     free(work);
     free(rwork);
-    free(B);
     OK
 }
 
@@ -298,8 +284,7 @@ int zgesdd_ (char *jobz, integer *m, integer *n,
     integer *ldu, doublecomplex *vt, integer *ldvt, doublecomplex *work,
     integer *lwork, doublereal *rwork, integer* iwork, integer *info);
 
-int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
-    //printf("entro\n");
+int svd_l_Cdd(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
     integer m = ar;
     integer n = ac;
     integer q = MIN(m,n);
@@ -319,9 +304,6 @@ int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
         }
     }
     DEBUGMSG("svd_l_Cdd");
-    doublecomplex *B = (doublecomplex*)malloc(m*n*sizeof(doublecomplex));
-    CHECK(!B,MEM);
-    memcpy(B,ap,m*n*sizeof(doublecomplex));
     integer* iwk = (integer*) malloc(8*q*sizeof(integer));
     CHECK(!iwk,MEM);
     int lrwk;
@@ -337,18 +319,17 @@ int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
     integer res;
     // ask for optimal lwk
     doublecomplex ans;
-    zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
+    zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,&ans,&lwk,rwk,iwk,&res);
     lwk = ans.r;
     //printf("lwk = %ld\n",lwk);
     doublecomplex * workv = (doublecomplex*)malloc(lwk*sizeof(doublecomplex));
     CHECK(!workv,MEM);
-    zgesdd_ (jobz,&m,&n,B,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);
+    zgesdd_ (jobz,&m,&n,ap,&m,sp,up,&m,vp,&ldvt,workv,&lwk,rwk,iwk,&res);
     //printf("res = %ld\n",res);
     CHECK(res,res);
     free(workv); // printf("freed workv\n");
     free(rwk);   // printf("freed rwk\n");
     free(iwk);   // printf("freed iwk\n");
-    free(B);     // printf("freed B, salgo\n");
     OK
 }
 
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index 124e353..c91cddd 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -102,25 +102,26 @@ foreign import ccall unsafe "svd_l_Cdd" zgesdd :: TSVD C
 
 -- | Full SVD of a real matrix using LAPACK's /dgesvd/.
 svdR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdR = svdAux dgesvd "svdR" . fmat
+svdR = svdAux dgesvd "svdR"
 
 -- | Full SVD of a real matrix using LAPACK's /dgesdd/.
 svdRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-svdRd = svdAux dgesdd "svdRdd" . fmat
+svdRd = svdAux dgesdd "svdRdd"
 
 -- | Full SVD of a complex matrix using LAPACK's /zgesvd/.
 svdC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
-svdC = svdAux zgesvd "svdC" . fmat
+svdC = svdAux zgesvd "svdC"
 
 -- | Full SVD of a complex matrix using LAPACK's /zgesdd/.
 svdCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
-svdCd = svdAux zgesdd "svdCdd" . fmat
+svdCd = svdAux zgesdd "svdCdd"
 
 svdAux f st x = unsafePerformIO $ do
+    a <- copy ColumnMajor x
     u <- createMatrix ColumnMajor r r
     s <- createVector (min r c)
     v <- createMatrix ColumnMajor c c
-    f # x # u # s # v #| st
+    f # a # u # s # v #| st
     return (u,s,v)
   where
     r = rows x
@@ -129,25 +130,26 @@ svdAux f st x = unsafePerformIO $ do
 
 -- | Thin SVD of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'S\'.
 thinSVDR :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-thinSVDR = thinSVDAux dgesvd "thinSVDR" . fmat
+thinSVDR = thinSVDAux dgesvd "thinSVDR"
 
 -- | Thin SVD of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'S\'.
 thinSVDC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
-thinSVDC = thinSVDAux zgesvd "thinSVDC" . fmat
+thinSVDC = thinSVDAux zgesvd "thinSVDC"
 
 -- | Thin SVD of a real matrix, using LAPACK's /dgesdd/ with jobz == \'S\'.
 thinSVDRd :: Matrix Double -> (Matrix Double, Vector Double, Matrix Double)
-thinSVDRd = thinSVDAux dgesdd "thinSVDRdd" . fmat
+thinSVDRd = thinSVDAux dgesdd "thinSVDRdd"
 
 -- | Thin SVD of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'S\'.
 thinSVDCd :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double, Matrix (Complex Double))
-thinSVDCd = thinSVDAux zgesdd "thinSVDCdd" . fmat
+thinSVDCd = thinSVDAux zgesdd "thinSVDCdd"
 
 thinSVDAux f st x = unsafePerformIO $ do
+    a <- copy ColumnMajor x
     u <- createMatrix ColumnMajor r q
     s <- createVector q
     v <- createMatrix ColumnMajor q c
-    f # x # u # s # v #| st
+    f # a # u # s # v #| st
     return (u,s,v)
   where
     r = rows x
@@ -157,23 +159,24 @@ thinSVDAux f st x = unsafePerformIO $ do
 
 -- | Singular values of a real matrix, using LAPACK's /dgesvd/ with jobu == jobvt == \'N\'.
 svR :: Matrix Double -> Vector Double
-svR = svAux dgesvd "svR" . fmat
+svR = svAux dgesvd "svR"
 
 -- | Singular values of a complex matrix, using LAPACK's /zgesvd/ with jobu == jobvt == \'N\'.
 svC :: Matrix (Complex Double) -> Vector Double
-svC = svAux zgesvd "svC" . fmat
+svC = svAux zgesvd "svC"
 
 -- | Singular values of a real matrix, using LAPACK's /dgesdd/ with jobz == \'N\'.
 svRd :: Matrix Double -> Vector Double
-svRd = svAux dgesdd "svRd" . fmat
+svRd = svAux dgesdd "svRd"
 
 -- | Singular values of a complex matrix, using LAPACK's /zgesdd/ with jobz == \'N\'.
 svCd :: Matrix (Complex Double) -> Vector Double
-svCd = svAux zgesdd "svCd" . fmat
+svCd = svAux zgesdd "svCd"
 
 svAux f st x = unsafePerformIO $ do
+    a <- copy ColumnMajor x
     s <- createVector q
-    g # x # s #| st
+    g # a # s #| st
     return s
   where
     r = rows x
@@ -184,16 +187,17 @@ svAux f st x = unsafePerformIO $ do
 
 -- | Singular values and all right singular vectors of a real matrix, using LAPACK's /dgesvd/ with jobu == \'N\' and jobvt == \'A\'.
 rightSVR :: Matrix Double -> (Vector Double, Matrix Double)
-rightSVR = rightSVAux dgesvd "rightSVR" . fmat
+rightSVR = rightSVAux dgesvd "rightSVR"
 
 -- | Singular values and all right singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'N\' and jobvt == \'A\'.
 rightSVC :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
-rightSVC = rightSVAux zgesvd "rightSVC" . fmat
+rightSVC = rightSVAux zgesvd "rightSVC"
 
 rightSVAux f st x = unsafePerformIO $ do
+    a <- copy ColumnMajor x
     s <- createVector q
     v <- createMatrix ColumnMajor c c
-    g # x # s # v #| st
+    g # a # s # v #| st
     return (s,v)
   where
     r = rows x
@@ -204,16 +208,17 @@ rightSVAux f st x = unsafePerformIO $ do
 
 -- | Singular values and all left singular vectors of a real matrix, using LAPACK's /dgesvd/  with jobu == \'A\' and jobvt == \'N\'.
 leftSVR :: Matrix Double -> (Matrix Double, Vector Double)
-leftSVR = leftSVAux dgesvd "leftSVR" . fmat
+leftSVR = leftSVAux dgesvd "leftSVR"
 
 -- | Singular values and all left singular vectors of a complex matrix, using LAPACK's /zgesvd/ with jobu == \'A\' and jobvt == \'N\'.
 leftSVC :: Matrix (Complex Double) -> (Matrix (Complex Double), Vector Double)
-leftSVC = leftSVAux zgesvd "leftSVC" . fmat
+leftSVC = leftSVAux zgesvd "leftSVC"
 
 leftSVAux f st x = unsafePerformIO $ do
+    a <- copy ColumnMajor x
     u <- createMatrix ColumnMajor r r
     s <- createVector q
-    g # x # u # s #| st
+    g # a # u # s #| st
     return (u,s)
   where
     r = rows x