/*--------------------------------------------------------------------------------
ハウスホルダー変換を用いたQR分解
QR decomposition using the Householder transformation
householder_QR_decomp.c
--------------------------------------------------------------------------------*/
#include "main.h"
void householder_QR_decomp( 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 alpha;
double norm;
double v[ROW_NUM][COL_NUM]; // Householder vector, ハウスホルダー変換ベクトル
double u[ROW_NUM][COL_NUM]; // temporary matrix
// initialize Q as an identity matrix
for (i = 0; i < ROW_NUM; i++) {
for (j = 0; j < ROW_NUM; j++) {
Q[i][j] = (i == j) ? 1.0 : 0.0;
}
}
for (i = 0; i < ROW_NUM; i++) {
for (j = 0; j < COL_NUM; j++) {
R[i][j] = X[i][j];
}
}
// perform QR decomposition
// QR分解を実行
for (k = 0; k < COL_NUM; k++) {
// compute the norm of the k-th column of R
// Rのk列のノルムを計算
norm = 0.0;
for (i = k; i < ROW_NUM; i++) {
norm += R[i][k] * R[i][k];
}
norm = sqrt( norm );
// compute alpha and update the k-th column of R
// alphaを計算し、Rのk列を更新
alpha = (R[k][k] > 0) ? -norm : norm;
u[k][k] = R[k][k] - alpha;
for (i = k+1; i < ROW_NUM; i++) {
u[i][k] = R[i][k];
}
// compute the Householder vector v
// ハウスホルダー変換ベクトルvを計算
norm = 0.0;
for (i = k; i < ROW_NUM; i++) {
norm += u[i][k] * u[i][k];
}
norm = sqrt(norm);
for (i = k; i < ROW_NUM; i++) {
v[i][k] = u[i][k] / norm;
}
// update R and Q using the Householder transformation
// ハウスホルダー変換を使ってRとQを更新
for (j = k; j < COL_NUM; j++) {
alpha = 0.0;
for (i = k; i < ROW_NUM; i++) {
alpha += v[i][k] * R[i][j];
}
alpha *= 2.0;
for (i = k; i < ROW_NUM; i++) {
R[i][j] -= alpha * v[i][k];
}
}
for (j = 0; j < ROW_NUM; j++) {
alpha = 0.0;
for (i = k; i < ROW_NUM; i++) {
alpha += v[i][k] * Q[j][i];
}
alpha *= 2.0;
for (i = k; i < ROW_NUM; i++) {
Q[j][i] -= alpha * v[i][k];
}
}
}
}