pass copied slice (eig)

2 files changed, 34 insertions, 53 deletions

diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 8524962..72c44cb 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -354,12 +354,12 @@ int svd_l_Cdd(KOCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) {
 
 //////////////////// general complex eigensystem ////////////
 
-/* Subroutine */ int zgeev_(char *jobvl, char *jobvr, integer *n,
+int zgeev_(char *jobvl, char *jobvr, integer *n,
         doublecomplex *a, integer *lda, doublecomplex *w, doublecomplex *vl,
         integer *ldvl, doublecomplex *vr, integer *ldvr, doublecomplex *work,
         integer *lwork, doublereal *rwork, integer *info);
 
-int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
+int eig_l_C(OCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
     integer n = ar;
     REQUIRES(ac==n && sn==n, BAD_SIZE);
     REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
@@ -367,18 +367,14 @@ int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
     REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
     char jobvr = vp==NULL?'N':'V';
     DEBUGMSG("eig_l_C");
-    doublecomplex *B = (doublecomplex*)malloc(n*n*sizeof(doublecomplex));
-    CHECK(!B,MEM);
-    memcpy(B,ap,n*n*sizeof(doublecomplex));
     double *rwork = (double*) malloc(2*n*sizeof(double));
     CHECK(!rwork,MEM);
     integer lwork = -1;
     integer res;
     // ask for optimal lwork
     doublecomplex ans;
-    //printf("ask zgeev\n");
     zgeev_  (&jobvl,&jobvr,
-             &n,B,&n,
+             &n,ap,&n,
              sp,
              up,&n,
              vp,&n,
@@ -386,12 +382,10 @@ int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
              rwork,
              &res);
     lwork = ceil(ans.r);
-    //printf("ans = %d\n",lwork);
     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
     CHECK(!work,MEM);
-    //printf("zgeev\n");
     zgeev_  (&jobvl,&jobvr,
-             &n,B,&n,
+             &n,ap,&n,
              sp,
              up,&n,
              vp,&n,
@@ -401,7 +395,6 @@ int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
     CHECK(res,res);
     free(work);
     free(rwork);
-    free(B);
     OK
 }
 
@@ -409,12 +402,12 @@ int eig_l_C(KOCMAT(a), OCMAT(u), CVEC(s),OCMAT(v)) {
 
 //////////////////// general real eigensystem ////////////
 
-/* Subroutine */ int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *
+int dgeev_(char *jobvl, char *jobvr, integer *n, doublereal *
         a, integer *lda, doublereal *wr, doublereal *wi, doublereal *vl,
         integer *ldvl, doublereal *vr, integer *ldvr, doublereal *work,
         integer *lwork, integer *info);
 
-int eig_l_R(KODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
+int eig_l_R(ODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
     integer n = ar;
     REQUIRES(ac==n && sn==n, BAD_SIZE);
     REQUIRES(up==NULL || (ur==n && uc==n), BAD_SIZE);
@@ -422,28 +415,22 @@ int eig_l_R(KODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
     REQUIRES(vp==NULL || (vr==n && vc==n), BAD_SIZE);
     char jobvr = vp==NULL?'N':'V';
     DEBUGMSG("eig_l_R");
-    double *B = (double*)malloc(n*n*sizeof(double));
-    CHECK(!B,MEM);
-    memcpy(B,ap,n*n*sizeof(double));
     integer lwork = -1;
     integer res;
     // ask for optimal lwork
     double ans;
-    //printf("ask dgeev\n");
     dgeev_  (&jobvl,&jobvr,
-             &n,B,&n,
+             &n,ap,&n,
              (double*)sp, (double*)sp+n,
              up,&n,
              vp,&n,
              &ans, &lwork,
              &res);
     lwork = ceil(ans);
-    //printf("ans = %d\n",lwork);
     double * work = (double*)malloc(lwork*sizeof(double));
     CHECK(!work,MEM);
-    //printf("dgeev\n");
     dgeev_  (&jobvl,&jobvr,
-             &n,B,&n,
+             &n,ap,&n,
              (double*)sp, (double*)sp+n,
              up,&n,
              vp,&n,
@@ -451,37 +438,32 @@ int eig_l_R(KODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
              &res);
     CHECK(res,res);
     free(work);
-    free(B);
     OK
 }
 
 
 //////////////////// symmetric real eigensystem ////////////
 
-/* Subroutine */ int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
+int dsyev_(char *jobz, char *uplo, integer *n, doublereal *a,
          integer *lda, doublereal *w, doublereal *work, integer *lwork,
         integer *info);
 
-int eig_l_S(int wantV,KODMAT(a),DVEC(s),ODMAT(v)) {
-    integer n = ar;
-    REQUIRES(ac==n && sn==n, BAD_SIZE);
+int eig_l_S(int wantV,DVEC(s),ODMAT(v)) {
+    integer n = sn;
     REQUIRES(vr==n && vc==n, BAD_SIZE);
     char jobz = wantV?'V':'N';
     DEBUGMSG("eig_l_S");
-    memcpy(vp,ap,n*n*sizeof(double));
     integer lwork = -1;
     char uplo = 'U';
     integer res;
     // ask for optimal lwork
     double ans;
-    //printf("ask dsyev\n");
     dsyev_  (&jobz,&uplo,
              &n,vp,&n,
              sp,
              &ans, &lwork,
              &res);
     lwork = ceil(ans);
-    //printf("ans = %d\n",lwork);
     double * work = (double*)malloc(lwork*sizeof(double));
     CHECK(!work,MEM);
     dsyev_  (&jobz,&uplo,
@@ -496,17 +478,15 @@ int eig_l_S(int wantV,KODMAT(a),DVEC(s),ODMAT(v)) {
 
 //////////////////// hermitian complex eigensystem ////////////
 
-/* Subroutine */ int zheev_(char *jobz, char *uplo, integer *n, doublecomplex
+int zheev_(char *jobz, char *uplo, integer *n, doublecomplex
         *a, integer *lda, doublereal *w, doublecomplex *work, integer *lwork,
         doublereal *rwork, integer *info);
 
-int eig_l_H(int wantV,KOCMAT(a),DVEC(s),OCMAT(v)) {
-    integer n = ar;
-    REQUIRES(ac==n && sn==n, BAD_SIZE);
+int eig_l_H(int wantV,DVEC(s),OCMAT(v)) {
+    integer n = sn;
     REQUIRES(vr==n && vc==n, BAD_SIZE);
     char jobz = wantV?'V':'N';
     DEBUGMSG("eig_l_H");
-    memcpy(vp,ap,2*n*n*sizeof(double));
     double *rwork = (double*) malloc((3*n-2)*sizeof(double));
     CHECK(!rwork,MEM);
     integer lwork = -1;
@@ -514,7 +494,6 @@ int eig_l_H(int wantV,KOCMAT(a),DVEC(s),OCMAT(v)) {
     integer res;
     // ask for optimal lwork
     doublecomplex ans;
-    //printf("ask zheev\n");
     zheev_  (&jobz,&uplo,
              &n,vp,&n,
              sp,
@@ -522,7 +501,6 @@ int eig_l_H(int wantV,KOCMAT(a),DVEC(s),OCMAT(v)) {
              rwork,
              &res);
     lwork = ceil(ans.r);
-    //printf("ans = %d\n",lwork);
     doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
     CHECK(!work,MEM);
     zheev_  (&jobz,&uplo,
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index 65deceb..049c11e 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -225,13 +225,14 @@ leftSVAux f st x = unsafePerformIO $ do
 
 foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok
 foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok
-foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R ::> R :> R ::> Ok
-foreign import ccall unsafe "eig_l_H" zheev :: CInt -> C ::> R :> C ::> Ok
+foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok
+foreign import ccall unsafe "eig_l_H" zheev :: CInt -> R :> C ::> Ok
 
 eigAux f st m = unsafePerformIO $ do
+    a <- copy ColumnMajor m
     l <- createVector r
     v <- createMatrix ColumnMajor r r
-    g # m # l # v #| st
+    g # a # l # v #| st
     return (l,v)
   where
     r = rows m
@@ -241,11 +242,12 @@ eigAux f st m = unsafePerformIO $ do
 -- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/.
 -- The eigenvectors are the columns of v. The eigenvalues are not sorted.
 eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double))
-eigC = eigAux zgeev "eigC" . fmat
+eigC = eigAux zgeev "eigC"
 
 eigOnlyAux f st m = unsafePerformIO $ do
+    a <- copy ColumnMajor m
     l <- createVector r
-    g # m # l #| st
+    g # a # l #| st
     return l
   where
     r = rows m
@@ -254,22 +256,23 @@ eigOnlyAux f st m = unsafePerformIO $ do
 -- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'.
 -- The eigenvalues are not sorted.
 eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double)
-eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat
+eigOnlyC = eigOnlyAux zgeev "eigOnlyC"
 
 -- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/.
 -- The eigenvectors are the columns of v. The eigenvalues are not sorted.
 eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double))
 eigR m = (s', v'')
-    where (s,v) = eigRaux (fmat m)
+    where (s,v) = eigRaux m
           s' = fixeig1 s
           v' = toRows $ trans v
           v'' = fromColumns $ fixeig (toList s') v'
 
 eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double)
 eigRaux m = unsafePerformIO $ do
+    a <- copy ColumnMajor m
     l <- createVector r
     v <- createMatrix ColumnMajor r r
-    g # m # l # v #| "eigR"
+    g # a # l # v #| "eigR"
     return (l,v)
   where
     r = rows m
@@ -289,15 +292,15 @@ fixeig _ _ = error "fixeig with impossible inputs"
 -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
 -- The eigenvalues are not sorted.
 eigOnlyR :: Matrix Double -> Vector (Complex Double)
-eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat
+eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR"
 
 
 -----------------------------------------------------------------------------
 
 eigSHAux f st m = unsafePerformIO $ do
     l <- createVector r
-    v <- createMatrix ColumnMajor r r
-    f # m # l # v #| st
+    v <- copy ColumnMajor m
+    f # l # v #| st
     return (l,v)
   where
     r = rows m
@@ -307,12 +310,12 @@ eigSHAux f st m = unsafePerformIO $ do
 -- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order).
 eigS :: Matrix Double -> (Vector Double, Matrix Double)
 eigS m = (s', fliprl v)
-    where (s,v) = eigS' (fmat m)
+    where (s,v) = eigS' m
           s' = fromList . reverse . toList $  s
 
 -- | 'eigS' in ascending order
 eigS' :: Matrix Double -> (Vector Double, Matrix Double)
-eigS' = eigSHAux (dsyev 1) "eigS'" . fmat
+eigS' = eigSHAux (dsyev 1) "eigS'"
 
 -- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/.
 -- The eigenvectors are the columns of v.
@@ -320,23 +323,23 @@ eigS' = eigSHAux (dsyev 1) "eigS'" . fmat
 eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
 eigH m = (s', fliprl v)
   where
-    (s,v) = eigH' (fmat m)
+    (s,v) = eigH' m
     s' = fromList . reverse . toList $  s
 
 -- | 'eigH' in ascending order
 eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double))
-eigH' = eigSHAux (zheev 1) "eigH'" . fmat
+eigH' = eigSHAux (zheev 1) "eigH'"
 
 
 -- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'.
 -- The eigenvalues are sorted in descending order.
 eigOnlyS :: Matrix Double -> Vector Double
-eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat
+eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'"
 
 -- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'.
 -- The eigenvalues are sorted in descending order.
 eigOnlyH :: Matrix (Complex Double) -> Vector Double
-eigOnlyH = vrev . fst. eigSHAux (zheev 0) "eigH'" . fmat
+eigOnlyH = vrev . fst. eigSHAux (zheev 0) "eigH'"
 
 vrev = flatten . flipud . reshape 1