1 2 #include <../src/mat/impls/aij/mpi/mpiaij.h> 3 4 extern PetscErrorCode MatCreateColmap_MPIAIJ_Private(Mat); 5 6 #undef __FUNCT__ 7 #define __FUNCT__ "MatFDColoringApply_MPIAIJ" 8 PetscErrorCode MatFDColoringApply_MPIAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 9 { 10 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void*))coloring->f; 11 PetscErrorCode ierr; 12 PetscInt k,start,end,l,row,col,colb,**columnsforrow=coloring->columnsforrow,nz; 13 PetscScalar dx,*y,*xx,*w3_array; 14 PetscScalar *vscale_array; 15 PetscReal epsilon = coloring->error_rel,umin = coloring->umin,unorm; 16 Vec w1 = coloring->w1,w2=coloring->w2,w3; 17 void *fctx = coloring->fctx; 18 PetscBool flg = PETSC_FALSE; 19 PetscInt ctype = coloring->ctype,nxloc; 20 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)J->data; 21 Mat A=aij->A,B=aij->B; 22 Mat_SeqAIJ *cspA=(Mat_SeqAIJ*)A->data, *cspB=(Mat_SeqAIJ*)B->data; 23 PetscMPIInt rank; 24 MatEntry *Jentry=coloring->matentry; 25 26 PetscFunctionBegin; 27 ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)J), &rank);CHKERRQ(ierr); 28 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 29 ierr = PetscOptionsGetBool(NULL,"-mat_fd_coloring_dont_rezero",&flg,NULL);CHKERRQ(ierr); 30 if (flg) { 31 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 32 } else { 33 PetscBool assembled; 34 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 35 if (assembled) { 36 ierr = MatZeroEntries(J);CHKERRQ(ierr); 37 } 38 } 39 40 if (!coloring->vscale) { //ctype == IS_COLORING_GHOSTED ? 41 printf("[%d] create a non-ghostedcoloring->vscale\n",rank); 42 ierr = VecDuplicate(x1,&coloring->vscale);CHKERRQ(ierr); 43 //SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_NULL,"Null Object: coloring->scale"); 44 } 45 46 /* Set w1 = F(x1) */ 47 if (!coloring->fset) { 48 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 49 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 50 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 51 } else { 52 coloring->fset = PETSC_FALSE; 53 } 54 55 if (!coloring->w3) { 56 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 57 ierr = PetscLogObjectParent((PetscObject)coloring,(PetscObject)coloring->w3);CHKERRQ(ierr); 58 } 59 w3 = coloring->w3; 60 61 /* Compute vscale = 1./dx - the local scale factors, including ghost points */ 62 ierr = VecGetOwnershipRange(w1,&start,&end);CHKERRQ(ierr); /* =J's row ownershipRange, used by ghosted x! */ 63 //printf("[%d] start end %d %d\n",rank,start,end); 64 65 if (coloring->htype[0] == 'w') { 66 /* tacky test; need to make systematic if we add other approaches to computing h*/ 67 ierr = VecNorm(x1,NORM_2,&unorm);CHKERRQ(ierr); 68 dx = 1.0/((1.0 + unorm)*epsilon); 69 ierr = VecSet(coloring->vscale,dx);CHKERRQ(ierr); 70 } else { 71 ierr = VecGetLocalSize(x1,&nxloc);CHKERRQ(ierr); 72 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr); 73 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 74 for (col=0; col<nxloc; col++) { 75 dx = xx[col]; 76 if (PetscAbsScalar(dx) < umin) { 77 if (PetscRealPart(dx) >= 0.0) dx = umin; 78 else if (PetscRealPart(dx) < 0.0 ) dx = -umin; 79 } 80 dx *= epsilon; 81 vscale_array[col] = 1.0/dx; 82 } 83 ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 84 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 85 } 86 if (ctype == IS_COLORING_GLOBAL) { 87 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 88 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 89 } 90 91 /* Loop over each color */ 92 nz = 0; 93 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 94 for (k=0; k<coloring->ncolors; k++) { 95 coloring->currentcolor = k; 96 97 /* 98 Loop over each column associated with color 99 adding the perturbation to the vector w3 = x1 + dx. 100 */ 101 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 102 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr); 103 if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array - start; 104 105 for (l=0; l<coloring->ncolumns[k]; l++) { 106 col = coloring->columns[k][l]; /* local column (in global index!) of the matrix we are probing for */ 107 /* convert global col to local colb */ 108 if (col >= start && col < end) { 109 colb = col-start; 110 } else { 111 #if defined(PETSC_USE_CTABLE) 112 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 113 colb--; 114 #else 115 colb = aij->colmap[col] - 1; /* local column index */ 116 #endif 117 colb += end - start; 118 } 119 w3_array[col] += 1.0/vscale_array[colb]; 120 } 121 if (ctype == IS_COLORING_GLOBAL) w3_array = w3_array + start; 122 ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 123 124 /* 125 Evaluate function at w3 = x1 + dx (here dx is a vector of perturbations) 126 w2 = F(x1 + dx) - F(x1) 127 */ 128 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 129 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 130 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 131 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 132 133 /* Loop over rows of vector, putting results into Jacobian matrix */ 134 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 135 for (l=0; l<coloring->nrows[k]; l++) { 136 row = Jentry[nz].row; /* local row index */ 137 *(Jentry[nz].valaddr) = y[row]*vscale_array[Jentry[nz].col]; 138 nz++; 139 } 140 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 141 } /* endof for each color */ 142 if (nz != cspA->nz + cspB->nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"nz %d != mat->nz %d\n",nz,cspA->nz+cspB->nz); 143 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 144 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 145 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 146 147 coloring->currentcolor = -1; 148 PetscFunctionReturn(0); 149 } 150 151 #undef __FUNCT__ 152 #define __FUNCT__ "MatFDColoringCreate_MPIAIJ_new" 153 PetscErrorCode MatFDColoringCreate_MPIAIJ_new(Mat mat,ISColoring iscoloring,MatFDColoring c) 154 { 155 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)mat->data; 156 PetscErrorCode ierr; 157 PetscMPIInt size,*ncolsonproc,*disp,nn; 158 PetscInt i,n,nrows,j,k,m,ncols,col; 159 const PetscInt *is,*A_ci,*A_cj,*B_ci,*B_cj,*row = NULL,*ltog=NULL; 160 PetscInt nis = iscoloring->n,nctot,*cols; 161 PetscInt *rowhit,cstart,cend,colb; 162 IS *isa; 163 ISLocalToGlobalMapping map = mat->cmap->mapping; 164 PetscInt ctype=c->ctype; 165 Mat A=aij->A,B=aij->B; 166 Mat_SeqAIJ *cspA=(Mat_SeqAIJ*)A->data, *cspB=(Mat_SeqAIJ*)B->data; 167 PetscScalar *A_val=cspA->a,*B_val=cspB->a; 168 PetscInt spidx; 169 PetscInt *spidxA,*spidxB,nz; 170 PetscScalar **valaddrhit; 171 MatEntry *Jentry; 172 173 PetscFunctionBegin; 174 if (ctype == IS_COLORING_GHOSTED) { 175 if(!map) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"When using ghosted differencing matrix must have local to global mapping provided with MatSetLocalToGlobalMapping"); 176 ierr = ISLocalToGlobalMappingGetIndices(map,<og);CHKERRQ(ierr); 177 } 178 179 m = mat->rmap->n; 180 cstart = mat->cmap->rstart; 181 cend = mat->cmap->rend; 182 c->M = mat->rmap->N; /* set the global rows and columns and local rows */ 183 c->N = mat->cmap->N; 184 c->m = m; 185 c->rstart = mat->rmap->rstart; 186 187 c->ncolors = nis; 188 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 189 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 190 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 191 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 192 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 193 ierr = PetscLogObjectMemory((PetscObject)c,5*nis*sizeof(PetscInt));CHKERRQ(ierr); 194 195 nz = cspA->nz + cspB->nz; /* total nonzero entries of mat */ 196 ierr = PetscMalloc(nz*sizeof(PetscScalar*),&c->valaddr);CHKERRQ(ierr); 197 ierr = PetscMalloc(nz*sizeof(MatEntry),&Jentry);CHKERRQ(ierr); 198 ierr = PetscLogObjectMemory((PetscObject)c,nz*sizeof(PetscScalar*));CHKERRQ(ierr); 199 200 /* Allow access to data structures of local part of matrix 201 - creates aij->colmap which maps global column number to local number in part B */ 202 if (!aij->colmap) { 203 ierr = MatCreateColmap_MPIAIJ_Private(mat);CHKERRQ(ierr); 204 } 205 ierr = MatGetColumnIJ_SeqAIJ_Color(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);CHKERRQ(ierr); 206 ierr = MatGetColumnIJ_SeqAIJ_Color(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);CHKERRQ(ierr); 207 208 ierr = PetscMalloc2(m+1,PetscInt,&rowhit,m+1,PetscScalar*,&valaddrhit);CHKERRQ(ierr); 209 nz = 0; 210 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 211 for (i=0; i<nis; i++) { /* for each local color */ 212 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 213 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 214 215 c->ncolumns[i] = n; /* local number of columns of this color on this process */ 216 if (n) { 217 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 218 ierr = PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));CHKERRQ(ierr); 219 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 220 } else { 221 c->columns[i] = 0; 222 } 223 224 if (ctype == IS_COLORING_GLOBAL) { 225 /* Determine nctot, the total (parallel) number of columns of this color */ 226 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr); 227 ierr = PetscMalloc2(size,PetscMPIInt,&ncolsonproc,size,PetscMPIInt,&disp);CHKERRQ(ierr); 228 229 /* ncolsonproc[j]: local ncolumns on proc[j] of this color */ 230 ierr = PetscMPIIntCast(n,&nn);CHKERRQ(ierr); 231 ierr = MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 232 nctot = 0; for (j=0; j<size; j++) nctot += ncolsonproc[j]; 233 if (!nctot) { 234 ierr = PetscInfo(mat,"Coloring of matrix has some unneeded colors with no corresponding rows\n");CHKERRQ(ierr); 235 } 236 237 disp[0] = 0; 238 for (j=1; j<size; j++) { 239 disp[j] = disp[j-1] + ncolsonproc[j-1]; 240 } 241 242 /* Get cols, the complete list of columns for this color on each process */ 243 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 244 ierr = MPI_Allgatherv((void*)is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 245 ierr = PetscFree2(ncolsonproc,disp);CHKERRQ(ierr); 246 } else if (ctype == IS_COLORING_GHOSTED) { 247 /* Determine local number of columns of this color on this process, including ghost points */ 248 nctot = n; 249 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 250 ierr = PetscMemcpy(cols,is,n*sizeof(PetscInt));CHKERRQ(ierr); 251 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this MatFDColoring type"); 252 253 /* Mark all rows affect by these columns */ 254 ierr = PetscMemzero(rowhit,m*sizeof(PetscInt));CHKERRQ(ierr); 255 for (j=0; j<nctot; j++) { /* loop over columns*/ 256 if (ctype == IS_COLORING_GHOSTED) { 257 col = ltog[cols[j]]; 258 } else { 259 col = cols[j]; 260 } 261 if (col >= cstart && col < cend) { /* column is in diagonal block of matrix A */ 262 row = A_cj + A_ci[col-cstart]; 263 nrows = A_ci[col-cstart+1] - A_ci[col-cstart]; 264 /* loop over columns of A marking them in rowhit */ 265 for (k=0; k<nrows; k++) { 266 /* set valaddrhit for part A */ 267 spidx = spidxA[A_ci[col-cstart] + k]; 268 valaddrhit[*row] = &A_val[spidx]; 269 rowhit[*row++] = col - cstart + 1; /* local column index */ 270 } 271 } else { /* column is in off-diagonal block of matrix B */ 272 #if defined(PETSC_USE_CTABLE) 273 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 274 colb--; 275 #else 276 colb = aij->colmap[col] - 1; /* local column index */ 277 #endif 278 if (colb == -1) { 279 nrows = 0; 280 } else { 281 row = B_cj + B_ci[colb]; 282 nrows = B_ci[colb+1] - B_ci[colb]; 283 } 284 /* loop over columns of B marking them in rowhit */ 285 for (k=0; k<nrows; k++) { 286 /* set valaddrhit for part B */ 287 spidx = spidxB[B_ci[colb] + k]; 288 valaddrhit[*row] = &B_val[spidx]; 289 rowhit[*row++] = colb + 1 + cend - cstart; /* local column index */ 290 } 291 } 292 } 293 294 /* count the number of hits */ 295 nrows = 0; 296 for (j=0; j<m; j++) { 297 if (rowhit[j]) nrows++; 298 } 299 c->nrows[i] = nrows; 300 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 301 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 302 ierr = PetscLogObjectMemory((PetscObject)c,2*(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr); 303 nrows = 0; 304 for (j=0; j<m; j++) { 305 if (rowhit[j]) { 306 Jentry[nz].row = j; /* local row index */ 307 Jentry[nz].col = rowhit[j] - 1; /* local column index */ 308 Jentry[nz].valaddr = valaddrhit[j]; /* address of mat value for this entry */ 309 310 c->rows[i][nrows] = j; /* local row index */ 311 c->columnsforrow[i][nrows] = rowhit[j] - 1; /* local column index */ 312 c->valaddr[nz] = valaddrhit[j]; /* address of mat value for this entry */ 313 nz++; 314 nrows++; 315 } 316 } 317 ierr = PetscFree(cols);CHKERRQ(ierr); 318 } 319 if (nz != cspA->nz + cspB->nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"nz %d != mat->nz %d\n",nz,cspA->nz+cspB->nz); 320 321 /* create vscale for storing dx */ 322 if (ctype == IS_COLORING_GLOBAL) { 323 ierr = VecCreateGhost(PetscObjectComm((PetscObject)mat),m,PETSC_DETERMINE,B->cmap->n,aij->garray,&c->vscale);CHKERRQ(ierr); 324 } 325 326 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 327 328 c->matentry = Jentry; 329 ierr = PetscFree2(rowhit,valaddrhit);CHKERRQ(ierr); 330 ierr = MatRestoreColumnIJ_SeqAIJ_Color(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);CHKERRQ(ierr); 331 ierr = MatRestoreColumnIJ_SeqAIJ_Color(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);CHKERRQ(ierr); 332 if (ctype == IS_COLORING_GHOSTED) { 333 ierr = ISLocalToGlobalMappingRestoreIndices(map,<og);CHKERRQ(ierr); 334 } 335 336 mat->ops->fdcoloringapply = MatFDColoringApply_MPIAIJ; 337 PetscFunctionReturn(0); 338 } 339 340 /*------------------------------------------------------*/ 341 #undef __FUNCT__ 342 #define __FUNCT__ "MatFDColoringCreate_MPIAIJ" 343 PetscErrorCode MatFDColoringCreate_MPIAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 344 { 345 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)mat->data; 346 PetscErrorCode ierr; 347 PetscMPIInt size,*ncolsonproc,*disp,nn; 348 PetscInt i,n,nrows,j,k,m,ncols,col; 349 const PetscInt *is,*A_ci,*A_cj,*B_ci,*B_cj,*rows = 0,*ltog; 350 PetscInt nis = iscoloring->n,nctot,*cols; 351 PetscInt *rowhit,M,cstart,cend,colb; 352 PetscInt *columnsforrow,l; 353 IS *isa; 354 PetscBool done,flg; 355 ISLocalToGlobalMapping map = mat->cmap->mapping; 356 PetscInt ctype=c->ctype; 357 PetscBool new_impl=PETSC_FALSE; 358 359 PetscFunctionBegin; 360 ierr = PetscOptionsName("-new","using new impls","",&new_impl);CHKERRQ(ierr); 361 if (new_impl){ 362 ierr = MatFDColoringCreate_MPIAIJ_new(mat,iscoloring,c);CHKERRQ(ierr); 363 PetscFunctionReturn(0); 364 } 365 if (ctype == IS_COLORING_GHOSTED && !map) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"When using ghosted differencing matrix must have local to global mapping provided with MatSetLocalToGlobalMapping"); 366 367 if (map) {ierr = ISLocalToGlobalMappingGetIndices(map,<og);CHKERRQ(ierr);} 368 else ltog = NULL; 369 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 370 371 M = mat->rmap->n; 372 cstart = mat->cmap->rstart; 373 cend = mat->cmap->rend; 374 c->M = mat->rmap->N; /* set the global rows and columns and local rows */ 375 c->N = mat->cmap->N; 376 c->m = mat->rmap->n; 377 c->rstart = mat->rmap->rstart; 378 379 c->ncolors = nis; 380 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 381 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 382 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 383 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 384 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 385 ierr = PetscLogObjectMemory((PetscObject)c,5*nis*sizeof(PetscInt));CHKERRQ(ierr); 386 387 /* Allow access to data structures of local part of matrix */ 388 if (!aij->colmap) { 389 ierr = MatCreateColmap_MPIAIJ_Private(mat);CHKERRQ(ierr); 390 } 391 ierr = MatGetColumnIJ(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr); 392 ierr = MatGetColumnIJ(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr); 393 394 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 395 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); 396 397 for (i=0; i<nis; i++) { 398 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 399 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 400 401 c->ncolumns[i] = n; /* local number of columns of this color on this process */ 402 if (n) { 403 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 404 ierr = PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));CHKERRQ(ierr); 405 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 406 } else { 407 c->columns[i] = 0; 408 } 409 410 if (ctype == IS_COLORING_GLOBAL) { 411 /* Determine the total (parallel) number of columns of this color */ 412 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr); 413 ierr = PetscMalloc2(size,PetscMPIInt,&ncolsonproc,size,PetscMPIInt,&disp);CHKERRQ(ierr); 414 415 ierr = PetscMPIIntCast(n,&nn);CHKERRQ(ierr); 416 ierr = MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 417 nctot = 0; for (j=0; j<size; j++) nctot += ncolsonproc[j]; 418 if (!nctot) { 419 ierr = PetscInfo(mat,"Coloring of matrix has some unneeded colors with no corresponding rows\n");CHKERRQ(ierr); 420 } 421 422 disp[0] = 0; 423 for (j=1; j<size; j++) { 424 disp[j] = disp[j-1] + ncolsonproc[j-1]; 425 } 426 427 /* Get complete list of columns for color on each processor */ 428 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 429 ierr = MPI_Allgatherv((void*)is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 430 ierr = PetscFree2(ncolsonproc,disp);CHKERRQ(ierr); 431 } else if (ctype == IS_COLORING_GHOSTED) { 432 /* Determine local number of columns of this color on this process, including ghost points */ 433 nctot = n; 434 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 435 ierr = PetscMemcpy(cols,is,n*sizeof(PetscInt));CHKERRQ(ierr); 436 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this MatFDColoring type"); 437 438 /* 439 Mark all rows affect by these columns 440 */ 441 /* Temporary option to allow for debugging/testing */ 442 flg = PETSC_FALSE; 443 ierr = PetscOptionsGetBool(NULL,"-matfdcoloring_slow",&flg,NULL);CHKERRQ(ierr); 444 if (!flg) { /*-----------------------------------------------------------------------------*/ 445 /* crude, fast version */ 446 ierr = PetscMemzero(rowhit,M*sizeof(PetscInt));CHKERRQ(ierr); 447 /* loop over columns*/ 448 for (j=0; j<nctot; j++) { 449 if (ctype == IS_COLORING_GHOSTED) { 450 col = ltog[cols[j]]; 451 } else { 452 col = cols[j]; 453 } 454 if (col >= cstart && col < cend) { 455 /* column is in diagonal block of matrix */ 456 rows = A_cj + A_ci[col-cstart]; 457 m = A_ci[col-cstart+1] - A_ci[col-cstart]; 458 } else { 459 #if defined(PETSC_USE_CTABLE) 460 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 461 colb--; 462 #else 463 colb = aij->colmap[col] - 1; 464 #endif 465 if (colb == -1) { 466 m = 0; 467 } else { 468 rows = B_cj + B_ci[colb]; 469 m = B_ci[colb+1] - B_ci[colb]; 470 } 471 } 472 /* loop over columns marking them in rowhit */ 473 for (k=0; k<m; k++) { 474 rowhit[*rows++] = col + 1; 475 } 476 } 477 478 /* count the number of hits */ 479 nrows = 0; 480 for (j=0; j<M; j++) { 481 if (rowhit[j]) nrows++; 482 } 483 c->nrows[i] = nrows; 484 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 485 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 486 ierr = PetscLogObjectMemory((PetscObject)c,2*(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr); 487 nrows = 0; 488 for (j=0; j<M; j++) { 489 if (rowhit[j]) { 490 c->rows[i][nrows] = j; /* local row index */ 491 c->columnsforrow[i][nrows] = rowhit[j] - 1; /* global column index */ 492 nrows++; 493 } 494 } 495 } else { /*-------------------------------------------------------------------------------*/ 496 /* slow version, using rowhit as a linked list */ 497 PetscInt currentcol,fm,mfm; 498 rowhit[M] = M; 499 nrows = 0; 500 /* loop over columns*/ 501 for (j=0; j<nctot; j++) { 502 if (ctype == IS_COLORING_GHOSTED) { 503 col = ltog[cols[j]]; 504 } else { 505 col = cols[j]; 506 } 507 if (col >= cstart && col < cend) { 508 /* column is in diagonal block of matrix */ 509 rows = A_cj + A_ci[col-cstart]; 510 m = A_ci[col-cstart+1] - A_ci[col-cstart]; 511 } else { 512 #if defined(PETSC_USE_CTABLE) 513 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 514 colb--; 515 #else 516 colb = aij->colmap[col] - 1; 517 #endif 518 if (colb == -1) { 519 m = 0; 520 } else { 521 rows = B_cj + B_ci[colb]; 522 m = B_ci[colb+1] - B_ci[colb]; 523 } 524 } 525 526 /* loop over columns marking them in rowhit */ 527 fm = M; /* fm points to first entry in linked list */ 528 for (k=0; k<m; k++) { 529 currentcol = *rows++; 530 /* is it already in the list? */ 531 do { 532 mfm = fm; 533 fm = rowhit[fm]; 534 } while (fm < currentcol); 535 /* not in list so add it */ 536 if (fm != currentcol) { 537 nrows++; 538 columnsforrow[currentcol] = col; 539 /* next three lines insert new entry into linked list */ 540 rowhit[mfm] = currentcol; 541 rowhit[currentcol] = fm; 542 fm = currentcol; 543 /* fm points to present position in list since we know the columns are sorted */ 544 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid coloring of matrix detected"); 545 } 546 } 547 c->nrows[i] = nrows; 548 549 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 550 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 551 ierr = PetscLogObjectMemory((PetscObject)c,(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr); 552 /* now store the linked list of rows into c->rows[i] */ 553 nrows = 0; 554 fm = rowhit[M]; 555 do { 556 c->rows[i][nrows] = fm; 557 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 558 fm = rowhit[fm]; 559 } while (fm < M); 560 } /* ---------------------------------------------------------------------------------------*/ 561 ierr = PetscFree(cols);CHKERRQ(ierr); 562 } 563 564 /* Optimize by adding the vscale, and scaleforrow[][] fields */ 565 /* 566 vscale will contain the "diagonal" on processor scalings followed by the off processor 567 */ 568 if (ctype == IS_COLORING_GLOBAL) { 569 ierr = VecCreateGhost(PetscObjectComm((PetscObject)mat),aij->A->rmap->n,PETSC_DETERMINE,aij->B->cmap->n,aij->garray,&c->vscale);CHKERRQ(ierr); 570 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 571 for (k=0; k<c->ncolors; k++) { 572 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 573 for (l=0; l<c->nrows[k]; l++) { 574 col = c->columnsforrow[k][l]; 575 if (col >= cstart && col < cend) { 576 /* column is in diagonal block of matrix */ 577 colb = col - cstart; 578 } else { 579 /* column is in "off-processor" part */ 580 #if defined(PETSC_USE_CTABLE) 581 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 582 colb--; 583 #else 584 colb = aij->colmap[col] - 1; 585 #endif 586 colb += cend - cstart; 587 } 588 c->vscaleforrow[k][l] = colb; 589 } 590 } 591 } else if (ctype == IS_COLORING_GHOSTED) { 592 /* Get gtol mapping */ 593 PetscInt N = mat->cmap->N,nlocal,*gtol; 594 ierr = PetscMalloc((N+1)*sizeof(PetscInt),>ol);CHKERRQ(ierr); 595 for (i=0; i<N; i++) gtol[i] = -1; 596 ierr = ISLocalToGlobalMappingGetSize(map,&nlocal);CHKERRQ(ierr); 597 for (i=0; i<nlocal; i++) gtol[ltog[i]] = i; 598 599 c->vscale = 0; /* will be created in MatFDColoringApply() */ 600 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 601 for (k=0; k<c->ncolors; k++) { 602 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 603 for (l=0; l<c->nrows[k]; l++) { 604 col = c->columnsforrow[k][l]; /* global column index */ 605 c->vscaleforrow[k][l] = gtol[col]; /* local column index */ 606 } 607 } 608 ierr = PetscFree(gtol);CHKERRQ(ierr); 609 } 610 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 611 612 ierr = PetscFree(rowhit);CHKERRQ(ierr); 613 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 614 ierr = MatRestoreColumnIJ(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr); 615 ierr = MatRestoreColumnIJ(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr); 616 if (map) {ierr = ISLocalToGlobalMappingRestoreIndices(map,<og);CHKERRQ(ierr);} 617 PetscFunctionReturn(0); 618 } 619 620 621 622 623 624 625