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,
53 Standard_RangeError_Raise_if ((LowerRow > UpperRow)
54 || (LowerCol > UpperCol), "math_Matrix() - invalid dimensions");
57 math_Matrix::math_Matrix (const Standard_Integer LowerRow,
58 const Standard_Integer UpperRow,
59 const Standard_Integer LowerCol,
60 const Standard_Integer UpperCol,
61 const Standard_Real InitialValue):
63 LowerRowIndex(LowerRow),
64 UpperRowIndex(UpperRow),
65 LowerColIndex(LowerCol),
66 UpperColIndex(UpperCol),
67 Array(LowerRow, UpperRow,
70 Standard_RangeError_Raise_if ((LowerRow > UpperRow)
71 || (LowerCol > UpperCol), "math_Matrix() - invalid dimensions");
72 Array.Init(InitialValue);
75 math_Matrix::math_Matrix (const Standard_Address Tab,
76 const Standard_Integer LowerRow,
77 const Standard_Integer UpperRow,
78 const Standard_Integer LowerCol,
79 const Standard_Integer UpperCol) :
81 LowerRowIndex(LowerRow),
82 UpperRowIndex(UpperRow),
83 LowerColIndex(LowerCol),
84 UpperColIndex(UpperCol),
85 Array(Tab, LowerRow, UpperRow, LowerCol, UpperCol)
87 Standard_RangeError_Raise_if ((LowerRow > UpperRow)
88 || (LowerCol > UpperCol), "math_Matrix() - invalid dimensions");
91 void math_Matrix::Init(const Standard_Real InitialValue)
93 Array.Init(InitialValue);
96 math_Matrix::math_Matrix (const math_Matrix& Other):
98 LowerRowIndex(Other.LowerRow()),
99 UpperRowIndex(Other.UpperRow()),
100 LowerColIndex(Other.LowerCol()),
101 UpperColIndex(Other.UpperCol()),
108 math_Matrix math_Matrix::Divided (const Standard_Real Right) const
110 Standard_DivideByZero_Raise_if (Abs(Right) <= RealEpsilon(),
111 "math_Matrix::Divided() - zero divisor");
112 math_Matrix temp = Multiplied(1./Right);
117 Standard_Real math_Matrix::Determinant() const
119 math_Gauss Sol(*this);
122 return Sol.Determinant();
129 void math_Matrix::Transpose()
131 math_NotSquare_Raise_if (RowNumber() != ColNumber(),
132 "math_Matrix::Transpose() - matrix is not square");
134 Standard_Integer Row = LowerRowIndex;
135 Standard_Integer Col = LowerColIndex;
136 SetLowerCol(LowerRowIndex);
138 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
139 for(Standard_Integer J = I; J <= UpperColIndex; J++) {
141 Array(I, J) = Array(J, I);
149 math_Matrix math_Matrix::Transposed() const
151 math_Matrix Result(LowerColIndex, UpperColIndex,
152 LowerRowIndex, UpperRowIndex);
154 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
155 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
156 Result.Array(J, I) = Array(I, J);
162 void math_Matrix::Invert()
164 math_NotSquare_Raise_if (RowNumber() != ColNumber(),
165 "math_Matrix::Transpose() - matrix is not square");
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(),
222 "math_Matrix::Divide() - zero divisor");
224 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
225 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
226 Array(I, J) = Array(I, J) / Right;
231 void math_Matrix::Add (const math_Matrix& Right)
233 Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
234 || (ColNumber() != Right.ColNumber()),
235 "math_Matrix::Add() - input matrix has different dimensions");
237 Standard_Integer I2 = Right.LowerRowIndex;
238 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
239 Standard_Integer J2 = Right.LowerColIndex;
240 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
241 Array(I, J) = Array(I, J) + Right.Array(I2, J2);
248 void math_Matrix::Subtract (const math_Matrix& Right)
250 Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
251 || (ColNumber() != Right.ColNumber()),
252 "math_Matrix::Subtract() - input matrix has different dimensions");
254 Standard_Integer I2 = Right.LowerRowIndex;
255 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
256 Standard_Integer J2 = Right.LowerColIndex;
257 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
258 Array(I, J) = Array(I, J) - Right.Array(I2, J2);
265 void math_Matrix::Set(const Standard_Integer I1,const Standard_Integer I2,
266 const Standard_Integer J1,const Standard_Integer J2,
267 const math_Matrix& M)
270 Standard_DimensionError_Raise_if((I1 < LowerRowIndex) ||
271 (I2 > UpperRowIndex) ||
272 (J1 < LowerColIndex) ||
273 (J2 > UpperColIndex) ||
274 (I1 > I2) || (J1 > J2) ||
275 (I2-I1+1 != M.RowNumber()) ||
276 (J2-J1+1 != M.ColNumber()),
277 "math_Matrix::Set() - invalid indices");
279 Standard_Integer II = M.LowerRow();
280 for(Standard_Integer I = I1; I <= I2; I++) {
281 Standard_Integer JJ = M.LowerCol();
282 for(Standard_Integer J = J1; J <= J2; J++) {
283 Array(I, J) = M.Array(II, JJ);
290 void math_Matrix::SetRow (const Standard_Integer Row,
291 const math_Vector& V)
293 Standard_RangeError_Raise_if ((Row < LowerRowIndex)
294 || (Row > UpperRowIndex),
295 "math_Matrix::SetRow() - invalid index");
297 Standard_DimensionError_Raise_if (ColNumber() != V.Length(),
298 "math_Matrix::SetRow() - input vector has wrong dimensions");
300 Standard_Integer I = V.Lower();
301 for(Standard_Integer Index = LowerColIndex; Index <= UpperColIndex; Index++) {
302 Array(Row, Index) = V.Array(I);
307 void math_Matrix::SetCol (const Standard_Integer Col,
308 const math_Vector& V)
310 Standard_RangeError_Raise_if ((Col < LowerColIndex)
311 || (Col > UpperColIndex),
312 "math_Matrix::SetCol() - invalid index");
314 Standard_DimensionError_Raise_if (RowNumber() != V.Length(),
315 "math_Matrix::SetCol() - input vector has wrong dimensions");
317 Standard_Integer I = V.Lower();
318 for(Standard_Integer Index = LowerRowIndex; Index <= UpperRowIndex; Index++) {
319 Array(Index, Col) = V.Array(I);
324 void math_Matrix::SetDiag(const Standard_Real Value)
326 math_NotSquare_Raise_if (RowNumber() != ColNumber(),
327 "math_Matrix::SetDiag() - matrix is not square");
329 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
333 math_Vector math_Matrix::Row (const Standard_Integer Row) const
336 math_Vector Result(LowerColIndex, UpperColIndex);
338 for(Standard_Integer Index = LowerColIndex; Index <= UpperColIndex; Index++) {
339 Result.Array(Index) = Array(Row, Index);
344 math_Vector math_Matrix::Col (const Standard_Integer Col) const
347 math_Vector Result(LowerRowIndex, UpperRowIndex);
349 for(Standard_Integer Index = LowerRowIndex; Index <= UpperRowIndex; Index++) {
350 Result.Array(Index) = Array(Index, Col);
355 void math_Matrix::SwapRow(const Standard_Integer Row1,
356 const Standard_Integer Row2)
358 Standard_RangeError_Raise_if ((Row1 < LowerRowIndex)
359 || (Row1 > UpperRowIndex)
360 || (Row2 < LowerRowIndex)
361 || (Row2 > UpperRowIndex),
362 "math_Matrix::SetCol() - invalid indices");
364 math_Vector V1 = Row(Row1);
365 math_Vector V2 = Row(Row2);
370 void math_Matrix::SwapCol(const Standard_Integer Col1,
371 const Standard_Integer Col2)
373 Standard_RangeError_Raise_if ((Col1 < LowerColIndex)
374 || (Col1 > UpperColIndex)
375 || (Col2 < LowerColIndex)
376 || (Col2 > UpperColIndex),
377 "math_Matrix::SetCol() - invalid indices");
379 math_Vector V1 = Col(Col1);
380 math_Vector V2 = Col(Col2);
387 math_Matrix math_Matrix::Multiplied (const math_Matrix& Right) const
389 Standard_DimensionError_Raise_if (ColNumber() != Right.RowNumber(),
390 "math_Matrix::Multiplied() - matrices have incompatible dimensions");
392 math_Matrix Result(LowerRowIndex, UpperRowIndex,
393 Right.LowerColIndex, Right.UpperColIndex);
396 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
397 for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
399 Standard_Integer I2 = Right.LowerRowIndex;
400 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
401 Som = Som + Array(I, J) * Right.Array(I2, J2);
404 Result.Array(I, J2) = Som;
410 math_Matrix math_Matrix::TMultiply (const math_Matrix& Right) const
412 Standard_DimensionError_Raise_if (RowNumber() != Right.RowNumber(),
413 "math_Matrix::TMultiply() - matrices have incompatible dimensions");
415 math_Matrix Result(LowerColIndex, UpperColIndex,
416 Right.LowerColIndex, Right.UpperColIndex);
419 for(Standard_Integer I = LowerColIndex; I <= UpperColIndex; I++) {
420 for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
422 Standard_Integer I2 = Right.LowerRowIndex;
423 for(Standard_Integer J = LowerRowIndex; J <= UpperRowIndex; J++) {
424 Som = Som + Array(J, I) * Right.Array(I2, J2);
427 Result.Array(I, J2) = Som;
433 math_Matrix math_Matrix::Added (const math_Matrix& Right) const
435 Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
436 || (ColNumber() != Right.ColNumber()),
437 "math_Matrix::Added() - input matrix has different dimensions");
439 math_Matrix Result(LowerRowIndex, UpperRowIndex,
440 LowerColIndex, UpperColIndex);
442 Standard_Integer I2 = Right.LowerRowIndex;
443 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
444 Standard_Integer J2 = Right.LowerColIndex;
445 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
446 Result.Array(I, J) = Array(I, J) + Right.Array(I2, J2);
454 math_Matrix math_Matrix::Opposite ()
457 math_Matrix Result(LowerRowIndex, UpperRowIndex,
458 LowerColIndex, UpperColIndex);
460 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
461 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
462 Result.Array(I, J) = - Array(I, J);
468 math_Matrix math_Matrix::Subtracted (const math_Matrix& Right) const
470 Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
471 || (ColNumber() != Right.ColNumber()),
472 "math_Matrix::Subtracted() - input matrix has different dimensions");
474 math_Matrix Result(LowerRowIndex, UpperRowIndex,
475 LowerColIndex, UpperColIndex);
477 Standard_Integer I2 = Right.LowerRowIndex;
478 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
479 Standard_Integer J2 = Right.LowerColIndex;
480 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
481 Result.Array(I, J) = Array(I, J) - Right.Array(I2, J2);
489 void math_Matrix::Multiply(const math_Vector& Left,
490 const math_Vector& Right)
492 Standard_DimensionError_Raise_if ((RowNumber() != Left.Length())
493 || (ColNumber() != Right.Length()),
494 "math_Matrix::Multiply() - input vectors have incompatible dimensions");
496 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
497 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
498 Array(I, J) = Left.Array(I) * Right.Array(J);
503 void math_Matrix::Multiply(const math_Matrix& Left,
504 const math_Matrix& Right)
506 Standard_DimensionError_Raise_if ((Left.ColNumber() != Right.RowNumber())
507 || (RowNumber() != Left.RowNumber())
508 || (ColNumber() != Right.ColNumber()),
509 "math_Matrix::Multiply() - matrices have incompatible dimensions");
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)
534 Standard_DimensionError_Raise_if ((TLeft.RowNumber() != Right.RowNumber())
535 || (RowNumber() != TLeft.ColNumber())
536 || (ColNumber() != Right.ColNumber()),
537 "math_Matrix::TMultiply() - matrices have incompatible dimensions");
540 Standard_Integer I1 = TLeft.LowerColIndex;
541 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
542 Standard_Integer J2 = Right.LowerColIndex;
543 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
545 Standard_Integer J1 = TLeft.LowerRowIndex;
546 Standard_Integer I2 = Right.LowerRowIndex;
547 for(Standard_Integer K = TLeft.LowerRowIndex; K <= TLeft.UpperRowIndex; K++) {
548 Som = Som + TLeft.Array(J1, I1) * Right.Array(I2, J2);
559 void math_Matrix::Add (const math_Matrix& Left,
560 const math_Matrix& Right)
562 Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
563 || (ColNumber() != Right.ColNumber())
564 || (Right.RowNumber() != Left.RowNumber())
565 || (Right.ColNumber() != Left.ColNumber()),
566 "math_Matrix::Add() - matrices have incompatible dimensions");
568 Standard_Integer I1 = Left.LowerRowIndex;
569 Standard_Integer I2 = Right.LowerRowIndex;
570 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
571 Standard_Integer J1 = Left.LowerColIndex;
572 Standard_Integer J2 = Right.LowerColIndex;
573 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
574 Array(I, J) = Left.Array(I1, J1) + Right.Array(I2, J2);
583 void math_Matrix::Subtract(const math_Matrix& Left,
584 const math_Matrix& Right)
586 Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
587 || (ColNumber() != Right.ColNumber())
588 || (Right.RowNumber() != Left.RowNumber())
589 || (Right.ColNumber() != Left.ColNumber()),
590 "math_Matrix::Subtract() - matrices have incompatible dimensions");
592 Standard_Integer I1 = Left.LowerRowIndex;
593 Standard_Integer I2 = Right.LowerRowIndex;
594 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
595 Standard_Integer J1 = Left.LowerColIndex;
596 Standard_Integer J2 = Right.LowerColIndex;
597 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
598 Array(I, J) = Left.Array(I1, J1) - Right.Array(I2, J2);
608 void math_Matrix::Multiply(const math_Matrix& Right)
610 Standard_DimensionError_Raise_if (ColNumber() != Right.RowNumber(),
611 "math_Matrix::Multiply() - input matrix has incompatible dimensions");
614 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
615 for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
617 Standard_Integer I2 = Right.LowerRowIndex;
618 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
619 Som = Som + Array(I, J) * Right.Array(I2, J2);
629 math_Vector math_Matrix::Multiplied(const math_Vector& Right)const
631 Standard_DimensionError_Raise_if (ColNumber() != Right.Length(),
632 "math_Matrix::Multiplied() - input vector has incompatible dimensions");
634 math_Vector Result(LowerRowIndex, UpperRowIndex);
636 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
637 Result.Array(I) = 0.0;
638 Standard_Integer II = Right.Lower();
639 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
640 Result.Array(I) = Result.Array(I) + Array(I, J) * Right.Array(II);
647 math_Matrix& math_Matrix::Initialized(const math_Matrix& Other)
649 Standard_DimensionError_Raise_if ((RowNumber() != Other.RowNumber())
650 || (ColNumber() != Other.ColNumber()),
651 "math_Matrix::Initialized() - input matrix has different dimensions");
653 (Other.Array).Copy(Array);
660 void math_Matrix::Dump(Standard_OStream& o)const
663 o << "math_Matrix of RowNumber = " << RowNumber();
664 o << " and ColNumber = " << ColNumber() << "\n";
666 for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
667 for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
668 o << "math_Matrix ( " << I << ", " << J << " ) = ";
669 o << Array(I, J) << "\n";