source: issm/branches/trunk-jpl-damage/src/c/toolkits/petsc/patches/MatPartition.cpp@ 11427

Last change on this file since 11427 was 11197, checked in by Mathieu Morlighem, 14 years ago

Fixed many warnings on pleiades: missing return statement, uninitialized
values, wrong print argument,... etc
File size: 3.7 KB
Line 
1/*!\file: MatPartition.cpp
2 * \brief partition matrix according to node sets
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 "../../../shared/shared.h"
17
18int MatPartition(Mat* poutmatrix,Mat matrixA,double* row_partition_vector,int row_partition_vector_size ,
19 double* col_partition_vector,int col_partition_vector_size){
20
21 int i;
22
23 /*Petsc matrix*/
24 int d_nz;
25 int o_nz;
26 int* node_rows=NULL;
27 int* node_cols=NULL;
28 int count;
29 IS col_index=NULL;
30 IS row_index=NULL;
31 int lower_row,upper_row,range;
32 int MA,NA; //matrixA dimensions
33 const char* type=NULL;
34 int csize;
35
36 /*output*/
37 Mat outmatrix=NULL;
38
39 /*get input matrix size: */
40 MatGetSize(matrixA,&MA,&NA);
41
42 /*If one of the partitioning row or col vectors has a 0 dimension, we return a NILL matrix. Return it with the same type as matrixA.*/
43
44 if ((row_partition_vector_size==0) || (col_partition_vector==0)){
45 MatGetType(matrixA,&type);
46 if (strcmp(type,"mpiaij")==0){
47 d_nz=0;
48 o_nz=0;
49 MatCreateMPIAIJ(MPI_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE, 0,0,d_nz,PETSC_NULL,o_nz,PETSC_NULL,&outmatrix);
50 }
51 else if (strcmp(type,"mpidense")==0){
52 MatCreateMPIDense(MPI_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE, 0,0,PETSC_NULL,&outmatrix);
53 }
54 /*Assemble*/
55 MatAssemblyBegin(outmatrix,MAT_FINAL_ASSEMBLY);
56 MatAssemblyEnd(outmatrix,MAT_FINAL_ASSEMBLY);
57 }
58 else{
59 /*Both vectors are non nill, use MatGetSubMatrix to condense out*/
60 /*Figure out which rows each node is going to get from matrix A.*/
61 MatGetOwnershipRange(matrixA,&lower_row,&upper_row);
62 upper_row--;
63 range=upper_row-lower_row+1;
64
65 count=0;
66 if (range){
67 node_rows=(int*)xmalloc(range*sizeof(int)); //this is the maximum number of rows one node can extract.
68
69 for (i=0;i<row_partition_vector_size;i++){
70 if ( ((int)(*(row_partition_vector+i))>=(lower_row+1)) && ((int)(*(row_partition_vector+i))<=(upper_row+1)) ){
71 *(node_rows+count)=(int)*(row_partition_vector+i)-1;
72 count++;
73 }
74 }
75 }
76 else{
77 count=0;
78 }
79
80 /*Now each node has a node_rows vectors holding which rows they should extract from matrixA. Create an Index Set from node_rows.*/
81 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
82 ISCreateGeneral(MPI_COMM_WORLD,count,node_rows,&row_index);
83 #else
84 ISCreateGeneral(MPI_COMM_WORLD,count,node_rows,PETSC_COPY_VALUES,&row_index);
85 #endif
86
87 /*Same deal for columns*/
88 node_cols=(int*)xmalloc(col_partition_vector_size*sizeof(int));
89 for (i=0;i<col_partition_vector_size;i++){
90 *(node_cols+i)=(int)*(col_partition_vector+i)-1;
91 }
92 #if _PETSC_MAJOR_ < 3 || (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ < 2)
93 ISCreateGeneral(MPI_COMM_WORLD,col_partition_vector_size,node_cols,&col_index);
94 #else
95 ISCreateGeneral(MPI_COMM_WORLD,col_partition_vector_size,node_cols,PETSC_COPY_VALUES,&col_index);
96 #endif
97
98 /*Call MatGetSubMatrix*/
99 csize=PetscDetermineLocalSize(col_partition_vector_size);
100 if(col_partition_vector_size==row_partition_vector_size){
101 #if _PETSC_MAJOR_ >= 3
102 MatGetSubMatrix(matrixA,row_index,col_index,MAT_INITIAL_MATRIX,&outmatrix);
103 #else
104 MatGetSubMatrix(matrixA,row_index,col_index,count,MAT_INITIAL_MATRIX,&outmatrix);
105 #endif
106 }
107 else{
108 #if _PETSC_MAJOR_ >= 3
109 MatGetSubMatrix(matrixA,row_index,col_index,MAT_INITIAL_MATRIX,&outmatrix);
110 #else
111 MatGetSubMatrix(matrixA,row_index,col_index,csize,MAT_INITIAL_MATRIX,&outmatrix);
112 #endif
113 }
114
115 }
116
117 /*Free ressources:*/
118 xfree((void**)&node_rows);
119 xfree((void**)&node_cols);
120 ISFree(&col_index);
121 ISFree(&row_index);
122
123 /*Assign output pointers:*/
124 *poutmatrix=outmatrix;
125 return 1;
126}
Note: See TracBrowser for help on using the repository browser.