1 // Created on: 2005-12-15
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.
17 #include <math_EigenValuesSearcher.hxx>
18 #include <StdFail_NotDone.hxx>
20 //==========================================================================
22 // Computation of sqrt(x*x + y*y).
23 //==========================================================================
24 static inline Standard_Real pythag(const Standard_Real x,
25 const Standard_Real y)
27 return Sqrt(x*x + y*y);
30 math_EigenValuesSearcher::math_EigenValuesSearcher(const TColStd_Array1OfReal& Diagonal,
31 const TColStd_Array1OfReal& Subdiagonal)
33 myIsDone = Standard_False;
35 Standard_Integer n = Diagonal.Length();
36 if (Subdiagonal.Length() != n)
37 Standard_Failure::Raise("math_EigenValuesSearcher : dimension mismatch");
39 myDiagonal = new TColStd_HArray1OfReal(1, n);
40 myDiagonal->ChangeArray1() = Diagonal;
41 mySubdiagonal = new TColStd_HArray1OfReal(1, n);
42 mySubdiagonal->ChangeArray1() = Subdiagonal;
44 myEigenValues = new TColStd_HArray1OfReal(1, n);
45 myEigenVectors = new TColStd_HArray2OfReal(1, n, 1, n);
47 Standard_Real* d = new Standard_Real [n+1];
48 Standard_Real* e = new Standard_Real [n+1];
49 Standard_Real** z = new Standard_Real* [n+1];
50 Standard_Integer i, j;
51 for (i = 1; i <= n; i++)
52 z[i] = new Standard_Real [n+1];
54 for (i = 1; i <= n; i++)
55 d[i] = myDiagonal->Value(i);
56 for (i = 2; i <= n; i++)
57 e[i] = mySubdiagonal->Value(i);
58 for (i = 1; i <= n; i++)
59 for (j = 1; j <= n; j++)
60 z[i][j] = (i == j)? 1. : 0.;
62 Standard_Boolean result;
65 Standard_Integer iter;
77 result = Standard_True;
82 for (i = 2; i <= n; i++)
87 for (l = 1; l <= n; l++) {
91 for (m = l; m <= n-1; m++) {
92 dd = Abs(d[m]) + Abs(d[m + 1]);
94 if (Abs(e[m]) + dd == dd)
100 result = Standard_False;
101 break; //return result;
104 g = (d[l + 1] - d[l])/(2.*e[l]);
108 g = d[m] - d[l] + e[l]/(g - r);
110 g = d[m] - d[l] + e[l]/(g + r);
116 for (i = m - 1; i >= l; i--) {
131 r = (d[i] - g)*s + 2.0*c*b;
136 for (k = 1; k <= n; k++) {
138 z[k][i + 1] = s*z[k][i] + c*f;
139 z[k][i] = c*z[k][i] - s*f;
143 if (r == 0 && i >= 1)
152 if (result == Standard_False)
154 } //end of for (l = 1; l <= n; l++)
155 } //end of if (n != 1)
159 for (i = 1; i <= n; i++)
160 myEigenValues->ChangeValue(i) = d[i];
161 for (i = 1; i <= n; i++)
162 for (j = 1; j <= n; j++)
163 myEigenVectors->ChangeValue(i, j) = z[i][j];
170 for (i = 1; i <= n; i++)
175 Standard_Boolean math_EigenValuesSearcher::IsDone() const
180 Standard_Integer math_EigenValuesSearcher::Dimension() const
185 Standard_Real math_EigenValuesSearcher::EigenValue(const Standard_Integer Index) const
187 return myEigenValues->Value(Index);
190 math_Vector math_EigenValuesSearcher::EigenVector(const Standard_Integer Index) const
192 math_Vector theVector(1, myN);
195 for (i = 1; i <= myN; i++)
196 theVector(i) = myEigenVectors->Value(i, Index);