//
//  MatrixEig.h
//
//  Qdos2 - Class for calculating eigenvalues/vectors of a symmetric numeric
//  matrix
//
//  Janet Barnes
//  Sept 1999
//

#ifndef QDOS_EIG_H
#define QDOS_EIG_H

#include <algorithm>
#include "Matrix.h"

using namespace std;

namespace QDOS {
	void eigen(double * a, double * v, int n);
// +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ class MatrixEig
// Class for calculating eigenvalues/egigenvectors for a rank 2 square
// array
// -----------------------------------------------------------------------------
	class MatrixEig
	{
		protected:
			size_t _n;                    // no. rows/cols
			ScalarSeq<double> _eigvals;   // _n elt sequence of _rank 0 matrices
			MatrixSeq<double> _eigvecs;   // _n elt sequence of rank 1 matrices
    
		public:
     // --- Constructors
     // Empty constructor
			MatrixEig(){;}
     // Construct eigenvalues/eigenvectors from Matrix
			template<class T>
					MatrixEig(const Matrix<T>& m):_n(m.dims()[0]),\
					_eigvecs(MatrixSeq<double>(m.dims()[0])),\
							_eigvals(ScalarSeq<double>(m.dims()[0]))
							{
								try
								{
									int res;
             // verify that matrix argument is rank 2
									if (m.rank()!=2)
										throw MatrixRankError();
             //  verify that matrix is symmetric
									if (!m.is_symmetric())
									{
										throw MatrixDomainError(ExceptionMsg()<<\
										"Matrix not symmetric");
									}   
									char jobz='V'; // Compute eigenvalues AND eigenvectors
									char uplo='L';
									long int n=_n;             
									Matrix<double> eigvecs(Dim(_n,_n,ROW_MAJOR_ORDER));
									eigvecs.copy(m);
									long int lda=_n;
									long int info;
									Matrix<double> eigvals(Dim(_n,ROW_MAJOR_ORDER));
             // call CLAPACK routine
#ifdef DIAG
									cout << "calling CLAPACK eigenvalue/eigenvector routine" << endl;
#endif
									eigen(eigvecs.ptr(), eigvals.ptr(), n);
             // convert to MatrixSeq
									for(int i=0;i<n;i++)
									{  
										_eigvecs[i]=eigvecs.col(i);
										_eigvals[i]=eigvals.row(i);
									}
                  
								}           
								catch(Exception& e)
								{
									e.msg() << "[MatrixEig<T>]";
									throw e;
								}
							}
     
     // --- Misc functions
							ScalarSeq<double>& eigenvalues() 
							{
								return _eigvals;
							}
     
							MatrixSeq<double>& eigenvectors() 
							{
								return _eigvecs;
							}
     
	};

}  // end namespace QDOS        
#endif
