--- mgcv/src/mat.c.orig Fri Feb 20 18:34:25 2009 +++ mgcv/src/mat.c Tue Dec 29 07:21:44 2009 @@ -30,7 +30,7 @@ { double *work,*p1,*p2,*p; int piv=1; work=(double *)calloc((size_t) *n,sizeof(double)); - F77_NAME(dchdc)(a,n,n,work,pivot,&piv,rank); + F77_CALL(dchdc)(a,n,n,work,pivot,&piv,rank); /* zero stuff below the leading diagonal */ for (p2=a+ *n,p1=a+1;p20.5) lwork++; - work=(double *)calloc((size_t)lwork,sizeof(double)); - /* actual call */ - F77_NAME(dgesvd)(&jobu,&jobvt, r, c, x, &lda, d, u, &ldu, vt,&ldvt, - work, &lwork, &info); - free(work); + SPL_CALL(dgesvd)(jobu, jobvt, *r, *c, x, lda, d, u, ldu, vt, ldvt, &info); } void mgcv_svd_full(double *x,double *vt,double *d,int *r,int *c) @@ -106,21 +96,11 @@ matrix(um[[2]],q,q);er$v */ { const char jobu='O',jobvt='A'; - int lda,ldu,ldvt,lwork; + int lda,ldu,ldvt; int info; - double work1,*work,*u=NULL; + double *u=NULL; ldu=lda= *r;ldvt = *c; - lwork=-1; - /* workspace query */ - F77_NAME(dgesvd)(&jobu,&jobvt, r, c, x, &lda, d, u, &ldu, vt,&ldvt, - &work1, &lwork, &info); - lwork=(int)floor(work1); - if (work1-lwork>0.5) lwork++; - work=(double *)calloc((size_t)lwork,sizeof(double)); - /* actual call */ - F77_NAME(dgesvd)(&jobu,&jobvt, r, c, x, &lda, d, u, &ldu, vt,&ldvt, - work, &lwork, &info); - free(work); + SPL_CALL(dgesvd)(jobu, jobvt, *r, *c, x, lda, d, u, ldu, vt, ldvt, &info); } @@ -133,17 +113,11 @@ Calls LAPACK routine dormtr */ { char trans='N',side='R',uplo='U'; - int nq,lwork=-1,info; - double *work,work1; + int nq,info; if (*left) { side = 'L';nq = *m;} else nq = *n; if (*transpose) trans = 'T'; /* workspace query ... */ - F77_NAME(dormtr)(&side,&uplo,&trans,m,n,S,&nq,tau,B,m,&work1,&lwork,&info); - lwork=(int)floor(work1);if (work1-lwork>0.5) lwork++; - work=(double *)calloc((size_t)lwork,sizeof(double)); - /* actual call ... */ - F77_NAME(dormtr)(&side,&uplo,&trans,m,n,S,&nq,tau,B,m,work,&lwork,&info); - free(work); + SPL_CALL(dormtr)(side, uplo, trans, *m, *n, S, nq, tau, B, *m, &info); } void mgcv_tri_diag(double *S,int *n,double *tau) @@ -164,18 +138,13 @@ The routine calls dsytrd from LAPACK. */ -{ int lwork=-1,info; - double *work,work1,*e,*d; +{ int info; + double *e,*d; char uplo='U'; d = (double *)calloc((size_t)*n,sizeof(double)); e = (double *)calloc((size_t)*n-1,sizeof(double)); - /* work space query... */ - F77_NAME(dsytrd)(&uplo,n,S,n,d,e,tau,&work1,&lwork,&info); - lwork=(int)floor(work1);if (work1-lwork>0.5) lwork++; - work=(double *)calloc((size_t)lwork,sizeof(double)); - /* Actual call... */ - F77_NAME(dsytrd)(&uplo,n,S,n,d,e,tau,work,&lwork,&info); - free(work);free(d);free(e); + SPL_CALL(dsytrd)(uplo, *n, S, *n, d, e, tau, &info); + free(d);free(e); } @@ -231,15 +200,8 @@ um<-.C("mgcv_qr",as.double(X),as.integer(r),as.integer(c),as.integer(pivot),as.double(tau)) qr.R(qr(X));matrix(um[[1]],r,c)[1:c,1:c] */ -{ int info,lwork=-1,*ip; - double work1,*work; - /* workspace query */ - F77_NAME(dgeqp3)(r,c,x,r,pivot,tau,&work1,&lwork,&info); - lwork=(int)floor(work1);if (work1-lwork>0.5) lwork++; - work=(double *)calloc((size_t)lwork,sizeof(double)); - /* actual call */ - F77_NAME(dgeqp3)(r,c,x,r,pivot,tau,work,&lwork,&info); - free(work); +{ int info,*ip; + SPL_CALL(dgeqp3)(*r, *c, x, *r, pivot, tau, &info); /*if (*r<*c) lwork= *r; else lwork= *c;*/ for (ip=pivot;ip < pivot + *c;ip++) (*ip)--; /* ... for 'tis C in which we work and not the 'cursed Fortran... */ @@ -267,19 +229,11 @@ */ { char side='L',trans='N'; - int lda,lwork=-1,info; - double *work,work1; + int lda,info; if (! *left) { side='R';lda = *c;} else lda= *r; if ( *tp) trans='T'; - /* workspace query */ - F77_NAME(dormqr)(&side,&trans,r,c,k,a,&lda,tau,b,r,&work1,&lwork,&info); - lwork=(int)floor(work1);if (work1-lwork>0.5) lwork++; - work=(double *)calloc((size_t)lwork,sizeof(double)); - /* actual call */ - F77_NAME(dormqr)(&side,&trans,r,c,k,a,&lda,tau,b,r,work,&lwork,&info); - free(work); - + SPL_CALL(dormqr)(side, trans, *r, *c, *k, a, lda, tau, b, *r, &info); } @@ -382,31 +336,17 @@ */ { char jobz='V',uplo='U',range='A'; - double work1,*work,dum1=0,abstol=0.0,*Z,*dum2,x,*p; - int lwork = -1,liwork = -1,iwork1,info,*iwork,dumi=0,n_eval=0,*isupZ,i; + double *work,dum1=0,abstol=0.0,*Z,*dum2,x,*p; + int info,dumi=0,n_eval=0,*isupZ,i; if (*get_vectors) jobz='V'; else jobz='N'; if (*use_dsyevd) - { F77_NAME(dsyevd)(&jobz,&uplo,n,A,n,ev,&work1,&lwork,&iwork1,&liwork,&info); - lwork=(int)floor(work1);if (work1-lwork>0.5) lwork++; - work=(double *)calloc((size_t)lwork,sizeof(double)); - liwork = iwork1;iwork= (int *)calloc((size_t)liwork,sizeof(int)); - F77_NAME(dsyevd)(&jobz,&uplo,n,A,n,ev,work,&lwork,iwork,&liwork,&info); - free(work);free(iwork); + { SPL_CALL(dsyevd)(jobz, uplo, *n, A, *n, ev, &info); } else { Z=(double *)calloc((size_t)(*n * *n),sizeof(double)); /* eigen-vector storage */ isupZ=(int *)calloc((size_t)(2 * *n),sizeof(int)); /* eigen-vector support */ - F77_NAME(dsyevr)(&jobz,&range,&uplo, - n,A,n,&dum1,&dum1,&dumi,&dumi, - &abstol,&n_eval,ev, - Z,n,isupZ, &work1,&lwork,&iwork1,&liwork,&info); - lwork=(int)floor(work1);if (work1-lwork>0.5) lwork++; - work=(double *)calloc((size_t)lwork,sizeof(double)); - liwork = iwork1;iwork= (int *)calloc((size_t)liwork,sizeof(int)); - F77_NAME(dsyevr)(&jobz,&range,&uplo, - n,A,n,&dum1,&dum1,&dumi,&dumi, - &abstol,&n_eval,ev, - Z,n,isupZ, work,&lwork,iwork,&liwork,&info); - free(work);free(iwork); + SPL_CALL(dsyevr)(jobz, range, uplo, + *n, A, *n, dum1, dum1, dumi, dumi, + abstol, &n_eval, ev, Z, *n, isupZ, &info); if (*descending) for (i=0;i<*n/2;i++) { /* reverse the eigenvalues */ x = ev[i]; ev[i] = ev[*n-i-1];ev[*n-i-1] = x;