3x3 complex matrix

[zeviah5]

I am new to C++. I am trying to read a 3x3 matrix from a file then needs to find eigenvalues and eigenvectors manually. My file is not opening. I am not sure of how to get the eigenvalues and eigenvectors.

#include <iostream>
#include <math.h>
#include <complex>
#include <vector>
#include <fstream>
#include <cstdlib>

using namespace std;

const int n = 3;

int main() {

vector<vector<complex<double>>> matrix;

ifstream inputMatrix;
inputMatrix.open("Matrix.dat");

int rows = 0;
int columns = 0;
complex<double> value;
vector<complex<double>> rowvector;
istringstream row;
string line;
do {
getline(inputMatrix, line, '\n');
row.str(line);
do {
row >> value;
rowvector.push_back(value);
columns++;
} while (!row.eof());
matrix.push_back(rowvector);
rowvector.clear();
rows++;
} while (!inputMatrix.eof());
cout << rows << columns << endl;

vector<complex<double>> printrow;
for (int i = 0; i < rows; i++) {
printrow = matrix[i];
for (int j = 0; j < printrow.size(); j++) {
cout << printrow[i] << ' ';
}
cout << endl;
} //*/

system("pause");
return 0;


}

[gunnerfunner]

for file problems always post some/sufficient sample text from the file

[jonnin]

I think there is a cheezy 3x3 ONLY (and a 2x2 only) way to get the e-vals which can then be used to get the evecs. If you want it more generally (any dimensions) it is more trouble. I think these methods use the determinates which are intractable to compute for bigger matrices, but I am a little rusty here.

Ive got code that can do it for a general matrix but it predates the STL so its C-ish. Havent' used it in a long, long time ... at some point we back-hacked it to accept the <complex> type but we didnt clean up much else.

[zeviah5]

Here is the sample text file that I am trying to read

(1,1) (2,0) (0,1)
(2,0) (0,0) (-1,0)
(0,-1)(-1,0) (2,0)

[lastchance]

Are you trying to read a data file with brackets and commas in? It might be better to read each line into a string, change all brackets and commas for spaces and then you will be able to read real and imaginary parts of each complex number directly.

Don't use .eof - you know exactly how many entries there are per line.

Your matrix is Hermitian - look up "Rayleigh quotient iteration" to find its eigenvalues and eigenvectors. For a Hermitian matrix the eigenvalues should be real. EDIT: Whoops; bad mistake on my part - the top-left component ruins that as it has a non-zero imaginary term on the diagonal. You might be stuck with thrashing through an algebraic solution.

[jonnin]

error handling (reading all as a string and dealing with that) is one way, but you can just read a simpler file format,
R I R I R I
R I ...
...

and read in a loop
file >> vec[i].real;
file >>vec[i].imaginary;

[zeviah5]

Thank you guys.

I have removed the brackets and commas, however I am still unable to read from the file.

Jonnin, should I have a for loop to read in each row and column?

I am not sure how to proceed.

[JLBorges]

> I have removed the brackets and commas, however I am still unable to read from the file.

There is no need to remove the parentheses and commas; std::complex<> is input and output streamable.

#include <iostream>
#include <complex>
#include <array>
#include <sstream>
#include <iomanip>

template < std::size_t N, typename T = double >
using matrix_type = std::array< std::array< std::complex<T>, N >, N > ;

template < std::size_t N, typename T >
std::istream& get_matrix( std::istream& stm, matrix_type<N,T>& matrix )
{
    for( auto& row : matrix )
        for( auto& value : row ) stm >> value ;

    if( stm.fail() ) matrix = {} ;
    return stm ;
}

template < std::size_t N, typename T >
std::ostream& put_matrix( std::ostream& stm, const matrix_type<N,T>& matrix )
{
    for( const auto& row : matrix )
    {
        for( const auto& value : row ) stm << value << ' ' ;
        stm << '\n' ;
    }

    return stm ;
}

int main()
{
    const std::size_t N = 3 ;
    matrix_type<N> matrix ;

    {
        std::istringstream test_file( "(1,1) (2,0) (0,1)\n"
                                      "(2,0) (0,0) (-1,0)\n"
                                      "(0,-1)(-1,0) (2,0)\n" ) ;
        if( get_matrix( test_file, matrix ) )
            std::cout << "matrix was read in successfully\n" ;
    }

    std::cout << std::fixed << std::setprecision(1) ;
    put_matrix( std::cout, matrix ) ;
}

