数理情報学４ pivot 枢軸選択つきガウスの消去法

#include

#include

#define N 4
main(){

	double a[N][N],b[N],ell,a_max;

	int pivot[N],i,j,k,m,m_opt,k_step_old;

    

	//初期設定

	a[0][0]= 1.0; a[0][1]= 1.0; a[0][2]= 1.0; a[0][3]= 1.0;

	a[1][0]= 1.0; a[1][1]= 1.0; a[1][2]= 2.0; a[1][3]= 2.0;

	a[2][0]= 1.0; a[2][1]= 2.0; a[2][2]= -1.0; a[2][3]= 2.0;

	a[3][0]= 1.0; a[3][1]= 2.0; a[3][2]= 2.0; a[3][3]= -1.0;
	b[0]= 1.0; b[1]= 2.0; b[2]= 1.0; b[3]= 2.0; 
    //仮のピボット選択

    for(i=0;ia_max){

                a_max = fabs(a[pivot[m]][k]);

                m_opt = m;

            }

        }

        /*m_opt != kのとき

          「k段のピポット」と「m_opt段のピポット」を交換*/

        if(m_opt != k){

            k_step_old = pivot[k];

            pivot[k] = pivot[m_opt];

            pivot[m_opt] = k_step_old;

        }

        //前進消去

        for(m=k+1;m=0;m--){

        for(j=N-1;j>m;j--){

            b[pivot[m]] =  b[pivot[m]] -a[pivot[m]][j]*b[pivot[j]];

        }

        b[pivot[m]] =  b[pivot[m]]/a[pivot[m]][m];

    }

    

    printf("%f\n",b[pivot[0]]);

    printf("%f\n",b[pivot[1]]);

    printf("%f\n",b[pivot[2]]);

    printf("%f\n",b[pivot[3]]);

}