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