← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/external/ccmath/svduv.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/svduv.c
 
1
/*  svduv.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 <stdlib.h>
 
9
#include "ccmath.h"
 
10
int svduv(double *d, double *a, double *u, int m, double *v, int n)
 
11
{
 
12
    double *p, *p1, *q, *pp, *w, *e;
 
13
 
 
14
    double s, h, r, t, sv;
 
15
 
 
16
    int i, j, k, mm, nm, ms;
 
17
 
 
18
    if (m < n)
 
19
        return -1;
 
20
    w = (double *)calloc(m + n, sizeof(double));
 
21
    e = w + m;
 
22
    for (i = 0, mm = m, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
 
23
        if (mm > 1) {
 
24
            sv = h = 0.;
 
25
            for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
 
26
                w[j] = *q;
 
27
                s += *q * *q;
 
28
            }
 
29
            if (s > 0.) {
 
30
                h = sqrt(s);
 
31
                if (*p < 0.)
 
32
                    h = -h;
 
33
                s += *p * h;
 
34
                s = 1. / s;
 
35
                t = 1. / (w[0] += h);
 
36
                sv = 1. + fabs(*p / h);
 
37
                for (k = 1, ms = n - i; k < ms; ++k) {
 
38
                    for (j = 0, q = p + k, r = 0.; j < mm; q += n)
 
39
                        r += w[j++] * *q;
 
40
                    r *= s;
 
41
                    for (j = 0, q = p + k; j < mm; q += n)
 
42
                        *q -= r * w[j++];
 
43
                }
 
44
                for (j = 1, q = p; j < mm;)
 
45
                    *(q += n) = t * w[j++];
 
46
            }
 
47
            *p = sv;
 
48
            d[i] = -h;
 
49
        }
 
50
        if (mm == 1)
 
51
            d[i] = *p;
 
52
        p1 = p + 1;
 
53
        sv = h = 0.;
 
54
        if (nm > 1) {
 
55
            for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
 
56
                s += *q * *q;
 
57
            if (s > 0.) {
 
58
                h = sqrt(s);
 
59
                if (*p1 < 0.)
 
60
                    h = -h;
 
61
                sv = 1. + fabs(*p1 / h);
 
62
                s += *p1 * h;
 
63
                s = 1. / s;
 
64
                t = 1. / (*p1 += h);
 
65
                for (k = n, ms = n * (m - i); k < ms; k += n) {
 
66
                    for (j = 0, q = p1, pp = p1 + k, r = 0.; j < nm; ++j)
 
67
                        r += *q++ * *pp++;
 
68
                    r *= s;
 
69
                    for (j = 0, q = p1, pp = p1 + k; j < nm; ++j)
 
70
                        *pp++ -= r * *q++;
 
71
                }
 
72
                for (j = 1, q = p1 + 1; j < nm; ++j)
 
73
                    *q++ *= t;
 
74
            }
 
75
            *p1 = sv;
 
76
            e[i] = -h;
 
77
        }
 
78
        if (nm == 1)
 
79
            e[i] = *p1;
 
80
    }
 
81
    ldvmat(a, v, n);
 
82
    ldumat(a, u, m, n);
 
83
    qrbdv(d, e, u, m, v, n);
 
84
    for (i = 0; i < n; ++i) {
 
85
        if (d[i] < 0.) {
 
86
            d[i] = -d[i];
 
87
            for (j = 0, p = v + i; j < n; ++j, p += n)
 
88
                *p = -*p;
 
89
        }
 
90
    }
 
91
    free(w);
 
92
    return 0;
 
93
}