1 // Created on: 2005-12-21
2 // Created by: Julia GERASIMOVA
3 // Copyright (c) 2005-2014 OPEN CASCADE SAS
5 // This file is part of Open CASCADE Technology software library.
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
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.
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
16 #include <math_ComputeKronrodPointsAndWeights.ixx>
17 #include <math_EigenValuesSearcher.hxx>
18 #include <math_Array1OfValueAndWeight.hxx>
19 #include <math_CompareOfValueAndWeight.hxx>
20 #include <math_QuickSortOfValueAndWeight.hxx>
21 #include <Standard_ErrorHandler.hxx>
23 math_ComputeKronrodPointsAndWeights::math_ComputeKronrodPointsAndWeights(const Standard_Integer Number)
25 myIsDone = Standard_False;
28 Standard_Integer i, j;
29 Standard_Integer a2NP1 = 2*Number + 1;
31 myPoints = new TColStd_HArray1OfReal(1, a2NP1);
32 myWeights = new TColStd_HArray1OfReal(1, a2NP1);
34 TColStd_Array1OfReal aDiag(1, a2NP1);
35 TColStd_Array1OfReal aSubDiag(1, a2NP1);
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++) {
47 Standard_Integer sqrIm1 = (i-1)*(i-1);
48 aSubDiag(i) = sqrIm1/(4.*sqrIm1 - 1);
52 for (i = a3KN2p1 + 1; i <= aKronrodN; i++) {
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++;
64 for (i = -1; i <= aNd2; i++) {
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++) {
76 Standard_Real *ptrtmp;
82 Standard_Real *a = aa+1;
83 Standard_Real *b = bb+1;
88 for (m = 0; m <= n - 2; m++) {
91 for (k = (m + 1)/2; k >= 0; k--) {
93 u += (a[k + n + 1] - a[l])*t[k] + b[k + n + 1]*s[k - 1] - b[l]*s[k];
102 for (j = aNd2; j >= 0; j--)
106 for (m = n - 1; m <= 2*n - 3; m++) {
109 for (k = m + 1 - n; k <= (m - 1)/2; k++) {
112 u += -(a[k + n + 1] - a[l])*t[j] - b[k + n + 1]*s[j] + b[l]*s[j + 1];
118 a[k + n + 1] = a[k] + (s[j] - b[k + n + 1]*s[j + 1])/ t[j + 1];
121 b[k + n + 1] = s[j]/s[j + 1];
129 // Termination phase.
130 a[2*Number] = a[n - 1] - b[2*Number]*s[0]/t[0];
134 for (i = 1; i <= a2NP1; i++) {
141 for (i = 1; i <= a2NP1; i++)
142 aSubDiag(i) = Sqrt(aSubDiag(i));
144 // Compute eigen values.
145 math_EigenValuesSearcher EVsearch(aDiag, aSubDiag);
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 );
157 math_CompareOfValueAndWeight theComparator;
158 math_QuickSortOfValueAndWeight::Sort(VWarray, theComparator);
160 for (i = 1; i <= a2NP1; i++) {
161 myPoints->ChangeValue(i) = VWarray(i).Value();
162 myWeights->ChangeValue(i) = VWarray(i).Weight();
164 myIsDone = Standard_True;
166 } catch (Standard_Failure) {
170 Standard_Boolean math_ComputeKronrodPointsAndWeights::IsDone() const
175 math_Vector math_ComputeKronrodPointsAndWeights::Points() const
177 Standard_Integer Number = myPoints->Length();
178 math_Vector thePoints(1, Number);
179 for (Standard_Integer i = 1; i <= Number; i++)
180 thePoints(i) = myPoints->Value(i);
185 math_Vector math_ComputeKronrodPointsAndWeights::Weights() const
187 Standard_Integer Number = myWeights->Length();
188 math_Vector theWeights(1, Number);
189 for (Standard_Integer i = 1; i <= Number; i++)
190 theWeights(i) = myWeights->Value(i);