Re: Problems with CLAPACK SVD routines on OS X

On Aug 9, 4:18 pm, "H.S." <hs.saDELETEME...@xxxxxxxxx> wrote:
H.S. wrote:

Ah, that's where your problem is. You should be saving your matrix a
column by column (that is what Fortran routines expect). So a_1,2 should
be followed by a_2,2. I think you are saving it row by row at present.

Once you save your matrices in column by column format, then the number
of rows in your matrix is the leading dimension.

Bye the way, I never checked your input format (the way you read your
input file). You are doing it in a convoluted way. May I suggest another
method:
Let the first line of your input file be two number: Rows Cols
Then you should read in the values of the matrix.

Here is an example input file:
6 4
4 1 2 1
3 1 2 3
4 1 3 4
1 3 1 2
2 2 5 1
3 2 2 1

Here is the code that reads this data file (if opened using an input
stream), it is much more elegant now:
infile >> M; infile >> N;//read in dimens. of matrix
double dTmp;
doublereal *a = new double[M*N];
//now read in the matrix values in
//column by column format for fortan routines
for (r = 0; r < M; r++){
for (c = 0; c < N; c++){
infile >> a[r + c*M];//read in the matrix value
}
}
infile.close();//don't forget to close the file

And, for your convenience, here is the source code that works now
(solution agrees with MATLAB):

The code I posted previously contained some unnecessary statements. Here
is the code again, but much cleaned up:
---------------------------------------------------------------------
//#include <vecLib/vecLib.h>
//#include "clapack.h"

#include <iostream>
#include <fstream>
#include <string>
#include <sstream>
#include <iomanip>

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

//#include "../taucs/external/src/f2c.h"

typedef int integer;
typedef double doublereal;

#ifdef __cplusplus
extern "C" {
#endif
extern 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);
#ifdef __cplusplus}

#endif // __cplusplus

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

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

int r,c;
ifstream infile(argv[1],ios::in);

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

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

cout << "Initializing a " << M << " x " << N << " matrix." << endl;
infile >> M; infile >> N;//read in dimens. of matrix
doublereal *a = new double[M*N];
//now read in the matrix values in
//column by column format for fortan routines
for (r = 0; r < M; r++){
for (c = 0; c < N; c++){
infile >> a[r + c*M];//read in the matrix value
}
}
infile.close();//don't forget to close the file
cout << "Finished initializing..." << endl;

char JOBU, JOBVT, JOBZ;

integer LDA = M;
integer LDU = M;
integer LDVT = N;
integer INFO = 0;

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

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

cout << "Initializing LWORK..." << endl;
integer LWORK = 3*M*N;
doublereal *IWORK = new doublereal[8*LWORK];

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

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

if (INFO != 0){ //svd routine returned an error code
std::cerr << "Error solving SVD. INFO = "
<<". See dgesvd.f documentation" << endl;
}
else {// success in dgesvd
cout << "INFO=" << INFO;

for (int i= 0; i< M; i++ ) {
cout << "s[" << std::setw(4) << i << "] = " << s[i] << endl;
}

cout << endl;
}

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

return 0;}

--------------------------------------------------------------------------



LDVT). When I use M as the leading dimension of A (for LDA), I get the
following run-time error:

Parameter 5 to routine DGESDD was incorrect
Mac OS BLAS parameter error in DGESDD, parameter #0, (unavailable), is
0

Additionally, does initialization of the WORK array based on LWORK
hold for dgesdd as well? I would assume yes.

I wouldn't assume, I would look at dgesdd.f documentation. One site is:http://www.netlib.org/lapack/double/dgesdd.f

regards,
->HS

I had a feeling Fortran's column-major ordering was an issue. Thanks
for all your help.

Best,
Robert

.

Relevant Pages

  • Re: Problems with CLAPACK SVD routines on OS X
    ... column by column (that is what Fortran routines expect). ... I think you are saving it row by row at present. ... Once you save your matrices in column by column format, ... Let the first line of your input file be two number: ...
    (sci.math.num-analysis)
  • Re: Proper use of scanf
    ... have published some routines which operate directly on the input ... Barring EOF it will NOT be a digit. ... This forbids overflow, ... int readxwd ...
    (comp.lang.c)
  • [PATCH 15/20] inflate: (arch) tidy user declarations
    ... bulk of the users are now _not_ a collection of routines copied from ... * Modified for ARM Linux by Russell King ... static int input_data_size; ... - * gzip declarations ...
    (Linux-Kernel)
  • Bulletproof reading of ints (was Re: string concatenation)
    ... this overflow is not caught. ... not happy with the EOF treatment on systems with non-sticky EOF ... return positive int for error, 0 for no error, and EOF for EOF. ... /* These stream input routines are written so that simple ...
    (comp.lang.c)
  • Bulletproof reading of ints (was Re: string concatenation)
    ... this overflow is not caught. ... not happy with the EOF treatment on systems with non-sticky EOF ... return positive int for error, 0 for no error, and EOF for EOF. ... /* These stream input routines are written so that simple ...
    (comp.unix.programmer)