#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "mpi.h"
#include "dmumps_c.h"
#include "string.h"
#define JOB_INIT -1
#define JOB_END -2
#define USE_COMM_WORLD -987654
#define ICNTL(I) icntl[(I)-1] /*macro s.t. indices match documentation*/

void mumps_set_options(DMUMPS_STRUC_C *id, int *options, int n){
  for(int i = 0; i < n; ++i){
    id->ICNTL(options[2*i+0]) = options[2*i+1]; 
  }
}

DMUMPS_STRUC_C *mumps_element_new(int n, int nelt, int *eltptr, int *eltvar, double *a_elt, double *rhs) {
  DMUMPS_STRUC_C *id = malloc(sizeof(DMUMPS_STRUC_C));
  int myid, ierr;
  ierr = MPI_Init(NULL, NULL);
  //ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  /*Initialize a MUMPS instance. Use MPI_COMM_WORLD.*/
  id->job=JOB_INIT;
  id->par=1;
  id->sym=0;
  id->comm_fortran=USE_COMM_WORLD;
  dmumps_c(id);
  /*Define the problem on the host*/
  id->n = n;
  id->nelt = nelt;
  id->eltptr = malloc(sizeof(int)*(nelt+1));
  memcpy(id->eltptr, eltptr, sizeof(int)*(nelt+1));
  id->eltvar = eltvar;
  id->eltvar = malloc(sizeof(int)*(eltptr[nelt]-1));
  memcpy(id->eltvar, eltvar, sizeof(int)*(eltptr[nelt]-1));
  /*No outputs*/
  /*id->ICNTL(1)=-1;
  id->ICNTL(2)=-1;
  id->ICNTL(3)=-1;*/
  /* output level */
  id->ICNTL(4)=-1;
  // id->ICNTL(4)=2;

  /* element format */
  id->ICNTL(5)=1;
  /* permutes the matrix to a zero-free diagonal */
  id->ICNTL(6)=2;
  /* reordering metis */
  id->ICNTL(7)=5;
  /* transpose matrix (local matrices c->f) */
  id->ICNTL(9)=0;
  /* no scaling */
  id->ICNTL(8)=0;
  /* When significant extra fill-in is caused by numerical pivoting, increasing ICNTL(14) may help */
  /* 50 makes all the validation functional */
  id->ICNTL(14)=50;
  /* num of openmp threads */
  id->ICNTL(16)=0;
  return id;
}

void mumps_element_solve(DMUMPS_STRUC_C *id, double *a_elt, double *rhs) {
  int init = id->a_elt == NULL;
  id->a_elt = a_elt;
  id->rhs = rhs;
  id->job = (init ? 6 : 5);
  dmumps_c(id);
}

void mumps_element_delete(DMUMPS_STRUC_C *id) {
  free(id->eltptr);
  free(id->eltvar);
  id->job=JOB_END;
  dmumps_c(id);
  free(id);
}