/* sqmatrix.h
K.Shiota
2006.1
*/

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define RANDOM(x)	(rand()%(x))		// 0..x-1 間の整数乱数
#define REALRANDOM(x)	(((x)*rand())/RAND_MAX)	// 0..x 間の実数乱数

/* 行列の最大サイズを設定 */
#define MAXSIZE 128

/* 実行対象の行列サイズを設定
   （デフォルト値は 2, 主プログラムで指定せよ ） */
int SIZE=2; 

/* MATRIX, VECTOR, ROWVECTOR という型を定義 */
typedef struct { double ent[MAXSIZE]; } VECTOR;
typedef struct { double ent[MAXSIZE]; } ROWVECTOR;
typedef struct { double ent[MAXSIZE][MAXSIZE]; } MATRIX;

/********************************************************

VECTOR v;
と宣言すると v.ent[i] が列ベクトル v の 第i成分を表す。 

ROWVECTOR v;
と宣言すると v.ent[i] が行ベクトル v の 第i成分を表す。 

MATRIX a;
と宣言すると a.ent[i][j] が行列 a の (i,j)-成分を表す。 

*********************************************************/

/* 列ベクトルを表示する関数 */
void vecprint(VECTOR v);

/* 行ベクトルを表示する関数 */
void rowvecprint(ROWVECTOR v);

/* 行列を表示する関数 */
void matprint(MATRIX a);

/* 列ベクトル v のスカラー倍 w = c v を返す関数 */
VECTOR vecscmul(double c, VECTOR v);

/* 行ベクトル v のスカラー倍 w = c v を返す関数 */
ROWVECTOR rowvecscmul(double c, ROWVECTOR v);

/* 行列 A のスカラー倍 B = c A を返す関数 */
MATRIX matscmul(double c, MATRIX a);

/* ふたつの列ベクトル u, v の和 w = u + v を返す関数 */
VECTOR vecadd(VECTOR u, VECTOR v);

/* ふたつの列ベクトル u, v の差 w = u - v を返す関数 */
VECTOR vecsub(VECTOR u, VECTOR v);

/* ふたつの行ベクトル u, v の和 w = u + v を返す関数 */
ROWVECTOR rowvecadd(ROWVECTOR u, ROWVECTOR v);

/* ふたつの行ベクトル u, v の差 w = u - v を返す関数 */
ROWVECTOR rowvecsub(ROWVECTOR u, ROWVECTOR v);

/* ふたつの行列 A, B の和 C = A + B を返す関数 */
MATRIX matadd(MATRIX a, MATRIX b);

/* ふたつの行列 A, B の差 C = A - B を返す関数 */
MATRIX matsub(MATRIX a, MATRIX b);

/* 行ベクトル u と列ベクトルv の積 w = u v を返す関数 */
double rowvecvecmul(ROWVECTOR u, VECTOR v);

/* 列ベクトル u と行ベクトル v の積 A = u v を返す関数 */
MATRIX vecrowvecmul(VECTOR u, ROWVECTOR v);

/* 行列 A と列ベクトル v の積 w = A v を返す関数 */
VECTOR matvecmul(MATRIX a, VECTOR v);

/* 行ベクトル v と行列 A の積 w = v A を返す関数 */
ROWVECTOR rowvecmatmul(ROWVECTOR b, MATRIX a);

/* 行列 A, B の積 C = A B を返す関数 */
MATRIX matmul(MATRIX a, MATRIX b);

/* ゼロベクトル（列ベクトル）を返す関数 */
VECTOR zerovec();

/* ゼロベクトル（行ベクトル）を返す関数 */
ROWVECTOR zerorowvec();

/* ゼロ行列を返す関数 */
MATRIX zeromat();

/* 単位行列を返す関数 */
MATRIX unitmat();

/* 第 i,j 座標の角度 x の回転を表す行列を返す関数 */
MATRIX rotmat(int i, int j, double x);

/* 列ベクトル v の転置ベクトルを返す関数 */
ROWVECTOR transvec(VECTOR v);

/* 行ベクトル v の転置ベクトルを返す関数 */
VECTOR transrowvec(ROWVECTOR v);

/* 行列 A の転置行列を返す関数 */
MATRIX transmat(MATRIX a);

/* ふたつのベクトル u, v の標準内積を返す関数 */
double innerprod(VECTOR u, VECTOR v);

/* 行ベクトル v の長さ（ユークリッドノルム）を返す関数 */
double vecnorm(VECTOR v);

