Re: Problems with CLAPACK SVD routines on OS X

So I figured out some of these problems yesterday. After much
wrestling with the dgesvd routine, I tried using the dgesdd routine
(the divide-and-conquer svd function) and after I learned that some of
my parameters were invalid, I was able to get some reasonable results.
I'm still not exactly clear on what they should be. The code below
produces the correct singular values but the u and v matrices are not
right. If you have any thoughts given the new code included below, let
me know.

Thanks again,
Robert

#include </System/Library/Frameworks/vecLib.framework/Versions/A/
Headers/clapack.h>
#include <iostream>
#include "f2c.h"
#include <fstream>
#include <string>
#include <sstream>
#include <stdio.h>
#include <algorithm>

using std::cout;
using std::endl;
using std::ifstream;
using std::istringstream;
using std::cerr;
using std::string;
using std::getline;
using std::ios;
using std::max;
using std::min;

int main(int argc, char **argv) {

if (argc < 2) {
cout << "Error. No input file" << endl;
exit(-1);
}

int number = 0;
ifstream infile(argv[1]);

cout << "Opening " << argv[1] << endl;
if (infile.fail()) {
cerr << "unable to open file " << argv[1] << " for reading" <<
endl;
exit(1);
}

// initialize matrix dimension variables
integer M = 0;
integer N = 0;

// pass through infile to determine dimensions
infile.clear();
infile.seekg(0, ios::beg);
string line = "";

while (getline(infile, line)) {
N = 0;
istringstream s(line);
while(s >> number) {
N++;
}
M++;
}

char JOBU;
char JOBVT;
char JOBZ;

integer LDA = N;
integer LDU = N;
integer LDVT = M;
integer INFO = 0;
integer *IWORK = new integer[8*min(M,N)];

doublereal *s = new doublereal[N];
doublereal *wk = new doublereal[3*M*N];
doublereal *uu = new doublereal[M*M];
doublereal *vt = new doublereal[N*N];

// reset file handle
infile.clear();
infile.seekg(0, ios::beg);

cout << "Initializing a " << M << " x " << N << " matrix." << endl;
doublereal *a = new doublereal[M*N];

cout << "Finished initializing..." << endl;

int row = 0;
int column = 0;

int index = 0;
while (getline(infile, line)) {
column = 0;
istringstream s(line);
while(s >> number && (index < (M*N)) && (row < M) && (column < N))
{
if (number > 0) {
//cout << "A[" << row << "][" << column << "] = " << number <<
endl;
a[index] = number;
} else {
a[index] = 0;
}
index++;
column++;
}
row++;
}
infile.close();

JOBU = 'O';
JOBVT = 'A';
JOBZ = 'O';

cout << "Initializing LWORK..." << endl;

//integer LWORK = max(3*min(M,N)+max(M,N),5*min(M,N));
integer SLWORK = (3 * min(M,N) * min(M,N)) + max( max(M,N), 4 *
min(M,N) * min(M,N) + 4 * min(M,N));

/* 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)
*/

cout << "Calculating SVD..." << endl;

//dgesvd_( &JOBU, &JOBVT, &M, &N, a, &LDA, s, uu, &LDU, vt, &LDVT,
wk, &LWORK, &INFO);
dgesdd_(&JOBZ, &N, &M, a, &LDA, s, uu, &LDU, vt, &LDVT, wk, &SLWORK,
IWORK, &INFO);

cout << "\nINFO=" << INFO << endl;

cout << "u =" << endl;
for (int i=0; i<(M*M); i++) {
if (i%(M)==0) {
printf("\n");
}
printf("%.4f ", uu[i]);
}
cout << endl;

for (int i= 0; i< M; i++ ) {
printf("\ns[ %d ] = %.4f", i, s[i]);
}
cout << endl;

cout << "\nvt =" << endl;
for (int i=0; i<(N*N); i++) {
if (i%(N)==0) {
printf("\n");
}
printf("%.4f ", vt[i]);
}

cout << endl;

// delete all of our dynamically allocated objects
delete [] IWORK;
delete [] s;
delete [] wk;
delete [] uu;
delete [] vt;

return 0;
}

OUTPUT
Opening sample.matrix
Initializing a 4 x 5 matrix.
Finished initializing...
Initializing LWORK...
Calculating SVD...

INFO=0
u =

0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000

s[ 0 ] = 4.0000
s[ 1 ] = 3.0000
s[ 2 ] = 2.2361
s[ 3 ] = -0.0000

vt =

-0.0000 -0.0000 -1.0000 -0.0000 -0.0000
-1.0000 -0.0000 -0.0000 -0.0000 -0.0000
-0.0000 -1.0000 -1.0000 -0.0000 -0.0000
-0.0000 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000 0.0000 0.0000 0.0000


Okay, after going over the documentation of dsevd.f
(http://www.netlib.org/lapack/double/dgesvd.f), I have noticed that some
of your parameters are invalid. For example, you are not allocating
proper space for uu and vt. Also, you are allocating the WORK array
prematurely. You should allocate it after defining LWORK and use LWORK
to allocate memory for WORK -- this reduces the changes of error.

Finally, you should be deleting the allocated memory using "delete []
var;" statement since you are allocating arrays. Also, this should be
done just before the return statement in main().

Please fix the above and see if you still get the segmentation fault. We
will then visit the problem of comparison of your results with MATLAB's.

As a side note, you are using functions from LAPACK as well as from BLAS
libraries.

regards,
->HS


.

Relevant Pages

  • Re: Send message to all client.
    ... int *fdes; ... per connection is a poor design, as is frequently discussed here -- ... that part is ontopic) statically allocating 1000 int's won't even be ... Even if your protocol is designed so no client voluntarily ...
    (comp.programming.threads)
  • [PATCH] ext2: Fix data corruption for racing writes
    ... If two writers allocating blocks to file race with each other (e.g. because ... writepages races with ordinary write or two writepages race with each other), ... int main{ ...
    (Linux-Kernel)
  • Re: declaration of variable in for loop
    ... int i,j; ... block,its allocating in the same address? ... That's what the C language requires. ... for "k" is allocated and deallocated on each iteration of the loop. ...
    (comp.lang.c)
  • Re: Problems with CLAPACK SVD routines on OS X
    ... you are allocating the WORK array ... Let the first line of your input file be two number: ... int main(int argc, char **argv) { ...
    (sci.math.num-analysis)
  • Re: Problems with CLAPACK SVD routines on OS X
    ... you are allocating the WORK array ... You should allocate it after defining LWORK and use LWORK ... Parameter 5 to routine DGESDD was incorrect ... does initialization of the WORK array based on LWORK ...
    (sci.math.num-analysis)