Re: Problems with CLAPACK SVD routines on OS X

Robert V. wrote:
On Aug 9, 1:57 pm, "H.S." <hs.saDELETEME...@xxxxxxxxx> wrote:
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

I am trying to incorporate your comments and I am realizing that I am
not exactly sure about all the leading dimension parameters (LDA, LDU,

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):
----------------------------------------------
#include <iostream>
#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;

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

typedef int integer;
typedef double doublereal;
#define SIZE 5



#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 number = 0;
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);
}

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


char JOBU;
char JOBVT;
char 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*N];
doublereal *uu = new doublereal[M*M];
//doublereal *vt = new doublereal[M*N];
doublereal *vt = new doublereal[N*N];

cout << "Initializing a " << M << " x " << N << " matrix." << endl;

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

int row = 0;
int column = 0;

int index = 0;

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, wk,
&LWORK, &INFO);

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

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

cout << endl;

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

return 0;
}
----------------------------------------------


The solution is:
INFO=0
s[ 0 ] = 11.485
s[ 1 ] = 3.26975
s[ 2 ] = 2.65336
s[ 3 ] = 2.08873
s[ 4 ] = 0
s[ 5 ] = 0


Solution given by matlab is:
[U,S,V]=svd([4 1 2 1;3 1 2 3;4 1 3 4;1 3 1 2;2 2 5 1;3
2 2 1])

U =

-0.3805 0.1197 -0.4391 -0.5654 0.0386 -0.5718
-0.4038 0.3451 0.0566 0.2148 0.8057 0.1394
-0.5451 0.4293 -0.0514 0.4321 -0.5732 0.0260
-0.2648 -0.0683 0.8839 -0.2153 -0.0549 -0.3077
-0.4463 -0.8168 -0.1419 0.3213 0.0641 -0.0784
-0.3546 -0.1021 0.0043 -0.5458 -0.1174 0.7430


S =

11.4850 0 0 0
0 3.2698 0 0
0 0 2.6534 0
0 0 0 2.0887
0 0 0 0
0 0 0 0


V =

-0.6212 0.3739 -0.4444 -0.5261
-0.3244 -0.3514 0.7321 -0.4851
-0.5581 -0.6543 -0.2775 0.4283
-0.4443 0.5555 0.4354 0.5518


regards,
->HS




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

->HS



Thanks,
Robert

.

Relevant Pages

  • Re: Help using a C dll function in VB
    ... bool function(char** data, int size); ... it is generally a "taboo" for allocating memory across boundaries - Allocating in the DLL and deallocating in the main process. ... The indirection char **data makes it harder with .NET interop. ... Declare Auto Function initdata Lib "mydll.dll" _ ...
    (microsoft.public.dotnet.languages.vb)
  • 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
    ... int main{ ... wk, &LWORK, &INFO); ... you are allocating the WORK array ...
    (sci.math.num-analysis)