source: Sophya/trunk/Poubelle/DPC:FitsIOServer/NTools/matrix.h

// 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 "anydataobj.h"
#include "ppersist.h"

namespace SOPHYA {
  class GeneralFit;
  template <class T> class TMatrix;
}
// <summary> Operations matricielles </summary>
// Operations matricielles, algebre lineaire...
class OMatrix : public PPersist, public AnyDataObj {
  friend class OVector;
  friend class OMatrixRC;
  friend class TMatrix<r_8>;
public:
  OMatrix();
  // Constructeur, matrice a zero
  OMatrix(int r, int c);
  // Constructeur avec des valeurs. Pas d'allocation, juste pointe.
  OMatrix(int r, int c, double* values);
  // Constructeur par copie
  OMatrix(const OMatrix& a);
  // Constructeur par copie a partir d'une TMatrix<r_8>
  OMatrix(const TMatrix<r_8>& a);
  // Destructeur
  virtual ~OMatrix();

  // Construction automatique pour PPF
  // <group>
  enum {classId = ClassId_Matrix};
  int_4                 ClassId() const         { return classId; }
  static PPersist*      Create() {return new OMatrix(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
  OMatrix& operator = (const OMatrix& a);
  // Affectation depuis scalaire : identite * scalaire
  OMatrix& operator = (double x);

  // Impression
  friend ostream& operator << (ostream& s, const OMatrix& 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 OMatrix operator * (const OMatrix& a, double b);
  friend OMatrix operator * (double b, const OMatrix& a);
  friend OMatrix operator + (const OMatrix& a, double b);
  friend OMatrix operator + (double b, const OMatrix& a);
  friend OMatrix operator - (const OMatrix& a, double b);
  friend OMatrix operator - (double b, const OMatrix& a);
  friend OMatrix operator / (const OMatrix& 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 * OVector(2), BUG !
  // <group>
  friend OMatrix operator * (const OMatrix& a, int b);
  friend OMatrix operator * (int b, const OMatrix& a);
  friend OMatrix operator + (const OMatrix& a, int b);
  friend OMatrix operator + (int b, const OMatrix& a);
  friend OMatrix operator - (const OMatrix& a, int b);
  friend OMatrix operator - (int b, const OMatrix& a);
  friend OMatrix operator / (const OMatrix& a, int b);
  // </group>

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

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

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

  // Matrice transposee
  OMatrix Transpose() const;
  // Matrice inverse
  // <thrown> singMatxErr </thrown>
  OMatrix 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>
  OMatrixRC Row(int r) const;
  OMatrixRC Col(int c) const;
  OMatrixRC Diag() const;
  // </group>
  /* Versions const ? Alors il faut constructeur const pour OMatrixRC,    */
  /* mais OMatrixRC contient un OMatrix* 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(OMatrix& A, OMatrix& B);

  // Residus et fonction fittees.
  OMatrix FitResidus(GeneralFit& gfit
            ,double xorg=0.,double yorg=0.,double dx=1.,double dy=1.);
  OMatrix FitFunction(GeneralFit& gfit
            ,double xorg=0.,double yorg=0.,double dx=1.,double dy=1.);

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


inline r_8& OMatrix::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& OMatrix::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 OMatrixRC EXC_AWARE {
friend class OVector;
friend class OMatrix;
public:
  OMatrixRC();

  virtual ~OMatrixRC() {}

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

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

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

  OMatrixRC& operator = (const OMatrixRC& rc);
  //operator OVector() const;
  OVector GetVect() const;

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

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

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

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

  int IMaxAbs(int first=0);

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

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


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


inline double* OMatrixRC::Org(const OMatrix& 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 OMatrixRC::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& OMatrixRC::operator()(int i)
{
  return data[i*step];
}


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


#endif /* MATRIX_SEEN */