source: issm/trunk-jpl/src/c/toolkits/petsc/patches/NewMat.cpp

Last change on this file was 27157, checked in by Mathieu Morlighem, 3 years ago

CHG: simplifying all petsc includes: only petscksp needs to be included now
File size: 2.8 KB
Line 
1/*!\file: NewMat.cpp
2 * \brief create matrix using the Petsc library
3 */
4
5#ifdef HAVE_CONFIG_H
6 #include <config.h>
7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
11/*Petsc includes: */
12#include <petscksp.h>
13
14#include "./petscpatches.h"
15#include "../../../shared/shared.h"
16#include "../../mpi/issmmpi.h"
17
18/*NewMat(int M,int N){{{*/
19Mat NewMat(int M,int N,ISSM_MPI_Comm comm){
20
21 /*output:*/
22 Mat outmatrix=NULL;
23
24 /*parameters: */
25 double sparsity=0.001; //default
26 int m,n;
27 int d_nz,o_nz,nnz;
28
29 /*Determine local sizes: */
30 m=DetermineLocalSize(M,comm);
31 n=DetermineLocalSize(N,comm);
32
33 nnz=(int)((double)M*(double)N*sparsity); //number of non zeros.
34 d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
35 o_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
36
37 #if PETSC_VERSION_GT(3,2,0)
38 MatCreateAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
39 #else
40 MatCreateMPIAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
41 #endif
42
43 return outmatrix;
44}
45/*}}}*/
46/*NewMat(int M,int N,double sparsity,ISSM_MPI_Comm comm){{{*/
47Mat NewMat(int M,int N,double sparsity,ISSM_MPI_Comm comm){
48
49 /*output:*/
50 Mat outmatrix=NULL;
51
52 /*parameters: */
53 int m,n;
54 int d_nz,o_nz;
55 int nnz;
56
57 /*Determine local sizes: */
58 m=DetermineLocalSize(M,comm);
59 n=DetermineLocalSize(N,comm);
60
61 nnz=(int)((double)M*(double)N*sparsity); //number of non zeros.
62 d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
63 o_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
64
65 #if PETSC_VERSION_GT(3,2,0)
66 if(sparsity==1){
67 MatCreateDense(comm,m,n,M,N,NULL,&outmatrix);
68 }
69 else{
70 MatCreateAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
71 }
72 #else
73 MatCreateMPIAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
74 #endif
75
76 return outmatrix;
77}
78/*}}}*/
79/*NewMat(int M,int N,int connectivity,int numberofdofspernode){{{*/
80Mat NewMat(int M,int N,int connectivity,int numberofdofspernode,ISSM_MPI_Comm comm){
81
82 /*output:*/
83 Mat outmatrix=NULL;
84
85 /*parameters: */
86 int m,n;
87 int d_nz,o_nz;
88
89 #if PETSC_VERSION_MAJOR >= 3
90 #if defined(_HAVE_PETSCDEV_) || PETSC_VERSION_MINOR >=4
91 MatType type;
92 #else
93 const MatType type;
94 #endif
95 #else
96 MatType type;
97 #endif
98
99 /*Determine local sizes: */
100 m=DetermineLocalSize(M,comm);
101 n=DetermineLocalSize(N,comm);
102
103 /*Figure out number of non zeros per row: */
104 d_nz=(int)connectivity*numberofdofspernode/2;
105 o_nz=(int)connectivity*numberofdofspernode/2;
106
107 MatCreate(comm,&outmatrix);
108 MatSetSizes(outmatrix,m,n,M,N);
109 MatSetFromOptions(outmatrix);
110 MatSetOption(outmatrix,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
111
112 /*preallocation according to type: */
113 MatGetType(outmatrix,&type);
114
115 if((strcmp(type,"mpiaij")==0) || (strcmp(type,"mpidense")==0)){
116 MatMPIAIJSetPreallocation(outmatrix,d_nz,NULL,o_nz,NULL);
117 }
118
119 return outmatrix;
120}
121/*}}}*/
Note: See TracBrowser for help on using the repository browser.