1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: fdmatrix.c,v 1.14 1997/08/24 23:27:24 curfman Exp bsmith $"; 3 #endif 4 5 /* 6 This is where the abstract matrix operations are defined that are 7 used for finite difference computations of Jacobians using coloring. 8 */ 9 10 #include "petsc.h" 11 #include "src/mat/matimpl.h" /*I "mat.h" I*/ 12 #include "src/vec/vecimpl.h" 13 #include "pinclude/pviewer.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 ierr = ViewerDrawGetDraw(viewer,&draw); CHKERRQ(ierr); 26 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) return 0; 27 ierr = DrawSyncClear(draw); CHKERRQ(ierr); 28 29 xr = fd->N; yr = fd->M; h = yr/10.0; w = xr/10.0; 30 xr += w; yr += h; xl = -w; yl = -h; 31 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 32 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 DrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1); 40 } 41 } 42 ierr = DrawSyncFlush(draw); CHKERRQ(ierr); 43 ierr = DrawGetPause(draw,&pause); CHKERRQ(ierr); 44 if (pause >= 0) { PetscSleep(pause); return 0;} 45 ierr = DrawCheckResizedWindow(draw); 46 ierr = DrawSyncGetMouseButton(draw,&button,&xc,&yc,0,0); 47 while (button != BUTTON_RIGHT) { 48 ierr = DrawSyncClear(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 DrawRectangle(draw,x,y,x+1,y+1,i+1,i+1,i+1,i+1); 63 } 64 } 65 ierr = DrawCheckResizedWindow(draw); 66 ierr = DrawSyncGetMouseButton(draw,&button,&xc,&yc,0,0); 67 } 68 69 return 0; 70 } 71 72 73 #undef __FUNC__ 74 #define __FUNC__ "MatFDColoringView" 75 /*@C 76 MatFDColoringView - Views a finite difference coloring context. 77 78 Input Parameter: 79 . color - the coloring context 80 81 .seealso: MatFDColoringCreate() 82 83 @*/ 84 int MatFDColoringView(MatFDColoring color,Viewer viewer) 85 { 86 ViewerType vtype; 87 int i,j,format,ierr; 88 89 if (!viewer) viewer = VIEWER_STDOUT_WORLD; 90 91 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 92 if (vtype == DRAW_VIEWER) { 93 ierr = MatFDColoringView_Draw(color,viewer); CHKERRQ(ierr); 94 return 0; 95 } 96 97 printf("MatFDColoring Object:\n"); 98 printf(" Error tolerance %g\n",color->error_rel); 99 printf(" Umin %g\n",color->umin); 100 printf(" Number of colors %d\n",color->ncolors); 101 102 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 103 if (format != VIEWER_FORMAT_ASCII_INFO) { 104 for ( i=0; i<color->ncolors; i++ ) { 105 printf("Information for color %d\n",i); 106 printf("Number of columns %d\n",color->ncolumns[i]); 107 for ( j=0; j<color->ncolumns[i]; j++ ) { 108 printf(" %d\n",color->columns[i][j]); 109 } 110 printf("Number of rows %d\n",color->nrows[i]); 111 for ( j=0; j<color->nrows[i]; j++ ) { 112 printf(" %d %d \n",color->rows[i][j],color->columnsforrow[i][j]); 113 } 114 } 115 } 116 return 0; 117 } 118 119 #undef __FUNC__ 120 #define __FUNC__ "MatFDColoringSetParameters" 121 /*@ 122 MatFDColoringSetParameters - Sets the parameters for the approximation of 123 Jacobian using finite differences. 124 125 $ J(u)_{:,i} = [J(u+h*dx_{i}) - J(u)]/h where 126 $ h = error_rel*u[i] if u[i] > umin 127 $ = error_rel*umin else 128 $ 129 $ dx_{i} = (0, ... 1, .... 0) 130 131 Input Parameters: 132 . coloring - the jacobian coloring context 133 . error_rel - relative error 134 . umin - minimum allowable u-value 135 136 .keywords: SNES, Jacobian, finite differences, parameters 137 @*/ 138 int MatFDColoringSetParameters(MatFDColoring matfd,double error,double umin) 139 { 140 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 141 142 if (error != PETSC_DEFAULT) matfd->error_rel = error; 143 if (umin != PETSC_DEFAULT) matfd->umin = umin; 144 return 0; 145 } 146 147 #undef __FUNC__ 148 #define __FUNC__ "MatFDColoringSetFrequency" 149 /*@ 150 MatFDColoringSetFrequency - Sets the frequency at which new Jacobians 151 are computed. 152 153 Input Parameters: 154 . coloring - the jacobian coloring context 155 . freq - frequency 156 157 .keywords: SNES, Jacobian, finite differences, parameters 158 @*/ 159 int MatFDColoringSetFrequency(MatFDColoring matfd,int freq) 160 { 161 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 162 163 matfd->freq = freq; 164 return 0; 165 } 166 167 #undef __FUNC__ 168 #define __FUNC__ "MatFDColoringSetFunction" 169 /*@ 170 MatFDColoringSetFunction - Sets the function to use for computing the Jacobian. 171 172 Input Parameters: 173 . coloring - the jacobian coloring context 174 . f - the function 175 . fctx - the function context 176 177 178 .keywords: SNES, Jacobian, finite differences, parameters, function 179 @*/ 180 int MatFDColoringSetFunction(MatFDColoring matfd,int (*f)(void *,Vec,Vec,void *),void *fctx) 181 { 182 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 183 184 matfd->f = f; 185 matfd->fctx = fctx; 186 187 return 0; 188 } 189 190 #undef __FUNC__ 191 #define __FUNC__ "MatFDColoringSetFromOptions" 192 /*@ 193 MatFDColoringSetFromOptions - Set coloring finite difference parameters from 194 the options database. 195 196 $ J(u)_{:,i} = [J(u+h*dx_{i}) - J(u)]/h where 197 $ h = error_rel*u[i] if u[i] > umin 198 $ = error_rel*umin else 199 $ 200 $ dx_{i} = (0, ... 1, .... 0) 201 202 Input Parameters: 203 . coloring - the jacobian coloring context 204 205 Options Database: 206 . -mat_fd_coloring_error square root of relative error in function 207 . -mat_fd_coloring_umin see above 208 . -mat_fd_coloring_freq frequency at which Jacobian is computed 209 210 .keywords: SNES, Jacobian, finite differences, parameters 211 @*/ 212 int MatFDColoringSetFromOptions(MatFDColoring matfd) 213 { 214 int ierr,flag,freq = 1; 215 double error = PETSC_DEFAULT,umin = PETSC_DEFAULT; 216 PetscValidHeaderSpecific(matfd,MAT_FDCOLORING_COOKIE); 217 218 ierr = OptionsGetDouble(matfd->prefix,"-mat_fd_coloring_err",&error,&flag);CHKERRQ(ierr); 219 ierr = OptionsGetDouble(matfd->prefix,"-mat_fd_coloring_umin",&umin,&flag);CHKERRQ(ierr); 220 ierr = MatFDColoringSetParameters(matfd,error,umin); CHKERRQ(ierr); 221 ierr = OptionsGetInt(matfd->prefix,"-mat_fd_coloring_freq",&freq,&flag);CHKERRQ(ierr); 222 ierr = MatFDColoringSetFrequency(matfd,freq);CHKERRQ(ierr); 223 ierr = OptionsHasName(PETSC_NULL,"-help",&flag); CHKERRQ(ierr); 224 if (flag) { 225 ierr = MatFDColoringPrintHelp(matfd); CHKERRQ(ierr); 226 } 227 return 0; 228 } 229 230 #undef __FUNC__ 231 #define __FUNC__ "MatFDColoringPrintHelp" 232 /*@ 233 MatFDColoringPrintHelp - Prints help message for matrix finite difference calculations 234 using coloring. 235 236 Input Parameter: 237 . fdcoloring - the MatFDColoring context 238 239 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringSetFromOptions() 240 @*/ 241 int MatFDColoringPrintHelp(MatFDColoring fd) 242 { 243 PetscValidHeaderSpecific(fd,MAT_FDCOLORING_COOKIE); 244 245 PetscPrintf(fd->comm,"-mat_fd_coloring_err <err>: set sqrt rel tol in function, defaults to 1.e-8\n"); 246 PetscPrintf(fd->comm,"-mat_fd_coloring_umin <umin>: see users manual, defaults to 1.e-8\n"); 247 PetscPrintf(fd->comm,"-mat_fd_coloring_freq <freq>: frequency that Jacobian is recomputed, defaults 1\n"); 248 PetscPrintf(fd->comm,"-mat_fd_coloring_view\n"); 249 PetscPrintf(fd->comm,"-mat_fd_coloring_view_draw\n"); 250 PetscPrintf(fd->comm,"-mat_fd_coloring_view_info\n"); 251 return 0; 252 } 253 254 int MatFDColoringView_Private(MatFDColoring fd) 255 { 256 int ierr,flg; 257 258 ierr = OptionsHasName(PETSC_NULL,"-mat_fd_coloring_view",&flg); CHKERRQ(ierr); 259 if (flg) { 260 Viewer viewer; 261 ierr = ViewerFileOpenASCII(fd->comm,"stdout",&viewer);CHKERRQ(ierr); 262 ierr = MatFDColoringView(fd,viewer); CHKERRQ(ierr); 263 ierr = ViewerDestroy(viewer); CHKERRQ(ierr); 264 } 265 ierr = OptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_info",&flg); CHKERRQ(ierr); 266 if (flg) { 267 Viewer viewer; 268 ierr = ViewerFileOpenASCII(fd->comm,"stdout",&viewer);CHKERRQ(ierr); 269 ierr = ViewerPushFormat(viewer,VIEWER_FORMAT_ASCII_INFO,PETSC_NULL);CHKERRQ(ierr); 270 ierr = MatFDColoringView(fd,viewer); CHKERRQ(ierr); 271 ierr = ViewerPopFormat(viewer);CHKERRQ(ierr); 272 ierr = ViewerDestroy(viewer); CHKERRQ(ierr); 273 } 274 ierr = OptionsHasName(PETSC_NULL,"-mat_fd_coloring_view_draw",&flg); CHKERRQ(ierr); 275 if (flg) { 276 ierr = MatFDColoringView(fd,VIEWER_DRAWX_(fd->comm)); CHKERRQ(ierr); 277 ierr = ViewerFlush(VIEWER_DRAWX_(fd->comm)); CHKERRQ(ierr); 278 } 279 return 0; 280 } 281 282 #undef __FUNC__ 283 #define __FUNC__ "MatFDColoringCreate" 284 /*@C 285 MatFDColoringCreate - Creates a matrix coloring context for finite difference 286 computation of Jacobians. 287 288 Input Parameters: 289 . mat - the matrix containing the nonzero structure of the Jacobian 290 . iscoloring - the coloring of the matrix 291 292 Output Parameter: 293 . color - the new coloring context 294 295 Options Database: 296 . -mat_fd_coloring_error square root of relative error in function 297 . -mat_fd_coloring_umin see above 298 . -mat_fd_coloring_view 299 . -mat_fd_coloring_view_draw 300 301 .seealso: MatFDColoringDestroy() 302 @*/ 303 int MatFDColoringCreate(Mat mat,ISColoring iscoloring,MatFDColoring *color) 304 { 305 MatFDColoring c; 306 MPI_Comm comm; 307 int ierr,M,N; 308 309 ierr = MatGetSize(mat,&M,&N); CHKERRQ(ierr); 310 if (M != N) SETERRQ(PETSC_ERR_SUP,0,"Only for square matrices"); 311 312 PetscObjectGetComm((PetscObject)mat,&comm); 313 PetscHeaderCreate(c,_p_MatFDColoring,MAT_FDCOLORING_COOKIE,0,comm,MatFDColoringDestroy,MatFDColoringView); 314 PLogObjectCreate(c); 315 316 if (mat->ops.fdcoloringcreate) { 317 ierr = (*mat->ops.fdcoloringcreate)(mat,iscoloring,c); CHKERRQ(ierr); 318 } else { 319 SETERRQ(PETSC_ERR_SUP,0,"Code not yet written for this matrix type"); 320 } 321 322 c->error_rel = 1.e-8; 323 c->umin = 1.e-6; 324 c->freq = 1; 325 326 ierr = MatFDColoringView_Private(c); CHKERRQ(ierr); 327 328 *color = c; 329 330 return 0; 331 } 332 333 #undef __FUNC__ 334 #define __FUNC__ "MatFDColoringDestroy" 335 /*@C 336 MatFDColoringDestroy - Destroys a matrix coloring context that was created 337 via MatFDColoringCreate(). 338 339 Input Paramter: 340 . c - coloring context 341 342 .seealso: MatFDColoringCreate() 343 @*/ 344 int MatFDColoringDestroy(MatFDColoring c) 345 { 346 int i,ierr,flag; 347 348 if (--c->refct > 0) return 0; 349 350 ierr = OptionsHasName(PETSC_NULL,"-matfdcoloring_view",&flag); 351 if (flag) { 352 Viewer viewer; 353 ierr = ViewerFileOpenASCII(c->comm,"stdout",&viewer);CHKERRQ(ierr); 354 ierr = MatFDColoringView(c,viewer);CHKERRQ(ierr); 355 ierr = ViewerDestroy(viewer); CHKERRQ(ierr); 356 } 357 ierr = OptionsHasName(PETSC_NULL,"-matfdcoloring_view_info",&flag); 358 if (flag) { 359 Viewer viewer; 360 ierr = ViewerFileOpenASCII(c->comm,"stdout",&viewer);CHKERRQ(ierr); 361 ierr = ViewerPushFormat(viewer,VIEWER_FORMAT_ASCII_INFO, PETSC_NULL); CHKERRQ(ierr); 362 ierr = MatFDColoringView(c,viewer);CHKERRQ(ierr); 363 ierr = ViewerPopFormat(viewer); 364 ierr = ViewerDestroy(viewer); CHKERRQ(ierr); 365 } 366 367 for ( i=0; i<c->ncolors; i++ ) { 368 if (c->columns[i]) PetscFree(c->columns[i]); 369 if (c->rows[i]) PetscFree(c->rows[i]); 370 if (c->columnsforrow[i]) PetscFree(c->columnsforrow[i]); 371 } 372 PetscFree(c->ncolumns); 373 PetscFree(c->columns); 374 PetscFree(c->nrows); 375 PetscFree(c->rows); 376 PetscFree(c->columnsforrow); 377 PetscFree(c->scale); 378 if (c->w1) { 379 ierr = VecDestroy(c->w1); CHKERRQ(ierr); 380 ierr = VecDestroy(c->w2); CHKERRQ(ierr); 381 ierr = VecDestroy(c->w3); CHKERRQ(ierr); 382 } 383 PLogObjectDestroy(c); 384 PetscHeaderDestroy(c); 385 return 0; 386 } 387 388 #include "snes.h" 389 390 #undef __FUNC__ 391 #define __FUNC__ "MatFDColoringApply" 392 /*@ 393 MatFDColoringApply - Given a matrix for which a MatFDColoring has been created, 394 computes the Jacobian for a function via finite differences. 395 396 Input Parameters: 397 . mat - location to store Jacobian 398 . coloring - coloring context created with MatFDColoringCreate() 399 . x1 - location at which Jacobian is to be computed 400 . sctx - optional context required by function (actually a SNES context) 401 402 .seealso: MatFDColoringCreate(), MatFDColoringDestroy(), MatFDColoringView() 403 404 .keywords: coloring, Jacobian, finite differences 405 @*/ 406 int MatFDColoringApply(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 407 { 408 int k, ierr,N,start,end,l,row,col,srow; 409 Scalar dx, mone = -1.0,*y,*scale = coloring->scale,*xx,*wscale = coloring->wscale; 410 double epsilon = coloring->error_rel, umin = coloring->umin; 411 MPI_Comm comm = coloring->comm; 412 Vec w1,w2,w3; 413 int (*f)(void *,Vec,Vec,void *) = coloring->f; 414 void *fctx = coloring->fctx; 415 416 /* 417 ierr = SNESGetIterationNumber((SNES)sctx,&it); CHKERRQ(ierr); 418 if ((freq > 1) && ((it % freq) != 1)) { 419 PLogInfo(coloring,"MatFDColoringApply:Skipping Jacobian, iteration %d, freq %d\n",it,freq); 420 *flag = SAME_PRECONDITIONER; 421 return 0; 422 } else { 423 PLogInfo(coloring,"MatFDColoringApply:Computing Jacobian, iteration %d, freq %d\n",it,freq); 424 *flag = SAME_NONZERO_PATTERN; 425 }*/ 426 427 if (!coloring->w1) { 428 ierr = VecDuplicate(x1,&coloring->w1); CHKERRQ(ierr); 429 PLogObjectParent(coloring,coloring->w1); 430 ierr = VecDuplicate(x1,&coloring->w2); CHKERRQ(ierr); 431 PLogObjectParent(coloring,coloring->w2); 432 ierr = VecDuplicate(x1,&coloring->w3); CHKERRQ(ierr); 433 PLogObjectParent(coloring,coloring->w3); 434 } 435 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 436 437 ierr = MatZeroEntries(J); CHKERRQ(ierr); 438 439 ierr = VecGetOwnershipRange(x1,&start,&end); CHKERRQ(ierr); 440 ierr = VecGetSize(x1,&N); CHKERRQ(ierr); 441 ierr = VecGetArray(x1,&xx); CHKERRQ(ierr); 442 ierr = (*f)(sctx,x1,w1,fctx); CHKERRQ(ierr); 443 444 PetscMemzero(wscale,N*sizeof(Scalar)); 445 /* 446 Loop over each color 447 */ 448 449 for (k=0; k<coloring->ncolors; k++) { 450 ierr = VecCopy(x1,w3); CHKERRQ(ierr); 451 /* 452 Loop over each column associated with color adding the 453 perturbation to the vector w3. 454 */ 455 for (l=0; l<coloring->ncolumns[k]; l++) { 456 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 457 dx = xx[col-start]; 458 if (dx == 0.0) dx = 1.0; 459 #if !defined(PETSC_COMPLEX) 460 if (dx < umin && dx >= 0.0) dx = umin; 461 else if (dx < 0.0 && dx > -umin) dx = -umin; 462 #else 463 if (abs(dx) < umin && real(dx) >= 0.0) dx = umin; 464 else if (real(dx) < 0.0 && abs(dx) < umin) dx = -umin; 465 #endif 466 dx *= epsilon; 467 wscale[col] = 1.0/dx; 468 VecSetValues(w3,1,&col,&dx,ADD_VALUES); 469 } 470 VecRestoreArray(x1,&xx); 471 /* 472 Evaluate function at x1 + dx (here dx is a vector, of perturbations) 473 */ 474 ierr = (*f)(sctx,w3,w2,fctx); CHKERRQ(ierr); 475 ierr = VecAXPY(&mone,w1,w2); CHKERRQ(ierr); 476 /* Communicate scale to all processors */ 477 #if !defined(PETSC_COMPLEX) 478 MPI_Allreduce(wscale,scale,N,MPI_DOUBLE,MPI_SUM,comm); 479 #else 480 MPI_Allreduce(wscale,scale,2*N,MPI_DOUBLE,MPI_SUM,comm); 481 #endif 482 /* 483 Loop over rows of vector putting results into Jacobian matrix 484 */ 485 VecGetArray(w2,&y); 486 for (l=0; l<coloring->nrows[k]; l++) { 487 row = coloring->rows[k][l]; 488 col = coloring->columnsforrow[k][l]; 489 y[row] *= scale[col]; 490 srow = row + start; 491 ierr = MatSetValues(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 492 } 493 VecRestoreArray(w2,&y); 494 } 495 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 496 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 497 return 0; 498 } 499