← Back to branch summary

~ubuntu-branches/debian/sid/octave-tisean/sid

« back to all changes in this revision

Viewing changes to src/routines_c/eigen.cc

  • Committer: Package Import Robot
  • Author(s): Rafael Laboissiere
  • Date: 2017-08-14 12:53:47 UTC
  • Revision ID: package-import@ubuntu.com-20170814125347-ju5owr4dggr53a2n
Tags: upstream-0.2.3
Import upstream version 0.2.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/routines_c/eigen.cc
 
1
/* Modified: Piotr Held <pjheld@gmail.com> (2015). 
 
2
 * These functions are based on eigen.c of TISEAN 3.0.1 https://github.com/heggus/Tisean"
 
3
 */
 
4
/********************************************************************/
 
5
/********************************************************************/
 
6
 
 
7
#include <octave/oct.h>
 
8
 
 
9
#include <cmath>
 
10
#include <cstdlib>
 
11
#include <cstdio>
 
12
 
 
13
typedef double doublereal;
 
14
typedef int integer;
 
15
 
 
16
#define abs(x) (((x)>=0.0)?(x):-(x))
 
17
#define min(x,y) (((x)<=(y))?(x):(y))
 
18
#define max(x,y) (((x)>=(y))?(x):(y))
 
19
 
 
20
static doublereal c_b10 = 1.;
 
21
 
 
22
extern void check_alloc(void*);
 
23
 
 
24
double d_sign(double *a,double *b)
 
25
{
 
26
  double x;
 
27
  x = (*a >= 0 ? *a : - *a);
 
28
  return ( *b >= 0 ? x : -x);
 
29
}
 
30
 
 
31
doublereal pythag(doublereal *a, doublereal *b)
 
32
{
 
33
    doublereal ret_val, d__1, d__2, d__3;
 
34
    static doublereal p, r__, s, t, u;
 
35
 
 
36
    d__1 = abs(*a), d__2 = abs(*b);
 
37
    p = max(d__1,d__2);
 
38
    if (p == 0.) {
 
39
        goto L20;
 
40
    }
 
41
    d__2 = abs(*a), d__3 = abs(*b);
 
42
    d__1 = min(d__2,d__3) / p;
 
43
    r__ = d__1 * d__1;
 
44
L10:
 
45
    t = r__ + 4.;
 
46
    if (t == 4.) {
 
47
        goto L20;
 
48
    }
 
49
    s = r__ / t;
 
50
    u = s * 2. + 1.;
 
51
    p = u * p;
 
52
    d__1 = s / u;
 
53
    r__ = d__1 * d__1 * r__;
 
54
    goto L10;
 
55
L20:
 
56
    ret_val = p;
 
57
    return ret_val;
 
58
}
 
59
 
 
60
 
 
61
int tred2(const integer *nm, const integer *n, doublereal *a, 
 
62
        doublereal *d__, doublereal *e, doublereal *z__)
 
