- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content
So I suspect something is amiss. Any one have any ideas?
I've extracted the problem to some sample code which follows. When I run this I get:
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content
See this post in the Lapack Forum. In the MKL documentation, remove the power 2 in the first term of the formula for lwork for the cases jobz = 'S' and 'O' . Correspondingly, replace the line
[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]
Note that you had incorrectly used minmn as the first argument to MAX, but it had no effect with the values for m and n that you used.
The source code at Netlib of dgesdd.f actually computes minlwork for its internal use, and has the correct formula both in the introductory comments and in the body of the code. Here are the comments in that file as to the size of lwork:
[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]
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content
This "documentation" error seems to be somewhat ubuquitous -- I find it not just in the MKL version, but in all the lapack versions I've found. Addtionally, I have found some interesting ways people have dealt with this. For example, in Python's numpy package this routine is used to provide an SVD. After making the workspace query, the souce code says:
[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]You'll notice that they additionlly add an extra 500 to the size without explanation.
Thanks again!
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content
If this were about boiler design rather than software, we'd call that the usual "corrosion allowance"!
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content
[bash]C:Documents and SettingsmorgentdMy DocumentswikipediadgesddTestDebug>dgesddTest.exe Optimum lwork size = 15190 Minimum lwork size = 20090[/bash]
This happens for a variety of other matrix dimensions too. So I'm leaning toward the approach used by the author of numpy. Thioughts?
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content
The calculation of the optimum LWORK is rather complicated and involves several calls to ILAENV with M,N and other arguments.
(new edits:) The formulae for the minimum LWORK needed are incorrect: those in the MKL documentation for DGESDD as well as those in the comments of dgesdd.f at Netlib.
By running a test driver (in Fortran, see below) with the public domain Lapack 3.2 sources as well as with MKL, for the following input parameter ranges,
JOBZ = 'S'; M, N = 1 to 5000; LDA = M and M + 20; LDU = M and M + 10; LDVT = N and N + 5; LWORK = -1
I arrived at the revised formulae for the minimum LWORK needed (lines 20 to 24).
[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] However, since you are making a dry run with LWORK=-1, and you are allocating the work space dynamically in your C code, it should be perfectly safe for you to CALLOC that amount of memory, and completely ignore the "minimum LWORK" estimate given by the formulae.
I do not think that you need to use ad-hoc adjustments such as those that you mentioned w.r.t. the Python adaptation.
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page