プログラマ専用SNS ミクプラ

コード:

#include <stdio.h>
#include <math.h>
#define N 4
#define N1 8
#define EPSILON 1.0E-5
#define TRUE 1
#define FALSE 0

float a[N][N1]={{1.0,0.75,0.0,-0.5,1.0,0.0,0.0,0.0},
{0.0,-0.25,0.0,0.5,0.0,1.0,0.0,0.0},
{0.0,0.75,0.0,-0.5,0.0,0.0,1.0,0.0},
{-2.0,-1.5,-1.0,1.0,0.0,0.0,0.0,1.0}};
float det;
int sw;
void main()
{
int i,j;
void sweep();
printf("\n\n===== 元の行列 ====\n\n");
for(i=0;i<N;i++)
{
printf("%d:",i+1);
for(j=0;j<N;j++) 
printf(" %10.4f",a[i][j]);
printf("\n");
}
printf("\n\n");
sw=TRUE;

sweep(&sw);

if(sw==TRUE)
{
printf("\n\n==== 逆行列 ====\n\n");
for(i=0;i<N;i++)
{
printf(" %d;",i+1);
for(j=N;i<N1;j++) printf(" %d10.4f ",a[i][j]);
printf("\n");
}
printf("\n 行列式=%8.4f \n",det);
}
else
{
printf("\n逆行列は求まりませんでした!\n");
}
}

void swap(wk1,wk2)
float *wk1,*wk2;
{
float w;
w=*wk1; *wk1=*wk2;*wk2=w;
}

void sweep(sw)
int *sw;
{
int i,j,k,ik;
float ak,aik,p;

det=1.0;
p=1.0;
for(k=0;k<N;k++)
{
ak=a[k][k];
if(fabs(ak)<=EPSILON)
{
ik=k+1;
while((ik<N) && (fabs(a[ik][k])<EPSILON)) ik++;
if(ik<N)
{
for(j=k;j<N1;j++){ swap(&a[k][j],&a[ik][j]);}
ak=a[k][k];
p=-p;
}
if(ik>=N)
{
printf(" ピボットが零です ! \n");
*sw=FALSE;
goto owari;
}
}
det=det*ak;
for(j=k;j<N1;j++) a[k][j]=a[k][j]/ak;
for(i=0;i<N;i++)
{
if(i!=k)
{
aik=a[i][k];
for(j=k;j<N1;j++) a[i][j]=a[i][j]-aik*a[k][j];
}
}
}
det=p*det;
owari:;
}