/*
*   DIRECT FACTORIZATION ALGORITHM 6.4
*
*   To factor the n by n matrix A = (A(I,J)) into the product of the
*   lower triangular matrix L = (L(I,J)) and U = (U(I,J)), that is
*   A = LU, where the main diagonal of either L or U consists of all ones:
*
*   INPUT:   dimension n; the entries A(I,J), 1<=I, J<=n, of A;
*            the diagonal L(1,1), ..., L(N,N) of L or the diagonal
*            U(1,1), ..., U(N,N) of U.
*
*   OUTPUT:  the entries L(I,J), 1<=J<=I, 1<=I<=n of L and the entries
*            U(I,J), I<=J<=n, 1<=I<=n of U.
*/

#include<stdio.h>
#include<math.h>
#include<iostream>

#define ZERO 1.0E-20
#define true 1
#define false 0

#define N   4

static double absval(double);
static void OUTPUT(int, double [][N], int);

using namespace std;

int  main()
{
      double A[N][N] = {
                       {1.0, 1.0, 0.0, 3.0},
                       {2.0, 1.0, -1.0, 1.0},
                       {3.0, -1.0, -1.0, 2.0},
                       {-1.0, 2.0, 3.0, -1.0} 
                       };
      double  XL[N];
      double S,SS;
      int M,I,J,ISW = 0,JJ,K,KK,OK = true;
                /* ISW = 0 specifies that entris on diagonal of L are all ones */

      cout.setf(ios::fixed,ios::floatfield);
      cout.precision(8);

      for (I=1; I<=N; I++) XL[I-1] = 1.0;
      /* STEP 1 */
      if (absval(A[0][0]) <= ZERO) OK = false;
      else {
         /* the entries of L below the main diagonal will be
            placed in the corresponding entries of A; the
            entries of U above the main diagonal will be
            placed in the corresponding entries of A; the
            main diagonal which was NOT input will become
            the main diagonal of A; the input main diagonal
            of L or U is, of course, placed in XL.           */
         A[0][0] = A[0][0] / XL[0];
         /* STEP 2 */
         for (J=2; J<=N; J++) {
            if (ISW == 0) {
               /* first row of U */
               A[0][J-1] = A[0][J-1] / XL[0];
               /* first column of L */
               A[J-1][0] = A[J-1][0] / A[0][0];
            }
            else {
               /* first row of U */
               A[0][J-1] = A[0][J-1] / A[0][0];
               /* first column of L */
               A[J-1][0] = A[J-1][0] / XL[0];
            }
         }
         /* STEP 3 */
         M = N - 1;
         I = 2;
         while ((I <= M) && OK) {
            /* STEP 4 */
            KK = I - 1;
            S = 0.0;
            for (K=1; K<=KK; K++) S = S - A[I-1][K-1] * A[K-1][I-1];
            A[I-1][I-1] = ( A[I-1][I-1] + S ) / XL[I-1];
            if (absval(A[I-1][I-1]) <= ZERO) OK = false;
            else { 
               /* STEP 5 */
               JJ = I + 1;
               for (J=JJ; J<=N; J++) {
                  SS = 0.0;
                  S = 0.0;
                  for (K=1; K<=KK; K++) {
                     SS = SS - A[I-1][K-1] * A[K-1][J-1];
                     S = S - A[J-1][K-1] * A[K-1][I-1];
                  }  
                  if (ISW == 0) {
                     /* Ith row of U */
                     A[I-1][J-1] = (A[I-1][J-1] + SS) / XL[I-1];
                     /* Ith column of L */
                     A[J-1][I-1] = (A[J-1][I-1] + S) / A[I-1][I-1];
                  }
                  else {
                     /* Ith row of U */
                     A[I-1][J-1] = (A[I-1][J-1] + SS) / A[I-1][I-1];
                     /* Ith column of L */
                     A[J-1][I-1] = (A[J-1][I-1] + S) / XL[I-1];
                  }  
               }  
            }  
            I++;
         }  
         if (OK) {
            /* STEP 6 */
            S = 0.0;
            for (K=1; K<=M; K++)  S = S - A[N-1][K-1] * A[K-1][N-1];
            A[N-1][N-1] = (A[N-1][N-1] + S) / XL[N-1];
            /* If A(N-1,N-1) = 0 then A = LU but the matrix is singular.
               Process is complete, all entries of A have been
               determined. */
            /* STEP 7 */
            OUTPUT(N, A, ISW);
         }  
      }  

      if (!OK) cout<<"System has no unique solution"<<endl;
      return 0;
}

static void OUTPUT(int n, double A[][N], int ISW)
{
   int I, J, FLAG;

   cout<<"GENERAL LU FACTORIZATION"<<endl;
   if (ISW == 0) 
      cout<<"The diagonal of L consists of all entries = 1.0"<<endl;
   else
      cout<<"The diagonal of U consists of all entries = 1.0"<<endl;
   cout<<"Entries of L below diagonal and entries of U above/on diagonal"<<endl;
   cout<<"- output by rows in overwrite format:"<<endl;
   for (I=1; I<=n; I++) {
      for (J=1; J<=n; J++) cout<<"  "<< A[I-1][J-1];
      cout<<endl;
   }
}

/* Absolute Value Function */
static double absval(double val)
{
   if (val >= 0) return val;
   else return -val;
}