GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
sv2val.c
Go to the documentation of this file.
1/* sv2val.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"
10int sv2val(double *d, double *a, int m, int n)
11{
12 double *p, *p1, *q, *w, *v;
13
14 double s, h, u;
15
16 int i, j, k, mm, nm, ms;
17
18 if (m < n)
19 return -1;
20 w = (double *)calloc(m, sizeof(double));
21 for (i = 0, mm = m, p = a; i < n && mm > 1; ++i, --mm, p += n + 1) {
22 for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
23 w[j] = *q;
24 s += *q * *q;
25 }
26 if (s > 0.) {
27 h = sqrt(s);
28 if (*p < 0.)
29 h = -h;
30 s += *p * h;
31 s = 1. / s;
32 w[0] += h;
33 for (k = 1, ms = n - i; k < ms; ++k) {
34 for (j = 0, q = p + k, u = 0.; j < mm; q += n)
35 u += w[j++] * *q;
36 u = u * s;
37 for (j = 0, q = p + k; j < mm; q += n)
38 *q -= u * w[j++];
39 }
40 *p = -h;
41 }
42 }
43 for (i = 0, p = a; i < n; ++i, p += n) {
44 for (j = 0, q = p; j < i; ++j)
45 *q++ = 0.;
46 }
47 for (i = 0, mm = n, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
48 if (i && mm > 1) {
49 for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
50 w[j] = *q;
51 s += *q * *q;
52 }
53 if (s > 0.) {
54 h = sqrt(s);
55 if (*p < 0.)
56 h = -h;
57 s += *p * h;
58 s = 1. / s;
59 w[0] += h;
60 for (k = 1, ms = n - i; k < ms; ++k) {
61 for (j = 0, q = p + k, u = 0.; j < mm; q += n)
62 u += w[j++] * *q;
63 u *= s;
64 for (j = 0, q = p + k; j < mm; q += n)
65 *q -= u * w[j++];
66 }
67 *p = -h;
68 }
69 }
70 p1 = p + 1;
71 if (nm > 1) {
72 for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
73 s += *q * *q;
74 if (s > 0.) {
75 h = sqrt(s);
76 if (*p1 < 0.)
77 h = -h;
78 s += *p1 * h;
79 s = 1. / s;
80 *p1 += h;
81 for (k = n, ms = n * (m - i); k < ms; k += n) {
82 for (j = 0, q = p1, v = p1 + k, u = 0.; j < nm; ++j)
83 u += *q++ * *v++;
84 u *= s;
85 for (j = 0, q = p1, v = p1 + k; j < nm; ++j)
86 *v++ -= u * *q++;
87 }
88 *p1 = -h;
89 }
90 }
91 }
92 for (j = 0, p = a; j < n; ++j, p += n + 1) {
93 d[j] = *p;
94 if (j < n - 1)
95 w[j] = *(p + 1);
96 else
97 w[j] = 0.;
98 }
99 qrbdi(d, w, n);
100 for (i = 0; i < n; ++i)
101 if (d[i] < 0.)
102 d[i] = -d[i];
103 free(w);
104 return 0;
105}
int qrbdi(double *d, double *e, int n)
Definition qrbdi.c:9
int sv2val(double *d, double *a, int m, int n)
Definition sv2val.c:10