src/math/math_EigenValuesSearcher.cxx

   1 // Created on: 2005-12-15
   2 // Created by: Julia GERASIMOVA
   3 // Copyright (c) 2005-2012 OPEN CASCADE SAS
   4 //
   5 // The content of this file is subject to the Open CASCADE Technology Public
   6 // License Version 6.5 (the "License"). You may not use the content of this file
   7 // except in compliance with the License. Please obtain a copy of the License
   8 // at http://www.opencascade.org and read it completely before using this file.
   9 //
  10 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  11 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  12 //
  13 // The Original Code and all software distributed under the License is
  14 // distributed on an "AS IS" basis, without warranty of any kind, and the
  15 // Initial Developer hereby disclaims all such warranties, including without
  16 // limitation, any warranties of merchantability, fitness for a particular
  17 // purpose or non-infringement. Please see the License for the specific terms
  18 // and conditions governing the rights and limitations under the License.
  19
  20
  21
  22 #include <math_EigenValuesSearcher.ixx>
  23
  24 //==========================================================================
  25 //function : pythag
  26 //           Computation of sqrt(x*x + y*y).
  27 //==========================================================================
  28
  29 static inline Standard_Real pythag(const Standard_Real x,
  30                                    const Standard_Real y)
  31 {
  32   return Sqrt(x*x + y*y);
  33 }
  34
  35 math_EigenValuesSearcher::math_EigenValuesSearcher(const TColStd_Array1OfReal& Diagonal,
  36                                                    const TColStd_Array1OfReal& Subdiagonal)
  37 {
  38   myIsDone = Standard_False;
  39
  40   Standard_Integer n = Diagonal.Length();
  41   if (Subdiagonal.Length() != n)
  42     Standard_Failure::Raise("math_EigenValuesSearcher : dimension mismatch");
  43
  44   myDiagonal    = new TColStd_HArray1OfReal(1, n);
  45   myDiagonal->ChangeArray1() = Diagonal;
  46   mySubdiagonal = new TColStd_HArray1OfReal(1, n);
  47   mySubdiagonal->ChangeArray1() = Subdiagonal;
  48   myN = n;
  49   myEigenValues  = new TColStd_HArray1OfReal(1, n);
  50   myEigenVectors = new TColStd_HArray2OfReal(1, n, 1, n);
  51
  52   Standard_Real* d  = new Standard_Real [n+1];
  53   Standard_Real* e  = new Standard_Real [n+1];
  54   Standard_Real** z = new Standard_Real* [n+1];
  55   Standard_Integer i, j;
  56   for (i = 1; i <= n; i++)
  57     z[i] = new Standard_Real [n+1];
  58
  59   for (i = 1; i <= n; i++)
  60     d[i] = myDiagonal->Value(i);
  61   for (i = 2; i <= n; i++)
  62     e[i] = mySubdiagonal->Value(i);
  63   for (i = 1; i <= n; i++)
  64     for (j = 1; j <= n; j++)
  65       z[i][j] = (i == j)? 1. : 0.;
  66
  67   Standard_Boolean result;
  68   Standard_Integer m;
  69   Standard_Integer l;
  70   Standard_Integer iter;
  71   //Standard_Integer i;
  72   Standard_Integer k;
  73   Standard_Real    s;
  74   Standard_Real    r;
  75   Standard_Real    p;
  76   Standard_Real    g;
  77   Standard_Real    f;
  78   Standard_Real    dd;
  79   Standard_Real    c;
  80   Standard_Real    b;
  81
  82   result = Standard_True;
  83
  84   if (n != 1)
  85     {
  86       // Shift e.
  87       for (i = 2; i <= n; i++)
  88         e[i - 1] = e[i];
  89
  90       e[n] = 0.0;
  91
  92       for (l = 1; l <= n; l++) {
  93         iter = 0;
  94
  95         do {
  96           for (m = l; m <= n-1; m++) {
  97             dd = Abs(d[m]) + Abs(d[m + 1]);
  98
  99             if (Abs(e[m]) + dd == dd)
 100               break;
 101           }
 102
 103           if (m != l) {
 104             if (iter++ == 30) {
 105               result = Standard_False;
 106               break; //return result;
 107             }
 108
 109             g = (d[l + 1] - d[l])/(2.*e[l]);
 110             r = pythag(1., g);
 111
 112             if (g < 0)
 113               g = d[m] - d[l] + e[l]/(g - r);
 114             else
 115               g = d[m] - d[l] + e[l]/(g + r);
 116
 117             s = 1.;
 118             c = 1.;
 119             p = 0.;
 120
 121             for (i = m - 1; i >= l; i--) {
 122               f        = s*e[i];
 123               b        = c*e[i];
 124               r        = pythag (f, g);
 125               e[i + 1] = r;
 126
 127               if (r == 0.) {
 128                 d[i + 1] -= p;
 129                 e[m]      = 0.;
 130                 break;
 131               }
 132
 133               s        = f/r;
 134               c        = g/r;
 135               g        = d[i + 1] - p;
 136               r        = (d[i] - g)*s + 2.0*c*b;
 137               p        = s*r;
 138               d[i + 1] = g + p;
 139               g        = c*r - b;
 140
 141               for (k = 1; k <= n; k++) {
 142                 f           = z[k][i + 1];
 143                 z[k][i + 1] = s*z[k][i] + c*f;
 144                 z[k][i]     = c*z[k][i] - s*f;
 145               }
 146             }
 147
 148             if (r == 0 && i >= 1)
 149               continue;
 150
 151             d[l] -= p;
 152             e[l]  = g;
 153             e[m]  = 0.;
 154           }
 155         }
 156         while (m != l);
 157         if (result == Standard_False)
 158           break;
 159       } //end of for (l = 1; l <= n; l++)
 160     } //end of if (n != 1)
 161
 162   if (result)
 163     {
 164       for (i = 1; i <= n; i++)
 165         myEigenValues->ChangeValue(i) = d[i];
 166       for (i = 1; i <= n; i++)
 167         for (j = 1; j <= n; j++)
 168           myEigenVectors->ChangeValue(i, j) = z[i][j];
 169     }
 170
 171   myIsDone = result;
 172
 173   delete [] d;
 174   delete [] e;
 175   for (i = 1; i <= n; i++)
 176     delete [] z[i];
 177   delete [] z;
 178 }
 179
 180 Standard_Boolean math_EigenValuesSearcher::IsDone() const
 181 {
 182   return myIsDone;
 183 }
 184
 185 Standard_Integer math_EigenValuesSearcher::Dimension() const
 186 {
 187   return myN;
 188 }
 189
 190 Standard_Real math_EigenValuesSearcher::EigenValue(const Standard_Integer Index) const
 191 {
 192   return myEigenValues->Value(Index);
 193 }
 194
 195 math_Vector math_EigenValuesSearcher::EigenVector(const Standard_Integer Index) const
 196 {
 197   math_Vector theVector(1, myN);
 198
 199   Standard_Integer i;
 200   for (i = 1; i <= myN; i++)
 201     theVector(i) = myEigenVectors->Value(i, Index);
 202
 203   return theVector;
 204 }