/**************************************************
* unit:    matrix          release 0.6            *
* purpose: general manipulation with matrices     *
***************************************************/
#ifndef __MATRIX__
#define __MATRIX__

#define __ATOMS_MATRIX

#include <math.h>
#include <stdio.h>

#define MatrixId 0x400
#ifndef No_Memory
 #define No_Memory	  0x1
#endif
#define Invalid_Size      0x2
#define Bad_Allocation    0x3
#define Wrong_Pointer     0x4
#define Bad_Index         0x5
#define Incompatible_Size 0x6
#define NotSquare         0x7
#define NoFullRank        0x8


#ifdef Streams
  #include <iostream.h>
#endif
#ifndef __Common_H_INCLUDED
  #include <common.h>
#endif


#ifdef __FixNonTemplateFriend
  #define FIX
#else
  #define FIX <>        //Fix for Non Template Friends
#endif

typedef unsigned TSIZE;         //variable used for indexing inside the matrix
#ifndef EPSILON
 #define EPSILON 1e-12
#endif

template <class FLOAT> class temp_matrix;

template <class FLOAT>
class matrix {
  protected:
      TSIZE size2D,size1D;
      FLOAT **data;

  public:
      matrix(void);
      matrix(FLOAT number);
      matrix(const matrix &);
      matrix(temp_matrix<FLOAT> &);
      matrix(TSIZE x, TSIZE y);
      matrix(TSIZE x, TSIZE y, FLOAT contents);
      matrix(TSIZE x, TSIZE y, FLOAT **DATA);
      ~matrix(void);
                //matrix indexing
      TSIZE rows(void) const    {return(size2D);};
      TSIZE cols(void) const	{return(size1D);};
      FLOAT &operator()(TSIZE x, TSIZE y) const {return member(x,y);}
      FLOAT *operator[](TSIZE y)           {return(data[y]);}
      FLOAT &member(TSIZE x, TSIZE y) const;
      void row(TSIZE y, FLOAT *DataRow);
      FLOAT *row(TSIZE y) const            {return(data[y]);}
      FLOAT **Data(void) const             {return(data);};

      matrix &operator=(const FLOAT & F);
      matrix &operator=(const matrix & m);
      matrix &operator=(temp_matrix<FLOAT> & m);
      friend int operator==(const matrix & m1,const matrix & m2) {return(m1.compare(m2));};
      friend int operator!=(const matrix & m1,const matrix & m2) {return(!(m1.compare(m2)));};

      matrix &operator += (const matrix &);
      matrix &operator += (const FLOAT &);
      matrix &operator -= (const matrix &);
      matrix &operator -= (const FLOAT &);
      matrix &operator *= (const FLOAT &);
      matrix &operator *= (const matrix & m);
      matrix &operator /= (const matrix &);
      matrix &operator /= (const FLOAT &);

      friend temp_matrix<FLOAT> operator- FIX(const matrix & m);
      //friend matrix operator+ (const matrix & m)      {return(m);}

      friend temp_matrix<FLOAT> operator+ FIX(const matrix &, const FLOAT &);
      friend temp_matrix<FLOAT> operator+ (const FLOAT &f, const matrix &m)  {return m + f;};
      friend temp_matrix<FLOAT> operator+ FIX(const matrix &, const matrix &);
      friend temp_matrix<FLOAT> operator- FIX(const matrix &, const FLOAT &);
      friend temp_matrix<FLOAT> operator- FIX(const FLOAT &, const matrix &);
      friend temp_matrix<FLOAT> operator- FIX(const matrix &, const matrix &);
      friend temp_matrix<FLOAT> operator* (const matrix &m, const FLOAT &f)     {return(m.multiply(f));};
      friend temp_matrix<FLOAT> operator* (const FLOAT &f, const matrix &m)     {return(m.multiply(f));};
      friend temp_matrix<FLOAT> operator* (const matrix &m1, const matrix &m2)  {return m1.multiply(m2);};
      friend temp_matrix<FLOAT> operator/ (const matrix &m, const FLOAT &f)     {return(m.divide(f));};
      friend temp_matrix<FLOAT> operator/ FIX(const FLOAT &f, const matrix &m);
      friend temp_matrix<FLOAT> operator/ FIX(const matrix &m1, const matrix &m2);

//    friend matrix minor(matrix & m,int r,int s);

