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

int main(void)
{

  int i,j,k,n=4,m=3;
  double x[4][3], q[4][3], r[3][3];

  FILE *fi1, *fi2;

  void mgs(double x[][m],int n,int m,double q[][m],double r[][m]); 
  double dot(double x[],double y[],int length);


  fi1=fopen("mgsdata.txt","r");
  fi2=fopen("mgsdata.txt","r");
  i=0;
  while (fscanf(fi2, "%lf %lf %lf", &x[i][0],&x[i][1],&x[i][2]) != EOF) {
 
    fscanf(fi1, "%lf %lf %lf", &x[i][0],&x[i][1],&x[i][2]);
    i=i+1;
  } 
  n=i;

  printf("X:\n");
  for (i=0; i<n; i++){
   printf("%f %f %f \n ", x[i][0],x[i][1],x[i][2]);
  }

  fclose(fi1);
  fclose(fi2);

  mgs(x,4,3,q,r);

  printf("\nQ:\n");
  for (i=0; i<4; i++){
    printf("%f %f %f \n ", q[i][0],q[i][1],q[i][2]);
  }
  printf("\nR:\n");
  for (i=0; i<3; i++){
    printf("%f %f %f \n ", r[i][0],r[i][1],r[i][2]);
  }

  return 0;

}

double dot(double x[],double y[],int length){

  int i;
  double z;

  z=0.0;
  for(i=0;i<length;i++){
    z+=x[i]*y[i];
  }
  /*  printf("z: %f \n",z);*/
  return z;

}

void mgs(double x[][3],int n,int m,double q[][3],double r[][3]){
int i,j,k,l;
double qk[4],qj[4];

/* initialize r and q */
for(i=0;i<m;i++){
  for(j=0;j<m;j++){
    r[i][j]=0.0;
  }
  r[i][i]=1.0;
}
for(i=0;i<n;i++){
  for(j=0;j<m;j++){
    q[i][j]=x[i][j];
  }
}

for(k=0;k<m;k++){

  for(i=0;i<n;i++){
    qk[i]=q[i][k];
  }
  r[k][k]=sqrt(dot(qk,qk,n));

  for(i=0;i<n;i++){
    q[i][k]=q[i][k]/r[k][k];
    qk[i]=q[i][k];
  }

  for(j=k+1;j<m;j++){
    for(i=0;i<n;i++){
      qj[i]=q[i][j];
    }
    r[k][j]=dot(qk,qj,n);

    for(l=0;l<n;l++){
      q[l][j]=q[l][j]-q[l][k]*r[k][j];
    }
  }
}
}