#include <math.h>
#include <stdio.h>
#define real double

real dot(real *a,real *b,int n){
int i; real s=0; for(i=0;i<n;i++) s+=a[i]*b[i]; return s; }

void QR_dec(real **a,real **r,int n,int m){
int i,j,k;
for(i=0;i<m;i++){
	r[i][i]=sqrt(dot(a[i],a[i],n));
	if(r[i][i]==0) printf("QR_dec: singular matrix\n");
	for(k=0;k<n;k++) a[i][k]/=r[i][i];
	for(j=i+1;j<m;j++) {
		r[j][i]=dot(a[i],a[j],n);
		for(k=0;k<n;k++){
				a[j][k]-=r[j][i]*a[i][k];
				}
		}
	}
}

void QR_back(real **q,real **r,int n,int m,real *bb){
int i,j; real *b,s;

b = (real *)malloc(n*sizeof(real));
for(i=0;i<m;i++){b[i]=0; for(j=0;j<n;j++) b[i]+=q[i][j]*bb[j];}
for(i=m-1;i>=0;i--){s=0; for(j=i+1;j<n;j++) s+=r[j][i]*bb[j];
				bb[i]=(b[i]-s)/r[i][i];}
free(b);
}