63
{
 
64
    integer a_dim1, a_offset, z_dim1, z_offset, i__1, i__2, i__3;
 
65
    doublereal d__1;
 
66
 
 
67
    double sqrt(doublereal), d_sign(doublereal *, doublereal *);
 
68
 
 
69
    static doublereal f, g, h__;
 
70
    static integer i__, j, k, l;
 
71
    static doublereal hh;
 
72
    static integer ii, jp1;
 
73
    static doublereal scale;
 
74
 
 
75
 
 
76
 
 
77
/*     this subroutine is a translation of the algol procedure tred2, */
 
78
/*     num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson. */
 
79
/*     handbook for auto. comp., vol.ii-linear algebra, 212-226(1971). */
 
80
 
 
81
/*     this subroutine reduces a real symmetric matrix to a */
 
82
/*     symmetric tridiagonal matrix using and accumulating */
 
83
/*     orthogonal similarity transformations. */
 
84
 
 
85
/*     on input */
 
86
 
 
87
/*        nm must be set to the row dimension of two-dimensional */
 
88
/*          array parameters as declared in the calling program */
 
89
/*          dimension statement. */
 
90
 
 
91
/*        n is the order of the matrix. */
 
92
 
 
93
/*        a contains the real symmetric input matrix.  only the */
 
94
/*          lower triangle of the matrix need be supplied. */
 
95
 
 
96
/*     on output */
 
97
 
 
98
/*        d contains the diagonal elements of the tridiagonal matrix. */
 
99
 
 
100
/*        e contains the subdiagonal elements of the tridiagonal */
 
101
/*          matrix in its last n-1 positions.  e(1) is set to zero. */
 
102
 
 
103
/*        z contains the orthogonal transformation matrix */
 
104
/*          produced in the reduction. */
 
105
 
 
106
/*        a and z may coincide.  if distinct, a is unaltered. */
 
107
 
 
108
/*     questions and comments should be directed to burton s. garbow, */
 
109
/*     mathematics and computer science div, argonne national laboratory */
 
110
 
 
111
/*     this version dated august 1983. */
 
112
 
 
113
/*     ------------------------------------------------------------------ */
 
114
 
 
115
    z_dim1 = *nm;
 
116
    z_offset = 1 + z_dim1 * 1;
 
117
    z__ -= z_offset;
 
118
    --e;
 
119
    --d__;
 
120
    a_dim1 = *nm;
 
121
    a_offset = 1 + a_dim1 * 1;
 
122
    a -= a_offset;
 
123
 
 
124
    i__1 = *n;
 
125
    for (i__ = 1; i__ <= i__1; ++i__) {
 
126
 
 
127
        i__2 = *n;
 
128
        for (j = i__; j <= i__2; ++j) {
 
129
            z__[j + i__ * z_dim1] = a[j + i__ * a_dim1];
 
130
        }
 
131
 
 
132
        d__[i__] = a[*n + i__ * a_dim1];
 
133
    }
 
134
 
 
135
    if (*n == 1) {
 
136
        goto L510;
 
137
    }
 
138
    i__1 = *n;
 
139
    for (ii = 2; ii <= i__1; ++ii) {
 
140
        i__ = *n + 2 - ii;
 
141
        l = i__ - 1;
 
142
        h__ = 0.;
 
143
        scale = 0.;
 
144
        if (l < 2) {
 
145
            goto L130;
 
146
        }
 
147
        i__2 = l;
 
148
        for (k = 1; k <= i__2; ++k) {
 
149
            scale += (d__1 = d__[k], abs(d__1));
 
150
        }
 
151
 
 
152
        if (scale != 0.) {
 
153
            goto L140;
 
154
        }
 
155
L130:
 
156
        e[i__] = d__[l];
 
157
 
 
158
        i__2 = l;
 
159
        for (j = 1; j <= i__2; ++j) {
 
160
            d__[j] = z__[l + j * z_dim1];
 
161
            z__[i__ + j * z_dim1] = 0.;
 
162
            z__[j + i__ * z_dim1] = 0.;
 
163
        }
 
164
 
 
165
        goto L290;
 
166
 
 
167
L140:
 
168
        i__2 = l;
 
169
        for (k = 1; k <= i__2; ++k) {
 
170
            d__[k] /= scale;
 
171
            h__ += d__[k] * d__[k];
 
172
        }
 
173
 
 
174
        f = d__[l];
 
175
        d__1 = sqrt(h__);
 
176
        g = -d_sign(&d__1, &f);
 
177
        e[i__] = scale * g;
 
178
        h__ -= f * g;
 
179
        d__[l] = f - g;
 
180
        i__2 = l;
 
181
        for (j = 1; j <= i__2; ++j) {
 
182
            e[j] = 0.;
 
183
        }
 
184
 
 
185
        i__2 = l;
 
186
        for (j = 1; j <= i__2; ++j) {
 
187
            f = d__[j];
 
188
            z__[j + i__ * z_dim1] = f;
 
189
            g = e[j] + z__[j + j * z_dim1] * f;
 
190
            jp1 = j + 1;
 
191
            if (l < jp1) {
 
192
                goto L220;
 
193
            }
 
194
 
 
195
            i__3 = l;
 
196
            for (k = jp1; k <= i__3; ++k) {
 
197
                g += z__[k + j * z_dim1] * d__[k];
 
198
                e[k] += z__[k + j * z_dim1] * f;
 
199
            }
 
200
 
 
201
L220:
 
202
            e[j] = g;
 
203
        }
 
204
        f = 0.;
 
205
 
 
206
        i__2 = l;
 
207
        for (j = 1; j <= i__2; ++j) {
 
208
            e[j] /= h__;
 
209
            f += e[j] * d__[j];
 
210
        }
 
211
 
 
212
        hh = f / (h__ + h__);
 
213
        i__2 = l;
 
214
        for (j = 1; j <= i__2; ++j) {
 
215
            e[j] -= hh * d__[j];
 
216
        }
 
217
        i__2 = l;
 
218
        for (j = 1; j <= i__2; ++j) {
 
219
            f = d__[j];
 
220
            g = e[j];
 
221
 
 
222
            i__3 = l;
 
223
            for (k = j; k <= i__3; ++k) {
 
224
                z__[k + j * z_dim1] = z__[k + j * z_dim1] - f * e[k] - g * 
 
225
                        d__[k];
 
226
            }
 
227
 
 
228
            d__[j] = z__[l + j * z_dim1];
 
229
            z__[i__ + j * z_dim1] = 0.;
 
230
        }
 
231
 
 
232
L290:
 
233
        d__[i__] = h__;
 
234
    }
 
235
    i__1 = *n;
 
236
    for (i__ = 2; i__ <= i__1; ++i__) {
 
237
        l = i__ - 1;
 
238
        z__[*n + l * z_dim1] = z__[l + l * z_dim1];
 
239
        z__[l + l * z_dim1] = 1.;
 
240
        h__ = d__[i__];
 
241
        if (h__ == 0.) {
 
242
            goto L380;
 
243
        }
 
244
 
 
245
        i__2 = l;
 
246
        for (k = 1; k <= i__2; ++k) {
 
247
            d__[k] = z__[k + i__ * z_dim1] / h__;
 
248
        }
 
249
 
 
250
        i__2 = l;
 
251
        for (j = 1; j <= i__2; ++j) {
 
252
            g = 0.;
 
253
 
 
254
            i__3 = l;
 
255
            for (k = 1; k <= i__3; ++k) {
 
256
                g += z__[k + i__ * z_dim1] * z__[k + j * z_dim1];
 
257
            }
 
258
 
 
259
            i__3 = l;
 
260
            for (k = 1; k <= i__3; ++k) {
 
261
                z__[k + j * z_dim1] -= g * d__[k];
 
262
            }
 
263
        }
 
264
 
 
265
L380:
 
266
        i__3 = l;
 
267
        for (k = 1; k <= i__3; ++k) {
 
268
            z__[k + i__ * z_dim1] = 0.;
 
269
        }
 
270
 
 
271
    }
 
272
 
 
273
L510:
 
274
    i__1 = *n;
 
275
    for (i__ = 1; i__ <= i__1; ++i__) {
 
276
        d__[i__] = z__[*n + i__ * z_dim1];
 
277
        z__[*n + i__ * z_dim1] = 0.;
 
278
    }
 
279
 
 
280
    z__[*n + *n * z_dim1] = 1.;
 
281
    e[1] = 0.;
 
282
    return 0;
 
283
} 
 
