← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to lib/external/ccmath/qrbdu1.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
lib/external/ccmath/qrbdu1.c
 
1
/*  qrbdu1.c    CCMATH mathematics library source code.
 
2
 *
 
3
 *  Copyright (C)  2000   Daniel A. Atkinson    All rights reserved.
 
4
 *  This code may be redistributed under the terms of the GNU library
 
5
 *  public license (LGPL). ( See the lgpl.license file for details.)
 
6
 * ------------------------------------------------------------------------
 
7
 */
 
8
#include "ccmath.h"
 
9
int qrbdu1(double *dm, double *em, double *um, int mm, double *vm, int m)
 
10
{
 
11
    int i, j, k, n, jj, nm;
 
12
 
 
13
    double u, x, y, a, b, c, s, t, w, *p, *q;
 
14
 
 
15
    for (j = 1, t = fabs(dm[0]); j < m; ++j)
 
16
        if ((s = fabs(dm[j]) + fabs(em[j - 1])) > t)
 
17
            t = s;
 
18
    t *= 1.e-15;
 
19
    n = 100 * m;
 
20
    nm = m;
 
21
    for (j = 0; m > 1 && j < n; ++j) {
 
22
        for (k = m - 1; k > 0; --k) {
 
23
            if (fabs(em[k - 1]) < t)
 
24
                break;
 
25
            if (fabs(dm[k - 1]) < t) {
 
26
                for (i = k, s = 1., c = 0.; i < m; ++i) {
 
27
                    a = s * em[i - 1];
 
28
                    b = dm[i];
 
29
                    em[i - 1] *= c;
 
30
                    dm[i] = u = sqrt(a * a + b * b);
 
31
                    s = -a / u;
 
32
                    c = b / u;
 
33
                    for (jj = 0, p = um + k - 1; jj < mm; ++jj, p += nm) {
 
34
                        q = p + i - k + 1;
 
35
                        w = c * *p + s * *q;
 
36
                        *q = c * *q - s * *p;
 
37
                        *p = w;
 
38
                    }
 
39
                }
 
40
                break;
 
41
            }
 
42
        }
 
43
        y = dm[k];
 
44
        x = dm[m - 1];
 
45
        u = em[m - 2];
 
46
        a = (y + x) * (y - x) - u * u;
 
47
        s = y * em[k];
 
48
        b = s + s;
 
49
        u = sqrt(a * a + b * b);
 
50
        if (u > 0.) {
 
51
            c = sqrt((u + a) / (u + u));
 
52
            if (c != 0.)
 
53
                s /= (c * u);
 
54
            else
 
55
                s = 1.;
 
56
            for (i = k; i < m - 1; ++i) {
 
57
                b = em[i];
 
58
                if (i > k) {
 
59
                    a = s * em[i];
 
60
                    b *= c;
 
61
                    em[i - 1] = u = sqrt(x * x + a * a);
 
62
                    c = x / u;
 
63
                    s = a / u;
 
64
                }
 
65
                a = c * y + s * b;
 
66
                b = c * b - s * y;
 
67
                for (jj = 0, p = vm + i; jj < nm; ++jj, p += nm) {
 
68
                    w = c * *p + s * *(p + 1);
 
69
                    *(p + 1) = c * *(p + 1) - s * *p;
 
70
                    *p = w;
 
71
                }
 
72
                s *= dm[i + 1];
 
73
                dm[i] = u = sqrt(a * a + s * s);
 
74
                y = c * dm[i + 1];
 
75
                c = a / u;
 
76
                s /= u;
 
77
                x = c * b + s * y;
 
78
                y = c * y - s * b;
 
79
                for (jj = 0, p = um + i; jj < mm; ++jj, p += nm) {
 
80
                    w = c * *p + s * *(p + 1);
 
81
                    *(p + 1) = c * *(p + 1) - s * *p;
 
82
                    *p = w;
 
83
                }
 
84
            }
 
85
        }
 
86
        em[m - 2] = x;
 
87
        dm[m - 1] = y;
 
88
        if (fabs(x) < t)
 
89
            --m;
 
90
        if (m == k + 1)
 
91
            --m;
 
92
    }
 
93
    return j;
 
94
}