source: Sophya/trunk/SophyaLib/TArray/sopemtx.cc@ 882

Last change on this file since 882 was 813, checked in by ansari, 25 years ago

Correction bug/amelioarions TArray,TMatrix,TVector - Reza 5/4/2000
File size: 17.6 KB
Line 
1// C.Magneville, R.Ansari 03/2000
2
3#include "machdefs.h"
4#include <stdio.h>
5#include <iostream.h>
6#include <complex>
7#include <math.h>
8#include "sopemtx.h"
9#include "smathconst.h"
10
11////////////////////////////////////////////////////////////////
12// -------------------------------------------------------------
13// La classe de gestion des lignes et colonnes d'une matrice
14// -------------------------------------------------------------
15////////////////////////////////////////////////////////////////
16/////////////////////////////////////////////////////////////////////////
17// Classe de lignes/colonnes de matrices
18enum TRCKind {TmatrixRow=0, TmatrixCol=1, TmatrixDiag=2};
19template <class T>
20class TMatrixRC {
21public:
22 TMatrixRC();
23 TMatrixRC(TMatrix<T>& m, TRCKind kind, uint_4 index=0);
24 virtual ~TMatrixRC() {}
25
26 // Acces aux rangees et colonnes de matrices
27 static TMatrixRC<T> Row(TMatrix<T> & m, uint_4 r);
28 static TMatrixRC<T> Col(TMatrix<T> & m, uint_4 c);
29 static TMatrixRC<T> Diag(TMatrix<T> & m);
30
31 // ---- A virer $CHECK$ Reza 03/2000
32 // int_4 Next();
33 // int_4 Prev();
34 int_4 SetCol(int_4 c);
35 int_4 SetRow(int_4 r);
36 int_4 SetDiag();
37
38 static uint_4 Step(const TMatrix<T>& m, TRCKind rckind);
39 static T* Org(const TMatrix<T>&, TRCKind rckind, uint_4 ind=0);
40
41 TRCKind Kind() const { return kind; }
42 uint_4 NElts() const;
43 T& operator()(uint_4 i);
44 T operator()(uint_4 i) const;
45
46
47 TMatrixRC<T>& operator = (const TMatrixRC<T>& rc);
48
49// ---- A virer $CHECK$ Reza 03/2000
50// TVector<T> GetVect() const;
51// TMatrixRC<T>& operator += (const TMatrixRC<T>& rc);
52// TMatrixRC<T>& operator -= (const TMatrixRC<T>& rc);
53
54 TMatrixRC<T>& operator *= (T x);
55 TMatrixRC<T>& operator /= (T x);
56
57// TMatrixRC<T>& operator -= (T x);
58// TMatrixRC<T>& operator += (T x);
59
60
61 TMatrixRC<T>& LinComb(T a, T b, const TMatrixRC& rc, uint_4 first=0);
62 TMatrixRC<T>& LinComb(T b, const TMatrixRC<T>& rc, uint_4 first=0);
63
64 uint_4 IMaxAbs(uint_4 first=0);
65 void Print(ostream & os) const ;
66
67 static void Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2);
68
69 inline static double Abs_Value(uint_1 v) {return (double) v;}
70 inline static double Abs_Value(uint_2 v) {return (double) v;}
71 inline static double Abs_Value(int_2 v) {return (v>0)? (double) v: (double) -v;}
72 inline static double Abs_Value(int_4 v) {return (v>0)? (double) v: (double) -v;}
73 inline static double Abs_Value(int_8 v) {return (v>0)? (double) v: (double) -v;}
74 inline static double Abs_Value(uint_4 v) {return (double) v;}
75 inline static double Abs_Value(uint_8 v) {return (double) v;}
76 inline static double Abs_Value(r_4 v) {return (double) fabsf(v);}
77 inline static double Abs_Value(r_8 v) {return fabs(v);}
78 inline static double Abs_Value(complex<float> v)
79 {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
80 inline static double Abs_Value(complex<double> v)
81 {return sqrt(v.real()*v.real()+v.imag()*v.imag());}
82
83protected:
84 TMatrix<T>* matrix;
85 T* data;
86 int_4 index;
87 uint_4 step;
88 TRCKind kind;
89};
90
91
92
93template <class T>
94inline T operator * (const TMatrixRC<T>& a, const TMatrixRC<T>& b)
95 {
96 if ( a.NElts() != b.NElts() )
97 throw(SzMismatchError("TMatrixRC::operator * size mismatch\n"));
98 if ( a.Kind() != b.Kind() )
99 throw(SzMismatchError("TMatrixRC::operator * type mismatch\n"));
100 T sum = 0;
101 for(uint_4 i=0; i<a.NElts(); i++) sum += a(i)*b(i);
102 return sum;
103 }
104
105
106template <class T>
107inline uint_4 TMatrixRC<T>::Step(const TMatrix<T>& m, TRCKind rckind)
108 { switch (rckind) { case TmatrixRow : return m.Step(m.ColsKA());
109 case TmatrixCol : return m.Step(m.RowsKA());
110 case TmatrixDiag : return (m.Step(m.RowsKA())+m.Step(m.ColsKA())); }
111 return 0; }
112
113template <class T>
114inline T* TMatrixRC<T>::Org(const TMatrix<T>& m, TRCKind rckind, uint_4 index)
115 { switch (rckind) { case TmatrixRow : return const_cast<T *>(&(m(index,0)));
116 case TmatrixCol : return const_cast<T *>(&(m(0,index)));
117 case TmatrixDiag : return const_cast<T *>(m.Data()); }
118 return NULL; }
119
120template <class T> inline uint_4 TMatrixRC<T>::NElts() const
121 { if (!matrix) return 0;
122 switch (kind) { case TmatrixRow : return matrix->NCols();
123 case TmatrixCol : return matrix->NRows();
124 case TmatrixDiag : return matrix->NCols(); }
125 return 0; }
126
127template <class T>
128inline T& TMatrixRC<T>::operator()(uint_4 i) {return data[i*step];}
129template <class T>
130inline T TMatrixRC<T>::operator()(uint_4 i) const {return data[i*step];}
131
132////////////////////////////////////////////////////////////////
133// Typedef pour simplifier et compatibilite Peida
134typedef TMatrixRC<r_8> MatrixRC;
135
136
137template <class T> TMatrixRC<T>::TMatrixRC()
138: matrix(NULL), data(NULL), index(0), step(0)
139{}
140
141template <class T> TMatrixRC<T>::TMatrixRC(TMatrix<T>& m,TRCKind rckind,uint_4 ind)
142: matrix(&m), data(Org(m,rckind,ind)),
143 index(ind), step(Step(m,rckind)), kind(rckind)
144{
145if (kind == TmatrixDiag && m.NCols() != m.NRows())
146 throw(SzMismatchError("TMatrixRC::TMatrixRC(...,TmatrixDiag,...) size mismatch\n"));
147}
148
149////////////////////////////////////////////////////////////////
150// Acces aux rangees et colonnes de matrices
151
152template <class T>
153TMatrixRC<T> TMatrixRC<T>::Row(TMatrix<T> & m, uint_4 r)
154{
155TMatrixRC<T> rc(m, TmatrixRow, r);
156return rc;
157}
158
159template <class T>
160TMatrixRC<T> TMatrixRC<T>::Col(TMatrix<T> & m, uint_4 c)
161{
162TMatrixRC<T> rc(m, TmatrixCol, c);
163return rc;
164}
165
166template <class T>
167TMatrixRC<T> TMatrixRC<T>::Diag(TMatrix<T> & m)
168{
169TMatrixRC<T> rc(m, TmatrixDiag);
170return rc;
171}
172
173
174// ---- A virer $CHECK$ Reza 03/2000
175// template <class T> int_4 TMatrixRC<T>::Next()
176// {
177// if (!matrix || kind==TmatrixDiag) return -1;
178// index++;
179// if(kind == TmatrixRow) {
180// if(index > (int_4)matrix->NRows()) {index = (int_4)matrix->NRows(); return -1;}
181// data += matrix->NCols();
182// } else {
183// if (index > (int_4)matrix->NCols()) {index = (int_4)matrix->NCols(); return -1;}
184// data++;
185// }
186// return index;
187// }
188
189// template <class T> int_4 TMatrixRC<T>::Prev()
190// {
191// if (!matrix || kind == TmatrixDiag) return -1;
192// index--;
193// if(index < 0) {index = 0; return -1;}
194// if(kind == TmatrixRow) data -= matrix->NCols();
195// else data--;
196// return index;
197// }
198
199
200template <class T> int_4 TMatrixRC<T>::SetCol(int_4 c)
201{
202if(!matrix) return -1;
203if(c<0 || c>(int_4)matrix->NCols()) return -1;
204kind = TmatrixCol;
205index = c;
206step = Step(*matrix, TmatrixCol);
207data = Org(*matrix, TmatrixCol, c);
208return c;
209}
210
211template <class T> int_4 TMatrixRC<T>::SetRow(int_4 r)
212{
213if(!matrix) return -1;
214if(r<0 && r>(int_4)matrix->NRows()) return -1;
215kind = TmatrixRow;
216index = r;
217step = Step(*matrix, TmatrixRow);
218data = Org(*matrix, TmatrixRow, r);
219return r;
220}
221
222template <class T> int_4 TMatrixRC<T>::SetDiag()
223{
224if (!matrix) return -1;
225if (matrix->NCols() != matrix->NRows())
226 throw(SzMismatchError("TMatrixRC::SetDiag size mismatch\n"));
227kind = TmatrixDiag;
228index = 0;
229step = Step(*matrix, TmatrixDiag);
230data = Org(*matrix, TmatrixDiag);
231return 0;
232}
233
234
235template <class T> TMatrixRC<T>& TMatrixRC<T>::operator = (const TMatrixRC<T>& rc)
236{
237matrix = rc.matrix;
238data = rc.data;
239index = rc.index;
240step = rc.step;
241kind = rc.kind;
242return *this;
243}
244
245// ---- A virer $CHECK$ Reza 03/2000
246// template <class T> TVector<T> TMatrixRC<T>::GetVect() const
247// {
248// TVector<T> v(NElts());
249// for (uint_4 i=0; i<NElts(); i++) v(i) = (*this)(i);
250// return v;
251// }
252
253// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (const TMatrixRC<T>& rc)
254// {
255// if ( NElts() != rc.NElts() )
256// throw(SzMismatchError("TMatrixRC::operator+= size mismatch\n"));
257// if ( kind != rc.kind )
258// throw(SzMismatchError("TMatrixRC::operator+= type mismatch\n"));
259// for (uint_4 i=0; i<NElts(); i++) (*this)(i) += rc(i);
260// return *this;
261// }
262
263// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (const TMatrixRC<T>& rc)
264// {
265// if( NElts() != rc.NElts() )
266// throw(SzMismatchError("TMatrixRC::operator-= size mismatch\n"));
267// if( kind != rc.kind )
268// throw(SzMismatchError("TMatrixRC::operator-= type mismatch\n"));
269// for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= rc(i);
270// return *this;
271// }
272
273
274template <class T> TMatrixRC<T>& TMatrixRC<T>::operator *= (T x)
275{
276for(uint_4 i=0; i<NElts(); i++) (*this)(i) *= x;
277return *this;
278}
279
280template <class T> TMatrixRC<T>& TMatrixRC<T>::operator /= (T x)
281{
282for(uint_4 i=0; i<NElts(); i++) (*this)(i) /= x;
283return *this;
284}
285
286
287// ---- A virer $CHECK$ Reza 03/2000
288// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator -= (T x)
289// {
290// for(uint_4 i=0; i<NElts(); i++) (*this)(i) -= x;
291// return *this;
292// }
293
294// template <class T> TMatrixRC<T>& TMatrixRC<T>::operator += (T x)
295// {
296// for(uint_4 i=0; i<NElts(); i++) (*this)(i) += x;
297// return *this;
298// }
299
300template <class T>
301TMatrixRC<T>& TMatrixRC<T>::LinComb(T a, T b, const TMatrixRC<T>& rc, uint_4 first)
302{
303if ( NElts() != rc.NElts() )
304 throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
305if ( kind != rc.kind )
306 throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
307for(uint_4 i=first; i<NElts(); i++) (*this)(i) = (*this)(i)*a + rc(i)*b;
308return *this;
309}
310
311template <class T>
312TMatrixRC<T>& TMatrixRC<T>::LinComb(T b, const TMatrixRC<T>& rc, uint_4 first)
313{
314if ( NElts() != rc.NElts() )
315 throw(SzMismatchError("TMatrixRC::LinComb size mismatch\n"));
316if ( kind != rc.kind )
317 throw(SzMismatchError("TMatrixRC::LinComb type mismatch\n"));
318for(uint_4 i=first; i<NElts(); i++) (*this)(i) += rc(i)*b;
319return *this;
320}
321
322template <class T> uint_4 TMatrixRC<T>::IMaxAbs(uint_4 first)
323{
324if (first>NElts())
325 throw(SzMismatchError("TMatrixRC::IMaxAbs size mismatch\n"));
326uint_4 imax = first;
327double vmax = Abs_Value((*this)(first));
328for(uint_4 i=first+1; i<NElts(); i++) {
329 double v = Abs_Value((*this)(i));
330 if(v > vmax) {vmax = v; imax = i;}
331}
332return imax;
333}
334
335template <class T>
336void TMatrixRC<T>::Print(ostream & os) const
337{
338 os << " TMatrixRC<T>::Print(ostream & os) " << NElts() << " Kind="
339 << kind << " Index=" << index << " Step= " << step << endl;
340 for(uint_4 i=0; i<NElts(); i++) {
341 os << (*this)(i);
342 if (kind == TmatrixCol) os << endl;
343 else os << ", ";
344 }
345 os << endl;
346}
347
348template <class T>
349void TMatrixRC<T>::Swap(TMatrixRC<T>& rc1, TMatrixRC<T>& rc2)
350{
351if(rc1.NElts() != rc2.NElts())
352 throw(SzMismatchError("TMatrixRC::Swap size mismatch\n"));
353if(rc1.kind != rc2.kind)
354 throw(SzMismatchError("TMatrixRC::Swap type mismatch\n"));
355if(rc1.data == rc2.data) return;
356for(uint_4 i=0; i<rc1.NElts(); i++)
357 {T tmp = rc1(i); rc1(i) = rc2(i); rc2(i) = tmp;}
358}
359
360
361
362
363////////////////////////////////////////////////////////////////
364//**** Pour inversion
365#ifndef M_LN2
366#define M_LN2 0.69314718055994530942
367#endif
368#ifndef LN_MINDOUBLE
369#define LN_MINDOUBLE (M_LN2 * (DMINEXP - 1))
370#endif
371#ifndef LN_MAXDOUBLE
372#define LN_MAXDOUBLE (M_LN2 * DMAXEXP)
373#endif
374
375template <class T>
376T SimpleMatrixOperation<T>::GausPiv(TMatrix<T>& a, TMatrix<T>& b)
377// Pivot de Gauss
378// * Attention: egcs impose que cette fonction soit mise dans le .cc
379// avant ::Inverse() (car Inverse() l'utilise)
380// {TMatrix A(a); TMatrix B(b); return (T) TMatrix::GausPiv(A,B);}
381{
382uint_4 n = a.NRows();
383if(n!=b.NRows())
384 throw(SzMismatchError("TMatrix::GausPiv size mismatch\n"));
385// On fait une normalisation un peu brutale...
386double vmin=MAXDOUBLE;
387double vmax=0;
388for(uint_4 iii=0; iii<a.NRows(); iii++)
389 for(uint_4 jjj=0; jjj<a.NCols(); jjj++) {
390 double v = TMatrixRC<T>::Abs_Value(a(iii,jjj));
391 if(v>vmax) vmax = v;
392 if(v<vmin && v>0) vmin = v;
393}
394double nrm = sqrt(vmin*vmax);
395if(nrm > 1.e5 || nrm < 1.e-5) {
396 a /= nrm;
397 b /= nrm;
398 //cout << "normalisation matrice " << nrm << endl;
399} else nrm=1;
400
401double det = 1.0;
402if(nrm != 1) {
403 double ld = a.NRows() * log(nrm);
404 if (ld <= LN_MINDOUBLE || ld >= LN_MAXDOUBLE) {
405 // cerr << "TMatrix warning, overflow for det" << endl;
406 } else {
407 det = exp(ld);
408 }
409}
410
411TMatrixRC<T> pivRowa(a,TmatrixRow);
412TMatrixRC<T> pivRowb(b,TmatrixRow);
413
414for(uint_4 k=0; k<n-1; k++) {
415 uint_4 iPiv = TMatrixRC<T>::Col(a, k).IMaxAbs(k);
416 if(iPiv != k) {
417 TMatrixRC<T> aIPiv(TMatrixRC<T>::Row(a,iPiv));
418 TMatrixRC<T> aK(TMatrixRC<T>::Row(a, k));
419 TMatrixRC<T>::Swap(aIPiv,aK);
420 TMatrixRC<T> bIPiv(TMatrixRC<T>::Row(b, iPiv));
421 TMatrixRC<T> bK(TMatrixRC<T>::Row(b, k));
422 TMatrixRC<T>::Swap(bIPiv,bK);
423 }
424 double pivot = a(k,k);
425 if (fabs(pivot) < 1.e-50) return 0.0;
426 //det *= pivot;
427 pivRowa.SetRow(k); // to avoid constructors
428 pivRowb.SetRow(k);
429 for (uint_4 i=k+1; i<n; i++) {
430 double r = -a(i,k)/pivot;
431 TMatrixRC<T>::Row(a, i).LinComb(r, pivRowa); // + rapide que -= r * pivRowa
432 TMatrixRC<T>::Row(b, i).LinComb(r, pivRowb);
433 }
434}
435det *= a(n-1, n-1);
436
437// on remonte
438for(uint_4 kk=n-1; kk>0; kk--) {
439 double pivot = a(kk,kk);
440 if (fabs(pivot) <= 1.e-50) return 0.0;
441 pivRowa.SetRow(kk); // to avoid constructors
442 pivRowb.SetRow(kk);
443 for(uint_4 jj=0; jj<kk; jj++) {
444 double r = -a(jj,kk)/pivot;
445 TMatrixRC<T>::Row(a, jj).LinComb(r, pivRowa);
446 TMatrixRC<T>::Row(b, jj).LinComb(r, pivRowb);
447 }
448}
449
450for(uint_4 l=0; l<n; l++) {
451 if (fabs((double)a(l,l)) <= 1.e-50) return 0.0;
452 TMatrixRC<T>::Row(b, l) /= a(l,l);
453}
454
455return det;
456}
457
458template <class T>
459TMatrix<T> SimpleMatrixOperation<T>::Inverse(TMatrix<T> const & A)
460// Inversion
461{
462TMatrix<T> a(A);
463TMatrix<T> b(a.NCols(),a.NRows()); b = IdentityMatrix(1.);
464if(fabs((double)GausPiv(a,b)) < 1.e-50)
465 throw(MathExc("TMatrix Inverse() Singular OMatrix"));
466return b;
467}
468
469
470LinFitter::LinFitter()
471{
472}
473
474LinFitter::~LinFitter()
475{
476}
477
478double LinFitter::LinFit(const Vector& x, const Vector& y, int nf,
479 double (*f)(int, double), Vector& c)
480{
481 int n = x.NElts();
482 if (n != y.NElts()) THROW(sizeMismatchErr);
483
484 Matrix fx(nf, n);
485 for (int i=0; i<nf; i++)
486 for (int j=0; j<n; j++)
487 fx(i,j) = f(i,x(j));
488
489 return LinFit(fx,y,c);
490}
491
492
493
494double LinFitter::LinFit(const Matrix& fx, const Vector& y, Vector& c)
495{
496 int n = y.NElts();
497 if ( n != fx.NCol()) THROW(sizeMismatchErr);
498
499 int nf = fx.NRows();
500
501 Matrix a(nf,nf);
502
503 for (int j=0; j<nf; j++)
504 for (int k=j; k<nf; k++)
505 a(j,k) = a(k,j) = TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), j)
506
507 * TMatrixRC<r_8>::Row(const_cast<Matrix &>(fx), k); /* $CHECK$ Reza 10/3/2000 */
508
509 c = fx * y;
510
511 if (SimpleMatrixOperation<r_8>::GausPiv(a,c) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
512
513 double xi2 = 0;
514
515 for (int k=0; k<n; k++) {
516 double x = 0;
517 for (int i=0; i<nf; i++)
518 x += c(i) * fx(i,k);
519 x -= y(k);
520 xi2 += x*x;
521 }
522 return xi2;
523}
524
525
526
527double LinFitter::LinFit(const Vector& x, const Vector& y, const Vector& errY2, int nf,
528 double (*f)(int, double), Vector& c, Vector& errC)
529{
530 int n = x.NElts();
531 if (n != y.NElts()) THROW(sizeMismatchErr);
532
533 Matrix fx(nf, n);
534 for (int i=0; i<nf; i++)
535 for (int j=0; j<n; j++)
536 fx(i,j) = f(i,x(j));
537
538 return LinFit(fx,y,errY2,c,errC);
539}
540
541
542double LinFitter::LinFit(const Matrix& fx, const Vector& y, const Vector& errY2,
543 Vector& c, Vector& errC)
544{
545 int n = y.NElts();
546 if ( n != errY2.NElts()) THROW(sizeMismatchErr);
547 if ( n != fx.NCol()) THROW(sizeMismatchErr);
548
549 int nf = fx.NRows();
550
551 Matrix a(nf,nf);
552
553 c.Realloc(nf);
554 errC.Realloc(nf);
555
556 for (int j=0; j<nf; j++)
557 for (int k=j; k<nf; k++) {
558 double x=0;
559 for (int l=0; l<n; l++)
560 x += fx(j,l) * fx(k,l) / errY2(l); // Matrice a inverser
561 a(j,k) = a(k,j) = x;
562 }
563
564 Matrix d(nf,nf+1);
565 for (int k=0; k<nf; k++) {
566 double x=0;
567 for (int l=0; l<n; l++)
568 x += fx(k,l) * y(l) / errY2(l); // Second membre 1ere colonne
569 d(k,0) = x;
570 for (int m=1; m<=nf; m++)
571 d(k,m) = (k==m) ? 1 : 0; // Reste second membre = Id.
572 }
573
574 if (SimpleMatrixOperation<r_8>::GausPiv(a,d) == 0) THROW(singMatxErr); /* $CHECK$ Reza 10/3/2000 */
575
576 for (int l=0; l<nf; l++) {
577 c(l) = d(l,0); // Parametres = 1ere colonne
578 errC(l) = d(l,l+1); // Erreurs = diag inverse.
579 }
580
581 double xi2 = 0;
582
583 for (int jj=0; jj<n; jj++) {
584 double x = 0;
585 for (int ii=0; ii<nf; ii++)
586 x += c(ii) * fx(ii,jj);
587 x -= y(jj);
588 xi2 += x*x/errY2(jj);
589 }
590 return xi2;
591}
592
593
594void _ZZ_TestTMatrixRC_YY_(TMatrix<r_8> & m)
595{
596 cout << " ::::: _ZZ_TestTMatrixRC_YY_ :::: M= \n" << m << endl;
597 TMatrixRC<r_8> l0 = TMatrixRC<r_8>::Row(m,0);
598 cout << "TMatrixRC<r_8>::Row(m,0) = \n" ;
599 l0.Print(cout);
600 TMatrixRC<r_8> l1 = TMatrixRC<r_8>::Row(m,1);
601 cout << "TMatrixRC<r_8>::Row(m,1) = \n" ;
602 l1.Print(cout);
603 l0.SetRow(2);
604 cout << "TMatrixRC<r_8>::l0.SetRow(2 = \n" ;
605 l0.Print(cout);
606
607 TMatrixRC<r_8> c0 = TMatrixRC<r_8>::Col(m,0);
608 cout << "TMatrixRC<r_8>::Col(m,0) = \n" ;
609 c0.Print(cout);
610 TMatrixRC<r_8> c1 = TMatrixRC<r_8>::Col(m,1);
611 cout << "TMatrixRC<r_8>::Col(m,1) = \n" ;
612 c1.Print(cout);
613 c0.SetCol(2);
614 cout << "TMatrixRC<r_8>::c0.SetCol(2) = \n" ;
615 c0.Print(cout);
616 TMatrixRC<r_8> c00 = TMatrixRC<r_8>::Col(m,0);
617 TMatrixRC<r_8>::Swap(c0, c00);
618 cout << " ::::: M Apres Swap = \n" << m << endl;
619
620}
621
622///////////////////////////////////////////////////////////////
623#ifdef __CXX_PRAGMA_TEMPLATES__
624// Instances gestion lignes/colonnes
625#pragma define_template TMatrixRC<int_4>
626#pragma define_template TMatrixRC<r_4>
627#pragma define_template TMatrixRC<r_8>
628#pragma define_template SimpleMatrixOperation<int_4>
629#pragma define_template SimpleMatrixOperation<r_4>
630#pragma define_template SimpleMatrixOperation<r_8>
631#endif
632
633#if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES)
634// Instances gestion lignes/colonnes
635template class TMatrixRC<int_4>;
636template class TMatrixRC<r_4>;
637template class TMatrixRC<r_8>;
638template class SimpleMatrixOperation<int_4>;
639template class SimpleMatrixOperation<r_4>;
640template class SimpleMatrixOperation<r_8>;
641#endif
Note: See TracBrowser for help on using the repository browser.