問題：n×n行列A、n×m行列Bに関して、次の連立方程式を、Gauss-Jordan法を用いて解け。
ただし、完全ピボット選択ではなく、部分ピボット選択を用いよ。
また、解があっているかをチェックする機能もつけよ。

　　　　　　　　　　　　　　　　　　　　AX＝Ｂ

/*2kadai1.c*/
#include <stdio.h>
#include<math.h>
#define EPS 1e-12
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}
double temp;

/*領域の確保*/
typedef double **mat;
mat matrix(m,n)
int m,n;
{
int i; mat x;

x=(double **) malloc(m*sizeof(double *));
if(!x) printf("allocation failure 1 in matrix()");
x[0] = (double *) malloc(m*n*sizeof(double));
if(!x[0]) printf("allocation failure 2 in matrix()");
for(i=1;i<m;i++) x[i] = x[i-1] + n;
return(x);
}
void freemat(x)
mat x;
{
free(x[0]);
free(x);
}

/*ピボットの定義*/
void pivot(a,n,k,l)
mat a;int n,k,*l;
{
int i; double big;

big=fabs(a[k][k]);
for(i=k+1;i<n;i++)
if(fabs(a[i][k])>big){
*l=i;
big=fabs(a[i][k]);
}
if(fabs(big)<EPS)printf("error in pivot()");
}

/*Gauss Jordanの定義*/
void gauss(a,b,n,m)
int n,m; mat a,b;
{
int i,j,k,h,piv; double p,q;

for(i=0;i<n;i++){
piv = i;
pivot(a,n,i,&piv);
if(i!=piv){
for(j=i;j<n;j++){SWAP(a[i][j],a[piv][j]);}
for(j=i;j<m;j++){SWAP(b[i][j],b[piv][j]);}}
p = a[i][i];
for(j=0;j<n;j++){a[i][j] = a[i][j]/p;}
for(j=0;j<m;j++){b[i][j] = b[i][j]/p;}
for(k=0;k<n;k++) {if(k!=i){q=a[k][i];
for(h=0;h<n;h++){a[k][h] -= q*a[i][h];}
for(h=0;h<m;h++){b[k][h] -= q*b[i][h];}}

}
}
}

main()
{
int n,m,i,j,k;
mat A,B,C,D,P;
double num,d,p;

/*行列A(行数n,列数n),行列B(行数n,列数m)のサイズを決める*/
printf("行列A(行数n,列数n),行列B(行数n,列数m)のサイズを決めます。\n");
printf("nを入力してください。");
scanf("%d",&n);
printf("mを入力してください。");
scanf("%d",&m);
printf("行列Aは%d×%d行列、行列Bは%d×%d行列です。\n",n,n,n,m);

A=matrix(n,n);
B=matrix(n,m);
C=matrix(n,n);
D=matrix(n,m);
P=matrix(n,m);

/*行列Aの各成分を入力*/
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("行列Aの%d行%d列の成分は%lfです。\n",i+1,j+1,A[i][j]);
}
}
/*行列Bの各成分を入力*/
for(i=0;i<n;i++){
for(j=0;j<m;j++){
printf("行列Bの%d行%d列の成分は？・・・",i+1,j+1);
scanf("%lf",&B[i][j]);
printf("行列Bの%d行%d列の成分は%lfです。\n",i+1,j+1,B[i][j]);
}
}

/* A,Bをコピーしておく(検算用 C,Dとする) */
for(i=0;i<n;i++){
for(j=0;j<n;j++){
C[i][j]=A[i][j];}
for(j=0;j<m;j++){
D[i][j]=B[i][j];}
}

/*Gauss-Jordan法*/
gauss(A,B,n,m);

/* 解の出力 */
for(j=0; j<n; j++){
for(k=0; k<m; k++){
printf("解Xの[%d]行[%d]列成分 = %lf\n", j+1,k+1,B[j][k]);
}
}

/* 検算 */
for(i=0; i<n; i++){
for(j=0; j<m; j++){
for(k=0; k<n; k++){
P[i][j] += C[i][k]*B[k][j];
}
printf("Bの%d行%d列成分は%lfです。\n",i+1,j+1,P[i][j]);

/* 結果と一致しなければエラー終了 */
if(fabs(P[i][j]-D[i][j])>EPS)
{printf("%d行%d列・・・Matching Error!\n",i+1,j+1);}
else{printf("%d行%d列・・・OK!\n",i+1,j+1);}
}
}

/*解放*/
freemat(A);
freemat(B);
freemat(C);
freemat(D);
freemat(P);
}