src/math/math_Crout.cxx

   1 // Copyright (c) 1997-1999 Matra Datavision
   2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   3 //
   4 // The content of this file is subject to the Open CASCADE Technology Public
   5 // License Version 6.5 (the "License"). You may not use the content of this file
   6 // except in compliance with the License. Please obtain a copy of the License
   7 // at http://www.opencascade.org and read it completely before using this file.
   8 //
   9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  11 //
  12 // The Original Code and all software distributed under the License is
  13 // distributed on an "AS IS" basis, without warranty of any kind, and the
  14 // Initial Developer hereby disclaims all such warranties, including without
  15 // limitation, any warranties of merchantability, fitness for a particular
  16 // purpose or non-infringement. Please see the License for the specific terms
  17 // and conditions governing the rights and limitations under the License.
  18
  19 // lpa le 20/08/91
  20
  21 //#ifndef DEB
  22 #define No_Standard_RangeError
  23 #define No_Standard_OutOfRange
  24 #define No_Standard_DimensionError
  25 //#endif
  26
  27 #include <math_Crout.ixx>
  28 #include <math_NotSquare.hxx>
  29 #include <StdFail_NotDone.hxx>
  30 #include <math_Vector.hxx>
  31
  32 math_Crout::math_Crout(const math_Matrix& A, const Standard_Real MinPivot):
  33                        InvA(1, A.RowNumber(), 1, A.ColNumber())
  34 {
  35   Standard_Integer i,j,k;
  36   Standard_Integer Nctl = A.RowNumber();
  37   Standard_Integer lowr = A.LowerRow(), lowc = A.LowerCol();
  38   Standard_Real scale;
  39
  40   math_Matrix L(1, Nctl, 1, Nctl);
  41   math_Vector Diag(1, Nctl);
  42
  43
  44
  45   math_NotSquare_Raise_if(Nctl != A.ColNumber(), " ");
  46
  47   Det = 1;
  48   for (i =1; i <= Nctl; i++) {
  49     for (j = 1; j <= i -1; j++) {
  50       scale = 0.0;
  51       for (k = 1; k <= j-1; k++) {
  52         scale += L(i,k)*L(j,k)*Diag(k);
  53       }
  54       L(i,j) = (A(i+lowr-1,j+lowc-1)-scale)/Diag(j);
  55     }
  56     scale = 0.0;
  57     for (k = 1; k <= i-1; k++) {
  58       scale += L(i,k)*L(i,k)*Diag(k);
  59     }
  60     Diag(i) = A(i+lowr-1,i+lowc-1)-scale;
  61     Det *= Diag(i);
  62     if (Abs(Diag(i)) <= MinPivot) {
  63       Done = Standard_False;
  64       return;
  65     }
  66     L(i,i) = 1.0;
  67   }
  68
  69 // Calcul de l inverse de L:
  70 //==========================
  71
  72   L(1,1) = 1./L(1,1);
  73   for (i = 2; i <= Nctl; i++) {
  74     for (k = 1; k <= i-1; k++) {
  75       scale = 0.0;
  76       for (j = k; j <= i-1; j++) {
  77         scale += L(i,j)*L(j,k);
  78       }
  79       L(i,k) = -scale/L(i,i);
  80     }
  81     L(i,i) = 1./L(i,i);
  82   }
  83
  84 // Calcul de l inverse de Mat:
  85 //============================
  86
  87   for (j = 1; j <= Nctl; j++) {
  88     scale = L(j,j)*L(j,j)/Diag(j);
  89     for (k = j+1; k <= Nctl; k++) {
  90       scale += L(k,j) *L(k,j)/Diag(k);
  91     }
  92     InvA(j,j) = scale;
  93     for (i = j+1; i <= Nctl; i++) {
  94       scale = L(i,j) *L(i,i)/Diag(i);
  95       for (k = i+1; k <= Nctl; k++) {
  96         scale += L(k,j)*L(k,i)/Diag(k);
  97       }
  98       InvA(i,j) = scale;
  99     }
 100   }
 101   Done = Standard_True;
 102 }
 103
 104
 105
 106 void math_Crout::Solve(const math_Vector& B, math_Vector& X) const
 107 {
 108   StdFail_NotDone_Raise_if(!Done, " ");
 109   Standard_DimensionError_Raise_if((B.Length() != InvA.RowNumber()) ||
 110                                    (X.Length() != B.Length()), " ");
 111
 112   Standard_Integer n = InvA.RowNumber();
 113   Standard_Integer lowb = B.Lower(), lowx = X.Lower();
 114   Standard_Integer i, j;
 115
 116   for (i = 1; i <= n; i++) {
 117     X(i+lowx-1) = InvA(i, 1)*B(1+lowb-1);
 118     for ( j = 2; j <= i; j++) {
 119       X(i+lowx-1) += InvA(i, j)*B(j+lowb-1);
 120     }
 121     for (j = i+1; j <= n; j++) {
 122       X(i+lowx-1) += InvA(j,i)*B(j+lowb-1);
 123     }
 124   }
 125 }
 126
 127 void math_Crout::Dump(Standard_OStream& o) const
 128 {
 129   o << "math_Crout ";
 130   if(Done) {
 131     o << " Status = Done \n";
 132   }
 133   else {
 134     o << " Status = not Done \n";
 135   }
 136 }
 137
 138
 139
 140
 141
 142
 143