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