Ice Sheet System Model  4.18
Code documentation
Functions
petscpatches.h File Reference
#include <petscmat.h>
#include <petscvec.h>
#include <petscksp.h>
#include <petscsys.h>
#include "./SolverEnum.h"
#include "../../toolkitsenums.h"
#include "../../../shared/io/Comm/IssmComm.h"

Go to the source code of this file.

Functions

Vec NewVec (int size, ISSM_MPI_Comm comm, bool fromlocalsize=false)
 
Mat NewMat (int M, int N, ISSM_MPI_Comm comm)
 
Mat NewMat (int M, int N, double sparsity, ISSM_MPI_Comm comm)
 
Mat NewMat (int M, int N, int connectivity, int numberofdofspernode, ISSM_MPI_Comm comm)
 
int VecToMPISerial (double **pgathered_vector, Vec vector, ISSM_MPI_Comm comm, bool broadcast=true)
 
void MatFree (Mat *pmat)
 
void ISFree (IS *pis)
 
void VecFree (Vec *pvec)
 
void KSPFree (KSP *pksp)
 
int MatPartition (Mat *poutmatrix, Mat matrixA, double *row_partition_vector, int row_partition_vector_size, double *col_partition_vector, int col_partition_vector_size)
 
void PetscOptionsDetermineSolverType (int *psolver_type)
 
void MatMultPatch (Mat A, Vec X, Vec AX, ISSM_MPI_Comm comm)
 
void MatToSerial (double **poutmatrix, Mat matrix, ISSM_MPI_Comm comm)
 
Vec SerialToVec (double *vector, int vector_size)
 
InsertMode ISSMToPetscInsertMode (InsMode mode)
 
NormType ISSMToPetscNormMode (NormMode mode)
 
MatType ISSMToPetscMatrixType (MatrixType type)
 

Function Documentation

◆ NewVec()

Vec NewVec ( int  size,
ISSM_MPI_Comm  comm,
bool  fromlocalsize = false 
)

Definition at line 19 of file NewVec.cpp.

19  {
20 
21  int local_size;
22 
23  /*output: */
24  Vec vector=NULL;
25 
26  /*determine local size of vector: */
27  if(fromlocalsize){
28  local_size=size;
29  }
30  else{
31  local_size=DetermineLocalSize(size,comm);
32  }
33 
34  VecCreate(comm,&vector);
35  VecSetSizes(vector,local_size,PETSC_DECIDE);
36  VecSetFromOptions(vector);
37 
38  return vector;
39 }

◆ NewMat() [1/3]

Mat NewMat ( int  M,
int  N,
ISSM_MPI_Comm  comm 
)

Definition at line 21 of file NewMat.cpp.

21  {
22 
23  /*output:*/
24  Mat outmatrix=NULL;
25 
26  /*parameters: */
27  double sparsity=0.001; //default
28  int m,n;
29  int d_nz,o_nz,nnz;
30 
31  /*Determine local sizes: */
32  m=DetermineLocalSize(M,comm);
33  n=DetermineLocalSize(N,comm);
34 
35  nnz=(int)((double)M*(double)N*sparsity); //number of non zeros.
36  d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
37  o_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
38 
39  #if _PETSC_MAJOR_ == 3 && _PETSC_MINOR_ > 2
40  MatCreateAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
41  #else
42  MatCreateMPIAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
43  #endif
44 
45  return outmatrix;
46 }

◆ NewMat() [2/3]

Mat NewMat ( int  M,
int  N,
double  sparsity,
ISSM_MPI_Comm  comm 
)

Definition at line 49 of file NewMat.cpp.

49  {
50 
51  /*output:*/
52  Mat outmatrix=NULL;
53 
54  /*parameters: */
55  int m,n;
56  int d_nz,o_nz;
57  int nnz;
58 
59  /*Determine local sizes: */
60  m=DetermineLocalSize(M,comm);
61  n=DetermineLocalSize(N,comm);
62 
63  nnz=(int)((double)M*(double)N*sparsity); //number of non zeros.
64  d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
65  o_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
66 
67  #if _PETSC_MAJOR_ == 3 && _PETSC_MINOR_ > 2
68  if(sparsity==1){
69  MatCreateDense(comm,m,n,M,N,NULL,&outmatrix);
70  }
71  else{
72  MatCreateAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
73  }
74  #else
75  MatCreateMPIAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
76  #endif
77 
78  return outmatrix;
79 }

