/* *************************************************************
   
    CSymMatrix class - Class that implements basic 
	               operations with sqare matrices in a
			userfriendly form. Also it is able to
			calculate eigenvalues and eigenvectors 
			of a matrix.

     Comment: Diagonalization procedure 
	             void eigenproblem(CSymMatrix&, vector<T>&)
		     is an simple adjustment of EISPACK code.
			  
        	     The emphasize was on code transparency. Efficiency was
		     traded off if it would considerably alter readability 
		     of the code

	Version: 0.5 (not optimized for speed and memory)
	Author: Rastko Sknepnek
	Date:   Feb, 2003

    *************************************************************   */



#include <iostream>
#include <vector>
#include <math.h>

using namespace std;

// CSymMatrix - Class for symmetric matrix

template <class T>
class CSymMatrix
{
	private:
		vector<T*> m; //contains matrix
		int size; // cintains matrix size
		CSymMatrix tridiag(CSymMatrix&); // Returns tridiagonal form of matrix
		bool isSymmetric(); // Terurns true if matrix is symmetric
	public:
		CSymMatrix(const int&); // Creates empty matrix
		CSymMatrix(const int, const T& t); // Creates a matrix, and sets diagonal elemets to t; off-diagonal is zero
		CSymMatrix(const int, vector<T *>); //Creates a matrix and initializes it to a given matix
		CSymMatrix(const CSymMatrix&); // Same as above
		~CSymMatrix(); // Destroys matrix

		void print(ostream&); // Prints matrix
		int getsize() const; //Returns matrix size
		T trace(); // Calculates trace
		CSymMatrix transpose(); // returns transposed matrix
		void eigenproblem(CSymMatrix&, vector<T>&); //Returns Eigenvalues and eigen vectors
		
		CSymMatrix& operator= (const CSymMatrix&); // Set value to the matrix
		CSymMatrix operator+ (const CSymMatrix&);  //Add two matrices
		CSymMatrix operator- (const CSymMatrix&);
		CSymMatrix operator/ (const T&);
		CSymMatrix operator+= (const CSymMatrix&);
		CSymMatrix operator-= (const CSymMatrix&);
		CSymMatrix operator/= (const T&);
		CSymMatrix operator* (const CSymMatrix&); // Multiply two matrices
		CSymMatrix& operator*= (const CSymMatrix&);
		vector<T> operator*(vector<T>&); // Matrix-vector multiplication
		T* operator[](const int i) const; // Returns index of a row
};

template<class T>
CSymMatrix<T> operator* (const vector<T>&, const vector<T>&); //outer vector product

template<class T>
T inner (const vector<T>&, const vector<T>&);

template<class T>
vector<T>& operator *= (vector<T>&, const T&);

// Used by EISPACK procedure for diagonalization
template <class T>
T SIGN(T a, T b);

template<class T>
void print_vector(const vector<T>&);

Usage:

CSymMatrix class provides four different constructors. Three of them take as
the frist agrument size of the matrix. Fourth constructor just makes a copy of 
a given matrix.

Eigen value is calulated using method eigneproblem. A matrix sent to eigenproblem as the
first argument will ountain in its columns eigenvectors. The vector in the second argument
will contain eigenvalues. Both, eigenvectors and eigen values are sorted.

Here is a sample code:

#include "matrix.h"
#include <stdlib.h>

const int N = 20;

int main()
{
   // .....

   CSymMatrix<double> A(N); // Creates empty NxN matrix of double precision elemets

   vector<double> EigVal(N); // EigVal will contain eigenvalues
   CSymMatrix<double> EigVec(N); // EigVec will contain eigenvectors

  // Prepare matrix A

  for(int i=0; i<N; i++)
	for(int j=0; j<N; j++)
             A[i][j] = double(rand())/RAND_MAX; // Sets elements of A to random numbers between 0 and 1

  A.eigenproblem(EigVec, EigVal);  // Diagonalizes A
 
  print_vector(cout,EigVal); // Prints eigenvalues

  EigVec.print(cout);  // Prints eigenvectors

  // .....

  return 0;
}