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