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

void getmat(int n, double l[n][n], double u[n][n]) {
  for( int i=0;i<n; i++) {
    for ( int j=0; j<=i; j++ ) {
      l[i][j] = 1.0/(1+fabs(i-j));
    }
    for ( int j=i; j<n; j++ ) {
      u[i][j] = 1.0 + (i+j)%5;
    }
  }
}

// Matrix multiplication LU = L*U, where L is lower triangular and
// U is upper triangular
//
void trimul(int n, double l[n][n], double u[n][n], double lu[n][n]) {
#pragma omp parallel for schedule (dynamic,10) 
  for ( int i=0; i<n; i++ ) {
    for ( int j=0; j<n; j++ )
      lu[i][j] = 0.0;
    for ( int k=0; k<=i; k++ ) {
      for ( int j=0; j<n; j++ )
	lu[i][j] += l[i][k]*u[k][j];
    }
  }
}

double residual( int n, double a[n][n], double b[n][n], double c[n][n] ) {
  double r, cij;
  r=0.0;
  for ( int i=0; i<n; i++ ) {
    for ( int j=0; j<n; j++ ) {
      cij = 0.0;
      for ( int k=0; k<n; k++ ) 
	cij += a[i][k]*b[k][j];
      r = fmax(r,fabs(c[i][j]-cij));
    }
  }
  return r;
}

void printmat(int n, double a[n][n]) {
  for (int i=0; i<n; i++ ) {
    for (int j=0;j<n; j++ ) 
      printf("%.5lf ",a[i][j]);
    printf("\n");
  }
}

#define N 1000

double l[N][N], u[N][N], lu[N][N];

int main() {
  double t;
  getmat(N,l,u);
  t = omp_get_wtime();
  trimul(N,l,u,lu);  
  printf("elapsed time = %lf ",omp_get_wtime()-t);
  printf("residual = %le\n",residual(N,l,u,lu));
}