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

Last change on this file since 13596 was 13596, checked in by Eric.Larour, 12 years ago

CHG: fixed segfault in NewMat call
File size: 2.9 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 "petscmat.h"
13#include "petscvec.h"
14#include "petscksp.h"
15
16#include "./petscpatches.h"
17
18#include "../../../shared/shared.h"
19#include "../../../include/include.h"
20#include "../../mpi/patches/mpipatches.h"
21
22/*NewMat(int M,int N){{{*/
23Mat NewMat(int M,int N,COMM comm){
24
25 /*output:*/
26 Mat outmatrix=NULL;
27
28 /*parameters: */
29 double sparsity=.001; //default
30 int m,n;
31 int d_nz,o_nz,nnz;
32
33 /*Determine local sizes: */
34 m=DetermineLocalSize(M);
35 n=DetermineLocalSize(N);
36
37 nnz=(int)((double)M*(double)N*sparsity); //number of non zeros.
38 d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
39 o_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
40
41 #if _PETSC_MAJOR_ == 3 && _PETSC_MINOR_ > 2
42 MatCreateAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
43 #else
44 MatCreateMPIAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
45 #endif
46
47 return outmatrix;
48}
49/*}}}*/
50/*NewMat(int M,int N,double sparsity,COMM comm){{{*/
51Mat NewMat(int M,int N,double sparsity,COMM comm){
52
53 /*output:*/
54 Mat outmatrix=NULL;
55
56 /*parameters: */
57 int m,n;
58 int d_nz,o_nz;
59 int nnz;
60
61 /*Determine local sizes: */
62 m=DetermineLocalSize(M);
63 n=DetermineLocalSize(N);
64
65 nnz=(int)((double)M*(double)N*sparsity); //number of non zeros.
66 d_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
67 o_nz=(int)((double)nnz/(double)M/2.0); //number of non zeros per row/2
68
69 #if _PETSC_MAJOR_ == 3 && _PETSC_MINOR_ > 2
70 if(sparsity==1){
71 MatCreateDense(comm,m,n,M,N,NULL,&outmatrix);
72 }
73 else{
74 MatCreateAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
75 }
76 #else
77 MatCreateMPIAIJ(comm,m,n,M,N,d_nz,NULL,o_nz,NULL,&outmatrix);
78 #endif
79
80 return outmatrix;
81}
82/*}}}*/
83/*NewMat(int M,int N,int connectivity,int numberofdofspernode){{{*/
84Mat NewMat(int M,int N,int connectivity,int numberofdofspernode,COMM comm){
85
86 /*output:*/
87 Mat outmatrix=NULL;
88
89 /*parameters: */
90 int m,n;
91 int d_nz,o_nz;
92 int nnz;
93
94 #if _PETSC_MAJOR_ >= 3
95 const MatType type;
96 #else
97 MatType type;
98 #endif
99
100 /*Determine local sizes: */
101 m=DetermineLocalSize(M);
102 n=DetermineLocalSize(N);
103
104 /*Figure out number of non zeros per row: */
105 d_nz=(int)connectivity*numberofdofspernode/2;
106 o_nz=(int)connectivity*numberofdofspernode/2;
107
108 MatCreate(comm,&outmatrix);
109 MatSetSizes(outmatrix,m,n,M,N);
110 MatSetFromOptions(outmatrix);
111
112 /*preallocation according to type: */
113 MatGetType(outmatrix,&type);
114
115 #if _PETSC_MAJOR_ == 2
116 if((strcmp(type,"mpiaij")==0) || (strcmp(type,"aijmumps")==0)){
117 MatMPIAIJSetPreallocation(outmatrix,d_nz,NULL,o_nz,NULL);
118 }
119 #else
120 if((strcmp(type,"mpiaij")==0) || (strcmp(type,"mpidense")==0)){
121 MatMPIAIJSetPreallocation(outmatrix,d_nz,NULL,o_nz,NULL);
122 }
123 #endif
124
125 return outmatrix;
126}
127/*}}}*/
Note: See TracBrowser for help on using the repository browser.