Changeset 299 in Sophya for trunk/SophyaLib

Timestamp:

May 17, 1999, 7:16:44 PM (26 years ago)

Author:

ansari

Message:

TMatrix et TVector suite cmv 17/5/99

Location:

trunk/SophyaLib/NTools

Files:

• tmatrix.cc (modified) (7 diffs)
• tmatrix.h (modified) (2 diffs)
• tvector.cc (added)
• tvector.h (added)

trunk/SophyaLib/NTools/tmatrix.cc

@@ -1 @@
-// $Id: tmatrix.cc,v 1.4 1999-05-05 05:51:54 ansari Exp $
+// $Id: tmatrix.cc,v 1.5 1999-05-17 17:16:43 ansari Exp $
 //                         C.Magneville          04/99
 #include "machdefs.h"
@@ -9 +9 @@
 #include "tmatrix.h"
 #include "objfio.h"
+#include "generalfit.h"
 
 using namespace PlanckDPC;
@@ -130 +131 @@
 {
 if(mNr!=b.mNr || mNc!=b.mNc)
-    throw(SzMismatchError("NDataBlock operator C=A+B size mismatch\n"));
+    throw(SzMismatchError("TMatrix operator C=A+B size mismatch\n"));
 TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
 result.mNDBlock = mNDBlock+b.mNDBlock;
@@ -139 +140 @@
 {
 if(mNr!=b.mNr || mNc!=b.mNc)
-    throw(SzMismatchError("NDataBlock operator C=A-B size mismatch\n"));
+    throw(SzMismatchError("TMatrix operator C=A-B size mismatch\n"));
 TMatrix<T> result; result.SetTemp(true); result.mNr=mNr; result.mNc=mNc;
 result.mNDBlock = mNDBlock-b.mNDBlock;
@@ -149 +150 @@
 {
 if(mNr==0 || mNc==0 || b.mNc==0 || mNc!=b.mNr)
-  throw(SzMismatchError("NDataBlock operator C=A*B size mismatch\n"));
+  throw(SzMismatchError("TMatrix operator C=A*B size mismatch\n"));
 TMatrix<T> r; r.SetTemp(true); r.ReSize(mNr,b.mNc);
 T *ai,*aik,*bj,*bkj,*ri,*rij;
@@ -160 +161 @@
 }
 return r;
+}
+
+#include "matrix.h"
+////////////////////////////////////////////////////////////////
+//**** Pour inversion
+TMatrix<r_8> TMatrix<r_8>::Inverse() const
+// Inversion
+{
+TMatrix<r_8> b(mNc,mNr);
+TMatrix<r_8> a(*this);
+b = 1.;
+if (fabs(TMatrix<r_8>::GausPiv(a,b)) < 1.e-50)
+  throw(MathExc("TMatrix Inverse() Singular Matrix"));
+return b;
+}
+
+double TMatrix<r_8>::GausPiv(TMatrix<r_8>& a, TMatrix<r_8>& b)
+// Pivot de Gauss
+{
+Matrix A(a);
+Matrix B(b);
+return Matrix::GausPiv(A,B);
+}
+
+//////////////////////////////////////////////////////////
+//**** Residus des fits
+TMatrix<r_8> TMatrix<r_8>::FitResidus(GeneralFit& gfit)
+//      Retourne une classe contenant les residus du fit ``gfit''.
+//      On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
+{
+if(NCols()<=0||NRows()<=0)
+  throw(SzMismatchError("TMatrix::FitResidus size mismatch\n"));
+GeneralFunction* f = gfit.GetFunction();
+if(f==NULL)
+  throw(NullPtrError("TMatrix::FitResidus GeneraFit==NULL\n"));
+int npar =  gfit.GetNPar();
+if(npar==0)
+  throw(SzMismatchError("TMatrix::FitResidus GeneraFit 0 parametre\n"));
+double* par = new double[npar];
+{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
+TMatrix<r_8> m(*this);
+for(int i=0;i<NRows();i++) for(int j=0;j<NCols();j++) {
+  double x[2] = {(double)j,(double)i}; // attention x=col y=row!!
+  m(i,j) -= f->Value(x,par);
+}
+delete [] par;
+return m;
+}
+
+TMatrix<r_8> TMatrix<r_8>::FitFunction(GeneralFit& gfit)
+//      Retourne une classe contenant la fonction du fit ``gfit''.
+//      On suppose que x=j (colonnes) et y=i (lignes) pour m(i,j).
+{
+if(NCols()<=0||NRows()<=0)
+  throw(SzMismatchError("TMatrix::FitFunction size mismatch\n"));
+GeneralFunction* f = gfit.GetFunction();
+if(f==NULL)
+  throw(NullPtrError("TMatrix::FitFunction GeneraFit==NULL\n"));
+int npar =  gfit.GetNPar();
+if(npar==0)
+  throw(SzMismatchError("TMatrix::FitFunction GeneraFit 0 parametre\n"));
+double* par = new double[npar];
+{for(int i=0;i<npar;i++) par[i] = gfit.GetParm(i);}
+TMatrix<r_8> m(*this);
+for(int i=0;i<NRows();i++) for(int j=0;j<NCols();j++) {
+  double x[2] = {(double)j,(double)i};
+  m(i,j) = f->Value(x,par);
+}
+delete [] par;
+return m;
 }
 
@@ -235 +306 @@
 os.Put(itab,3);
 // On ecrit le NDataBlock
-//cmv   encore des problemes avec les consteries
 FIO_NDataBlock<T> fio_nd(&dobj->DataBlock());
 fio_nd.Write(os);

trunk/SophyaLib/NTools/tmatrix.h

@@ -10 +10 @@
 #include "anydataobj.h"
 #include "ndatablock.h"
+class GeneralFit;
 
 namespace PlanckDPC {
@@ -93 +94 @@
   TMatrix<T> Mul(const TMatrix<T>& b) const;
 
+  // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes
+  // operations sur la matrice B
+  TMatrix<T> Inverse() const;
+  static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);
+
+  // Residus et fonction fittees.
+  TMatrix<T> FitResidus(GeneralFit& gfit);
+  TMatrix<T> FitFunction(GeneralFit& gfit);
+
 protected:
   // partage les donnees si "a" temporaire, clone sinon.