← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

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

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

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
}