=====================================================================
			Source Code
=====================================================================

#include "stdafx.h"
#include <math.h>
#include <stdio.h>

void householder(double* a, double* q, double* r, int row, int column, int maxi, int maxj);


int _tmain(int argc, _TCHAR* argv[])
{

	int i, m, n, j, k;

	double a[8][4]={7.0605, 9.5022, 4.4559, 1.19,
			0.3183, 0.3445, 6.4631, 4.9836,
			2.7692, 4.3874, 7.0936, 9.5974,
			0.4617, 3.8156, 7.5469, 3.4039,
			0.9713, 7.6552, 2.7603, 5.8527,
			8.2346, 7.952,  6.797,  2.2381,
			6.9483, 1.8687, 6.551,  7.5127,
			3.171,  4.8976, 1.6261, 2.551};

	double q[8][8];
	double r[8][4];
	double qr[8][4];
	double x[4];
	double y[8];
	double xx[5];
	double yy[10];

	double a2[10][5];
	double q2[10][10];
	double r2[10][5];
	double qr2[10][5];
	double a2x[10];
	double res[10];

	double b[8] =  {0.7585,
			0.5395,
			5.3080,
			7.7917,
			9.3401,
			1.2991,
			5.6882,
			4.6939};

	double bb[10] ={0.1125,
			0.2227,
			0.3286,
			0.4284,
			0.5205,
			0.6039,
			0.6778,
			0.7421,
			0.7969,
			0.8427};

	for(i=0; i<8; i++){
		for(j=0; j<8; j++){
			q[i][j] = 0.0;
		}
	}

	for(i=0; i<8; i++){
		for(j=0; j<4; j++){
			r[i][j] = 0.0;
		}
	}

/* 1 */
	printf("1\n");

	householder((double*)a, (double*)q, (double*)r, 8, 4, 8, 4);

	printf(" matrix Q\n");
	for(i=0; i<8; i++){
		for(j=0; j<8; j++){
			printf(" %8.4f", q[i][j]);
		}
		printf("\n");
	}
	printf(" matrix R\n");
	for(i=0; i<8; i++){
		for(j=0; j<4; j++){
			printf(" %8.4f", r[i][j]);
		}
		printf("\n");
	}
	// Q*R
	for(n=0; n<8; n++){
		for(m=0; m<4; m++){
			qr[n][m] = 0;
		}
	}
	for(i=0; i<4; i++){
		for(n=0; n<8; n++){
			for(m=0; m<8; m++){
				qr[n][i] += q[n][m] * r[m][i];
			}
		}
	}
	printf(" Q*R\n");
	for(n=0; n<8; n++){
		for(m=0; m<4; m++){
			printf(" %8.4f", qr[n][m]);
		}
		printf("\n");
	}

	// y=Qt*b
	for(i=0; i<8; i++){
		y[i] = 0;
	}
	for(i=0; i<8; i++){
		for(j=0; j<8; j++){
			y[i] += q[j][i] * b[j];
		}
	}

	// solve for x (back substitution)
	for (i = 0; i <4; i++){
		x[i]=y[i];
	}
	for (i=4-1; i>=0; i--) {
		for (j=i+1; j< 4; j++){
			x[i] -= r[i][j]*x[j];
		}
		x[i] /= r[i][i];
	}

	printf("x=\n");
	for(i=0; i<4; i++){
			printf(" %7.4f", x[i]);
		printf("\n");
	}
	printf("\n");

/* 2.1 */
	printf("2.1\n");

	for (i=0; i<10; i++){
		for (j=4; j>=0; j--){
			if(j == 0){
				if(i != 9){
					a2[i][4-j] = i+1;
				}
				else{
					// row10 is 1s
					a2[i][4-j] = 1;
				}
			}
			else{
				a2[i][4-j] = pow((i+1.0)*0.1, j);
			}
		}
	}
	printf("coefficient matrix\n");
	for(i=0; i<10; i++){
		for(j=0; j<5; j++){
			printf(" %8.4f", a2[i][j]);
		}
		printf("\n");
	}

	for(i=0; i<10; i++){
		for(j=0; j<10; j++){
			q2[i][j] = 0.0;
		}
	}
	for(i=0; i<10; i++){
		for(j=0; j<5; j++){
			r2[i][j] = 0.0;
		}
	}
	

	householder((double*)a2, (double*)q2, (double*)r2, 10, 5, 10, 5);

	printf(" matrix Q\n");
	for(i=0; i<10; i++){
		for(j=0; j<10; j++){
			printf(" %8.4f", q2[i][j]);
		}
		printf("\n");
	}
	printf(" matrix R\n");
	for(i=0; i<10; i++){
		for(j=0; j<5; j++){
			printf(" %8.4f", r2[i][j]);
		}
		printf("\n");
	}

	// Q*R
	for(n=0; n<10; n++){
		for(m=0; m<5; m++){
			qr2[n][m] = 0;
		}
	}
	for(i=0; i<5; i++){
		for(n=0; n<10; n++){
			for(m=0; m<10; m++){
				qr2[n][i] += q2[n][m] * r2[m][i];
			}
		}
	}
	printf(" Q*R\n");
	for(n=0; n<10; n++){
		for(m=0; m<5; m++){
			printf(" %8.4f", qr2[n][m]);
		}
		printf("\n");
	}

	// y
	for(i=0; i<10; i++){
		yy[i] = 0;
	}
	for(i=0; i<10; i++){
		for(j=0; j<10; j++){
			yy[i] += q2[j][i] * bb[j];
		}
	}

	// x
	for (i = 0; i <5; i++){
		xx[i]=yy[i];
	}
	for (i=5-1; i>=0; i--) {
		for (j=i+1; j< 5; j++){
			xx[i] -= r2[i][j]*xx[j];
		}
		xx[i] /= r2[i][i];
	}
	
	printf("coefficents are (descending order)\n");
	for(i=0; i<5; i++){
			printf(" %7.4f", xx[i]);
	}
	printf("\n\n");

/* 2.2 */
	printf("2.2\n");
	for(i=0; i<10; i++){
		a2x[i]=0;
		for (k = 0; k < 5; k++) {
			a2x[i]+=a2[i][k]*xx[k];
		}
	}
	for(i=0; i<10; i++){
		res[i]=bb[i]-a2x[i];
	}
	printf("residual\n");
	for(i=0; i<10; i++){
			printf(" %7.4f", res[i]);
			printf("\n");
	}
	return 0;
}

