Hi,

Why isn't "MatComplex" defined in "nr3.h"?

In "nr3.h" I find several type definitions, e.g. "VecDoub" and "MatDoub". I also find "VecComplex", but "MatComplex" is missing.

Despite this, in the ebook (third edition), it says: "Matrices:
MatInt, MatUint, MatChar, MatUchar, MatLlong, MatUllong, MatDoub,MatComplex, and MatBool."

I have the source code "NR_C302", which I interpret as version 3.02.

Anyone know anything about this?

Best regards, perr

...
Anyone know anything about this? I will make an observation (but I wouldn't dream of trying to put words into the Authors' mouths).

Since none of the NR routines use MatComplex, leaving its definition out of nr3.h doesn't seem to be a big shortcoming.

However...

If you want to experiment with it, try putting the following into your source code somewhere:

#include "nr3.h"
typedef const NRmatrix<Complex> MatComplex_I;
typedef NRmatrix<Complex> MatComplex, MatComplex_O, MatComplex_IO;


For example:

#include "../nr3.h"
typedef const NRmatrix<Complex> MatComplex_I;
typedef NRmatrix<Complex> MatComplex, MatComplex_O, MatComplex_IO;

int main()
{
// The "empty" constructor
MatComplex y;

// Make an array that can be used to initialize a MatComplex:
Complex z[] = {
Complex(1.1,.1), Complex(.2,.3), Complex(.4,.5), Complex(.6,.7),
Complex(1.2,.2), Complex(.3,.4), Complex(.5,.6), Complex(.7,.8),
Complex(1.3,.3), Complex(.4,.5), Complex(.6,.7), Complex(.8,.9)
};

// Create an object that is initialized from the array
MatComplex mc (3, 4, z);

// Overloaded assignment operator
y = mc;

int i, j;
cout << "mc.nrows() = " << mc.nrows()
<< ", mc.ncols() = " << mc.ncols() << endl;

//
// Add 42 to each element and print using real() and imag()
// member functions of the std::complex<double> elements
//
cout << "mc[0][0].real() = " << mc[0][0].real()
<< ", mc[0][0].imag() = " << mc[0][0].imag()
<< endl << endl;
for (i = 0; i < mc.nrows(); i++) {
for (j = 0; j < mc.ncols(); j++) {
mc[i][j] += 42.0;
cout << " (" << mc[i][j].real() << "," << mc[i][j].imag()
<< ")";
}
cout << endl;
}
cout << endl;


cout << "y.nrows() = " << y.nrows()
<< ", y.ncols() = " << y.ncols() << endl;

//
// Print the copy of the original matrix,
// using the std::complex overloaded << operator
//
cout << "y[0][0] = " << y[0][0] << endl << endl;
for (i = 0; i < y.nrows(); i++) {
for (j = 0; j < y.ncols(); j++) {
cout << " " << y[i][j];
}
cout << endl;
}

return 0;
}



Output:

mc.nrows() = 3, mc.ncols() = 4
mc[0][0].real() = 1.1, mc[0][0].imag() = 0.1

(43.1,0.1) (42.2,0.3) (42.4,0.5) (42.6,0.7)
(43.2,0.2) (42.3,0.4) (42.5,0.6) (42.7,0.8)
(43.3,0.3) (42.4,0.5) (42.6,0.7) (42.8,0.9)

y.nrows() = 3, y.ncols() = 4
y[0][0] = (1.1,0.1)

(1.1,0.1) (0.2,0.3) (0.4,0.5) (0.6,0.7)
(1.2,0.2) (0.3,0.4) (0.5,0.6) (0.7,0.8)
(1.3,0.3) (0.4,0.5) (0.6,0.7) (0.8,0.9)


Notes:

Arithmetic and I/O operations on the individual elements are defined by the std::complex stuff in the C++ Standard Library.

Matrix operations (arithmetic or I/O) are not defined in nr3.h for any of the matrix or vector typedefs.

Regards,

Dave

First off, thank you very much for your post Dave, it was a big help for me. After using your approach to creating NRmatrix<Complex> objects I successfully modified gaussj.h to invert complex matrixes.

typedef const NRmatrix<Complex> MatComplex_I;
typedef NRmatrix<Complex> MatComplex, MatComplex_O, MatComplex_IO;

