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