src/math/math_Gauss.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 //#ifndef OCCT_DEBUG
  16 #define No_Standard_RangeError
  17 #define No_Standard_OutOfRange
  18 #define No_Standard_DimensionError
  19
  20 //#endif
  21
  22 #include <math_Gauss.hxx>
  23 #include <math_Matrix.hxx>
  24 #include <math_NotSquare.hxx>
  25 #include <math_Recipes.hxx>
  26 #include <Standard_DimensionError.hxx>
  27 #include <Standard_NotImplemented.hxx>
  28 #include <StdFail_NotDone.hxx>
  29
  30 math_Gauss::math_Gauss(const math_Matrix& A,
  31                            const Standard_Real MinPivot)
  32                            : LU   (1, A.RowNumber(), 1, A.ColNumber()),
  33                              Index(1, A.RowNumber()) {
  34
  35       math_NotSquare_Raise_if(A.RowNumber() != A.ColNumber(), " ");
  36       LU = A;
  37       Standard_Integer Error = LU_Decompose(LU,
  38                                   Index,
  39                                   D,
  40                                   MinPivot);
  41       if(!Error) {
  42         Done = Standard_True;
  43       }
  44       else {
  45         Done = Standard_False;
  46       }
  47     }
  48
  49     void  math_Gauss::Solve(const math_Vector& B, math_Vector& X) const{
  50
  51        StdFail_NotDone_Raise_if(!Done, " ");
  52
  53        X = B;
  54        LU_Solve(LU,
  55                 Index,
  56                 X);
  57     }
  58
  59     void math_Gauss::Solve (math_Vector& X) const{
  60
  61        StdFail_NotDone_Raise_if(!Done, " ");
  62
  63        if(X.Length() != LU.RowNumber()) {
  64          Standard_DimensionError::Raise();
  65        }
  66        LU_Solve(LU,
  67                 Index,
  68                 X);
  69     }
  70
  71     Standard_Real math_Gauss::Determinant() const{
  72
  73        StdFail_NotDone_Raise_if(!Done, " ");
  74
  75        Standard_Real Result = D;
  76        for(Standard_Integer J = 1; J <= LU.UpperRow(); J++) {
  77          Result *= LU(J,J);
  78        }
  79        return Result;
  80     }
  81
  82     void math_Gauss::Invert(math_Matrix& Inv) const{
  83
  84        StdFail_NotDone_Raise_if(!Done, " ");
  85
  86        Standard_DimensionError_Raise_if((Inv.RowNumber() != LU.RowNumber()) ||
  87                                         (Inv.ColNumber() != LU.ColNumber()),
  88                                         " ");
  89
  90        Standard_Integer LowerRow = Inv.LowerRow();
  91        Standard_Integer LowerCol = Inv.LowerCol();
  92        math_Vector Column(1, LU.UpperRow());
  93
  94        Standard_Integer I, J;
  95        for(J = 1; J <= LU.UpperRow(); J++) {
  96          for(I = 1; I <= LU.UpperRow(); I++) {
  97            Column(I) = 0.0;
  98          }
  99          Column(J) = 1.0;
 100          LU_Solve(LU, Index, Column);
 101          for(I = 1; I <= LU.RowNumber(); I++) {
 102            Inv(I+LowerRow-1,J+LowerCol-1) = Column(I);
 103          }
 104        }
 105
 106     }
 107
 108
 109 void math_Gauss::Dump(Standard_OStream& o) const {
 110        o << "math_Gauss ";
 111        if(Done) {
 112          o<< " Status = Done \n";
 113          o << " Determinant of A = " << D << endl;
 114        }
 115        else {
 116          o << " Status = not Done \n";
 117        }
 118     }
 119
 120
 121
 122