/*--------------------------------------------------------------------------------
ハウスホルダー変換を用いたQR分解
QR decomposition using the Householder transformation
main.c
--------------------------------------------------------------------------------*/
#include "main.h"
int main(void)
{
double X[ROW_NUM][COL_NUM];
// XをQR分解によってQとRに分解
// decompose X into Q and R by QR decomposition
double Q[ROW_NUM][ROW_NUM]; // X = Q * R
double R[ROW_NUM][COL_NUM]; // X = Q * R
/*----------------------------------------------------------------------*/
X[0][0] = 1.0; X[0][1] = 2.0;
X[1][0] = 3.0; X[1][1] = 4.0;
X[2][0] = 5.0; X[2][1] = 6.0;
/*----------------------------------------------------------------------*/
// QR分解の実行
// perform QR decomposition
householder_QR_decomp( X, Q, R );
// QR分解の結果の表示
// show X, Q and R (Q and R are the result of QR decomposition)
show_XQR( X, Q, R );
return 0;
}
/*----------------------------------------------------------------------*/
void show_XQR( double X[ROW_NUM][COL_NUM], double Q[ROW_NUM][ROW_NUM], double R[ROW_NUM][COL_NUM] )
{
int i;
int j;
int k;
double s;
double t;
double temp;
// X[ROW_NUM][COL_NUM]
printf("Matrix X:\n");
for (i = 0; i < ROW_NUM; i++) {
for (j = 0; j < COL_NUM; j++) {
printf("%lf\t", X[i][j]);
// printf("%.16g\t", X[i][j]);
}
printf("\n");
}
printf("\n");
// Q[ROW_NUM][ROW_NUM]
printf("Matrix Q:\n");
for (i = 0; i < ROW_NUM; i++) {
for (j = 0; j < ROW_NUM; j++) {
printf("%lf\t", Q[i][j]);
// printf("%.16g\t", Q[i][j]);
}
printf("\n");
}
printf("\n");
// R[ROW_NUM][COL_NUM]
printf("Matrix R:\n");
for (i = 0; i < ROW_NUM; i++) {
for (j = 0; j < COL_NUM; j++) {
printf("%lf\t", R[i][j]);
// printf("%.16g\t", R[i][j]);
}
printf("\n");
}
printf("\n");
// Q * R
printf("Matrix Q * R:\n");
t = 0.0;
for (int i = 0; i < ROW_NUM; i++) {
for (int j = 0; j < COL_NUM; j++) {
s = 0.0;
for (k = 0; k < ROW_NUM; k++) {
s += Q[i][k] * R[k][j];
}
printf("%lf\t", s);
// printf("%.16g\t", s);
s -= X[i][j];
t += s * s;
}
printf("\n");
}
temp = t/(double)(ROW_NUM * COL_NUM);
printf("# 二乗平均平方根誤差: %g\n", sqrt(temp) );
printf("# root-mean-square error: %g\n", sqrt(temp) );
}