/* 출처는 모대학에서 리포트용으로 작성된 소스인 것같습니다.
알고리즘에 관한한 저도 문외한이니 저에게 질문을 하지는 말아주세요 ^^ */


#include <stdio.h>
#include <math.h>  // 수학적 계산이 있기때문에 사용. 이 프로그램에서는 사용 안해도 됨.

void gauss(int);    // 가우스 소거법을 실행하는 서브루틴을 설정
double A[50][51];  // 입력행렬을 광역변수로 선언

void main(void)
{
    int i, j, n;

    // 정사각행렬 A의 크기를 입력
    printf("Enter the dimension of matrix : ");
    scanf("%d", &n);
   
    // 행렬 A의 원소를 입력
    printf("\nEnter the elements of matrix A: \n");
    for(i=0;i<n;i++)
    {
        for(j=0;j<n;j++)
    {
        printf("A[%d][%d] = ",i+1, j+1);
               scanf("%lf", &A[i][j]);
        }
          printf("\n");
    }

    // 벡터 b를 입력, 편의상 b[i]=A[i][n]으로 놓는다.
    printf("Enter the elements of b: \n");
    for(i=0;i<n;i++)
    {
        printf("b[%d]= ",i+1);
        scanf("%lf",&A[i][n]);
    }

    gauss(n);  // 가우스 소거법을 실행 할 수 있는 서브루틴을 호출

}

// 가우스 소거법
void gauss(int n)
{
    int i, j, k, h;
    double sum, m;
   
    for(j=0; j<n; j++)
    {
    // 행렬의 대각성분의 원소가 0이 아니도록 만듬.
    //if(A[j][j]==0)
    //{
    //    for(i=j+1; i<n; i++)
    //    {
    //        if(A[i][j]!=0)
    //        {    // i행과 j행을 바꾸어 준다.
    //            for(k=j; k<n+1; k++)
    //            {
    //                m=A[j][k];  
    //                A[j][k]=A[i][k];
    //                A[i][k]=m;
    //            }
    //            break;
    //        }
    //    }
    //}
    // 피벗 Gauss 소거법
        // 열벡터에서 피벗행을 찾음
    m=fabs(A[j][j]);
    for(i=j+1; i<n; i++)
    {
        if(fabs(A[i][j]) > m)
        {
            m=fabs(A[i][j]);    h=i;
        }
    }
    // h행과 j행을 바꾸어 준다.
    for(k=j; k<n+1; k++)
    {
        m=A[j][k];  
        A[j][k]=A[h][k];
        A[h][k]=m;
    }
    // 피벗 Gauss 소거법 끝
    // 대각성분의 원소중 0이 있으면 실행 중단.
    if(A[j][j]==0)
    {
        printf("\n-----------------------------\n");
        printf("     Matrix is singular!        ");   // 정칙행렬이 아니면 출력.
        printf("\n-----------------------------\n\n");
        return;  // 실행 종료
        }
        // 소거 시작
        for(i=j+1; i<n; i++)
        for(k=j+1; k<n+1; k++)
            A[i][k]=A[i][k]-(A[i][j]/A[j][j])*A[j][k];
    }   

    // 역 대입법
    for(i=n-1; i>=0; i--)
    {
        for(k=i+1; k<n; k++)
        {
            A[i][n]=A[i][n]-A[i][k]*A[k][n];
        }
        A[i][n]=A[i][n]/A[i][i];
    }

    printf("\nSolution : \n");
    for(i=0; i<n; i++)
    printf("x[%d] = %lf \n", i+1, A[i][n]);

    return;

}