C語言 例題6-8利用LU分解法求線性方程式的解

C  code program for LU Decomposition for solving linear equations
C語言 例題6-8利用LU分解法求線性方程式的解

/* ex6-8.c based on Gaussian Elimination Method
 * to reduce a determinant to upper triangular
 * determinant, then get the result form it.
 */
#include <stdio.h>
#include <math.h>

void printfmatrix(int n,double a[n][n])
{
    //print matrix
    for(int i=0;i<n;i++)
    {
        for(int j=0;j<n;j++)
        printf("%7.3lf\t\t",a[i][j]);
        printf("\n");
    }
}

double dett(int n,double a[n][n])
{
    int i,j,k,m,flag=0; /* check number of row interchange */
    double temp,bb,acc=1.0;
 
    for(k=0;k<n;k++)
    {
        /* chech if a[k][k]=0 is true,then
        * interchange E(k) and E(k+1)
        */
        if(a[k][k]==0)
        {
        flag++;
        for(m=1;m<n;m++)
        {
            temp=a[k][m];
            a[k][m]=a[k+1][m];
            a[k+1][m]=temp;
        }
        }
        /* To reduce the matrix to triangular form */
        for(i=k;i<n-1;i++)
        {
        bb=a[i+1][k]/a[k][k];
        printf("factor=%0.3lf\n",bb);
        for(j=k;j<n;j++)
        {
            a[i+1][j]=a[i+1][j]-bb*a[k][j];
            printf("a[%1d][%1d]=%0.3lf\t\t",i+1,j,a[i+1][j]);
        }
        printf("\n");
        printfmatrix(n,a);
        printf("\n");
       
        }
       
        //print matrix
    }
   
    if(fabs(a[n-1][n-1]==0.0))
    {
        printf("The determinant=0\n");
        return 0;
    }
 
    printf("Upper triangular determinant\n");
    printfmatrix(n,a);
    /* multiplying the diagonal element  */
    /* and modifying the answer based on */
    /* flag's odd or even number ....... */
    for(i=0;i<n;i++)
        acc=acc*a[i][i];
    if(flag%2==0)
        return (acc);
    else
        return (-acc);
}



int main()
{
    int i,j,k,m,n;
    n=3;
    double a[3][3] = {  {4,-2,1},
                        {3,0,-5},
                        {1,-3,-4}};
    printf("n=%d\n",n);
    printf("The Matrix = \n");
    printfmatrix(n,a);
    printf("\n");
    /* call function dett() and output the result */
    printf("Ans=%8.3lf\n",dett(n, a));
    return 0;
}

輸出畫面
n=3
The Matrix =
  4.000 -2.000   1.000
  3.000   0.000 -5.000
  1.000 -3.000 -4.000

factor=0.750
a[1][0]=0.000 a[1][1]=1.500 a[1][2]=-5.750
  4.000 -2.000   1.000
  0.000   1.500 -5.750
  1.000 -3.000 -4.000

factor=0.250
a[2][0]=0.000 a[2][1]=-2.500 a[2][2]=-4.250
  4.000 -2.000   1.000
  0.000   1.500 -5.750
  0.000 -2.500 -4.250

factor=-1.667
a[2][1]=0.000 a[2][2]=-13.833
  4.000 -2.000   1.000
  0.000   1.500 -5.750
  0.000   0.000 -13.833

Upper triangular determinant
  4.000 -2.000   1.000
  0.000   1.500 -5.750
  0.000   0.000 -13.833
Ans= -83.000