← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/external/ccmath/svdval.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/svdval.c
 
1
/*  svdval.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 svdval(double *d, double *a, int m, int n)
 
11
{
 
12
    double *p, *p1, *q, *w, *v;
 
13
 
 
14
    double s, h, u;
 
15
 
 
16
    int i, j, k, mm, nm, ms;
 
17
 
 
18
    if (m < n)
 
19
        return -1;
 
20
    w = (double *)calloc(m, sizeof(double));
 
21
    for (i = 0, mm = m, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
 
22
        if (mm > 1) {
 
23
            for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
 
24
                w[j] = *q;
 
25
                s += *q * *q;
 
26
            }
 
27
            if (s > 0.) {
 
28
                h = sqrt(s);
 
29
                if (*p < 0.)
 
30
                    h = -h;
 
31
                s += *p * h;
 
32
                s = 1. / s;
 
33
                w[0] += h;
 
34
                for (k = 1, ms = n - i; k < ms; ++k) {
 
35
                    for (j = 0, q = p + k, u = 0.; j < mm; q += n)
 
36
                        u += w[j++] * *q;
 
37
                    u *= s;
 
38
                    for (j = 0, q = p + k; j < mm; q += n)
 
39
                        *q -= u * w[j++];
 
40
                }
 
41
                *p = -h;
 
42
            }
 
43
        }
 
44
        p1 = p + 1;
 
45
        if (nm > 1) {
 
46
            for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
 
47
                s += *q * *q;
 
48
            if (s > 0.) {
 
49
                h = sqrt(s);
 
50
                if (*p1 < 0.)
 
51
                    h = -h;
 
52
                s += *p1 * h;
 
53
                s = 1. / s;
 
54
                *p1 += h;
 
55
                for (k = n, ms = n * (m - i); k < ms; k += n) {
 
56
                    for (j = 0, q = p1, v = p1 + k, u = 0.; j < nm; ++j)
 
57
                        u += *q++ * *v++;
 
58
                    u *= s;
 
59
                    for (j = 0, q = p1, v = p1 + k; j < nm; ++j)
 
60
                        *v++ -= u * *q++;
 
61
                }
 
62
                *p1 = -h;
 
63
            }
 
64
        }
 
65
    }
 
66
 
 
67
    for (j = 0, p = a; j < n; ++j, p += n + 1) {
 
68
        d[j] = *p;
 
69
        if (j != n - 1)
 
70
            w[j] = *(p + 1);
 
71
        else
 
72
            w[j] = 0.;
 
73
    }
 
74
    qrbdi(d, w, n);
 
75
    for (i = 0; i < n; ++i)
 
76
        if (d[i] < 0.)
 
77
            d[i] = -d[i];
 
78
    free(w);
 
79
    return 0;
 
80
}