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; 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 && !map) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_INCOMP,"When using ghosted differencing matrix must have local to global mapping provided with MatSetLocalToGlobalMapping"); 182 183 if (map) {ierr = ISLocalToGlobalMappingGetIndices(map,<og);CHKERRQ(ierr);} 184 else ltog = NULL; 185 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 186 187 m = mat->rmap->n; 188 cstart = mat->cmap->rstart; 189 cend = mat->cmap->rend; 190 c->M = mat->rmap->N; /* set the global rows and columns and local rows */ 191 c->N = mat->cmap->N; 192 c->m = m; 193 c->rstart = mat->rmap->rstart; 194 195 c->ncolors = nis; 196 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 197 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 198 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 199 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 200 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 201 ierr = PetscLogObjectMemory((PetscObject)c,5*nis*sizeof(PetscInt));CHKERRQ(ierr); 202 203 /* Allow access to data structures of local part of matrix 204 - creates aij->colmap which maps global column number to local number in part B */ 205 if (!aij->colmap) { 206 ierr = MatCreateColmap_MPIAIJ_Private(mat);CHKERRQ(ierr); 207 } 208 ierr = MatGetColumnIJ_SeqAIJ_Color(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&spidxA,NULL);CHKERRQ(ierr); 209 ierr = MatGetColumnIJ_SeqAIJ_Color(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&spidxB,NULL);CHKERRQ(ierr); 210 211 ierr = PetscMalloc2(m+1,PetscInt,&rowhit,m+1,PetscScalar*,&valaddrhit);CHKERRQ(ierr); 212 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); 213 nz = cspA->nz + cspB->nz; /* total nonzero entries of mat */ 214 ierr = PetscMalloc(nz*sizeof(PetscScalar*),&c->valaddr);CHKERRQ(ierr); 215 216 nz = 0; 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 (map) {ierr = ISLocalToGlobalMappingRestoreIndices(map,<og);CHKERRQ(ierr);} 373 374 mat->ops->fdcoloringapply = MatFDColoringApply_MPIAIJ; 375 PetscFunctionReturn(0); 376 } 377 378 /*------------------------------------------------------*/ 379 #undef __FUNCT__ 380 #define __FUNCT__ "MatFDColoringCreate_MPIAIJ" 381 PetscErrorCode MatFDColoringCreate_MPIAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 382 { 383 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)mat->data; 384 PetscErrorCode ierr; 385 PetscMPIInt size,*ncolsonproc,*disp,nn; 386 PetscInt i,n,nrows,j,k,m,ncols,col; 387 const PetscInt *is,*A_ci,*A_cj,*B_ci,*B_cj,*rows = 0,*ltog; 388 PetscInt nis = iscoloring->n,nctot,*cols; 389 PetscInt *rowhit,M,cstart,cend,colb; 390 PetscInt *columnsforrow,l; 391 IS *isa; 392 PetscBool done,flg; 393 ISLocalToGlobalMapping map = mat->cmap->mapping; 394 PetscInt ctype=c->ctype; 395 PetscBool new_impl=PETSC_FALSE; 396 397 PetscFunctionBegin; 398 ierr = PetscOptionsName("-new","using new impls","",&new_impl);CHKERRQ(ierr); 399 if (new_impl){ 400 ierr = MatFDColoringCreate_MPIAIJ_new(mat,iscoloring,c);CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 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"); 404 405 if (map) {ierr = ISLocalToGlobalMappingGetIndices(map,<og);CHKERRQ(ierr);} 406 else ltog = NULL; 407 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 408 409 M = mat->rmap->n; 410 cstart = mat->cmap->rstart; 411 cend = mat->cmap->rend; 412 c->M = mat->rmap->N; /* set the global rows and columns and local rows */ 413 c->N = mat->cmap->N; 414 c->m = mat->rmap->n; 415 c->rstart = mat->rmap->rstart; 416 417 c->ncolors = nis; 418 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 419 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 420 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 421 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 422 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 423 ierr = PetscLogObjectMemory((PetscObject)c,5*nis*sizeof(PetscInt));CHKERRQ(ierr); 424 425 /* Allow access to data structures of local part of matrix */ 426 if (!aij->colmap) { 427 ierr = MatCreateColmap_MPIAIJ_Private(mat);CHKERRQ(ierr); 428 } 429 ierr = MatGetColumnIJ(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr); 430 ierr = MatGetColumnIJ(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr); 431 432 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 433 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); 434 435 for (i=0; i<nis; i++) { 436 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 437 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 438 439 c->ncolumns[i] = n; /* local number of columns of this color on this process */ 440 if (n) { 441 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 442 ierr = PetscLogObjectMemory((PetscObject)c,n*sizeof(PetscInt));CHKERRQ(ierr); 443 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 444 } else { 445 c->columns[i] = 0; 446 } 447 448 if (ctype == IS_COLORING_GLOBAL) { 449 /* Determine the total (parallel) number of columns of this color */ 450 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);CHKERRQ(ierr); 451 ierr = PetscMalloc2(size,PetscMPIInt,&ncolsonproc,size,PetscMPIInt,&disp);CHKERRQ(ierr); 452 453 ierr = PetscMPIIntCast(n,&nn);CHKERRQ(ierr); 454 ierr = MPI_Allgather(&nn,1,MPI_INT,ncolsonproc,1,MPI_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 455 nctot = 0; for (j=0; j<size; j++) nctot += ncolsonproc[j]; 456 if (!nctot) { 457 ierr = PetscInfo(mat,"Coloring of matrix has some unneeded colors with no corresponding rows\n");CHKERRQ(ierr); 458 } 459 460 disp[0] = 0; 461 for (j=1; j<size; j++) { 462 disp[j] = disp[j-1] + ncolsonproc[j-1]; 463 } 464 465 /* Get complete list of columns for color on each processor */ 466 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 467 ierr = MPI_Allgatherv((void*)is,n,MPIU_INT,cols,ncolsonproc,disp,MPIU_INT,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr); 468 ierr = PetscFree2(ncolsonproc,disp);CHKERRQ(ierr); 469 } else if (ctype == IS_COLORING_GHOSTED) { 470 /* Determine local number of columns of this color on this process, including ghost points */ 471 nctot = n; 472 ierr = PetscMalloc((nctot+1)*sizeof(PetscInt),&cols);CHKERRQ(ierr); 473 ierr = PetscMemcpy(cols,is,n*sizeof(PetscInt));CHKERRQ(ierr); 474 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not provided for this MatFDColoring type"); 475 476 /* 477 Mark all rows affect by these columns 478 */ 479 /* Temporary option to allow for debugging/testing */ 480 flg = PETSC_FALSE; 481 ierr = PetscOptionsGetBool(NULL,"-matfdcoloring_slow",&flg,NULL);CHKERRQ(ierr); 482 if (!flg) { /*-----------------------------------------------------------------------------*/ 483 /* crude, fast version */ 484 ierr = PetscMemzero(rowhit,M*sizeof(PetscInt));CHKERRQ(ierr); 485 /* loop over columns*/ 486 for (j=0; j<nctot; j++) { 487 if (ctype == IS_COLORING_GHOSTED) { 488 col = ltog[cols[j]]; 489 } else { 490 col = cols[j]; 491 } 492 if (col >= cstart && col < cend) { 493 /* column is in diagonal block of matrix */ 494 rows = A_cj + A_ci[col-cstart]; 495 m = A_ci[col-cstart+1] - A_ci[col-cstart]; 496 } else { 497 #if defined(PETSC_USE_CTABLE) 498 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 499 colb--; 500 #else 501 colb = aij->colmap[col] - 1; 502 #endif 503 if (colb == -1) { 504 m = 0; 505 } else { 506 rows = B_cj + B_ci[colb]; 507 m = B_ci[colb+1] - B_ci[colb]; 508 } 509 } 510 /* loop over columns marking them in rowhit */ 511 for (k=0; k<m; k++) { 512 rowhit[*rows++] = col + 1; 513 } 514 } 515 516 /* count the number of hits */ 517 nrows = 0; 518 for (j=0; j<M; j++) { 519 if (rowhit[j]) nrows++; 520 } 521 c->nrows[i] = nrows; 522 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 523 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 524 ierr = PetscLogObjectMemory((PetscObject)c,2*(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr); 525 nrows = 0; 526 for (j=0; j<M; j++) { 527 if (rowhit[j]) { 528 c->rows[i][nrows] = j; /* local row index */ 529 c->columnsforrow[i][nrows] = rowhit[j] - 1; /* global column index */ 530 nrows++; 531 } 532 } 533 } else { /*-------------------------------------------------------------------------------*/ 534 /* slow version, using rowhit as a linked list */ 535 PetscInt currentcol,fm,mfm; 536 rowhit[M] = M; 537 nrows = 0; 538 /* loop over columns*/ 539 for (j=0; j<nctot; j++) { 540 if (ctype == IS_COLORING_GHOSTED) { 541 col = ltog[cols[j]]; 542 } else { 543 col = cols[j]; 544 } 545 if (col >= cstart && col < cend) { 546 /* column is in diagonal block of matrix */ 547 rows = A_cj + A_ci[col-cstart]; 548 m = A_ci[col-cstart+1] - A_ci[col-cstart]; 549 } else { 550 #if defined(PETSC_USE_CTABLE) 551 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 552 colb--; 553 #else 554 colb = aij->colmap[col] - 1; 555 #endif 556 if (colb == -1) { 557 m = 0; 558 } else { 559 rows = B_cj + B_ci[colb]; 560 m = B_ci[colb+1] - B_ci[colb]; 561 } 562 } 563 564 /* loop over columns marking them in rowhit */ 565 fm = M; /* fm points to first entry in linked list */ 566 for (k=0; k<m; k++) { 567 currentcol = *rows++; 568 /* is it already in the list? */ 569 do { 570 mfm = fm; 571 fm = rowhit[fm]; 572 } while (fm < currentcol); 573 /* not in list so add it */ 574 if (fm != currentcol) { 575 nrows++; 576 columnsforrow[currentcol] = col; 577 /* next three lines insert new entry into linked list */ 578 rowhit[mfm] = currentcol; 579 rowhit[currentcol] = fm; 580 fm = currentcol; 581 /* fm points to present position in list since we know the columns are sorted */ 582 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Invalid coloring of matrix detected"); 583 } 584 } 585 c->nrows[i] = nrows; 586 587 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 588 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 589 ierr = PetscLogObjectMemory((PetscObject)c,(nrows+1)*sizeof(PetscInt));CHKERRQ(ierr); 590 /* now store the linked list of rows into c->rows[i] */ 591 nrows = 0; 592 fm = rowhit[M]; 593 do { 594 c->rows[i][nrows] = fm; 595 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 596 fm = rowhit[fm]; 597 } while (fm < M); 598 } /* ---------------------------------------------------------------------------------------*/ 599 ierr = PetscFree(cols);CHKERRQ(ierr); 600 } 601 602 /* Optimize by adding the vscale, and scaleforrow[][] fields */ 603 /* 604 vscale will contain the "diagonal" on processor scalings followed by the off processor 605 */ 606 if (ctype == IS_COLORING_GLOBAL) { 607 ierr = VecCreateGhost(PetscObjectComm((PetscObject)mat),aij->A->rmap->n,PETSC_DETERMINE,aij->B->cmap->n,aij->garray,&c->vscale);CHKERRQ(ierr); 608 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 609 for (k=0; k<c->ncolors; k++) { 610 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 611 for (l=0; l<c->nrows[k]; l++) { 612 col = c->columnsforrow[k][l]; 613 if (col >= cstart && col < cend) { 614 /* column is in diagonal block of matrix */ 615 colb = col - cstart; 616 } else { 617 /* column is in "off-processor" part */ 618 #if defined(PETSC_USE_CTABLE) 619 ierr = PetscTableFind(aij->colmap,col+1,&colb);CHKERRQ(ierr); 620 colb--; 621 #else 622 colb = aij->colmap[col] - 1; 623 #endif 624 colb += cend - cstart; 625 } 626 c->vscaleforrow[k][l] = colb; 627 } 628 } 629 } else if (ctype == IS_COLORING_GHOSTED) { 630 /* Get gtol mapping */ 631 PetscInt N = mat->cmap->N,nlocal,*gtol; 632 ierr = PetscMalloc((N+1)*sizeof(PetscInt),>ol);CHKERRQ(ierr); 633 for (i=0; i<N; i++) gtol[i] = -1; 634 ierr = ISLocalToGlobalMappingGetSize(map,&nlocal);CHKERRQ(ierr); 635 for (i=0; i<nlocal; i++) gtol[ltog[i]] = i; 636 637 c->vscale = 0; /* will be created in MatFDColoringApply() */ 638 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 639 for (k=0; k<c->ncolors; k++) { 640 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 641 for (l=0; l<c->nrows[k]; l++) { 642 col = c->columnsforrow[k][l]; /* global column index */ 643 644 c->vscaleforrow[k][l] = gtol[col]; /* local column index */ 645 } 646 } 647 ierr = PetscFree(gtol);CHKERRQ(ierr); 648 } 649 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 650 651 ierr = PetscFree(rowhit);CHKERRQ(ierr); 652 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 653 ierr = MatRestoreColumnIJ(aij->A,0,PETSC_FALSE,PETSC_FALSE,&ncols,&A_ci,&A_cj,&done);CHKERRQ(ierr); 654 ierr = MatRestoreColumnIJ(aij->B,0,PETSC_FALSE,PETSC_FALSE,&ncols,&B_ci,&B_cj,&done);CHKERRQ(ierr); 655 if (map) {ierr = ISLocalToGlobalMappingRestoreIndices(map,<og);CHKERRQ(ierr);} 656 PetscFunctionReturn(0); 657 } 658 659 660 661 662 663 664