ピボット選択がどのように行われてるか確かめたいので
関数pivotを呼び出す前と呼び出して式を入れ替えた後のakk,ak+1k,...,ankを調べたいのですが
どうすればいいのでしょうか？
int gauss(double *x, double *a, double *b, int n)
{
int i,j,k;
double tmp,p,sum;
/* step 1: 前進消去 */
for(k=0; k<n-1;k++){
printf("---- Step %d ----\n",k+1);
/* ピボットの選択 */
/*** 追加(2)：ピボット選択前，選択して式を入れ替えた後それぞれの
* a(k,k), a(k+1,k), ... , a(n-1,k)
* を表示して，ピボット選択がどのように行なわれたか調べる．
***/
/*--- ピボット選択前の位置 ---*/
printf("-- before --\n");
/* a(k,k), a(k+1,k), ... , a(n-1,k) を表示させる．*/
/*----------------------------*/
j = pivot(a,n,k); /* j: ピボットとして選ばれたa(j,k)の行番号 */
if(j == ERROR) {
return ERROR;
} else {
if(j != k) { /* ピボットにはa(k,k)ではなくa(j,k)が選ばれた．式の入れ替えが必要 */
/* Aのk行とj行の入れ替え．*/
for(i=0; i<n; i++){
tmp = a[n*k+i];
a[n*k+i] = a[n*j+i];
a[n*j+i] = tmp;
}
/* b[k] と b[j] の入れ替え */
tmp=b[j];
b[j]=b[k];
b[k]=tmp;
}
}
/*--- ピボット選択をし，式の入れ替えをした直後の位置 ---*/
printf("-- after --\n");
/* a(k,k), a(k+1,k), ... , a(n-1,k) を表示させる．
/*------------------------------------------------------*/
/* x[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]);
}
/*--- 追加(1)：第k段によって x[k] を消去した後の，各式の状態を表示する． ---*/
/* a(1,1),a(1,2),...,a(1,n),b(1); a(2,1),a(2,2),... a(2,n),b(2), ... */
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;
}
}
お礼
ありがとうございます。 一回のみならず二回も... 理解することができました。