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