◆ NewMat() [3/3]

Mat NewMat ( int  M,
int  N,
int  connectivity,
int  numberofdofspernode,
ISSM_MPI_Comm  comm 
)

Definition at line 82 of file NewMat.cpp.

82  {
83 
84  /*output:*/
85  Mat outmatrix=NULL;
86 
87  /*parameters: */
88  int m,n;
89  int d_nz,o_nz;
90 
91  #if _PETSC_MAJOR_ >= 3
92  #if defined(_HAVE_PETSCDEV_) || _PETSC_MINOR_ >=4
93  MatType type;
94  #else
95  const MatType type;
96  #endif
97  #else
98  MatType type;
99  #endif
100 
101  /*Determine local sizes: */
102  m=DetermineLocalSize(M,comm);
103  n=DetermineLocalSize(N,comm);
104 
105  /*Figure out number of non zeros per row: */
106  d_nz=(int)connectivity*numberofdofspernode/2;
107  o_nz=(int)connectivity*numberofdofspernode/2;
108 
109  MatCreate(comm,&outmatrix);
110  MatSetSizes(outmatrix,m,n,M,N);
111  MatSetFromOptions(outmatrix);
112 
113  /*preallocation according to type: */
114  MatGetType(outmatrix,&type);
115 
116  if((strcmp(type,"mpiaij")==0) || (strcmp(type,"mpidense")==0)){
117  MatMPIAIJSetPreallocation(outmatrix,d_nz,NULL,o_nz,NULL);
118  }
119 
120  return outmatrix;
121 }

◆ VecToMPISerial()

int VecToMPISerial ( double **  pgathered_vector,
Vec  vector,
ISSM_MPI_Comm  comm,
bool  broadcast = true 
)

Definition at line 14 of file VecToMPISerial.cpp.

14  {
15 
16  int i;
17  int num_procs;
18  int my_rank;
19 
20  /*Petsc*/
21  ISSM_MPI_Status status;
22  PetscInt lower_row,upper_row;
23  int range;
24  int * idxn=NULL;
25  int buffer[3];
26 
27  /*intermediary results*/
28  double* local_vector=NULL;
29 
30  /*input*/
31  int vector_size;
32 
33  /*Output*/
34  double* gathered_vector=NULL; //Global vector holding the final assembled vector on all nodes.
35 
36  /*recover my_rank and num_procs*/
37  ISSM_MPI_Comm_size(comm,&num_procs);
38  ISSM_MPI_Comm_rank(comm,&my_rank);
39 
40  VecGetSize(vector,&vector_size);
41  if(vector_size==0){
42  *pgathered_vector=NULL;
43  return 1;
44  }
45 
46  /*Allocate gathered vector on all nodes .*/
47  if(broadcast || my_rank==0){
48  gathered_vector=xNew<double>(vector_size);
49  }
50 
51  /*Allocate local vectors*/
52  VecGetOwnershipRange(vector,&lower_row,&upper_row);
53  upper_row--;
54  range=upper_row-lower_row+1;
55 
56  if (range){
57  idxn=xNew<int>(range);
58  for (i=0;i<range;i++){
59  *(idxn+i)=lower_row+i;
60  }
61  local_vector=xNew<double>(range);
62  /*Extract values from MPI vector to serial local_vector on each node*/
63  VecGetValues(vector,range,idxn,local_vector);
64  }
65 
66  /*Now each node holds its local_vector containing range rows.
67  * We send this local_vector to the gathered_vector on node 0*/
68  for (i=1;i<num_procs;i++){
69  if (my_rank==i){
70  buffer[0]=my_rank;
71  buffer[1]=lower_row;
72  buffer[2]=range;
73  ISSM_MPI_Send(buffer,3,ISSM_MPI_INT,0,1,comm);
74  if (range)ISSM_MPI_Send(local_vector,range,ISSM_MPI_PDOUBLE,0,1,comm);
75  }
76  if (my_rank==0){
77  ISSM_MPI_Recv(buffer,3,ISSM_MPI_INT,i,1,comm,&status);
78  if (buffer[2])ISSM_MPI_Recv(gathered_vector+buffer[1],buffer[2],ISSM_MPI_PDOUBLE,i,1,comm,&status);
79  }
80  }
81 
82  if (my_rank==0){
83  //Still have the local_vector on node 0 to take care of.
84  if (range)memcpy(gathered_vector+lower_row,local_vector,range*sizeof(double));
85  }
86 
87  if(broadcast){
88  /*Now, broadcast gathered_vector from node 0 to other nodes: */
89  ISSM_MPI_Bcast(gathered_vector,vector_size,ISSM_MPI_PDOUBLE,0,comm);
90  }
91 
92  /*Assign output pointers: */
93  *pgathered_vector=gathered_vector;
94 
95  /*free ressources: */
96  xDelete<int>(idxn);
97  xDelete<double>(local_vector);
98 
99  return 1;
100 }