void householder(double* a, double* q, double* r, int row, int column, int maxi, int maxj)
{
	int i, j, k;
	double aa1;	// |A1|
	double uu;	// |u1|^2

	double* u1 = new double[row];
	double* e1 = new double[row];
	double* uut = new double[row*row];
	double* h1 = new double[row*row];
	double* I = new double[row*row];
	double* h1a = new double[row*column];
	double* a2 = new double[(row-1)*(column-1)];
	double* h = new double[maxi*maxi];
	double* q0 = new double[maxi*maxi];

	for(i=0; i<row; i++){
		e1[i] = 0.0;
	}
	e1[0] = 1.0;

	for(i=0; i<row; i++){
		for(j=0; j<row; j++){
			if(i == j){
				I[i*row+j] = 1.0;
			}
			else{
				I[i*row+j] = 0.0;
			}
		}
	}
	for(i=0; i<maxi; i++){
		for(j=0; j<maxi; j++){
			if(i == j){
				h[i*maxi+j] = 1.0;
			}
			else{
				h[i*maxi+j] = 0.0;
			}
		}
	}
	for(i=0; i<maxi; i++){
		for(j=0; j<maxi; j++){
			q0[i*maxi+j] = q[i*maxi+j];
		}
	}
	for(i=0; i<maxi; i++){
		for(j=0; j<maxi; j++){
			q[i*maxi+j] = 0;
		}
	}

	// norm
	aa1 = 0;
	for(i=0; i<row; i++){
		aa1 += pow(a[i*column], 2.0);
	}
	aa1 = sqrt(aa1);

	// u
	for(i=0; i<row; i++){
		u1[i] = a[i*column] + aa1 * e1[i];
	}

	// H1
	// |u1|^2
	uu = 0;
	for(i=0; i<row; i++){
		uu += pow(u1[i], 2.0);
	}

	// u1*u1t
	for(i=0; i<row; i++){
		for(j=0; j<row; j++){
			uut[i*row+j] = u1[i] * u1[j];
		}
	}

	// H1
	for(i=0; i<row; i++){
		for(j=0; j<row; j++){
			h1[i*row+j] = I[i*row+j] - (2.0/uu)*uut[i*row+j];
		}
	}

	// h1->h
	k = 0;
	for(i=(maxi-row); i<maxi; i++){
		for(j=(maxi-row); j<maxi; j++){
			h[i*maxi+j] = h1[k];
			k++;
		}
	}

	// q
	if(row != maxi){
		for(i=0; i<maxi; i++){
			for(j=0; j<maxi; j++){
				for(k=0; k<maxi; k++){
					q[j*maxi+i] += q0[j*maxi+k] * h[k*maxi+i];
				}
			}
		}
	}
	else{
		for(i=0; i<maxi; i++){
			for(j=0; j<maxi; j++){
				q[i*maxi+j] = h[i*maxi+j];
			}
		}
	}

	// H1*A
	for(i=0; i<row; i++){
		for(j=0; j<column; j++){
			h1a[i*column+j] = 0;
		}
	}
	for(i=0; i<column; i++){
		for(j=0; j<row; j++){
			for(k=0; k<row; k++){
				h1a[j*column+i] += h1[j*row+k] * a[k*column+i];
			}
		}
	}

	// r
	for(i=0; i<row; i++){
		for(j=0; j<column; j++){
			r[((maxi-row)+i)*maxj+(maxj-column)+j] = h1a[i*column+j];
		}
	}

	if(column > 1){
		// h -> a2
		for(i=0; i<(row-1); i++){
			for(j=0; j<(column-1); j++){
				a2[i*(column-1)+j] = h1a[(i+1)*column+j+1];
			}
		}
		householder(a2, q, r, row-1, column-1, maxi, maxj);
	}

	delete [] u1;
	delete [] e1;
	delete [] uut;
	delete [] h1;
	delete [] I;
	delete [] h1a;
	delete [] a2;
	delete [] h;
	delete [] q0;

	return ;
}