/* 列ベクトル v の長さ（ユークリッドノルム）を返す関数 */
double rowvecnorm(ROWVECTOR v);

/* ランダムな列ベクトルを返す関数（成分の最大絶対値 < x） */
VECTOR randomvec(double x);

/* ランダムな整数成分列ベクトルを返す関数（成分の最大絶対値 < x） */
VECTOR randomintvec(int x);

/* ランダムな行ベクトルを返す関数（成分の最大絶対値 < x） */
ROWVECTOR randomrowvec(double x);

/* ランダムな整数成分行ベクトルを返す関数（成分の最大絶対値 < x） */
ROWVECTOR randomintrowvec(int x);

/* ランダムな行列を返す関数（成分の最大絶対値 < x）*/
MATRIX randommat(double x);

/* ランダムな整数成分行列を返す関数（成分の最大絶対値 < x）*/
MATRIX randomintmat(int x);

/* 行列 A の LU 分解を計算する関数
   ２つの行列を返すために、L と U はポインタで受け渡しする。 */
void ludecomp(MATRIX a, MATRIX *l, MATRIX *u);

/* LU 分解法の前進代入 : L y = b の解 y を返す関数 */
VECTOR lustep2(MATRIX l, VECTOR b);

/* LU 分解法の後退代入 : U x = y の解 x を求める関数 */
VECTOR lustep3(MATRIX u, VECTOR y);


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

void vecprint(VECTOR v)
{
	int i;
	for(i=0;i<SIZE;i++) printf("%12.7lf\n",v.ent[i]);
}

void rowvecprint(ROWVECTOR v)
{
	int i;
	for(i=0;i<SIZE;i++) printf("%12.7lf  ",v.ent[i]);
	printf("\n");
}

void matprint(MATRIX a)
{
	int i,j;
	for(i=0;i<SIZE;i++){
		for(j=0;j<SIZE;j++) printf("%12.7lf",a.ent[i][j]);
		printf("\n");
	}
}

VECTOR vecscmul(double c, VECTOR v)
{
	VECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=c*v.ent[i];
	return w;
}

ROWVECTOR rowvecscmul(double c, ROWVECTOR v)
{
	ROWVECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=c*v.ent[i];
	return w;
}

MATRIX matscmul(double c, MATRIX a)
{
	MATRIX b;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			b.ent[i][j]=c*a.ent[i][j];
	return b;
}

VECTOR vecadd(VECTOR u, VECTOR v)
{
	VECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=u.ent[i]+v.ent[i];
	return w;
}

VECTOR vecsub(VECTOR u, VECTOR v)
{
	VECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=u.ent[i]-v.ent[i];
	return w;
}

ROWVECTOR rowvecadd(ROWVECTOR u, ROWVECTOR v)
{
	ROWVECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=u.ent[i]+v.ent[i];
	return w;
}

ROWVECTOR rowvecsub(ROWVECTOR u, ROWVECTOR v)
{
	ROWVECTOR w;
	int i;
	for(i=0;i<SIZE;i++) w.ent[i]=u.ent[i]-v.ent[i];
	return w;
}

MATRIX matadd(MATRIX a, MATRIX b)
{
	MATRIX c;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			c.ent[i][j]=a.ent[i][j]+b.ent[i][j];
	return c;
}

MATRIX matsub(MATRIX a, MATRIX b)
{
	MATRIX c;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			c.ent[i][j]=a.ent[i][j]-b.ent[i][j];
	return c;
}

double rowvecvecmul(ROWVECTOR u, VECTOR v)
{
	double x=0.0;
	int k;
	for(k=0;k<SIZE;k++) x+=u.ent[k]*v.ent[k];
	return x;
}

MATRIX vecrowvecmul(VECTOR u, ROWVECTOR v)
{
	MATRIX a;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			a.ent[i][j]=u.ent[i]*v.ent[j];
	return a;
}

VECTOR matvecmul(MATRIX a, VECTOR v)
{
	VECTOR w;
	int i,k;
	double x;
	for(i=0;i<SIZE;i++){
		x=0.0;
		for(k=0;k<SIZE;k++) x+=a.ent[i][k]*v.ent[k];
		w.ent[i]=x;
	}
	return w;
}

