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

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

int qrdec(double** a,double** r,int n,int m){
// QR-decomposition routine.
// a : input : n times m matrix to be decomposed
// a : output : the n times m orthogonal Q matrix
// r : output : the m times m right triangular R matrix
   for(int i=0;i<m;i++){
      r[i][i]=sqrt(dot(a[i],a[i],n));
      if(r[i][i]==0) {fprintf(stderr,"qrdec: singular matrix\n"); return 1;}
      for(int k=0;k<n;k++) a[i][k]/=r[i][i];
      for(int j=i+1;j<m;j++) {
         r[j][i]=dot(a[i],a[j],n);
         for(int k=0;k<n;k++){
            a[j][k]-=r[j][i]*a[i][k]; } } }
   return 0; }