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