From c5ed204b8d6a36681c7ec6b227c634bfae501435 Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Sun, 28 Jun 2015 20:04:02 +0200 Subject: pass copied slice (qr, hess,schur,lu) --- packages/base/src/Internal/C/lapack-aux.c | 136 +++++++++++++----------------- 1 file changed, 58 insertions(+), 78 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 ab78dac..8524962 100644 --- a/packages/base/src/Internal/C/lapack-aux.c +++ b/packages/base/src/Internal/C/lapack-aux.c @@ -901,18 +901,17 @@ int chol_l_S(ODMAT(l)) { //////////////////// QR factorization ///////////////////////// -/* Subroutine */ int dgeqr2_(integer *m, integer *n, doublereal *a, integer * +int dgeqr2_(integer *m, integer *n, doublereal *a, integer * lda, doublereal *tau, doublereal *work, integer *info); -int qr_l_R(KODMAT(a), DVEC(tau), ODMAT(r)) { - integer m = ar; - integer n = ac; +int qr_l_R(DVEC(tau), ODMAT(r)) { + integer m = rr; + integer n = rc; integer mn = MIN(m,n); - REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); + REQUIRES(m>=1 && n >=1 && taun == mn, BAD_SIZE); DEBUGMSG("qr_l_R"); double *WORK = (double*)malloc(n*sizeof(double)); CHECK(!WORK,MEM); - memcpy(rp,ap,m*n*sizeof(double)); integer res; dgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); CHECK(res,res); @@ -920,18 +919,17 @@ int qr_l_R(KODMAT(a), DVEC(tau), ODMAT(r)) { OK } -/* Subroutine */ int zgeqr2_(integer *m, integer *n, doublecomplex *a, +int zgeqr2_(integer *m, integer *n, doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex *work, integer *info); -int qr_l_C(KOCMAT(a), CVEC(tau), OCMAT(r)) { - integer m = ar; - integer n = ac; +int qr_l_C(CVEC(tau), OCMAT(r)) { + integer m = rr; + integer n = rc; integer mn = MIN(m,n); - REQUIRES(m>=1 && n >=1 && rr== m && rc == n && taun == mn, BAD_SIZE); + REQUIRES(m>=1 && n >=1 && taun == mn, BAD_SIZE); DEBUGMSG("qr_l_C"); doublecomplex *WORK = (doublecomplex*)malloc(n*sizeof(doublecomplex)); CHECK(!WORK,MEM); - memcpy(rp,ap,m*n*sizeof(doublecomplex)); integer res; zgeqr2_ (&m,&n,rp,&m,taup,WORK,&res); CHECK(res,res); @@ -939,19 +937,18 @@ int qr_l_C(KOCMAT(a), CVEC(tau), OCMAT(r)) { OK } -/* Subroutine */ int dorgqr_(integer *m, integer *n, integer *k, doublereal * +int dorgqr_(integer *m, integer *n, integer *k, doublereal * a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, integer *info); -int c_dorgqr(KODMAT(a), KDVEC(tau), ODMAT(r)) { - integer m = ar; - integer n = MIN(ac,ar); +int c_dorgqr(KDVEC(tau), ODMAT(r)) { + integer m = rr; + integer n = MIN(rc,rr); integer k = taun; DEBUGMSG("c_dorgqr"); integer lwork = 8*n; // FIXME double *WORK = (double*)malloc(lwork*sizeof(double)); CHECK(!WORK,MEM); - memcpy(rp,ap,m*k*sizeof(double)); integer res; dorgqr_ (&m,&n,&k,rp,&m,(double*)taup,WORK,&lwork,&res); CHECK(res,res); @@ -959,19 +956,18 @@ int c_dorgqr(KODMAT(a), KDVEC(tau), ODMAT(r)) { OK } -/* Subroutine */ int zungqr_(integer *m, integer *n, integer *k, +int zungqr_(integer *m, integer *n, integer *k, doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * work, integer *lwork, integer *info); -int c_zungqr(KOCMAT(a), KCVEC(tau), OCMAT(r)) { - integer m = ar; - integer n = MIN(ac,ar); +int c_zungqr(KCVEC(tau), OCMAT(r)) { + integer m = rr; + integer n = MIN(rc,rr); integer k = taun; DEBUGMSG("z_ungqr"); integer lwork = 8*n; // FIXME doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); CHECK(!WORK,MEM); - memcpy(rp,ap,m*k*sizeof(doublecomplex)); integer res; zungqr_ (&m,&n,&k,rp,&m,(doublecomplex*)taup,WORK,&lwork,&res); CHECK(res,res); @@ -982,20 +978,19 @@ int c_zungqr(KOCMAT(a), KCVEC(tau), OCMAT(r)) { //////////////////// Hessenberg factorization ///////////////////////// -/* Subroutine */ int dgehrd_(integer *n, integer *ilo, integer *ihi, +int dgehrd_(integer *n, integer *ilo, integer *ihi, doublereal *a, integer *lda, doublereal *tau, doublereal *work, integer *lwork, integer *info); -int hess_l_R(KODMAT(a), DVEC(tau), ODMAT(r)) { - integer m = ar; - integer n = ac; +int hess_l_R(DVEC(tau), ODMAT(r)) { + integer m = rr; + integer n = rc; integer mn = MIN(m,n); - REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE); + REQUIRES(m>=1 && n == m && taun == mn-1, BAD_SIZE); DEBUGMSG("hess_l_R"); - integer lwork = 5*n; // fixme + integer lwork = 5*n; // FIXME double *WORK = (double*)malloc(lwork*sizeof(double)); CHECK(!WORK,MEM); - memcpy(rp,ap,m*n*sizeof(double)); integer res; integer one = 1; dgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); @@ -1005,20 +1000,19 @@ int hess_l_R(KODMAT(a), DVEC(tau), ODMAT(r)) { } -/* Subroutine */ int zgehrd_(integer *n, integer *ilo, integer *ihi, +int zgehrd_(integer *n, integer *ilo, integer *ihi, doublecomplex *a, integer *lda, doublecomplex *tau, doublecomplex * work, integer *lwork, integer *info); -int hess_l_C(KOCMAT(a), CVEC(tau), OCMAT(r)) { - integer m = ar; - integer n = ac; +int hess_l_C(CVEC(tau), OCMAT(r)) { + integer m = rr; + integer n = rc; integer mn = MIN(m,n); - REQUIRES(m>=1 && n == m && rr== m && rc == n && taun == mn-1, BAD_SIZE); + REQUIRES(m>=1 && n == m && taun == mn-1, BAD_SIZE); DEBUGMSG("hess_l_C"); - integer lwork = 5*n; // fixme + integer lwork = 5*n; // FIXME doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); CHECK(!WORK,MEM); - memcpy(rp,ap,m*n*sizeof(doublecomplex)); integer res; integer one = 1; zgehrd_ (&n,&one,&n,rp,&n,taup,WORK,&lwork,&res); @@ -1029,23 +1023,17 @@ int hess_l_C(KOCMAT(a), CVEC(tau), OCMAT(r)) { //////////////////// Schur factorization ///////////////////////// -/* Subroutine */ int dgees_(char *jobvs, char *sort, L_fp select, integer *n, +int dgees_(char *jobvs, char *sort, L_fp select, integer *n, doublereal *a, integer *lda, integer *sdim, doublereal *wr, doublereal *wi, doublereal *vs, integer *ldvs, doublereal *work, integer *lwork, logical *bwork, integer *info); -int schur_l_R(KODMAT(a), ODMAT(u), ODMAT(s)) { - integer m = ar; - integer n = ac; - REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE); +int schur_l_R(ODMAT(u), ODMAT(s)) { + integer m = sr; + integer n = sc; + REQUIRES(m>=1 && n==m && ur==n && uc==n, BAD_SIZE); DEBUGMSG("schur_l_R"); - //int k; - //printf("---------------------------\n"); - //printf("%p: ",ap); for(k=0;k0) { return NOCONVER; } @@ -1069,18 +1054,17 @@ int schur_l_R(KODMAT(a), ODMAT(u), ODMAT(s)) { } -/* Subroutine */ int zgees_(char *jobvs, char *sort, L_fp select, integer *n, +int zgees_(char *jobvs, char *sort, L_fp select, integer *n, doublecomplex *a, integer *lda, integer *sdim, doublecomplex *w, doublecomplex *vs, integer *ldvs, doublecomplex *work, integer *lwork, doublereal *rwork, logical *bwork, integer *info); -int schur_l_C(KOCMAT(a), OCMAT(u), OCMAT(s)) { - integer m = ar; - integer n = ac; - REQUIRES(m>=1 && n==m && ur==n && uc==n && sr==n && sc==n, BAD_SIZE); +int schur_l_C(OCMAT(u), OCMAT(s)) { + integer m = sr; + integer n = sc; + REQUIRES(m>=1 && n==m && ur==n && uc==n, BAD_SIZE); DEBUGMSG("schur_l_C"); - memcpy(sp,ap,n*n*sizeof(doublecomplex)); - integer lwork = 6*n; // fixme + integer lwork = 6*n; // FIXME doublecomplex *WORK = (doublecomplex*)malloc(lwork*sizeof(doublecomplex)); doublecomplex *W = (doublecomplex*)malloc(n*sizeof(doublecomplex)); // W not really required in this call @@ -1103,21 +1087,20 @@ int schur_l_C(KOCMAT(a), OCMAT(u), OCMAT(s)) { //////////////////// LU factorization ///////////////////////// -/* Subroutine */ int dgetrf_(integer *m, integer *n, doublereal *a, integer * +int dgetrf_(integer *m, integer *n, doublereal *a, integer * lda, integer *ipiv, integer *info); -int lu_l_R(KODMAT(a), DVEC(ipiv), ODMAT(r)) { - integer m = ar; - integer n = ac; +int lu_l_R(DVEC(ipiv), ODMAT(r)) { + integer m = rr; + integer n = rc; integer mn = MIN(m,n); REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE); DEBUGMSG("lu_l_R"); integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); - memcpy(rp,ap,m*n*sizeof(double)); integer res; dgetrf_ (&m,&n,rp,&m,auxipiv,&res); if(res>0) { - res = 0; // fixme + res = 0; // FIXME } CHECK(res,res); int k; @@ -1129,21 +1112,20 @@ int lu_l_R(KODMAT(a), DVEC(ipiv), ODMAT(r)) { } -/* Subroutine */ int zgetrf_(integer *m, integer *n, doublecomplex *a, +int zgetrf_(integer *m, integer *n, doublecomplex *a, integer *lda, integer *ipiv, integer *info); -int lu_l_C(KOCMAT(a), DVEC(ipiv), OCMAT(r)) { - integer m = ar; - integer n = ac; +int lu_l_C(DVEC(ipiv), OCMAT(r)) { + integer m = rr; + integer n = rc; integer mn = MIN(m,n); REQUIRES(m>=1 && n >=1 && ipivn == mn, BAD_SIZE); DEBUGMSG("lu_l_C"); integer* auxipiv = (integer*)malloc(mn*sizeof(integer)); - memcpy(rp,ap,m*n*sizeof(doublecomplex)); integer res; zgetrf_ (&m,&n,rp,&m,auxipiv,&res); if(res>0) { - res = 0; // fixme + res = 0; // FIXME } CHECK(res,res); int k; @@ -1157,11 +1139,11 @@ int lu_l_C(KOCMAT(a), DVEC(ipiv), OCMAT(r)) { //////////////////// LU substitution ///////////////////////// -/* Subroutine */ int dgetrs_(char *trans, integer *n, integer *nrhs, +int dgetrs_(char *trans, integer *n, integer *nrhs, doublereal *a, integer *lda, integer *ipiv, doublereal *b, integer * ldb, integer *info); -int luS_l_R(KODMAT(a), KDVEC(ipiv), KODMAT(b), ODMAT(x)) { +int luS_l_R(KODMAT(a), KDVEC(ipiv), ODMAT(b)) { integer m = ar; integer n = ac; integer mrhs = br; @@ -1174,19 +1156,18 @@ int luS_l_R(KODMAT(a), KDVEC(ipiv), KODMAT(b), ODMAT(x)) { auxipiv[k] = (integer)ipivp[k]; } integer res; - memcpy(xp,bp,mrhs*nrhs*sizeof(double)); - dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&m,auxipiv,xp,&mrhs,&res); + dgetrs_ ("N",&n,&nrhs,(/*no const (!?)*/ double*)ap,&m,auxipiv,bp,&mrhs,&res); CHECK(res,res); free(auxipiv); OK } -/* Subroutine */ int zgetrs_(char *trans, integer *n, integer *nrhs, +int zgetrs_(char *trans, integer *n, integer *nrhs, doublecomplex *a, integer *lda, integer *ipiv, doublecomplex *b, integer *ldb, integer *info); -int luS_l_C(KOCMAT(a), KDVEC(ipiv), KOCMAT(b), OCMAT(x)) { +int luS_l_C(KOCMAT(a), KDVEC(ipiv), OCMAT(b)) { integer m = ar; integer n = ac; integer mrhs = br; @@ -1199,8 +1180,7 @@ int luS_l_C(KOCMAT(a), KDVEC(ipiv), KOCMAT(b), OCMAT(x)) { auxipiv[k] = (integer)ipivp[k]; } integer res; - memcpy(xp,bp,mrhs*nrhs*sizeof(doublecomplex)); - zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,xp,&mrhs,&res); + zgetrs_ ("N",&n,&nrhs,(doublecomplex*)ap,&m,auxipiv,bp,&mrhs,&res); CHECK(res,res); free(auxipiv); OK -- cgit v1.2.3