1 2 #include "src/mat/impls/aij/mpi/mpiaij.h" 3 4 EXTERN PetscErrorCode CreateColmap_MPIAIJ_Private(Mat); 5 EXTERN PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat,PetscInt,PetscTruth,PetscInt*,PetscInt*[],PetscInt*[],PetscTruth*); 6 EXTERN PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat,PetscInt,PetscTruth,PetscInt*,PetscInt*[],PetscInt*[],PetscTruth*); 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "MatFDColoringCreate_MPIAIJ" 10 PetscErrorCode MatFDColoringCreate_MPIAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 11 { 12 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)mat->data; 13 PetscErrorCode ierr; 14 PetscMPIInt size,*ncolsonproc,*disp,nn; 15 PetscInt i,*is,n,nrows,j,k,m,*rows = 0,*A_ci,*A_cj,ncols,col; 16 PetscInt nis = iscoloring->n,nctot,*cols,*B_ci,*B_cj; 17 PetscInt *rowhit,M = mat->m,cstart = aij->cstart,cend = aij->cend,colb; 18 PetscInt *columnsforrow,l; 19 IS *isa; 20 PetscTruth done,flg; 21 22 PetscFunctionBegin; 23 if (!mat->assembled) { 24 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled first; MatAssemblyBegin/End();"); 25 } 26 27 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 28 c->M = mat->M; /* set the global rows and columns and local rows */ 29 c->N = mat->N; 30 c->m = mat->m; 31 c->rstart = aij->rstart; 32 33 c->ncolors = nis; 34 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 35 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 36 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 37 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 38 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 39 PetscLogObjectMemory(c,5*nis*sizeof(PetscInt)); 40 41 /* Allow access to data structures of local part of matrix */ 42 if (!aij->colmap) { 43 ierr = CreateColmap_MPIAIJ_Private(mat);CHKERRQ(ierr); 44 } 45 /* 46 Calls the _SeqAIJ() version of these routines to make sure it does not 47 get the reduced (by inodes) version of I and J 48 */ 49 ierr = MatGetColumnIJ_SeqAIJ(aij->A,0,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr); 50 ierr = MatGetColumnIJ_SeqAIJ(aij->B,0,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr); 51 52 ierr = MPI_Comm_size(mat->comm,&size);CHKERRQ(ierr); 53 ierr = PetscMalloc(2*size*sizeof(PetscInt*),&ncolsonproc);CHKERRQ(ierr); 54 disp = ncolsonproc + size; 55 56 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 57 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); 58 59 /* 60 Temporary option to allow for debugging/testing 61 */ 62 ierr = PetscOptionsHasName(PETSC_NULL,"-matfdcoloring_slow",&flg);CHKERRQ(ierr); 63 64 for (i=0; i<nis; i++) { 65 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 66 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 67 c->ncolumns[i] = n; 68 c->ncolumns[i] = n; 69 if (n) { 70 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 71 PetscLogObjectMemory(c,n*sizeof(PetscInt)); 72 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 73 } else { 74 c->columns[i] = 0; 75 } 76 77 /* Determine the total (parallel) number of columns of this color */ 78 nn = (PetscMPIInt)n; 79 ierr = MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,mat->comm);CHKERRQ(ierr); 80 nctot = 0; for (j=0; j<size; j++) {nctot += ncolsonproc[j];} 81 if (!nctot) { 82 PetscLogInfo((PetscObject)mat,"MatFDColoringCreate_MPIAIJ: Coloring of matrix has some unneeded colors with no corresponding rows\n"); 83 } 84 85 disp[0] = 0; 86 for (j=1; j<size; j++) { 87 disp[j] = disp[j-1] + ncolsonproc[j-1]; 88 } 89 90 /* Get complete list of columns for color on each processor */ 91 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 92 ierr = MPI_Allgatherv(is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,mat->comm);CHKERRQ(ierr); 93 94 /* 95 Mark all rows affect by these columns 96 */ 97 if (!flg) {/*-----------------------------------------------------------------------------*/ 98 /* crude, fast version */ 99 ierr = PetscMemzero(rowhit,M*sizeof(PetscInt));CHKERRQ(ierr); 100 /* loop over columns*/ 101 for (j=0; j<nctot; j++) { 102 col = cols[j]; 103 if (col >= cstart && col < cend) { 104 /* column is in diagonal block of matrix */ 105 rows = A_cj + A_ci[col-cstart]; 106 m = A_ci[col-cstart+1] - A_ci[col-cstart]; 107 } else { 108 #if defined (PETSC_USE_CTABLE) 109 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr) 110 colb --; 111 #else 112 colb = aij->colmap[col] - 1; 113 #endif 114 if (colb == -1) { 115 m = 0; 116 } else { 117 rows = B_cj + B_ci[colb]; 118 m = B_ci[colb+1] - B_ci[colb]; 119 } 120 } 121 /* loop over columns marking them in rowhit */ 122 for (k=0; k<m; k++) { 123 rowhit[*rows++] = col + 1; 124 } 125 } 126 127 /* count the number of hits */ 128 nrows = 0; 129 for (j=0; j<M; j++) { 130 if (rowhit[j]) nrows++; 131 } 132 c->nrows[i] = nrows; 133 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 134 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 135 PetscLogObjectMemory(c,2*(nrows+1)*sizeof(PetscInt)); 136 nrows = 0; 137 for (j=0; j<M; j++) { 138 if (rowhit[j]) { 139 c->rows[i][nrows] = j; 140 c->columnsforrow[i][nrows] = rowhit[j] - 1; 141 nrows++; 142 } 143 } 144 } else {/*-------------------------------------------------------------------------------*/ 145 /* slow version, using rowhit as a linked list */ 146 PetscInt currentcol,fm,mfm; 147 rowhit[M] = M; 148 nrows = 0; 149 /* loop over columns*/ 150 for (j=0; j<nctot; j++) { 151 col = cols[j]; 152 if (col >= cstart && col < cend) { 153 /* column is in diagonal block of matrix */ 154 rows = A_cj + A_ci[col-cstart]; 155 m = A_ci[col-cstart+1] - A_ci[col-cstart]; 156 } else { 157 #if defined (PETSC_USE_CTABLE) 158 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 159 colb --; 160 #else 161 colb = aij->colmap[col] - 1; 162 #endif 163 if (colb == -1) { 164 m = 0; 165 } else { 166 rows = B_cj + B_ci[colb]; 167 m = B_ci[colb+1] - B_ci[colb]; 168 } 169 } 170 /* loop over columns marking them in rowhit */ 171 fm = M; /* fm points to first entry in linked list */ 172 for (k=0; k<m; k++) { 173 currentcol = *rows++; 174 /* is it already in the list? */ 175 do { 176 mfm = fm; 177 fm = rowhit[fm]; 178 } while (fm < currentcol); 179 /* not in list so add it */ 180 if (fm != currentcol) { 181 nrows++; 182 columnsforrow[currentcol] = col; 183 /* next three lines insert new entry into linked list */ 184 rowhit[mfm] = currentcol; 185 rowhit[currentcol] = fm; 186 fm = currentcol; 187 /* fm points to present position in list since we know the columns are sorted */ 188 } else { 189 SETERRQ(PETSC_ERR_PLIB,"Invalid coloring of matrix detected"); 190 } 191 } 192 } 193 c->nrows[i] = nrows; 194 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 195 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 196 PetscLogObjectMemory(c,(nrows+1)*sizeof(PetscInt)); 197 /* now store the linked list of rows into c->rows[i] */ 198 nrows = 0; 199 fm = rowhit[M]; 200 do { 201 c->rows[i][nrows] = fm; 202 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 203 fm = rowhit[fm]; 204 } while (fm < M); 205 } /* ---------------------------------------------------------------------------------------*/ 206 ierr = PetscFree(cols);CHKERRQ(ierr); 207 } 208 209 /* Optimize by adding the vscale, and scaleforrow[][] fields */ 210 /* 211 vscale will contain the "diagonal" on processor scalings followed by the off processor 212 */ 213 ierr = VecCreateGhost(mat->comm,aij->A->m,PETSC_DETERMINE,aij->B->n,aij->garray,&c->vscale);CHKERRQ(ierr) 214 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 215 for (k=0; k<c->ncolors; k++) { 216 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 217 for (l=0; l<c->nrows[k]; l++) { 218 col = c->columnsforrow[k][l]; 219 if (col >= cstart && col < cend) { 220 /* column is in diagonal block of matrix */ 221 colb = col - cstart; 222 } else { 223 /* column is in "off-processor" part */ 224 #if defined (PETSC_USE_CTABLE) 225 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 226 colb --; 227 #else 228 colb = aij->colmap[col] - 1; 229 #endif 230 colb += cend - cstart; 231 } 232 c->vscaleforrow[k][l] = colb; 233 } 234 } 235 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 236 237 ierr = PetscFree(rowhit);CHKERRQ(ierr); 238 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 239 ierr = PetscFree(ncolsonproc);CHKERRQ(ierr); 240 ierr = MatRestoreColumnIJ_SeqAIJ(aij->A,0,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr); 241 ierr = MatRestoreColumnIJ_SeqAIJ(aij->B,0,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr); 242 PetscFunctionReturn(0); 243 } 244 245 246 247 248 249 250