LU_unopt.c

#include <math.h>
#include "LU.h"

double LU_num_flops(int N)
{
    /* rougly 2/3*N^3 */

    double Nd = (double)N;

    return (2.0 * Nd * Nd * Nd / 3.0);
}

void LU_copy_matrix(int M, int N, double **lu, double **A)
{
    int i;
    int j;

    for (i = 0; i < M; i++)
        for (j = 0; j < N; j++)
            lu[i][j] = A[i][j];
}

int LU_factor(int M, int N, double **A, int *pivot)
{

    int minMN = M < N ? M : N;
    int j = 0;

    for (j = 0; j < minMN; j++)
    {
        /* find pivot in column j and  test for singularity. */

        int jp = j;
        int i;

        double t = fabs(A[j][j]);
        for (i = j + 1; i < M; i++)
        {
            double ab = fabs(A[i][j]);
            if (ab > t)
            {
                jp = i;
                t = ab;
            }
        }

        pivot[j] = jp;

        /* jp now has the index of maximum element  */
        /* of column j, below the diagonal          */

        if (A[jp][j] == 0)
            return 1; /* factorization failed because of zero pivot */

        if (jp != j)
        {
            /* swap rows j and jp */
            double *tA = A[j];
            A[j] = A[jp];
            A[jp] = tA;
        }

        if (j < M - 1) /* compute elements j+1:M of jth column  */
        {
            /* note A(j,j), was A(jp,p) previously which was */
            /* guarranteed not to be zero (Label #1)         */

            double recp = 1.0 / A[j][j];
            int k;
            for (k = j + 1; k < M; k++)
                A[k][j] *= recp;
        }

        if (j < minMN - 1)
        {
            /* rank-1 update to trailing submatrix:   E = E - x*y; */
            /* E is the region A(j+1:M, j+1:N) */
            /* x is the column vector A(j+1:M,j) */
            /* y is row vector A(j,j+1:N)        */

            int ii;
            for (ii = j + 1; ii < M; ii++)
            {
                double *Aii = A[ii];
                double *Aj = A[j];
                double AiiJ = Aii[j];
                int jj;
                for (jj = j + 1; jj < N; jj++)
                    Aii[jj] -= AiiJ * Aj[jj];
            }
        }
    }

    return 0;
}