C語言 例題6-1 利用高斯消去法 求行列式的解

C語言 例題6-1 利用高斯消去法 求行列式的解

參考  https://www.geeksforgeeks.org/gaussian-elimination/

// C program to demostrate working of Guassian Elimination
// method
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define N 3        // Number of unknowns

// function to reduce matrix to r.e.f.  Returns a value to
// indicate whether matrix is singular or not
int forwardElim(double mat[N][N+1]);

// function to calculate the values of the unknowns
void backSub(double mat[N][N+1]);

// function to get matrix content
void gaussianElimination(double mat[N][N+1])
{
    /* reduction into r.e.f. */
    int singular_flag = forwardElim(mat);

    /* if matrix is singular */
    if (singular_flag != -1)
    {
        printf("Singular Matrix.\n");

        /* if the RHS of equation corresponding to
           zero row  is 0, * system has infinitely
           many solutions, else inconsistent*/
        if (mat[singular_flag][N])
            printf("Inconsistent System.");
        else
            printf("May have infinitely many "
                   "solutions.");

        return;
    }

    /* get solution to system and print it using
       backward substitution */
    backSub(mat);
}

// function for elemntary operation of swapping two rows
void swap_row(double mat[N][N+1], int i, int j)
{
    //printf("Swapped rows %d and %d\n", i, j);

    for (int k=0; k<=N; k++)
    {
        double temp = mat[i][k];
        mat[i][k] = mat[j][k];
        mat[j][k] = temp;
    }
}

// function to print matrix content at any stage
void print(double mat[N][N+1])
{
    for (int i=0; i<N; i++, printf("\n"))
        for (int j=0; j<=N; j++)
            printf("%lf ", mat[i][j]);

    printf("\n");
}

// function to reduce matrix to r.e.f.
int forwardElim(double mat[N][N+1])
{
    for (int k=0; k<N; k++)
    {
        // Initialize maximum value and index for pivot
        int i_max = k;
        int v_max = mat[i_max][k];

        /* find greater amplitude for pivot if any */
        for (int i = k+1; i < N; i++)
            if (abs(mat[i][k]) > v_max)
                v_max = mat[i][k], i_max = i;

        /* if a prinicipal diagonal element  is zero,
         * it denotes that matrix is singular, and
         * will lead to a division-by-zero later. */
        if (!mat[k][i_max])
            return k; // Matrix is singular

        /* Swap the greatest value row with current row */
        if (i_max != k)
            swap_row(mat, k, i_max);
         
        printf("第 %1d 次的行列式 \n",k);
        print(mat);
     
        for (int i=k+1; i<N; i++)
        {
            /* factor f to set current row kth elemnt to 0,
             * and subsequently remaining kth column to 0 */
            double f = mat[i][k]/mat[k][k];

            /* subtract fth multiple of corresponding kth
               row element*/
            for (int j=k+1; j<=N; j++)
                mat[i][j] -= mat[k][j]*f;

            /* filling lower triangular matrix with zeros*/
            mat[i][k] = 0;
        }

        print(mat);        //for matrix state
    }
    printf("上三角矩陣\n");
    print(mat);            //for matrix state
    return -1;
}

// function to calculate the values of the unknowns
void backSub(double mat[N][N+1])
{
    double x[N];  // An array to store solution

    /* Start calculating from last equation up to the
       first */
    for (int i = N-1; i >= 0; i--)
    {
        /* start with the RHS of the equation */
        x[i] = mat[i][N];

        /* Initialize j to i+1 since matrix is upper
           triangular*/
        for (int j=i+1; j<N; j++)
        {
            /* subtract all the lhs values
             * except the coefficient of the variable
             * whose value is being calculated */
            x[i] -= mat[i][j]*x[j];
        }

        /* divide the RHS by the coefficient of the
           unknown being calculated */
        x[i] = x[i]/mat[i][i];
    }

    printf("\nSolution for the system:\n");
    for (int i=0; i<N; i++)
        printf("%lf\n", x[i]);
}

// Driver program
int main()
{
    /* input matrix */
    double mat[N][N+1] = {{3.0, 2.0,-4.0, 3.0},
                          {2.0, 3.0, 3.0, 15.0},
                          {5.0, -3, 1.0, 14.0}
                         };
 
    printf("原始行列式\n");
    print(mat);            //for matrix state
    gaussianElimination(mat);

    return 0;
}

輸出畫面
原始行列式
3.000000 2.000000 -4.000000 3.000000
2.000000 3.000000 3.000000 15.000000
5.000000 -3.000000 1.000000 14.000000

第 0 次的行列式
5.000000 -3.000000 1.000000 14.000000
2.000000 3.000000 3.000000 15.000000
3.000000 2.000000 -4.000000 3.000000

5.000000 -3.000000 1.000000 14.000000
0.000000 4.200000 2.600000 9.400000
0.000000 3.800000 -4.600000 -5.400000

第 1 次的行列式
5.000000 -3.000000 1.000000 14.000000
0.000000 4.200000 2.600000 9.400000
0.000000 3.800000 -4.600000 -5.400000

5.000000 -3.000000 1.000000 14.000000
0.000000 4.200000 2.600000 9.400000
0.000000 0.000000 -6.952381 -13.904762

第 2 次的行列式
5.000000 -3.000000 1.000000 14.000000
0.000000 4.200000 2.600000 9.400000
0.000000 0.000000 -6.952381 -13.904762

5.000000 -3.000000 1.000000 14.000000
0.000000 4.200000 2.600000 9.400000
0.000000 0.000000 -6.952381 -13.904762

上三角矩陣
5.000000 -3.000000 1.000000 14.000000
0.000000 4.200000 2.600000 9.400000
0.000000 0.000000 -6.952381 -13.904762


Solution for the system:
3.000000
1.000000
2.000000