From 9c05df0cd663bafaf0b69eafee53fce45569dc95 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Mon, 29 Jun 2015 16:37:51 +0200 Subject: pass copied slice in linearSolve --- packages/base/src/Internal/C/lapack-aux.c | 165 ++++++++++-------------------- 1 file changed, 52 insertions(+), 113 deletions(-) (limited to 'packages/base/src/Internal/C') diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c index 1018126..ca60846 100644 --- a/packages/base/src/Internal/C/lapack-aux.c +++ b/packages/base/src/Internal/C/lapack-aux.c @@ -314,22 +314,19 @@ int svd_l_Cdd(OCMAT(a),OCMAT(u), DVEC(s),OCMAT(v)) { } double *rwk = (double*)malloc(lrwk*sizeof(double));; CHECK(!rwk,MEM); - //printf("%s %ld %d\n",jobz,q,lrwk); integer lwk = -1; integer res; // ask for optimal lwk doublecomplex ans; 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,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(workv); + free(rwk); + free(iwk); OK } @@ -498,80 +495,72 @@ int eig_l_H(int wantV,DVEC(s),OCMAT(v)) { //////////////////// general real linear system //////////// -/* Subroutine */ int dgesv_(integer *n, integer *nrhs, doublereal *a, integer +int dgesv_(integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer *ldb, integer *info); -int linearSolveR_l(KODMAT(a),KODMAT(b),ODMAT(x)) { +int linearSolveR_l(ODMAT(a),ODMAT(b)) { integer n = ar; integer nhrs = bc; REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); DEBUGMSG("linearSolveR_l"); - double*AC = (double*)malloc(n*n*sizeof(double)); - memcpy(AC,ap,n*n*sizeof(double)); - memcpy(xp,bp,n*nhrs*sizeof(double)); integer * ipiv = (integer*)malloc(n*sizeof(integer)); integer res; dgesv_ (&n,&nhrs, - AC, &n, + ap, &n, ipiv, - xp, &n, + bp, &n, &res); if(res>0) { return SINGULAR; } CHECK(res,res); free(ipiv); - free(AC); OK } //////////////////// general complex linear system //////////// -/* Subroutine */ int zgesv_(integer *n, integer *nrhs, doublecomplex *a, +int zgesv_(integer *n, integer *nrhs, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer * info); -int linearSolveC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) { +int linearSolveC_l(OCMAT(a),OCMAT(b)) { integer n = ar; integer nhrs = bc; REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); DEBUGMSG("linearSolveC_l"); - doublecomplex*AC = (doublecomplex*)malloc(n*n*sizeof(doublecomplex)); - memcpy(AC,ap,n*n*sizeof(doublecomplex)); - memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); integer * ipiv = (integer*)malloc(n*sizeof(integer)); integer res; zgesv_ (&n,&nhrs, - AC, &n, + ap, &n, ipiv, - xp, &n, + bp, &n, &res); if(res>0) { return SINGULAR; } CHECK(res,res); free(ipiv); - free(AC); OK } //////// symmetric positive definite real linear system using Cholesky //////////// -/* Subroutine */ int dpotrs_(char *uplo, integer *n, integer *nrhs, +int dpotrs_(char *uplo, integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, integer * info); -int cholSolveR_l(KODMAT(a),KODMAT(b),ODMAT(x)) { +int cholSolveR_l(KODMAT(a),ODMAT(b)) { integer n = ar; + integer lda = aXc; integer nhrs = bc; REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); DEBUGMSG("cholSolveR_l"); - memcpy(xp,bp,n*nhrs*sizeof(double)); integer res; dpotrs_ ("U", &n,&nhrs, - (double*)ap, &n, - xp, &n, + (double*)ap, &lda, + bp, &n, &res); CHECK(res,res); OK @@ -579,21 +568,21 @@ int cholSolveR_l(KODMAT(a),KODMAT(b),ODMAT(x)) { //////// Hermitian positive definite real linear system using Cholesky //////////// -/* Subroutine */ int zpotrs_(char *uplo, integer *n, integer *nrhs, +int zpotrs_(char *uplo, integer *n, integer *nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, integer *info); -int cholSolveC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) { +int cholSolveC_l(KOCMAT(a),OCMAT(b)) { integer n = ar; + integer lda = aXc; integer nhrs = bc; REQUIRES(n>=1 && ar==ac && ar==br,BAD_SIZE); DEBUGMSG("cholSolveC_l"); - memcpy(xp,bp,n*nhrs*sizeof(doublecomplex)); integer res; zpotrs_ ("U", &n,&nhrs, - (doublecomplex*)ap, &n, - xp, &n, + (doublecomplex*)ap, &lda, + bp, &n, &res); CHECK(res,res); OK @@ -601,41 +590,30 @@ int cholSolveC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) { //////////////////// least squares real linear system //////////// -/* Subroutine */ int dgels_(char *trans, integer *m, integer *n, integer * +int dgels_(char *trans, integer *m, integer *n, integer * nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal *work, integer *lwork, integer *info); -int linearSolveLSR_l(KODMAT(a),KODMAT(b),ODMAT(x)) { +int linearSolveLSR_l(ODMAT(a),ODMAT(b)) { integer m = ar; integer n = ac; integer nrhs = bc; - integer ldb = xr; - REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); + integer ldb = bXc; + REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE); DEBUGMSG("linearSolveLSR_l"); - double*AC = (double*)malloc(m*n*sizeof(double)); - memcpy(AC,ap,m*n*sizeof(double)); - if (m>=n) { - memcpy(xp,bp,m*nrhs*sizeof(double)); - } else { - int k; - for(k = 0; k0) { @@ -643,47 +621,35 @@ int linearSolveLSR_l(KODMAT(a),KODMAT(b),ODMAT(x)) { } CHECK(res,res); free(work); - free(AC); OK } //////////////////// least squares complex linear system //////////// -/* Subroutine */ int zgels_(char *trans, integer *m, integer *n, integer * +int zgels_(char *trans, integer *m, integer *n, integer * nrhs, doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb, doublecomplex *work, integer *lwork, integer *info); -int linearSolveLSC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) { +int linearSolveLSC_l(OCMAT(a),OCMAT(b)) { integer m = ar; integer n = ac; integer nrhs = bc; - integer ldb = xr; - REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); + integer ldb = bXc; + REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE); DEBUGMSG("linearSolveLSC_l"); - doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); - memcpy(AC,ap,m*n*sizeof(doublecomplex)); - if (m>=n) { - memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); - } else { - int k; - for(k = 0; k0) { @@ -691,52 +657,40 @@ int linearSolveLSC_l(KOCMAT(a),KOCMAT(b),OCMAT(x)) { } CHECK(res,res); free(work); - free(AC); OK } //////////////////// least squares real linear system using SVD //////////// -/* Subroutine */ int dgelss_(integer *m, integer *n, integer *nrhs, +int dgelss_(integer *m, integer *n, integer *nrhs, doublereal *a, integer *lda, doublereal *b, integer *ldb, doublereal * s, doublereal *rcond, integer *rank, doublereal *work, integer *lwork, integer *info); -int linearSolveSVDR_l(double rcond,KODMAT(a),KODMAT(b),ODMAT(x)) { +int linearSolveSVDR_l(double rcond,ODMAT(a),ODMAT(b)) { integer m = ar; integer n = ac; integer nrhs = bc; - integer ldb = xr; - REQUIRES(m>=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); + integer ldb = bXc; + REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE); DEBUGMSG("linearSolveSVDR_l"); - double*AC = (double*)malloc(m*n*sizeof(double)); double*S = (double*)malloc(MIN(m,n)*sizeof(double)); - memcpy(AC,ap,m*n*sizeof(double)); - if (m>=n) { - memcpy(xp,bp,m*nrhs*sizeof(double)); - } else { - int k; - for(k = 0; k=1 && n>=1 && ar==br && xr==MAX(m,n) && xc == bc, BAD_SIZE); + integer ldb = bXc; + REQUIRES(m>=1 && n>=1 && br==MAX(m,n), BAD_SIZE); DEBUGMSG("linearSolveSVDC_l"); - doublecomplex*AC = (doublecomplex*)malloc(m*n*sizeof(doublecomplex)); double*S = (double*)malloc(MIN(m,n)*sizeof(double)); double*RWORK = (double*)malloc(5*MIN(m,n)*sizeof(double)); - memcpy(AC,ap,m*n*sizeof(doublecomplex)); - if (m>=n) { - memcpy(xp,bp,m*nrhs*sizeof(doublecomplex)); - } else { - int k; - for(k = 0; k