ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/solveur/src/matpleine.cpp
Revision: 1158
Committed: Thu Jun 13 22:18:49 2024 UTC (11 months ago) by francois
File size: 22443 byte(s)
Error occurred while calculating annotation data.
Log Message:
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

# Content
1 //####//------------------------------------------------------------
2 //####//------------------------------------------------------------
3 //####// MAGiC
4 //####// Jean Christophe Cuilliere et Vincent FRANCOIS
5 //####// Departement de Genie Mecanique - UQTR
6 //####//------------------------------------------------------------
7 //####// MAGIC est un projet de recherche de l equipe ERICCA
8 //####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9 //####// http://www.uqtr.ca/ericca
10 //####// http://www.uqtr.ca/
11 //####//------------------------------------------------------------
12 //####//------------------------------------------------------------
13 //####//
14 //####// matpleine.cpp
15 //####//
16 //####//------------------------------------------------------------
17 //####//------------------------------------------------------------
18 //####// COPYRIGHT 2000-2024
19 //####// jeu 13 jun 2024 11:58:57 EDT
20 //####//------------------------------------------------------------
21 //####//------------------------------------------------------------
22 #include "gestionversion.h"
23
24
25 #pragma hdrstop
26 #include <iomanip>
27 #include "matpleine.h"
28 #pragma package(smart_init)
29
30 #include <math>
31
32
33
34
35 #include <fstream>
36 #include <cmath>
37
38
39 Matrice::Matrice()
40 {
41 Data = NULL;
42 NbCols = NbRows = 0;
43 }
44
45
46 Matrice::Matrice( int Rows, int Cols,double InitVal)
47 {
48 NbRows = Rows;
49 NbCols = Cols;
50
51 Data = new double*[NbRows];
52 for(int i=0; i<NbRows; ++i) {
53 Data[i] = new double[NbCols];
54 for(int j=0; j<NbCols; ++j) {
55 Data[i][j] = InitVal;
56 }
57 }
58 }
59
60
61 Matrice::Matrice(double* m_pData, int Rows, int Cols)
62 {
63 NbRows = Rows;
64 NbCols = Cols;
65
66 Data = new double*[NbRows];
67 for(int i=0; i<NbRows; ++i) {
68 Data[i] = new double[NbCols];
69 for(int j=0; j<NbCols; ++j) {
70 Data[i][j] = m_pData[i*NbCols+j];
71 }
72 }
73 }
74
75 Matrice::Matrice(double** m_nData, int Rows, int Cols)
76 {
77 NbRows = Rows;
78 NbCols = Cols;
79
80 Data = new double*[NbRows];
81 for(int i=0; i<NbRows; ++i) {
82 Data[i] = new double[NbCols];
83 for(int j=0; j<NbCols; ++j) {
84 Data[i][j] = m_nData[i][j];
85 }
86 }
87 }
88
89
90 Matrice::Matrice(const Matrice& obj)
91 {
92 NbRows = obj.NbRows;
93 NbCols = obj.NbCols;
94
95 Data = new double*[NbRows];
96 for(int i=0; i<NbRows; ++i) {
97 Data[i] = new double[NbCols];
98 for(int j=0; j<NbCols; ++j) {
99 Data[i][j] = obj.Data[i][j];
100 }
101 }
102 }
103
104
105 Matrice::~Matrice()
106 {
107 if(Data){
108 for(int i=0; i<NbRows; ++i) {
109 delete[] Data[i];
110 Data[i] = NULL;
111 }
112 delete[] Data;
113 }
114 Data = NULL;
115 }
116
117 Matrice& Matrice::operator =(const Matrice& obj)
118 {
119 if(&obj == this) return *this;
120
121 if(NbRows || NbCols) {
122 for(int i=0; i<NbRows; ++i) {
123 delete[] Data[i];
124 }
125 delete[] Data;
126 }
127
128 NbRows = obj.NbRows;
129 NbCols = obj.NbCols;
130
131 Data = new double*[NbRows];
132 for(int i=0; i<NbRows; ++i) {
133 Data[i] = new double[NbCols];
134 for(int j=0; j<NbCols; ++j) {
135 Data[i][j] = obj.Data[i][j];
136 }
137 }
138
139 return *this;
140 }
141
142 Matrice& Matrice::operator +(const Matrice& obj) const
143 {
144 assert((NbRows == obj.NbRows)&& (NbCols == obj.NbCols)) ;
145
146
147 static Matrice temp;
148 temp = Matrice( NbRows, NbCols,0.0);
149
150 for(int i=0; i<NbRows; ++i) {
151 for(int j=0; j<NbCols; ++j) {
152 temp.Data[i][j] = Data[i][j] + obj.Data[i][j];
153 }
154 }
155
156 return temp;
157 }
158
159
160 Matrice& Matrice::operator -(const Matrice& obj) const
161 {
162 assert((NbRows == obj.NbRows)&& (NbCols == obj.NbCols)) ;
163
164 static Matrice temp;
165 temp = Matrice(NbRows, NbCols,0.);
166
167 for(int i=0; i<NbRows; ++i) {
168 for(int j=0; j<NbCols; ++j) {
169 temp.Data[i][j] = Data[i][j] - obj.Data[i][j];
170 }
171 }
172
173 return temp;
174 }
175
176
177 Matrice& Matrice::operator *(const Matrice& obj) const
178 {
179 assert(NbCols == obj.NbRows);
180
181 double sum = double(0);
182 double prod = double(1);
183
184 static Matrice temp;
185 temp = Matrice( NbRows, obj.NbCols,0.0);
186
187 for(int i=0; i<temp.NbRows; ++i) {
188 for(int j=0; j<temp.NbCols; ++j) {
189 sum = double(0);
190 for(int q=0; q<NbCols; ++q) {
191 prod = Data[i][q] * obj.Data[q][j];
192 sum += prod;
193 }
194 temp.Data[i][j] = sum;
195 }
196 }
197
198 return temp;
199 }
200
201 Matrice& Matrice::operator *(const double _d) const
202 {
203 static Matrice temp;
204 temp = Matrice(Data, NbRows,NbCols);
205
206 for(int i=0; i<NbRows; ++i) {
207 for(int j=0; j<NbCols; ++j) {
208 temp.Data[i][j] *= _d;
209 }
210 }
211
212 return temp;
213 }
214
215
216
217
218
219 Matrice& Matrice::operator /(const double _d) const
220 {
221 assert(_d != 0) ;
222
223 static Matrice temp;
224 temp = Matrice(Data, NbCols, NbRows);
225
226 for(int i=0; i<temp.NbRows; ++i) {
227 for(int j=0; j<temp.NbCols; ++j) {
228 temp.Data[i][j] /= _d;
229 }
230 }
231
232 return temp;
233 }
234
235
236
237 Matrice& Matrice::operator +=(const Matrice& obj)
238 {
239 return(*this = *this + obj);
240 }
241
242 Matrice& Matrice::operator -=(const Matrice& obj)
243 {
244 return(*this = *this - obj);
245 }
246
247 Matrice& Matrice::operator *=(const Matrice& obj)
248 {
249 return(*this = *this * obj);
250 }
251
252 Matrice& Matrice::operator *=(const double _d)
253 {
254 return(*this = *this * _d);
255 }
256
257
258
259 Matrice& Matrice::operator /=(const double _d)
260 {
261 return(*this = *this / _d);
262 }
263
264
265
266
267
268 bool Matrice::operator ==(const Matrice& obj) const
269 {
270 if(NbRows != obj.NbRows) return false;
271 if(NbCols != obj.NbCols) return false;
272
273 for(int i=0; i<NbRows; ++i) {
274 for(int j=0; j<NbCols; ++j) {
275 if(Data[i][j] != obj.Data[i][j]) return false;
276 }
277 }
278
279 return true;
280 }
281
282 bool Matrice ::operator !=(const Matrice& obj) const
283 {
284 if(NbRows != obj.NbRows) return true;
285 if(NbCols != obj.NbCols) return true;
286
287 for(int i=0; i<NbRows; ++i) {
288 for(int j=0; j<NbCols; ++j) {
289 if(Data[i][j] != obj.Data[i][j]) return true;
290 }
291 }
292
293 return false;
294 }
295
296 double* Matrice ::operator [](const int _i) const
297 {
298 assert((_i < NbRows)&& (_i >= 0));
299
300 return Data[_i];
301 }
302
303
304 double& Matrice ::operator ()(const int _i, const int _j) const
305 {
306 assert((_i < NbRows) && (_j < NbCols));
307 assert((_i >= 0) && (_j >= 0));
308
309 return Data[_i][_j];
310 }
311
312 bool Matrice::IsIdentity() const
313 {
314 if(NbCols != NbRows) return false;
315
316 for(int i=0; i<NbCols; ++i) {
317 for(int j=0; j<NbRows; ++j) {
318 if(i == j) {
319 if(Data[i][j] != 1) return false;
320 }
321 else if(Data[i][j] != 0) return false;
322 }
323 }
324
325 return true;
326 }
327
328 bool Matrice::IsZeroMatrice() const
329 {
330 for(int i=0; i<NbRows; ++i) {
331 for(int j=0; j<NbCols; ++j) {
332 if(Data[i][j] != 0) return false;
333 }
334 }
335
336 return true;
337 }
338
339
340
341 double Matrice::Get_SumAll() const
342 {
343 double sum = 0;
344
345 for(int i=0; i<NbRows; ++i) {
346 for(int j=0; j<NbCols; ++j) {
347 sum += Data[i][j];
348 }
349 }
350
351 return sum;
352 }
353
354
355
356
357 double Matrice::Get_SumRow(const int Row) const
358 {
359 double sum = double(0.);
360
361 for(int j=0; j<NbCols; ++j) {
362 sum += Data[Row][j];
363 }
364
365 return sum;
366 }
367
368
369 double Matrice::Get_SumColumn(const int Col) const
370 {
371 double sum = double(0);
372
373 for(int i=0; i<NbRows; ++i) {
374 sum += Data[i][Col];
375 }
376
377 return sum;
378 }
379
380
381
382
383
384
385 double Matrice::Get_Max() const
386 {
387 double max = Data[0][0];
388
389 for(int i=0; i<NbRows; ++i) {
390 for(int j=0; j<NbCols; ++j) {
391 if(Data[i][j] > max) max = Data[i][j];
392 }
393 }
394
395 return max;
396 }
397
398 double Matrice::Get_AbMax() const
399 {
400 double max = fabs(Data[0][0]);
401
402 for(int i=0; i<NbRows; ++i) {
403 for(int j=0; j<NbCols; ++j) {
404 if(fabs(Data[i][j]) > max) max = fabs(Data[i][j]);
405 }
406 }
407
408 return max;
409 }
410
411
412 double Matrice ::Get_Min() const
413 {
414 double min = Data[0][0];
415
416 for(int i=0; i<NbRows; ++i) {
417 for(int j=0; j<NbCols; ++j) {
418 if(Data[i][j] < min) min = Data[i][j];
419 }
420 }
421
422 return min;
423 }
424
425 double Matrice::Get_RowMax(const int Row) const
426 {
427 double max = Data[Row][0];
428
429 for(int j=0; j<NbCols; ++j) {
430 if(Data[Row][j] > max) max = Data[Row][j];
431 }
432
433 return max;
434 }
435
436 double Matrice::Get_RowMin(const int Row) const
437 {
438 double min = Data[Row][0];
439
440 for(int j=0; j<NbCols; ++j) {
441 if(Data[Row][j] < min) min = Data[Row][j];
442 }
443
444 return min;
445 }
446
447 double Matrice::Get_ColumnMax(const int Col) const
448 {
449 double max = Data[0][Col];
450
451 for(int i=0; i<NbRows; ++i) {
452 if(Data[i][Col] > max) max = Data[i][Col];
453 }
454
455 return max;
456 }
457
458 double Matrice::Get_ColumnMin(const int Col) const
459 {
460 double min = Data[0][Col];
461
462 for(int i=0; i<NbRows; ++i) {
463 if(Data[i][Col] < min) min = Data[i][Col];
464 }
465
466 return min;
467 }
468
469
470
471 double* Matrice ::Get_DataOneDimen() const
472 {
473 double* newData;
474
475 newData = new double[NbRows * NbCols];
476
477 for(int i=0; i<NbRows; ++i) {
478 for(int j=0; j<NbCols; ++j) {
479 newData[(i*NbRows)+j] = Data[i][j];
480 }
481 }
482
483 return newData;
484 }
485
486
487 double** Matrice ::Get_DataTwoDimen() const
488 {
489 double** newData;
490
491 newData = new double*[NbRows];
492 for(int i=0; i<NbRows; ++i) {
493 newData[i] = new double[NbCols];
494 for(int j=0; j<NbCols; ++j) {
495 newData[i][j] = Data[i][j];
496 }
497 }
498
499 return newData;
500 }
501
502 int Matrice ::Get_NbRows() const
503 {
504 return NbRows;
505 }
506
507 int Matrice::Get_NbCols() const
508 {
509 return NbCols;
510 }
511
512 Matrice& Matrice::Clear()
513 {
514 for(int i=0; i<NbRows; ++i) {
515 for(int j=0; j<NbCols; ++j) {
516 Data[i][j] = 0;
517 }
518 }
519
520 return *this;
521 }
522
523
524 Matrice& Matrice::ClearRow(const int Row)
525 {
526 for(int j=0; j<NbCols; ++j) {
527 Data[Row][j] = 0;
528 }
529
530 return *this;
531 }
532
533
534 Matrice& Matrice::ClearColumn(const int Col)
535 {
536 for(int i=0; i<NbRows; ++i) {
537 Data[i][Col] = 0;
538 }
539
540 return *this;
541 }
542
543 Matrice& Matrice::Fill(const double _d)
544 {
545 for(int i=0; i<NbRows; ++i) {
546 for(int j=0; j<NbCols; ++j) {
547 Data[i][j] = _d;
548 }
549 }
550
551 return *this;
552 }
553
554 Matrice& Matrice::FillRow(const int Row, const double _d)
555 {
556 for(int j=0; j<NbCols; ++j) {
557 Data[Row][j] = _d;
558 }
559
560 return *this;
561 }
562
563 Matrice& Matrice::FillColumn(const int Col, const double _d)
564 {
565 for(int i=0; i<NbRows; ++i) {
566 Data[i][Col] = _d;
567 }
568
569 return *this;
570 }
571
572
573
574
575
576 Matrice& Matrice::Get_SubMatrice(const int RowSpot, const int ColSpot, const int RowLen, const int ColLen) const
577 {
578 static Matrice temp;
579 temp = Matrice( RowLen, ColLen,0.0);
580
581 for(int i=RowSpot, k=0; i<(RowLen+RowSpot); i++, k++) {
582 for(int j=ColSpot, l=0; j<(ColLen+ColSpot); j++, l++) {
583 temp(k,l) = Data[i][j];
584 }
585 }
586
587 return temp;
588 }
589
590 Matrice& Matrice::Set_SubMatrice(Matrice& sub,const int RowSpot, const int ColSpot, const int RowLen, const int ColLen)
591 {
592
593 for(int i=RowSpot, k=0;i<RowSpot+RowLen;i++,++k)
594 {
595 for(int j=ColSpot, l=0;j< ColSpot+ColLen;j++,++l)
596 Data[i][j]=sub(k,l);
597 }
598
599 return *this;
600 }
601
602 Matrice Matrice::Get_SubMatrice(vecteur& lgn,vecteur& col)
603 {
604 int t_tmp=lgn.getDim();
605 Matrice tmp(t_tmp,t_tmp,0.);
606 for(int i=0;i<t_tmp;i++)
607 {
608 int I= lgn[i];
609 for(int j=0;j<t_tmp;j++)
610 {
611 int J= col[j];
612 tmp.Data[i][j] = Data[I][J];
613 }
614 }
615 return tmp;
616 }
617
618 Matrice& Matrice::Get_SwapRows(const int Row1, const int Row2)
619 {
620 double* temp;
621
622 temp = new double[NbCols];
623
624 for(int j=0; j<NbCols; ++j) {
625 temp[j] = Data[Row1][j];
626 Data[Row1][j] = Data[Row2][j];
627 Data[Row2][j] = temp[j];
628 }
629
630 delete[] temp;
631
632 return *this;
633 }
634
635 Matrice& Matrice::Get_SwapCols(const int Col1, const int Col2)
636 {
637 double* temp;
638
639 temp = new double[NbRows];
640
641 for(int i=0; i<NbRows; ++i) {
642 temp[i] = Data[i][Col1];
643 Data[i][Col1] = Data[i][Col2];
644 Data[i][Col2] = temp[i];
645 }
646
647 delete[] temp;
648
649 return *this;
650 }
651
652 Matrice& Matrice::Get_Transpose() const
653 {
654 static Matrice temp;
655 temp = Matrice(NbCols, NbRows,0.0);
656
657 for(int i=0; i<NbRows; ++i) {
658 for(int j=0; j<NbCols; ++j) {
659 temp.Data[j][i] = Data[i][j];
660 }
661 }
662
663 return temp;
664 }
665
666 Matrice& Matrice::Transpose()
667 {
668 *this = this->Get_Transpose();
669
670 return *this;
671 }
672
673
674
675
676
677 Matrice& Matrice::Get_ConMatAsRow(const Matrice& obj)
678 {
679 assert(NbCols == obj.NbCols);// ErrorMsg("mismatched matrices in row concatenation", true);
680
681 Matrice temp((NbRows + obj.NbRows), NbCols,0.0);
682 int i;
683 for(int i=0; i<NbRows; ++i) {
684 for(int j=0; j<NbCols; ++j) {
685 temp.Data[i][j] = Data[i][j];
686 }
687 }
688 for(int k=0; i<temp.NbRows; ++i, ++k) {
689 for(int j=0; j<NbCols; ++j) {
690 temp.Data[i][j] = obj.Data[k][j];
691 }
692 }
693
694 *this = temp;
695
696 return *this;
697 }
698
699 Matrice& Matrice::Get_ConMatAsCol(const Matrice& obj)
700 {
701 assert(NbRows == obj.NbRows);// ErrorMsg("mismatched matrices in column concatenation", true);
702
703 Matrice temp(NbRows, (NbCols + obj.NbCols),0.0);
704 int i,j;
705 for(i=0; i<NbRows; ++i) {
706 for( j=0; j<NbCols; ++j) {
707 temp.Data[i][j] = Data[i][j];
708 }
709 for(int l=0; l<obj.NbCols; ++l, ++j) {
710 temp.Data[i][j] = obj.Data[i][l];
711 }
712 }
713
714 *this = temp;
715
716 return *this;
717 }
718
719 Matrice& Matrice::GetCMAR(const Matrice& obj) const
720 {
721 static Matrice temp;
722
723 temp = *this;
724 temp.Get_ConMatAsRow(obj);
725
726 return temp;
727 }
728
729 Matrice& Matrice::GetCMAC(const Matrice& obj) const
730 {
731 static Matrice temp;
732
733 temp = *this;
734 temp.Get_ConMatAsCol(obj);
735
736 return temp;
737 }
738
739 Matrice& Matrice::ConcatenateRow(const double* RowData)
740 {
741 Matrice temp( NbRows+1, NbCols,0.0);
742 int i;
743 for(i=0; i<NbRows; ++i) {
744 for(int j=0; j<NbCols; ++j) {
745 temp.Data[i][j] = Data[i][j];
746 }
747 }
748 int j;
749 for(i=NbRows, j=0; j<NbCols; ++j) {
750 temp.Data[i][j] = RowData[j];
751 }
752
753 *this = temp;
754
755 return *this;
756 }
757
758 Matrice& Matrice::ConcatenateColumn(const double* ColumnData)
759 {
760
761 Matrice temp( NbRows, NbCols+1,0.0);
762 int i;
763 for(i=0; i<NbRows; ++i) {
764 for(int j=0; j<NbCols; ++j) {
765 temp.Data[i][j] = Data[i][j];
766 }
767 }
768 int j;
769 for(j=NbCols, i=0; i<NbRows; ++i) {
770 temp.Data[i][j] = ColumnData[i];
771 }
772
773 *this = temp;
774
775 return *this;
776 }
777
778 Matrice& Matrice::Add_Row(const vecteur& RowData, const int RowSpot)
779 {
780 int New_NbRows=NbRows;
781 int New_NbCols=NbCols;
782
783 Matrice temp( NbRows+1, NbCols,0.0);
784
785 //if(NbRows==0)
786 New_NbRows=NbRows+1;
787
788
789
790 for(int i=0, k=0; i<New_NbRows; ++i, ++k) {
791 if(i == RowSpot) {
792 for(int j=0; j<New_NbCols; ++j) {
793 temp.Data[i][j] += RowData[j];
794 }
795 ++k;
796 }
797
798 if(k<New_NbRows){
799 for(int j=0; j<New_NbCols; ++j)
800 temp.Data[k][j] += Data[i][j];
801
802 }
803 }
804
805 *this = temp;
806
807 return *this;
808 }
809
810
811 Matrice& Matrice::Add_Column(const vecteur& ColumnData, const int ColumnSpot)
812 {
813 int New_NbRows=NbRows;
814 int New_NbCols=NbCols;
815
816 Matrice temp(NbRows, NbCols+1,0.0);
817
818 // if(NbCols==0)
819 New_NbCols=NbCols+1;
820
821 for(int i=0; i<NbRows; ++i) {
822 for(int j=0, k=0; j<New_NbCols; ++j, ++k) {
823 if(j == ColumnSpot) {
824 temp.Data[i][k] += ColumnData[i];
825 ++k;
826 }
827 if(k<New_NbCols)
828 temp.Data[i][k] += Data[i][j];
829 }
830 }
831
832 *this = temp;
833
834 return *this;
835 }
836
837
838 vecteur Matrice::Get_Row(int RowSpot)
839 {
840 vecteur RowData(NbCols);
841
842 for(int j=0; j<NbCols; ++j) {
843 // if(i == RowSpot) {
844 // for(int j=0; j<NbCols; ++j) {
845 RowData[j]=Data[RowSpot][j];
846 // }
847 // break;
848 // }
849
850 }
851 return RowData;
852 }
853
854 vecteur Matrice::Get_Col(int ColSpot)
855 {
856 vecteur ColData(NbRows);
857
858 for(int i=0; i<NbRows; ++i) {
859 //for(int j=0; j<NbCols; ++j) {
860 //if(j == ColSpot) {
861 ColData[i]=Data[i][ColSpot];
862 // }
863 // break;
864 // }
865
866 }
867 return ColData;
868
869 }
870
871
872 Matrice& Matrice::RemoveRow(const int Row)
873 {
874 Matrice temp(NbRows-1, NbCols,0.0);
875
876 for(int i=0, k=0; i<NbRows; ++i, ++k) {
877 if(i == Row) ++i;
878 for(int j=0; j<NbCols; ++j) {
879 temp.Data[k][j] = Data[i][j];
880 }
881 }
882
883 *this = temp;
884
885 return *this;
886 }
887
888 Matrice& Matrice::RemoveColumn(const int Column)
889 {
890 Matrice temp(0.0, NbRows, NbCols-1);
891
892 for(int i=0; i<NbRows; ++i) {
893 for(int j=0, l=0; j<NbCols; ++j, ++l) {
894 if(j == Column) ++j;
895 temp.Data[i][l] = Data[i][j];
896 }
897 }
898
899 *this = temp;
900
901 return *this;
902 }
903
904
905 Matrice& Matrice::Get_Covariant() const
906 {
907 Matrice temp;
908
909 temp = *this;
910 temp.Transpose();
911
912 return(*this * temp);
913 }
914
915 Matrice& Matrice::Set_Covariant()
916 {
917 *this = this->Get_Covariant();
918
919 return *this;
920 }
921
922
923 Matrice& Matrice::IdentityMatrice(int Diagonal)
924 {
925 static Matrice temp;
926 temp = Matrice( Diagonal, Diagonal,0.0);
927
928 for(int q=0; q<Diagonal; ++q) {
929 temp.Data[q][q] = 1;
930 }
931
932 return temp;
933 }
934
935
936 void Matrice:: affiche()
937 {
938 cout << "[";
939 for(int i=0; i<NbRows; ++i) {
940 for(int j=0; j<NbCols; ++j) {
941 cout << "[" << Data[i][j] << "]";
942 }
943 if(NbRows-i-1) cout << "\n ";
944 }
945 cout << "]\n";
946
947
948 }
949
950
951
952 Matrice& Matrice::Get_LU(vecteur& perm,double&d)
953 {
954 static const double petit=double(1.0e-20);
955 int taille=NbRows;
956 int imax,i,j,k;
957 vector<double> v1(taille);
958 double maxi,temp,som,dum;
959 d=1;
960 if(taille&&taille==NbCols)
961 {
962 for(i=0;i<taille;i++)
963 {
964 maxi=double(0);
965 for(j=0;j<taille;j++)
966 if((temp=(double)fabs(Data[i][j]))>maxi)
967 maxi=temp;
968 if(maxi==0)
969 {
970 this->Clear();
971 return *this;
972 }
973 v1[i]=1./maxi;
974 }
975
976 for(j=0;j<taille;j++)
977 {
978 //cherche les colonnnes de lamatrice L (triangulaire inferieur
979 for(i=0;i<j;i++)
980 {
981 som=Data[i][j];
982 for(k=0;k<i;k++)
983 som-=Data[i][k]*Data[k][j];
984 Data[i][j]=som;
985 }
986 maxi=0;
987
988 //cherches les lignes de la matrice triangulaire superieur
989 for(i=j;i<taille;i++)
990 {
991 som=Data[i][j];
992 for(k=0;k<j;k++)
993 som-=Data[i][k]*Data[k][j];
994 Data[i][j]=som;
995 if((dum=v1[i]*fabs(som))>=maxi)
996 {
997 maxi=dum;
998 imax=i;
999 }
1000 }
1001 //permuation des lignes imax et j
1002 if(j!=imax)
1003 {
1004 for(k=0;k<taille;k++)
1005 {
1006 dum=Data[imax][k];
1007 Data[imax][k]=Data[j][k];
1008 Data[j][k]=dum;
1009 }
1010 d*=-1;
1011 v1[imax]=v1[j];
1012 }
1013 // stoker dans perm la permuattion
1014 perm[j]=imax;
1015 // si le pivot est null
1016 if(Data[j][j]==0.0)
1017 Data[j][j]=petit;
1018 if(j!=taille-1)
1019 {
1020 dum=1./Data[j][j];
1021 for(i=j+1;i<taille;i++)
1022 Data[i][j]*=dum;
1023 }
1024 }
1025 }
1026 return *this;
1027 }
1028
1029
1030 vecteur Matrice::Get_Solution(vecteur& perm,const vecteur&b)
1031 {
1032 vecteur r1;
1033 int taille=NbRows,i,j,ip,ii=-1;
1034 double som;
1035
1036 if(taille&&taille==NbCols&&taille==(int)b.getDim())// Get_Size())
1037 {
1038 r1=b;
1039 for(i=0;i<taille;i++)
1040 {
1041 ip=perm[i];
1042 som=r1[ip];
1043 r1[ip]=r1[i];
1044 if(ii>=0)
1045 for(j=ii;j<=i-1;j++)
1046 som-=Data[i][j]*r1[j];
1047 else if(som) ii=i;
1048 r1[i]=som;
1049 }
1050 for(i=taille-1;i>=0;i--)
1051 {
1052 som=r1[i];
1053 for(j=i+1;j<taille;j++)
1054 som-=Data[i][j]*r1[j];
1055 r1[i]=som/Data[i][i];
1056 }
1057 }
1058 return r1;
1059 }
1060
1061
1062
1063 Matrice Matrice ::Get_Inverse(Matrice& K)
1064 {
1065 double d;
1066 Matrice inv(NbRows,0,0.);
1067 vecteur perm(NbRows);
1068 K= K. Get_LU(perm,d);
1069
1070 for(int i=0;i<NbRows;i++)
1071 {
1072 vecteur ei(NbRows);
1073 ei[i]=1.;
1074 if (i>0) ei[i-1]=0.;
1075 vecteur& xi=K.Get_Solution(perm,ei);
1076 inv.Add_Column(xi,i);
1077
1078 }
1079 return inv;
1080
1081 }
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092 void Matrice::Get_Solution(vecteur& perm, vecteur& b, vecteur& x)
1093 {
1094
1095 int taille=NbRows,i,j,ip,ii=-1;
1096 double som;
1097
1098
1099 for(i=0;i<taille;i++) x[i]=b[i];
1100
1101 for(i=0;i<taille;i++)
1102 {
1103 ip=perm[i];
1104 som=x[ip];
1105 x[ip]=x[i];
1106 if(ii>=0)
1107 for(j=ii;j<=i-1;j++)
1108 som-=Data[i][j]*x[j];
1109 else if(som) ii=i;
1110 x[i]=som;
1111 }
1112 for(i=taille-1;i>=0;i--)
1113 {
1114 som=x[i];
1115 for(j=i+1;j<taille;j++)
1116 som-=Data[i][j]*x[j];
1117 x[i]=som/Data[i][i];
1118 }
1119
1120
1121 }
1122
1123
1124 vecteur Matrice:: Multiply(const vecteur& v)
1125 {
1126 assert (NbCols==v.getDim());
1127 vecteur z(NbRows);
1128 for(int i=0; i<NbRows; ++i) {
1129 for(int j=0; j<NbCols; ++j) {
1130 z[i] += Data[i][j]*v[j];
1131 }
1132 }
1133
1134 return z;
1135
1136 }
1137
1138
1139
1140 //extraire une matrice de la ligne 1 jusqua iigne 2 et de la colone 1 jusqua la colonne2
1141 /*
1142 Matrice Matrice ::Get_SubMatrice(int ln1,int ln2, int cln1,int cln2)
1143 {
1144 int nb_lgne=ln2-ln1+1;
1145 int nb_colne=cln2-cln1+1;
1146 Matrice m(nb_lgne, nb_colne,0.);
1147
1148 for(int i=ln1;i<=ln2;i++)
1149 {
1150 int I=i-ln1;
1151 for(int j=cln1;j<=cln2;j++)
1152 {
1153 int J=j-cln1;
1154 m(I,J)=Data[i][j] ;
1155 }
1156
1157
1158 }
1159
1160 return m;
1161 }
1162 */
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175 std::ostream& operator <<(std::ostream& os, const Matrice& obj)
1176 {
1177 os<<"------------------------------------------------------------------------------------"<<endl;
1178 os<<"# MATRICE: DIM1= "<< obj.Get_NbRows()<<" , DIM2: "<<obj.Get_NbCols()<<" #"<<endl;
1179 os<<"------------------------------------------------------------------------------------"<<endl;
1180 for(int i=0;i<obj.Get_NbRows();i++) {
1181 for(int j=0;j<obj.Get_NbCols();j++) {
1182 os<<setw(18)<<obj(i,j);
1183 }
1184 os<<endl;
1185 }
1186 os<<"------------------------------------------------------------------------------------"<<endl;
1187 return os;
1188 /*
1189 for(int i=0; i<obj.Get_NbRows(); ++i) {
1190 for(int j=0; j<obj.Get_NbCols(); ++j) {
1191 if (i==0)ostr<<"["<<setw(9);
1192 ostr <<setw(10)<<obj[i][j];
1193 }
1194 if(obj.Get_NbRows()-i-1) ostr << "\n ";
1195 }
1196 ostr<<"]"<<endl;
1197
1198 return ostr; */
1199 }
1200
1201 vecteur operator * (const Matrice& m,const vecteur& v)
1202 {
1203 assert (m.Get_NbCols()==v.getDim()) ;
1204 vecteur z(m.Get_NbRows());
1205 for(int i=0; i<m.Get_NbRows(); ++i) {
1206 for(int j=0; j<m.Get_NbCols(); ++j) {
1207 z[i] += m(i,j)*v[j];
1208 }
1209 }
1210
1211 return z;
1212
1213 }
1214 /*
1215 void Matrice::convert_Mat(const Matrice &mat,double valmin,MatCr &mat1)
1216 {
1217 int i,j,k,taille=mat.Get_NbRows();
1218
1219 if(mat.Get_NbCols()==taille)
1220 {
1221 mat1._index.resize((unsigned int)(taille+1));
1222 mat1._values.resize((unsigned int)(taille+1));
1223 for(i=0;i<taille;mat1._values[(unsigned int)i]=mat(i,i),i++);
1224 mat1._index[0]=taille+1;
1225 k=taille;
1226 for(i=0;i<taille;i++)
1227 {
1228 for(j=0;j<taille;j++)
1229 if(fabs(mat(i,j))>valmin&&i!=j)
1230 {
1231 k++;
1232 mat1._index.push_back(j);
1233 mat1._values.push_back(mat(i,j));
1234 }
1235 mat1._index[i+1]=k+1;
1236 }
1237 }
1238 else
1239 {
1240 mat1._index.resize(0);
1241 mat1._values.resize(0);
1242 }
1243 }
1244
1245
1246 */
1247
1248