/********************************************************
*   SOLVE A LINEAR SYSTEM BY TRIANGULARIZATION METHOD   *
* ----------------------------------------------------- *
* SAMPLE RUN:                                           *
*                                                       *
* Linear system AX = B                                  *
*                                                       *
* Matrix A:                                             *
*   2.00   1.00  -2.00                                  *
*   3.00  -1.00   1.00                                  *
*   7.00   5.00  -3.00                                  *
*                                                       *
* Right side:                                           *
*   B(1) =   1.00                                       *
*   B(2) =   0.00                                       *
*   B(3) =   0.00                                       *
*                                                       *
* Solution of AX=B:                                     *
*   X(1) =   0.0625                                     *
*   X(2) =  -0.5000                                     *
*   X(3) =  -0.6875                                     *
*                                                       *
* ----------------------------------------------------- *
* Reference: "Méthodes de calcul numérique - Tome 1 By  *
*             Claude Nowakowski, PS1 1981" [BIBLI 07].  *
*                                                       *
*                    C++ Release By J-P Moreau, Paris.  *
*                          (www.jpmoreau.fr)            *
********************************************************/
#include <stdio.h>

#define  SIZE  25

typedef double MAT[SIZE][SIZE];  // index 0 not used here
typedef double VEC[SIZE];

 MAT A;
 VEC B, X;

 int I, J, K, N;
 double S;


void main() {

  N = 3;

  A[1][1] = 2.0; A[1][2] =  1.0; A[1][3] = -2.0;
  A[2][1] = 3.0; A[2][2] = -1.0; A[2][3] =  1.0;
  A[3][1] = 7.0; A[3][2] =  5.0; A[3][3] = -3.0;

  B[1] = 1.0;
  B[2] = 0.0;
  B[3] = 0.0;

  printf("\n Linear system AX = B\n\n");
  printf(" Matrix A:\n");
  for (I=1; I<=N; I++) {
    for (J=1; J<=N; J++) printf(" %6.2f", A[I][J]);
    printf("\n");
  }

  printf("\n Right Side:\n");
  for (I=1; I<=N; I++) printf("   B(%d) = %6.2f\n", I, B[I]);

// Transform A into triangular matrix
  for (K=1; K<N; K++) {
    for (I=K+1; I<=N; I++) {
      B[I] -= A[I][K] / A[K][K] * B[K];
      for (J=K+1; J<=N; J++)
        A[I][J] -= A[I][K] / A[K][K] * A[K][J];
    }
  }

// Solve triangular system
  X[N] = B[N] / A[N][N];
  for (I=N-1; I>0; I--) {
    S = 0.0;
    for (K=I+1; K<=N; K++)  S += A[I][K] * X[K];
    X[I] = (B[I] - S) / A[I][I];
  }

// Print results
  printf("\n Solution of AX=B:\n");
  for (I=1; I<=N; I++)  printf("   X(%d) = %8.4f\n", I, X[I]);
  printf("\n");

}

// end of file Tlinear.cpp