matrix was read in successfully

(1.0,1.0) (2.0,0.0) (0.0,1.0) 
(2.0,0.0) (0.0,0.0) (-1.0,0.0) 
(0.0,-1.0) (-1.0,0.0) (2.0,0.0)

http://coliru.stacked-crooked.com/a/b2e5b6910bdf81c5

[lastchance]

There is no need to remove the parentheses and commas

Thanks @JLBorges; I didn't realise that - I should just have tried it.

@zeviah5,
Should your program work for arbitrary NxN complex matrices, or just for 3x3?

The code below gets the eigenvalues by brute force from the characteristic equation
det(A-lambda.I) = 0
which is cubic for 3x3 matrices, and so just about solvable.

If you need something for NxN matrices that would be much much harder and you have to either attack the characteristic equation or use some iterative technique.

#include <iostream>
#include <fstream>
#include <cmath>
#include <complex>
#include <cassert>
using namespace std;


using cmplx = complex<double>;

const int N = 3;


// Prototypes
void eigenvalues3x3( cmplx A[3][3], cmplx lambda[3] );
cmplx det3x3( cmplx A[3][3] );


//======================================================================


int main()
{
   cmplx A[N][N];      // matrix
   cmplx lambda[N];    // eigenvalues

   // Read matrix
   ifstream in( "matrix.dat" );   assert( in );
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) in >> A[i][j];
   }
   in.close();

   cout << "\nComplex matrix is:\n";
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) cout << A[i][j] << '\t';
      cout << '\n';
   }


   // Compute eigenvalues (WARNING: brute-force method for 3x3 matrices ONLY)
   if ( N == 3 )
   {
      eigenvalues3x3( A, lambda );
      cout << "\nEigenvalues are:\n";
      for ( int i = 0; i < N; i++ ) cout << lambda[i] << '\n';
   }
}


//======================================================================


void eigenvalues3x3( cmplx A[3][3], cmplx lambda[3] )
{
   // Set up characteristic equation:   det( A - lambda I ) = 0
   //    as a cubic in lambda:  a.lambda^3 + b.lambda^2 + c.lambda + d = 0
   cmplx a = -1.0;                                     // -1
   cmplx b = A[0][0] + A[1][1] + A[2][2];              // trace
   cmplx c = A[1][2] * A[2][1] - A[1][1] * A[2][2]     // -sum of diagonal minors
           + A[2][0] * A[0][2] - A[2][2] * A[0][0]
           + A[0][1] * A[1][0] - A[0][0] * A[1][1];
   cmplx d = det3x3( A );                              // determinant

   // Solve cubic by Cardano's method (easier in complex numbers!)
   cmplx p = ( b * b - 3.0 * a * c ) / ( 9.0 * a * a );
   cmplx q = ( 9.0 * a * b * c - 27.0 * a * a * d - 2.0 * b * b * b ) / ( 54.0 * a * a * a );
   cmplx delta = q * q - p * p * p;
   cmplx g1 = pow( q + sqrt( delta ), 1.0 / 3.0 );     // warning: exponents and sqrt of complex
   cmplx g2 = pow( q - sqrt( delta ), 1.0 / 3.0 );     
   cmplx offset = -b / ( 3.0 * a );
   cmplx omega = cmplx( -0.5, 0.5 * sqrt( 3.0 ) );     // complex cube root of unity
   cmplx omega2 = omega * omega;

   lambda[0] = g1          + g2          + offset;
   lambda[1] = g1 * omega  + g2 * omega2 + offset;
   lambda[2] = g1 * omega2 + g2 * omega  + offset;
}


//======================================================================


cmplx det3x3( cmplx A[3][3] )                                           // determinant of 3x3 matrix
{
   return A[0][0] * ( A[1][1] * A[2][2] - A[1][2] * A[2][1] )
        - A[0][1] * ( A[1][0] * A[2][2] - A[1][2] * A[2][0] )
        + A[0][2] * ( A[1][0] * A[2][1] - A[1][1] * A[2][0] );
}


