[772] | 1 | // This may look like C code, but it is really -*- C++ -*-
|
---|
| 2 | #ifndef SOpeMatrix_SEEN
|
---|
| 3 | #define SOpeMatrix_SEEN
|
---|
| 4 |
|
---|
| 5 | #include "machdefs.h"
|
---|
| 6 | #include "tmatrix.h"
|
---|
| 7 | #include "tvector.h"
|
---|
| 8 |
|
---|
| 9 | namespace SOPHYA {
|
---|
| 10 |
|
---|
| 11 |
|
---|
| 12 |
|
---|
| 13 | ////////////////////////////////////////////////////////////////
|
---|
| 14 | template <class T>
|
---|
| 15 | class SimpleMatrixOperation {
|
---|
| 16 | public:
|
---|
| 17 | // Pivot de Gauss : diagonalise la matrice A, en effectuant les memes
|
---|
| 18 | // operations sur la matrice B
|
---|
| 19 | static TMatrix<T> Inverse(TMatrix<T> const & A);
|
---|
| 20 | static T GausPiv(TMatrix<T>& A, TMatrix<T>& B);
|
---|
| 21 | };
|
---|
| 22 |
|
---|
| 23 | // Resolution du systeme A*C = B
|
---|
| 24 | inline r_8 LinSolveInPlace(TMatrix<r_8>& a, TVector<r_8>& b)
|
---|
| 25 | {
|
---|
| 26 | if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
|
---|
| 27 | throw(SzMismatchError("LinSolveInPlace(TMatrix<r_8>,TVector<r_8>) size mismatch"));
|
---|
| 28 | return SimpleMatrixOperation<r_8>::GausPiv(a,b);
|
---|
| 29 | }
|
---|
| 30 |
|
---|
| 31 | // Resolution du systeme A*C = B, avec C retourne dans B
|
---|
| 32 | inline r_8 LinSolve(const TMatrix<r_8>& a, const TVector<r_8>& b, TVector<r_8>& c)
|
---|
| 33 | {
|
---|
| 34 | if(a.NCols() != b.NRows() || a.NCols() != a.NRows())
|
---|
| 35 | throw(SzMismatchError("LinSolve(TMatrix<r_8>,TVector<r_8>) size mismatch"));
|
---|
| 36 | c = b;
|
---|
| 37 | TMatrix<r_8> a1(a);
|
---|
| 38 | return SimpleMatrixOperation<r_8>::GausPiv(a1,c);
|
---|
| 39 | }
|
---|
| 40 |
|
---|
[804] | 41 | //--------------------------------------
|
---|
| 42 | // Linear fitting
|
---|
| 43 | //--------------------------------------
|
---|
[772] | 44 |
|
---|
[804] | 45 | class LinFitter {
|
---|
| 46 | public :
|
---|
| 47 | LinFitter();
|
---|
| 48 | virtual ~LinFitter();
|
---|
| 49 |
|
---|
| 50 | double LinFit(const Vector& x, const Vector& y, int nf,
|
---|
| 51 | double (*f)(int, double), Vector& c);
|
---|
| 52 | // fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1;
|
---|
| 53 |
|
---|
| 54 | double LinFit(const Matrix& fx, const Vector& y, Vector& c);
|
---|
| 55 | // fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
|
---|
| 56 | // la matrice fx contient les valeurs des f:
|
---|
| 57 | // fx(i,j) = f(i, x(j)).
|
---|
| 58 |
|
---|
| 59 | double LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
|
---|
| 60 | double (*f)(int, double), Vector& c, Vector& errC);
|
---|
| 61 | // fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
|
---|
| 62 | // errY2 contient les carres des erreurs sur les Y.
|
---|
| 63 | // au retour, errC contient les erreurs sur les coefs.
|
---|
| 64 |
|
---|
| 65 | double LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
|
---|
| 66 | Vector& c, Vector& errC);
|
---|
| 67 | // fit lineaire des y en tant que somme de c(i)f(i,x), i=0..nf-1,
|
---|
| 68 | // la matrice fx contient les valeurs des f:
|
---|
| 69 | // fx(i,j) = f(i, x(j)).
|
---|
| 70 | // errY2 contient les carres des erreurs sur les Y.
|
---|
| 71 | // au retour, errC contient les erreurs sur les coefs.
|
---|
| 72 | };
|
---|
| 73 |
|
---|
| 74 |
|
---|
[772] | 75 | } // Fin du namespace
|
---|
| 76 |
|
---|
| 77 | #endif
|
---|