SVD C example for m < n

[Intel_C_Intel]

Hi,

Does anyone have a SVD example in C when m < n? I am using the code below (it works fine for m >= n), but I get incorrect results for S, U, and V when m < n. I know the section where I am copying p/q back to u/v is incorrect, but p, q, and s are not even correct at that point.

Thanks,
Marcus

int DLLEXPORT dna_lapack_dsvd( int block_size, int m, int n, double* a, double* s, double* u, double* v ){
double* e = new double[n-1];
double* tauq = new double[max(1,min(m,n))];
double* taup = new double[max(1,min(m,n))];
int info;
char Q = 'Q';
char P = 'P';
char U = 'U';

int lwork = (m + n)*block_size;
double* work = new double[lwork];
dgebrd(&m, &n, a, &m, s, e, tauq, taup, work, &lwork, &info);

if( info != 0 ){
delete[] work;
delete[] e;
delete[] taup;
delete[] tauq;
return info;
}

int len = m * n;
double* q = new double[len];
double* p = new double[len];
for( int i = 0; i < len; i++ ){
q = a;
p = a;
}

lwork = min(m,n)*block_size;
delete[] work;
work = new double[lwork];

if( m > n ){
dorgbr(&Q, &m, &n, &n, q, &m, tauq, work, &lwork, &info);
} else{
dorgbr(&Q, &m, &m, &n, q, &m, tauq, work, &lwork, &info);
}

if( info != 0 ){
delete[] work;
delete[] e;
delete[] p;
delete[] q;
delete[] taup;
delete[] tauq;
return info;
}

if( m < n ){
dorgbr(&P, &m, &n, &m, p, &m, taup, work, &lwork, &info);
} else{
dorgbr(&P, &n, &n, &m, p, &m, taup, work, &lwork, &info);
}

if( info != 0 ){
delete[] work;
delete[] e;
delete[] p;
delete[] q;
delete[] taup;
delete[] tauq;
return info;
}
lwork = max(1,4*(n-1));
delete[] work;
work = new double[lwork];

int ncc = 0;
int ldc = 1;
double* ee = new double;
for(int i = 0; i < n; i++ ){
ee = e;
}
delete[] e;

dbdsqr(&U, &n, &n, &m, &ncc, s, ee, p, &m, q, &m, NULL, &ldc, work, &info);

//code is wrong for m < n
if( m < n ){
for( int i = 0; i < m*m; i++){
u = q;
}
for( int i = 0; i < n*n; i++){
v = p;
}
} else{
for( int i = 0; i < m*n; i++){
u = q;
}
for( int i = 0; i < n; i++){
for( int j = 0; j < n; j++){
v[i*n+j] = p[j*m+i];
}
}
}
delete[] work;
delete[] ee;
delete[] p;
delete[] q;
delete[] taup;
delete[] tauq;
return info;
}

Message Edited by Marcus-Cuda2 on 08-23-2004 06:19 AM

Message Edited by Marcus-Cuda2 on 08-23-2004 06:19 AM