◆ MatFree()

void MatFree ( Mat *  pmat)

Definition at line 16 of file MatFree.cpp.

16  {
17 
18  #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
19  if(*pmat)MatDestroy(*pmat);
20  *pmat=NULL;
21  #else
22  if(*pmat)MatDestroy(pmat);
23  *pmat=NULL;
24  #endif
25 
26 }

◆ ISFree()

void ISFree ( IS *  pis)

◆ VecFree()

void VecFree ( Vec *  pvec)

Definition at line 16 of file VecFree.cpp.

16  {
17 
18  #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
19  if(*pvec)VecDestroy(*pvec);
20  #else
21  if(*pvec)VecDestroy(pvec);
22  #endif
23  *pvec=NULL;
24 
25 }

◆ KSPFree()

void KSPFree ( KSP *  pksp)

Definition at line 16 of file KSPFree.cpp.

16  {
17 
18  #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
19  if(*pksp)KSPDestroy(*pksp);
20  *pksp=NULL;
21  #else
22  if(*pksp)KSPDestroy(pksp);
23  *pksp=NULL;
24  #endif
25 
26 }

◆ MatPartition()

int MatPartition ( Mat *  poutmatrix,
Mat  matrixA,
double *  row_partition_vector,
int  row_partition_vector_size,
double *  col_partition_vector,
int  col_partition_vector_size 
)

◆ PetscOptionsDetermineSolverType()

void PetscOptionsDetermineSolverType ( int *  psolver_type)

Definition at line 20 of file PetscOptionsDetermineSolverType.cpp.

20  {
21 
22  char option[100];
23  #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
24  PetscTruth flag;
25  #else
26  PetscBool flag;
27  #endif
28 
29  /*output: */
30  int solver_type=PETSCPACKAGE;
31 
32  /*retrieve mat_type option: */
33  #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=7)
34  PetscOptionsGetString(NULL,PETSC_NULL,"-mat_type",&option[0],100,&flag);
35  #else
36  PetscOptionsGetString(PETSC_NULL,"-mat_type",&option[0],100,&flag);
37  #endif
38 
39  if (strcmp(option,"aijmumps")==0){
40  solver_type=MUMPSPACKAGE_LU;
41  }
42  if (strcmp(option,"sbaijmumps")==0){
43  solver_type=MUMPSPACKAGE_CHOL;
44  }
45  if (strcmp(option,"aijspooles")==0){
46  solver_type=SPOOLESPACKAGE_LU;
47  }
48  if (strcmp(option,"sbaijspooles")==0){
49  solver_type=SPOOLESPACKAGE_CHOL;
50  }
51  if (strcmp(option,"superlu_dist")==0){
52  solver_type=SUPERLUDISTPACKAGE;
53  }
54  if (strcmp(option,"")==0){
55  solver_type=SUPERLUDISTPACKAGE;
56  }
57 
58  #if _PETSC_MAJOR_ >= 3
59  #if (_PETSC_MINOR_>=7)
60  PetscOptionsGetString(NULL,PETSC_NULL,"-pc_factor_mat_solver_package",&option[0],100,&flag);
61  #else
62  PetscOptionsGetString(PETSC_NULL,"-pc_factor_mat_solver_package",&option[0],100,&flag);
63  #endif
64  if (strcmp(option,"mumps")==0){
65  solver_type=MUMPSPACKAGE_LU;
66  }
67  #endif
68 
69  #if (_PETSC_MAJOR_==3) && (_PETSC_MINOR_>=7)
70  PetscOptionsGetString(NULL,PETSC_NULL,"-issm_option_solver",&option[0],100,&flag);
71  #else
72  PetscOptionsGetString(PETSC_NULL,"-issm_option_solver",&option[0],100,&flag);
73  #endif
74  if(strcmp(option,"FS")==0 || strcmp(option,"stokes")==0){
75  solver_type=FSSolverEnum;
76  }
77 
78  *psolver_type=solver_type;
79 }

