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


 main() {

   double sumatoria;
   int i,j,k,tam,s;
   double a[4][4];
   double u[4][4];
   double l[4][4];
   
   
   
   
   
   tam=4;
   
   
   
   
   a[0][0]=6.0;
   a[0][1]=2.0;
   a[0][2]=1.0;
   a[0][3]=-1.0;
   a[1][0]=2.0;
   a[1][1]=4.0;
   a[1][2]=1.0;
   a[1][3]=0.0;
   a[2][0]=1.0;
   a[2][1]=1.0;
   a[2][2]=4.0;
   a[2][3]=-1.0;
   a[3][0]=-1.0;
   a[3][1]=0.0;
   a[3][2]=-1.0;
   a[3][3]=3.0;

   


   l[0][0]=0.0;
   l[0][1]=0.0;
   l[0][2]=0.0;
   l[0][3]=0.0;
   l[1][0]=0.0;
   l[1][1]=0.0;
   l[1][2]=0.0;
   l[1][3]=0.0;
   l[2][0]=0.0;
   l[2][1]=0.0;
   l[2][2]=0.0;
   l[2][3]=0.0;
   l[3][0]=0.0;
   l[3][1]=0.0;
   l[3][2]=0.0;
   l[3][3]=0.0;


   u[0][0]=0.0;
   u[0][1]=0.0;
   u[0][2]=0.0;
   u[0][3]=0.0;
   u[1][0]=0.0;
   u[1][1]=0.0;
   u[1][2]=0.0;
   u[1][3]=0.0;
   u[2][0]=0.0;
   u[2][1]=0.0;
   u[2][2]=0.0;
   u[2][3]=0.0;
   u[3][0]=0.0;
   u[3][1]=0.0;
   u[3][2]=0.0;
   u[3][3]=0.0;





   for (i = 0; i < tam; i++)  l[i][i] = 1.0;

  

   for ( k = 0; k < tam; k++) {
  
   
     for ( j = 0; j < tam; j++) {
 	  
     
                      sumatoria = 0.0;
                      for ( s = 0; s < tam; s++) {
 
	                             sumatoria += l[k][s] * u[s][j];
				     
				     
	                            if (k==tam-1 && j==1) {

printf("cout << ERROR   << sumatoria << endl;\n");
printf("cout << l[ << k << , << s << ] X u[ << s << , << j << ] =  << l[k][s]*u[s][j] << endl;cout << ERROR   << sumatoria << endl;\n");
     /* 
     cout << "ERROR  " << sumatoria << endl;
	  cout << "l[" << k << "," << s << "] X u[" << s << "," << j << "] = " << l[k][s]*u[s][j] << endl;
       */
	                                                  }

                                                    }
                    u[k][j] = (a[k][j] - sumatoria) / l[k][k];
		    
		    printf("u[%d][%d]=%e\n" ,k,j,u[k][j]);
                                   }  /* fin de j */

     for ( i = k+1; i < tam; i++) {
        sumatoria = 0.0;
 	          for ( s = 0; s < tam; s++) {
	          sumatoria += l[i][s] * u[s][k];
 	        			        }
       l[i][k] = (a[i][k] - sumatoria) / u[k][k];
       printf("l[%d][%d]=%e\n" ,i,k,l[i][k]);

                                    }

                                 }   /* fin de i */

 }   /* fin de main */





el_choya@pc-joel:~# gcc bimbo.c -o bimbo
el_choya@pc-joel:~# ./bimbo
u[0][0]=6.000000e+00
u[0][1]=2.000000e+00
u[0][2]=1.000000e+00
u[0][3]=-1.000000e+00
l[1][0]=3.333333e-01
l[2][0]=1.666667e-01
l[3][0]=-1.666667e-01
u[1][0]=0.000000e+00
u[1][1]=3.333333e+00
u[1][2]=6.666667e-01
u[1][3]=3.333333e-01
l[2][1]=2.000000e-01
l[3][1]=1.000000e-01
u[2][0]=0.000000e+00
u[2][1]=0.000000e+00
u[2][2]=3.700000e+00
u[2][3]=-9.000000e-01
l[3][2]=-2.432432e-01
u[3][0]=0.000000e+00
cout << ERROR   << sumatoria << endl;
cout << l[ << k << , << s << ] X u[ << s << , << j << ] =  << l[k][s]*u[s][j] << endl;cout << ERROR   << sumatoria << endl;
cout << ERROR   << sumatoria << endl;
cout << l[ << k << , << s << ] X u[ << s << , << j << ] =  << l[k][s]*u[s][j] << endl;cout << ERROR   << sumatoria << endl;
cout << ERROR   << sumatoria << endl;
cout << l[ << k << , << s << ] X u[ << s << , << j << ] =  << l[k][s]*u[s][j] << endl;cout << ERROR   << sumatoria << endl;
cout << ERROR   << sumatoria << endl;
cout << l[ << k << , << s << ] X u[ << s << , << j << ] =  << l[k][s]*u[s][j] << endl;cout << ERROR   << sumatoria << endl;
u[3][1]=-5.556536e-18
u[3][2]=0.000000e+00
u[3][3]=2.581081e+00
el_choya@pc-joel:~#                                   


:)