src/math/math_Crout.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 // lpa le 20/08/91
  16
  17 //#ifndef OCCT_DEBUG
  18 #define No_Standard_RangeError
  19 #define No_Standard_OutOfRange
  20 #define No_Standard_DimensionError
  21
  22 //#endif
  23
  24 #include <math_Crout.hxx>
  25 #include <math_Matrix.hxx>
  26 #include <math_NotSquare.hxx>
  27 #include <math_Vector.hxx>
  28 #include <Standard_DimensionError.hxx>
  29 #include <StdFail_NotDone.hxx>
  30
  31 math_Crout::math_Crout(const math_Matrix& A, const Standard_Real MinPivot):
  32                        InvA(1, A.RowNumber(), 1, A.ColNumber())
  33 {
  34   Standard_Integer i,j,k;
  35   Standard_Integer Nctl = A.RowNumber();
  36   Standard_Integer lowr = A.LowerRow(), lowc = A.LowerCol();
  37   Standard_Real scale;
  38
  39   math_Matrix L(1, Nctl, 1, Nctl);
  40   math_Vector Diag(1, Nctl);
  41
  42
  43
  44   math_NotSquare_Raise_if(Nctl != A.ColNumber(), " ");
  45
  46   Det = 1;
  47   for (i =1; i <= Nctl; i++) {
  48     for (j = 1; j <= i -1; j++) {
  49       scale = 0.0;
  50       for (k = 1; k <= j-1; k++) {
  51         scale += L(i,k)*L(j,k)*Diag(k);
  52       }
  53       L(i,j) = (A(i+lowr-1,j+lowc-1)-scale)/Diag(j);
  54     }
  55     scale = 0.0;
  56     for (k = 1; k <= i-1; k++) {
  57       scale += L(i,k)*L(i,k)*Diag(k);
  58     }
  59     Diag(i) = A(i+lowr-1,i+lowc-1)-scale;
  60     Det *= Diag(i);
  61     if (Abs(Diag(i)) <= MinPivot) {
  62       Done = Standard_False;
  63       return;
  64     }
  65     L(i,i) = 1.0;
  66   }
  67
  68 // Calcul de l inverse de L:
  69 //==========================
  70
  71   L(1,1) = 1./L(1,1);
  72   for (i = 2; i <= Nctl; i++) {
  73     for (k = 1; k <= i-1; k++) {
  74       scale = 0.0;
  75       for (j = k; j <= i-1; j++) {
  76         scale += L(i,j)*L(j,k);
  77       }
  78       L(i,k) = -scale/L(i,i);
  79     }
  80     L(i,i) = 1./L(i,i);
  81   }
  82
  83 // Calcul de l inverse de Mat:
  84 //============================
  85
  86   for (j = 1; j <= Nctl; j++) {
  87     scale = L(j,j)*L(j,j)/Diag(j);
  88     for (k = j+1; k <= Nctl; k++) {
  89       scale += L(k,j) *L(k,j)/Diag(k);
  90     }
  91     InvA(j,j) = scale;
  92     for (i = j+1; i <= Nctl; i++) {
  93       scale = L(i,j) *L(i,i)/Diag(i);
  94       for (k = i+1; k <= Nctl; k++) {
  95         scale += L(k,j)*L(k,i)/Diag(k);
  96       }
  97       InvA(i,j) = scale;
  98     }
  99   }
 100   Done = Standard_True;
 101 }
 102
 103
 104
 105 void math_Crout::Solve(const math_Vector& B, math_Vector& X) const
 106 {
 107   StdFail_NotDone_Raise_if(!Done, " ");
 108   Standard_DimensionError_Raise_if((B.Length() != InvA.RowNumber()) ||
 109                                    (X.Length() != B.Length()), " ");
 110
 111   Standard_Integer n = InvA.RowNumber();
 112   Standard_Integer lowb = B.Lower(), lowx = X.Lower();
 113   Standard_Integer i, j;
 114
 115   for (i = 1; i <= n; i++) {
 116     X(i+lowx-1) = InvA(i, 1)*B(1+lowb-1);
 117     for ( j = 2; j <= i; j++) {
 118       X(i+lowx-1) += InvA(i, j)*B(j+lowb-1);
 119     }
 120     for (j = i+1; j <= n; j++) {
 121       X(i+lowx-1) += InvA(j,i)*B(j+lowb-1);
 122     }
 123   }
 124 }
 125
 126 void math_Crout::Dump(Standard_OStream& o) const
 127 {
 128   o << "math_Crout ";
 129   if(Done) {
 130     o << " Status = Done \n";
 131   }
 132   else {
 133     o << " Status = not Done \n";
 134   }
 135 }
 136
 137
 138
 139
 140
 141
 142