cholesky1.c

来自「这是数值方法中各种插值算法」· C语言 代码 · 共 141 行

#include <stdio.h>
#include <math.h>
#define N 3
   int main()
       {
          int n1=3;
          double a1[3][3]={3,2,3,2,2,0,3,0,12};
          double b[3]={5,3,7};
          void Doolittle(int n,double a[N][N],double b[N]);
          Doolittle(n1,a1,b);
          getchar();

        }



  void  Doolittle(int n,double a[N][N],double b[N])
       {
          int i,j,k;
          double x[N],y[N];
          double sum,sum2;
          double l[N][N];
          double u[N][N];
          for(i=0;i<n;i++)
            for(j=0;j<n;j++)
              {
               l[i][j]=0;
               u[i][j]=0;
              }
          for(i=0;i<n;i++)
           {
            if(i==0)
             {
               l[i][i]=sqrt(a[i][i]);
              }

               l[i][0]=a[i][0]/l[0][0];

            }

           for(i=1;i<n;i++)
              {
                    sum=0;

                    for(k=0;k<i;k++)
                      {

                        sum=sum+l[i][k]*l[i][k];

                      }

                        l[i][i]=sqrt(a[i][i]-sum);

                for(j=i+1;j<n;j++)
                  {
                     sum2=0;

                    for(k=0;k<i;k++)
                      {

                        sum2=sum2+l[j][k]*l[i][k];

                       }


                       l[j][i]=(a[j][i]-sum2)/l[i][i];

                   }

                   }



              for(i=0;i<n;i++)
                for(j=0;j<n;j++)
                  {

                     u[j][i]=l[i][j];

                   }

              for(i=0;i<n;i++)
                {
                  sum=0;
                  if(i==0)
                    {
                      y[i]=b[i]/l[i][i];

                     }
                   else
                     {
                       for(k=0;k<i;k++)
                         {

                           sum=sum+l[i][k]*y[k];


                         }

                         y[i]=(b[i]-sum)/l[i][i];
                      }

                     }
          printf("The result of x is:\n");

               for(i=n-1;i>=0;i--)
                 {
                   sum=0;

                   if(i==n-1)
                     {
                       x[i]=y[i]/l[i][i];

                      }
                   else
                       {
                         for(k=i;k<n;k++)
                          {

                             sum=sum+u[i][k]*x[k];


                           }

                           x[i]=(y[i]-sum)/u[i][i];





                         }

                        printf("%f\n",x[i]);

                   }



        }