GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-6596d49e63
housev.c
Go to the documentation of this file.
1 /* housev.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 void housev(double *a, double *d, double *dp, int n)
11 {
12  double sc, x, y, h;
13 
14  int i, j, k, m, e;
15 
16  double *qw, *qs, *pc, *p;
17 
18  qs = (double *)calloc(n, sizeof(double));
19  for (j = 0, pc = a; j < n - 2; ++j, pc += n + 1) {
20  m = n - j - 1;
21  for (i = 1, sc = 0.; i <= m; ++i)
22  sc += pc[i] * pc[i];
23  if (sc > 0.) {
24  sc = sqrt(sc);
25  if ((x = *(pc + 1)) < 0.) {
26  y = x - sc;
27  h = 1. / sqrt(-2. * sc * y);
28  }
29  else {
30  y = x + sc;
31  h = 1. / sqrt(2. * sc * y);
32  sc = -sc;
33  }
34  for (i = 0, qw = pc + 1; i < m; ++i) {
35  qs[i] = 0.;
36  if (i)
37  qw[i] *= h;
38  else
39  qw[i] = y * h;
40  }
41  for (i = 0, e = j + 2, p = pc + n + 1, h = 0.; i < m;
42  ++i, p += e++) {
43  qs[i] += (y = qw[i]) * *p++;
44  for (k = i + 1; k < m; ++k) {
45  qs[i] += qw[k] * *p;
46  qs[k] += y * *p++;
47  }
48  h += y * qs[i];
49  }
50  for (i = 0; i < m; ++i) {
51  qs[i] -= h * qw[i];
52  qs[i] += qs[i];
53  }
54  for (i = 0, e = j + 2, p = pc + n + 1; i < m; ++i, p += e++) {
55  for (k = i; k < m; ++k)
56  *p++ -= qw[i] * qs[k] + qs[i] * qw[k];
57  }
58  }
59  d[j] = *pc;
60  dp[j] = sc;
61  }
62  d[j] = *pc;
63  dp[j] = *(pc + 1);
64  d[j + 1] = *(pc += n + 1);
65  free(qs);
66  for (i = 0, m = n + n, p = pc; i < m; ++i)
67  *p-- = 0.;
68  *pc = 1.;
69  *(pc -= n + 1) = 1.;
70  qw = pc - n;
71  for (m = 2; m < n; ++m, qw -= n + 1) {
72  for (j = 0, p = pc, *pc = 1.; j < m; ++j, p += n) {
73  for (i = 0, qs = p, h = 0.; i < m;)
74  h += qw[i++] * *qs++;
75  for (i = 0, qs = p, h += h; i < m;)
76  *qs++ -= h * qw[i++];
77  }
78  for (i = 0, p = qw + m; i < n; ++i)
79  *(--p) = 0.;
80  *(pc -= n + 1) = 1.;
81  }
82 }
void housev(double *a, double *d, double *dp, int n)
Definition: housev.c:10
void free(void *)
#define x