GRASS GIS 8 Programmer's Manual  8.2.2dev(2023)-3d2c704037
solvers_direct_cholesky_band.c
Go to the documentation of this file.
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <math.h>
4 #include <grass/gis.h>
5 #include <grass/gmath.h>
6 #include <grass/glocale.h>
7 
8 /**
9  * \brief Cholesky decomposition of a symmetric band matrix
10  *
11  * \param A (double**) the input symmetric band matrix
12  * \param T (double**) the resulting lower tringular symmetric band matrix
13  * \param rows (int) number of rows
14  * \param bandwidth (int) the bandwidth of the symmetric band matrix
15  *
16  * */
17 
18 void G_math_cholesky_sband_decomposition(double **A, double **T, int rows, int bandwidth)
19 {
20  int i, j, k, end;
21  double sum;
22 
23  G_debug(2, "G_math_cholesky_sband_decomposition(): n=%d bandwidth=%d", rows, bandwidth);
24 
25  for (i = 0; i < rows; i++) {
26  G_percent(i, rows, 9);
27  /* For j = 0 */
28  sum = A[i][0];
29  end = ((bandwidth - 0) < (i + 1) ? (bandwidth - 0) : (i + 1));
30  for (k = 1; k < end; k++)
31  sum -= T[i - k][k] * T[i - k][0 + k];
32  if (sum <= 0.0)
33  G_fatal_error(_("Decomposition failed at row %i and col %i"), i, 0);
34  T[i][0] = sqrt(sum);
35 
36  #pragma omp parallel for schedule (static) private(j, k, end, sum) shared(A, T, i, bandwidth)
37  for (j = 1; j < bandwidth; j++) {
38  sum = A[i][j];
39  end = ((bandwidth - j) < (i + 1) ? (bandwidth - j) : (i + 1));
40  for (k = 1; k < end; k++)
41  sum -= T[i - k][k] * T[i - k][j + k];
42  T[i][j] = sum / T[i][0];
43  }
44  }
45 
46  G_percent(i, rows, 2);
47  return;
48 }
49 
50 /**
51  * \brief Cholesky symmetric band matrix solver for linear equation systems of type Ax = b
52  *
53  * \param A (double**) the input symmetric band matrix
54  * \param x (double*) the resulting vector, result is written in here
55  * \param b (double*) the right hand side of Ax = b
56  * \param rows (int) number of rows
57  * \param bandwidth (int) the bandwidth of the symmetric band matrix
58  *
59  * */
60 
61 void G_math_solver_cholesky_sband(double **A, double *x, double *b, int rows, int bandwidth)
62 {
63 
64  double **T;
65 
66  T = G_alloc_matrix(rows, bandwidth);
67 
68  G_math_cholesky_sband_decomposition(A, T, rows, bandwidth); /* T computation */
69  G_math_cholesky_sband_substitution(T, x, b, rows, bandwidth);
70 
71  G_free_matrix(T);
72 
73  return;
74 }
75 
76 /**
77  * \brief Forward and backward substitution of a lower tringular symmetric band matrix of A from system Ax = b
78  *
79  * \param T (double**) the lower triangle symmetric band matrix
80  * \param x (double*) the resulting vector
81  * \param b (double*) the right hand side of Ax = b
82  * \param rows (int) number of rows
83  * \param bandwidth (int) the bandwidth of the symmetric band matrix
84  *
85  * */
86 
87 void G_math_cholesky_sband_substitution(double **T, double *x, double *b, int rows, int bandwidth)
88 {
89 
90  int i, j, start, end;
91 
92  /* Forward substitution */
93  x[0] = b[0] / T[0][0];
94  for (i = 1; i < rows; i++) {
95  x[i] = b[i];
96  /* start = 0 or i - bandwidth + 1 */
97  start = ((i - bandwidth + 1) < 0 ? 0 : (i - bandwidth + 1));
98  /* end = i */
99  for (j = start; j < i; j++)
100  x[i] -= T[j][i - j] * x[j];
101  x[i] = x[i] / T[i][0];
102  }
103 
104  /* Backward substitution */
105  x[rows - 1] = x[rows - 1] / T[rows - 1][0];
106  for (i = rows - 2; i >= 0; i--) {
107  /* start = i + 1 */
108  /* end = rows or bandwidth + i */
109  end = (rows < (bandwidth + i) ? rows : (bandwidth + i));
110  for (j = i + 1; j < end; j++)
111  x[i] -= T[i][j - i] * x[j];
112  x[i] = x[i] / T[i][0];
113  }
114 
115  return;
116 }
117 
118 /*--------------------------------------------------------------------------------------*/
119 /* Tcholetsky matrix invertion */
120 
121 void G_math_cholesky_sband_invert(double **A, double *invAdiag, int rows, int bandwidth)
122 {
123  double **T = NULL;
124  double *vect = NULL;
125  int i, j, k, start;
126  double sum;
127 
128  T = G_alloc_matrix(rows, bandwidth);
129  vect = G_alloc_vector(rows);
130 
131  /* T computation */
132  G_math_cholesky_sband_decomposition(A, T, rows, bandwidth);
133 
134  /* T Diagonal invertion */
135  for (i = 0; i < rows; i++) {
136  T[i][0] = 1.0 / T[i][0];
137  }
138 
139  /* A Diagonal invertion */
140  for (i = 0; i < rows; i++) {
141  vect[0] = T[i][0];
142  invAdiag[i] = vect[0] * vect[0];
143  for (j = i + 1; j < rows; j++) {
144  sum = 0.0;
145  /* start = i or j - bandwidth + 1 */
146  start = ((j - bandwidth + 1) < i ? i : (j - bandwidth + 1));
147  /* end = j */
148  for (k = start; k < j; k++) {
149  sum -= vect[k - i] * T[k][j - k];
150  }
151  vect[j - i] = sum * T[j][0];
152  invAdiag[i] += vect[j - i] * vect[j - i];
153  }
154  }
155 
156  G_free_matrix(T);
157  G_free_vector(vect);
158 
159  return;
160 }
161 
162 /*--------------------------------------------------------------------------------------*/
163 /* Tcholetsky matrix solution and invertion */
164 
165 void G_math_solver_cholesky_sband_invert(double **A, double *x, double *b, double *invAdiag,
166  int rows, int bandwidth)
167 {
168 
169  double **T = NULL;
170  double *vect = NULL;
171  int i, j, k, start;
172  double sum;
173 
174  T = G_alloc_matrix(rows, bandwidth);
175  vect = G_alloc_vector(rows);
176 
177  /* T computation */
178  G_math_cholesky_sband_decomposition(A, T, rows, bandwidth);
179  G_math_cholesky_sband_substitution(T, x, b, rows, bandwidth);
180 
181  /* T Diagonal invertion */
182  for (i = 0; i < rows; i++) {
183  T[i][0] = 1.0 / T[i][0];
184  }
185 
186  /* A Diagonal invertion */
187  for (i = 0; i < rows; i++) {
188  vect[0] = T[i][0];
189  invAdiag[i] = vect[0] * vect[0];
190  for (j = i + 1; j < rows; j++) {
191  sum = 0.0;
192  /* start = i or j - bandwidth + 1 */
193  start = ((j - bandwidth + 1) < i ? i : (j - bandwidth + 1));
194  /* end = j */
195  for (k = start; k < j; k++) {
196  sum -= vect[k - i] * T[k][j - k];
197  }
198  vect[j - i] = sum * T[j][0];
199  invAdiag[i] += vect[j - i] * vect[j - i];
200  }
201  }
202 
203  G_free_matrix(T);
204  G_free_vector(vect);
205 
206  return;
207 }
208 
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
void G_free_vector(double *)
Vector memory deallocation.
Definition: dalloc.c:129
void G_math_cholesky_sband_invert(double **A, double *invAdiag, int rows, int bandwidth)
double ** G_alloc_matrix(int, int)
Matrix memory allocation.
Definition: dalloc.c:60
#define NULL
Definition: ccmath.h:32
#define x
double b
Definition: r_raster.c:39
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
Definition: dalloc.c:41
void G_math_cholesky_sband_decomposition(double **A, double **T, int rows, int bandwidth)
Cholesky decomposition of a symmetric band matrix.
void G_math_cholesky_sband_substitution(double **T, double *x, double *b, int rows, int bandwidth)
Forward and backward substitution of a lower tringular symmetric band matrix of A from system Ax = b...
void G_math_solver_cholesky_sband_invert(double **A, double *x, double *b, double *invAdiag, int rows, int bandwidth)
void G_percent(long, long, int)
Print percent complete messages.
Definition: percent.c:62
#define _(str)
Definition: glocale.h:10
void G_free_matrix(double **)
Matrix memory deallocation.
Definition: dalloc.c:169
void G_math_solver_cholesky_sband(double **A, double *x, double *b, int rows, int bandwidth)
Cholesky symmetric band matrix solver for linear equation systems of type Ax = b. ...
int G_debug(int, const char *,...) __attribute__((format(printf