#include <stdio.h>
#include <math.h>

#include "nr.h"

int main()
{
	int i,j;
	double d;			// 行交換回数：1なら偶数回，-1なら奇数回
	double **input;		// 入力行列
	int *index;			// 部分ピボット選択の際の行交換を記録したもの
						// ピボット選択をすれば，LU分解されるのはinput行列そのものではなく，
						// input行列を変換したものである．
	double col[DIM],**y;

	/*  メモリ確保   */
	input = matrix(DIM,DIM);
	y = matrix(DIM,DIM);
	index = ivector(DIM);


	/*  入力行列を作成  */
	for (i=0;i<DIM;i++){
		for (j=0;j<DIM;j++){
			input[i][j] = 2;
		}
	}
	input[0][0] = 4;
	input[1][1] = 3;
	input[2][2] = 3;

	/*  LU分解   */
	ludcmp(input, DIM, index, &d);
	/* 　結果表示   */
	for (i=0;i<DIM;i++){
		printf("ピボット：%d\n",index[i]);	// ピボット表示，分からなければ無視してよい．
	}
	printf("行変換回数は奇数？偶数？：%lf\n",d);	// 偶数：1,奇数-1
	for (i=0;i<DIM;i++){
		for (j=0;j<DIM;j++){
			printf("%10lf",input[i][j]);
		}
		printf("\n");
	}
	printf("\n");
	/////////////////////////////////////////////

	/* 逆行列求める */
	for (j=0;j<DIM;j++){
		for (i=0;i<DIM;i++)	col[i] = 0.0;
		col[j] = 1.0;
		lubksb(input, DIM, index, col);

		for (i=0;i<DIM;i++)	y[i][j] = col[i];
	}
	/* 結果表示 */
	for (i=0;i<DIM;i++){
		for (j=0;j<DIM;j++){
			printf("%10lf",y[i][j]);
		}
		printf("\n");
	}
	printf("\n");

	/* 解を求める */
	col[0] = 0;
	col[1] = 3;
	col[2] = 4;
	lubksb(input, DIM, index, col);
	/* 結果表示 */
	for (i=0;i<DIM;i++){
			printf("%10lf",col[i]);
	}
	printf("\n");

	/*   メモリ開放   */
	free_matrix(input);
	free_matrix(y);
	free_ivector(index);

	return 0;
}