← Back to branch summary

~ubuntu-branches/debian/squeeze/gmsh/squeeze

« back to all changes in this revision

Viewing changes to contrib/Chaco/inertial/eigenvec2.c

  • Committer: Bazaar Package Importer
  • Author(s): Christophe Prud'homme
  • Date: 2009-01-07 16:02:08 UTC
  • mfrom: (5.1.2 sid)
  • Revision ID: james.westby@ubuntu.com-20090107160208-vklhtj69br5yw5bh
Tags: 2.2.6.dfsg-2
* debian/control: fixed lintian warning "debhelper-but-no-misc-depends"
* debian/watch: fixed lintian warning
  "debian-watch-file-should-mangle-version"

Show diffs side-by-side

added added

removed removed

Lines of Context:
contrib/Chaco/inertial/eigenvec2.c
 
1
/* This software was developed by Bruce Hendrickson and Robert Leland   *
 
2
 * at Sandia National Laboratories under US Department of Energy        *
 
3
 * contract DE-AC04-76DP00789 and is copyrighted by Sandia Corporation. */
 
4
 
 
5
#include        <math.h>
 
6
#include        "defs.h"
 
7
 
 
8
 
 
9
/* Find eigenvalues of 2x2 symmetric system by solving quadratic. */
 
10
void      evals2(H, eval1, eval2)
 
11
double    H[2][2];              /* symmetric matrix for eigenvalues */
 
12
double   *eval1;                /* smallest eigenvalue */
 
13
double   *eval2;                /* middle eigenvalue */
 
14
{
 
15
    double    M[2][2];          /* normalized version of matrix */
 
16
    double    b, c;             /* coefficents of cubic equation */
 
17
    double    root1, root2;     /* roots of quadratic */
 
18
    double    xmax;             /* largest matrix element */
 
19
    int       i, j;             /* loop counters */
 
20
 
 
21
    xmax = 0.0;
 
22
    for (i = 0; i < 2; i++) {
 
23
        for (j = i; j < 2; j++) {
 
24
            if (fabs(H[i][j]) > xmax)
 
25
                xmax = fabs(H[i][j]);
 
26
        }
 
27
    }
 
28
    if (xmax != 0) {
 
29
        for (i = 0; i < 2; i++) {
 
30
            for (j = 0; j < 2; j++)
 
31
                M[i][j] = H[i][j] / xmax;
 
32
        }
 
33
    }
 
34
 
 
35
 
 
36
    b = -M[0][0] - M[1][1];
 
37
    c = M[0][0] * M[1][1] - M[1][0] * M[1][0];
 
38
    root1 = -.5 * (b + sign(b) * sqrt(b * b - 4 * c));
 
39
    root2 = c / root1;
 
40
 
 
41
    root1 *= xmax;
 
42
    root2 *= xmax;
 
43
    *eval1 = min(root1, root2);
 
44
    *eval2 = max(root1, root2);
 
45
}
 
46
 
 
47
 
 
48
/* Solve for eigenvector of SPD 2x2 matrix, with given eigenvalue. */
 
49
void      eigenvec2(A, eval, evec, res)
 
50
double    A[2][2];              /* matrix */
 
51
double    eval;                 /* eigenvalue */
 
52
double    evec[2];              /* eigenvector returned */
 
53
double   *res;                  /* normalized residual */
 
54
{
 
55
    double    norm;             /* norm of eigenvector */
 
56
    double    res1, res2;       /* components of residual vector */
 
57
    int       i;                /* loop counter */
 
58
 
 
59
    if (fabs(A[0][0] - eval) > fabs(A[1][1] - eval)) {
 
60
        evec[0] = -A[1][0];
 
61
        evec[1] = A[0][0] - eval;
 
62
    }
 
63
    else {
 
64
        evec[0] = A[1][1] - eval;
 
65
        evec[1] = -A[1][0];
 
66
    }
 
67
 
 
68
    /* Normalize eigenvector and calculate a normalized eigen-residual. */
 
69
    norm = sqrt(evec[0] * evec[0] + evec[1] * evec[1]);
 
70
    if (norm == 0) {
 
71
        evec[0] = 1;
 
72
        evec[1] = 0;
 
73
        norm = 1;
 
74
    }
 
75
    for (i = 0; i < 2; i++)
 
76
        evec[i] /= norm;
 
77
    res1 = (A[0][0] - eval) * evec[0] + A[1][0] * evec[1];
 
78
    res2 = A[1][0] * evec[0] + (A[1][1] - eval) * evec[1];
 
79
    *res = sqrt(res1 * res1 + res2 * res2);
 
80
 
 
81
    res1 = fabs(A[0][0]) + fabs(A[1][0]);
 
82
    res2 = fabs(A[1][1]) + fabs(A[1][0]);
 
83
    *res /= max(res1, res2);
 
84
}