hello everyone,

I'm having troubles in obtaining the SVD solution of a matix, does anyone have a working driver program and data set I can use to compare mine with?

Much appreciated!

hi
you can see Numerical Recipes Forum > Obsolete Editions Forum > Methods: Chapters 2, 11, and 18
regards
amin

...driver program and data set...


//
// Driver to test svdcmp function
//
//
// First line of data file is whatever comment you want
//
// Second line has integer values for number of rows
// and number of columns (the rest of the second line
// is thrown away
//
// Matrix values are read from successive lines
//
//
// davekw7x
// April, 2008

// Whatever you need for the nr3 include files
#include "../code/nr3.h"
#include "../code/svd.h"


int main(void)
{
int i, j, k, m, n;
string txt;

char *inname = "svdtest.dat";
ifstream fp(inname);

if (!fp) {
cerr << "*****File " << inname
<< " could not be opened for reading*****" << endl;
throw("Exiting program"); // Or could simply exit(EXIT_FAILURE)
}
cout << "Opened file " << inname << " for reading." << endl;


getline(fp, txt); // First line is any kind of comment
cout << txt << endl; // Echo the first line

fp >> m >> n; // Read rows, columns from second line
getline(fp, txt); // Eat the rest of this line
if (!fp) {
cerr << "*****Could not read rows and columns from second line*****"
<< endl;
throw("Exiting program"); // Or could simply exit(EXIT_FAILURE)
}

cout << "m = " << m << ", n = " << n << endl;
MatDoub a(m, n);
MatDoub uvw(m, n);

// Get matrix coefficients
for (i = 0; i < m && fp; i++) {
for (j = 0; j < n && fp; j++) {
fp >> a[i][j];
if (!fp) {
cerr << "***** Error reading coefficient["
<< i << "][" << j << "] *****" << endl;
throw("Exiting program"); // Or exit(EXIT_FAILURE);
}
}
}
fp.close();

// The constructor performs the decomposition
SVD mySVD(a);
cout << fixed << setprecision(6);
cout << "**********Decomposition Matrices***********" << endl;
cout << "Matrix U" << endl;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << mySVD.u[i][j];
}
cout << endl;
}
cout << "Diagonal of matrix W" << endl;
for (i = 0; i < n; i++) {
cout << setw(12) << mySVD.w[i];
}

cout << endl << "Matrix V-Transpose" << endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << mySVD.v[j][i];
}
cout << endl;
}
cout << endl;
cout << "***Check product against original matrix***" << endl;
cout << "Original Matrix" << endl;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << a[i][j];
}
cout << endl;
}
cout << "Product [U] x [W] x [V-transpose]" << endl;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
uvw[i][j] = 0.0;
for (k = 0; k < n; k++) {
uvw[i][j] += mySVD.u[i][k] * mySVD.w[k] * mySVD.v[j][k];
}
cout << setw(12) << uvw[i][j];
}
cout << endl;
}
cout << endl;

//This gives a small number for an ill-conditioned system.
cout << "Inverse Condition number = "
<< scientific << mySVD.inv_condition() << endl;

return 0;
}


Test input file:

Test file for xsvd program
3 3 Number of rows, columns
0.299000 0.587000 0.114000
-0.168636 -0.331068 0.499704
0.499813 -0.418531 -0.081282


Output:

Opened file svdtest.dat for reading.
Test file for xsvd program
m = 3, n = 3
**********Decomposition Matrices***********
Matrix U
-0.732559 -0.119020 0.670218
0.553771 0.468367 0.688455
0.395848 -0.875481 0.277196
Diagonal of matrix W
0.803700 0.640264 0.458726
Matrix V-Transpose
-0.142554 -0.969294 0.200366
-0.862374 0.220986 0.455495
0.485787 0.107858 0.867397

***Check product against original matrix***
Original Matrix
0.299000 0.587000 0.114000
-0.168636 -0.331068 0.499704
0.499813 -0.418531 -0.081282
Product [U] x [W] x [V-transpose]
0.299000 0.587000 0.114000
-0.168636 -0.331068 0.499704
0.499813 -0.418531 -0.081282

Inverse Condition number = 5.707672e-01


Regards,

Dave

Thanks for the reply davekw7x, i really appreciate your help!

A problem occurs when I implement your example. The code compiles fine, the output is very odd. It seems the SVD mySVD(a) command does not pass the matrix correctly. The rank of the matrix is 0, mySVD.rank(), and mySVD.nullity is 3. Why would the matrix not be passed?

Opened file svdtest.dat for reading.
Test file for xsvd program
m = 3, n = 3
**********Decomposition Matrices***********
Matrix U
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000
Diagonal of matrix W
0.000000 0.000000 0.000000
Matrix V-Transpose
1.000000 0.000000 0.000000
0.000000 1.000000 0.000000
0.000000 0.000000 1.000000

***Check product against original matrix***
Original Matrix
0.299000 0.587000 0.114000
-0.168636 -0.331068 0.499704
0.499813 -0.418531 -0.081282
Product [u] x [W] x [V-transpose]
0.000000 0.000000 0.000000
0.000000 0.000000 0.000000
0.000000 0.000000 0.000000

Inverse Condition number = 0.000000e+000

Cheers,

...
A problem occurs when I implement your example...

Well, since your output shows the correct input matrix, you are obviously reading it OK. (That's always a potential problem area when using other people's input: maybe DOS or UNIX line endings can confuse certain compilers---anyhow, that part is working.)


Just to be sure that I posted the right stuff, here's what I did:

I pasted the code and example from my previous post into new files on my Linux system.
I compiled and executed with GNU g++ version 4.1.2. Results were the same as I posted.

Then I went to my Windows platform.

