MatComplex

[perr]

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

[davekw7x]

Originally Posted by 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

[SteveV45]

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?

[davekw7x]

Originally Posted by SteveV45

...
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_for_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!

Originally Posted by SteveV45

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.)

[SteveV45]

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