ludecomp.c

ludecomp.c

#include<stdio.h>

main(){
	long i,j,N=4,NRHS=1,IPIV[N],LDA=N,LDB=N,INFO;
	double A[N*N],b[N];
	/* 変数の説明 N:変数の数 NRHS:右辺の列数 IPIV[N]:Pivot選択の情報
	LDA:左辺の行数 LDB:右辺の行数 INFO:ルーチンが成功したかどうかを返す */
	
	/* 配列の引数は列に対して連続になるようにする(例:1行1列の次は2行1列) */
	
	A[0]=1;A[4]=-1;A[8]=0;A[12]=2;
	A[1]=0;A[5]=1;A[9]=2;A[13]=3;
	A[2]=2;A[6]=0;A[10]=-1;A[14]=3;
	A[3]=2;A[7]=-3;A[11]=1;A[15]=4;
	
	b[0]=0;
	b[1]=2;
	b[2]=1;
	b[3]=1;
	
	/* dgetrf_で行列をLU分解する。AにはLU分解後の行列が入る
	(Lの対角成分は1なのでそれ以外をUの下三角部分に格納したもの)*/
	
	dgetrf_(&N,&N,A,&N,IPIV,&INFO);
	
	/*
	for(i=0;i<N;i++){
		for(j=0;j<N;j++)
			printf("%lf ",A[i+j*N]);
		printf("\n");
	}
	*/
	
	/* dgetrs_によって、LU分解後の行列Aを用いて行列方程式を解く*/
	
	dgetrs_("N",&N,&NRHS,A,&N,IPIV,b,&LDB,&INFO);
	
	for(i=0;i<N;i++)
		printf("%lf\n",b[i]);
}