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