b311480e |
1 | // Created on: 2005-12-21 |
2 | // Created by: Julia GERASIMOVA |
973c2be1 |
3 | // Copyright (c) 2005-2014 OPEN CASCADE SAS |
b311480e |
4 | // |
973c2be1 |
5 | // This file is part of Open CASCADE Technology software library. |
b311480e |
6 | // |
d5f74e42 |
7 | // This library is free software; you can redistribute it and/or modify it under |
8 | // the terms of the GNU Lesser General Public License version 2.1 as published |
973c2be1 |
9 | // by the Free Software Foundation, with special exception defined in the file |
10 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
11 | // distribution for complete text of the license and disclaimer of any warranty. |
b311480e |
12 | // |
973c2be1 |
13 | // Alternatively, this file may be used under the terms of Open CASCADE |
14 | // commercial license or contractual agreement. |
7fd59977 |
15 | |
42cf5bc1 |
16 | |
7fd59977 |
17 | #include <math_Array1OfValueAndWeight.hxx> |
42cf5bc1 |
18 | #include <math_ComputeKronrodPointsAndWeights.hxx> |
19 | #include <math_EigenValuesSearcher.hxx> |
7fd59977 |
20 | #include <Standard_ErrorHandler.hxx> |
21 | |
e35db416 |
22 | #include <algorithm> |
7fd59977 |
23 | math_ComputeKronrodPointsAndWeights::math_ComputeKronrodPointsAndWeights(const Standard_Integer Number) |
24 | { |
25 | myIsDone = Standard_False; |
26 | |
27 | try { |
28 | Standard_Integer i, j; |
29 | Standard_Integer a2NP1 = 2*Number + 1; |
30 | |
31 | myPoints = new TColStd_HArray1OfReal(1, a2NP1); |
32 | myWeights = new TColStd_HArray1OfReal(1, a2NP1); |
33 | |
34 | TColStd_Array1OfReal aDiag(1, a2NP1); |
35 | TColStd_Array1OfReal aSubDiag(1, a2NP1); |
36 | |
37 | // Initialize symmetric tridiagonal matrix. |
38 | Standard_Integer n = Number; |
39 | Standard_Integer aKronrodN = 2*Number + 1; |
40 | Standard_Integer a3KN2p1 = Min(3*(Number + 1)/2 + 1, aKronrodN); |
41 | for (i = 1; i <= a3KN2p1; i++) { |
42 | aDiag(i) = 0.; |
43 | |
44 | if (i == 1) |
45 | aSubDiag(i) = 0.; |
46 | else { |
47 | Standard_Integer sqrIm1 = (i-1)*(i-1); |
48 | aSubDiag(i) = sqrIm1/(4.*sqrIm1 - 1); |
49 | } |
50 | } |
51 | |
52 | for (i = a3KN2p1 + 1; i <= aKronrodN; i++) { |
53 | aDiag(i) = 0.; |
54 | aSubDiag(i) = 0.; |
55 | } |
56 | |
57 | // Initialization of temporary data structures. |
58 | Standard_Integer aNd2 = Number/2; |
59 | Standard_Real *s = new Standard_Real[aNd2 + 2]; |
60 | Standard_Real *t = new Standard_Real[aNd2 + 2]; |
61 | Standard_Real *ss = s++; |
62 | Standard_Real *tt = t++; |
63 | |
64 | for (i = -1; i <= aNd2; i++) { |
65 | s[i] = 0.; |
66 | t[i] = 0.; |
67 | } |
68 | |
69 | // Generation of Jacobi-Kronrod matrix. |
70 | Standard_Real *aa = new Standard_Real [a2NP1+1]; |
71 | Standard_Real *bb = new Standard_Real [a2NP1+1]; |
72 | for (i = 1; i <= a2NP1; i++) { |
73 | aa[i] = aDiag(i); |
74 | bb[i] = aSubDiag(i); |
75 | } |
76 | Standard_Real *ptrtmp; |
77 | Standard_Real u; |
78 | Standard_Integer m; |
79 | Standard_Integer k; |
80 | Standard_Integer l; |
81 | |
82 | Standard_Real *a = aa+1; |
83 | Standard_Real *b = bb+1; |
84 | |
85 | // Eastward phase. |
86 | t[0] = b[Number + 1]; |
87 | |
88 | for (m = 0; m <= n - 2; m++) { |
89 | u = 0; |
90 | |
91 | for (k = (m + 1)/2; k >= 0; k--) { |
92 | l = m - k; |
93 | u += (a[k + n + 1] - a[l])*t[k] + b[k + n + 1]*s[k - 1] - b[l]*s[k]; |
94 | s[k] = u; |
95 | } |
96 | |
97 | ptrtmp = t; |
98 | t = s; |
99 | s = ptrtmp; |
100 | } |
101 | |
102 | for (j = aNd2; j >= 0; j--) |
103 | s[j] = s[j - 1]; |
104 | |
105 | // Southward phase. |
106 | for (m = n - 1; m <= 2*n - 3; m++) { |
107 | u = 0; |
108 | |
109 | for (k = m + 1 - n; k <= (m - 1)/2; k++) { |
110 | l = m - k; |
111 | j = n - 1 - l; |
112 | u += -(a[k + n + 1] - a[l])*t[j] - b[k + n + 1]*s[j] + b[l]*s[j + 1]; |
113 | s[j] = u; |
114 | } |
115 | |
116 | if (m % 2 == 0) { |
117 | k = m/2; |
118 | a[k + n + 1] = a[k] + (s[j] - b[k + n + 1]*s[j + 1])/ t[j + 1]; |
119 | } else { |
120 | k = (m + 1)/2; |
121 | b[k + n + 1] = s[j]/s[j + 1]; |
122 | } |
123 | |
124 | ptrtmp = t; |
125 | t = s; |
126 | s = ptrtmp; |
127 | } |
128 | |
129 | // Termination phase. |
130 | a[2*Number] = a[n - 1] - b[2*Number]*s[0]/t[0]; |
131 | |
132 | delete [] ss; |
133 | delete [] tt; |
134 | for (i = 1; i <= a2NP1; i++) { |
135 | aDiag(i) = aa[i]; |
136 | aSubDiag(i) = bb[i]; |
137 | } |
138 | delete [] aa; |
139 | delete [] bb; |
140 | |
141 | for (i = 1; i <= a2NP1; i++) |
142 | aSubDiag(i) = Sqrt(aSubDiag(i)); |
143 | |
144 | // Compute eigen values. |
145 | math_EigenValuesSearcher EVsearch(aDiag, aSubDiag); |
146 | |
147 | if (EVsearch.IsDone()) { |
148 | math_Array1OfValueAndWeight VWarray(1, a2NP1); |
149 | for (i = 1; i <= a2NP1; i++) { |
150 | math_Vector anEigenVector = EVsearch.EigenVector(i); |
151 | Standard_Real aWeight = anEigenVector(1); |
152 | aWeight = 2. * aWeight * aWeight; |
153 | math_ValueAndWeight EVW( EVsearch.EigenValue(i), aWeight ); |
154 | VWarray(i) = EVW; |
155 | } |
156 | |
e35db416 |
157 | std::sort (VWarray.begin(), VWarray.end()); |
7fd59977 |
158 | |
159 | for (i = 1; i <= a2NP1; i++) { |
160 | myPoints->ChangeValue(i) = VWarray(i).Value(); |
161 | myWeights->ChangeValue(i) = VWarray(i).Weight(); |
162 | } |
163 | myIsDone = Standard_True; |
164 | } |
a738b534 |
165 | } catch (Standard_Failure const&) { |
7fd59977 |
166 | } |
167 | } |
168 | |
169 | Standard_Boolean math_ComputeKronrodPointsAndWeights::IsDone() const |
170 | { |
171 | return myIsDone; |
172 | } |
173 | |
174 | math_Vector math_ComputeKronrodPointsAndWeights::Points() const |
175 | { |
176 | Standard_Integer Number = myPoints->Length(); |
177 | math_Vector thePoints(1, Number); |
178 | for (Standard_Integer i = 1; i <= Number; i++) |
179 | thePoints(i) = myPoints->Value(i); |
180 | |
181 | return thePoints; |
182 | } |
183 | |
184 | math_Vector math_ComputeKronrodPointsAndWeights::Weights() const |
185 | { |
186 | Standard_Integer Number = myWeights->Length(); |
187 | math_Vector theWeights(1, Number); |
188 | for (Standard_Integer i = 1; i <= Number; i++) |
189 | theWeights(i) = myWeights->Value(i); |
190 | |
191 | return theWeights; |
192 | } |