src/math/math_EigenValuesSearcher.cxx

   1 // Created on: 2005-12-15
   2 // Created by: Julia GERASIMOVA
   3 // Copyright (c) 2005-2014 OPEN CASCADE SAS
   4 //
   5 // This file is part of Open CASCADE Technology software library.
   6 //
   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.
  12 //
  13 // Alternatively, this file may be used under the terms of Open CASCADE
  14 // commercial license or contractual agreement.
  15
  16
  17 #include <math_EigenValuesSearcher.hxx>
  18 #include <StdFail_NotDone.hxx>
  19
  20 //==========================================================================
  21 //function : pythag
  22 //           Computation of sqrt(x*x + y*y).
  23 //==========================================================================
  24 static inline Standard_Real pythag(const Standard_Real x,
  25                                    const Standard_Real y)
  26 {
  27   return Sqrt(x*x + y*y);
  28 }
  29
  30 math_EigenValuesSearcher::math_EigenValuesSearcher(const TColStd_Array1OfReal& Diagonal,
  31                                                    const TColStd_Array1OfReal& Subdiagonal)
  32 {
  33   myIsDone = Standard_False;
  34
  35   Standard_Integer n = Diagonal.Length();
  36   if (Subdiagonal.Length() != n)
  37     Standard_Failure::Raise("math_EigenValuesSearcher : dimension mismatch");
  38
  39   myDiagonal    = new TColStd_HArray1OfReal(1, n);
  40   myDiagonal->ChangeArray1() = Diagonal;
  41   mySubdiagonal = new TColStd_HArray1OfReal(1, n);
  42   mySubdiagonal->ChangeArray1() = Subdiagonal;
  43   myN = n;
  44   myEigenValues  = new TColStd_HArray1OfReal(1, n);
  45   myEigenVectors = new TColStd_HArray2OfReal(1, n, 1, n);
  46
  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];
  53
  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.;
  61
  62   Standard_Boolean result;
  63   Standard_Integer m;
  64   Standard_Integer l;
  65   Standard_Integer iter;
  66   //Standard_Integer i;
  67   Standard_Integer k;
  68   Standard_Real    s;
  69   Standard_Real    r;
  70   Standard_Real    p;
  71   Standard_Real    g;
  72   Standard_Real    f;
  73   Standard_Real    dd;
  74   Standard_Real    c;
  75   Standard_Real    b;
  76
  77   result = Standard_True;
  78
  79   if (n != 1)
  80     {
  81       // Shift e.
  82       for (i = 2; i <= n; i++)
  83         e[i - 1] = e[i];
  84
  85       e[n] = 0.0;
  86
  87       for (l = 1; l <= n; l++) {
  88         iter = 0;
  89
  90         do {
  91           for (m = l; m <= n-1; m++) {
  92             dd = Abs(d[m]) + Abs(d[m + 1]);
  93
  94             if (Abs(e[m]) + dd == dd)
  95               break;
  96           }
  97
  98           if (m != l) {
  99             if (iter++ == 30) {
 100               result = Standard_False;
 101               break; //return result;
 102             }
 103
 104             g = (d[l + 1] - d[l])/(2.*e[l]);
 105             r = pythag(1., g);
 106
 107             if (g < 0)
 108               g = d[m] - d[l] + e[l]/(g - r);
 109             else
 110               g = d[m] - d[l] + e[l]/(g + r);
 111
 112             s = 1.;
 113             c = 1.;
 114             p = 0.;
 115
 116             for (i = m - 1; i >= l; i--) {
 117               f        = s*e[i];
 118               b        = c*e[i];
 119               r        = pythag (f, g);
 120               e[i + 1] = r;
 121
 122               if (r == 0.) {
 123                 d[i + 1] -= p;
 124                 e[m]      = 0.;
 125                 break;
 126               }
 127
 128               s        = f/r;
 129               c        = g/r;
 130               g        = d[i + 1] - p;
 131               r        = (d[i] - g)*s + 2.0*c*b;
 132               p        = s*r;
 133               d[i + 1] = g + p;
 134               g        = c*r - b;
 135
 136               for (k = 1; k <= n; k++) {
 137                 f           = z[k][i + 1];
 138                 z[k][i + 1] = s*z[k][i] + c*f;
 139                 z[k][i]     = c*z[k][i] - s*f;
 140               }
 141             }
 142
 143             if (r == 0 && i >= 1)
 144               continue;
 145
 146             d[l] -= p;
 147             e[l]  = g;
 148             e[m]  = 0.;
 149           }
 150         }
 151         while (m != l);
 152         if (result == Standard_False)
 153           break;
 154       } //end of for (l = 1; l <= n; l++)
 155     } //end of if (n != 1)
 156
 157   if (result)
 158     {
 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];
 164     }
 165
 166   myIsDone = result;
 167
 168   delete [] d;
 169   delete [] e;
 170   for (i = 1; i <= n; i++)
 171     delete [] z[i];
 172   delete [] z;
 173 }
 174
 175 Standard_Boolean math_EigenValuesSearcher::IsDone() const
 176 {
 177   return myIsDone;
 178 }
 179
 180 Standard_Integer math_EigenValuesSearcher::Dimension() const
 181 {
 182   return myN;
 183 }
 184
 185 Standard_Real math_EigenValuesSearcher::EigenValue(const Standard_Integer Index) const
 186 {
 187   return myEigenValues->Value(Index);
 188 }
 189
 190 math_Vector math_EigenValuesSearcher::EigenVector(const Standard_Integer Index) const
 191 {
 192   math_Vector theVector(1, myN);
 193
 194   Standard_Integer i;
 195   for (i = 1; i <= myN; i++)
 196     theVector(i) = myEigenVectors->Value(i, Index);
 197
 198   return theVector;
 199 }