#include "header.h"

int main(void){
  int i;
  int const Dim=4;
  gsl_matrix *A, *Q, *R, *Out;
  gsl_vector *tau;
  FILE *fp, *fpOut;
  fp =fopen("A.mat","r");
  fpOut = fopen("test","w");
  
  //Initialization 
  A = gsl_matrix_alloc(Dim,Dim);
  Q = gsl_matrix_alloc(Dim,Dim);
  R = gsl_matrix_alloc(Dim,Dim);
  Out = gsl_matrix_alloc(Dim,Dim);
  tau = gsl_vector_alloc(Dim);
  gsl_matrix_fscanf(fp, A);
  
  //QR Decomp
  for(i=0; i< 120; i++){
    gsl_linalg_QR_decomp (A, tau);
    gsl_linalg_QR_unpack (A, tau, Q, R);
    gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, R, Q, 0.0, A);
    }
  
  // Checking Output
  gsl_matrix_fprintf (fpOut, A, "%12g");
  fprintf(fpOut, "\n\n");
  
  return 1;
}