◆ MatMultPatch()

void MatMultPatch ( Mat  A,
Vec  X,
Vec  AX,
ISSM_MPI_Comm  comm 
)

Definition at line 23 of file MatMultPatch.cpp.

23  { //same prototype as MatMult in Petsc
24 
25  int m,n;
26  Vec X_rel=NULL;
27 
28  _assert_(A); _assert_(X);
29 
30  if (MatMultCompatible(A,X,comm)){
31  MatMult(A,X,AX);
32  }
33  else{
34  MatGetLocalSize(A,&m,&n);;
35  VecRelocalize(&X_rel,X,n,comm);
36  MatMult(A,X_rel,AX); ;
37  #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
38  VecDestroy(X_rel);
39  #else
40  VecDestroy(&X_rel);
41  #endif
42  }
43 }

◆ MatToSerial()

void MatToSerial ( double **  poutmatrix,
Mat  matrix,
ISSM_MPI_Comm  comm 
)

Definition at line 14 of file MatToSerial.cpp.

14  {
15 
16  int i;
17  int my_rank;
18  int num_procs;
19 
20  /*Petsc variables*/
21  PetscInt lower_row,upper_row;
22  int range;
23  int M,N; //size of matrix
24  ISSM_MPI_Status status;
25  int* idxm=NULL;
26  int* idxn=NULL;
27  double* local_matrix=NULL; /*matrix local to each node used for temporary holding matrix values*/
28  int buffer[3];
29 
30  /*recover my_rank and num_procs:*/
31  ISSM_MPI_Comm_size(comm,&num_procs);
32  ISSM_MPI_Comm_rank(comm,&my_rank);
33 
34  /*Output*/
35  double* outmatrix=NULL;
36 
37  /*get matrix size: */
38  MatGetSize(matrix,&M,&N);
39 
40  /*partition: */
41  MatGetOwnershipRange(matrix,&lower_row,&upper_row);
42  upper_row--;
43  range=upper_row-lower_row+1;
44 
45  /*Local and global allocation*/
46  if (my_rank==0)outmatrix=xNew<double>(M*N);
47 
48  if (range){
49  local_matrix=xNew<double>(N*range);
50  idxm=xNew<int>(range);
51  idxn=xNew<int>(N);
52 
53  for (i=0;i<N;i++){
54  *(idxn+i)=i;
55  }
56  for (i=0;i<range;i++){
57  *(idxm+i)=lower_row+i;
58  }
59 
60  MatGetValues(matrix,range,idxm,N,idxn,local_matrix);
61  }
62 
63  /*Now each node holds its local_matrix containing range rows.
64  * We send these rows to the matrix on node 0*/
65 
66  for (i=1;i<num_procs;i++){
67  if (my_rank==i){
68  buffer[0]=my_rank;
69  buffer[1]=lower_row;
70  buffer[2]=range;
71  ISSM_MPI_Send(buffer,3,ISSM_MPI_INT,0,1,comm);
72  if (range)ISSM_MPI_Send(local_matrix,N*range,ISSM_MPI_PDOUBLE,0,1,comm);
73  }
74  if (my_rank==0){
75  ISSM_MPI_Recv(buffer,3,ISSM_MPI_INT,i,1,comm,&status);
76  if (buffer[2])ISSM_MPI_Recv(outmatrix+(buffer[1]*N),N*buffer[2],ISSM_MPI_PDOUBLE,i,1,comm,&status);
77  }
78  }
79  if (my_rank==0){
80  //Still have the local_matrix on node 0 to take care of.
81  memcpy(outmatrix,local_matrix,N*range*sizeof(double));
82  }
83 
84  /*Assign output pointer: */
85  *poutmatrix=outmatrix;
86  xDelete<int>(idxm);
87  xDelete<int>(idxn);
88  xDelete<double>(local_matrix);
89 }

◆ SerialToVec()

Vec SerialToVec ( double *  vector,
int  vector_size 
)

◆ ISSMToPetscInsertMode()

InsertMode ISSMToPetscInsertMode ( InsMode  mode)

Definition at line 20 of file ISSMToPetscInsertMode.cpp.

