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