householder_QR_decomp.c

/*--------------------------------------------------------------------------------

ハウスホルダー変換を用いたQR分解

QR decomposition using the Householder transformation

householder_QR_decomp.c

--------------------------------------------------------------------------------*/

#include "lorenz_lyap.h"

void householder_QR_decomp( double X[DIM][DIM], double Q[DIM][DIM], double R[DIM][DIM] )

{

int i;

int j;

int k;

double alpha;

double norm;

double v[DIM][DIM];  // Householder vector, ハウスホルダー変換ベクトル

double u[DIM][DIM];  // temporary matrix

// initialize Q as an identity matrix

for (i = 0; i < DIM; i++) {

for (j = 0; j < DIM; j++) {

Q[i][j] = (i == j) ? 1.0 : 0.0;

}

}

for (i = 0; i < DIM; i++) {

for (j = 0; j < DIM; j++) {

R[i][j] = X[i][j];

}

}

// perform QR decomposition

// QR分解を実行

for (k = 0; k < DIM; k++) {

// compute the norm of the k-th column of R

// Rのk列のノルムを計算

norm = 0.0;

for (i = k; i < DIM; 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 < DIM; i++) {

u[i][k] = R[i][k];

}

// compute the Householder vector v

// ハウスホルダー変換ベクトルvを計算

norm = 0.0;

for (i = k; i < DIM; i++) {

norm += u[i][k] * u[i][k];

}

norm = sqrt(norm);

for (i = k; i < DIM; i++) {

v[i][k] = u[i][k] / norm;

}

// update R and Q using the Householder transformation

// ハウスホルダー変換を使ってRとQを更新

for (j = k; j < DIM; j++) {

alpha = 0.0;

for (i = k; i < DIM; i++) {

alpha += v[i][k] * R[i][j];

}

alpha *= 2.0;

for (i = k; i < DIM; i++) {

R[i][j] -= alpha * v[i][k];

}

}

for (j = 0; j < DIM; j++) {

alpha = 0.0;

for (i = k; i < DIM; i++) {

alpha += v[i][k] * Q[j][i];

}

alpha *= 2.0;

for (i = k; i < DIM; i++) {

Q[j][i] -= alpha * v[i][k];

}

}

}

}