닫기

/**
Project Name : Find NxN inverse Matrix(A) ( P01.c )
use : Report.

Version : 1.0
Copyright (c) 2007 : Ohyung ( ohyung@ohyung.com )
Last modified at : 2007.11.07

Base  : linux (ubuntu 7.10 - Gutsy Gilbbon, kernel 2.6.22-12-generic)
compiler : gcc (GCC) 4.1.3
Tool : Code::Blocks svn4596

사용방법 : $ gcc -o P01 P01.c
실행 : $ ./P01
*/

#include < stdio.h>
#include < stdlib.h> // malloc

static double
Determinant(double** matrixA, int uesrN);

static int
GaussElemination(double** matrixA, double** matrixB, int userN);

static int
Pivoting(double** matrixA, double** matrixB, int count, int userN);

static double
Absolute(double value);

int
main(void)
{
    int userN = 0;
    int userM = 0;
    int row = 0;
    int column = 0;

    double** matrixA;
    double** matrixI=0;
    double determinantResult = 0;

    printf("n by n = ");
    scanf("%d", &userN);

    userM = userN - 1;

    matrixA = (double**)malloc(sizeof(double*) * userN);
    matrixI = (double**)malloc(sizeof(double*) * userN);

    for (row = 0; row < userN; row++)
    {
        matrixA[row] = (double*)malloc(sizeof(double) * userN);
        matrixI[row] = (double*)malloc(sizeof(double) * userN);
    } /* for (row = 0; row < userN; row++) */

    for (row = 0; row < userN; row++)
    {
        matrixI[row][row] = 1;
    }

    for (row = 0; row < userN; row++)
    {
        for (column = 0; column < userN; column++)
        {
            printf("%d , %d =", row, column);
            scanf("%lf", &matrixA[row][column]);
        } /* for (column = 0; column < userN; column++) */
    } /* for (row = 0; row < userN; row++) */
    printf("\nA =\n");

    for (row = 0; row < userN; row++)
    {
        for (column = 0; column < userN; column++)
        {
            printf("%10.5lf ", matrixA[row][column]);
            if(column % userN == userM)
            {
                printf("\n");
            } /* if(column % userN == userM) */
        } /* for (column = 0; column < userN; column++) */
    } /* for (row = 0; row < userN; row++) */

    printf("\nI =\n");
    for (row = 0; row < userN; row++)
    {
        for (column = 0; column < userN; column++)
        {
            printf("%10.5lf ", matrixI[row][column]);
            if(column % userN == userM)
            {
                printf("\n");
            } /* if(column % userN == userM) */
        } /* for (column = 0; column < userN; column++) */
    } /* for (row = 0; row < userN; row++) */


    determinantResult = Determinant(matrixA, userN);
    printf("\n|A| = %.5g\n", determinantResult);

    if (!determinantResult)
    {
        printf("\nNo determinant(A)\n");
    } /* if (!determinantResult) */
    else
    {
        GaussElemination(matrixA, matrixI, userN );

        // origin I to inversed A
        printf("\nA' =\n");
        for (row = 0; row < userN; row++)
        {
            for (column = 0; column < userN; column++)
            {
                printf("%10.5lf ", matrixI[row][column]);
                if(column % userN == userM)
                {
                    printf("\n");
                } /* if(column % userN == userM) */
            } /* for (column = 0; column < userN; column++) */
        } /* for (row = 0; row < userN; row++) */

        // origin A to I
        printf("\nAA' =\n");
        for (row = 0; row < userN; row++)
        {
            for (column = 0; column < userN; column++)
            {
                printf("%10.5lf ", matrixA[row][column]);
                if(column % userN == userM)
                {
                    printf("\n");
                } /* if(column % userN == userM) */
            } /* for (column = 0; column < userN; column++) */
        } /* for (row = 0; row < userN; row++) */

    } /* else if(determinantResult)*/

    for(row = 0; row < userN; row++)
    {
        free(matrixA[row]);
        free(matrixI[row]);
    } /* for(row = 0; row < userN; row++) */
    free(matrixA);
    free(matrixI);

    return 0;
}

/**
 * Function Name: Determinant
 *
 * Function Description:
 *        This function solve the determinant from a matrix A
 *
 * Input: matrixA , userN
 * Output: determinantResult
 * Version: 1.0
 */
