src/math/math_ComputeKronrodPointsAndWeights.cxx

   1 // Created on: 2005-12-21
   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 #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>
  22
  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
 157       math_CompareOfValueAndWeight theComparator;
 158       math_QuickSortOfValueAndWeight::Sort(VWarray, theComparator);
 159
 160       for (i = 1; i <= a2NP1; i++) {
 161         myPoints->ChangeValue(i)  = VWarray(i).Value();
 162         myWeights->ChangeValue(i) = VWarray(i).Weight();
 163       }
 164       myIsDone = Standard_True;
 165     }
 166   } catch (Standard_Failure) {
 167   }
 168 }
 169
 170 Standard_Boolean math_ComputeKronrodPointsAndWeights::IsDone() const
 171 {
 172   return myIsDone;
 173 }
 174
 175 math_Vector math_ComputeKronrodPointsAndWeights::Points() const
 176 {
 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);
 181
 182   return thePoints;
 183 }
 184
 185 math_Vector math_ComputeKronrodPointsAndWeights::Weights() const
 186 {
 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);
 191
 192   return theWeights;
 193 }