/*************************************************************
*       SOLVE A LINEAR SYSTEM BY DIRECT FACTORIZATION        *
* ---------------------------------------------------------- *
* SAMPLE RUN:                                                *
* (Solve linear system A X = B, where:                       *
*                                                            *
*       1  0  0  0  0  1              1                      *
*       1  1  0  0  0 -1              0                      *
*  A = -1  1  1  0  0  1          B = 1                      *
*       1 -1  1  1  0 -1              0                      *
*      -1  1 -1  1  1  1              1                      *
*       1 -1  1 -1  1 -1              0 )                    *
*                                                            *
* LINEAR SYSTEM AX = B:                                      *
*                                                            *
*  1.0000  0.0000  0.0000  0.0000  0.0000  1.0000    1.0000  *
*  1.0000  1.0000  0.0000  0.0000  0.0000 -1.0000    0.0000  *
* -1.0000  1.0000  1.0000  0.0000  0.0000  1.0000    1.0000  *
*  1.0000 -1.0000  1.0000  1.0000  0.0000 -1.0000    0.0000  *
* -1.0000  1.0000 -1.0000  1.0000  1.0000  1.0000    1.0000  *
*  1.0000 -1.0000  1.0000 -1.0000  1.0000 -1.0000    0.0000  *
*                                                            *
* SOLUTION:                                                  *
*  0.343750                                                  *
*  0.312500                                                  *
*  0.375000                                                  *
*  0.250000                                                  *
*  0.500000                                                  *
*  0.656250                                                  *
*                                                            *
* ---------------------------------------------------------- *
* Ref.: From Numath Library By Tuan Dang Trong in fortran 77 *
*       [BIBLI 18].                                          * 
*                                                            *
*                           C++ Release By J-P Moreau, Paris *
*                                 (www.jpmoreau.fr)          *
*************************************************************/
#include <stdio.h>
#include <math.h>

#define  NMAX  50