I pasted the stuff from my previous post there.
I compiled and executed with GNU g++ 3.4.4 (cygwin port). Got the same results as I posted.

I compiled and executed with Borland bcc32.exe (version 5.82 from the free download of TurboC++ builder). Got the same results as I posted.

I compiled and executed with Microsoft cl.exe (version 15.00.20706.01 from the free download of Visual C++ Express 2008).
I got the same results as I posted. (Except the Microsoft printout for the condition number has three digits in the exponent, which is their normal output. All values were the same.)

What compiler are you using? Where did you get the nr3.h and svd.h routines? Mine are from the Numerical Recipes CD, with source file dates 05/10/2007.

Regards,

Dave

Hi there,

I'm using MS Visual C++ 6.0.

The NR code was downloaded from the NR website last week.

How strange. I will try a different compiler....

...
I'm using MS Visual C++ 6.0...

Yoweee.

I just tried it and the results were the same as yours. Totally bad on my Windows XP platform.

I just about never use that old bewhiskered version except for comparison with other people's installations and problems, and I didn't think of testing this particular example. There are lots and lots of bugs in Visual C++ version 6. The compiler version in my installation is 12.00.8804, and I have installed service packs up to and including service pack 5. Still buggy, and will be forever, I'm thinking.

At least later versions have some chance of getting bugs fixed. (And, as I mentioned, there are free downloads for Microsoft and Borland compilers, but I usually use gcc/cygwin for my Windows command-line stuff.)

Anyhow, that's why I go to great pains to let people know what I used when testing the code.

Regards,

Dave

OK, I just installed MSV C++ 2008 Express and I am now seeing the correct result....so many thanks to you! I have learned a lot!

Cheers,

...I have learned a lot!

Me too.

Regards,

Dave

Hi again. I'm gonna be cheeky...sorry. Could you help?

I have tried making your above code into a separate header file with a structure which contains the reading in operations. However, when I try and display an element of the matrix defined in the structure from within the driver i get a crash. Why does this happen? A similar method works fine when using the svd.h header (i.e. can display matrix u from within driver)...help!

Compiler = MS Visual C++ 2008 Express

Header

// header based on code written by davekw7x, April, 2008
//
//http://www.nr.com/forum/showthread.php?t=1437
//

struct simpleIO {
Int m,n;
void displaymat(const MatDoub a);
void read(char* inname);
string txt;
int rows, cols;
int i, j, k ;
MatDoub a;
MatDoub geta();
};


void simpleIO::displaymat(MatDoub a){

m = a.nrows();
n = a.ncols();
cout << "Original Matrix" << endl;

for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << a[i][j];
}
cout << endl;
}
}
void simpleIO::read(char* inname)
{
ifstream fp(inname);

if (!fp) {
cerr << "*****File " << inname
<< " could not be opened for reading*****" << endl;
throw("Exiting program"); // Or could simply exit(EXIT_FAILURE)
}
cout << "Opened file " << inname << " for reading." << endl;

getline(fp, txt); // First line is any kind of comment
cout << txt << endl; // Echo the first line

fp >> m >> n; // Read rows, columns from second line
getline(fp, txt); // Eat the rest of this line
if (!fp) {
cerr << "*****Could not read rows and columns from second line*****"
<< endl;
throw("Exiting program"); // Or could simply exit(EXIT_FAILURE)
}

cout << "m = " << m << ", n = " << n << endl;
MatDoub a(m, n);

// Get matrix coefficients
for (i = 0; i < m && fp; i++) {
for (j = 0; j < n && fp; j++) {
fp >> a[i][j];
if (!fp) {
cerr << "***** Error reading coefficient["
<< i << "][" << j << "] *****" << endl;
throw("Exiting program"); // Or exit(EXIT_FAILURE);
}
}
}
fp.close();

displaymat(a);

}

Driver


#include "nr3.h"
#include "matrixIO.h"

int main(void)
{

char *inname = "svdtest.dat";
simpleIO myIO;

myIO.read(inname);

cout << myIO.a[1][1];

return 0
}

Hi again. I'm gonna be cheeky...sorry. Could you help?

The struct definition has a MatDoub member called a. When you create myIO, the default constructor creates an "empty" matrix (rows = 0, columns = zero, vector is empty). (The constructor definition for NRMatrix is in nr3.h")

In your read() function, you create a brand-new, local matrix with m rows and n columns:

MatDoub a(m, n);

Everything inside that function uses this new matrix. You can put stuff in it and print it out, etc. When the program leaves read() that matrix goes out of scope.

In the main function you try to access myIO.a[1][1], which doesn't exist, since the matrix in myIO is still empty.

The [] operator to access vector and matrix elements doesn't test to see whether they are within the defined range, and the behavior is undefined. A "memory access violation" commonly causes a program to abort, but even if it didn't, you still wouldn't be able to get to the good stuff, since it isn't there.

Bottom line:
change the above line to

a.resize(m,n);


Now, you will be working with and changing the member matrix, and the results will be available for that object.


Try the following main() function with your original version, then make the change that I showed.

int main(void)
{

char *inname = "svdtest.dat";
simpleIO myIO;

myIO.read(inname);

cout << "After my.IO.read(): myIO.a.nrows() = " << myIO.a.nrows()
<< ", myIO.a.ncols() = " << myIO.a.ncols()
<< endl;
for (int i = 0; i < myIO.a.nrows(); i++) {
for (int j = 0; j < myIO.a.ncols(); j++) {
cout << setw(12) << myIO.a[i][j] << " ";
}
cout << endl;
}

return 0;
}


Regards,

Dave

Cheers Dave, if you're ever in Manchester, UK, i'll buy you a pint!