7fd59977 |
1 | // File: FEmTool_ProfileMatrix.cxx |
2 | // Created: Fri Oct 31 10:02:27 1997 |
3 | // Author: Roman BORISOV |
4 | // <rbv@velox.nnov.matra-dtv.fr> |
5 | |
6 | #define No_Standard_RangeError |
7 | #define No_Standard_OutOfRange |
8 | #define No_Standard_DimensionError |
9 | |
10 | #include <FEmTool_ProfileMatrix.ixx> |
11 | #include <gp.hxx> |
12 | |
13 | |
14 | //======================================================================= |
15 | //function : :FEmTool_ProfileMatrix |
16 | //purpose : |
17 | //======================================================================= |
18 | FEmTool_ProfileMatrix::FEmTool_ProfileMatrix(const TColStd_Array1OfInteger& FirstIndexes) |
19 | : profile(1, 2, 1, FirstIndexes.Length()) |
20 | { |
21 | Standard_Integer i, j, k, l; |
22 | profile(1, 1) = 0; |
23 | profile(2, 1) = 1; |
24 | for(i = 2; i <= FirstIndexes.Length(); i++) { |
25 | profile(1, i) = i - FirstIndexes(i); |
26 | profile(2, i) = profile(2, i-1) + profile(1, i) + 1; |
27 | } |
28 | NextCoeff = new TColStd_HArray1OfInteger(1, profile(2, FirstIndexes.Length())); |
29 | |
30 | for(i = 1, k = 1; i <= FirstIndexes.Length(); i++) |
31 | for(j = FirstIndexes(i); j <= i; j++) { |
32 | for(l = i+1; l <= FirstIndexes.Length() && j < FirstIndexes(l); l++); |
33 | |
34 | if(l > FirstIndexes.Length()) NextCoeff->SetValue(k, 0); |
35 | else NextCoeff->SetValue(k, l); |
36 | k++; |
37 | } |
38 | |
39 | |
40 | ProfileMatrix = new TColStd_HArray1OfReal(1, profile(2, FirstIndexes.Length())); |
41 | SMatrix = new TColStd_HArray1OfReal(1, profile(2, FirstIndexes.Length())); |
42 | IsDecomp = Standard_False; |
43 | } |
44 | |
45 | //======================================================================= |
46 | //function : Init |
47 | //purpose : |
48 | //======================================================================= |
49 | |
50 | void FEmTool_ProfileMatrix::Init(const Standard_Real Value) |
51 | { |
52 | ProfileMatrix->Init(Value); |
53 | IsDecomp = Standard_False; |
54 | } |
55 | |
56 | //======================================================================= |
57 | //function : ChangeValue |
58 | //purpose : |
59 | //======================================================================= |
60 | |
61 | Standard_Real& FEmTool_ProfileMatrix::ChangeValue(const Standard_Integer I, |
62 | const Standard_Integer J) |
63 | { |
64 | Standard_Integer Ind; |
65 | Ind = I-J; |
66 | if (Ind < 0) { |
67 | Ind = -Ind; |
68 | Standard_OutOfRange_Raise_if(Ind>profile(1, J), |
69 | "FEmTool_ProfileMatrix::ChangeValue"); |
70 | Ind = profile(2, J) - Ind; |
71 | } |
72 | else { |
73 | Standard_OutOfRange_Raise_if(Ind>profile(1, I), |
74 | "FEmTool_ProfileMatrix::ChangeValue"); |
75 | Ind = profile(2, I)-Ind; |
76 | } |
77 | return ProfileMatrix->ChangeValue(Ind); |
78 | } |
79 | |
80 | //======================================================================= |
81 | //function : Decompose |
82 | //purpose : Choleski's decomposition |
83 | //======================================================================= |
84 | Standard_Boolean FEmTool_ProfileMatrix::Decompose() |
85 | { |
86 | Standard_Integer i, j, k, ik, jk, DiagAddr, CurrAddr, Kmin, Kj; |
87 | Standard_Real Sum, a, Eps = 1.e-32; |
88 | |
89 | SMatrix->Init(0.); |
90 | Standard_Real * SMA = &SMatrix->ChangeValue(1); |
91 | SMA--; |
92 | const Standard_Real * PM = &ProfileMatrix->Value(1); |
93 | PM--; |
94 | |
95 | for(j = 1; j <= RowNumber(); j++) { |
96 | DiagAddr = profile(2, j); |
97 | Kj = j - profile(1, j); |
98 | Sum = 0; |
99 | for(k = DiagAddr - profile(1, j); k < DiagAddr; k++) |
100 | Sum += SMA[k] * SMA[k]; |
101 | |
102 | a = PM[DiagAddr] - Sum; |
103 | if(a < Eps) { |
104 | return IsDecomp = Standard_False;// Matrix is not positive defined |
105 | } |
106 | a = Sqrt(a); |
107 | SMA[DiagAddr] = a; |
108 | |
109 | CurrAddr = DiagAddr; |
110 | while ((i = NextCoeff->Value(CurrAddr)) > 0) { |
111 | CurrAddr = profile(2, i) - (i - j); |
112 | |
113 | // Computation of Sum of S .S for k = 1,..,j-1 |
114 | // ik jk |
115 | Sum = 0; |
116 | Kmin = Max((i - profile(1, i)), Kj); |
117 | ik = profile(2, i) - i + Kmin; |
118 | jk = DiagAddr - j + Kmin; |
119 | for(k = Kmin; k <j; k++,ik++,jk++) { |
120 | Sum += SMA[ik]*SMA[jk]; |
121 | } |
122 | SMA[CurrAddr] = (PM[CurrAddr] - Sum)/a; |
123 | } |
124 | } |
125 | return IsDecomp = Standard_True; |
126 | } |
127 | |
128 | //======================================================================= |
129 | //function : Solve |
130 | //purpose : Resolution of the system S*t(S)X = B |
131 | //======================================================================= |
132 | void FEmTool_ProfileMatrix::Solve(const math_Vector& B,math_Vector& X) const |
133 | { |
134 | if (!IsDecomp) StdFail_NotDone::Raise("Decomposition must be done"); |
135 | |
136 | Standard_Integer i, j, jj,DiagAddr, CurrAddr; |
137 | Standard_Real Sum; |
138 | |
139 | Standard_Real * x = &X(X.Lower()); |
140 | x--; |
141 | Standard_Real * b = &B(B.Lower()); |
142 | b--; |
143 | const Standard_Real * SMA = &SMatrix->Value(1); |
144 | SMA --; |
145 | const Standard_Integer * NC = &NextCoeff->Value(1); |
146 | NC--; |
147 | |
148 | // Resolution of Sw = B; |
149 | for(i = 1; i <= RowNumber(); i++) { |
150 | DiagAddr = profile(2, i); |
151 | Sum = 0; |
152 | for(j = i - profile(1, i), jj = DiagAddr - (i - j); |
153 | j < i; j++, jj++) |
154 | Sum += SMA[jj]* x[j]; |
155 | x[i] = (b[i] - Sum)/SMA[DiagAddr]; |
156 | } |
157 | |
158 | // Resolution of t(S)X = w; |
159 | for(i = ColNumber(); i >= 1; i--) { |
160 | DiagAddr = profile(2, i); |
161 | j = NC[DiagAddr]; |
162 | Sum = 0; |
163 | while(j > 0) { |
164 | CurrAddr = profile(2, j) - (j-i); |
165 | Sum += SMA[CurrAddr]*x[j]; |
166 | j = NC[CurrAddr]; |
167 | } |
168 | x[i] = (x[i] - Sum)/SMA[DiagAddr]; |
169 | } |
170 | } |
171 | |
172 | Standard_Boolean FEmTool_ProfileMatrix::Prepare() |
173 | { |
174 | Standard_NotImplemented::Raise("FEmTool_ProfileMatrix::Prepare"); |
175 | return Standard_False; |
176 | } |
177 | |
178 | // void FEmTool_ProfileMatrix::Solve(const math_Vector& B,const math_Vector& Init,math_Vector& X,math_Vector& Residual,const Standard_Real Tolerance,const Standard_Integer NbIterations) const |
179 | void FEmTool_ProfileMatrix::Solve(const math_Vector& ,const math_Vector& ,math_Vector& ,math_Vector& ,const Standard_Real ,const Standard_Integer ) const |
180 | { |
181 | Standard_NotImplemented::Raise("FEmTool_ProfileMatrix::Solve"); |
182 | } |
183 | |
184 | |
185 | //======================================================================= |
186 | //function : Multiplied |
187 | //purpose : MX = H*X |
188 | //======================================================================= |
189 | void FEmTool_ProfileMatrix::Multiplied(const math_Vector& X,math_Vector& MX) const |
190 | { |
191 | Standard_Integer i, j, jj, DiagAddr, CurrAddr; |
192 | Standard_Real * m = &MX(MX.Lower()); |
193 | m--; |
194 | Standard_Real * x = &X(X.Lower()); |
195 | x--; |
196 | const Standard_Real * PM = &ProfileMatrix->Value(1); |
197 | PM--; |
198 | const Standard_Integer * NC = &NextCoeff->Value(1); |
199 | NC--; |
200 | |
201 | for(i = 1; i <= RowNumber(); i++) { |
202 | DiagAddr = profile(2, i); |
203 | m[i] = 0; |
204 | for(j = i - profile(1, i), jj = DiagAddr - (i - j); |
205 | j <= i; j++, jj++) |
206 | m[i] += PM[jj]*x[j]; |
207 | |
208 | CurrAddr = DiagAddr; |
209 | for(j = NC[CurrAddr]; j > 0; j = NC[CurrAddr]) { |
210 | CurrAddr = profile(2, j) - (j-i); |
211 | m[i] += PM[CurrAddr]*x[j]; |
212 | } |
213 | } |
214 | } |
215 | |
216 | Standard_Integer FEmTool_ProfileMatrix::RowNumber() const |
217 | { |
218 | return profile.RowLength(); |
219 | } |
220 | |
221 | Standard_Integer FEmTool_ProfileMatrix::ColNumber() const |
222 | { |
223 | return profile.RowLength(); |
224 | } |
225 | |
226 | Standard_Boolean FEmTool_ProfileMatrix::IsInProfile(const Standard_Integer i, |
227 | const Standard_Integer j) const |
228 | { |
229 | if (j <= i) { |
230 | if ((i - j) <= profile(1, i)) return Standard_True; |
231 | else return Standard_False; |
232 | } |
233 | else if ((j - i) <= profile(1, j)) return Standard_True; |
234 | |
235 | return Standard_False; |
236 | } |
237 | |
238 | void FEmTool_ProfileMatrix::OutM() const |
239 | { |
240 | Standard_Integer i, j; |
241 | cout<<"Matrix A"<<endl; |
242 | for(i = 1; i <= RowNumber(); i++) { |
243 | for(j = 1; j < i - profile(1, i); j++) |
244 | cout<<"0 "; |
245 | |
246 | for(j = profile(2, i) - profile(1, i); j <= profile(2, i); j++) |
247 | cout<<ProfileMatrix->Value(j)<<" "; |
248 | cout<<endl; |
249 | } |
250 | |
251 | cout<<"NextCoeff"<<endl; |
252 | for(i = 1; i <= profile(2, RowNumber()); i++) |
253 | cout<<NextCoeff->Value(i)<<" "; |
254 | cout<<endl; |
255 | } |
256 | |
257 | void FEmTool_ProfileMatrix::OutS() const |
258 | { |
259 | Standard_Integer i, j; |
260 | cout<<"Matrix S"<<endl; |
261 | for(i = 1; i <= RowNumber(); i++) { |
262 | for(j = 1; j < i - profile(1, i); j++) |
263 | cout<<"0 "; |
264 | |
265 | for(j = profile(2, i) - profile(1, i); j <= profile(2, i); j++) |
266 | cout<<SMatrix->Value(j)<<" "; |
267 | cout<<endl; |
268 | } |
269 | } |