← Back to branch summary

~ubuntu-branches/ubuntu/raring/rheolef/raring-proposed

« back to all changes in this revision

Viewing changes to nfem/geo_element/tensor-eig3.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito, Pierre Saramito, Sylvestre Ledru
  • Date: 2012-05-14 14:02:09 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20120514140209-dzbdlidkotyflf9e
Tags: 6.1-1
http://bugs.debian.org/671996
[ Pierre Saramito ]
* New upstream release 6.1 (minor changes):
  - support arbitrarily polynomial order Pk
  - source code supports g++-4.7 (closes: #671996)

[ Sylvestre Ledru ]
* update of the watch file

Show diffs side-by-side

added added

removed removed

Lines of Context:
nfem/geo_element/tensor-eig3.cc
36
36
// The source files in this directory have been copied from the public domain
37
37
// Java matrix library JAMA.  The derived source code is in the public domain
38
38
// as well.
39
 
#include <math.h>
 
39
#include "rheolef/tensor.h"
40
40
namespace rheolef { 
41
41
 
42
 
#ifdef MAX
43
 
#undef MAX
44
 
#endif
45
 
 
46
 
#define MAX(a, b) ((a)>(b)?(a):(b))
47
 
 
48
 
#define n 3
49
 
 
50
 
static double hypot2(double x, double y) {
51
 
  return sqrt(x*x+y*y);
52
 
}
 
42
template<class T>
 
43
static
 
44
inline
 
45
T hypot2(const T& x, const T& y) { return sqrt(x*x+y*y); }
53
46
 
54
47
// Symmetric Householder reduction to tridiagonal form.
55
 
 
56
 
static void tred2(double V[n][n], double d[n], double e[n]) {
 
48
template<class T>
 
49
static
 
50
void
 
51
tred2 (tensor_basic<T>& V, point_basic<T>& d, point_basic<T>& e)
 
52
{
 
53
  const int n = 3;
57
54
 
58
55
//  This is derived from the Algol procedures tred2 by
59
56
//  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
61
58
//  Fortran subroutine in EISPACK.
62
59
 
63
60
  for (int j = 0; j < n; j++) {
64
 
    d[j] = V[n-1][j];
 
61
    d[j] = V(n-1,j);
65
62
  }
66
63
 
67
64
  // Householder reduction to tridiagonal form.
70
67
 
71
68
    // Scale to avoid under/overflow.
72
69
 
73
 
    double scale = 0.0;
74
 
    double h = 0.0;
 
70
    T scale = 0.0;
 
71
    T h = 0.0;
75
72
    for (int k = 0; k < i; k++) {
76
73
      scale = scale + fabs(d[k]);
77
74
    }
78
75
    if (scale == 0.0) {
79
76
      e[i] = d[i-1];
80
77
      for (int j = 0; j < i; j++) {
81
 
        d[j] = V[i-1][j];
82
 
        V[i][j] = 0.0;
83
 
        V[j][i] = 0.0;
 
78
        d[j] = V(i-1,j);
 
79
        V(i,j) = 0.0;
 
80
        V(j,i) = 0.0;
84
81
      }
85
82
    } else {
86
83
 
90
87
        d[k] /= scale;
91
88
        h += d[k] * d[k];
92
89
      }
93
 
      double f = d[i-1];
94
 
      double g = sqrt(h);
 
90
      T f = d[i-1];
 
91
      T g = sqrt(h);
95
92
      if (f > 0) {
96
93
        g = -g;
97
94
      }
106
103
 
107
104
      for (int j = 0; j < i; j++) {
108
105
        f = d[j];
109
 
        V[j][i] = f;
110
 
        g = e[j] + V[j][j] * f;
 
106
        V(j,i) = f;
 
107
        g = e[j] + V(j,j) * f;
111
108
        for (int k = j+1; k <= i-1; k++) {
112
 
          g += V[k][j] * d[k];
113
 
          e[k] += V[k][j] * f;
 
109
          g += V(k,j) * d[k];
 
110
          e[k] += V(k,j) * f;
114
111
        }
115
112
        e[j] = g;
116
113
      }
119
116
        e[j] /= h;
120
117
        f += e[j] * d[j];
121
118
      }
122
 
      double hh = f / (h + h);
 
119
      T hh = f / (h + h);
123
120
      for (int j = 0; j < i; j++) {
124
121
        e[j] -= hh * d[j];
125
122
      }
127
124
        f = d[j];
128
125
        g = e[j];
129
126
        for (int k = j; k <= i-1; k++) {
130
 
          V[k][j] -= (f * e[k] + g * d[k]);
 
127
          V(k,j) -= (f * e[k] + g * d[k]);
131
128
        }
132
 
        d[j] = V[i-1][j];
133
 
        V[i][j] = 0.0;
 
129
        d[j] = V(i-1,j);
 
130
        V(i,j) = 0.0;
134
131
      }
135
132
    }
136
133
    d[i] = h;
139
136
  // Accumulate transformations.
140
137
 
141
138
  for (int i = 0; i < n-1; i++) {
142
 
    V[n-1][i] = V[i][i];
143
 
    V[i][i] = 1.0;
144
 
    double h = d[i+1];
 
139
    V(n-1,i) = V(i,i);
 
140
    V(i,i) = 1.0;
 
141
    T h = d[i+1];
145
142
    if (h != 0.0) {
146
143
      for (int k = 0; k <= i; k++) {
147
 
        d[k] = V[k][i+1] / h;
 
144
        d[k] = V(k,i+1) / h;
148
145
      }
149
146
      for (int j = 0; j <= i; j++) {
150
 
        double g = 0.0;
 
147
        T g = 0.0;
151
148
        for (int k = 0; k <= i; k++) {
152
 
          g += V[k][i+1] * V[k][j];
 
149
          g += V(k,i+1) * V(k,j);
153
150
        }
154
151
        for (int k = 0; k <= i; k++) {
155
 
          V[k][j] -= g * d[k];
 
152
          V(k,j) -= g * d[k];
156
153
        }
157
154
      }
158
155
    }
159
156
    for (int k = 0; k <= i; k++) {
160
 
      V[k][i+1] = 0.0;
 
157
      V(k,i+1) = 0.0;
161
158
    }
162
159
  }
163
160
  for (int j = 0; j < n; j++) {
164
 
    d[j] = V[n-1][j];
165
 
    V[n-1][j] = 0.0;
 
161
    d[j] = V(n-1,j);
 
162
    V(n-1,j) = 0.0;
166
163
  }
167
 
  V[n-1][n-1] = 1.0;
 
164
  V(n-1,n-1) = 1.0;
168
165
  e[0] = 0.0;
169
166
} 
170
 
 
171
167
// Symmetric tridiagonal QL algorithm.
172
 
 
173
 