=====================================================================
			Results
=====================================================================

1
 matrix Q
  -0.5192  -0.2621   0.2785   0.2631  -0.0529  -0.6925  -0.0232  -0.1790
  -0.0234  -0.0031  -0.6006   0.0827   0.3138  -0.3223   0.5595   0.3414
  -0.2036  -0.1712  -0.3867  -0.4860   0.2535  -0.1534  -0.6624   0.1300
  -0.0340  -0.3478  -0.5684   0.4393  -0.5125   0.1889  -0.1671  -0.1886
  -0.0714  -0.6931   0.0679  -0.4275   0.1074   0.2404   0.3989  -0.3140
  -0.6056   0.0183   0.0512   0.3580   0.4410   0.5374  -0.0410   0.1312
  -0.5110   0.5165  -0.2183  -0.3847  -0.3356   0.0571   0.2267  -0.3304
  -0.2332  -0.1829   0.1853  -0.1893  -0.5032   0.0864   0.0858   0.7576
 matrix R
 -13.5985 -13.4237 -12.2052  -9.0114
  -0.0000  -9.6611  -3.7297  -3.7561
   0.0000   0.0000 -10.2666  -8.9629
  -0.0000  -0.0000   0.0000  -7.5176
  -0.0000  -0.0000   0.0000   0.0000
   0.0000  -0.0000   0.0000  -0.0000
   0.0000   0.0000   0.0000   0.0000
   0.0000  -0.0000   0.0000   0.0000
 Q*R
   7.0605   9.5022   4.4559   1.1900
   0.3183   0.3445   6.4631   4.9836
   2.7692   4.3874   7.0936   9.5974
   0.4617   3.8156   7.5469   3.4039
   0.9713   7.6552   2.7603   5.8527
   8.2346   7.9520   6.7970   2.2381
   6.9483   1.8687   6.5510   7.5127
   3.1710   4.8976   1.6261   2.5510
