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