/* matrix.h for Numerical Analysis C
by K.Shiota
2003.11.25
*/

#define RANDOM(x)	(rand()%(x))

/***** 扱うベクトル・行列の最大サイズと、型定義 *****/

#define MATRIXSIZE 128
typedef struct { int row; double ent[MATRIXSIZE]; } VECTOR;
typedef struct { int row,col; double ent[MATRIXSIZE][MATRIXSIZE]; } MATRIX;

/* VECTOR v;
   は v が縦ベクトルであることを表し、サイズは v.row, 第 i 要素は v.ent[i] 

   MATRIX a;
   は a が行列であることを表し、サイズは a.row×a.col, (i,j) 成分は a.ent[i][j]  */


/***** 関数の形式と機能 *****/

MATRIX vectomat(VECTOR a);
/* m 項縦ベクトルを m×1 行列に変換して返り値とする */

VECTOR vecadd(VECTOR a, VECTOR b);
/* ふたつのベクトル a, b の和 a+b を返り値とする */

VECTOR vecsub(VECTOR a, VECTOR b);
/* ふたつのベクトル a, b の差 a-b を返り値とする */

MATRIX matadd(MATRIX a, MATRIX b);
/* ふたつの行列 A, B の和 A+B を返り値とする */

MATRIX matsub(MATRIX a, MATRIX b);
/* ふたつの行列 A, B の差 A-B を返り値とする */

VECTOR matvecmul(MATRIX a, VECTOR b);
/* 行列 A と縦ベクトル b の積 Ab を返り値とする */

MATRIX matmul(MATRIX a, MATRIX b);
/* ふたつの行列 A, B の積 AB を返り値とする */

MATRIX matscalarmul(double x, MATRIX a);
/* 行列 A のスカラー x 倍 xA を返り値とする */

VECTOR vecscalarmul(double x, VECTOR a);
/* ベクトル a のスカラー x 倍 xa を返り値とする */

MATRIX zero_mat(int m, int n);
/* m×n のゼロ行列を返り値とする */

MATRIX unit_mat(int n);
/* n 次の単位行列を返り値とする */

MATRIX rot_mat(int n, int i, int j, double x);
/* 第 i,j 座標の角度 x の回転を表す行列を返り値とする */

double innerproduct(VECTOR a, VECTOR b);
/* ふたつのベクトル a, b の標準内積を返り値とする */

double vecnorm(VECTOR a);
/* ベクトル a の長さ（ユークリッドノルム）を返り値とする */

void print_vec(VECTOR a);
/* ベクトル a を出力する */

void print_mat(MATRIX a);
/* 行列 A を出力する */

VECTOR randomvec(int m);
/* ランダムな m 項縦ベクトルを返り値とする */


/***** 以下関数定義部 *****/

MATRIX vectomat(VECTOR a)
{
	int m=a.row,i;
	MATRIX c;
	
	c.row=m;
	c.col=1;
	for(i=1; i<=m; i++)
		c.ent[i][1]=a.ent[i];
	return c;
}

VECTOR vecadd(VECTOR a, VECTOR b)
{
	int m=a.row,i;
	VECTOR c;
	
	if(m!=b.row){
		printf("Invalid Matrix Size\n");
		exit(1);
	}
	c.row=m;
	for(i=1; i<=m; i++)
		c.ent[i]=a.ent[i]+b.ent[i];
	return c;
}

VECTOR vecsub(VECTOR a, VECTOR b)
{
	int m=a.row,i;
	VECTOR c;
	
	if(m!=b.row){
		printf("Invalid Matrix Size\n");
		exit(1);
	}
	c.row=m;
	for(i=1; i<=m; i++)
		c.ent[i]=a.ent[i]-b.ent[i];
	return c;
}

MATRIX matadd(MATRIX a, MATRIX b)
{
	int m=a.row,n=a.col,i,j;
	MATRIX c;
	
	if(m!=b.row || n!=b.col){
		printf("Invalid Matrix Size\n");
		exit(1);
	}
	c.row=m;
	c.col=n;
	for(i=1; i<=m; i++)
		for(j=1; j<=n; j++)
			c.ent[i][j]=a.ent[i][j]+b.ent[i][j];
	return c;
}

