=============================================================
		C code
=============================================================

#include "stdafx.h"
#include "math.h"

int _tmain(int argc, _TCHAR* argv[])
{
   int i, j, k, p;

   double pmax, s;

   double a[8][8] ={2.0373, 0.2961, 1.4826, 0.0915, 2.4165, 2.9361, 1.5648, 2.919,
		    1.1865, 0.7857, 2.3373, 2.2323, 1.7301, 2.1381, 0.2901, 1.947,
		    1.1022, 1.0062, 2.145,  1.5,    0.5487, 1.5015, 2.4543, 2.4009,
		    2.964,  2.0391, 2.7111, 1.4397, 0.7197, 1.4133, 2.4525, 1.3614,
		    0.1131, 0.4098, 2.6727, 2.7141, 2.6595, 0.1788, 2.1672, 1.2972,
		    2.6556, 2.1636, 1.0026, 1.8297, 0.0861, 2.046,  0.4497, 2.4759,
		    2.7399, 0.3204, 2.0961, 1.8531, 1.4697, 0.1272, 1.9788, 0.2505,
		    2.3886, 1.9614, 0.5934, 2.5782, 0.5037, 0.2142, 1.5558, 0.3996};

   double b[8] = {9.5974,
		  3.4039,
		  5.8527,
		  2.2381,
		  7.5127,
		  2.5510,
		  5.0596,
		  6.9908};

   double w[8];
   double y[8];
   double deta, detu;
   double A[8][8];
   double AC[8][8];
   double u1[8][8];
   double u2[8][8];
   double l1[8][8];
   double l2[8][8];
   double bC[8];


/* ===================== Gaussian Elimination Method ===================== */


/* AC, A = copy of A */
   	for(i=0; i<8; i++){
		  for(j=0; j<8; j++){
			  AC[i][j]=0;
			  A[i][j]=0;
			  AC[i][j]=a[i][j];
			  A[i][j]=a[i][j];
		  }
	}

/* bC = copy of b */
   	for(i=0; i<8; i++){
			  bC[i]=0;
			  bC[i]=b[i];
	}

/* 1(1.1) */

/* find the pivot */
   for(k = 0; k < 8-1; k++){ 
      p = k;
      pmax = fabs( a[k][k] );
      for(i = k+1; i < 8; i++){
         if(fabs( a[i][k] ) < pmax){
            p = i;
            pmax = fabs( a[i][k] );
         }
      }

/* if the pivot is too small */
      if(fabs( pmax ) < 1.0e-12){
         printf("too small pivot! \n");
         return(0);
      }

/* swap the rows */
      if(p != k){
         for(i = k; i < 8; i++){
            s = a[k][i];
            a[k][i] = a[p][i];
            a[p][i] = s;
         }
         s = b[k];
         b[k] = b[p];
         b[p] = s;
      }

      for(i = k +1; i < 8; i++){
         w[i] = a[i][k] / a[k][k];
         a[i][k] = 0;
         for(j = k + 1; j < 8; j++){
            a[i][j] = a[i][j] - a[k][j] * w[i];
         }
         b[i] = b[i] - b[k] * w[i];
      }
   }

/* backward substitution */
      for(i = 8 - 1; i >= 0; i--){
         for(j = i + 1; j < 8; j++){
            b[i] = b[i] - a[i][j] * b[j];
         }
         b[i] = b[i] / a[i][i];
      }

	printf("1(1.1) \n");
	printf("x= \n");
      for(i = 0; i < 8; i++){
         printf("%12.8lf\n",b[i]);
      }

/* 1(1.2) */
	for(i=0; i<8; i++){
		y[i]=0;
		for (k = 0; k < 8; k++) {
			y[i]+=AC[i][k]*b[k];
		}
	}

	printf("1(1.2) \n");
	printf("Ax= \n");
	for(i=0;i<8;i++) {
        printf("%8.4f",y[i]);
                printf("\n");       
        }
	
/* 1(1.3) */		
	deta=1;
	for(i=0; i<8; i++){
		deta *= a[i][i];
	}

	printf("1(1.3) \n");
	printf("detA(diagonal) = %13.10f\n", deta);

/* 2(2.1) */
	for( i = 0; i < 8; i++ ){
		for( j = 0; j < 8; j++ ){
			if( i == j ){ l1[i][j] = 1; }
			else{ l1[i][j] = 0; }
			u1[i][j] = AC[i][j];
		}
	}

	for( i = 0; i < 8; i++ ){
		for( j = i+1; j < 8; j++ ){
			l1[j][i] = u1[j][i] / u1[i][i];
			for( k = i; k < 8; k++ ){
				u1[j][k] = u1[j][k] - l1[j][i]*u1[i][k];
			}
		}
	}

	printf("2(2.1) \n");
	printf("Matrix L: \n");
	for(i=0; i<8; i++){
		for(j=0; j<8; j++){
			printf(" %7.4f", l1[i][j]);
		}
		printf("\n");
	}

	printf("Matrix U: \n");
	for(i=0; i<8; i++){
		for(j=0; j<8; j++){
			printf(" %7.4f", u1[i][j]);
		}
		printf("\n");
	}

/* 2(2.2) */
	for( i = 0; i < 8; i++ ){
		for( j = i+1; j <8; j++ ){
			bC[j] -= bC[i] * l1[j][i];
		}
	}
	for( i = 8-1; i >= 0; i-- ){
		for( j = i+1; j<8; j++ ){
			bC[i] -= u1[i][j] * bC[j];
		}
		bC[i] /= u1[i][i];
	} 

	printf("2(2.2) \n");
	printf("x= \n");
	for(i=0; i<8; i++){
		printf(" %12.8lf", bC[i]);
		printf("\n");
	}


/* ===================== LU Factorization ===================== */


/* 3(3.1) */
	for (int i = 0; i < 8; ++i) {
		for(int j = i+1; j < 8; ++j) {
			A[j][i] /= A[i][i];
		for(int k = i+1; k < 8; ++k)
			A[j][k] -= A[i][k] * A[j][i];
		}
	}

	for( i = 0; i < 8; i++ ){
		for( j = 0; j < 8; j++ ){
			if( j > i-1 ){ u2[i][j] = A[i][j];}
			else{ u2[i][j] = 0;}
		}
	}

	printf("3(3.1) \n");
	printf("Matrix U \n");
	for (i=0 ; i<8; i++ ){
		for (j=0; j<8 ; j++ ){
			printf(" %7.4f", u2[i][j]);
       		}
		printf("\n");
	}

	 for( i = 0; i < 8; i++ ){
		for( j = 0; j < 8; j++ ){
			if( i == j ){ l2[i][j] = 1; }
			else if( j < i ){ l2[i][j] = A[i][j];}
			else{ l2[i][j] = 0;}
		}
	}
	printf("Matrix L \n");
	for (i=0 ; i<8; i++ ){
		for (j=0; j<8 ; j++ ){
			printf(" %7.4f", l2[i][j]);
        	}
		printf("\n");
	}

/* 3(3.2) */
	detu=1;
	for(i=0; i<8; i++){
		detu *= u2[i][i];
	}

	printf("3(3.2) \n");
	printf("detU(diagonal) = %13.10f\n", detu);
	
    return(0);
}


