GeneralRealMatrix.hh

Go to the documentation of this file.

#ifndef _GENERAL_REAL_MATRIX_H
#define _GENERAL_REAL_MATRIX_H

#ifndef _OCS_SYSTEM_HDRS_
#define _OCS_SYSTEM_HDRS_

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

#include <fstream.h>
#include <complex.h>
#include <stdlib.h>
#include <stdio.h>
#include <strstream>

#include <string>
  using namespace std;


#endif /* _OCS_SYSTEM_HDRS_ */

#include "ocsConst.hh"

// For most compliers Fortran77 subroutines like myfun are recognized
// as myfun_ Therefore:

#define F77NAME(x) x##_

extern "C" void F77NAME(dgetrf) (int*, int*, double*,int*,int*,int*);

extern "C" void F77NAME(dgetrs)(char*,int*,int*,double*,int*,int*,double*,
                                int*,int*);
                
extern "C" void F77NAME(dsyev)(char*,char*,int*,double*,int*,double*,double*,
                               int*,int*);

// The last argument is an int which is the length of the
// character array name.

extern "C" int F77NAME(ilaenv)(int*,char*,char*,int*,int*,int*,int*,int);

class GeneralRealMatrix
{

public:

  GeneralRealMatrix(int NumRows2, int NumCol2);      
  GeneralRealMatrix(const GeneralRealMatrix & Copy);       

  ~GeneralRealMatrix(void);

  void WriteFileMatrix(string OutFileName);


  inline double Get(int i, int j) const
   {
     BoundsCheck(i,j);
     return Data[ (i-1) + (j-1) * NumRows ];
   };


  inline void Set(int i, int j, double Val)
   {
    BoundsCheck(i,j);
    Data[ (i-1) + (j-1) * NumRows ] = Val;
   };


  inline double & Element(int i, int j)
  {
    BoundsCheck(i,j);
    return Data[ (i-1) + (j-1) * NumRows ];
  }


  inline int GetNumRows(void) const { return NumRows; };
  inline int GetNumCols(void) const { return NumCols; };


  void SetAll(double Val);


  void Zero(void);


  void Multiply(const GeneralRealMatrix & A, const GeneralRealMatrix & B);


  void MultiplyTrans(const GeneralRealMatrix & A, const GeneralRealMatrix & B);


  void Image(const GeneralRealMatrix & X, GeneralRealMatrix & Y) const;


  void Symmetrize(void);
 

  void Transpose(GeneralRealMatrix & Atrans) const;


  int Solve( GeneralRealMatrix& X, const GeneralRealMatrix& B) const;


  void SymmetricEig(GeneralRealMatrix & EigValue,
                    GeneralRealMatrix & EigVector) const;
 
  double* Data;

private:  

  inline void BoundsCheck(int i, int j) const
   {
     if( i > NumRows || i < 1)
       {  
         cerr << "GeneralRealMatrix::BoundsCheck: Row index out of bounds."
              << endl
              << "ABORT" << endl;
         exit(1);
       }

     if( j > NumCols || j < 1)
       {  
        cerr << "GeneralRealMatrix::BoundsCheck: Col index out of bounds."
              << endl
              << "ABORT" << endl;
         exit(1);
       }
   }


  int NumRows;
  int NumCols;

};

#endif // ## _GENERAL_REAL_MATRIX_H