src/math/math_Matrix.cxx

   1 // Copyright (c) 1997-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   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.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15
  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>
  24
  25 void math_Matrix::SetLowerRow(const Standard_Integer LowerRow) {
  26
  27   Array.SetLowerRow(LowerRow);
  28   Standard_Integer Rows = RowNumber();
  29   LowerRowIndex = LowerRow;
  30   UpperRowIndex = LowerRowIndex + Rows - 1;
  31 }
  32
  33 void math_Matrix::SetLowerCol(const Standard_Integer LowerCol) {
  34
  35   Array.SetLowerCol(LowerCol);
  36   Standard_Integer Cols = ColNumber();
  37   LowerColIndex = LowerCol;
  38   UpperColIndex = LowerColIndex + Cols - 1;
  39 }
  40
  41 math_Matrix::math_Matrix (const Standard_Integer LowerRow,
  42                           const Standard_Integer UpperRow,
  43                           const Standard_Integer LowerCol,
  44                           const Standard_Integer UpperCol):
  45
  46                           LowerRowIndex(LowerRow),
  47                           UpperRowIndex(UpperRow),
  48                           LowerColIndex(LowerCol),
  49                           UpperColIndex(UpperCol),
  50                           Array(LowerRow, UpperRow,
  51                                 LowerCol, UpperCol)
  52 {
  53   Standard_RangeError_Raise_if ((LowerRow > UpperRow)
  54                              || (LowerCol > UpperCol), "math_Matrix() - invalid dimensions");
  55 }
  56
  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):
  62
  63                           LowerRowIndex(LowerRow),
  64                           UpperRowIndex(UpperRow),
  65                           LowerColIndex(LowerCol),
  66                           UpperColIndex(UpperCol),
  67                           Array(LowerRow, UpperRow,
  68                                 LowerCol, UpperCol)
  69 {
  70   Standard_RangeError_Raise_if ((LowerRow > UpperRow)
  71                              || (LowerCol > UpperCol), "math_Matrix() - invalid dimensions");
  72   Array.Init(InitialValue);
  73 }
  74
  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) :
  80
  81                           LowerRowIndex(LowerRow),
  82                           UpperRowIndex(UpperRow),
  83                           LowerColIndex(LowerCol),
  84                           UpperColIndex(UpperCol),
  85                           Array(Tab, LowerRow, UpperRow, LowerCol, UpperCol)
  86 {
  87   Standard_RangeError_Raise_if ((LowerRow > UpperRow)
  88                              || (LowerCol > UpperCol), "math_Matrix() - invalid dimensions");
  89 }
  90
  91 void math_Matrix::Init(const Standard_Real InitialValue)
  92 {
  93   Array.Init(InitialValue);
  94 }
  95
  96 math_Matrix::math_Matrix (const math_Matrix& Other):
  97
  98 LowerRowIndex(Other.LowerRow()),
  99 UpperRowIndex(Other.UpperRow()),
 100 LowerColIndex(Other.LowerCol()),
 101 UpperColIndex(Other.UpperCol()),
 102 Array(Other.Array)
 103 {
 104 }
 105
 106
 107
 108 math_Matrix math_Matrix::Divided (const Standard_Real Right) const
 109 {
 110   Standard_DivideByZero_Raise_if (Abs(Right) <= RealEpsilon(),
 111                                   "math_Matrix::Divided() - zero divisor");
 112   math_Matrix temp = Multiplied(1./Right);
 113   return temp;
 114 }
 115
 116
 117 Standard_Real math_Matrix::Determinant() const
 118 {
 119   math_Gauss Sol(*this);
 120
 121   if(Sol.IsDone()) {
 122     return Sol.Determinant();
 123   }
 124   else {
 125     return 0.0;
 126   }
 127 }
 128
 129 void math_Matrix::Transpose()
 130 {
 131   math_NotSquare_Raise_if (RowNumber() != ColNumber(),
 132                            "math_Matrix::Transpose() - matrix is not square");
 133
 134   Standard_Integer Row = LowerRowIndex;
 135   Standard_Integer Col = LowerColIndex;
 136   SetLowerCol(LowerRowIndex);
 137   Standard_Real Temp;
 138   for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
 139     for(Standard_Integer J = I; J <= UpperColIndex; J++) {
 140       Temp = Array(I, J);
 141       Array(I, J) = Array(J, I);
 142       Array(J, I) = Temp;
 143     }
 144   }
 145   SetLowerRow(Col);
 146   SetLowerCol(Row);
 147 }
 148
 149 math_Matrix math_Matrix::Transposed() const
 150 {
 151   math_Matrix Result(LowerColIndex, UpperColIndex,
 152                      LowerRowIndex, UpperRowIndex);
 153
 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);
 157     }
 158   }
 159   return Result;
 160 }
 161
 162 void math_Matrix::Invert()
 163 {
 164   math_NotSquare_Raise_if (RowNumber() != ColNumber(),
 165                            "math_Matrix::Transpose() - matrix is not square");
 166
 167   math_Gauss Sol(*this);
 168   if(Sol.IsDone()) {
 169     Sol.Invert(*this);
 170   }
 171   else {
 172     throw math_SingularMatrix(); // SingularMatrix Exception;
 173   }
 174 }
 175
 176 math_Matrix math_Matrix::Inverse() const
 177 {
 178
 179   math_Matrix Result = *this;
 180   Result.Invert();
 181   return Result;
 182 }
 183
 184 void math_Matrix::Multiply (const Standard_Real Right)
 185 {
 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;
 189     }
 190   }
 191 }
 192
 193 math_Matrix math_Matrix::Multiplied (const Standard_Real Right) const
 194 {
 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;
 200     }
 201   }
 202   return Result;
 203 }
 204
 205 math_Matrix math_Matrix::TMultiplied (const Standard_Real Right) const
 206 {
 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;
 212     }
 213   }
 214   return Result;
 215 }
 216
 217
 218
 219 void math_Matrix::Divide (const Standard_Real Right)
 220 {
 221   Standard_DivideByZero_Raise_if (Abs(Right) <= RealEpsilon(),
 222                                   "math_Matrix::Divide() - zero divisor");
 223
 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;
 227     }
 228   }
 229 }
 230
 231 void math_Matrix::Add (const math_Matrix& Right)
 232 {
 233   Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
 234                                  || (ColNumber() != Right.ColNumber()),
 235                                     "math_Matrix::Add() - input matrix has different dimensions");
 236
 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);
 242       J2++;
 243     }
 244     I2++;
 245   }
 246 }
 247
 248 void math_Matrix::Subtract (const math_Matrix& Right)
 249 {
 250   Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
 251                                  || (ColNumber() != Right.ColNumber()),
 252                                     "math_Matrix::Subtract() - input matrix has different dimensions");
 253
 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);
 259       J2++;
 260     }
 261     I2++;
 262   }
 263 }
 264
 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)
 268 {
 269
 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");
 278
 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);
 284       JJ++;
 285     }
 286     II++;
 287   }
 288 }
 289
 290 void math_Matrix::SetRow (const Standard_Integer Row,
 291                           const math_Vector& V)
 292 {
 293   Standard_RangeError_Raise_if ((Row < LowerRowIndex)
 294                              || (Row > UpperRowIndex),
 295                                  "math_Matrix::SetRow() - invalid index");
 296
 297   Standard_DimensionError_Raise_if (ColNumber() != V.Length(),
 298                                     "math_Matrix::SetRow() - input vector has wrong dimensions");
 299
 300   Standard_Integer I = V.LowerIndex;
 301   for(Standard_Integer Index = LowerColIndex; Index <= UpperColIndex; Index++) {
 302     Array(Row, Index) = V.Array(I);
 303     I++;
 304   }
 305 }
 306
 307 void math_Matrix::SetCol (const Standard_Integer Col,
 308                           const math_Vector& V)
 309 {
 310   Standard_RangeError_Raise_if ((Col < LowerColIndex)
 311                              || (Col > UpperColIndex),
 312                                  "math_Matrix::SetCol() - invalid index");
 313
 314   Standard_DimensionError_Raise_if (RowNumber() != V.Length(),
 315                                     "math_Matrix::SetCol() - input vector has wrong dimensions");
 316
 317   Standard_Integer I = V.LowerIndex;
 318   for(Standard_Integer Index = LowerRowIndex; Index <= UpperRowIndex; Index++) {
 319     Array(Index, Col) = V.Array(I);
 320     I++;
 321   }
 322 }
 323
 324 void math_Matrix::SetDiag(const Standard_Real Value)
 325 {
 326   math_NotSquare_Raise_if (RowNumber() != ColNumber(),
 327                            "math_Matrix::SetDiag() - matrix is not square");
 328
 329   for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
 330     Array(I, I) = Value;
 331   }
 332 }
 333 math_Vector math_Matrix::Row (const Standard_Integer Row)  const
 334 {
 335
 336   math_Vector Result(LowerColIndex, UpperColIndex);
 337
 338   for(Standard_Integer Index = LowerColIndex; Index <= UpperColIndex; Index++) {
 339     Result.Array(Index) = Array(Row, Index);
 340   }
 341   return Result;
 342 }
 343
 344 math_Vector math_Matrix::Col (const Standard_Integer Col) const
 345 {
 346
 347   math_Vector Result(LowerRowIndex, UpperRowIndex);
 348
 349   for(Standard_Integer Index = LowerRowIndex; Index <= UpperRowIndex; Index++) {
 350     Result.Array(Index) = Array(Index, Col);
 351   }
 352   return Result;
 353 }
 354
 355 void math_Matrix::SwapRow(const Standard_Integer Row1,
 356                           const Standard_Integer Row2)
 357 {
 358   Standard_RangeError_Raise_if ((Row1 < LowerRowIndex)
 359                              || (Row1 > UpperRowIndex)
 360                              || (Row2 < LowerRowIndex)
 361                              || (Row2 > UpperRowIndex),
 362                                 "math_Matrix::SetCol() - invalid indices");
 363
 364   math_Vector V1 = Row(Row1);
 365   math_Vector V2 = Row(Row2);
 366   SetRow(Row1,V2);
 367   SetRow(Row2,V1);
 368 }
 369
 370 void math_Matrix::SwapCol(const Standard_Integer Col1,
 371                           const Standard_Integer Col2)
 372 {
 373   Standard_RangeError_Raise_if ((Col1 < LowerColIndex)
 374                              || (Col1 > UpperColIndex)
 375                              || (Col2 < LowerColIndex)
 376                              || (Col2 > UpperColIndex),
 377                                 "math_Matrix::SetCol() - invalid indices");
 378
 379   math_Vector V1 = Col(Col1);
 380   math_Vector V2 = Col(Col2);
 381   SetCol(Col1,V2);
 382   SetCol(Col2,V1);
 383 }
 384
 385
 386
 387 math_Matrix  math_Matrix::Multiplied (const math_Matrix& Right) const
 388 {
 389   Standard_DimensionError_Raise_if (ColNumber() != Right.RowNumber(),
 390                                     "math_Matrix::Multiplied() - matrices have incompatible dimensions");
 391
 392   math_Matrix Result(LowerRowIndex,       UpperRowIndex,
 393                      Right.LowerColIndex, Right.UpperColIndex);
 394
 395   Standard_Real Som;
 396   for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
 397     for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
 398       Som = 0.0;
 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);
 402         I2++;
 403       }
 404       Result.Array(I, J2) = Som;
 405     }
 406   }
 407   return Result;
 408 }
 409
 410 math_Matrix  math_Matrix::TMultiply (const math_Matrix& Right) const
 411 {
 412   Standard_DimensionError_Raise_if (RowNumber() != Right.RowNumber(),
 413                                     "math_Matrix::TMultiply() - matrices have incompatible dimensions");
 414
 415   math_Matrix Result(LowerColIndex,       UpperColIndex,
 416                      Right.LowerColIndex, Right.UpperColIndex);
 417
 418   Standard_Real Som;
 419   for(Standard_Integer I = LowerColIndex; I <= UpperColIndex; I++) {
 420     for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
 421       Som = 0.0;
 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);
 425         I2++;
 426       }
 427       Result.Array(I, J2) = Som;
 428     }
 429   }
 430   return Result;
 431 }
 432
 433 math_Matrix  math_Matrix::Added (const math_Matrix& Right) const
 434 {
 435   Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
 436                                  || (ColNumber() != Right.ColNumber()),
 437                                      "math_Matrix::Added() - input matrix has different dimensions");
 438
 439   math_Matrix Result(LowerRowIndex, UpperRowIndex,
 440                      LowerColIndex, UpperColIndex);
 441
 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);
 447       J2++;
 448     }
 449     I2++;
 450   }
 451   return Result;
 452 }
 453
 454 math_Matrix  math_Matrix::Opposite ()
 455 {
 456
 457   math_Matrix Result(LowerRowIndex, UpperRowIndex,
 458                      LowerColIndex, UpperColIndex);
 459
 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);
 463     }
 464   }
 465   return Result;
 466 }
 467
 468 math_Matrix  math_Matrix::Subtracted (const math_Matrix& Right) const
 469 {
 470   Standard_DimensionError_Raise_if ((RowNumber() != Right.RowNumber())
 471                                  || (ColNumber() != Right.ColNumber()),
 472                                      "math_Matrix::Subtracted() - input matrix has different dimensions");
 473
 474   math_Matrix Result(LowerRowIndex, UpperRowIndex,
 475                      LowerColIndex, UpperColIndex);
 476
 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);
 482       J2++;
 483     }
 484     I2++;
 485   }
 486   return Result;
 487 }
 488
 489 void  math_Matrix::Multiply(const math_Vector&  Left,
 490                             const math_Vector&  Right)
 491 {
 492   Standard_DimensionError_Raise_if ((RowNumber() != Left.Length())
 493                                  || (ColNumber() != Right.Length()),
 494                                     "math_Matrix::Multiply() - input vectors have incompatible dimensions");
 495
 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);
 499     }
 500   }
 501 }
 502
 503 void  math_Matrix::Multiply(const math_Matrix&  Left,
 504                             const math_Matrix&  Right)
 505 {
 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");
 510
 511   Standard_Real Som;
 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++) {
 516       Som = 0.0;
 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);
 521         J1++;
 522         I2++;
 523       }
 524       Array(I, J) = Som;
 525       J2++;
 526     }
 527     I1++;
 528   }
 529 }
 530
 531 void math_Matrix::TMultiply(const math_Matrix& TLeft,
 532                             const math_Matrix&  Right)
 533 {
 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");
 538
 539   Standard_Real Som;
 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++) {
 544       Som = 0.0;
 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);
 549         J1++;
 550         I2++;
 551       }
 552       Array(I, J) = Som;
 553       J2++;
 554     }
 555     I1++;
 556   }
 557 }
 558
 559 void  math_Matrix::Add (const math_Matrix&  Left,
 560                         const math_Matrix&  Right)
 561 {
 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");
 567
 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);
 575       J1++;
 576       J2++;
 577     }
 578     I1++;
 579     I2++;
 580   }
 581 }
 582
 583 void  math_Matrix::Subtract(const math_Matrix&  Left,
 584                             const math_Matrix&  Right)
 585 {
 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");
 591
 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);
 599       J1++;
 600       J2++;
 601     }
 602     I1++;
 603     I2++;
 604   }
 605 }
 606
 607
 608 void math_Matrix::Multiply(const math_Matrix& Right)
 609 {
 610   Standard_DimensionError_Raise_if (ColNumber() != Right.RowNumber(),
 611                                     "math_Matrix::Multiply() - input matrix has incompatible dimensions");
 612
 613   Standard_Real Som;
 614   for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
 615     for(Standard_Integer J2 = Right.LowerColIndex; J2 <= Right.UpperColIndex; J2++) {
 616       Som = 0.0;
 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);
 620         I2++;
 621       }
 622       Array(I, J2) = Som;
 623     }
 624   }
 625 }
 626
 627
 628
 629 math_Vector math_Matrix::Multiplied(const math_Vector& Right)const
 630 {
 631   Standard_DimensionError_Raise_if (ColNumber() != Right.Length(),
 632                                     "math_Matrix::Multiplied() - input vector has incompatible dimensions");
 633
 634   math_Vector Result(LowerRowIndex, UpperRowIndex);
 635
 636   for(Standard_Integer I = LowerRowIndex; I <= UpperRowIndex; I++) {
 637     Result.Array(I) = 0.0;
 638     Standard_Integer II = Right.LowerIndex;
 639     for(Standard_Integer J = LowerColIndex; J <= UpperColIndex; J++) {
 640       Result.Array(I) = Result.Array(I) + Array(I, J) * Right.Array(II);
 641       II++;
 642     }
 643   }
 644   return Result;
 645 }
 646
 647 math_Matrix& math_Matrix::Initialized(const math_Matrix& Other)
 648 {
 649   Standard_DimensionError_Raise_if ((RowNumber() != Other.RowNumber())
 650                                  || (ColNumber() != Other.ColNumber()),
 651                                     "math_Matrix::Initialized() - input matrix has different dimensions");
 652
 653   (Other.Array).Copy(Array);
 654   return *this;
 655 }
 656
 657
 658
 659
 660 void math_Matrix::Dump(Standard_OStream& o)const
 661
 662 {
 663   o << "math_Matrix of RowNumber = " << RowNumber();
 664   o << " and ColNumber = " << ColNumber() << "\n";
 665
 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";
 670     }
 671   }
 672 }
 673