#include "libscl.h"
#include <omp.h>

using namespace std;
using namespace scl;

int main(int argc, char** argp, char** envp)
{
  realmat a, b;
  if(!vecread("a.dat",a) || !vecread("b.dat",b)) error("Read failed");
  if (a.ncol() != b.nrow()) error("Not conformable");

  const INTEGER arows = a.nrow();
  const INTEGER acols = b.nrow();
  const INTEGER bcols = b.ncol();

  realmat r(arows,bcols,0.0);

  INTEGER chunk = 2;

  #pragma omp parallel shared(chunk) 
  {
    #pragma omp for schedule(static, chunk)

    for (INTEGER j=1; j<=bcols; ++j) {
      for (INTEGER k=1; k<=acols; ++k) {
        for (INTEGER i=1; i<=arows; ++i) {
          r(i,j) += a(i,k)*b(k,j);
        }
      }
    }

  }

  std::cout << a << b << r << '\n';

  return 0;
}