crout_solve.c

Solve Ax=B with Crout s method

#include <stdio.h>

main()
{
      int i,j,k;
      int A[3][3]={{1,-1,0},{2,2,3},{-1,3,2}},B[3]={2,-1,4};
      double temp;
      static double L[3][3];
      static double U[3][3];
      static double X[3],Y[3];
      for (i=0;i<3;i++)
      {
          L[i][0]=A[i][0];   
          U[0][i]=A[0][i]/L[0][0];      
      }
      for (j=1;j<3;j++)
      {
          U[j][j]=1;
          for (i=j;i<3;i++)
          {
              temp=0.0;
              for (k=0;k<=j-1;k++)
              {
                  temp+=L[i][k]*U[k][j];
              }
              L[i][j]=A[i][j]-temp;
          }
          for (i=j+1;i<3;i++)
          {
              temp=0.0;
              for (k=0;k<=j+1;k++)
              {
                  temp+=L[j][k]*U[k][i];
              }
              U[j][i]=(A[j][i]-temp)/L[j][j];
          }
      }
      Y[0]=B[0]/L[0][0];
      for (i=1;i<3;i++)
      {
         temp=0.0;
         Y[i]=B[i]/L[i][i];
         for (j=0;j<i;j++)
         {
             temp+=L[i][j]*Y[j];
         }
         Y[i]=(B[i]-temp)/L[i][i];
      }    
      X[2]=Y[2]/U[2][2];
      for (i=1;i>=0;i--)
      {
         temp=0.0;
         X[i]=Y[i];
         for (j=i+1;j<3;j++)
         {   
             temp+=U[i][j]*X[j];
         }
         X[i]=(Y[i]-temp)/U[i][i];
      }
      printf("The solution of A is as follows:\n");
      printf("\n");
      for (i=0;i<3;i++)
      {
          if (X[i]>=0)
             printf(" %c %c%f %c \n",'[','+',X[i],']');
          else
             printf(" %c %f %c \n",'[',X[i],']');
      }
      printf("\n");
      getch();
}