      void fill(FLOAT f);
      int compare(const matrix &m) const;
      temp_matrix<FLOAT> divide(const FLOAT &f) const;
      temp_matrix<FLOAT> multiply(const FLOAT &f) const;
      temp_matrix<FLOAT> multiply(const matrix &m) const;
      void ones(void);
      void resize(TSIZE, TSIZE);
      void zeros(void)          {fill(0);};
      void fliprows(const TSIZE row1,const TSIZE row2)
                   {FLOAT*p=data[row1];data[row1]=data[row2];data[row2]=p;};

      FLOAT det(void) const;
      matrix exp(int MaxIterations=1000) const;
      temp_matrix<FLOAT> inv(void) const;
      FLOAT norm(void) const;
      temp_matrix<FLOAT> triangle(const matrix& M, int *flips=NULL) const;
      temp_matrix<FLOAT> solve(const matrix& M) const;

      friend temp_matrix<FLOAT> transpose FIX( matrix  & m );
      friend temp_matrix<FLOAT> inv(const matrix & m) {return m.inv();};
      friend FLOAT det(const matrix & m) {return m.det();};
      friend matrix exp(const matrix & m) {return m.exp();};
      friend FLOAT norm(const matrix & m) {return m.norm();};

      friend void fill( matrix & m, FLOAT f) {m.fill(f);};
      friend void ones( matrix & m) {m.ones();};
      friend void zeros( matrix & m) {m.fill(0);};
      friend void resize( matrix & m, TSIZE x, TSIZE y) {m.resize(x,y);};

      int check(void);
      void erase(void);
                        //Functions for administration a class matrix
      friend void erase( matrix & m ) {m.erase();};
      friend int check( matrix & m ) {return(m.check());};
#ifdef Streams
//      friend istream& operator>> FIX(istream&, matrix&);
      friend ostream& operator<< FIX(ostream&, const matrix&);
#endif
      };


// Specializations of matrices - matrix types
typedef matrix<int>         int_matrix;
typedef matrix<long int>    long_matrix;
typedef matrix<float>       float_matrix;
typedef matrix<double>      double_matrix;
typedef matrix<long double> long_double_matrix;
#ifdef __ATOMS_INTERVAL
typedef matrix<float_interval>       float_interval_matrix;
typedef matrix<double_interval>      double_interval_matrix;
typedef matrix<long_double_interval> long_double_interval_matrix;
#endif

//Internal speedup class - do not use it!
template <class FLOAT>
class temp_matrix:public matrix<FLOAT>
        {
public: temp_matrix(void)                                       {};
        temp_matrix(FLOAT number):matrix<FLOAT>(number)         {};
        temp_matrix(const matrix<FLOAT> &m):matrix<FLOAT>(m)    {};
        temp_matrix(temp_matrix<FLOAT> &m):matrix<FLOAT>(m)     {};
        temp_matrix(TSIZE y, TSIZE x):matrix<FLOAT>(y,x)        {};
        temp_matrix(TSIZE y, TSIZE x, FLOAT contents):matrix<FLOAT>(y,x,contents) {};
        };

#ifdef FIX
 #undef FIX
#endif

#ifndef __BORLANDC__
template<class FLOAT> temp_matrix<FLOAT> operator - (const matrix<FLOAT> & m);

template<class FLOAT> temp_matrix<FLOAT> operator+ (const matrix<FLOAT> &, const FLOAT &);
template<class FLOAT> temp_matrix<FLOAT> operator+ (const matrix<FLOAT> &, const matrix<FLOAT> &);
template<class FLOAT> temp_matrix<FLOAT> operator- (const matrix<FLOAT> &, const FLOAT &);
template<class FLOAT> temp_matrix<FLOAT> operator- (const FLOAT &, const matrix<FLOAT> &);
template<class FLOAT> temp_matrix<FLOAT> operator- (const matrix<FLOAT> &, const matrix<FLOAT> &);
template<class FLOAT> temp_matrix<FLOAT> operator/ (const FLOAT &, const matrix<FLOAT> &);
template<class FLOAT> temp_matrix<FLOAT> operator/ (const matrix<FLOAT> &, const matrix<FLOAT> &);

template<class FLOAT> temp_matrix<FLOAT> transpose( matrix<FLOAT>  & m );

#ifdef Streams
//template<class FLOAT> istream& operator>>(istream&, matrix<FLOAT>&);
template<class FLOAT> ostream& operator<<(ostream&, const matrix<FLOAT>&);
#endif

#endif


template<class FLOAT> inline FLOAT Abs(FLOAT f) {return ((f>=0)?f:-f);}
inline int Abs(int f)           {return (abs(f));}
inline double Abs(double f)     {return (fabs(f));}


#endif
