Re: Problems with CLAPACK SVD routines on OS X

Hello,
I'm not using LAPACK, but for delete arrays you must use:

delete [] s;
delete [] wk;
etc.

I don't think that the problem is this, but tell me if now works!

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

delete IWORK;
delete s;
delete wk;
delete uu;
delete vt;


On Aug 7, 5:42 pm, "Robert V." <robert.vo...@xxxxxxxxx> wrote:
Hello,

I apologize if these questions have been addressed before, but I am
new to both C++ and LAPACK.

I am trying to write a program to calculate the SVD of a term-document
matrix. Since I have no experience with Fortran, I thought it best to
do this in C/C++. I am doing my development on an Intel Mac running OS
X 10.4.10.

I am able to get results using Lapack's dgesvd routine, but only when
allocating my array statically. The following is based on some sample
code from the netlib site and it works:

####
#include <vecLib/vecLib.h>
#include "f2c.h"
#include <iostream>

#define SIZE 4

using std::cout;
using std::endl;

int main() {
char JOBU;
char JOBVT;

int i;

integer M = SIZE;
integer N = SIZE;

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

double a[SIZE*SIZE] = { 16.0, 5.0, 9.0 , 4.0, 2.0, 11.0, 7.0 , 14.0,
3.0, 10.0, 6.0, 15.0, 13.0, 8.0, 12.0, 1.0};
double s[SIZE];
double wk[201];
double uu[SIZE*SIZE];
double vt[SIZE*SIZE];

JOBU = 'A';

JOBVT = 'A';

LWORK = 201;

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

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

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

for ( i= 0; i< SIZE; i++ ) {
cout << "\ns[ " << i << " ] = " << s[i];
}

cout << endl;

return 0;}

####

The following program is based on the sample above, but it differs in
that I have to allocate my arrays dynamically. Unfortunately, I can't
seem to get it to run without seg faulting during the call to dgesvd:

####
#include <vecLib/vecLib.h>
#include <iostream>
#include "f2c.h"
#include <fstream>
#include <string>
#include <sstream>

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

#define SIZE 4

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 = M;
integer LDU = M;
integer LDVT = N;
integer INFO = 0;
integer *IWORK = new integer[8*M];

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

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

cout << "Initializing a " << M << " x " << N << " matrix." << endl;
doublereal *a = new double[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 = 'A';
JOBVT = 'A';
JOBZ = 'A';

cout << "Initializing LWORK..." << endl;
integer LWORK = 3*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);

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

//cout << "\n INFO=" << INFO;

//for ( i= 0; i< SIZE; i++ ) {
//cout << "\ns[ " << i << " ] = " << s[i];
//}

//cout << endl;

return 0;}

####

g++ svd_test.cpp -I/usr/local/include -framework vecLib -o svd
./svd small.matrix
Opening small.matrix
Initializing a 5 x 15513 matrix
Finished initializing...
Initializing LWORK...
Calculating SVD...
Segmentation fault

I suspect that my problem stems from my confusion about row-major vs
column-major matrices, but this is a shot in the dark. Any insight
would be greatly appreciated.

Thanks!
Robert




.

Relevant Pages

  • Problems with CLAPACK SVD routines on OS X
    ... int main{ ... char JOBVT; ... integer LWORK; ... integer INFO; ...
    (sci.math.num-analysis)
  • Re: Two dimensional array Question
    ... int main ... Don't treat these arrays of char as if they were ... The easiest way to comprehend multidimensional arrays is to think of them ... In your second loop, you point p to the first ...
    (comp.lang.c)
  • Re: Is it good to use char instead of int to save memory?
    ... Is it recommended to use char instead of int for variables that you ... structs/unions). ... arrays when the arrays are large ...
    (comp.lang.c)
  • Re: Doubt
    ... int mainis better. ... void change (char *b) ... This is automatically done in the case of arrays, ...
    (comp.lang.c)
  • Re: C Strings
    ... I thought string literals were arrays until you operated on them, ... int main ... char *p1; ...
    (comp.lang.c)