1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
|
/* ldumat.c CCMATH mathematics library source code.
*
* Copyright (C) 2000 Daniel A. Atkinson All rights reserved.
* This code may be redistributed under the terms of the GNU library
* public license (LGPL). ( See the lgpl.license file for details.)
* ------------------------------------------------------------------------
*/
#include <stdlib.h>
void ldumat(double *a, double *u, int m, int n)
{
double *p0, *q0, *p, *q, *w;
int i, j, k, mm;
double s, h;
w = (double *)calloc(m, sizeof(double));
for (i = 0, mm = m * m, q = u; i < mm; ++i)
*q++ = 0.;
p0 = a + n * n - 1;
q0 = u + m * m - 1;
mm = m - n;
i = n - 1;
for (j = 0; j < mm; ++j, q0 -= m + 1)
*q0 = 1.;
if (mm == 0) {
p0 -= n + 1;
*q0 = 1.;
q0 -= m + 1;
--i;
++mm;
}
for (; i >= 0; --i, ++mm, p0 -= n + 1, q0 -= m + 1) {
if (*p0 != 0.) {
for (j = 0, p = p0 + n, h = 1.; j < mm; p += n)
w[j++] = *p;
h = *p0;
*q0 = 1. - h;
for (j = 0, q = q0 + m; j < mm; q += m)
*q = -h * w[j++];
for (k = i + 1, q = q0 + 1; k < m; ++k) {
for (j = 0, p = q + m, s = 0.; j < mm; p += m)
s += w[j++] * *p;
s *= h;
for (j = 0, p = q + m; j < mm; p += m)
*p -= s * w[j++];
*q++ = -s;
}
}
else {
*q0 = 1.;
for (j = 0, p = q0 + 1, q = q0 + m; j < mm; ++j, q += m)
*q = *p++ = 0.;
}
}
free(w);
}
|
