GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
svdval.c
Go to the documentation of this file.
1/* svdval.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 svdval(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, nm = n - 1, p = a; i < n; ++i, --mm, --nm, p += n + 1) {
22 if (mm > 1) {
23 for (j = 0, q = p, s = 0.; j < mm; ++j, q += n) {
24 w[j] = *q;
25 s += *q * *q;
26 }
27 if (s > 0.) {
28 h = sqrt(s);
29 if (*p < 0.)
30 h = -h;
31 s += *p * h;
32 s = 1. / s;
33 w[0] += h;
34 for (k = 1, ms = n - i; k < ms; ++k) {
35 for (j = 0, q = p + k, u = 0.; j < mm; q += n)
36 u += w[j++] * *q;
37 u *= s;
38 for (j = 0, q = p + k; j < mm; q += n)
39 *q -= u * w[j++];
40 }
41 *p = -h;
42 }
43 }
44 p1 = p + 1;
45 if (nm > 1) {
46 for (j = 0, q = p1, s = 0.; j < nm; ++j, ++q)
47 s += *q * *q;
48 if (s > 0.) {
49 h = sqrt(s);
50 if (*p1 < 0.)
51 h = -h;
52 s += *p1 * h;
53 s = 1. / s;
54 *p1 += h;
55 for (k = n, ms = n * (m - i); k < ms; k += n) {
56 for (j = 0, q = p1, v = p1 + k, u = 0.; j < nm; ++j)
57 u += *q++ * *v++;
58 u *= s;
59 for (j = 0, q = p1, v = p1 + k; j < nm; ++j)
60 *v++ -= u * *q++;
61 }
62 *p1 = -h;
63 }
64 }
65 }
66
67 for (j = 0, p = a; j < n; ++j, p += n + 1) {
68 d[j] = *p;
69 if (j != n - 1)
70 w[j] = *(p + 1);
71 else
72 w[j] = 0.;
73 }
74 qrbdi(d, w, n);
75 for (i = 0; i < n; ++i)
76 if (d[i] < 0.)
77 d[i] = -d[i];
78 free(w);
79 return 0;
80}
int qrbdi(double *d, double *e, int n)
Definition qrbdi.c:9
int svdval(double *d, double *a, int m, int n)
Definition svdval.c:10