MATRIX matsub(MATRIX a, MATRIX b)
{
	int m=a.row,n=a.col,i,j;
	MATRIX c;
	
	if(m!=b.row || n!=b.col){
		printf("Invalid Matrix Size\n");
		exit(1);
	}
	c.row=m;
	c.col=n;
	for(i=1; i<=m; i++)
		for(j=1; j<=n; j++)
			c.ent[i][j]=a.ent[i][j]-b.ent[i][j];
	return c;
}

MATRIX matmul(MATRIX a, MATRIX b)
{
	int l=a.row,m=a.col,n=b.col,i,j,k;
	MATRIX c;
	double x;
	
	if(m!=b.row){
		printf("Invalid Matrix Size\n");
		exit(1);
	}
	c.row=l;
	c.col=n;
	for(i=1; i<=l; i++)
		for(j=1; j<=n; j++){
			x=0;
			for(k=1; k<=m; k++) x+=a.ent[i][k]*b.ent[k][j];
			c.ent[i][j]=x;
		}
	return c;
}

VECTOR matvecmul(MATRIX a, VECTOR b)
{
	int m=a.row,n=a.col,i,j,k;
	VECTOR c;
	double x;
	
	if(n!=b.row){
		printf("Invalid Matrix Size\n");
		exit(1);
	}
	c.row=m;
	for(i=1; i<=m; i++){
		x=0;
		for(k=1; k<=n; k++) x+=a.ent[i][k]*b.ent[k];
		c.ent[i]=x;
	}
	return c;
}

MATRIX matscalarmul(double x, MATRIX a)
{
	int m=a.row,n=a.col,i,j;
	MATRIX c;
	
	c.row=m;
	c.col=n;
	for(i=1; i<=m; i++)
		for(j=1; j<=n; j++)
			c.ent[i][j]=x*a.ent[i][j];
	return c;
}

VECTOR vecscalarmul(double x, VECTOR a)
{
	int m=a.row,i;
	VECTOR c;
	
	c.row=m;
	for(i=1; i<=m; i++)
		c.ent[i]=x*a.ent[i];
	return c;
}

MATRIX tr_mat(MATRIX a)
{
	int i,j,m=a.row,n=a.col;
	MATRIX b;
	
	b.row=n;
	b.col=m;
	for(i=1; i<=m; i++)
		for(j=1; j<=n; j++)
			b.ent[j][i]=a.ent[i][j];
	return b;
}

MATRIX zero_mat(int m, int n)
{
	int i,j;
	MATRIX a;
	a.row=m;
	a.col=n;
	for(i=1; i<=m; i++)
		for(j=1; j<=n; j++)
			a.ent[i][j]=0.0;
	return a;
}

MATRIX unit_mat(int n)
{
	int i;
	MATRIX a=zero_mat(n,n);
	for(i=1; i<=n; i++) a.ent[i][i]=1.0;
	return a;
}

MATRIX rot_mat(int n, int i, int j, double x)
{
	MATRIX a=unit_mat(n);
	a.ent[i][i]=cos(x);
	a.ent[i][j]=-sin(x);
	a.ent[j][i]=sin(x);
	a.ent[j][j]=cos(x);
	return a;
}

double innerproduct(VECTOR a, VECTOR b)
{
	int i,m=a.row;
	double x=0.0;
	
	if(m!=b.row){
		printf("Invalid Matrix Size\n");
		exit(1);
	}
	for(i=1; i<=m; i++)
		x+=a.ent[i]*b.ent[i];
	return x;
}

double vecnorm(VECTOR a)
{
	double x=innerproduct(a,a);
	return sqrt(x);
}

void print_vec(VECTOR a)
{
	int i,m=a.row;
	for(i=1; i<=m; i++)
		printf("%12.7f\n",a.ent[i]);
}

void print_mat(MATRIX a)
{
	int i,j,m=a.row,n=a.col;
	for(i=1; i<=m; i++){
		for(j=1; j<=n; j++) printf("%12.7f",a.ent[i][j]);
		printf("\n");
	}
}

VECTOR randomvec(int m)
{
	int i;
	VECTOR a;
	
	a.row=m;
	for(i=1; i<=m; i++)
		a.ent[i]=RANDOM(10);
	return a;
}