//====================================================================== 

Complex matrix is:
(1,1)	(2,0)	(0,1)	
(2,0)	(0,0)	(-1,0)	
(0,-1)	(-1,0)	(2,0)	

Eigenvalues are:
(3.31839,0.360747)
(1.46466,0.290273)
(-1.78305,0.34898)

[zeviah5]

@JLBorgess, thank you so much.

@lastchance, yes I will have to extend to NxN matrix.
There are some portion of the code that I do not understand. Can I PM you?

[lastchance]

Best if you post your requests in the forum, @zeviah5. More input, more pairs of eyes and some checking of answers.

[zeviah5]

@lastchance, you are correct.

In line 9, using cmplx = complex<double>; are you defining the complex matrix?

line 15: eigenvalues is a void function but I don't see where it is been passed by reference to.

Is there some declarations missing for eg line 31:A , cmplx?

[lastchance]

Line 9 - no, I think that complex<double> is quite a long thing to type; I'm just creating an alias to, hopefully, make my program less cluttered. Every time you see cmplx in the rest of the code it is simply shorthand for complex<double>.

Line 15 - a void function is one that doesn't return a value via the name of the function. Any changed values will be returned from it via its arguments - in this case the array lambda[] with the eigenvalues in.

Matrix A[][] is declared on line 24 as type complex<double> - see my first answer above!

[48T7M4Gy]

Needs a check to make sure matrix size selected is compatible with number of complex numbers in text file:

#include <iostream>
#include <complex>
#include <vector>
#include <fstream>

using namespace std;

int main()
{
    int SIZE = 0;
    cout << "How big is the (square) matrix? ";
    cin >> SIZE;
    
    typedef complex<double> Complex;
    vector<vector<Complex>> matrix(SIZE, vector<Complex>(SIZE));
    
    ifstream inputMatrix;
    inputMatrix.open("complex_matrix.txt");
    
    Complex temp;
    
    if (inputMatrix.is_open())
    {
        int row = 0;
        int col = 0;
        
        // READ IN VALUES
        while(row < SIZE)
        {
            while(col < SIZE and inputMatrix >> matrix[row][col] )
            {
                ++col;
            }
            col = 0;
            ++row;
        }
        
        // PRINT OUT MATRIX
        for(auto row = 0; row < SIZE; ++row)
        {
            for(auto col = 0; col < SIZE; ++col)
            {
                cout << matrix[row][col] << '\t';
            }
            cout << '\n';
        }
        
        inputMatrix.close();
    }
    else cout << "Unable to open file";
    
    return 0;
}

[lastchance]

@zeviah5,

If you are allowed to use libraries, then both Eigen and Lapack have routines for computing eigenvalues and eigenvectors of arbitrary-size complex matrices. Writing such routines from scratch doesn't seem like a pleasant introduction to C++ ... !


Since everybody else seems to want to use vector<> or array<> here's another version of my program. It should be able to identify how big the matrix is (N) from the first line of the file, but it doesn't do input checking on the rest of the rows.

It will still only find eigenvalues for 3x3 matrices.

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <cmath>
#include <complex>
#include <vector>
#include <cassert>
using namespace std;


// Shorthand for some common types, 'cos I'm too lazy to write the whole lot out each time!
using cmplx  = complex<double>;        // complex number
using vec    = vector<cmplx>;          // vector
using matrix = vector<vec>;            // matrix (=collection of (row) vectors)


// Prototypes
int howBig( string filename );
void eigenvalues3x3( matrix A, vec &lambda );
cmplx det3x3( matrix A );


//======================================================================


int main()
{
   string filename = "matrix.dat";

   // Determine how large the matrix is and set array sizes accordingly
   int N = howBig( filename );         
   matrix A( N, vec(N) );             
   vec lambda(N);                     


   // Read matrix
   ifstream in( filename );   assert( in );
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) in >> A[i][j];
   }
   in.close();

   // Output check
   cout << "\nComplex matrix with N = " << N << ":\n";
   for ( int i = 0; i < N; i++ )
   {
      for ( int j = 0; j < N; j++ ) cout << A[i][j] << '\t';
      cout << '\n';
   }


   // Compute eigenvalues (WARNING: brute-force method for 3x3 matrices ONLY)
   if ( N == 3 )
   {
      eigenvalues3x3( A, lambda );
      cout << "\nEigenvalues are:\n";
      for ( int i = 0; i < N; i++ ) cout << lambda[i] << '\n';
   }
}


