C:\\Documents and Settings\\morgentd\\My Documents\\wikipedia\\dgesddTest\\Debug>dgesddTest.exe
Optimum lwork size = 47412
Minimum lwork size = 82080
Here's my code:
========================================================
// dgesddTest.c
#include
#include
#include "mkl_cblas.h"
#include "mkl_lapack.h"
#include "mkl.h"
#define MIN(a,b) ((a)<(b)? (a):(b))
#define MAX(a,b) ((a)>(b)? (a):(b))
int main(void) {
// set up the problem
int m=10000;
int n=108;
double *a = malloc(m*n*sizeof(double)); // just let array be filled with
// random junk for this test
// set up for return values (job type 'S')
int minmn = MIN(m,n);
double *u = malloc(m*n*sizeof(double));
double *vt = malloc(n*m*sizeof(double));
double *s = malloc(MAX(1,minmn)*sizeof(double));
// set up for dgesdd
int lwork,info,minlwork;
double *work;
double work0[1]={0.0};
int *iwork;
char jobz[1]={'S'};
iwork=calloc(MAX(1,8*minmn),sizeof(int));
// Workspace query
lwork = -1; info=0;
// signature (FORTRAN) of dgesdd is:
// dgesdd(jobz, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, iwork, info)
dgesdd(jobz,&m,&n,a,&m,s,u,&m,vt,&n,work0,&lwork,iwork,&info);
if ( info<0 ) {
printf("SVD error at %d: Parameter %d had an illegal value.\\n",__LINE__,-info);
return 1;
} else if ( info >0 ) {
printf("SVD error at %d: dbdsdc did not converge, updating process failed\\n",__LINE__);
return 1;
}
// Allocate work area
lwork = (int)(work0[0]);
minlwork = 3*minmn*minmn+MAX(minmn,4*minmn*minmn+4*minmn);
printf("Optimum lwork size = %d\\nMinimum lwork size = %d\\n",lwork,minlwork);
lwork = MAX(lwork,minlwork);
work = malloc(lwork*sizeof(double));
info = 0;
// Compute SVD
dgesdd(jobz,&m,&n,a,&m,s,u,&m,vt,&n,work,&lwork,iwork,&info);
if ( info<0 ) {
printf("SVD error at %d: Parameter %d had an illegal value.\\n",__LINE__,-info);
return 1;
} else if ( info >0 ) {
printf("SVD error at %d: dbdsdc did not converge, updating process failed\\n",__LINE__);
return 1;
}
free(iwork);iwork=0;
free(work);work=0;
return 0;
}
Link Copied
[cpp]minlwork = 3*minmn*minmn+MAX(minmn,(4*minmn*minmn+4*minmn));[/cpp]by
[cpp]minlwork = 3*minmn + MAX(maxmn,(4*minmn*minmn + 4*minmn));in your C program.
[/cpp]
[fxfortran]* The dimension of the array WORK. LWORK >= 1. * If JOBZ = 'N', * LWORK >= 3*min(M,N) + max(max(M,N),7*min(M,N)). * If JOBZ = 'O', * LWORK >= 3*min(M,N) + * max(max(M,N),5*min(M,N)*min(M,N)+4*min(M,N)). * If JOBZ = 'S' or 'A' * LWORK >= 3*min(M,N) + * max(max(M,N),4*min(M,N)*min(M,N)+4*min(M,N)). * For good performance, LWORK should generally be larger. * If LWORK = -1 but other input arguments are legal, WORK(1) * returns the optimal LWORK. [/fxfortran]
[bash] /* We need to check the result because sometimes the "optimal" value is actually too small. Change it to the maximum of the minimum and the optimal. */ long work0 = (long) *DDATA(work); int mn = MIN(m,n); int mx = MAX(m,n); switch(jobz){ case 'N': work0 = MAX(work0,3*mn + MAX(mx,6*mn)+500); break; case 'O': work0 = MAX(work0,3*mn*mn + MAX(mx,5*mn*mn+4*mn+500)); break; case 'S': case 'A': work0 = MAX(work0,3*mn*mn + MAX(mx,4*mn*(mn+1))+500); break; } [/bash]
[bash]C:Documents and SettingsmorgentdMy DocumentswikipediadgesddTestDebug>dgesddTest.exe Optimum lwork size = 15190 Minimum lwork size = 20090[/bash]
[fxfortran] program tdgesdd character jobz integer info,m,n,mxmn,mnmn,lwork,ifor,mnthr,lda,ldu,ldvt,iwork(2) Double Precision work(2) !dummy array double precision A(1,1),S(1,1),U(1,1),VT(1,1) !dummy arrays c do m=1,5000 if(mod(m,100).eq.0)write(*,*)' M = ',m do n=1,5000 JOBZ='S' lda=m+15 ldu=m+10 ldvt=n+5 lwork=-1 Call dgesdd( JOBZ, M, N, A, lda, S, U, ldu, + VT,ldvt, WoRK, LWORK, IWORK, INFO ) MXMN = max(M,N) MNMN = min(M,N) MNTHR = INT( MNMN*11.0D0 / 6.0D0 ) IF(MXMN.LT.MNTHR)THEN IFOR = 3*mnmn + max(mxMN,3*mnMN*mnMN+4*mnMN) ELSE IFOR = 4*mnMN*mnMN+7*mnMN ENDIF if(IFOR.gt.WoRK(1).OR.INFO.NE.0) + write(*,10)m,n,WORK(1),IFOR,INFO 10 format(1x,2I5,2x,F8.0,I11,2x,I2) end do end do end program tdgesdd [/fxfortran]
For more complete information about compiler optimizations, see our Optimization Notice.