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