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