0026377: Passing Handle objects as arguments to functions as non-const reference...
[occt.git] / src / FEmTool / FEmTool_ProfileMatrix.cxx
CommitLineData
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 29IMPLEMENT_STANDARD_RTTIEXT(FEmTool_ProfileMatrix,FEmTool_SparseMatrix)
30
7fd59977 31//=======================================================================
32//function : :FEmTool_ProfileMatrix
33//purpose :
34//=======================================================================
b311480e 35FEmTool_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
243Standard_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}