Loading...
Searching...
No Matches
solvers_direct_cholesky_band.c
Go to the documentation of this file.
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/gmath.h>
#include <grass/glocale.h>
/**
* \brief Cholesky decomposition of a symmetric band matrix
*
* \param A (double**) the input symmetric band matrix
* \param T (double**) the resulting lower tringular symmetric band matrix
* \param rows (int) number of rows
* \param bandwidth (int) the bandwidth of the symmetric band matrix
*
* */
void G_math_cholesky_sband_decomposition(double **A, double **T, int rows,
int bandwidth)
{
int i, j, k, end;
double sum;
G_debug(2, "G_math_cholesky_sband_decomposition(): n=%d bandwidth=%d",
rows, bandwidth);
for (i = 0; i < rows; i++) {
G_percent(i, rows, 9);
/* For j = 0 */
sum = A[i][0];
end = ((bandwidth - 0) < (i + 1) ? (bandwidth - 0) : (i + 1));
for (k = 1; k < end; k++)
sum -= T[i - k][k] * T[i - k][0 + k];
if (sum <= 0.0)
G_fatal_error(_("Decomposition failed at row %i and col %i"), i, 0);
T[i][0] = sqrt(sum);
#pragma omp parallel for schedule(static) private(j, k, end, sum) \
shared(A, T, i, bandwidth)
for (j = 1; j < bandwidth; j++) {
sum = A[i][j];
end = ((bandwidth - j) < (i + 1) ? (bandwidth - j) : (i + 1));
for (k = 1; k < end; k++)
sum -= T[i - k][k] * T[i - k][j + k];
T[i][j] = sum / T[i][0];
}
}
G_percent(i, rows, 2);
return;
}
/**
* \brief Cholesky symmetric band matrix solver for linear equation systems of
* type Ax = b
*
* \param A (double**) the input symmetric band matrix
* \param x (double*) the resulting vector, result is written in here
* \param b (double*) the right hand side of Ax = b
* \param rows (int) number of rows
* \param bandwidth (int) the bandwidth of the symmetric band matrix
*
* */
void G_math_solver_cholesky_sband(double **A, double *x, double *b, int rows,
int bandwidth)
{
double **T;
T = G_alloc_matrix(rows, bandwidth);
G_math_cholesky_sband_decomposition(A, T, rows,
bandwidth); /* T computation */
G_math_cholesky_sband_substitution(T, x, b, rows, bandwidth);
G_free_matrix(T);
return;
}
/**
* \brief Forward and backward substitution of a lower tringular symmetric band
* matrix of A from system Ax = b
*
* \param T (double**) the lower triangle symmetric band matrix
* \param x (double*) the resulting vector
* \param b (double*) the right hand side of Ax = b
* \param rows (int) number of rows
* \param bandwidth (int) the bandwidth of the symmetric band matrix
*
* */
void G_math_cholesky_sband_substitution(double **T, double *x, double *b,
int rows, int bandwidth)
{
int i, j, start, end;
/* Forward substitution */
x[0] = b[0] / T[0][0];
for (i = 1; i < rows; i++) {
x[i] = b[i];
/* start = 0 or i - bandwidth + 1 */
start = ((i - bandwidth + 1) < 0 ? 0 : (i - bandwidth + 1));
/* end = i */
for (j = start; j < i; j++)
x[i] -= T[j][i - j] * x[j];
x[i] = x[i] / T[i][0];
}
/* Backward substitution */
x[rows - 1] = x[rows - 1] / T[rows - 1][0];
for (i = rows - 2; i >= 0; i--) {
/* start = i + 1 */
/* end = rows or bandwidth + i */
end = (rows < (bandwidth + i) ? rows : (bandwidth + i));
for (j = i + 1; j < end; j++)
x[i] -= T[i][j - i] * x[j];
x[i] = x[i] / T[i][0];
}
return;
}
/*--------------------------------------------------------------------------------------*/
/* Tcholetsky matrix inversion */
void G_math_cholesky_sband_invert(double **A, double *invAdiag, int rows,
int bandwidth)
{
double **T = NULL;
double *vect = NULL;
int i, j, k, start;
double sum;
T = G_alloc_matrix(rows, bandwidth);
vect = G_alloc_vector(rows);
/* T computation */
G_math_cholesky_sband_decomposition(A, T, rows, bandwidth);
/* T Diagonal inversion */
for (i = 0; i < rows; i++) {
T[i][0] = 1.0 / T[i][0];
}
/* A Diagonal inversion */
for (i = 0; i < rows; i++) {
vect[0] = T[i][0];
invAdiag[i] = vect[0] * vect[0];
for (j = i + 1; j < rows; j++) {
sum = 0.0;
/* start = i or j - bandwidth + 1 */
start = ((j - bandwidth + 1) < i ? i : (j - bandwidth + 1));
/* end = j */
for (k = start; k < j; k++) {
sum -= vect[k - i] * T[k][j - k];
}
vect[j - i] = sum * T[j][0];
invAdiag[i] += vect[j - i] * vect[j - i];
}
}
G_free_matrix(T);
G_free_vector(vect);
return;
}
/*--------------------------------------------------------------------------------------*/
/* Tcholetsky matrix solution and inversion */
void G_math_solver_cholesky_sband_invert(double **A, double *x, double *b,
double *invAdiag, int rows,
int bandwidth)
{
double **T = NULL;
double *vect = NULL;
int i, j, k, start;
double sum;
T = G_alloc_matrix(rows, bandwidth);
vect = G_alloc_vector(rows);
/* T computation */
G_math_cholesky_sband_decomposition(A, T, rows, bandwidth);
G_math_cholesky_sband_substitution(T, x, b, rows, bandwidth);
/* T Diagonal inversion */
for (i = 0; i < rows; i++) {
T[i][0] = 1.0 / T[i][0];
}
/* A Diagonal inversion */
for (i = 0; i < rows; i++) {
vect[0] = T[i][0];
invAdiag[i] = vect[0] * vect[0];
for (j = i + 1; j < rows; j++) {
sum = 0.0;
/* start = i or j - bandwidth + 1 */
start = ((j - bandwidth + 1) < i ? i : (j - bandwidth + 1));
/* end = j */
for (k = start; k < j; k++) {
sum -= vect[k - i] * T[k][j - k];
}
vect[j - i] = sum * T[j][0];
invAdiag[i] += vect[j - i] * vect[j - i];
}
}
G_free_matrix(T);
G_free_vector(vect);
return;
}
NULL
#define NULL
Definition ccmath.h:32
G_free_vector
void G_free_vector(double *v)
Vector memory deallocation.
Definition dalloc.c:118
G_alloc_matrix
double ** G_alloc_matrix(int rows, int cols)
Matrix memory allocation.
Definition dalloc.c:55
G_free_matrix
void G_free_matrix(double **m)
Matrix memory deallocation.
Definition dalloc.c:154
G_alloc_vector
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:38
G_debug
int G_debug(int level, const char *msg,...)
Print debugging message.
Definition debug.c:66
b
double b
Definition driver/set_window.c:5
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
G_percent
void G_percent(long n, long d, int s)
Print percent complete messages.
Definition percent.c:61
G_math_solver_cholesky_sband
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.
Definition solvers_direct_cholesky_band.c:65
G_math_cholesky_sband_substitution
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.
Definition solvers_direct_cholesky_band.c:94
G_math_cholesky_sband_decomposition
void G_math_cholesky_sband_decomposition(double **A, double **T, int rows, int bandwidth)
Cholesky decomposition of a symmetric band matrix.
Definition solvers_direct_cholesky_band.c:18
G_math_solver_cholesky_sband_invert
void G_math_solver_cholesky_sband_invert(double **A, double *x, double *b, double *invAdiag, int rows, int bandwidth)
Definition solvers_direct_cholesky_band.c:174
G_math_cholesky_sband_invert
void G_math_cholesky_sband_invert(double **A, double *invAdiag, int rows, int bandwidth)
Definition solvers_direct_cholesky_band.c:129
x
#define x