ROWVECTOR rowmatvecmul(ROWVECTOR v, MATRIX a)
{
	ROWVECTOR w;
	int j,k;
	double x;
	for(j=0;j<SIZE;j++){
		x=0.0;
		for(k=0;k<SIZE;k++) x+=v.ent[k]*a.ent[k][j];
		w.ent[j]=x;
	}
	return w;
}

MATRIX matmul(MATRIX a, MATRIX b)
{
	MATRIX c;
	int i,j,k;
	double x;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++){
			x=0.0;
			for(k=0;k<SIZE;k++) x+=a.ent[i][k]*b.ent[k][j];
			c.ent[i][j]=x;
		}
	return c;
}

VECTOR zerovec()
{
	VECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=0;
	return v;
}

ROWVECTOR zerorowvec()
{
	ROWVECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=0;
	return v;
}

MATRIX zeromat()
{
	MATRIX a;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			a.ent[i][j]=0.0;
	return a;
}

MATRIX unitmat()
{
	int i;
	MATRIX a=zeromat();
	for(i=0;i<SIZE;i++) a.ent[i][i]=1.0;
	return a;
}

MATRIX rotmat(int i, int j, double x)
{
	MATRIX a=unitmat();
	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;
}

ROWVECTOR transvec(VECTOR v)
{
	ROWVECTOR w;
	int k;
	for(k=0;k<SIZE;k++) w.ent[k]=v.ent[k];
	return w;
}

VECTOR transrowvec(ROWVECTOR v)
{
	VECTOR w;
	int k;
	for(k=0;k<SIZE;k++) w.ent[k]=v.ent[k];
	return w;
}

MATRIX transmat(MATRIX a)
{
	int i,j;
	MATRIX b;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			b.ent[j][i]=a.ent[i][j];
	return b;
}

double innerprod(VECTOR u, VECTOR v)
{
	double x=0.0;
	int k;
	for(k=0;k<SIZE;k++) x+=u.ent[k]*v.ent[k];
	return x;
}

double vecnorm(VECTOR v)
{
	return sqrt(innerprod(v,v));
}

double rowvecnorm(ROWVECTOR v)
{
	return vecnorm(transrowvec(v));
}

VECTOR randomvec(double x)
{
	VECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=2*REALRANDOM(x)-x;
	return v;
}

VECTOR randomintvec(int x)
{
	VECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=RANDOM(2*x-1)-x+1;
	return v;
}

ROWVECTOR randomrowvec(double x)
{
	ROWVECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=2*REALRANDOM(x)-x;
	return v;
}

ROWVECTOR randomintrowvec(int x)
{
	ROWVECTOR v;
	int i;
	for(i=0;i<SIZE;i++) v.ent[i]=RANDOM(2*x-1)-x+1;
	return v;
}

MATRIX randommat(double x)
{
	MATRIX a;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			a.ent[i][j]=2*REALRANDOM(x)-x;
	return a;
}

MATRIX randomintmat(int x)
{
	MATRIX a;
	int i,j;
	for(i=0;i<SIZE;i++)
		for(j=0;j<SIZE;j++)
			a.ent[i][j]=RANDOM(2*x-1)-x+1;
	return a;
}

/********** 以下は２次行列限定 **********/

/* 正方行列 A の LU 分解を計算する関数 */
void lu(MATRIX a, MATRIX *l, MATRIX *u)
{
	l->ent[0][0]=a.ent[0][0];
	u->ent[0][1]=a.ent[0][1]/l->ent[0][0];
	l->ent[1][0]=a.ent[1][0];
	l->ent[1][1]=a.ent[1][1]-l->ent[1][0]*u->ent[0][1];

	l->ent[0][1]=0.0;
	u->ent[0][0]=1.0;
	u->ent[1][0]=0.0;
	u->ent[1][1]=1.0;
}
/* （注：pivot 選択は使っていない） */

/* LU 分解法の前進代入 : L y = b の解 y を返す関数 */
VECTOR step2(MATRIX l, VECTOR b)
{
	VECTOR y;
	y.ent[0]=b.ent[0]/l.ent[0][0];
	y.ent[1]=(b.ent[1]-l.ent[1][0]*y.ent[0])/l.ent[1][1];
	return y;
}

/* LU 分解法の後退代入 : U x = y の解 x を求める関数 */
VECTOR step3(MATRIX u, VECTOR y)
{
	VECTOR x;
	x.ent[1]=y.ent[1];
	x.ent[0]=y.ent[0]-u.ent[0][1]*x.ent[1];
	return x;
}