/*-----------------------------------------------------
グラム・シュミットの正規直交化法
Gram-Schmidt orhonormalisation
main.c
-----------------------------------------------------*/
#include "main.h"
int main(void)
{
int i, j, k;
double d;
double u[DIM][DIM];
double uT[DIM][DIM]; /* uの転置行列, transposed matrix of matrix u */
double e[DIM][DIM];
/*
u[i][j] =
｜ u[0][0] u[0][1] u[0][2] ｜
｜ ｜
｜ u[1][0] u[1][1] u[1][2] ｜
｜ ｜
｜ u[2][0] u[2][1] u[2][2] ｜
*/
/*----------------------------------------------------------------------*/
/*
Matlabを用いた正規直交基底の計算
Calculation of orthonormal bases using Matlab
https://jp.mathworks.com/help/symbolic/orth.html
A = sym([2 -3 -1; 1 1 -1; 0 1 -1]);
double( A )
B = double( orth(A) )
2 -3 -1
1 1 -1
0 1 -1
*/
/*
u[0][0] = 2.0; u[0][1] = -3.0; u[0][2] = -1.0;
u[1][0] = 1.0; u[1][1] = 1.0; u[1][2] = -1.0;
u[2][0] = 0.0; u[2][1] = 1.0; u[2][2] = -1.0;
*/
/*----------------------------------------------------------------------*/
/*
other example1
https://www.momoyama-usagi.com/entry/math-linear-algebra10
1 2 1
-1 -1 -1
0 -2 -2
*/
u[0][0] = 1.0; u[0][1] = 2.0; u[0][2] = 1.0;
u[1][0] = -1.0; u[1][1] = -1.0; u[1][2] = -1.0;
u[2][0] = 0.0; u[2][1] = -2.0; u[2][2] = -2.0;
/*----------------------------------------------------------------------*/
/*
other example2
-1 < u[i][j] < 1
*/
/*
srand((unsigned int)time(NULL));
for ( i = 0; i < DIM; i++ ) {
for ( j = 0; j < DIM; j++ ) {
u[i][j] = 2.0 * rand()/(double)RAND_MAX - 1.0;
}
}
*/
/*----------------------------------------------------------------------*/
// 行列 u
// matrix u
printf("# org u\n");
for ( i = 0; i < DIM; i++ ) {
for ( j = 0; j < DIM; j++ ) {
printf("%lf ", u[i][j]);
}
printf("\n");
}
printf("\n");
// 行列uの各列
// each column of matrix u
for ( j = 0; j < DIM; j++ ) {
printf("column %d: u%d\n", j, j);
for ( i = 0; i < DIM; i++ ) {
printf(" %lf\n", u[i][j]);
}
}
printf("#-------------------------------------------\n");
/*----------------------------------------------------------------------*/
// グラム・シュミットの正規直交化法の計算
// applying Gram-Schmidt orhonormalisation
// 任意の行列uを入力する。正規直交化された行列uと、直交化だけされた行列eが得られる。
// Input an arbitrary matrix u. Then, the orthonormalised matrix u
// and the only orthonormalised matrix e are obtained.
gram_schmidt_orth( u, e );
/*----------------------------------------------------------------------*/
// 正規直交化された行列u
// orthonormalised matrix u
printf("# orthonormalised u = e/|e|\n");
for ( i = 0; i < DIM; i++ ) {
for ( j = 0; j < DIM; j++ ) {
printf(" %lf ", u[i][j]);
}
printf("\n");
}
printf("#-------------------------------------------\n");
printf("# orthogonalised matrix e (not normalised)\n");
for ( i = 0; i < DIM; i++ ) {
for ( j = 0; j < DIM; j++) {
printf(" %lf ", e[i][j]);
}
printf("\n");
}
printf("#-------------------------------------------\n");
printf("# u = e/|e|\n");
for ( j = 0; j < DIM; j++ ) {
d = 0.0;
for ( i = 0; i < DIM; i++ ) {
d += e[i][j] * e[i][j];
}
printf("column %d: e%d/|e%d|\n", j, j, j);
printf(" |e|: %lf\n", sqrt(d) );
for ( i = 0; i < DIM; i++ ) {
printf(" %lf/%lf = %lf\n", e[i][j], sqrt(d), e[i][j]/(double)sqrt(d));
}
}
printf("#---------------------------------------------------------------------\n");
// 正規直交化された行列uの転置行列uT
// Transposed matrix uT of the orthonormalised matrix u
for ( i = 0; i < DIM; i++ ) {
for ( j = 0; j < DIM; j++ ) {
uT[j][i] = u[i][j];
}
}
printf("# transposed matrix uT of u\n");
for ( i = 0; i < DIM; i++ ) {
for ( j = 0; j < DIM; j++ ) {
printf(" %lf ", uT[i][j]);
}
printf("\n");
}
printf("#-------------------------------------------\n");
// 正規直交化したuとその転置行列uTの積が単位行列になることを確認
// check that the product of the orthonormalised matrix u
// and its transposed matrix uT is the unit matrix
printf("# check uT * u = unit matrix\n");
for ( i = 0; i < DIM; i++ ) {
for ( j = 0; j < DIM; j++ ) {
d = 0.0;
for ( k = 0; k < DIM; k++ ) {
d += uT[i][k] * u[k][j];
}
printf(" %lf ", d);
}
printf("\n");
}
printf("#-------------------------------------------\n");
// 正規直交化したuの各列の2ノルムが1であることを確認
// check the L2 norm of each column of the orthonormalised matrix u is 1
printf("# L2 norm of each column of orthonormalised u\n");
for ( j = 0; j < DIM; j++ ) {
d = 0.0;
for ( i = 0; i < DIM; i++ ) {
d += u[i][j] * u[i][j];
}
printf("column %d: %lf\n", j, d);
}
printf("#-------------------------------------------\n");
// 正規直交化したuの2つの列の内積を計算し、ゼロであることを確認
// calculate the inner product of two columns of the orthonormalised u
// and check that the product is zero
printf("# product of two columns of u\n");
for ( j = 0; j < DIM-1; j++ ) {
for ( k = j; k < DIM-1; k++ ) {
d = 0.0;
for ( i = 0; i < DIM; i++ ) {
d += u[i][j] * u[i][k+1];
}
printf("column %d * column %d: %lf\n", j, k+1, d);
}
}
return 0;
}