← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/external/ccmath/csolv.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/csolv.c
 
1
/*  csolv.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 csolv(Cpx * a, Cpx * b, int n)
 
11
{
 
12
    int i, j, k, lc;
 
13
 
 
14
    Cpx *ps, *p, *q, *pa, *pd;
 
15
 
 
16
    Cpx z, h, *q0;
 
17
 
 
18
    double s, t, tq = 0., zr = 1.e-15;
 
19
 
 
20
    q0 = (Cpx *) calloc(n, sizeof(Cpx));
 
21
    pa = a;
 
22
    pd = a;
 
23
    for (j = 0; j < n; ++j, ++pa, pd += n + 1) {
 
24
        if (j > 0) {
 
25
            for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
 
26
                *q++ = *p;
 
27
            for (i = 1; i < n; ++i) {
 
28
                lc = i < j ? i : j;
 
29
                z.re = z.im = 0.;
 
30
                for (k = 0, p = pa + i * n - j, q = q0; k < lc; ++k, ++q, ++p) {
 
31
                    z.re += p->re * q->re - p->im * q->im;
 
32
                    z.im += p->im * q->re + p->re * q->im;
 
33
                }
 
34
                q0[i].re -= z.re;
 
35
                q0[i].im -= z.im;
 
36
            }
 
37
            for (i = 0, p = pa, q = q0; i < n; ++i, p += n)
 
38
                *p = *q++;
 
39
        }
 
40
        s = fabs(pd->re) + fabs(pd->im);
 
41
        lc = j;
 
42
        for (k = j + 1, ps = pd; k < n; ++k) {
 
43
            ps += n;
 
44
            if ((t = fabs(ps->re) + fabs(ps->im)) > s) {
 
45
                s = t;
 
46
                lc = k;
 
47
            }
 
48
        }
 
49
        tq = tq > s ? tq : s;
 
50
        if (s < zr * tq) {
 
51
            free(q0);
 
52
            return -1;
 
53
        }
 
54
        if (lc != j) {
 
55
            h = b[j];
 
56
            b[j] = b[lc];
 
57
            b[lc] = h;
 
58
            p = a + n * j;
 
59
            q = a + n * lc;
 
60
            for (k = 0; k < n; ++k) {
 
61
                h = *p;
 
62
                *p++ = *q;
 
63
                *q++ = h;
 
64
            }
 
65
        }
 
66
        t = pd->re * pd->re + pd->im * pd->im;
 
67
        h.re = pd->re / t;
 
68
        h.im = -(pd->im) / t;
 
69
        for (k = j + 1, ps = pd + n; k < n; ++k, ps += n) {
 
70
            z.re = ps->re * h.re - ps->im * h.im;
 
71
            z.im = ps->im * h.re + ps->re * h.im;
 
72
            *ps = z;
 
73
        }
 
74
    }
 
75
    for (j = 1, ps = b + 1; j < n; ++j, ++ps) {
 
76
        for (k = 0, p = a + n * j, q = b, z.re = z.im = 0.; k < j; ++k) {
 
77
            z.re += p->re * q->re - p->im * q->im;
 
78
            z.im += p->im * q->re + p->re * q->im;
 
79
            ++p;
 
80
            ++q;
 
81
        }
 
82
        ps->re -= z.re;
 
83
        ps->im -= z.im;
 
84
    }
 
85
    for (j = n - 1, --ps, pd = a + n * n - 1; j >= 0; --j, pd -= n + 1) {
 
86
        for (k = j + 1, p = pd + 1, q = b + j + 1, z.re = z.im = 0.; k < n;
 
87
             ++k) {
 
88
            z.re += p->re * q->re - p->im * q->im;
 
89
            z.im += p->im * q->re + p->re * q->im;
 
90
            ++p;
 
91
            ++q;
 
92
        }
 
93
        h.re = ps->re - z.re;
 
94
        h.im = ps->im - z.im;
 
95
        t = pd->re * pd->re + pd->im * pd->im;
 
96
        ps->re = (h.re * pd->re + h.im * pd->im) / t;
 
97
        ps->im = (h.im * pd->re - h.re * pd->im) / t;
 
98
        --ps;
 
99
    }
 
100
    free(q0);
 
101
    return 0;
 
102
}