#include <stdio.h>
#define real double

void QR_dec(real**,real**,int,int);
void QR_back(real**,real**,int,int,real*);

void print_matrix(real** a,int n,int m){
  int i,j;
  for(i=0;i<n;i++){ for(j=0;j<m;j++) printf("%g\t",a[j][i]); printf("\n");}
}
void print_vector(real* a,int n){
  int i; for(i=0;i<n;i++) printf("%g\n",a[i]);
}
void matrix_mult(real** a,int n,int m,real** b,int l,real** c){
  int i,j,k;
  for(i=0;i<n;i++) { for(j=0;j<l;j++) { c[j][i]=0;
					for(k=0;k<m;k++) c[j][i]+=a[k][i]*b[j][k]; }}
}
void mult_mat_vec(real ** a,int n,int m,real * b,real * c){
  int i,j;
  for(i=0;i<n;i++) { c[i]=0; for(j=0;j<m;j++) c[i]+=a[j][i]*b[j]; }
}

int main(void){
int i,n=3,j,k;
real **r,**a,**q,*b,*x,d;

a = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) a[i] = (real *)malloc(n*sizeof(real)); 

r = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) r[i] = (real *)malloc(n*sizeof(real)); 

q = (real **)malloc(n*sizeof(real *)); 
for (i=0; i<n ; i++) q[i] = (real *)malloc(n*sizeof(real)); 

q[0][0]=-1; q[1][0]=1; q[2][0]=1;
q[0][1]=2; q[1][1]=-2; q[2][1]=2;
q[0][2]=3; q[1][2]=3; q[2][2]=-3;

printf("matrix A:\n"); print_matrix(q,n,n);
QR_dec(q,r,n,n);
printf("triangular matrix R:\n"); print_matrix(r,n,n);
d=1; for(i=0; i<n ; i++) d*=r[i][i];
printf("determinant:\t%g\n",d);
matrix_mult(q,n,n,r,n,a);
printf("matrix product QR:\n"); print_matrix(a,n,n);

b = (real *)malloc(n*sizeof(real)); 
x = (real *)malloc(n*sizeof(real)); 

b[0]=1 ; b[1]=2 ; b[2]=3;
printf("vector b:\n"); print_vector(b,n);
QR_back(q,r,n,n,b);
printf("solution x to Ax=b:\n"); print_vector(b,n);
mult_mat_vec(a,n,n,b,x);
printf("check: vector Ax:\n"); print_vector(x,n);

}