1 #ifdef HAVE_PETSC
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <assert.h>
5 #include <math.h>
6
7 #include <mpi.h>
8 #include <petsc.h>
9
10 #include "common_c.h"
11 #include "FCMangle.h"
12
13 #define fillsparsecpetscs FortranCInterface_GLOBAL_(fillsparsecpetscs, FILLSPARSECPETSCS)
14 #define fillsparsecpetscc FortranCInterface_GLOBAL_(fillsparsecpetscc, FILLSPARSECPETSCC)
15
16 #define COLMAJ2D(row,col,numrow) (row-1)+(col-1)*numrow
17 #define COLMAJ3D(a,b,c,amax,bmax,cmax) (a-1)+amax*((b-1)+bmax*(c-1))
18 #define ROWMAJ2D_ONE(row,col,numcol) (row-1)*numcol+(col-1)
19 typedef long long int gcorp_t;
20
fillsparsecpetscs(gcorp_t * ieng,double * EGmass,Mat * lhsP)21 void fillsparsecpetscs(gcorp_t* ieng, double* EGmass, Mat* lhsP)
22 {
23 int npro = propar.npro;
24 int nshl = shpdat.nshl;
25 double* mb = (double*) malloc(sizeof(double)*nshl*nshl); //block to insert
26 int e,i,j,aa; //following along with fillsparse.f
27 PetscInt* locat = (PetscInt*) malloc(sizeof(PetscInt)*nshl);
28 for(e=0;e<npro;e++)
29 {
30 for(aa=0;aa<nshl;aa++) locat[aa]=ieng[e+npro*aa]-1;
31 // for(aa=0;aa<nshl;aa++) assert(locat[aa]>=0);
32 for (i=0; i<nshl; i++) { // fill up Ke with respective egmass
33 for (j=0; j<nshl; j++) {
34 mb[nshl*i + j] = EGmass[e + npro*(i + nshl*j)];
35 }
36 }
37 //MatSetValuesBlocked(*lhsP, nshl , locat, nshl, locat, mb, ADD_VALUES);
38 PetscInt petsc_nshl;
39 petsc_nshl = (PetscInt) nshl;
40 MatSetValues(*lhsP, petsc_nshl , locat, petsc_nshl, locat, mb, ADD_VALUES);
41 }
42 free(mb);
43 free(locat);
44 }
fillsparsecpetscc(gcorp_t * ieng,double * EGmass,Mat * lhsP)45 void fillsparsecpetscc(gcorp_t* ieng, double* EGmass, Mat* lhsP)
46 {
47 int npro = propar.npro;
48 int nshl = shpdat.nshl;
49 int nedof = conpar.nedof;
50 double* mb = (double*) malloc(sizeof(double)*nedof*nedof); //block to insert
51 int e,i,j,aa; //following along with fillsparse.f
52 //int* locat = (int*) malloc(sizeof(int)*nshl);
53 PetscInt* locat = (PetscInt*) malloc(sizeof(PetscInt)*nshl);
54 for(e=0;e<npro;e++)
55 {
56 for(aa=0;aa<nshl;aa++) locat[aa]=ieng[e+npro*aa]-1;
57 // for(aa=0;aa<nshl;aa++) assert(locat[aa]>=0);
58 for (i=0; i<nedof; i++) { /* fill up Ke with respective egmass */
59 for (j=0; j<nedof; j++) {
60 mb[nedof*i + j] = EGmass[e + npro*(i + nedof*j)];
61 }
62 }
63 //MatSetValuesBlocked(*lhsP, nshl , locat, nshl, locat, mb, ADD_VALUES);
64 PetscInt petsc_nshl;
65 petsc_nshl = (PetscInt) nshl;
66 MatSetValuesBlocked(*lhsP, petsc_nshl , locat, petsc_nshl, locat, mb, ADD_VALUES);
67 }
68 free(mb);
69 free(locat);
70 }
71 #endif
72
73