mumps_solver.c 2.04 KB
Newer Older
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
1
2
#include <stdio.h>
#include <stdlib.h>
3
#include <string.h>
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#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 */
Jonathan Lambrechts's avatar
clean    
Jonathan Lambrechts committed
59
  id->ICNTL(16)=0;
Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
60
61
  return id;
}
62

Jonathan Lambrechts's avatar
Jonathan Lambrechts committed
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
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);
}