x=
 -0.4958
  0.5741
 -0.0304
  0.7339

2.1
coefficient matrix
   0.0001   0.0010   0.0100   0.1000   1.0000
   0.0016   0.0080   0.0400   0.2000   2.0000
   0.0081   0.0270   0.0900   0.3000   3.0000
   0.0256   0.0640   0.1600   0.4000   4.0000
   0.0625   0.1250   0.2500   0.5000   5.0000
   0.1296   0.2160   0.3600   0.6000   6.0000
   0.2401   0.3430   0.4900   0.7000   7.0000
   0.4096   0.5120   0.6400   0.8000   8.0000
   0.6561   0.7290   0.8100   0.9000   9.0000
   1.0000   1.0000   1.0000   1.0000   1.0000
 matrix Q
  -0.0001  -0.0052  -0.1033  -0.5948   0.2715   0.2968   0.1227  -0.1672  -0.4558  -0.4721
  -0.0012  -0.0366  -0.2993  -0.5624   0.1084  -0.2069  -0.3208  -0.2825   0.1523   0.5764
  -0.0063  -0.1065  -0.4562  -0.2541  -0.1357  -0.1828   0.1067   0.3974   0.5601  -0.4256
  -0.0198  -0.2122  -0.4938   0.0742  -0.2442   0.3997   0.4325   0.2315  -0.2677   0.4167
  -0.0483  -0.3358  -0.3841   0.2622  -0.1367  -0.5293  -0.2709  -0.1218  -0.4867  -0.2113
  -0.1001  -0.4446  -0.1512   0.2452   0.1303   0.5718  -0.4496  -0.2541   0.2829  -0.1365
  -0.1854  -0.4906   0.1292   0.0541   0.3635  -0.2562   0.5941  -0.3379   0.1967   0.0411
  -0.3163  -0.4107   0.3293  -0.1846   0.2332  -0.0374  -0.2333   0.6563  -0.1657   0.1425
  -0.5066  -0.1268   0.2694  -0.2489  -0.7270   0.0417   0.0254  -0.2375   0.0252  -0.0693
  -0.7721   0.4542  -0.2823   0.1789   0.2930   0.0000   0.0000   0.0000   0.0000  -0.0000
 matrix R
  -1.2951  -1.3961  -1.5276  -1.7051 -10.1014
   0.0000  -0.1716  -0.4409  -0.8913 -13.0008
  -0.0000  -0.0000  -0.0735  -0.3732  -1.1916
   0.0000  -0.0000   0.0000  -0.0953  -2.5626
   0.0000  -0.0000   0.0000  -0.0000  -2.6368
  -0.0000  -0.0000   0.0000  -0.0000   0.0000
  -0.0000   0.0000  -0.0000   0.0000   0.0000
  -0.0000   0.0000  -0.0000   0.0000   0.0000
  -0.0000  -0.0000  -0.0000  -0.0000  -0.0000
  -0.0000   0.0000   0.0000  -0.0000   0.0000
 Q*R
   0.0001   0.0010   0.0100   0.1000   1.0000
   0.0016   0.0080   0.0400   0.2000   2.0000
   0.0081   0.0270   0.0900   0.3000   3.0000
   0.0256   0.0640   0.1600   0.4000   4.0000
   0.0625   0.1250   0.2500   0.5000   5.0000
   0.1296   0.2160   0.3600   0.6000   6.0000
   0.2401   0.3430   0.4900   0.7000   7.0000
   0.4096   0.5120   0.6400   0.8000   8.0000
   0.6561   0.7290   0.8100   0.9000   9.0000
   1.0000   1.0000   1.0000   1.0000   1.0000
coefficents are (descending order)
  0.1627 -0.4699  0.0233  1.1266 -0.0000

2.2
residual
  0.0001
 -0.0000
 -0.0001
  0.0000
  0.0000
  0.0001
 -0.0000
 -0.0000
  0.0000
  0.0000