← Back to branch summary

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

« back to all changes in this revision

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