Gaussの消去法のプログラムなんですがこれを利用して、消去法による行

質問

Gaussの消去法のプログラムなんですがこれを利用して、消去法による行列式の計算プログラムをつくりたいのですが難しくてよくわかりません。。。
教えていただきたいです。
困ってるのでよろしくお願いします。
int gauss(double *x, double *a, double *b, int n)
{
int i,j,k,m;
double tmp,p,sum;
for(k=0; k<n-1;k++){
printf("---- Step %d ----\n",k+1);
printf("-- before --\n");
for( m = k;m < n;m++){
printf("%lf\n",a[n * m + k]);
}
printf("-- --\n");

j = pivot(a,n,k);
if(j == ERROR) {
return ERROR;
} else {
if(j != k) {
for(i=0; i<n; i++){
tmp = a[n*k+i];
a[n*k+i] = a[n*j+i];
a[n*j+i] = tmp;
}
tmp=b[j];
b[j]=b[k];
b[k]=tmp;
}
} printf("-- after --\n");
for( m = k;m < n;m++){
printf("%lf \n",a[n * m + k]);
}
for(i=k+1; i<n; i++){
p=a[n*i+k]/a[n*k+k];

for(j=0; j<n; j++){
a[n*i+j]=a[n*i+j]-p*a[n*k+j];
printf("a[%d %d]=%lf",i,j,a[n*i+j]);
}
b[i]=b[i]-p*b[k];
printf("b[%d]=%lf\n",i,b[i]);
}
printf("k=%d\n",k);
/*--------------------------------------------------------------------------*/
}
/* step 2: 後退代入 */
for(k=n-1; k>=0; k--){
if(fabs(a[n*k+k]) < EPS){
return(ERROR);
}
sum=0.0;
for(j=k+1; j<n; j++){
sum+=a[n*k+j]*x[j];
}
x[k]=(b[k] - sum)/a[k*n+k];
}

return 0;
}

int pivot(double *a, int n, int k)
{
int i,m;
double d;

/* ピボットの探索 */
m = k;
d = fabs(a[k*n+k]);
for(i=k+1; i<n; i++){
if(fabs(a[n*i+k]) > d){
m = i;
d = fabs(a[n*i+k]);
}
}

if(fabs(d) < EPS) {
return ERROR;
} else {
return m;
}
}

回答No.1

magicalpass · ベストアンサー

ガウスの消去法って一次連立方程式の解を求めるのでよかったっけ？

// コメントに「**」とある個所を追加・変更した
// 画面出力とかは省略

int gauss(double *x, double *a, double *b, int n)
{
int i,j,k,m;
double tmp,p,sum;

for(k=0; k<n;k++){ // **
j = pivot(a,n,k);
if(j == ERROR) {
return ERROR;
} else {
if(j != k) {
for(i=0; i<n; i++){
tmp = a[n*k+i];
a[n*k+i] = a[n*j+i];
a[n*j+i] = tmp;
}
tmp=b[j];
b[j]=b[k];
b[k]=tmp;
}
}
for(i=0; i<n; i++){ // **
if(i == k) { // **
p=a[n*k+k]; // **
for(j=0; j<n; j++){ // **
a[n*i+j]/=p; // **
} // **
b[i]/=p; // **
} else { // **
p=a[n*i+k]/a[n*k+k];
for(j=0; j<n; j++){
a[n*i+j]=a[n*i+j]-p*a[n*k+j];
}
b[i]=b[i]-p*b[k];
} // **
}
}
/* step 2: 後退代入 */
for(k=n-1; k>=0; k--){
if(fabs(a[n*k+k]) < EPS){
return(ERROR);
}
sum=0.0;
for(j=k+1; j<n; j++){
sum+=a[n*k+j]*x[j];
}
x[k]=(b[k] - sum)/a[k*n+k];
}
return 0;
}