/*
*   LDL^t ALGORITHM 6.5
*
*   To factor the positive definite n by n matrix A into LDL^t,
*   where L is a lower triangular matrix with ones along the diagonal
*   and D is a diagonal matrix with positive entries on the
*   diagonal.
*   Example 6.6.3
*
*   INPUT:   the dimension n; entries A(I,J), 1<=I, J<=n of A.
*
*   OUTPUT:  the entries L(I,J), 1<=J<I, 1<=I<=N of L and D(I),
*            1<=I<=n of D.
*/

#include<stdio.h>
#include<math.h>
#define true 1
#define false 0
#define N  3

void OUTPUT(double *, double [][N]);

main()
{
   double A[N][N] = {{4.0, -1.0, 1.0}, 
                     {-1.0, 4.25, 2.75},
                     {1.0, 2.75, 3.5}};
   double D[N],V[N];
   double S;
   int I,J,K,NN,JJ,KK,OK;

   /* STEP 1 */
   for (I=1; I<=N; I++) { 
      /* STEP 2 */
      for (J=1; J<=I-1; J++) 
          V[J-1] = A[I-1][J-1] * D[J-1];
      /* STEP 3 */
      D[I-1] = A[I-1][I-1];
      for (J=1; J<=I-1; J++) 
          D[I-1] = D[I-1] - A[I-1][J-1] * V[J-1];
      /* STEP 4 */
      for (J=I+1; J<=N; J++) {
         for (K=1; K<=I-1; K++)
            A[J-1][I-1] = A[J-1][I-1] - A[J-1][K-1] * V[K-1];
         A[J-1][I-1] = A[J-1][I-1] / D[I-1];
      }  
   }  
      /* STEP 5 */
   OUTPUT(D, A);
   return 0;
}

void OUTPUT(double *D, double A[][N])
{
   int I, J, FLAG;
   char NAME[30];
   FILE *OUP;

   printf("Choice of output method:\n");
   printf("1. Output to screen\n");
   printf("2. Output to text file\n");
   printf("Please enter 1 or 2.\n");
   scanf("%d", &FLAG);
   if (FLAG == 2) {
      printf("Input the file name in the form - drive:name.ext\n");
      printf("for example:   A:OUTPUT.DTA\n");
      scanf("%s", NAME);
      OUP = fopen(NAME, "w");
   }
   else OUP = stdout;
   fprintf(OUP, "LDL^t FACTORIZATION\n\n");
   fprintf(OUP, "The matrix L output by rows:\n");
   for (I=1; I<=N; I++) {
      for (J=1; J<=I-1; J++) fprintf(OUP, "  %12.8f", A[I-1][J-1]);
      fprintf(OUP, "\n");
   }
   fprintf(OUP, "The diagonal of D:\n");
   for (I=1; I<=N; I++) fprintf(OUP, "  %12.8f", D[I-1]);
   fprintf(OUP, "\n");
   fclose(OUP);
}