GRASS GIS 8 Programmer's Manual 8.2.0(2022)-exported
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
18void 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
61void 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
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
87void 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
121void 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
165void 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
#define NULL
Definition: ccmath.h:32
void G_free_vector(double *v)
Vector memory deallocation.
Definition: dalloc.c:129
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition: dalloc.c:60
void G_free_matrix(double **m)
Matrix memory deallocation.
Definition: dalloc.c:169
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition: debug.c:65
double b
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition: gis/error.c:160
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition: percent.c:62
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.
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_cholesky_sband_decomposition(double **A, double **T, int rows, int bandwidth)
Cholesky decomposition of a symmetric band matrix.
void G_math_solver_cholesky_sband_invert(double **A, double *x, double *b, double *invAdiag, int rows, int bandwidth)
void G_math_cholesky_sband_invert(double **A, double *invAdiag, int rows, int bandwidth)
#define x