← Back to branch summary

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

« back to all changes in this revision

Viewing changes to lib/external/ccmath/cminv.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/cminv.c
 
1
/*  cminv.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 cminv(Cpx * a, int n)
 
11
{
 
12
    int i, j, k, m, lc, *le;
 
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
    le = (int *)calloc(n, sizeof(int));
 
21
    q0 = (Cpx *) calloc(n, sizeof(Cpx));
 
22
    pa = 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(le - j);
 
52
            free(q0);
 
53
            return -1;
 
54
        }
 
55
        *le++ = lc;
 
56
        if (lc != j) {
 
57
            p = a + n * j;
 
58
            q = a + n * lc;
 
59
            for (k = 0; k < n; ++k, ++p, ++q) {
 
60
                h = *p;
 
61
                *p = *q;
 
62
                *q = h;
 
63
            }
 
64
        }
 
65
        t = pd->re * pd->re + pd->im * pd->im;
 
66
        h.re = pd->re / t;
 
67
        h.im = -(pd->im) / t;
 
68
        for (k = j + 1, ps = pd + n; k < n; ++k, ps += n) {
 
69
            z.re = ps->re * h.re - ps->im * h.im;
 
70
            z.im = ps->im * h.re + ps->re * h.im;
 
71
            *ps = z;
 
72
        }
 
73
        *pd = h;
 
74
    }
 
75
    for (j = 1, pd = ps = a; j < n; ++j) {
 
76
        for (k = 0, pd += n + 1, q = ++ps; k < j; ++k, q += n) {
 
77
            z.re = q->re * pd->re - q->im * pd->im;
 
78
            z.im = q->im * pd->re + q->re * pd->im;
 
79
            *q = z;
 
80
        }
 
81
    }
 
82
    for (j = 1, pa = a; j < n; ++j) {
 
83
        ++pa;
 
84
        for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
 
85
            *q++ = *p;
 
86
        for (k = 0; k < j; ++k) {
 
87
            h.re = h.im = 0.;
 
88
            for (i = k, p = pa + k * n + k - j, q = q0 + k; i < j; ++i) {
 
89
                h.re -= p->re * q->re - p->im * q->im;
 
90
                h.im -= p->im * q->re + p->re * q->im;
 
91
                ++p;
 
92
                ++q;
 
93
            }
 
94
            q0[k] = h;
 
95
        }
 
96
        for (i = 0, q = q0, p = pa; i < j; ++i, p += n)
 
97
            *p = *q++;
 
98
    }
 
99
    for (j = n - 2, pd = pa = a + n * n - 1; j >= 0; --j) {
 
100
        --pa;
 
101
        pd -= n + 1;
 
102
        for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
 
103
            *q++ = *p;
 
104
        for (k = n - 1, ps = pa; k > j; --k, ps -= n) {
 
105
            z.re = -ps->re;
 
106
            z.im = -ps->im;
 
107
            for (i = j + 1, p = ps + 1, q = q0; i < k; ++i, ++p, ++q) {
 
108
                z.re -= p->re * q->re - p->im * q->im;
 
109
                z.im -= p->im * q->re + p->re * q->im;
 
110
            }
 
111
            q0[--m] = z;
 
112
        }
 
113
        for (i = 0, m = n - j - 1, q = q0, p = pd + n; i < m; ++i, p += n)
 
114
            *p = *q++;
 
115
    }
 
116
    for (k = 0, pa = a; k < n - 1; ++k, ++pa) {
 
117
        for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
 
118
            *q++ = *p;
 
119
        for (j = 0, ps = a; j < n; ++j, ps += n) {
 
120
            if (j > k) {
 
121
                h.re = h.im = 0.;
 
122
                p = ps + j;
 
123
                i = j;
 
124
            }
 
125
            else {
 
126
                h = q0[j];
 
127
                p = ps + k + 1;
 
128
                i = k + 1;
 
129
            }
 
130
            for (; i < n; ++i, ++p) {
 
131
                h.re += p->re * q0[i].re - p->im * q0[i].im;
 
132
                h.im += p->im * q0[i].re + p->re * q0[i].im;
 
133
            }
 
134
            q0[j] = h;
 
135
        }
 
136
        for (i = 0, q = q0, p = pa; i < n; ++i, p += n)
 
137
            *p = *q++;
 
138
    }
 
139
    for (j = n - 2, le--; j >= 0; --j) {
 
140
        for (k = 0, p = a + j, q = a + *(--le); k < n; ++k, p += n, q += n) {
 
141
            h = *p;
 
142
            *p = *q;
 
143
            *q = h;
 
144
        }
 
145
    }
 
146
    free(le);
 
147
    free(q0);
 
148
    return 0;
 
149
}