行列計算プログラム
/* matrix.h*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define RANDOMIZE() srand(time(NULL))
#define RANDOM(x) (rand()%(x)) // 0..x-1 間の整数乱数
#define REALRANDOM(x) (((x)*rand())/RAND_MAX) // 0..x 間の実数乱数
/* 行列の最大サイズを設定 */
#define MAXSIZE 128
/* 実行対象の行列サイズを設定
* （デフォルト値は 3, 主プログラムで指定せよ ） */
int SIZE=3;
/* 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);
/********** 以下関数定義 **********/
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;
}
MATRIX array2mat(int a[SIZE][SIZE])
{
int i,j;
MATRIX b;
for(i=0;i<SIZE;i++)
for(j=0;j<SIZE;j++)
b.ent[i][j]=a[i][j];
return b;
}