PLUTO
4.4-patch2
|
Functions for LU decomposition and matrix inversion. More...
#include "pluto.h"
Functions | |
int | LUDecompose (double **a, int n, int *indx, double *d) |
void | LUBackSubst (double **a, int n, int *indx, double b[]) |
void | MatrixInverse (double **A, double **Ainv, int n) |
void | MatrixMultiply (double **A, double **B, double **C, int n) |
void | TridiagonalSolve (double *am, double *a0, double *ap, double *b, double *y, int n) |
void LUBackSubst | ( | double ** | a, |
int | n, | ||
int * | indx, | ||
double | b[] | ||
) |
Solve a linear system of the type Ax = b
, where A
is a square matrix, x
is the array of unknowns and b
is given. This function must be called after LUDecompose() (adapted from Numerical Recipes). For the 1st time, use:
The answer x
will be given back in b
and the original matrix A
will be destroyed. For all subsequent time, to solve the same system wth a different right-hand side b
, use
where A
is the matrix that has been decomposed by LUDecompose() and not the original matrix A
.
int LUDecompose | ( | double ** | a, |
int | n, | ||
int * | indx, | ||
double * | d | ||
) |
Perform LU decomposition. Adapted from Numerical Recipes. On output, replaces the matrix a[0..n-1][0..n-1] with its LU decomposition. indx[0..n-1] is a row vector that records the row permutation while d = 1 (-1) if the number of raw interchanges was even (odd). This function should be used in combination with LUBackSubst() to solve linear system of equations.
void MatrixInverse | ( | double ** | A, |
double ** | Ainv, | ||
int | n | ||
) |
Find the inverse of a matrix. ! Important: the matrix A is destroyed on output !
void MatrixMultiply | ( | double ** | A, |
double ** | B, | ||
double ** | C, | ||
int | n | ||
) |
Multiply two matrices, C = A.B
void TridiagonalSolve | ( | double * | am, |
double * | a0, | ||
double * | ap, | ||
double * | b, | ||
double * | y, | ||
int | n | ||
) |
Use recursion formula to solve a tridiagonal system of the type
am[i]*y[i-1] + a0[i]*y[i] + ap[i]*y[i+1] = b[i]
from i = 1..n-1. Boundary coefficients must have been imposed at y[0] and y[n] On output y[] is replaced with the updated solution.