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

Last change on this file since 11708 was 11708, checked in by cborstad, 13 years ago

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