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 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 } 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