=============================================================
			Results
=============================================================

1(1.1)
x=
 -0.12993705
 -1.49890981
 -3.78770710
  2.27071756
  1.74748010
  2.43346205
  3.45707911
 -0.36432996
______________________________________________

1(1.2)
Ax=
  9.5974
  3.4039
  5.8527
  2.2381
  7.5127
  2.5510
  5.0596
  6.9908

Error:  A*x-b = 

  1.0e-014 *

    0.1776
         0
         0
         0
   -0.2665
    0.0888
   -0.0888
         0

______________________________________________

1(1.3)
detA(diagonal) = 377.7375504427
______________________________________________

2(2.1)
Matrix L:
  1.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.5824  1.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.5410  1.3795  1.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  1.4549  2.6226  4.7965  1.0000  0.0000  0.0000  0.0000  0.0000
  0.0555  0.6414 -2.3829 -0.7846  1.0000  0.0000  0.0000  0.0000
  1.3035  2.8987  7.5358  2.3309  0.0921  1.0000  0.0000  0.0000
  1.3449 -0.1269 -0.4189  0.4438 -2.8421 -2.6615  1.0000  0.0000
  1.1724  2.6323  7.2782  2.6389 -0.0385  0.6009  0.3482  1.0000
Matrix U:
  2.0373  0.2961  1.4826  0.0915  2.4165  2.9361  1.5648  2.9190
  0.0000  0.6133  1.4739  2.1790  0.3228  0.4281 -0.6212  0.2470
  0.0000  0.0000 -0.6903 -1.5555 -1.2039 -0.6776  2.4647  0.4809
  0.0000  0.0000  0.0000  3.0530  2.1321 -0.7310 -10.0170 -5.8400
  0.0000  0.0000  0.0000  0.0000  1.1223 -2.4471  0.4929 -2.4592
  0.0000  0.0000  0.0000  0.0000  0.0000  4.0133  4.9400  8.1695
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 19.8220 13.9030
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -0.0000 -1.6067
______________________________________________

2(2.2)
x=
  -0.12993705
  -1.49890981
  -3.78770710
   2.27071756
   1.74748010
   2.43346205
   3.45707911
  -0.36432996
______________________________________________

3(3.1)
Matrix U
  2.0373  0.2961  1.4826  0.0915  2.4165  2.9361  1.5648  2.9190
  0.0000  0.6133  1.4739  2.1790  0.3228  0.4281 -0.6212  0.2470
  0.0000  0.0000 -0.6903 -1.5555 -1.2039 -0.6776  2.4647  0.4809
  0.0000  0.0000  0.0000  3.0530  2.1321 -0.7310 -10.0170 -5.8400
  0.0000  0.0000  0.0000  0.0000  1.1223 -2.4471  0.4929 -2.4592
  0.0000  0.0000  0.0000  0.0000  0.0000  4.0133  4.9400  8.1695
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 19.8220 13.9030
  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000 -1.6067
Matrix L
  1.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.5824  1.0000  0.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  0.5410  1.3795  1.0000  0.0000  0.0000  0.0000  0.0000  0.0000
  1.4549  2.6226  4.7965  1.0000  0.0000  0.0000  0.0000  0.0000
  0.0555  0.6414 -2.3829 -0.7846  1.0000  0.0000  0.0000  0.0000
  1.3035  2.8987  7.5358  2.3309  0.0921  1.0000  0.0000  0.0000
  1.3449 -0.1269 -0.4189  0.4438 -2.8421 -2.6615  1.0000  0.0000
  1.1724  2.6323  7.2782  2.6389 -0.0385  0.6009  0.3482  1.0000
______________________________________________

3(3.2)
detU(diagonal) = 377.7375504427