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