source: Sophya/trunk/SophyaLib/NTools/matrix.h@ 282

// This may look like C code, but it is really -*- C++ -*-
// Eric Aubourg, Septembre 94
#ifndef MATRIX_SEEN
#define MATRIX_SEEN

#define RGCHECK

#include <stdio.h>
#include <iostream.h>
#include "peida.h"
#include "ppersist.h"

class GeneralFit;

// <summary> Operations matricielles </summary>
// Operations matricielles, algebre lineaire...
class Matrix : public PPersist {
  friend class Vector;
  friend class MatrixRC;
public:
  Matrix();
  // Constructeur, matrice a zero
  Matrix(int r, int c);
  // Constructeur avec des valeurs. Pas d'allocation, juste pointe.
  Matrix(int r, int c, double* values);
  // Constructeur par copie
  Matrix(const Matrix& a);
  // Destructeur
  virtual ~Matrix();

  // Construction automatique pour PPF
  // <group>
  enum {classId = ClassId_Matrix};
  int_4                 ClassId() const         { return classId; }
  static PPersist*      Create() {return new Matrix(1,1);}
  // </group>

  // Remise a zero de tous les elements
  void Zero();
  // Change la taille de la matrice. Reallocation physique seulement si
  // pas assez de place, ou forcee si force=true
  void Realloc(int r, int c, bool force=false);

  // Operateur d'affectation
  Matrix& operator = (const Matrix& a);
  // Affectation depuis scalaire : identite * scalaire
  Matrix& operator = (double x);

  // Impression
  friend ostream& operator << (ostream& s, const Matrix& a);

  // Acces aux elements
  // <group>
  r_8&          operator()(int r, int c);
  r_8 const&    operator()(int r, int c) const;
  // </group>

  // Acces au tableau des elements
  // <group>
  inline r_8* Data() {return data;}
  inline const r_8* Data() const {return data;}
  // </group>
  // Norme 1 : somme des valeurs absolues
  double  Norm1();
  // Norme 2, euclidienne
  double  Norm2();
  
  // Operations matricielles avec scalaire, allocation d'une nouvelle matrice
  // <group>
  friend Matrix operator * (const Matrix& a, double b);
  friend Matrix operator * (double b, const Matrix& a);
  friend Matrix operator + (const Matrix& a, double b);
  friend Matrix operator + (double b, const Matrix& a);
  friend Matrix operator - (const Matrix& a, double b);
  friend Matrix operator - (double b, const Matrix& a);
  friend Matrix operator / (const Matrix& a, double b);
  // </group>

  // Operations matricielles avec scalaires entiers, pour eviter l'appel
  // du constructeur de vecteur a partir d'un entier...
  // sinon v*2 = v * Vector(2), BUG !
  // <group>
  friend Matrix operator * (const Matrix& a, int b);
  friend Matrix operator * (int b, const Matrix& a);
  friend Matrix operator + (const Matrix& a, int b);
  friend Matrix operator + (int b, const Matrix& a);
  friend Matrix operator - (const Matrix& a, int b);
  friend Matrix operator - (int b, const Matrix& a);
  friend Matrix operator / (const Matrix& a, int b);
  // </group>

  // Operations matricielles "en place", avec scalaire
  // <group>
  Matrix& operator *= (double b);
  Matrix& operator += (double b);
  Matrix& operator -= (double b);
  Matrix& operator /= (double b);
  // </group>

  // Operations matricielles
  // <group>
  friend Matrix operator * (const Matrix& a, const Matrix& b);
  friend Matrix operator + (const Matrix& a, const Matrix& b);
  friend Matrix operator - (const Matrix& a, const Matrix& b);

  Matrix& operator *= (const Matrix& a);
  Matrix& operator += (const Matrix& b);
  Matrix& operator -= (const Matrix& b);
  // </group>

  // Matrice transposee
  Matrix Transpose() const;
  // Matrice inverse
  // <thrown> singMatxErr </thrown>
  Matrix Inverse() const;

  // Nombre de lignes
  inline int NRows() const {return nr;}
  // Nombre de colonnes
  inline int NCol()  const {return nc;}

  // Entrees-sorties PPF
  // <group>
  void ReadSelf(PInPersist&);
  void WriteSelf(POutPersist&) const;
  // </group>

