#include <stdio.h> 
#include <math.h>
/* hess.c : 1999 3 27 by T.Oguni */
#define N 3
main(){
  int n=3, i, j;
  double A[N][N]={1.,4.,3.,2.,5.,4.,1.,-3.,-2.}, Q[N][N];
  hess(A,Q,n);
  printf("Matrix Q\n");
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("  %f", Q[i][j]); }
    printf("\n"); }
  printf("Matrix R\n");
  for (i=0; i<n; i++){
    for (j=0; j<n; j++){
      printf("  %f", A[i][j]); }
    printf("\n"); }
}
/* function */
hess(A,Q,n)
double A[N][N], Q[N][N];
int n;
{
  double W[N][N],E[N][N],Qi[N][N],Qj[N][N],AA[N][N],QQ[N][N];
  double wi[N], ai[N];
  double ss, s, alpha, t, si;
  int i, j, k;
  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-2; i++){
    for (k=0; k<n; k++){ ai[k]=0.0; wi[k]=0.0; }
    for (k=i+1; k<n; k++){ ai[k]=A[k][i]; wi[k]=ai[k]; }
    ss=0.0;
    for (k=i+1; k<n; k++){ ss+=ai[k]*ai[k]; }
    s=sqrt(ss);
    si=1.0;
    if (ai[i+1]<0) si=-1.0;
    wi[i+1]=ai[i+1]+si*s;
    t=0.0;
    for (k=0; k<n; k++){ t+=wi[k]*wi[k]; }
    alpha=2.0/t;
    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++){
        Qj[k][j]=Qi[j][k]; Q[k][j]=QQ[k][j]; }}
    matmul(AA,Qj,A,n);
  }
  return;
}
/* function */
matmul(A,B,C,n)
double A[N][N], B[N][N], C[N][N];
int n;
{
/* ŒvŽZ•”•ª */
  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;
    }
  }
}