#include "libscl.h"
#include <pthread.h>   // Header for pthread
#include <unistd.h>    // Header for sysconf
using namespace scl;

namespace {

  realmat a;
  REAL fnorm;
  pthread_mutex_t mutexsum;

  void* mult(void* arg)
  {
     INTEGER* jptr = (INTEGER*)(arg);
     INTEGER n = a.nrow();
     REAL* t = &a(1,*jptr);
     REAL* top = t + n;
     REAL sum = 0.0;
     while (t<top) sum += pow(*t++,2);
     pthread_mutex_lock(&mutexsum);
     fnorm += sum;
     pthread_mutex_unlock(&mutexsum);
     pthread_exit(NULL);
  }
}

int main (int argc, char** argp, char** envp)
{
  if(!vecread("a.dat",a)) error("Read failed");
  
  INTEGER num_threads = a.ncol();

  pthread_t threads[num_threads];

  pthread_attr_t attr;
  pthread_attr_init(&attr);
  pthread_attr_setdetachstate(&attr,PTHREAD_CREATE_JOINABLE);

  pthread_mutex_init(&mutexsum, NULL);

  INTEGER arg[num_threads];
  
  fnorm = 0.0;

  for (INTEGER t=0; t<num_threads; ++t) {
    arg[t] = t+1;
    int rc = pthread_create(&threads[t], &attr, &mult, (void*)(&arg[t]));
    if (rc) error("Cannot create threads");
  }
  
  pthread_attr_destroy(&attr);

  void* status;
  for (INTEGER t=0; t<num_threads; ++t) {
    int rc = pthread_join(threads[t], &status);
    if (rc) error("Cannot join threads");
  }

  fnorm = sqrt(fnorm);

  std::cout<<"The Frobenius norm of"<< a <<"\nis "<< fnorm <<'\n';

  pthread_mutex_destroy(&mutexsum);
  pthread_exit(NULL);
}