1 // Copyright (c) 1997-1999 Matra Datavision
2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
4 // This file is part of Open CASCADE Technology software library.
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License version 2.1 as published
8 // by the Free Software Foundation, with special exception defined in the file
9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
10 // distribution for complete text of the license and disclaimer of any warranty.
12 // Alternatively, this file may be used under the terms of Open CASCADE
13 // commercial license or contractual agreement.
16 #include <math_Gauss.hxx>
17 #include <math_Matrix.hxx>
18 #include <math_NotSquare.hxx>
19 #include <math_SingularMatrix.hxx>
20 #include <math_Vector.hxx>
21 #include <Standard_DimensionError.hxx>
22 #include <Standard_DivideByZero.hxx>
23 #include <Standard_RangeError.hxx>
25 void math_Matrix::SetLowerRow(const Standard_Integer LowerRow) {
27 Array.SetLowerRow(LowerRow);
28 Standard_Integer Rows = RowNumber();
29 LowerRowIndex = LowerRow;
30 UpperRowIndex = LowerRowIndex + Rows - 1;
33 void math_Matrix::SetLowerCol(const Standard_Integer LowerCol) {
35 Array.SetLowerCol(LowerCol);
36 Standard_Integer Cols = ColNumber();
37 LowerColIndex = LowerCol;
38 UpperColIndex = LowerColIndex + Cols - 1;
41 math_Matrix::math_Matrix (const Standard_Integer LowerRow,
42 const Standard_Integer UpperRow,
43 const Standard_Integer LowerCol,
44 const Standard_Integer UpperCol):
46 LowerRowIndex(LowerRow),
47 UpperRowIndex(UpperRow),
48 LowerColIndex(LowerCol),
49 UpperColIndex(UpperCol),
50 Array(LowerRow, UpperRow,
54 Standard_RangeError_Raise_if((LowerRow > UpperRow) ||
55 (LowerCol > UpperCol), "");
58 math_Matrix::math_Matrix (const Standard_Integer LowerRow,
59 const Standard_Integer UpperRow,
60 const Standard_Integer LowerCol,
61 const Standard_Integer UpperCol,
62 const Standard_Real InitialValue):
64 LowerRowIndex(LowerRow),
65 UpperRowIndex(UpperRow),
66 LowerColIndex(LowerCol),
67 UpperColIndex(UpperCol),
68 Array(LowerRow, UpperRow,
72 Standard_RangeError_Raise_if((LowerRow > UpperRow) ||
73 (LowerCol > UpperCol), "");
74 Array.Init(InitialValue);
77 math_Matrix::math_Matrix (const Standard_Address Tab,
78 const Standard_Integer LowerRow,
79 const Standard_Integer UpperRow,
80 const Standard_Integer LowerCol,
81 const Standard_Integer UpperCol) :
83 LowerRowIndex(LowerRow),
84 UpperRowIndex(UpperRow),
85 LowerColIndex(LowerCol),
86 UpperColIndex(UpperCol),
87 Array(Tab, LowerRow, UpperRow, LowerCol, UpperCol)
90 Standard_RangeError_Raise_if((LowerRow > UpperRow) ||
91 (LowerCol > UpperCol), "");
94 void math_Matrix::Init(const Standard_Real InitialValue)
96 Array.Init(InitialValue);
99 math_Matrix::math_Matrix (const math_Matrix& Other):
101 LowerRowIndex(Other.LowerRow()),
102 UpperRowIndex(Other.UpperRow()),
103 LowerColIndex(Other.LowerCol()),
104 UpperColIndex(Other.UpperCol()),
111 math_Matrix math_Matrix::Divided (const Standard_Real Right) const
113 Standard_DivideByZero_Raise_if(Abs(Right) <= RealEpsilon(), "");
114 math_Matrix temp = Multiplied(1./Right);
119 Standard_Real math_Matrix::Determinant() const
121 math_Gauss Sol(*this);
124 return Sol.Determinant();
131 void math_Matrix::Transpose()
133 math_NotSquare_Raise_if(RowNumber() != ColNumber(), "");
135 Standard_Integer Row = LowerRowIndex;
136 Standard_Integer Col = LowerColIndex;
137 SetLowerCol(LowerRowIndex);
139 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
140 for(Standard_Integer J = I; J <= UpperColIndex; J++) {
142 Array(I, J) = Array(J, I);
150 math_Matrix math_Matrix::Transposed() const
152 math_Matrix Result(LowerColIndex, UpperColIndex,
153 LowerRowIndex, UpperRowIndex);
155 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
156 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
157 Result.Array(J, I) = Array(I, J);
163 void math_Matrix::Invert()
165 math_NotSquare_Raise_if(RowNumber() != ColNumber(), "");
167 math_Gauss Sol(*this);
172 throw math_SingularMatrix(); // SingularMatrix Exception;
176 math_Matrix math_Matrix::Inverse() const
179 math_Matrix Result = *this;
184 void math_Matrix::Multiply (const Standard_Real Right)
186 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
187 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
188 Array(I, J) = Array(I, J) * Right;
193 math_Matrix math_Matrix::Multiplied (const Standard_Real Right) const
195 math_Matrix Result(LowerRowIndex, UpperRowIndex,
196 LowerColIndex, UpperColIndex);
197 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
198 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
199 Result.Array(I, J) = Array(I, J) * Right;
205 math_Matrix math_Matrix::TMultiplied (const Standard_Real Right) const
207 math_Matrix Result(LowerRowIndex, UpperRowIndex,
208 LowerColIndex, UpperColIndex);
209 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
210 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
211 Result.Array(I, J) = Array(I, J) * Right;
219 void math_Matrix::Divide (const Standard_Real Right)
221 Standard_DivideByZero_Raise_if(Abs(Right) <= RealEpsilon(), "");
223 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
224 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
225 Array(I, J) = Array(I, J) / Right;
230 void math_Matrix::Add (const math_Matrix& Right)
232 Standard_DimensionError_Raise_if((RowNumber() != Right.RowNumber()) ||
233 (ColNumber() != Right.ColNumber()),
236 Standard_Integer I2 = Right.LowerRowIndex;
237 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
238 Standard_Integer J2 = Right.LowerColIndex;
239 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
240 Array(I, J) = Array(I, J) + Right.Array(I2, J2);
247 void math_Matrix::Subtract (const math_Matrix& Right)
249 Standard_DimensionError_Raise_if((RowNumber() != Right.RowNumber()) ||
250 (ColNumber() != Right.ColNumber()),
253 Standard_Integer I2 = Right.LowerRowIndex;
254 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
255 Standard_Integer J2 = Right.LowerColIndex;
256 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
257 Array(I, J) = Array(I, J) - Right.Array(I2, J2);
264 void math_Matrix::Set(const Standard_Integer I1,const Standard_Integer I2,
265 const Standard_Integer J1,const Standard_Integer J2,
266 const math_Matrix& M)
269 Standard_DimensionError_Raise_if((I1 < LowerRowIndex) ||
270 (I2 > UpperRowIndex) ||
271 (J1 < LowerColIndex) ||
272 (J2 > UpperColIndex) ||
273 (I1 > I2) || (J1 > J2) ||
274 (I2-I1+1 != M.RowNumber()) ||
275 (J2-J1+1 != M.ColNumber()), "");
277 Standard_Integer II = M.LowerRow();
278 for(Standard_Integer I = I1; I <= I2; I++) {
279 Standard_Integer JJ = M.LowerCol();
280 for(Standard_Integer J = J1; J <= J2; J++) {
281 Array(I, J) = M.Array(II, JJ);
288 void math_Matrix::SetRow (const Standard_Integer Row,
289 const math_Vector& V)
292 Standard_RangeError_Raise_if((Row < LowerRowIndex) ||
293 (Row > UpperRowIndex) , "");
295 Standard_DimensionError_Raise_if(ColNumber() != V.Length(), "");
297 Standard_Integer I = V.LowerIndex;
298 for(Standard_Integer Index = LowerColIndex; Index <= UpperColIndex; Index++) {
299 Array(Row, Index) = V.Array(I);
304 void math_Matrix::SetCol (const Standard_Integer Col,
305 const math_Vector& V)
308 Standard_RangeError_Raise_if((Col < LowerColIndex) ||
309 (Col > UpperColIndex) , "");
311 Standard_DimensionError_Raise_if(RowNumber() != V.Length(), "");
313 Standard_Integer I = V.LowerIndex;
314 for(Standard_Integer Index = LowerRowIndex; Index <= UpperRowIndex; Index++) {
315 Array(Index, Col) = V.Array(I);
320 void math_Matrix::SetDiag(const Standard_Real Value)
323 math_NotSquare_Raise_if(RowNumber() != ColNumber(), "");
325 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
329 math_Vector math_Matrix::Row (const Standard_Integer Row) const
332 math_Vector Result(LowerColIndex, UpperColIndex);
334 for(Standard_Integer Index = LowerColIndex; Index <= UpperColIndex; Index++) {
335 Result.Array(Index) = Array(Row, Index);
340 math_Vector math_Matrix::Col (const Standard_Integer Col) const
343 math_Vector Result(LowerRowIndex, UpperRowIndex);
345 for(Standard_Integer Index = LowerRowIndex; Index <= UpperRowIndex; Index++) {
346 Result.Array(Index) = Array(Index, Col);
351 void math_Matrix::SwapRow(const Standard_Integer Row1,
352 const Standard_Integer Row2)
355 Standard_RangeError_Raise_if((Row1 < LowerRowIndex) ||
356 (Row1 > UpperRowIndex) ||
357 (Row2 < LowerRowIndex) ||
358 (Row2 > UpperRowIndex), "");
360 math_Vector V1 = Row(Row1);
361 math_Vector V2 = Row(Row2);
366 void math_Matrix::SwapCol(const Standard_Integer Col1,
367 const Standard_Integer Col2)
370 Standard_RangeError_Raise_if((Col1 < LowerColIndex) ||
371 (Col1 > UpperColIndex) ||
372 (Col2 < LowerColIndex) ||
373 (Col2 > UpperColIndex), "");
375 math_Vector V1 = Col(Col1);
376 math_Vector V2 = Col(Col2);
383 math_Matrix math_Matrix::Multiplied (const math_Matrix& Right) const
386 Standard_DimensionError_Raise_if(ColNumber() != Right.RowNumber(), "");
388 math_Matrix Result(LowerRowIndex, UpperRowIndex,
389 Right.LowerColIndex, Right.UpperColIndex);
392 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
393 for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
395 Standard_Integer I2 = Right.LowerRowIndex;
396 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
397 Som = Som + Array(I, J) * Right.Array(I2, J2);
400 Result.Array(I, J2) = Som;
406 math_Matrix math_Matrix::TMultiply (const math_Matrix& Right) const
409 Standard_DimensionError_Raise_if(RowNumber() != Right.RowNumber(), "");
411 math_Matrix Result(LowerColIndex, UpperColIndex,
412 Right.LowerColIndex, Right.UpperColIndex);
415 for(Standard_Integer I = LowerColIndex; I <= UpperColIndex; I++) {
416 for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
418 Standard_Integer I2 = Right.LowerRowIndex;
419 for(Standard_Integer J = LowerRowIndex; J <= UpperRowIndex; J++) {
420 Som = Som + Array(J, I) * Right.Array(I2, J2);
423 Result.Array(I, J2) = Som;
429 math_Matrix math_Matrix::Added (const math_Matrix& Right) const
432 Standard_DimensionError_Raise_if((RowNumber() != Right.RowNumber()) ||
433 (ColNumber() != Right.ColNumber()),
436 math_Matrix Result(LowerRowIndex, UpperRowIndex,
437 LowerColIndex, UpperColIndex);
439 Standard_Integer I2 = Right.LowerRowIndex;
440 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
441 Standard_Integer J2 = Right.LowerColIndex;
442 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
443 Result.Array(I, J) = Array(I, J) + Right.Array(I2, J2);
451 math_Matrix math_Matrix::Opposite ()
454 math_Matrix Result(LowerRowIndex, UpperRowIndex,
455 LowerColIndex, UpperColIndex);
457 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
458 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
459 Result.Array(I, J) = - Array(I, J);
465 math_Matrix math_Matrix::Subtracted (const math_Matrix& Right) const
468 Standard_DimensionError_Raise_if((RowNumber() != Right.RowNumber()) ||
469 (ColNumber() != Right.ColNumber()),
472 math_Matrix Result(LowerRowIndex, UpperRowIndex,
473 LowerColIndex, UpperColIndex);
475 Standard_Integer I2 = Right.LowerRowIndex;
476 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
477 Standard_Integer J2 = Right.LowerColIndex;
478 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
479 Result.Array(I, J) = Array(I, J) - Right.Array(I2, J2);
487 void math_Matrix::Multiply(const math_Vector& Left,
488 const math_Vector& Right)
491 Standard_DimensionError_Raise_if((RowNumber() != Left.Length()) ||
492 (ColNumber() != Right.Length()),
495 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
496 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
497 Array(I, J) = Left.Array(I) * Right.Array(J);
502 void math_Matrix::Multiply(const math_Matrix& Left,
503 const math_Matrix& Right)
506 Standard_DimensionError_Raise_if((Left.ColNumber() != Right.RowNumber()) ||
507 (RowNumber() != Left.RowNumber()) ||
508 (ColNumber() != Right.ColNumber()),
512 Standard_Integer I1 = Left.LowerRowIndex;
513 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
514 Standard_Integer J2 = Right.LowerColIndex;
515 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
517 Standard_Integer J1 = Left.LowerColIndex;
518 Standard_Integer I2 = Right.LowerRowIndex;
519 for(Standard_Integer K = Left.LowerColIndex; K <= Left.UpperColIndex; K++) {
520 Som = Som + Left.Array(I1, J1) * Right.Array(I2, J2);
531 void math_Matrix::TMultiply(const math_Matrix& TLeft,
532 const math_Matrix& Right)
535 Standard_DimensionError_Raise_if((TLeft.RowNumber() != Right.RowNumber()) ||
536 (RowNumber() != TLeft.ColNumber()) ||
537 (ColNumber() != Right.ColNumber()),
541 Standard_Integer I1 = TLeft.LowerColIndex;
542 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
543 Standard_Integer J2 = Right.LowerColIndex;
544 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
546 Standard_Integer J1 = TLeft.LowerRowIndex;
547 Standard_Integer I2 = Right.LowerRowIndex;
548 for(Standard_Integer K = TLeft.LowerRowIndex; K <= TLeft.UpperRowIndex; K++) {
549 Som = Som + TLeft.Array(J1, I1) * Right.Array(I2, J2);
560 void math_Matrix::Add (const math_Matrix& Left,
561 const math_Matrix& Right)
564 Standard_DimensionError_Raise_if((RowNumber() != Right.RowNumber()) ||
565 (ColNumber() != Right.ColNumber()) ||
566 (Right.RowNumber() != Left.RowNumber()) ||
567 (Right.ColNumber() != Left.ColNumber()),
570 Standard_Integer I1 = Left.LowerRowIndex;
571 Standard_Integer I2 = Right.LowerRowIndex;
572 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
573 Standard_Integer J1 = Left.LowerColIndex;
574 Standard_Integer J2 = Right.LowerColIndex;
575 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
576 Array(I, J) = Left.Array(I1, J1) + Right.Array(I2, J2);
585 void math_Matrix::Subtract(const math_Matrix& Left,
586 const math_Matrix& Right)
589 Standard_DimensionError_Raise_if((RowNumber() != Right.RowNumber()) ||
590 (ColNumber() != Right.ColNumber()) ||
591 (Right.RowNumber() != Left.RowNumber()) ||
592 (Right.ColNumber() != Left.ColNumber()),
595 Standard_Integer I1 = Left.LowerRowIndex;
596 Standard_Integer I2 = Right.LowerRowIndex;
597 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
598 Standard_Integer J1 = Left.LowerColIndex;
599 Standard_Integer J2 = Right.LowerColIndex;
600 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
601 Array(I, J) = Left.Array(I1, J1) - Right.Array(I2, J2);
611 void math_Matrix::Multiply(const math_Matrix& Right)
614 Standard_DimensionError_Raise_if(ColNumber() != Right.RowNumber(), "");
617 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
618 for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
620 Standard_Integer I2 = Right.LowerRowIndex;
621 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
622 Som = Som + Array(I, J) * Right.Array(I2, J2);
632 math_Vector math_Matrix::Multiplied(const math_Vector& Right)const
635 Standard_DimensionError_Raise_if(ColNumber() != Right.Length(), "");
637 math_Vector Result(LowerRowIndex, UpperRowIndex);
639 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
640 Result.Array(I) = 0.0;
641 Standard_Integer II = Right.LowerIndex;
642 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
643 Result.Array(I) = Result.Array(I) + Array(I, J) * Right.Array(II);
650 math_Matrix& math_Matrix::Initialized(const math_Matrix& Other)
653 Standard_DimensionError_Raise_if((RowNumber() != Other.RowNumber()) ||
654 (ColNumber() != Other.ColNumber()), "");
656 (Other.Array).Copy(Array);
663 void math_Matrix::Dump(Standard_OStream& o)const
666 o << "math_Matrix of RowNumber = " << RowNumber();
667 o << " and ColNumber = " << ColNumber() << "\n";
669 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
670 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
671 o << "math_Matrix ( " << I << ", " << J << " ) = ";
672 o << Array(I, J) << "\n";