Changeset 926 in Sophya for trunk/SophyaLib/TArray/sopemtx.cc


Ignore:
Timestamp:
Apr 13, 2000, 8:39:39 PM (25 years ago)
Author:
ansari
Message:

documentation cmv 13/4/00

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/SophyaLib/TArray/sopemtx.cc

    r813 r926  
    1414// -------------------------------------------------------------
    1515////////////////////////////////////////////////////////////////
    16 /////////////////////////////////////////////////////////////////////////
    17 // Classe de lignes/colonnes de matrices
    18 enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
     16
     17//! Class of line, column or diagonal of a TMatrix
     18/*!
     19  A TMatrixRC represents a line, a column or the diagonal of a TMatrix
     20 */
    1921template <class T>
    2022class TMatrixRC {
    2123public:
     24  //! Define type of TMatrixRC
     25  enum TRCKind {
     26    TmatrixRow=0,  //!< TMatrixRC ligne
     27    TmatrixCol=1,  //!< TMatrixRC column
     28    TmatrixDiag=2  //!< TMatrixRC diagonal
     29  };
    2230  TMatrixRC();
    2331  TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
    … …  
    3947  static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
    4048
     49  //! Return the kind of TMatrix (line,column,diagonal)
    4150  TRCKind Kind() const { return kind; }
    4251  uint_4 NElts() const;
    … …  
    6776  static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
    6877
     78  //! Define Absolute value for uint_1
    6979  inline static double Abs_Value(uint_1 v) {return (double) v;}
     80  //! Define Absolute value for uint_2
    7081  inline static double Abs_Value(uint_2 v) {return (double) v;}
     82  //! Define Absolute value for int_2
    7183  inline static double Abs_Value(int_2 v)  {return (v>0)? (double) v: (double) -v;}
     84  //! Define Absolute value for int_4
    7285  inline static double Abs_Value(int_4 v)  {return (v>0)? (double) v: (double) -v;}
     86  //! Define Absolute value for int_8
    7387  inline static double Abs_Value(int_8 v)  {return (v>0)? (double) v: (double) -v;}
     88  //! Define Absolute value for uint_4
    7489  inline static double Abs_Value(uint_4 v) {return (double) v;}
     90  //! Define Absolute value for uint_8
    7591  inline static double Abs_Value(uint_8 v) {return (double) v;}
     92  //! Define Absolute value for r_4
    7693  inline static double Abs_Value(r_4 v)    {return (double) fabsf(v);}
     94  //! Define Absolute value for r_8
    7795  inline static double Abs_Value(r_8 v)    {return fabs(v);}
    78   inline static double Abs_Value(complex<float> v)
     96  //! Define Absolute value for complex r_4
     97  inline static double Abs_Value(complex<r_4> v)
    7998                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
    80   inline static double Abs_Value(complex<double> v)
     99  //! Define Absolute value for complex r_8
     100  inline static double Abs_Value(complex<r_8> v)
    81101                {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
    82102
    83103protected:
    84   TMatrix<T>* matrix;
    85   T*          data;
    86   int_4       index;
    87   uint_4      step;
    88   TRCKind     kind;
     104  TMatrix<T>* matrix;  //!< pointer to the TMatrix
     105  T*          data;    //!< pointer to the beginnig of interesting datas
     106  int_4       index;   //!< index of the line/column
     107  uint_4      step;    //!< step of the line/column
     108  TRCKind     kind;    //!< type: line, column or diagonal
    89109};
    90110
    91111
    92 
     112//! Multiply two TMatrixRC
    93113template <class T>
    94114inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
    … …  
    103123  }
    104124
    105 
     125//! Get the step in datas for a TMatrix for type rckind (line/col/diag)
    106126template <class T>
    107127inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
    … …  
    111131    return 0; }
    112132
     133/*! Get the origin of datas for a TMatrix for type rckind and
     134  number index (line/col/diag). ex: origine for line "index". */
    113135template <class T>
    114136inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
    … …  
    118140    return NULL; }
    119141
     142//! return number of elements for a TMatrixRC
    120143template <class T> inline uint_4 TMatrixRC<T>::NElts() const
    121144  { if (!matrix) return 0;
    … …  
    125148    return 0; }
    126149
     150//! access of element \b i
    127151template <class T>
    128152inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
     153//! access of element \b i
    129154template <class T>
    130155inline T  TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
    131156
    132157////////////////////////////////////////////////////////////////
    133 // Typedef pour simplifier et compatibilite Peida
     158//! Typedef to simplifier TMatrixRC<r_8>
    134159typedef TMatrixRC<r_8> MatrixRC;
    135160
    136161
     162//! Default constructor
    137163template <class T> TMatrixRC<T>::TMatrixRC()
    138164: matrix(NULL), data(NULL), index(0), step(0)
    139165{}
    140166
     167//! Constructor
     168/*!
     169  \param m : matrix
     170  \param rckind : select line, column or diagonal
     171  \param ind : number of the line or column
     172*/
    141173template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,uint_4 ind)
    142174: matrix(&m), data(Org(m,rckind,ind)),
    … …  
    150182// Acces aux rangees et colonnes de matrices
    151183
     184//! Return TMatrixRC for line \b r of matrix \b m
    152185template <class T>
    153186TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, uint_4 r)
    … …  
    157190}
    158191
     192//! Return TMatrixRC for column \b r of matrix \b m
    159193template <class T>
    160194TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, uint_4 c)
    … …  
    164198}
    165199
     200//! Return TMatrixRC for diagonal of matrix \b m
    166201template <class T>
    167202TMatrixRC<T> TMatrixRC<T>::Diag(TMatrix<T> & m)
    … …  
    197232// }
    198233
    199  
     234//! Set column \b c for this TMatrixRC
    200235template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
    201236{
    … …  
    209244}
    210245
     246//! Set line \b r for this TMatrixRC
    211247template <class T> int_4 TMatrixRC<T>::SetRow(int_4 r)
    212248{
    … …  
    220256}
    221257
     258//! Set line diaginal for this TMatrixRC
    222259template <class T> int_4 TMatrixRC<T>::SetDiag()
    223260{
    … …  
    232269}
    233270
    234 
     271//! Operator =
    235272template <class T> TMatrixRC<T>& TMatrixRC<T>::operator = (const TMatrixRC<T>& rc)
    236273{
    … …  
    271308// }
    272309
    273 
     310//! Operator to multiply by constant \b x
    274311template <class T> TMatrixRC<T>& TMatrixRC<T>::operator *= (T x)
    275312{
    … …  
    278315}
    279316
     317//! Operator to divide by constant \b x
    280318template <class T> TMatrixRC<T>& TMatrixRC<T>::operator /= (T x)
    281319{
    … …  
    298336// }
    299337
     338//! Linear combination
     339/*!
     340  Do : \f$ MRC(i) = MRC(i)*a + rc(i)*b \f$
     341  \return *this
     342 */
    300343template <class T>
    301344TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4 first)
    … …  
    309352}
    310353
     354//! Linear combination
     355/*!
     356  Do : \f$ MRC(i) = MRC(i) + rc(i)*b \f$
     357 */
    311358template <class T>
    312359TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4 first)
    … …  
    320367}
    321368
     369//! Find maximum absolute value in TMatrixRC, search begin at \b first
    322370template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4 first)
    323371{
    … …  
    333381}
    334382
     383//! Print on stream \b os
    335384template <class T>
    336385void TMatrixRC<T>::Print(ostream & os) const
    … …  
    346395}
    347396
     397//! Swap two TMatrixRC of the same kind
    348398template <class T>
    349399void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
    … …  
    359409
    360410
    361 
    362 
    363411////////////////////////////////////////////////////////////////
     412// -------------------------------------------------------------
     413//   La classe de calcul simple sur les TMatrix
     414// -------------------------------------------------------------
     415////////////////////////////////////////////////////////////////
     416
    364417//**** Pour inversion
    365418#ifndef M_LN2
    … …  
    373426#endif
    374427
     428//! Gaussian pivoting
     429/*!
     430  Diagonalize matrix \b a, doing the same opreations on matrix \b b
     431  \return determinat of \b a
     432 */ 
    375433template <class T>
    376434T SimpleMatrixOperation<T>::GausPiv(TMatrix<T>& a, TMatrix<T>& b)
    … …  
    394452double nrm = sqrt(vmin*vmax);
    395453if(nrm > 1.e5 || nrm < 1.e-5) {
    396   a /= nrm;
    397   b /= nrm;
     454  a /= (T) nrm;
     455  b /= (T) nrm;
    398456  //cout << "normalisation matrice " << nrm << endl;
    399457} else nrm=1;
    400458
    401 double det = 1.0;
     459T det = 1;
    402460if(nrm != 1) {
    403461  double ld = a.NRows() * log(nrm);
    … …  
    405463   // cerr << "TMatrix warning, overflow for det" << endl;
    406464  } else {
    407     det = exp(ld);
    408   }
    409 }
    410 
    411 TMatrixRC<T> pivRowa(a,TmatrixRow);
    412 TMatrixRC<T> pivRowb(b,TmatrixRow);
     465    det = (T) exp(ld);
     466  }
     467}
     468
     469TMatrixRC<T> pivRowa(a,TMatrixRC<T>::TmatrixRow);
     470TMatrixRC<T> pivRowb(b,TMatrixRC<T>::TmatrixRow);
    413471
    414472for(uint_4 k=0; k<n-1; k++) {
    … …  
    422480    TMatrixRC<T>::Swap(bIPiv,bK);
    423481  }
    424   double pivot = a(k,k);
    425   if (fabs(pivot) < 1.e-50) return 0.0;
     482  T pivot = a(k,k);
     483  if( TMatrixRC<T>::Abs_Value(pivot) < 1.e-50 ) return (T) 0;
    426484  //det *= pivot;
    427485  pivRowa.SetRow(k); // to avoid constructors
    428486  pivRowb.SetRow(k);
    429487  for (uint_4 i=k+1; i<n; i++) {
    430     double r = -a(i,k)/pivot;
     488    T r = -a(i,k)/pivot;
    431489    TMatrixRC<T>::Row(a, i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa
    432490    TMatrixRC<T>::Row(b, i).LinComb(r, pivRowb);
    … …  
    437495// on remonte
    438496for(uint_4 kk=n-1; kk>0; kk--) {
    439   double pivot = a(kk,kk);
    440   if (fabs(pivot) <= 1.e-50) return 0.0;
     497  T pivot = a(kk,kk);
     498  if( TMatrixRC<T>::Abs_Value(pivot) <= 1.e-50 ) return (T) 0;
    441499  pivRowa.SetRow(kk); // to avoid constructors
    442500  pivRowb.SetRow(kk);
    443501  for(uint_4 jj=0; jj<kk; jj++) {
    444     double r = -a(jj,kk)/pivot;
     502    T r = -a(jj,kk)/pivot;
    445503    TMatrixRC<T>::Row(a, jj).LinComb(r, pivRowa);
    446504    TMatrixRC<T>::Row(b, jj).LinComb(r, pivRowb);
    … …  
    449507
    450508for(uint_4 l=0; l<n; l++) {
    451   if (fabs((double)a(l,l)) <= 1.e-50) return 0.0;
     509  if( TMatrixRC<T>::Abs_Value(a(l,l)) <= 1.e-50 ) return (T) 0;
    452510  TMatrixRC<T>::Row(b, l) /= a(l,l);
    453511}
    … …  
    456514}
    457515
     516//! Return the inverse matrix of \b A
    458517template <class T>
    459518TMatrix<T> SimpleMatrixOperation<T>::Inverse(TMatrix<T> const & A)
    460 // Inversion
    461519{
    462520TMatrix<T> a(A);
    463521TMatrix<T> b(a.NCols(),a.NRows());  b = IdentityMatrix(1.);
    464 if(fabs((double)GausPiv(a,b)) < 1.e-50)
    465   throw(MathExc("TMatrix Inverse() Singular OMatrix"));
     522if( TMatrixRC<T>::Abs_Value(GausPiv(a,b)) < 1.e-50)
     523  throw(MathExc("TMatrix Inverse() Singular Matrix"));
    466524return b;
    467525}
    468526
     527
     528////////////////////////////////////////////////////////////////
     529// -------------------------------------------------------------
     530//   La classe fit lineaire
     531// -------------------------------------------------------------
     532////////////////////////////////////////////////////////////////
    469533
    470534LinFitter::LinFitter()
    … …  
    626690#pragma define_template TMatrixRC<r_4>
    627691#pragma define_template TMatrixRC<r_8>
     692#pragma define_template TMatrixRC< complex<r_4> >
     693#pragma define_template TMatrixRC< complex<r_8> >
    628694#pragma define_template SimpleMatrixOperation<int_4>
    629695#pragma define_template SimpleMatrixOperation<r_4>
    630696#pragma define_template SimpleMatrixOperation<r_8>
     697#pragma define_template SimpleMatrixOperation< complex<r_4> >
     698#pragma define_template SimpleMatrixOperation< complex<r_8> >
    631699#endif
    632700
    … …  
    636704template class TMatrixRC<r_4>;
    637705template class TMatrixRC<r_8>;
     706template class TMatrixRC< complex<r_4> >;
     707template class TMatrixRC< complex<r_8> >;
    638708template class SimpleMatrixOperation<int_4>;
    639709template class SimpleMatrixOperation<r_4>;
    640710template class SimpleMatrixOperation<r_8>;
     711template class SimpleMatrixOperation< complex<r_4> >;
     712template class SimpleMatrixOperation< complex<r_8> >;
    641713#endif
Note: See TracChangeset for help on using the changeset viewer.