typedef  double MAT[NMAX][NMAX];

    MAT A, W;
    double B[NMAX], X[NMAX], Z[NMAX];
    int I,J,N;


    void DLITTL(MAT A, double *B, int N, double *X, MAT W, double *Z)  {
/*----------------------------------------------------------------------
!     SOLVE A LINEAR SYSTEM BY DIRECT FACTORIZATION
!     A*X = L*U*X = B
!     L(LOWER TRIANGULAR MATRIX WITH L(I,I)=1)
!     U(UPPER TRIANGULAR MATRIX),THE SUB-DIAGONAL TERMS OF L AND THE
!     TERMS U(I,J) ARE STORED IN THE SAME MATRIX, W. THE RESOLUTION
!     IS OBTAINED BY SOLVING TWO TRIANGULAR SYSTEMS:
!     L*Z = B  AND U*X = Z.
!     THE PARTIAL PIVOTING IS MADE BY CHOOSING THE BIGGEST ELEMENT OF
!     EACH COLUMN OF THE TRANSfor (MATION MATRIX.
!     INPUTS:
!     A = TABLE OF SIZE (N,N)                                      R*8
!     B = SECOND MEMBER VECTOR (N)                                 R*8
!     N = ORDER OF LINEAR SYSTEM                                   I*4
!     OUTPUTS:
!     X = SOLUTION VECTOR (N)                                      R*8
!     WORKING ZONE:
!     W = TABLE OF SIZE (LDA,N)                                    R*8
!     Z = AUXILIARY VECTOR (N)                                     R*8
!     NOTE:
!     MESSAGE '** DLITTL ** NO UNIQUE SOLUTION' IS GIVEN WHEN A 
!     NULL PIVOT IS FOUND.
!----------------------------------------------------------------------*/
//  Labels: e2,e10,e32
    double SUM,AMAX,T;
    int I,J,K,P;
      P=1;
      AMAX=fabs(A[1][1]);
//    SEEK MAX. ELEMENT OF FIRST COLUMN
      for (J=2; J<=N; J++) {
        if (fabs(A[J][1]) < AMAX) goto e2;
        AMAX=A[J][1];
        P=J;
e2:   ;}
      if (AMAX == 0.0) goto e32;
      if (P != 1) {
//    EXCHANGE LINES OF MATRIX A AND ORDER OF UNKNOWNS
//    LINKED TO B
		for (J=1; J<=N; J++) {
          T=A[1][J];
          A[1][J]=A[P][J];
          A[P][J]=T;
        }
        T=B[1];
        B[1]=B[P];
        B[P]=T;
      }
//    FIRST LINE OF U AND FIRST COLUMN OF L
      W[1][1]=A[1][1];
      for (J=2; J<=N; J++) {
        W[1][J]=A[1][J];
        W[J][1]=A[J][1]/W[1][1];
      }
//    SECOND LINE OF U AND SECOND COLUMN OF L
//    (N-1)TH LINE OF U AND (N-1)TH COLUMN OF L
      for (I=2; I<N; I++) {
		for (J=I; J<=N; J++) {
//    SEEK PIVOT
          AMAX=0.0;
          SUM=0.0;
          for (K=1; K<I; K++)
            SUM += W[J][K]*W[K][I];
          T=fabs(A[J][I]-SUM);
          if (T < AMAX) goto e10;
          AMAX=T;
          P=J;
e10:    ;}
        if (AMAX == 0.0) goto e32;
        if (P != I) {
//    EXCHANGE LINES OF MATRICES A , W AND ORDER OF UNKNOWNS
//    LINKED TO B
		  for (J=1; J<=N; J++) {
            T=A[I][J];
            A[I][J]=A[P][J];
            A[P][J]=T;
            T=W[I][J];
            W[I][J]=W[P][J];
            W[P][J]=T;
          }
          T=B[I];
          B[I]=B[P];
          B[P]=T;
        }
//    CALCULATE THE U(I,I)
        SUM=0.0;
        for (K=1; K<I; K++)
          SUM += W[I][K]*W[K][I];
        W[I][I]=A[I][I]-SUM;
        if (W[I][I] == 0.0) goto e32;
        for (J=I+1; J<=N; J++) {
          SUM=0.0;
          for (K=1; K<I; K++)
            SUM += W[I][K]*W[K][J];
//    CALCULATE THE U(I,J)
          W[I][J]=A[I][J]-SUM;
          SUM=0.0;
          for (K=1; K<I; K++)
            SUM += W[J][K]*W[K][I];
//    CALCULATE THE L(I,J)
          W[J][I]=(A[J][I]-SUM)/W[I][I];
        }
      }
      SUM=0.0;
      for (K=1; K<N; K++)
        SUM += W[N][K]*W[K][N];
//    CALCULATE  U(N,N)
      W[N][N]=A[N][N]-SUM;
      if (W[N][N] == 0.0) goto e32;
//    SOLVE SYSTEM  Z = L(-1)*B
      Z[1]=B[1];
      for (I=2; I<=N; I++) {
        SUM=0.0;
        for (K=1; K<I; K++)
          SUM += W[I][K]*Z[K];
        Z[I]=B[I]-SUM;
      }
//    SOLVE SYSTEM  X = U(-1)*Z
      X[N]=Z[N]/W[N][N];
      for (I=N-1; I>0; I--) {
        SUM=0.0;
        for (K=I+1; K<=N; K++)
          SUM += W[I][K]*X[K];
        X[I]=(Z[I]-SUM)/W[I][I];
      }
//    END OF RESOLUTION OF  A*X = L*(U*X) = B
      return;
e32:  printf("\n ** DLITTL ** NO UNIQUE SOLUTION.\n\n");
    }


void main()  {

  N=6;     //size of system

  for (I=1; I<=N; I++) {
    for (J=1; J<=N; J++) A[I][J]=0.0;
    B[I]=0.0;
  }

  // define elements of A <> 0
  A[1][1]= 1.0; A[1][6]= 1.0;
  A[2][1]= 1.0; A[2][2]= 1.0; A[2][6]=-1.0;
  A[3][1]=-1.0; A[3][2]= 1.0; A[3][3]= 1.0; A[3][6]= 1.0;
  A[4][1]= 1.0; A[4][2]=-1.0; A[4][3]= 1.0; A[4][4]= 1.0; A[4][6]=-1.0;
  A[5][1]=-1.0; A[5][2]= 1.0; A[5][3]=-1.0; A[5][4]= 1.0; A[5][5]= 1.0; A[5][6]= 1.0;
  A[6][1]= 1.0; A[6][2]=-1.0; A[6][3]= 1.0; A[6][4]=-1.0; A[6][5]= 1.0; A[6][6]=-1.0;

  //define elements of B <> 0
  B[1]=1.0; B[3]=1.0; B[5]=1.0;

  // print data
  printf("\n LINEAR SYSTEM AX=B:\n\n");
  for (I=1; I<=N; I++) {
    for (J=1; J<=N; J++) printf("%8.4f", A[I][J]);
    printf("   %8.4f\n", B[I]);
  }

  // call main function  
  DLITTL(A,B,N,X,W,Z);

  //print results
  printf("\n SOLUTION:\n\n");
  for (I=1; I<=N; I++) printf("  %f\n", X[I]);
  printf("\n");
  
}

// End of file tdlittl.cpp