SVD difficulties

[clusty]

Hey,

I am experiencing some small troubles with dgesvd_ when I am trying to compute the "economical" version of the SVD (JOBU=S,JOBVT=S).
The code works just fine when I allocate the fullsize arrays (NxN and MxM ) but fails when the arrays are allocated just to their needed size (Visual studio complains of a heap corruption).

I was wondering if you could give me some pointers as to what I am doing wrong.
I am also attaching the xcode in question
=========================================
void svdtest()
{
long m=3,n=8;
////////////////
char JOBU[1]={'S'}, JOBVT[1]={'S'};
double *u=new double[m*min(m,n)];
double *v=new double[min(m,n)*n];
double *s=new double[min(m,n)];
double *work;
long info,lwork=100000;
work=new double[lwork];
/////////////////
double *mat=new double[m*n];
for (int i=0;i<m;i++)
for (int j=0;j<n;j++)
mat[i*n+j]=i+j+1;


/////////////////////

dgesvd_(JOBU, JOBVT, &m, &n, mat, &m, s,u, &m,v,&n, work,&lwork, &info);


delete[] u;
delete[] v;
delete[] s;
delete[] work;
///////////////////


delete[] mat;
}

[admin]

You need first to query the routine with lwork=-1 to know the workspace then call it again for the calculations.
Something like this works:

Code: Select all

int m, n, i, j, lwork;
double *u;
double *v;
double *s;
double *mat;
double *work;
char JOBU[1],JOBVT[1];
static int INFO;

m=3;
n=8;

u   = (double*) malloc(m*min(m,n)*sizeof(double));
v   = (double*) malloc(min(m,n)*n*sizeof(double));
s   = (double*) malloc(min(m,n)*sizeof(double));
mat = (double*) malloc(m*n*sizeof(double));

JOBU[0]='S';
JOBVT[0]='S';

for (i=0;i<m;i++)
   for (j=0;j<n;j++)
   {
     mat[i*n+j]=i+j+1;
   }
lwork=-1;
work = (double*) malloc(2*sizeof(double));

 /* Subroutine  int dgesvd_(char *jobu, char *jobvt, integer *m, integer *n,
   doublereal *a, integer *lda, doublereal *s, doublereal *u, integer *
   ldu, doublereal *vt, integer *ldvt, doublereal *work, integer *lwork,
   integer *info); */
printf("Coucou\n");
dgesvd_(JOBU, JOBVT, &m, &n, mat, &m, s,u, &m,v,&n, work,&lwork, &INFO);
printf("info = %d \n",INFO);
printf("lwork size = %f \n",work[0]);
lwork=(int) work[0];
work = (double*) malloc(lwork*sizeof(double));
dgesvd_(JOBU, JOBVT, &m, &n, mat, &m, s,u, &m,v,&n, work,&lwork, &INFO);
printf("info = %d \n",INFO);

free(u);
free(v);
free(s);
free(work);
free(mat);
}

[clusty]

Hey,

Thanks for help. It is fix it halfway.
Now it work just fine when the matrix is "tall", but will not work for "wide" matrices (works for m>=n).
Seems that when I call free on v something bad happens.

Also had to change every int to long in order for it to compile.