inverse.c

Go to the documentation of this file.

 
#include "inverse.h"
 
// TODO : Bespoke code generation
 
void ludcmp(int n, UCFD_FLOAT A[n][n])
{
    int row, col, kdx;
    UCFD_FLOAT val;
 
    if (n == 1) {               // 1-equation RANS model
        A[0][0] = 1.0/A[0][0];
    }
 
    else {
        for (row=1; row<n; row++) {     // Skip first row
            A[row][0] /= A[0][0];
            for (col=1; col<n; col++) {
                // Lower triangular matrix
                if (row > col) {
                    val = 0.0;
                    for (kdx=0; kdx<col; kdx++)
                        val += A[row][kdx] * A[kdx][col];
                    A[row][col] = (A[row][col] - val)/A[col][col];
                }
 
                // Upper triangular matrix
                else {
                    val = 0.0;
                    for (kdx=0; kdx<row; kdx++)
                        val += A[row][kdx]*A[kdx][col];
                    A[row][col] -= val;
                }
            }
        }
    }
}
 
 
void lusub(int n, UCFD_FLOAT LU[n][n], UCFD_FLOAT *b)
{
    int row, col;
    UCFD_FLOAT val;
 
    if (n == 1) {                       // 1-equation RANS model
        b[0] *= LU[0][0];
    }
 
    else {
        // Forward substitution
        for (row=1; row<n; row++) {
            val = 0.0;
            for (col=0; col<row; col++)
                val += LU[row][col]*b[col];
            b[row] -= val;
        }
 
        // Backward substitution
        b[n-1] /= LU[n-1][n-1];
        for (row=n-2; row>-1; row--) {
            val = 0.0;
            for (col=row+1; col<n; col++)
                val += LU[row][col]*b[col];
            b[row] = (b[row] - val)/LU[row][row];
        }
    }
}