PLUTO  4.4-patch2
Functions
math_lu_decomp.c File Reference

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)
 

Detailed Description

Author
A. Mignone (migno.nosp@m.ne@t.nosp@m.o.inf.nosp@m.n.it)
Date
June 24, 2019

Function Documentation

◆ LUBackSubst()

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:

LUDecompose (A,n,indx,&d);
LUBackSubst (A,n,indx,b);

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

LUBackSubst(A,n,indx,b);

where A is the matrix that has been decomposed by LUDecompose() and not the original matrix A.

◆ LUDecompose()

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.

◆ MatrixInverse()

void MatrixInverse ( double **  A,
double **  Ainv,
int  n 
)

Find the inverse of a matrix. ! Important: the matrix A is destroyed on output !

◆ MatrixMultiply()

void MatrixMultiply ( double **  A,
double **  B,
double **  C,
int  n 
)

Multiply two matrices, C = A.B

◆ TridiagonalSolve()

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.