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