#include <stdio.h>
#include <math.h>
/* house.c : 1999 3 26 by T.Oguni */
#define N 3
main(){
  double A[N][N]={3.,2.,1.,2.,3.,2.,1.,2.,3.};
  double Q[N][N], B[N][N], C[N][N], QQ[N][N], AA[N][N];
  int i, j, n=3;
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      AA[i][j]=A[i][j];}}
  house(A,Q,n);
  printf("Matrix Q\n");
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("  %f", Q[j][i]); }
    printf("\n"); }
  printf("Matrix R\n");
  matmul(Q,AA,A,n);
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("  %f", A[i][j]); }
    printf("\n"); }
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      B[i][j]=Q[j][i]; }}
  matmul(B,Q,C,n);
  printf("check of orthogonality\n");
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("  %f",C[i][j]); }
    printf("\n"); }
}
/* function */
house(A,Q,n)
double A[N][N], Q[N][N];
int n;
{
  int i, j, k;
  double W[N][N], E[N][N], Qi[N][N], B[N][N], C[N][N];
  double ai[N], wi[N], QQ[N][N], AA[N][N];
  double tol=0.0000001, s, ss, si, alpha;
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      Q[i][j]=0.0; E[i][j]=0.0; }
    Q[i][i]=1.0; E[i][i]=1.0; }
  
  for (i=0; i<n-1; i++){
    for (k=0; k<i; k++){
      ai[k]=0.0; }
    for (k=i; k<n; k++){
      ai[k]=A[k][i]; }
    for (k=0; k<n; k++){
      wi[k]=ai[k]; }
    ss=0.0;
    for (k=i; k<n; k++){
      ss+=ai[k]*ai[k]; }
    s=sqrt(ss); si=1.0;
    if (ai[i]<0.0) si=-1.0;
    wi[i]=ai[i]+si*s;
    ss=0.0;
    for (k=0; k<n; k++){
      ss+=wi[k]*wi[k]; }
    alpha=2.0/ss;
    for (j=0; j<n; j++){
      for (k=0; k<n; k++){
        W[k][j]=wi[k]*wi[j]; }
    }
    for (k=0; k<n; k++){
      for (j=0; j<n; j++){
        Qi[k][j]=E[k][j]-alpha*W[k][j]; }
    }
    matmul(Qi,Q,QQ,n);
    matmul(Qi,A,AA,n);
    for (k=0; k<n; k++){
      for (j=0; j<n; j++){
        Q[k][j]=QQ[k][j]; }
      for (j=0; j<n; j++){
        A[k][j]=AA[k][j]; }
    }
  }
  return;
}
/* matmul */
matmul(A,B,C,n)
double A[N][N], B[N][N], C[N][N];
int n;
{
  int i, j, k;
  double s;
  for (k=0; k<n; k++){
    for (i=0; i<n; i++){
      s=0.0;
      for (j=0; j<n; j++){
        s+=A[i][j]*B[j][k];
      }
      C[i][k]=s;
    }
  }
  return;
}