GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
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"
10 int 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 }
void free(void *)
int sv2val(double *d, double *a, int m, int n)
Definition: sv2val.c:10
int qrbdi(double *d, double *e, int n)
Definition: qrbdi.c:9