void comgaussj(MatComplex_IO &a, MatComplex_IO &b)
{
Int i,icol,irow,j,k,l,ll,n=a.nrows(),m=b.ncols();
double big;
complex<double> dum,pivinv;
VecInt indxc(n),indxr(n),ipiv(n);
for (j=0;j<n;j++) ipiv[j]=0;
for (i=0;i<n;i++) {
big=0.0;
for (j=0;j<n;j++)
if (ipiv[j] != 1)
for (k=0;k<n;k++) {
if (ipiv[k] == 0) {
if (abs(a[j][k]) >= big) {
big=abs(a[j][k]);
irow=j;
icol=k;
}
}
}
++(ipiv[icol]);
if (irow != icol) {
for (l=0;l<n;l++) SWAP(a[irow][l],a[icol][l]);
for (l=0;l<m;l++) SWAP(b[irow][l],b[icol][l]);
}
indxr[i]=irow;
indxc[i]=icol;
if (a[icol][icol] == 0.0) throw("gaussj: Singular Matrix");
pivinv=1.0/a[icol][icol];
a[icol][icol]=1.0;
for (l=0;l<n;l++) a[icol][l] *= pivinv;
for (l=0;l<m;l++) b[icol][l] *= pivinv;
for (ll=0;ll<n;ll++)
if (ll != icol) {
dum=a[ll][icol];
a[ll][icol]=0.0;
for (l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
for (l=0;l<m;l++) b[ll][l] -= b[icol][l]*dum;
}
}
for (l=n-1;l>=0;l--) {
if (indxr[l] != indxc[l])
for (k=0;k<n;k++)
SWAP(a[k][indxr[l]],a[k][indxc[l]]);
}
}

void comgaussj(MatComplex_IO &a)
{
MatComplex_IO b(a.nrows(),0);
comgaussj(a,b);
}



However, I would like to be able to create a NRmatrix object by passing in a vector of complex objects like this.

#include "nr3.h"

typedef const NRmatrix<Complex> MatComplex_I;
typedef NRmatrix<Complex> MatComplex, MatComplex_O, MatComplex_IO;

int num_cav;
vector<complex<double>> M_con(num_cav*num_cav);
MatComplex a (num_cav, num_cav, M_con);

It seems like it should work if I could just overload NRmatrix correctly. Any pointers?

...
However, I would like to be able to create a NRmatrix object by passing in a vector of complex objects like this.

...
MatComplex a (num_cav, num_cav, M_con);...

It seems like it should work if I could just overload NRmatrix correctly. Any pointers?

No need to overload NRmatrix; you can use it as-is!

Here's the thing: It is guaranteed that a C++ std::vector stores its data elements in contiguous memory. (That is, internally, the vector object stores its elements as an array.)

So, just feed your MatComplex constructor a pointer whose value is the address of the first element in your vector.

You can try either of the following:

Create an "empty" vector and use push_back() to populate it. Make sure you give it enough values so that the result has the right size.


int nrows = 3; // Or whatever
int ncols = 4; // Or whatever

vector <complex<double> > M_con;

for (int i = 0; i < nrows*ncols; i++) {
M_con.push_back(whatever_complex_value_you_have_fo r_this_element);
}

// Create a complex matrix object that is initialized from the vector
MatComplex mc(nrows, ncols, &M_con[0]);


Or...

You can create the vector object with the desired size and populate it with assignment statements:


int nrows = 3;
int ncols = 4;

vector <complex<double> > M_con(nrows*ncols);

for (int i = 0; i < nrows*ncols; i++) {
M_con.[i] = (whatever_complex_value_you_have_for_this_element) ;
}

// Create a complex matrix object that is initialized from the vector
MatComplex mc(nrows, ncols, &M_con[0]);



Taa-daa!

Any pointers?Well, yeah! It's all done with pointers. (Or were you not intending a little pun there? Oh, well, I like it anyhow.)



Regards,

Dave

Footnote: NRVector objects also store their elements in arrays (guaranteed contiguous memory), so you could use my second snippet with a VecComplex if that suits your fancy. (NRVectors don't have push_back(), so my first example will work with a std::vector<complex<double> > but not a VecComplex.)

I went with option 2 and it worked perfectly. You're the man Dave, and thanks for the help.