inverse.c

Go to the documentation of this file.

 
void ludcmp(int n, double *A)
{
    int row, col, kdx;
    int nrow;
    double val;
 
    if (n == 1) {               // 1-equation RANS model
        A[0] = 1.0/A[0];
    }
 
    else {
        for (row=1; row<n; row++) {     // Skip first row
            nrow = n*row;
            A[nrow] /= A[0];
            for (col=1; col<n; col++) {
                // Lower triangular matrix
                if (row > col) {
                    val = 0.0;
                    for (kdx=0; kdx<col; kdx++)
                        val += A[nrow+kdx] * A[n*kdx+col];
                    A[nrow+col] = (A[nrow+col] - val)/A[(n+1)*col];
                }
 
                // Upper triangular matrix
                else {
                    val = 0.0;
                    for (kdx=0; kdx<row; kdx++)
                        val += A[nrow+kdx]*A[n*kdx+col];
                    A[nrow+col] -= val;
                }
            }
        }
    }
}
 
 
void lusubst(int n, double *LU, double *b)
{
    int row, col, nrow;
    double val;
 
    if (n == 1) {                       // 1-equation RANS model
        b[0] *= LU[0];
    }
 
    else {
        // Forward substitution
        for (row=1; row<n; row++) {
            nrow = n*row;
            val = 0.0;
            for (col=0; col<row; col++)
                val += LU[nrow+col]*b[col];
            b[row] -= val;
        }
 
        // Backward substitution
        b[n-1] /= LU[n*n-1];
        for (row=n-2; row>-1; row--) {
            nrow = n*row;
            val = 0.0;
            for (col=row+1; col<n; col++)
                val += LU[nrow+col]*b[col];
            b[row] = (b[row] - val)/LU[nrow+row];
        }
    }
}