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 |
|