284
 
 
285
int tql2(const integer *nm, const integer *n, doublereal *d__, 
 
286
        doublereal *e, doublereal *z__, integer *ierr)
 
287
{
 
288
    integer z_dim1, z_offset, i__1, i__2, i__3;
 
289
    doublereal d__1, d__2;
 
290
 
 
291
    double d_sign(doublereal *, doublereal *);
 
292
 
 
293
    static doublereal c__, f, g, h__;
 
294
    static integer i__, j, k, l, m;
 
295
    static doublereal p, r__, s, c2, c3;
 
296
    static integer l1, l2;
 
297
    static doublereal s2;
 
298
    static integer ii;
 
299
    static doublereal dl1, el1;
 
300
    static integer mml;
 
301
    static doublereal tst1, tst2;
 
302
    extern doublereal pythag_(doublereal *, doublereal *);
 
303
 
 
304
 
 
305
 
 
306
/*     this subroutine is a translation of the algol procedure tql2, */
 
307
/*     num. math. 11, 293-306(1968) by bowdler, martin, reinsch, and */
 
308
/*     wilkinson. */
 
309
/*     handbook for auto. comp., vol.ii-linear algebra, 227-240(1971). */
 
310
 
 
311
/*     this subroutine finds the eigenvalues and eigenvectors */
 
312
/*     of a symmetric tridiagonal matrix by the ql method. */
 
313
/*     the eigenvectors of a full symmetric matrix can also */
 
314
/*     be found if  tred2  has been used to reduce this */
 
315
/*     full matrix to tridiagonal form. */
 
316
 
 
317
/*     on input */
 
318
 
 
319
/*        nm must be set to the row dimension of two-dimensional */
 
320
/*          array parameters as declared in the calling program */
 
321
/*          dimension statement. */
 
322
 
 
323
/*        n is the order of the matrix. */
 
324
 
 
325
/*        d contains the diagonal elements of the input matrix. */
 
326
 
 
327
/*        e contains the subdiagonal elements of the input matrix */
 
328
/*          in its last n-1 positions.  e(1) is arbitrary. */
 
329
 
 
330
/*        z contains the transformation matrix produced in the */
 
331
/*          reduction by  tred2, if performed.  if the eigenvectors */
 
332
/*          of the tridiagonal matrix are desired, z must contain */
 
333
/*          the identity matrix. */
 
334
 
 
335
/*      on output */
 
336
 
 
337
/*        d contains the eigenvalues in ascending order.  if an */
 
338
/*          error exit is made, the eigenvalues are correct but */
 
339
/*          unordered for indices 1,2,...,ierr-1. */
 
340
 
 
341
/*        e has been destroyed. */
 
342
 
 
343
/*        z contains orthonormal eigenvectors of the symmetric */
 
344
/*          tridiagonal (or full) matrix.  if an error exit is made, */
 
345
/*          z contains the eigenvectors associated with the stored */
 
346
/*          eigenvalues. */
 
347
 
 
348
/*        ierr is set to */
 
349
/*          zero       for normal return, */
 
350
/*          j          if the j-th eigenvalue has not been */
 
351
/*                     determined after 30 iterations. */
 
352
 
 
353
/*     calls pythag for  dsqrt(a*a + b*b) . */
 
354
 
 
355
/*     questions and comments should be directed to burton s. garbow, */
 
356
/*     mathematics and computer science div, argonne national laboratory */
 
357
 
 
358
/*     this version dated august 1983. */
 
359
 
 
360
/*     ------------------------------------------------------------------ */
 
361
 
 
362
    z_dim1 = *nm;
 
363
    z_offset = 1 + z_dim1 * 1;
 
364
    z__ -= z_offset;
 
365
    --e;
 
366
    --d__;
 
367
 
 
368
    *ierr = 0;
 
369
    if (*n == 1) {
 
370
        goto L1001;
 
371
    }
 
372
 
 
373
    i__1 = *n;
 
374
    for (i__ = 2; i__ <= i__1; ++i__) {
 
375
        e[i__ - 1] = e[i__];
 
376
    }
 
377
 
 
378
    f = 0.;
 
379
    tst1 = 0.;
 
380
    e[*n] = 0.;
 
381
 
 
382
    i__1 = *n;
 
383
    for (l = 1; l <= i__1; ++l) {
 
384
        j = 0;
 
385
        h__ = (d__1 = d__[l], abs(d__1)) + (d__2 = e[l], abs(d__2));
 
386
        if (tst1 < h__) {
 
387
            tst1 = h__;
 
388
        }
 
389
        i__2 = *n;
 
390
        for (m = l; m <= i__2; ++m) {
 
391
            tst2 = tst1 + (d__1 = e[m], abs(d__1));
 
392
            if (tst2 == tst1) {
 
393
                goto L120;
 
394
            }
 
395
        }
 
396
 
 
397
L120:
 
398
        if (m == l) {
 
399
            goto L220;
 
400
        }
 
401
L130:
 
402
        if (j == 30) {
 
403
            goto L1000;
 
404
        }
 
405
        ++j;
 
406
        l1 = l + 1;
 
407
        l2 = l1 + 1;
 
408
        g = d__[l];
 
409
        p = (d__[l1] - g) / (e[l] * 2.);
 
410
        r__ = pythag(&p, &c_b10);
 
411
        d__[l] = e[l] / (p + d_sign(&r__, &p));
 
412
        d__[l1] = e[l] * (p + d_sign(&r__, &p));
 
413
        dl1 = d__[l1];
 
414
        h__ = g - d__[l];
 
415
        if (l2 > *n) {
 
416
            goto L145;
 
417
        }
 
418
 
 
419
        i__2 = *n;
 
420
        for (i__ = l2; i__ <= i__2; ++i__) {
 
421
            d__[i__] -= h__;
 
422
        }
 
423
 
 
424
L145:
 
425
        f += h__;
 
426
        p = d__[m];
 
427
        c__ = 1.;
 
428
        c2 = c__;
 
429
        el1 = e[l1];
 
430
        s = 0.;
 
431
        mml = m - l;
 
432
        i__2 = mml;
 
433
        for (ii = 1; ii <= i__2; ++ii) {
 
434
            c3 = c2;
 
435
            c2 = c__;
 
436
            s2 = s;
 
437
            i__ = m - ii;
 
438
            g = c__ * e[i__];
 
439
            h__ = c__ * p;
 
440
            r__ = pythag(&p, &e[i__]);
 
441
            e[i__ + 1] = s * r__;
 
442
            s = e[i__] / r__;
 
443
            c__ = p / r__;
 
444
            p = c__ * d__[i__] - s * g;
 
445
            d__[i__ + 1] = h__ + s * (c__ * g + s * d__[i__]);
 
446
            i__3 = *n;
 
447
            for (k = 1; k <= i__3; ++k) {
 
448
                h__ = z__[k + (i__ + 1) * z_dim1];
 
449
                z__[k + (i__ + 1) * z_dim1] = s * z__[k + i__ * z_dim1] + c__ 
 
450
                        * h__;
 
451
                z__[k + i__ * z_dim1] = c__ * z__[k + i__ * z_dim1] - s * h__;
 
452
            }
 
453
 
 
454
        }
 
455
 
 
456
        p = -s * s2 * c3 * el1 * e[l] / dl1;
 
457
        e[l] = s * p;
 
458
        d__[l] = c__ * p;
 
459
        tst2 = tst1 + (d__1 = e[l], abs(d__1));
 
460
        if (tst2 > tst1) {
 
461
            goto L130;
 
462
        }
 
463
L220:
 
464
        d__[l] += f;
 
465
    }
 
466
    i__1 = *n;
 
467
    for (ii = 2; ii <= i__1; ++ii) {
 
468
        i__ = ii - 1;
 
469
        k = i__;
 
470
        p = d__[i__];
 
471
 
 
472
        i__2 = *n;
 
473
        for (j = ii; j <= i__2; ++j) {
 
474
            if (d__[j] >= p) {
 
475
                goto L260;
 
476
            }
 
477
            k = j;
 
478
            p = d__[j];
 
479
L260:
 
480
            ;
 
481
        }
 
482
 
 
483
        if (k == i__) {
 
484
            goto L300;
 
485
        }
 
486
        d__[k] = d__[i__];
 
487
        d__[i__] = p;
 
488
 
 
489
        i__2 = *n;
 
490
        for (j = 1; j <= i__2; ++j) {
 
491
            p = z__[j + i__ * z_dim1];
 
492
            z__[j + i__ * z_dim1] = z__[j + k * z_dim1];
 
493
            z__[j + k * z_dim1] = p;
 
494
        }
 
495
 
 
496
L300:
 
497
        ;
 
498
    }
 
499
 
 
500
    goto L1001;
 
501
L1000:
 
502
    *ierr = l;
 
503
L1001:
 
504
    return 0;
 
505
}
 
506
 
 
507
void eigen(double **mat,octave_idx_type n,double *eig)
 
508
{
 
509
  int ierr,nm=(int)n;
 
510
  
 
511
  OCTAVE_LOCAL_BUFFER (double, trans, (octave_idx_type)nm*nm);
 
512
  OCTAVE_LOCAL_BUFFER (double, off, (octave_idx_type)nm);
 
513
 
 
514
  tred2(&nm,&nm,&mat[0][0],eig,off,trans);
 
515
  tql2(&nm,&nm,eig,off,trans,&ierr);
 
516
 
 
517
 
 
518
  if (ierr != 0)
 
519
    error ("Non converging eigenvalues!");
 
520
 
 
521
 
 
522
  for (octave_idx_type i=0;i<nm;i++)
 
523
    for (octave_idx_type j=0;j<nm;j++)
 
524
      mat[i][j]=trans[i+nm*j];
 
525
 
 
526
}