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