//======================================================================


int howBig( string filename )          // Reads one line to count the elements in it
{
   string s;
   ifstream in( filename );   assert( in );
   getline( in, s );
   in.close();

   stringstream ss( s );               // Creates a stream from this line
   int N = 0;
   cmplx dummy;
   while( ss >> dummy ) N++;           // Increments N for as many values as present

   return N;
}


//======================================================================


void eigenvalues3x3( matrix A, vec &lambda )
{
   // Set up characteristic equation:   det( A - lambda I ) = 0
   //    as a cubic in lambda:  a.lambda^3 + b.lambda^2 + c.lambda + d = 0
   cmplx a = -1.0;                                     // -1
   cmplx b = A[0][0] + A[1][1] + A[2][2];              // trace
   cmplx c = A[1][2] * A[2][1] - A[1][1] * A[2][2]     // -sum of diagonal minors
           + A[2][0] * A[0][2] - A[2][2] * A[0][0]
           + A[0][1] * A[1][0] - A[0][0] * A[1][1];
   cmplx d = det3x3( A );                              // determinant

   // Solve cubic by Cardano's method (easier in complex numbers!)
   cmplx p = ( b * b - 3.0 * a * c ) / ( 9.0 * a * a );
   cmplx q = ( 9.0 * a * b * c - 27.0 * a * a * d - 2.0 * b * b * b ) / ( 54.0 * a * a * a );
   cmplx delta = q * q - p * p * p;
   cmplx g1 = pow( q + sqrt( delta ), 1.0 / 3.0 );     // warning: complex exponents and sqrt
   cmplx g2 = pow( q - sqrt( delta ), 1.0 / 3.0 );     
   cmplx offset = -b / ( 3.0 * a );
   cmplx omega( -0.5, 0.5 * sqrt( 3.0 ) );             // complex cube root of unity
   cmplx omega2 = omega * omega;

   lambda[0] = g1          + g2          + offset;
   lambda[1] = g1 * omega  + g2 * omega2 + offset;
   lambda[2] = g1 * omega2 + g2 * omega  + offset;
}


//======================================================================


cmplx det3x3( matrix A )                                             // determinant of 3x3 matrix (unchecked)
{
   return A[0][0] * ( A[1][1] * A[2][2] - A[1][2] * A[2][1] )
        - A[0][1] * ( A[1][0] * A[2][2] - A[1][2] * A[2][0] )
        + A[0][2] * ( A[1][0] * A[2][1] - A[1][1] * A[2][0] );
}


//====================================================================== 

Complex matrix with N = 3:
(1,1)	(2,0)	(0,1)	
(2,0)	(0,0)	(-1,0)	
(0,-1)	(-1,0)	(2,0)	

Eigenvalues are:
(3.31839,0.360747)
(1.46466,0.290273)
(-1.78305,0.34898)

[jonnin]

I found lapack a bit frustrating to use. To call some of the advanced routines, you find that you need to call 10 other routines to get the matrix into some factored specialized form before you can solve the problem at hand. Its way overkill for this 3x3, but you will want to do it that way if you get to unstable large matrices.

[JLBorges]

Armadillo wraps LAPACK quite nicely, and it does not infect your code (Apache license)
http://arma.sourceforge.net/

[zeviah5]

I thank you guys for this lovely discussion. I am learning so much.

@JLBorges, I am going to try Armadillo. Thank you.

@lastchance, I tried to use LAPACk didn't work I may not have been understanding procedures correctly. Did you run this code in terminal?

[zeviah5]

@lastchance, I tried using xcode but getting error Assertion failed.

[lastchance]

I'm afraid that means that xcode is not finding your file matrix.dat. Please make sure that it is within whichever folder xcode goes looking for it. I don't use xcode (or any other IDE for that matter), so difficult for me to advise. Is there any option for "Current working directory" or the like? Or maybe put it in the same folder as source code or executable.

On a non-unix system I doubt that it makes much difference, but note that my code is looking for a matrix.dat with a lower-case m.