From 982b24442018c14510ce9bcf4d0e402613fcbea2 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sun, 28 Jun 2015 20:41:09 +0200 Subject: pass copied slice (eig) --- packages/base/src/Internal/C/lapack-aux.c | 50 +++++++++---------------------- packages/base/src/Internal/LAPACK.hs | 37 ++++++++++++----------- 2 files changed, 34 insertions(+), 53 deletions(-) (limited to 'packages/base') 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 -- cgit v1.2.3