1 2 #ifdef PETSC_RCS_HEADER 3 static char vcid[] = "$Id: fdmatrix.c,v 1.40 1998/12/21 00:59:47 bsmith Exp bsmith $"; 4 #endif 5 6 /* 7 This is where the abstract matrix operations are defined that are 8 used for finite difference computations of Jacobians using coloring. 9 */ 10 11 #include "petsc.h" 12 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 13 #include "src/vec/vecimpl.h" 14 15 #undef __FUNC__ 16 #define __FUNC__ "MatFDColoringView_Draw" 17 static int MatFDColoringView_Draw(MatFDColoring fd,Viewer viewer) 18 { 19 int ierr,i,j,pause; 20 PetscTruth isnull; 21 Draw draw; 22 double xr,yr,xl,yl,h,w,x,y,xc,yc,scale = 0.0; 23 DrawButton button; 24 25 PetscFunctionBegin; 26 ierr = ViewerDrawGetDraw(viewer,0,&draw); CHKERRQ(ierr); 27 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 28 ierr = DrawSynchronizedClear(draw); CHKERRQ(ierr); 29 30 xr = fd->N; yr = fd->M; h = yr/10.0; w = xr/10.0; 31 xr += w; yr += h; xl = -w; yl = -h; 32 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 33 34 /* loop over colors */ 35 for (i=0; i<fd->ncolors; i++ ) { 36 for ( j=0; j<fd->nrows[i]; j++ ) { 37 y = fd->M - fd->rows[i][j] - fd->rstart; 38 x = fd->columnsforrow[i][j]; 39 ierr = DrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1); CHKERRQ(ierr); 40 } 41 } 42 ierr = DrawSynchronizedFlush(draw); CHKERRQ(ierr); 43 ierr = DrawGetPause(draw,&pause); CHKERRQ(ierr); 44 if (pause >= 0) { PetscSleep(pause); PetscFunctionReturn(0);} 45 ierr = DrawCheckResizedWindow(draw); CHKERRQ(ierr); 46 ierr = DrawSynchronizedGetMouseButton(draw,&button,&xc,&yc,0,0); CHKERRQ(ierr); 47 while (button != BUTTON_RIGHT) { 48 ierr = DrawSynchronizedClear(draw); CHKERRQ(ierr); 49 if (button == BUTTON_LEFT) scale = .5; 50 else if (button == BUTTON_CENTER) scale = 2.; 51 xl = scale*(xl + w - xc) + xc - w*scale; 52 xr = scale*(xr - w - xc) + xc + w*scale; 53 yl = scale*(yl + h - yc) + yc - h*scale; 54 yr = scale*(yr - h - yc) + yc + h*scale; 55 w *= scale; h *= scale; 56 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 57 /* loop over colors */ 58 for (i=0; i<fd->ncolors; i++ ) { 59 for ( j=0; j<fd->nrows[i]; j++ ) { 60 y = fd->M - fd->rows[i][j] - fd->rstart; 61 x = fd->columnsforrow[i][j]; 62 ierr = DrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1); CHKERRQ(ierr); 63 } 64 } 65 ierr = DrawCheckResizedWindow(draw); CHKERRQ(ierr); 66 ierr = DrawSynchronizedGetMouseButton(draw,&button,&xc,&yc,0,0); CHKERRQ(ierr); 67 } 68 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNC__ 73 #define __FUNC__ "MatFDColoringView" 74 /*@C 75 MatFDColoringView - Views a finite difference coloring context. 76 77 Collective on MatFDColoring unless Viewer is VIEWER_STDOUT_SELF 78 79 Input Parameters: 80 + c - the coloring context 81 - viewer - visualization context 82 83 Notes: 84 The available visualization contexts include 85 + VIEWER_STDOUT_SELF - standard output (default) 86 . VIEWER_STDOUT_WORLD - synchronized standard 87 output where only the first processor opens 88 the file. All other processors send their 89 data to the first processor to print. 90 - VIEWER_DRAW_WORLD - graphical display of nonzero structure 91 92 .seealso: MatFDColoringCreate() 93 94 .keywords: Mat, finite differences, coloring, view 95 @*/ 96 int MatFDColoringView(MatFDColoring c,Viewer viewer) 97 { 98 ViewerType vtype; 99 int i,j,format,ierr; 100 101 PetscFunctionBegin; 102 PetscValidHeaderSpecific(c,MAT_FDCOLORING_COOKIE); 103 if (viewer) {PetscValidHeader(viewer);} 104 else {viewer = VIEWER_STDOUT_SELF;} 105 106 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 107 if (PetscTypeCompare(vtype,DRAW_VIEWER)) { 108 ierr = MatFDColoringView_Draw(c,viewer); CHKERRQ(ierr); 109 PetscFunctionReturn(0); 110 } else if (PetscTypeCompare(vtype,ASCII_VIEWER)) { 111 ViewerASCIIPrintf(viewer,"MatFDColoring Object:\n"); 112 ViewerASCIIPrintf(viewer," Error tolerance=%g\n",c->error_rel); 113 ViewerASCIIPrintf(viewer," Umin=%g\n",c->umin); 114 ViewerASCIIPrintf(viewer," Number of colors=%d\n",c->ncolors); 115 116 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 117 if (format != VIEWER_FORMAT_ASCII_INFO) { 118 for ( i=0; i<c->ncolors; i++ ) { 119 ViewerASCIIPrintf(viewer," Information for color %d\n",i); 120 ViewerASCIIPrintf(viewer," Number of columns %d\n",c->ncolumns[i]); 121 for ( j=0; j<c->ncolumns[i]; j++ ) { 122 ViewerASCIIPrintf(viewer," %d\n",c->columns[i][j]); 123 } 124 ViewerASCIIPrintf(viewer," Number of rows %d\n",c->nrows[i]); 125 for ( j=0; j<c->nrows[i]; j++ ) { 126 ViewerASCIIPrintf(viewer," %d %d \n",c->rows[i][j],c->columnsforrow[i][j]); 127 } 128 } 129 } 130 } else { 131 SETERRQ(1,1,"Viewer type not supported for this object"); 132 } 133 PetscFunctionReturn(0); 134 } 135 136 #undef __FUNC__ 137 #define __FUNC__ "MatFDColoringSetParameters" 138 /*@ 139 MatFDColoringSetParameters - Sets the parameters for the sparse approximation of 140 a Jacobian matrix using finite differences. 141 142 Collective on MatFDColoring 143 144 The Jacobian is estimated with the differencing approximation 145 .vb 146 J(u)_{:,i} = [J(u+h*dx_{i}) - J(u)]/h where 147 h = error_rel*u[i] if abs(u[i]) > umin 148 = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i] 149 dx_{i} = (0, ... 1, .... 0) 150 .ve 151 152 Input Parameters: 153 + coloring - the coloring context 154 . error_rel - relative error 155 - umin - minimum allowable u-value magnitude 156 157 .keywords: Mat, finite differences, coloring, set, parameters 158 159 .seealso: MatFDColoringCreate() 160 @*/ 161 int MatFDColoringSetParameters(MatFDColoring matfd,double error,double umin) 162 { 163 PetscFunctionBegin; 164 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 165 166 if (error != PETSC_DEFAULT) matfd->error_rel = error; 167 if (umin != PETSC_DEFAULT) matfd->umin = umin; 168 PetscFunctionReturn(0); 169 } 170 171 #undef __FUNC__ 172 #define __FUNC__ "MatFDColoringSetFrequency" 173 /*@ 174 MatFDColoringSetFrequency - Sets the frequency for computing new Jacobian 175 matrices. 176 177 Collective on MatFDColoring 178 179 Input Parameters: 180 + coloring - the coloring context 181 - freq - frequency (default is 1) 182 183 Notes: 184 Using a modified Newton strategy, where the Jacobian remains fixed for several 185 iterations, can be cost effective in terms of overall nonlinear solution 186 efficiency. This parameter indicates that a new Jacobian will be computed every 187 <freq> nonlinear iterations. 188 189 Options Database Keys: 190 . -mat_fd_coloring_freq <freq> - Sets coloring frequency 191 192 .keywords: Mat, finite differences, coloring, set, frequency 193 194 .seealso: MatFDColoringCreate(), MatFDColoringGetFrequency() 195 @*/ 196 int MatFDColoringSetFrequency(MatFDColoring matfd,int freq) 197 { 198 PetscFunctionBegin; 199 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 200 201 matfd->freq = freq; 202 PetscFunctionReturn(0); 203 } 204 205 #undef __FUNC__ 206 #define __FUNC__ "MatFDColoringGetFrequency" 207 /*@ 208 MatFDColoringGetFrequency - Gets the frequency for computing new Jacobian 209 matrices. 210 211 Not Collective 212 213 Input Parameters: 214 . coloring - the coloring context 215 216 Output Parameters: 217 . freq - frequency (default is 1) 218 219 Notes: 220 Using a modified Newton strategy, where the Jacobian remains fixed for several 221 iterations, can be cost effective in terms of overall nonlinear solution 222 efficiency. This parameter indicates that a new Jacobian will be computed every 223 <freq> nonlinear iterations. 224 225 Options Database Keys: 226 . -mat_fd_coloring_freq <freq> - Sets coloring frequency 227 228 .keywords: Mat, finite differences, coloring, get, frequency 229 230 .seealso: MatFDColoringSetFrequency() 231 @*/ 232 int MatFDColoringGetFrequency(MatFDColoring matfd,int *freq) 233 { 234 PetscFunctionBegin; 235 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 236 237 *freq = matfd->freq; 238 PetscFunctionReturn(0); 239 } 240 241 #undef __FUNC__ 242 #define __FUNC__ "MatFDColoringSetFunction" 243 /*@C 244 MatFDColoringSetFunction - Sets the function to use for computing the Jacobian. 245 246 Collective on MatFDColoring 247 248 Input Parameters: 249 + coloring - the coloring context 250 . f - the function 251 - fctx - the optional user-defined function context 252 253 .keywords: Mat, Jacobian, finite differences, set, function 254 @*/ 255 int MatFDColoringSetFunction(MatFDColoring matfd,int (*f)(void),void *fctx) 256 { 257 PetscFunctionBegin; 258 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 259 260 matfd->f = f; 261 matfd->fctx = fctx; 262 263 PetscFunctionReturn(0); 264 } 265 266 #undef __FUNC__ 267 #define __FUNC__ "MatFDColoringSetFromOptions" 268 /*@ 269 MatFDColoringSetFromOptions - Sets coloring finite difference parameters from 270 the options database. 271 272 Collective on MatFDColoring 273 274 The Jacobian is estimated with the differencing approximation 275 .vb 276 J(u)_{:,i} = [J(u+h*dx_{i}) - J(u)]/h where 277 h = error_rel*u[i] if abs(u[i]) > umin 278 = +/- error_rel*umin otherwise, with +/- determined by the sign of u[i] 279 dx_{i} = (0, ... 1, .... 0) 280 .ve 281 282 Input Parameter: 283 . coloring - the coloring context 284 285 Options Database Keys: 286 + -mat_fd_coloring_error <err> - Sets <err> (square root 287 of relative error in the function) 288 . -mat_fd_coloring_umin <umin> - Sets umin, the minimum allowable u-value magnitude 289 . -mat_fd_coloring_freq <freq> - Sets frequency of computing a new Jacobian 290 . -mat_fd_coloring_view - Activates basic viewing 291 . -mat_fd_coloring_view_info - Activates viewing info 292 - -mat_fd_coloring_view_draw - Activates drawing 293 294 .keywords: Mat, finite differences, parameters 295 @*/ 296 int MatFDColoringSetFromOptions(MatFDColoring matfd) 297 { 298 int ierr,flag,freq = 1; 299 double error = PETSC_DEFAULT,umin = PETSC_DEFAULT; 300 301 PetscFunctionBegin; 302 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 303 304 ierr = OptionsGetDouble(matfd->prefix,"-mat_fd_coloring_err",&error,&flag);CHKERRQ(ierr); 305 ierr = OptionsGetDouble(matfd->prefix,"-mat_fd_coloring_umin",&umin,&flag);CHKERRQ(ierr); 306 ierr = MatFDColoringSetParameters(matfd,error,umin); CHKERRQ(ierr); 307 ierr = OptionsGetInt(matfd->prefix,"-mat_fd_coloring_freq",&freq,&flag);CHKERRQ(ierr); 308 ierr = MatFDColoringSetFrequency(matfd,freq);CHKERRQ(ierr); 309 ierr = OptionsHasName(PETSC_NULL,"-help",&flag); CHKERRQ(ierr); 310 if (flag) { 311 ierr = MatFDColoringPrintHelp(matfd); CHKERRQ(ierr); 312 } 313 PetscFunctionReturn(0); 314 } 315 316 #undef __FUNC__ 317 #define __FUNC__ "MatFDColoringPrintHelp" 318 /*@ 319 MatFDColoringPrintHelp - Prints help message for matrix finite difference calculations 320 using coloring. 321 322 Collective on MatFDColoring 323 324 Input Parameter: 325 . fdcoloring - the MatFDColoring context 326 327 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringSetFromOptions() 328 @*/ 329 int MatFDColoringPrintHelp(MatFDColoring fd) 330 { 331 PetscFunctionBegin; 332 PetscValidHeaderSpecific(fd,MAT_FDCOLORING_COOKIE); 333 334 (*PetscHelpPrintf)(fd->comm,"-mat_fd_coloring_err <err>: set sqrt rel tol in function, defaults to %g\n",fd->error_rel); 335 (*PetscHelpPrintf)(fd->comm,"-mat_fd_coloring_umin <umin>: see users manual, defaults to %d\n",fd->umin); 336 (*PetscHelpPrintf)(fd->comm,"-mat_fd_coloring_freq <freq>: frequency that Jacobian is recomputed, defaults to %d\n",fd->freq); 337 (*PetscHelpPrintf)(fd->comm,"-mat_fd_coloring_view\n"); 338 (*PetscHelpPrintf)(fd->comm,"-mat_fd_coloring_view_draw\n"); 339 (*PetscHelpPrintf)(fd->comm,"-mat_fd_coloring_view_info\n"); 340 PetscFunctionReturn(0); 341 } 342 343 int MatFDColoringView_Private(MatFDColoring fd) 344 { 345 int ierr,flg; 346 347 PetscFunctionBegin; 348 ierr = OptionsHasName(PETSC_NULL,"-mat_fd_coloring_view",&flg); CHKERRQ(ierr); 349 if (flg) { 350 ierr = MatFDColoringView(fd,VIEWER_STDOUT_(fd->comm)); CHKERRQ(ierr); 351 } 352 ierr = OptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_info",&flg); CHKERRQ(ierr); 353 if (flg) { 354 ierr = ViewerPushFormat(VIEWER_STDOUT_(fd->comm),VIEWER_FORMAT_ASCII_INFO,PETSC_NULL);CHKERRQ(ierr); 355 ierr = MatFDColoringView(fd,VIEWER_STDOUT_(fd->comm)); CHKERRQ(ierr); 356 ierr = ViewerPopFormat(VIEWER_STDOUT_(fd->comm));CHKERRQ(ierr); 357 } 358 ierr = OptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_draw",&flg); CHKERRQ(ierr); 359 if (flg) { 360 ierr = MatFDColoringView(fd,VIEWER_DRAW_(fd->comm)); CHKERRQ(ierr); 361 ierr = ViewerFlush(VIEWER_DRAW_(fd->comm)); CHKERRQ(ierr); 362 } 363 PetscFunctionReturn(0); 364 } 365 366 #undef __FUNC__ 367 #define __FUNC__ "MatFDColoringCreate" 368 /*@C 369 MatFDColoringCreate - Creates a matrix coloring context for finite difference 370 computation of Jacobians. 371 372 Collective on Mat 373 374 Input Parameters: 375 + mat - the matrix containing the nonzero structure of the Jacobian 376 - iscoloring - the coloring of the matrix 377 378 Output Parameter: 379 . color - the new coloring context 380 381 Options Database Keys: 382 + -mat_fd_coloring_view - Activates basic viewing or coloring 383 . -mat_fd_coloring_view_draw - Activates drawing of coloring 384 - -mat_fd_coloring_view_info - Activates viewing of coloring info 385 386 .seealso: MatFDColoringDestroy() 387 @*/ 388 int MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color) 389 { 390 MatFDColoring c; 391 MPI_Comm comm; 392 int ierr,M,N; 393 394 PetscFunctionBegin; 395 ierr = MatGetSize(mat,&M,&N); CHKERRQ(ierr); 396 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Only for square matrices"); 397 398 PetscObjectGetComm((PetscObject)mat,&comm); 399 PetscHeaderCreate(c,_p_MatFDColoring,int,MAT_FDCOLORING_COOKIE,0,"MatFDColoring",comm, 400 MatFDColoringDestroy,MatFDColoringView); 401 PLogObjectCreate(c); 402 403 if (mat->ops->fdcoloringcreate) { 404 ierr = (*mat->ops->fdcoloringcreate)(mat,iscoloring,c); CHKERRQ(ierr); 405 } else { 406 SETERRQ(PETSC_ERR_SUP,0,"Code not yet written for this matrix type"); 407 } 408 409 c->error_rel = 1.e-8; 410 c->umin = 1.e-6; 411 c->freq = 1; 412 413 ierr = MatFDColoringView_Private(c); CHKERRQ(ierr); 414 415 *color = c; 416 417 PetscFunctionReturn(0); 418 } 419 420 #undef __FUNC__ 421 #define __FUNC__ "MatFDColoringDestroy" 422 /*@C 423 MatFDColoringDestroy - Destroys a matrix coloring context that was created 424 via MatFDColoringCreate(). 425 426 Collective on MatFDColoring 427 428 Input Parameter: 429 . c - coloring context 430 431 .seealso: MatFDColoringCreate() 432 @*/ 433 int MatFDColoringDestroy(MatFDColoring c) 434 { 435 int i,ierr; 436 437 PetscFunctionBegin; 438 if (--c->refct > 0) PetscFunctionReturn(0); 439 440 441 for ( i=0; i<c->ncolors; i++ ) { 442 if (c->columns[i]) PetscFree(c->columns[i]); 443 if (c->rows[i]) PetscFree(c->rows[i]); 444 if (c->columnsforrow[i]) PetscFree(c->columnsforrow[i]); 445 } 446 PetscFree(c->ncolumns); 447 PetscFree(c->columns); 448 PetscFree(c->nrows); 449 PetscFree(c->rows); 450 PetscFree(c->columnsforrow); 451 PetscFree(c->scale); 452 if (c->w1) { 453 ierr = VecDestroy(c->w1); CHKERRQ(ierr); 454 ierr = VecDestroy(c->w2); CHKERRQ(ierr); 455 ierr = VecDestroy(c->w3); CHKERRQ(ierr); 456 } 457 PLogObjectDestroy(c); 458 PetscHeaderDestroy(c); 459 PetscFunctionReturn(0); 460 } 461 462 #include "snes.h" 463 464 #undef __FUNC__ 465 #define __FUNC__ "MatFDColoringApply" 466 /*@ 467 MatFDColoringApply - Given a matrix for which a MatFDColoring context 468 has been created, computes the Jacobian for a function via finite differences. 469 470 Collective on MatFDColoring 471 472 Input Parameters: 473 + mat - location to store Jacobian 474 . coloring - coloring context created with MatFDColoringCreate() 475 . x1 - location at which Jacobian is to be computed 476 - sctx - optional context required by function (actually a SNES context) 477 478 Options Database Keys: 479 . -mat_fd_coloring_freq <freq> - Sets coloring frequency 480 481 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView() 482 483 .keywords: coloring, Jacobian, finite differences 484 @*/ 485 int MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 486 { 487 int k,fg,ierr,N,start,end,l,row,col,srow; 488 Scalar dx, mone = -1.0,*y,*scale = coloring->scale,*xx,*wscale = coloring->wscale; 489 double epsilon = coloring->error_rel, umin = coloring->umin; 490 MPI_Comm comm = coloring->comm; 491 Vec w1,w2,w3; 492 int (*f)(void *,Vec,Vec,void *) = ( int (*)(void *,Vec,Vec,void *))coloring->f; 493 void *fctx = coloring->fctx; 494 495 PetscFunctionBegin; 496 PetscValidHeaderSpecific(J,MAT_COOKIE); 497 PetscValidHeaderSpecific(coloring,MAT_FDCOLORING_COOKIE); 498 PetscValidHeaderSpecific(x1,VEC_COOKIE); 499 500 501 if (!coloring->w1) { 502 ierr = VecDuplicate(x1,&coloring->w1); CHKERRQ(ierr); 503 PLogObjectParent(coloring,coloring->w1); 504 ierr = VecDuplicate(x1,&coloring->w2); CHKERRQ(ierr); 505 PLogObjectParent(coloring,coloring->w2); 506 ierr = VecDuplicate(x1,&coloring->w3); CHKERRQ(ierr); 507 PLogObjectParent(coloring,coloring->w3); 508 } 509 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 510 511 ierr = OptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&fg); CHKERRQ(ierr); 512 if (fg) { 513 PLogInfo(coloring,"MatFDColoringApply: Not calling MatZeroEntries()\n"); 514 } else { 515 ierr = MatZeroEntries(J); CHKERRQ(ierr); 516 } 517 518 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERRQ(ierr); 519 ierr = VecGetSize(x1,&N); CHKERRQ(ierr); 520 ierr = (*f)(sctx,x1,w1,fctx); CHKERRQ(ierr); 521 522 PetscMemzero(wscale,N*sizeof(Scalar)); 523 /* 524 Loop over each color 525 */ 526 527 ierr = VecGetArray(x1,&xx); CHKERRQ(ierr); 528 for (k=0; k<coloring->ncolors; k++) { 529 ierr = VecCopy(x1,w3); CHKERRQ(ierr); 530 /* 531 Loop over each column associated with color adding the 532 perturbation to the vector w3. 533 */ 534 for (l=0; l<coloring->ncolumns[k]; l++) { 535 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 536 dx = xx[col-start]; 537 if (dx == 0.0) dx = 1.0; 538 #if !defined(USE_PETSC_COMPLEX) 539 if (dx < umin && dx >= 0.0) dx = umin; 540 else if (dx < 0.0 && dx > -umin) dx = -umin; 541 #else 542 if (PetscAbsScalar(dx) < umin && PetscReal(dx) >= 0.0) dx = umin; 543 else if (PetscReal(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 544 #endif 545 dx *= epsilon; 546 wscale[col] = 1.0/dx; 547 ierr = VecSetValues(w3,1,&col,&dx,ADD_VALUES);CHKERRQ(ierr); 548 } 549 550 /* 551 Evaluate function at x1 + dx (here dx is a vector of perturbations) 552 */ 553 ierr = (*f)(sctx,w3,w2,fctx); CHKERRQ(ierr); 554 ierr = VecAXPY(&mone,w1,w2); CHKERRQ(ierr); 555 /* Communicate scale to all processors */ 556 #if !defined(USE_PETSC_COMPLEX) 557 ierr = MPI_Allreduce(wscale,scale,N,MPI_DOUBLE,MPI_SUM,comm);CHKERRQ(ierr); 558 #else 559 ierr = MPI_Allreduce(wscale,scale,2*N,MPI_DOUBLE,MPI_SUM,comm);CHKERRQ(ierr); 560 #endif 561 /* 562 Loop over rows of vector, putting results into Jacobian matrix 563 */ 564 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 565 for (l=0; l<coloring->nrows[k]; l++) { 566 row = coloring->rows[k][l]; 567 col = coloring->columnsforrow[k][l]; 568 y[row] *= scale[col]; 569 srow = row + start; 570 ierr = MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 571 } 572 ierr = VecRestoreArray(w2,&y); CHKERRQ(ierr); 573 } 574 ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 575 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 576 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 577 PetscFunctionReturn(0); 578 } 579 580 #include "ts.h" 581 582 #undef __FUNC__ 583 #define __FUNC__ "MatFDColoringApplyTS" 584 /*@ 585 MatFDColoringApplyTS - Given a matrix for which a MatFDColoring context 586 has been created, computes the Jacobian for a function via finite differences. 587 588 Collective on Mat, MatFDColoring, and Vec 589 590 Input Parameters: 591 + mat - location to store Jacobian 592 . coloring - coloring context created with MatFDColoringCreate() 593 . x1 - location at which Jacobian is to be computed 594 - sctx - optional context required by function (actually a SNES context) 595 596 Options Database Keys: 597 . -mat_fd_coloring_freq <freq> - Sets coloring frequency 598 599 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView() 600 601 .keywords: coloring, Jacobian, finite differences 602 @*/ 603 int MatFDColoringApplyTS(Mat J,MatFDColoring coloring,double t,Vec x1,MatStructure *flag,void *sctx) 604 { 605 int k,fg,ierr,N,start,end,l,row,col,srow; 606 Scalar dx, mone = -1.0,*y,*scale = coloring->scale,*xx,*wscale = coloring->wscale; 607 double epsilon = coloring->error_rel, umin = coloring->umin; 608 MPI_Comm comm = coloring->comm; 609 Vec w1,w2,w3; 610 int (*f)(void *,double,Vec,Vec,void *) = ( int (*)(void *,double,Vec,Vec,void *))coloring->f; 611 void *fctx = coloring->fctx; 612 613 PetscFunctionBegin; 614 PetscValidHeaderSpecific(J,MAT_COOKIE); 615 PetscValidHeaderSpecific(coloring,MAT_FDCOLORING_COOKIE); 616 PetscValidHeaderSpecific(x1,VEC_COOKIE); 617 618 if (!coloring->w1) { 619 ierr = VecDuplicate(x1,&coloring->w1); CHKERRQ(ierr); 620 PLogObjectParent(coloring,coloring->w1); 621 ierr = VecDuplicate(x1,&coloring->w2); CHKERRQ(ierr); 622 PLogObjectParent(coloring,coloring->w2); 623 ierr = VecDuplicate(x1,&coloring->w3); CHKERRQ(ierr); 624 PLogObjectParent(coloring,coloring->w3); 625 } 626 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 627 628 ierr = OptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&fg); CHKERRQ(ierr); 629 if (fg) { 630 PLogInfo(coloring,"MatFDColoringApplyTS: Not calling MatZeroEntries()\n"); 631 } else { 632 ierr = MatZeroEntries(J); CHKERRQ(ierr); 633 } 634 635 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERRQ(ierr); 636 ierr = VecGetSize(x1,&N); CHKERRQ(ierr); 637 ierr = (*f)(sctx,t,x1,w1,fctx); CHKERRQ(ierr); 638 639 PetscMemzero(wscale,N*sizeof(Scalar)); 640 /* 641 Loop over each color 642 */ 643 644 ierr = VecGetArray(x1,&xx); CHKERRQ(ierr); 645 for (k=0; k<coloring->ncolors; k++) { 646 ierr = VecCopy(x1,w3); CHKERRQ(ierr); 647 /* 648 Loop over each column associated with color adding the 649 perturbation to the vector w3. 650 */ 651 for (l=0; l<coloring->ncolumns[k]; l++) { 652 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 653 dx = xx[col-start]; 654 if (dx == 0.0) dx = 1.0; 655 #if !defined(USE_PETSC_COMPLEX) 656 if (dx < umin && dx >= 0.0) dx = umin; 657 else if (dx < 0.0 && dx > -umin) dx = -umin; 658 #else 659 if (PetscAbsScalar(dx) < umin && PetscReal(dx) >= 0.0) dx = umin; 660 else if (PetscReal(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 661 #endif 662 dx *= epsilon; 663 wscale[col] = 1.0/dx; 664 ierr = VecSetValues(w3,1,&col,&dx,ADD_VALUES); CHKERRQ(ierr); 665 } 666 /* 667 Evaluate function at x1 + dx (here dx is a vector of perturbations) 668 */ 669 ierr = (*f)(sctx,t,w3,w2,fctx); CHKERRQ(ierr); 670 ierr = VecAXPY(&mone,w1,w2); CHKERRQ(ierr); 671 /* Communicate scale to all processors */ 672 #if !defined(USE_PETSC_COMPLEX) 673 ierr = MPI_Allreduce(wscale,scale,N,MPI_DOUBLE,MPI_SUM,comm);CHKERRQ(ierr); 674 #else 675 ierr = MPI_Allreduce(wscale,scale,2*N,MPI_DOUBLE,MPI_SUM,comm);CHKERRQ(ierr); 676 #endif 677 /* 678 Loop over rows of vector, putting results into Jacobian matrix 679 */ 680 ierr = VecGetArray(w2,&y); CHKERRQ(ierr); 681 for (l=0; l<coloring->nrows[k]; l++) { 682 row = coloring->rows[k][l]; 683 col = coloring->columnsforrow[k][l]; 684 y[row] *= scale[col]; 685 srow = row + start; 686 ierr = MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 687 } 688 ierr = VecRestoreArray(w2,&y); CHKERRQ(ierr); 689 } 690 ierr = VecRestoreArray(x1,&xx); CHKERRQ(ierr); 691 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 692 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 693 PetscFunctionReturn(0); 694 } 695 696 697 698