#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

using namespace std;
#ifndef ROWS
#define ROWS 1024
#endif
#ifndef COLS
#define COLS 1024
#endif

void gramSchmidt(float Q[][COLS], const int rows, const int cols)
{
    for(int k=0; k < cols; k++)
    {
        double tmp = 0.;
        #pragma omp parallel for reduction(+:tmp)
        for(int i=0; i < rows; i++) tmp +=  (Q[i][k] * Q[i][k]);
        tmp = sqrt(tmp);

        #pragma omp parallel for
        for(int i=0; i < rows; i++) Q[i][k] /= tmp;

        #pragma omp parallel for reduction(+:tmp)
        for(int j=k+1; j < cols; j++)
        {
            tmp=0.;
            for(int i=0; i < rows; i++) tmp += Q[i][k] * Q[i][j];
            for(int i=0; i < rows; i++) Q[i][j] -= tmp * Q[i][k];
        }
  }
}

void gramSchmidtSeq(float Q[][COLS], const int rows, const int cols)
{
    for(int k=0; k < cols; k++)
    {
        double tmp = 0.;

        for(int i=0; i < rows; i++) tmp +=  (Q[i][k] * Q[i][k]);
        tmp = sqrt(tmp);


        for(int i=0; i < rows; i++) Q[i][k] /= tmp;


        for(int j=k+1; j < cols; j++)
        {
            tmp=0.;
            for(int i=0; i < rows; i++) tmp += Q[i][k] * Q[i][j];
            for(int i=0; i < rows; i++) Q[i][j] -= tmp * Q[i][k];
        }
  }
}
void checkOctave(float A[][COLS], int rows, int cols)
{
  int found_error=0;
  for(int c1=0; c1 < cols; c1++)
    for(int c2=c1; c2 < cols; c2++) {
      double sum=0.;
      for(int i=0; i < rows; i++) sum += A[i][c1] * A[i][c2];
      if(c1 == c2) { // should be near 1 (unit length)
    if(sum < 0.9) {
      printf("Failed unit length: %d %d %g\n", c1,c2,sum);
      found_error = 1;
      exit(1);
    }
      } else { // should be very small (orthogonal)
    if(sum > 0.1) {
      printf("Failed orthogonal  %d %d %g\n", c1,c2,sum);
      found_error = 1;
      exit(1);
    }
      }
  }
  if(!found_error) printf("Check OK!\n");
}

int main(int argc, char *argv[])
{
  int rows=ROWS;
  int cols=COLS;
  float (*A)[COLS] = (float(*)[COLS])malloc(sizeof(float)*rows*cols);
  float (*B)[COLS] = (float(*)[COLS])malloc(sizeof(float)*rows*cols);
  // fill matrix A with random numbers
  for(int i=0; i < rows; i++)
    for(int j=0; j < cols; j++){
      A[i][j] = (float)rand()/ (float)RAND_MAX;
      B[i][j] = A[i][j];
    }
  printf("Done with init!\n");
  double startTime;
  double endTime;

 /* startTime = omp_get_wtime();
  gramSchmidtSeq(B, rows, cols);
  endTime = omp_get_wtime();

  printf("runtime %d %d matrix %g\n",rows, cols, (endTime-startTime));*/
  startTime = omp_get_wtime();
  gramSchmidt(A, rows, cols);
  endTime = omp_get_wtime();

  printf("runtime %d %d matrix %g\n",rows, cols, (endTime-startTime));

  startTime = omp_get_wtime();
  gramSchmidtSeq(B, rows, cols);
  endTime = omp_get_wtime();

  printf("runtime %d %d matrix %g\n",rows, cols, (endTime-startTime));
  checkOctave(A, rows, cols);
  free(A);
}