b311480e |
1 | // Created on: 1997-10-31 |
2 | // Created by: Roman BORISOV |
3 | // Copyright (c) 1997-1999 Matra Datavision |
973c2be1 |
4 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
b311480e |
5 | // |
973c2be1 |
6 | // This file is part of Open CASCADE Technology software library. |
b311480e |
7 | // |
d5f74e42 |
8 | // This library is free software; you can redistribute it and/or modify it under |
9 | // the terms of the GNU Lesser General Public License version 2.1 as published |
973c2be1 |
10 | // by the Free Software Foundation, with special exception defined in the file |
11 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
12 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
13 | // |
973c2be1 |
14 | // Alternatively, this file may be used under the terms of Open CASCADE |
15 | // commercial license or contractual agreement. |
7fd59977 |
16 | |
17 | #define No_Standard_RangeError |
18 | #define No_Standard_OutOfRange |
19 | #define No_Standard_DimensionError |
20 | |
7fd59977 |
21 | |
42cf5bc1 |
22 | #include <FEmTool_ProfileMatrix.hxx> |
23 | #include <gp.hxx> |
24 | #include <Standard_NotImplemented.hxx> |
25 | #include <Standard_OutOfRange.hxx> |
26 | #include <Standard_Type.hxx> |
27 | #include <StdFail_NotDone.hxx> |
7fd59977 |
28 | |
92efcf78 |
29 | IMPLEMENT_STANDARD_RTTIEXT(FEmTool_ProfileMatrix,FEmTool_SparseMatrix) |
30 | |
7fd59977 |
31 | //======================================================================= |
32 | //function : :FEmTool_ProfileMatrix |
33 | //purpose : |
34 | //======================================================================= |
b311480e |
35 | FEmTool_ProfileMatrix::FEmTool_ProfileMatrix(const TColStd_Array1OfInteger& FirstIndexes) |
7fd59977 |
36 | : profile(1, 2, 1, FirstIndexes.Length()) |
37 | { |
38 | Standard_Integer i, j, k, l; |
39 | profile(1, 1) = 0; |
40 | profile(2, 1) = 1; |
41 | for(i = 2; i <= FirstIndexes.Length(); i++) { |
42 | profile(1, i) = i - FirstIndexes(i); |
43 | profile(2, i) = profile(2, i-1) + profile(1, i) + 1; |
44 | } |
45 | NextCoeff = new TColStd_HArray1OfInteger(1, profile(2, FirstIndexes.Length())); |
46 | |
47 | for(i = 1, k = 1; i <= FirstIndexes.Length(); i++) |
48 | for(j = FirstIndexes(i); j <= i; j++) { |
49 | for(l = i+1; l <= FirstIndexes.Length() && j < FirstIndexes(l); l++); |
50 | |
51 | if(l > FirstIndexes.Length()) NextCoeff->SetValue(k, 0); |
52 | else NextCoeff->SetValue(k, l); |
53 | k++; |
54 | } |
55 | |
56 | |
57 | ProfileMatrix = new TColStd_HArray1OfReal(1, profile(2, FirstIndexes.Length())); |
58 | SMatrix = new TColStd_HArray1OfReal(1, profile(2, FirstIndexes.Length())); |
59 | IsDecomp = Standard_False; |
60 | } |
61 | |
62 | //======================================================================= |
63 | //function : Init |
64 | //purpose : |
65 | //======================================================================= |
66 | |
67 | void FEmTool_ProfileMatrix::Init(const Standard_Real Value) |
68 | { |
69 | ProfileMatrix->Init(Value); |
70 | IsDecomp = Standard_False; |
71 | } |
72 | |
73 | //======================================================================= |
74 | //function : ChangeValue |
75 | //purpose : |
76 | //======================================================================= |
77 | |
78 | Standard_Real& FEmTool_ProfileMatrix::ChangeValue(const Standard_Integer I, |
79 | const Standard_Integer J) |
80 | { |
81 | Standard_Integer Ind; |
82 | Ind = I-J; |
83 | if (Ind < 0) { |
84 | Ind = -Ind; |
85 | Standard_OutOfRange_Raise_if(Ind>profile(1, J), |
86 | "FEmTool_ProfileMatrix::ChangeValue"); |
87 | Ind = profile(2, J) - Ind; |
88 | } |
89 | else { |
90 | Standard_OutOfRange_Raise_if(Ind>profile(1, I), |
91 | "FEmTool_ProfileMatrix::ChangeValue"); |
92 | Ind = profile(2, I)-Ind; |
93 | } |
94 | return ProfileMatrix->ChangeValue(Ind); |
95 | } |
96 | |
97 | //======================================================================= |
98 | //function : Decompose |
99 | //purpose : Choleski's decomposition |
100 | //======================================================================= |
101 | Standard_Boolean FEmTool_ProfileMatrix::Decompose() |
102 | { |
103 | Standard_Integer i, j, k, ik, jk, DiagAddr, CurrAddr, Kmin, Kj; |
104 | Standard_Real Sum, a, Eps = 1.e-32; |
105 | |
106 | SMatrix->Init(0.); |
107 | Standard_Real * SMA = &SMatrix->ChangeValue(1); |
108 | SMA--; |
109 | const Standard_Real * PM = &ProfileMatrix->Value(1); |
110 | PM--; |
111 | |
112 | for(j = 1; j <= RowNumber(); j++) { |
113 | DiagAddr = profile(2, j); |
114 | Kj = j - profile(1, j); |
115 | Sum = 0; |
116 | for(k = DiagAddr - profile(1, j); k < DiagAddr; k++) |
117 | Sum += SMA[k] * SMA[k]; |
118 | |
119 | a = PM[DiagAddr] - Sum; |
120 | if(a < Eps) { |
121 | return IsDecomp = Standard_False;// Matrix is not positive defined |
122 | } |
123 | a = Sqrt(a); |
124 | SMA[DiagAddr] = a; |
125 | |
126 | CurrAddr = DiagAddr; |
127 | while ((i = NextCoeff->Value(CurrAddr)) > 0) { |
128 | CurrAddr = profile(2, i) - (i - j); |
129 | |
130 | // Computation of Sum of S .S for k = 1,..,j-1 |
131 | // ik jk |
132 | Sum = 0; |
133 | Kmin = Max((i - profile(1, i)), Kj); |
134 | ik = profile(2, i) - i + Kmin; |
135 | jk = DiagAddr - j + Kmin; |
136 | for(k = Kmin; k <j; k++,ik++,jk++) { |
137 | Sum += SMA[ik]*SMA[jk]; |
138 | } |
139 | SMA[CurrAddr] = (PM[CurrAddr] - Sum)/a; |
140 | } |
141 | } |
142 | return IsDecomp = Standard_True; |
143 | } |
144 | |
145 | //======================================================================= |
146 | //function : Solve |
147 | //purpose : Resolution of the system S*t(S)X = B |
148 | //======================================================================= |
149 | void FEmTool_ProfileMatrix::Solve(const math_Vector& B,math_Vector& X) const |
150 | { |
151 | if (!IsDecomp) StdFail_NotDone::Raise("Decomposition must be done"); |
152 | |
153 | Standard_Integer i, j, jj,DiagAddr, CurrAddr; |
154 | Standard_Real Sum; |
155 | |
156 | Standard_Real * x = &X(X.Lower()); |
157 | x--; |
158 | Standard_Real * b = &B(B.Lower()); |
159 | b--; |
160 | const Standard_Real * SMA = &SMatrix->Value(1); |
161 | SMA --; |
162 | const Standard_Integer * NC = &NextCoeff->Value(1); |
163 | NC--; |
164 | |
165 | // Resolution of Sw = B; |
166 | for(i = 1; i <= RowNumber(); i++) { |
167 | DiagAddr = profile(2, i); |
168 | Sum = 0; |
169 | for(j = i - profile(1, i), jj = DiagAddr - (i - j); |
170 | j < i; j++, jj++) |
171 | Sum += SMA[jj]* x[j]; |
172 | x[i] = (b[i] - Sum)/SMA[DiagAddr]; |
173 | } |
174 | |
175 | // Resolution of t(S)X = w; |
176 | for(i = ColNumber(); i >= 1; i--) { |
177 | DiagAddr = profile(2, i); |
178 | j = NC[DiagAddr]; |
179 | Sum = 0; |
180 | while(j > 0) { |
181 | CurrAddr = profile(2, j) - (j-i); |
182 | Sum += SMA[CurrAddr]*x[j]; |
183 | j = NC[CurrAddr]; |
184 | } |
185 | x[i] = (x[i] - Sum)/SMA[DiagAddr]; |
186 | } |
187 | } |
188 | |
189 | Standard_Boolean FEmTool_ProfileMatrix::Prepare() |
190 | { |
191 | Standard_NotImplemented::Raise("FEmTool_ProfileMatrix::Prepare"); |
192 | return Standard_False; |
193 | } |
194 | |
195 | // 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 |
196 | void FEmTool_ProfileMatrix::Solve(const math_Vector& ,const math_Vector& ,math_Vector& ,math_Vector& ,const Standard_Real ,const Standard_Integer ) const |
197 | { |
198 | Standard_NotImplemented::Raise("FEmTool_ProfileMatrix::Solve"); |
199 | } |
200 | |
201 | |
202 | //======================================================================= |
203 | //function : Multiplied |
204 | //purpose : MX = H*X |
205 | //======================================================================= |
206 | void FEmTool_ProfileMatrix::Multiplied(const math_Vector& X,math_Vector& MX) const |
207 | { |
208 | Standard_Integer i, j, jj, DiagAddr, CurrAddr; |
209 | Standard_Real * m = &MX(MX.Lower()); |
210 | m--; |
211 | Standard_Real * x = &X(X.Lower()); |
212 | x--; |
213 | const Standard_Real * PM = &ProfileMatrix->Value(1); |
214 | PM--; |
215 | const Standard_Integer * NC = &NextCoeff->Value(1); |
216 | NC--; |
217 | |
218 | for(i = 1; i <= RowNumber(); i++) { |
219 | DiagAddr = profile(2, i); |
220 | m[i] = 0; |
221 | for(j = i - profile(1, i), jj = DiagAddr - (i - j); |
222 | j <= i; j++, jj++) |
223 | m[i] += PM[jj]*x[j]; |
224 | |
225 | CurrAddr = DiagAddr; |
226 | for(j = NC[CurrAddr]; j > 0; j = NC[CurrAddr]) { |
227 | CurrAddr = profile(2, j) - (j-i); |
228 | m[i] += PM[CurrAddr]*x[j]; |
229 | } |
230 | } |
231 | } |
232 | |
233 | Standard_Integer FEmTool_ProfileMatrix::RowNumber() const |
234 | { |
235 | return profile.RowLength(); |
236 | } |
237 | |
238 | Standard_Integer FEmTool_ProfileMatrix::ColNumber() const |
239 | { |
240 | return profile.RowLength(); |
241 | } |
242 | |
243 | Standard_Boolean FEmTool_ProfileMatrix::IsInProfile(const Standard_Integer i, |
244 | const Standard_Integer j) const |
245 | { |
246 | if (j <= i) { |
247 | if ((i - j) <= profile(1, i)) return Standard_True; |
248 | else return Standard_False; |
249 | } |
250 | else if ((j - i) <= profile(1, j)) return Standard_True; |
251 | |
252 | return Standard_False; |
253 | } |
254 | |
255 | void FEmTool_ProfileMatrix::OutM() const |
256 | { |
257 | Standard_Integer i, j; |
258 | cout<<"Matrix A"<<endl; |
259 | for(i = 1; i <= RowNumber(); i++) { |
260 | for(j = 1; j < i - profile(1, i); j++) |
261 | cout<<"0 "; |
262 | |
263 | for(j = profile(2, i) - profile(1, i); j <= profile(2, i); j++) |
264 | cout<<ProfileMatrix->Value(j)<<" "; |
265 | cout<<endl; |
266 | } |
267 | |
268 | cout<<"NextCoeff"<<endl; |
269 | for(i = 1; i <= profile(2, RowNumber()); i++) |
270 | cout<<NextCoeff->Value(i)<<" "; |
271 | cout<<endl; |
272 | } |
273 | |
274 | void FEmTool_ProfileMatrix::OutS() const |
275 | { |
276 | Standard_Integer i, j; |
277 | cout<<"Matrix S"<<endl; |
278 | for(i = 1; i <= RowNumber(); i++) { |
279 | for(j = 1; j < i - profile(1, i); j++) |
280 | cout<<"0 "; |
281 | |
282 | for(j = profile(2, i) - profile(1, i); j <= profile(2, i); j++) |
283 | cout<<SMatrix->Value(j)<<" "; |
284 | cout<<endl; |
285 | } |
286 | } |