← Back to branch summary

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

« back to all changes in this revision

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