1 #define PETSCMAT_DLL 2 3 #include "../src/mat/impls/aij/mpi/mpiaij.h" 4 5 EXTERN PetscErrorCode CreateColmap_MPIAIJ_Private(Mat); 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "MatFDColoringCreate_MPIAIJ" 9 /* 10 This routine is shared by AIJ and BAIJ matrices, since it operators only on the nonzero structure of the elements or blocks. 11 This is why it has the ugly code with the MatGetBlockSize() 12 */ 13 PetscErrorCode MatFDColoringCreate_MPIAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 14 { 15 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)mat->data; 16 PetscErrorCode ierr; 17 PetscMPIInt size,*ncolsonproc,*disp,nn; 18 PetscInt bs = 1,i,n,nrows,j,k,m,*rows = 0,*A_ci,*A_cj,ncols,col; 19 const PetscInt *is; 20 PetscInt nis = iscoloring->n,nctot,*cols,*B_ci,*B_cj; 21 PetscInt *rowhit,M = mat->rmap->n,cstart = mat->cmap->rstart,cend = mat->cmap->rend,colb; 22 PetscInt *columnsforrow,l; 23 IS *isa; 24 PetscTruth done,flg; 25 ISLocalToGlobalMapping map = mat->mapping; 26 PetscInt *ltog = (map ? map->indices : (PetscInt*) PETSC_NULL) ,ctype=c->ctype; 27 PetscTruth flg1,flg2; 28 29 PetscFunctionBegin; 30 if (!mat->assembled) { 31 SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled first; MatAssemblyBegin/End();"); 32 } 33 if (ctype == IS_COLORING_GHOSTED && !map) SETERRQ(PETSC_ERR_ARG_INCOMP,"When using ghosted differencing matrix must have local to global mapping provided with MatSetLocalToGlobalMapping"); 34 35 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 36 37 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */ 38 ierr = PetscTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);CHKERRQ(ierr); 39 ierr = PetscTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);CHKERRQ(ierr); 40 if (flg1 || flg2) { 41 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 42 } 43 44 c->M = mat->rmap->N/bs; /* set the global rows and columns and local rows */ 45 c->N = mat->cmap->N/bs; 46 c->m = mat->rmap->n/bs; 47 c->rstart = mat->rmap->rstart/bs; 48 49 c->ncolors = nis; 50 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 51 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 52 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 53 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 54 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 55 ierr = PetscLogObjectMemory(c,5*nis*sizeof(PetscInt));CHKERRQ(ierr); 56 57 /* Allow access to data structures of local part of matrix */ 58 if (!aij->colmap) { 59 ierr = CreateColmap_MPIAIJ_Private(mat);CHKERRQ(ierr); 60 } 61 ierr = MatGetColumnIJ(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr); 62 ierr = MatGetColumnIJ(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr); 63 64 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 65 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); 66 67 for (i=0; i<nis; i++) { 68 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 69 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 70 c->ncolumns[i] = n; 71 if (n) { 72 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 73 ierr = PetscLogObjectMemory(c,n*sizeof(PetscInt));CHKERRQ(ierr); 74 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 75 } else { 76 c->columns[i] = 0; 77 } 78 79 if (ctype == IS_COLORING_GLOBAL){ 80 /* Determine the total (parallel) number of columns of this color */ 81 ierr = MPI_Comm_size(((PetscObject)mat)->comm,&size);CHKERRQ(ierr); 82 ierr = PetscMalloc(2*size*sizeof(PetscInt*),&ncolsonproc);CHKERRQ(ierr); 83 disp = ncolsonproc + size; 84 85 nn = PetscMPIIntCast(n); 86 ierr = MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,((PetscObject)mat)->comm);CHKERRQ(ierr); 87 nctot = 0; for (j=0; j<size; j++) {nctot += ncolsonproc[j];} 88 if (!nctot) { 89 ierr = PetscInfo(mat,"Coloring of matrix has some unneeded colors with no corresponding rows\n");CHKERRQ(ierr); 90 } 91 92 disp[0] = 0; 93 for (j=1; j<size; j++) { 94 disp[j] = disp[j-1] + ncolsonproc[j-1]; 95 } 96 97 /* Get complete list of columns for color on each processor */ 98 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 99 ierr = MPI_Allgatherv((void*)is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,((PetscObject)mat)->comm);CHKERRQ(ierr); 100 ierr = PetscFree(ncolsonproc);CHKERRQ(ierr); 101 } else if (ctype == IS_COLORING_GHOSTED){ 102 /* Determine local number of columns of this color on this process, including ghost points */ 103 nctot = n; 104 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 105 ierr = PetscMemcpy(cols,is,n*sizeof(PetscInt));CHKERRQ(ierr); 106 } else { 107 SETERRQ(PETSC_ERR_SUP,"Not provided for this MatFDColoring type"); 108 } 109 110 /* 111 Mark all rows affect by these columns 112 */ 113 /* Temporary option to allow for debugging/testing */ 114 flg = PETSC_FALSE; 115 ierr = PetscOptionsGetTruth(PETSC_NULL,"-matfdcoloring_slow",&flg,PETSC_NULL);CHKERRQ(ierr); 116 if (!flg) {/*-----------------------------------------------------------------------------*/ 117 /* crude, fast version */ 118 ierr = PetscMemzero(rowhit,M*sizeof(PetscInt));CHKERRQ(ierr); 119 /* loop over columns*/ 120 for (j=0; j<nctot; j++) { 121 if (ctype == IS_COLORING_GHOSTED) { 122 col = ltog[cols[j]]; 123 } else { 124 col = cols[j]; 125 } 126 if (col >= cstart && col < cend) { 127 /* column is in diagonal block of matrix */ 128 rows = A_cj + A_ci[col-cstart]; 129 m = A_ci[col-cstart+1] - A_ci[col-cstart]; 130 } else { 131 #if defined (PETSC_USE_CTABLE) 132 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr) 133 colb --; 134 #else 135 colb = aij->colmap[col] - 1; 136 #endif 137 if (colb == -1) { 138 m = 0; 139 } else { 140 rows = B_cj + B_ci[colb]; 141 m = B_ci[colb+1] - B_ci[colb]; 142 } 143 } 144 /* loop over columns marking them in rowhit */ 145 for (k=0; k<m; k++) { 146 rowhit[*rows++] = col + 1; 147 } 148 } 149 150 /* count the number of hits */ 151 nrows = 0; 152 for (j=0; j<M; j++) { 153 if (rowhit[j]) nrows++; 154 } 155 c->nrows[i] = nrows; 156 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 157 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 158 ierr = PetscLogObjectMemory(c,2*(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr); 159 nrows = 0; 160 for (j=0; j<M; j++) { 161 if (rowhit[j]) { 162 c->rows[i][nrows] = j; 163 c->columnsforrow[i][nrows] = rowhit[j] - 1; 164 nrows++; 165 } 166 } 167 } else {/*-------------------------------------------------------------------------------*/ 168 /* slow version, using rowhit as a linked list */ 169 PetscInt currentcol,fm,mfm; 170 rowhit[M] = M; 171 nrows = 0; 172 /* loop over columns*/ 173 for (j=0; j<nctot; j++) { 174 if (ctype == IS_COLORING_GHOSTED) { 175 col = ltog[cols[j]]; 176 } else { 177 col = cols[j]; 178 } 179 if (col >= cstart && col < cend) { 180 /* column is in diagonal block of matrix */ 181 rows = A_cj + A_ci[col-cstart]; 182 m = A_ci[col-cstart+1] - A_ci[col-cstart]; 183 } else { 184 #if defined (PETSC_USE_CTABLE) 185 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 186 colb --; 187 #else 188 colb = aij->colmap[col] - 1; 189 #endif 190 if (colb == -1) { 191 m = 0; 192 } else { 193 rows = B_cj + B_ci[colb]; 194 m = B_ci[colb+1] - B_ci[colb]; 195 } 196 } 197 198 /* loop over columns marking them in rowhit */ 199 fm = M; /* fm points to first entry in linked list */ 200 for (k=0; k<m; k++) { 201 currentcol = *rows++; 202 /* is it already in the list? */ 203 do { 204 mfm = fm; 205 fm = rowhit[fm]; 206 } while (fm < currentcol); 207 /* not in list so add it */ 208 if (fm != currentcol) { 209 nrows++; 210 columnsforrow[currentcol] = col; 211 /* next three lines insert new entry into linked list */ 212 rowhit[mfm] = currentcol; 213 rowhit[currentcol] = fm; 214 fm = currentcol; 215 /* fm points to present position in list since we know the columns are sorted */ 216 } else { 217 SETERRQ(PETSC_ERR_PLIB,"Invalid coloring of matrix detected"); 218 } 219 } 220 } 221 c->nrows[i] = nrows; 222 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 223 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 224 ierr = PetscLogObjectMemory(c,(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr); 225 /* now store the linked list of rows into c->rows[i] */ 226 nrows = 0; 227 fm = rowhit[M]; 228 do { 229 c->rows[i][nrows] = fm; 230 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 231 fm = rowhit[fm]; 232 } while (fm < M); 233 } /* ---------------------------------------------------------------------------------------*/ 234 ierr = PetscFree(cols);CHKERRQ(ierr); 235 } 236 237 /* Optimize by adding the vscale, and scaleforrow[][] fields */ 238 /* 239 vscale will contain the "diagonal" on processor scalings followed by the off processor 240 */ 241 if (ctype == IS_COLORING_GLOBAL) { 242 ierr = VecCreateGhost(((PetscObject)mat)->comm,aij->A->rmap->n,PETSC_DETERMINE,aij->B->cmap->n,aij->garray,&c->vscale);CHKERRQ(ierr); 243 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 244 for (k=0; k<c->ncolors; k++) { 245 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 246 for (l=0; l<c->nrows[k]; l++) { 247 col = c->columnsforrow[k][l]; 248 if (col >= cstart && col < cend) { 249 /* column is in diagonal block of matrix */ 250 colb = col - cstart; 251 } else { 252 /* column is in "off-processor" part */ 253 #if defined (PETSC_USE_CTABLE) 254 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 255 colb --; 256 #else 257 colb = aij->colmap[col] - 1; 258 #endif 259 colb += cend - cstart; 260 } 261 c->vscaleforrow[k][l] = colb; 262 } 263 } 264 } else if (ctype == IS_COLORING_GHOSTED) { 265 /* Get gtol mapping */ 266 PetscInt N = mat->cmap->N, *gtol; 267 ierr = PetscMalloc((N+1)*sizeof(PetscInt),>ol);CHKERRQ(ierr); 268 for (i=0; i<N; i++) gtol[i] = -1; 269 for (i=0; i<map->n; i++) gtol[ltog[i]] = i; 270 271 c->vscale = 0; /* will be created in MatFDColoringApply() */ 272 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 273 for (k=0; k<c->ncolors; k++) { 274 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 275 for (l=0; l<c->nrows[k]; l++) { 276 col = c->columnsforrow[k][l]; /* global column index */ 277 c->vscaleforrow[k][l] = gtol[col]; /* local column index */ 278 } 279 } 280 ierr = PetscFree(gtol);CHKERRQ(ierr); 281 } 282 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 283 284 ierr = PetscFree(rowhit);CHKERRQ(ierr); 285 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 286 ierr = MatRestoreColumnIJ(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr); 287 ierr = MatRestoreColumnIJ(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290 291 292 293 294 295 296