PDA

View Full Version : Ch 11 - Jacobi Struct Question


csallred
08-14-2009, 10:40 AM
I'm starting this with a disclaimer: my degree is in laser control, not numerical simulations so sorry for the dumb questions and plea for help.

The question: is there anything fundamentally wrong with the following code?
//define Floquet Hamiltonian
double HAM[size*size];
double OmegaR = omegaR(ampR, z, chi);
double OmegaB = omegaB(ampB, z, chi);
writeH(size, HAM, OmegaR, OmegaB);
MatDoub ham(size, size, HAM);
//finish defining Floquet Hamiltonian
At this point ham is a tri-diagonal matrix. In comparing to old code I would like to find the eigenvals and eigenvecs using different methods. The Symmeig struc doesn't complain, but the Jacobi struct (code below) does.

Jacobi res(ham);

cout << res.d;//eigenvalues
cout << endl;
cout << res.v;//eigenvectors

When I try to compile and run this it complains saying: "error: no match for 'operator<<' in 'std::cout << res.Jacobi::d'". I've #included everything I think I need, and (don't groan) Im running this in XCode on a mac in the C++ environment.

Thank you in advance for any brilliant ideas.

MPD78
08-14-2009, 11:19 AM
I believe the problem is this;

Your lines should read

cout << res.d << endl; //eigenvalues

cout << res.v << endl; //eigenvectors

Try that and let me know if the problem is cured.

Thanks
Matt

csallred
08-14-2009, 12:50 PM
Matt,

Thanks for weighing in, I really appreciate it.

Unfortunately it still spits out the same error. I'm wondering if there is something wrong with the way that the code is talking to eigen_sym.h b/c when I create a Symeig object (probably wrong word) the text "Symeig" is a different color in my editor to indicate that it is a type of structure, the text "Jacobi" doesn't imitate that. I'm using the code from the NR3.02 dvd. Do you know if there is any trick asscoiated with the type of complier for this?

Thanks

Claire

MPD78
08-14-2009, 02:23 PM
Okay,

Upon looking more deeply, your problem is that you are not passing your arguments into the struct correctly.

Here is a sample program that I wrote to help you out.

The program uses this matrix as the test.

1.0 -2.0 2.0
-2.0 1.0 -2.0
2.0 -2.0 1.0

The resulting eigenvalues are

lambda1 = lambda2 = -1
lambda3 = 5

#include "eigen_sym.h"

int main()
{

string filename;
ifstream matrixtest;

int n,i,j,k,d,pause;

// obtain the name of the input file
cout << "Enter the name of the matrix input file " << endl;
cin >> filename;
cout << endl;
matrixtest.open(filename.c_str());
if(matrixtest.fail())
cout << "Error opening matrix input file\n";
cout << "Enter the size of the matrix " << endl;
cin >> n;
cout << endl;

// obtain the matrix values
MatDoub a(n,n);
for(i=0;i<n;i++) {
for(j=0;j<n;j++) {
matrixtest >> a[i][j];
}
}

// print out the input matrix and vector
cout << fixed << setprecision(3);
cout << "The matrix read in is of size " << n << " x " << n <<" and contains the following elements" << endl << endl;
// print out the matrix that was read into the program
for(i=0;i<n;i++) {
for(j=0;j<n;j++) {
cout << setw(14) << a[i][j];
}
cout << endl;
}
cout << endl;

// instantiate the struct
Jacobi jac(a);
VecDoub yy(n);
yy = jac.d;
Doub nrt = jac.nrot;

cout << endl << "Results for jac.d" << endl << endl;
for(i=0;i<n;i++) {
cout << setw(14) << yy[i];
cout << endl;
}
cout << endl << "The number of rotations is " << nrt << endl;

// this is a dummy input used to keep the output screen visible. (There are much better ways to
// do this but this is quick and does the job)

cin >> pause;

return 0;
}

The output from the program is shown below.

Enter the name of the matrix input file
symmatrix.dat

Enter the size of the matrix
3

The matrix read in is of size 3 x 3 and contains the following elements

1.000 -2.000 2.000
-2.000 1.000 -2.000
2.000 -2.000 1.000


Results for jac.d

5.000
-1.000
-1.000

The number of rotations is 9.000

The eigenvalues are listed in reverse order but they are correct.

For reference the sample matrix was taken from "A First Course In Differential Equations" Eight Edition by Dennis G. Zill. The example is located on page 343.

EDIT: Additional information.

"Symeig" is a different color in my editor to indicate that it is a type of structure, the text "Jacobi" doesn't imitate that. I'm using the code from the NR3.02 dvd. Do you know if there is any trick asscoiated with the type of complier for this?


The word, struct, is associated with the C++ compilier and turns blue in color when typed in. The word Jacobi is simply the name of the struct so it does not change color.

Also, the method of ouput I showed above for the jac.d results will work for a variable of typdef VecDoub. The method to obtain the results from jac.v is different because it is of typedef MatDoub so an additional loop is required.

Here is how you can obtain the results from jac.v

cout << fixed << setprecision(3);
cout << endl << "Results for jac.v" << endl << endl;
for(i=0;i<n;i++) {
for(j=0;j<n;j++) {
cout << setw(14) << xx[i][j];
}
cout << endl;
}

Here is how you define the MatDoub variable

MatDoub xx(n,m);
xx = jac.v;

Here are my results for jac.v

Results for jac.v

0.577 0.789 0.211
-0.577 0.211 0.789
0.577 -0.577 0.577

END EDIT

Hope this helps you out.

Thanks
Matt

Disclaimer - Use my program at you own risk.