static void tql2(double V[n][n], double d[n], double e[n]) {
 
168
template<class T>
 
169
static
 
170
void
 
171
tql2 (tensor_basic<T>& V, point_basic<T>& d, point_basic<T>& e)
 
172
{
 
173
  const int n = 3;
174
174
 
175
175
//  This is derived from the Algol procedures tql2, by
176
176
//  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
182
182
  }
183
183
  e[n-1] = 0.0;
184
184
 
185
 
  double f = 0.0;
186
 
  double tst1 = 0.0;
187
 
  double eps = pow(2.0,-52.0);
 
185
  T f = 0.0;
 
186
  T tst1 = 0.0;
 
187
  T eps = pow(2.0,-52.0);
188
188
  for (int l = 0; l < n; l++) {
189
189
 
190
190
    // Find small subdiagonal element
191
191
 
192
 
    tst1 = MAX(tst1,fabs(d[l]) + fabs(e[l]));
 
192
    tst1 = std::max (tst1, fabs(d[l]) + fabs(e[l]));
193
193
    int m = l;
194
194
    while (m < n) {
195
195
      if (fabs(e[m]) <= eps*tst1) {
208
208
 
209
209
        // Compute implicit shift
210
210
 
211
 
        double g = d[l];
212
 
        double p = (d[l+1] - g) / (2.0 * e[l]);
213
 
        double r = hypot2(p,1.0);
 
211
        T g = d[l];
 
212
        T p = (d[l+1] - g) / (2.0 * e[l]);
 
213
        T r = hypot2(p,1.0);
214
214
        if (p < 0) {
215
215
          r = -r;
216
216
        }
217
217
        d[l] = e[l] / (p + r);
218
218
        d[l+1] = e[l] * (p + r);
219
 
        double dl1 = d[l+1];
220
 
        double h = g - d[l];
 
219
        T dl1 = d[l+1];
 
220
        T h = g - d[l];
221
221
        for (int i = l+2; i < n; i++) {
222
222
          d[i] -= h;
223
223
        }
226
226
        // Implicit QL transformation.
227
227
 
228
228
        p = d[m];
229
 
        double c = 1.0;
230
 
        double c2 = c;
231
 
        double c3 = c;
232
 
        double el1 = e[l+1];
233
 
        double s = 0.0;
234
 
        double s2 = 0.0;
 
229
        T c = 1.0;
 
230
        T c2 = c;
 
231
        T c3 = c;
 
232
        T el1 = e[l+1];
 
233
        T s = 0.0;
 
234
        T s2 = 0.0;
235
235
        for (int i = m-1; i >= l; i--) {
236
236
          c3 = c2;
237
237
          c2 = c;
248
248
          // Accumulate transformation.
249
249
 
250
250
          for (int k = 0; k < n; k++) {
251
 
            h = V[k][i+1];
252
 
            V[k][i+1] = s * V[k][i] + c * h;
253
 
            V[k][i] = c * V[k][i] - s * h;
 
251
            h = V(k,i+1);
 
252
            V(k,i+1) = s * V(k,i) + c * h;
 
253
            V(k,i) = c * V(k,i) - s * h;
254
254
          }
255
255
        }
256
256
        p = -s * s2 * c3 * el1 * e[l] / dl1;
269
269
 
270
270
  for (int i = 0; i < n-1; i++) {
271
271
    int k = i;
272
 
    double p = d[i];
 
272
    T p = d[i];
273
273
    for (int j = i+1; j < n; j++) {
274
274
      if (d[j] > p) {
275
275
        k = j;
280
280
      d[k] = d[i];
281
281
      d[i] = p;
282
282
      for (int j = 0; j < n; j++) {
283
 
        p = V[j][i];
284
 
        V[j][i] = V[j][k];
285
 
        V[j][k] = p;
 
283
        p = V(j,i);
 
284
        V(j,i) = V(j,k);
 
285
        V(j,k) = p;
286
286
      }
287
287
    }
288
288
  }
289
289
}
290
 
 
291
 
void eigen_decomposition(double A[n][n], double V[n][n], double d[n]) {
292
 
  double e[n];
293
 
  for (int i = 0; i < n; i++) {
294
 
    for (int j = 0; j < n; j++) {
295
 
      V[i][j] = A[i][j];
296
 
    }
297
 
  }
298
 
  tred2(V, d, e);
299
 
  tql2(V, d, e);
 
290
template<class T>
 
291
void
 
292
eigen_decomposition (const tensor_basic<T>& A, tensor_basic<T>& V, point_basic<T>& d)
 
293
{
 
294
  V = A;
 
295
  point_basic<T> e;
 
296
  tred2 (V, d, e);
 
297
  tql2  (V, d, e);
300
298
}
 
299
// ----------------------------------------------------------------------------
 
300
// instanciation in library
 
301
// ----------------------------------------------------------------------------
 
302
#define _RHEOLEF_instanciation(T)                                               \
 
303
template void eigen_decomposition (const tensor_basic<T>& A, tensor_basic<T>& V, point_basic<T>& d);
 
304
 
 
305
_RHEOLEF_instanciation(Float)
 
306
 
301
307
}// namespace rheolef