20  {
21 
22  switch(mode){
23  case ADD_VAL:
24  return ADD_VALUES;
25  break;
26  case INS_VAL:
27  return INSERT_VALUES;
28  break;
29  default:
30  _error_("unknown insert mode!");
31  break;
32  }
33 }

◆ ISSMToPetscNormMode()

NormType ISSMToPetscNormMode ( NormMode  mode)

Definition at line 20 of file ISSMToPetscNormMode.cpp.

20  {
21 
22  switch(mode){
23  case NORM_INF:
24  return NORM_INFINITY;
25  break;
26  case NORM_TWO:
27  return NORM_2;
28  break;
29  default:
30  _error_("unknown norm !");
31  break;
32  }
33 }

◆ ISSMToPetscMatrixType()

MatType ISSMToPetscMatrixType ( MatrixType  type)

Definition at line 20 of file ISSMToPetscMatrixType.cpp.

20  {
21 
22  switch(type){
23  case DENSE_SEQUENTIAL:
24  return MATSEQDENSE;
25  break;
26  case SPARSE_SEQUENTIAL:
27  return MATSEQAIJ;
28  break;
29  default:
30  _error_("unknown matrix type !");
31  break;
32  }
33 }
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
SPOOLESPACKAGE_CHOL
@ SPOOLESPACKAGE_CHOL
Definition: SolverEnum.h:13
DENSE_SEQUENTIAL
@ DENSE_SEQUENTIAL
Definition: toolkitsenums.h:16
VecRelocalize
void VecRelocalize(Vec *outvector, Vec vector, int m, ISSM_MPI_Comm comm)
Definition: MatMultPatch.cpp:75
SUPERLUDISTPACKAGE
@ SUPERLUDISTPACKAGE
Definition: SolverEnum.h:14
ADD_VAL
@ ADD_VAL
Definition: toolkitsenums.h:14
SPOOLESPACKAGE_LU
@ SPOOLESPACKAGE_LU
Definition: SolverEnum.h:12
NORM_INF
@ NORM_INF
Definition: toolkitsenums.h:15
FSSolverEnum
@ FSSolverEnum
Definition: EnumDefinitions.h:1061
NORM_TWO
@ NORM_TWO
Definition: toolkitsenums.h:15
SPARSE_SEQUENTIAL
@ SPARSE_SEQUENTIAL
Definition: toolkitsenums.h:16
ISSM_MPI_Comm_rank
int ISSM_MPI_Comm_rank(ISSM_MPI_Comm comm, int *rank)
Definition: issmmpi.cpp:198
DetermineLocalSize
int DetermineLocalSize(int global_size, ISSM_MPI_Comm comm)
Definition: DetermineLocalSize.cpp:9
ISSM_MPI_INT
#define ISSM_MPI_INT
Definition: issmmpi.h:127
ISSM_MPI_Send
int ISSM_MPI_Send(void *buf, int count, ISSM_MPI_Datatype datatype, int dest, int tag, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:484
INS_VAL
@ INS_VAL
Definition: toolkitsenums.h:14
ISSM_MPI_Status
int ISSM_MPI_Status
Definition: issmmpi.h:121
ISSM_MPI_PDOUBLE
#define ISSM_MPI_PDOUBLE
Definition: issmmpi.h:126
ISSM_MPI_Bcast
int ISSM_MPI_Bcast(void *buffer, int count, ISSM_MPI_Datatype datatype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:162
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
ISSM_MPI_Recv
int ISSM_MPI_Recv(void *buf, int count, ISSM_MPI_Datatype datatype, int source, int tag, ISSM_MPI_Comm comm, ISSM_MPI_Status *status)
Definition: issmmpi.cpp:342
MatMultCompatible
int MatMultCompatible(Mat A, Vec x, ISSM_MPI_Comm comm)
Definition: MatMultPatch.cpp:45
MUMPSPACKAGE_CHOL
@ MUMPSPACKAGE_CHOL
Definition: SolverEnum.h:11
PETSCPACKAGE
@ PETSCPACKAGE
Definition: SolverEnum.h:9
ISSM_MPI_Comm_size
int ISSM_MPI_Comm_size(ISSM_MPI_Comm comm, int *size)
Definition: issmmpi.cpp:209
MUMPSPACKAGE_LU
@ MUMPSPACKAGE_LU
Definition: SolverEnum.h:10
Generated on Thu Jul 2 2020 08:09:25 for Ice Sheet System Model by   doxygen 1.8.19