C語言 例題6-1 利用Guassian Elimination 解行列式

C語言 例題6-1 利用Guassian Elimination 解行列式
      input matrix
    -x1 + x2 +  2x3 = 2
    3x1 - x2 +   x3 = 6
    -x1 + 3x2 + 4x3 = 4

/* ex6-1.c based on Gaussian Elimination method
 * for solving the n x n linear algebra system
 * a11 x1+a12 x2+...+a1n xn=b1
 * a21 x1+a22 x2+...+a2n xn=b2
 * .       .         .       .
 * .       .         .       .
 * an1 x1+an2 x2+...+ann xn=bn
 * Input number of unknowns and equations n
 * with coefficent a11,a12,...,ann and b1,b2,
 * ...bn. Output solution x1,x2,x3,...,xn.
 *
      input matrix
    -x1 + x2 +  2x3 = 2
    3x1 - x2 +   x3 = 6
    -x1 + 3x2 + 4x3 = 4
*/

#include <stdio.h>
#include <math.h>
#define MAX 3
void gaussh(int n,double a[MAX][MAX+1],double x[]);
void printmatrix(int n , double a[MAX][MAX+1]) ;
void backSub(double a[MAX][MAX+1]) ;

int main()
{   
    int i,j,k,m,n;
    n=MAX;
    double x[MAX];
    double a[MAX][MAX+1]= { {-1.0 ,   1.0  , 2.0 ,  2.0},
                            { 3.0 ,  -1.0  , 1.0 ,  6.0},
                            {-1.0 ,   3.0  , 4.0 ,  4.0}
                        };

    printf("原始行列式\n");
    printf("================================\n");
    printmatrix(n,a);
    gaussh(n,a,x); /* call the function gaussh() */
    return 0;
}


void gaussh(int n,double a[MAX][MAX+1],double x[])
{
    int i,j,k,m;
    double temp,bb,cc;
    for(k=0;k<=n-1;k++)
    {
       /* check if a[k][k]=0 is true then interchange */
       /* E(k) and E(k+1).............................*/
       if(a[k][k]==0)
       {
        for(m=0;m<=n;m++)
        {
            temp=a[k][m];
            a[k][m]=a[k+1][m];
            a[k+1][m]=temp;
        }
       }
     
       printf("\n");
       printf("================================\n");
       printmatrix(n,a);
     
       /* To reduce the matrix to triangular form */
       for(i=k;i<n-1;i++)
       {
        bb=a[i+1][k]/a[k][k];
        printf("i=%1d--: %3.2lf , ",i,bb);
        for(j=k;j<=n+1;j++)
            a[i+1][j]=a[i+1][j]-bb*a[k][j];
       }
    }
   
   
    if(fabs(a[n-1][n-1])==0.0)
    {
        printf("NO UNIQUE SOLUTION!!!\n");
        return;
    }
   

   
    printf("\n\n上三角矩陣\n");
    printf("================================\n");
    printmatrix(n,a);
   
    backSub(a);
    return;
}

void printmatrix(int n , double a[MAX][MAX+1])
{
    for (int i=0; i<n; i++)
    {
        for (int j=0; j<=n; j++)
            printf("%4.2lf ", a[i][j]);
        printf("\n");
    }       
}


// function to calculate the values of the unknowns
void backSub(double a[MAX][MAX+1])
{
    double x[MAX];  // An array to store solution
 
    /* Start calculating from last equation up to the
       first */
    for (int i = MAX-1; i >= 0; i--)
    {
        /* start with the RHS of the equation */
        x[i] = a[i][MAX];
 
        /* Initialize j to i+1 since matrix is upper
           triangular*/
        for (int j=i+1; j<MAX; j++)
        {
            /* subtract all the lhs values
             * except the coefficient of the variable
             * whose value is being calculated */
            x[i] -= a[i][j]*x[j];
        }
 
        /* divide the RHS by the coefficient of the
           unknown being calculated */
        x[i] = x[i]/a[i][i];
    }
   
    printf("\n行列式 解 \n");
    printf("================================\n");
    for (int i=0; i<MAX; i++)
        printf("x%1d = %3.3lf\n", i, x[i]);
}

輸出畫面
原始行列式
================================
-1.00 1.00 2.00 2.00
3.00 -1.00 1.00 6.00
-1.00 3.00 4.00 4.00

================================
-1.00 1.00 2.00 2.00
3.00 -1.00 1.00 6.00
-1.00 3.00 4.00 4.00
i=0--: -3.00 , i=1--: 1.00 ,
================================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 2.00 2.00 2.00
i=1--: 1.00 ,
================================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 0.00 -5.00 -10.00


上三角矩陣
================================
-1.00 1.00 2.00 2.00
0.00 2.00 7.00 12.00
0.00 0.00 -5.00 -10.00

行列式 解
================================
x0 = 1.000
x1 = -1.000
x2 = 2.000