static double
Determinant (double** matrixA, int userN)
{
    int userM = 0;
    int counter = 0;
    int determinantACheck = 0;
    int determinantARow = 0;
    int determinantAColumn = 0;
    int sign = 1;

    double determinantResult = 0;
    double** matrixDeterminantA;

    userM = userN-1;

    matrixDeterminantA = (double**)malloc(sizeof(double*) * userM);
    for (counter = 0; counter < userM; counter++)
    {
        matrixDeterminantA[counter] = (double*)malloc(sizeof(double) * userM);
    } /* for (counter = 0; counter < userM; counter++) */

    if (userN == 2)
    {
        determinantResult = (matrixA[0][0] * matrixA[1][1]) -
                                (matrixA[0][1] * matrixA[1][0]);
    } /* if (userN == 2) */

    else
    {
        for (counter = 0; counter < userN; counter++)
        {
            for (determinantARow = 0; determinantARow < userM;
                    determinantARow++)
            {
                determinantACheck = 0;
                for (determinantAColumn = 0; determinantAColumn < userM;
                        determinantAColumn++)
                {
                    if (determinantAColumn == counter)
                    {
                        determinantACheck++;
                    } /* if (determinantAColumn == counter) */
                    matrixDeterminantA[determinantARow][determinantAColumn] =
                            matrixA[determinantARow + 1][determinantACheck];
                    determinantACheck++;
                } /* for (determinantAColumn = 0; determinantAColumn < userM;
                        determinantAColumn++) */
            } /* for (determinantARow = 0; determinantARow < userM;
                    determinantARow++) */

            if(counter % 2 == 1)
            {
                sign = -1;
            } /* if(counter % 2 == 1) */
            else
            {
                sign = 1;
            } /* else */

            determinantResult += sign * matrixA[0][counter] *
                                    Determinant(matrixDeterminantA,userM);
        } /* for (counter = 0; counter < userN; counter++) */
    } /* else */


    for(counter = 0; counter < userM; counter++)
    {
        free(matrixDeterminantA[counter]);
    } /* for(counter = 0; counter < userM; counter++) */
    free(matrixDeterminantA);

    return determinantResult;
}

/**
 * Function Name: GaussElemination
 *
 * Function Description:
 *        This function solve the Inverse and I matrix from a matrix A
 *
 * Input: matrixA, matrixB, count, userN
 * Output: pivoted matrixA,B
 * Version: 1.0
 */
static int
GaussElemination(double** matrixA, double** matrixB, int userN)
{
    int row = 0;
    int column = 0;
    int count = 0;
    double division=0;

    // forward elemination
    for (count = 0; count < (userN - 1); count++)
    {
        Pivoting(matrixA, matrixB, count, userN);
        for (row = count; row < userN -1; row++)
        {
            division = matrixA[row+1][count] / matrixA[count][count];
            for (column = 0; column < userN; column++)
            {
                matrixA[row+1][column] -= division * matrixA[count][column];
                matrixB[row+1][column] -= division * matrixB[count][column];
            } /* for (column = 0; column < userN; column++) */
        } /* for (row = count; row < userN -1; row++) */
    } /* for (count = 0; count < (userN - 1); count++) */

    // backward elemination
    for (count = userN - 1; count > 0; count--)
    {
        for (row = count; row > 0; row--)
        {
            division = matrixA[row-1][count] / matrixA[count][count];
            for (column = 0; column < userN; column++)
            {
                matrixA[row-1][column] -= division * matrixA[count][column];
                matrixB[row-1][column] -= division * matrixB[count][column];
            } /* for (column = 0; column < userN; column++) */
        } /* for (row = count; row > 0; row--) */
    } /* for (count = userN - 1; count > 0 ;count--) */

    // adjust I
    for (row = 0; row < userN; row++)
    {
        division = matrixA[row][row];
        for (column = 0; column < userN ; column++)
        {
                matrixA[row][column] /= division;
                matrixB[row][column] /= division;
        } /* for (column = 0; column < userN ; column++) */
    } /* for (row = 0; row < userN; row++) */
    return 0;

}

/**
 * Function Name: Pivoting
 *
 * Function Description:
 *        This function sort matrix A ( Partial Pivoting )
 *
 * Input: matrixA, matrixB, count, userN
 * Output: pivoted matrixA,B
 * Version: 1.0
 */
static int
Pivoting(double** matrixA, double** matrixB, int count, int userN)
{
    int check;
    int where;

    double pivotTemporary;
    double big;

    where = count;
    big = Absolute(matrixA[count][count]);

    for(check = count ;check < (userN - 1);check++)
    {
        pivotTemporary = Absolute(matrixA[count][check]);
        if(pivotTemporary > big)
            {
                big = pivotTemporary;
                where = check;
            } /* if(pivotTemporary > big) */
    } /* for(check = count ;check < (userN - 1);check++) */

    if( where != count )
    {
        for(check = 0; check < userN; check++)
        {
            pivotTemporary = matrixA[where][check];
            matrixA[where][check] = matrixA[count][check];
            matrixA[count][check] = pivotTemporary;
            pivotTemporary = matrixB[where][check];
            matrixB[where][check] = matrixB[count][check];
            matrixB[count][check] = pivotTemporary;
        } /* for(check = 0 ; check < userN ; check++) */
    } /* if( where != count ) */

    return 0;
}

/**
 * Function Name: Absolute
 *
 * Function Description:
 *        This function return the abs(value) without math headerfile
 *
 * Input: value
 * Output: abs(value)
 * Version: 1.0
 */
static double
Absolute(double value)
{
    return (value > 0) ? value : -value;
}

닫기