/* プログラム内で与えられた行列に対して */

/* 掃出法で逆行列と行列式の値を求める。 */

/* (画面表示し、同時にファイルに保存。) */

#include<stdio.h>

#define N 4 /* N次正方行列 */

#define M 8 /* M=N+N */

main(){

        float ma=1,b,c;

        int   i,j,k;

        float a[N][M]={{ 2.0, 5.0,-1.0, 4.5, 1.0, 0.0, 0.0, 0.0},

                       { 1.0, 3.0, 1.0, 2.0, 0.0, 1.0, 0.0, 0.0}

{ 3.0,-1.0,-2.0,10.0, 0.0, 0.0, 1.0, 0.0},

                       { 5.0,-2.0, 8.0,-4.0, 0.0, 0.0, 0.0, 1.0}};

        FILE *fp;

        fp=fopen("E:/Cpp/40102A.dat","w");

  /* 与えられた行列をファイルに保存 */

        fprintf(fp,"\n a[][]=\n");

        for(j=0;j<N;j++){

                for(k=0;k<N;k++){

                        fprintf(fp,"7.2f",a[j][k]);

                }

                fprintf(fp,"\n");

        }

  /* step i :i行i列は(i,i)要素以外は0にする操作 */

        for(i=0;i<N;i++){ 

                for(j=0;j<N;j++){

                        if(i!=j){

                                b=a[j][i]/a[i][i];

                                for(k=i;k<M;k++){

                                        a[j][k]=a[j][k]-a[i][k]*b;

                                }

                        }

                }

                ma=ma*a[i][i];

                printf("STEP %d\n",i+1);

                fprintf(fp,"\nSTEP %d\n",i+1);

                for(j=0;j<N;j++){

                        for(k=0;k<M;k++){

                                printf("%9.4f",a[j][k]);

                                fprintf(fp,"%9.4f",a[j][k]);

                        }

                        printf("\n");

fprintf(fp,"\n");

                }

        }

  /* end */

  /* 左半分を単位行列にする操作 */

        for(j=0;j<N;j++){

                c=a[j][j];

                for(k=j;k<M;k++){

                        a[j][k]=a[j][k]/c;

                }

        }

        printf("STEP %d\n",N+1);

        fprintf(fp,"\nSTEP %d\n",N+1);

        for(j=0;j<N;j++){

                for(k=0;k<M;k++){

                        printf("%9.4f",a[j][k]);

                        fprintf(fp,"%9.4f",a[j][k]);

                }

                printf("\n");

                fprintf(fp,"\n");

        }

  /* end */

  /* 逆行列の表示 */

        printf("\nInverse of a[][]=\n");

        fprintf(fp,"\nInverse of a[][]=\n");

        for(j=0;j<N;j++){

                for(k=0;k<M;k++){

                        printf("%9.4f",a[j][k]);

                        fprintf(fp,"%9.4f",a[j][k]);

                }

                printf("\n");

                fprintf(fp,"\n");

        }

  /* end */

  /* 行列式の値を表示 */

        printf("\nDeterminant of a[][] =%9.3f         [End]\n\n",ma);

        fprintf(fp,"\nDeterminant of a[][] =%9.3f         [End]\n\n",ma);

  /* end */

        fclose(fp);

        return(1);

}