  // Acces aux rangees et colonnes
  // <group>
  MatrixRC Row(int r) const;
  MatrixRC Col(int c) const;
  MatrixRC Diag() const;
  // </group>
  /* Versions const ? Alors il faut constructeur const pour MatrixRC,    */
  /* mais MatrixRC contient un Matrix* qui detruit const....             */
  /* mais en eux-meme, ils ne modifient pas la matrice. Ils permettent   */
  /* de le faire...                                                      */

  // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes
  // operations sur la matrice B
  static double GausPiv(Matrix& A, Matrix& B);

  // Residus et fonction fittees.
  Matrix* FitResidus(GeneralFit& gfit);
  Matrix* FitFunction(GeneralFit& gfit);

protected:
  int_4 nr;
  int_4 nc;
  int ndata;
  int nalloc;
  r_8* data;
};


inline r_8& Matrix::operator()(int r, int c)
{
#ifdef RGCHECK
  if (r<0 || r>=nr || c<0 || c>=nc) THROW(rangeCheckErr);
#endif
  return(data[r*nc + c]);
}


inline r_8 const& Matrix::operator()(int r, int c) const
{
#ifdef RGCHECK
  if (r<0 || r>=nr || c<0 || c>=nc) THROW(rangeCheckErr);
#endif
  return(data[r*nc + c]);
}

enum RCKind {matrixRow=0, matrixCol=1, matrixDiag=2};


// <summary> Acces aux lignes et colonnes des matrices </summary>
// Permet d'acceder a une ligne ou colonne d'une matrice, en place, comme
// avec un vecteur.
class MatrixRC EXC_AWARE {
friend class Vector;
friend class Matrix;
public:
  MatrixRC();

  virtual ~MatrixRC() {}

  int Next();
  int Prev();
  int SetCol(int c);
  int SetRow(int r);
  int SetDiag();

  static int Step(const Matrix& m, RCKind rckind);
  static double* Org(const Matrix&, RCKind rckind, int ind=0);

  int NElts() const;
  r_8& operator()(int i);
  r_8  operator()(int i) const;

  MatrixRC& operator = (const MatrixRC& rc);
  //operator Vector() const;
  Vector GetVect() const;

  MatrixRC& operator += (const MatrixRC& rc);
  MatrixRC& operator -= (const MatrixRC& rc);

  MatrixRC& operator *= (double x);
  MatrixRC& operator /= (double x);
  MatrixRC& operator -= (double x);
  MatrixRC& operator += (double x);

  friend double operator * (const MatrixRC& a, const MatrixRC& b);

  MatrixRC& LinComb(double a, double b, const MatrixRC& rc, int first=0);
  MatrixRC& LinComb(double b, const MatrixRC& rc, int first=0);

  int IMaxAbs(int first=0);

  static void Swap(MatrixRC& rc1, MatrixRC& rc2);
  // comme c'est symetrique, soit c'est static MatrixRC, soit static Matrix,
  // soit Matrix avec verification qu'on parle bien de la meme...

protected:
  MatrixRC(Matrix& m, RCKind kind, int index=0);
  Matrix* matrix;
  r_8*    data;
  int     index;
  int     step;
  RCKind  kind;
};


inline int MatrixRC::Step(const Matrix& m, RCKind rckind)
{
  switch (rckind)
    {
    case matrixRow  :
      return 1;
    case matrixCol  :
      return m.nc;
    case matrixDiag :
      return m.nc+1;
    }
  return -1;
}


inline double* MatrixRC::Org(const Matrix& m, RCKind rckind, int index)
{
  switch (rckind)
    {
    case matrixRow :
      return m.data + index * m.nc;
    case matrixCol :
      return m.data + index;
    case matrixDiag :
      return m.data;
    }
  return (double*) -1;
}


inline int MatrixRC::NElts() const
{
  if (!matrix) return -1;             // Failure ?
  switch (kind)
    {
    case matrixRow :
      return matrix->nc;
    case matrixCol :
      return matrix->nr;
    case matrixDiag :
      return matrix->nc;
    }
  return -1;
}



inline  double& MatrixRC::operator()(int i)
{
  return data[i*step];
}


inline  double MatrixRC::operator()(int i) const
{
  return data[i*step];
}

#endif /* MATRIX_SEEN */