1! 2! Description: This example solves a nonlinear system on 1 processor with SNES. 3! We solve the Bratu (SFI - solid fuel ignition) problem in a 2D rectangular 4! domain. The command line options include: 5! -par <parameter>, where <parameter> indicates the nonlinearity of the problem 6! problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81) 7! -mx <xg>, where <xg> = number of grid points in the x-direction 8! -my <yg>, where <yg> = number of grid points in the y-direction 9! 10 11! 12! -------------------------------------------------------------------------- 13! 14! Solid Fuel Ignition (SFI) problem. This problem is modeled by 15! the partial differential equation 16! 17! -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 18! 19! with boundary conditions 20! 21! u = 0 for x = 0, x = 1, y = 0, y = 1. 22! 23! A finite difference approximation with the usual 5-point stencil 24! is used to discretize the boundary value problem to obtain a nonlinear 25! system of equations. 26! 27! The parallel version of this code is snes/tutorials/ex5f.F 28! 29! -------------------------------------------------------------------------- 30 subroutine postcheck(snes, x, y, w, changed_y, changed_w, ctx, ierr) 31#include <petsc/finclude/petscsnes.h> 32 use petscsnes 33 implicit none 34 SNES snes 35 PetscReal norm 36 Vec tmp, x, y, w 37 PetscBool changed_w, changed_y 38 PetscErrorCode ierr 39 PetscInt ctx 40 PetscScalar mone 41 MPI_Comm comm 42 43 character(len=PETSC_MAX_PATH_LEN) :: outputString 44 45 PetscCallA(PetscObjectGetComm(snes, comm, ierr)) 46 PetscCallA(VecDuplicate(x, tmp, ierr)) 47 mone = -1.0 48 PetscCallA(VecWAXPY(tmp, mone, x, w, ierr)) 49 PetscCallA(VecNorm(tmp, NORM_2, norm, ierr)) 50 PetscCallA(VecDestroy(tmp, ierr)) 51 write (outputString, *) norm 52 PetscCallA(PetscPrintf(comm, 'Norm of search step '//trim(outputString)//'\n', ierr)) 53 end 54 55 program main 56#include <petsc/finclude/petscdraw.h> 57 use petscdraw 58 use petscsnes 59 implicit none 60! 61! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62! Variable declarations 63! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 64! 65! Variables: 66! snes - nonlinear solver 67! x, r - solution, residual vectors 68! J - Jacobian matrix 69! its - iterations for convergence 70! matrix_free - flag - 1 indicates matrix-free version 71! lambda - nonlinearity parameter 72! draw - drawing context 73! 74 SNES snes 75 MatColoring mc 76 Vec x, r 77 PetscDraw draw 78 Mat J 79 PetscBool matrix_free, flg, fd_coloring 80 PetscErrorCode ierr 81 PetscInt its, N, mx, my, i5 82 PetscMPIInt size, rank 83 PetscReal lambda_max, lambda_min, lambda 84 MatFDColoring fdcoloring 85 ISColoring iscoloring 86 PetscBool pc 87 integer4 imx, imy 88 external postcheck 89 character(len=PETSC_MAX_PATH_LEN) :: outputString 90 PetscScalar, pointer :: lx_v(:) 91 integer4 xl, yl, width, height 92 93! Store parameters in common block 94 95 common/params/lambda, mx, my, fd_coloring 96 97! Note: Any user-defined Fortran routines (such as FormJacobian) 98! MUST be declared as external. 99 100 external FormFunction, FormInitialGuess, FormJacobian 101 102! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 103! Initialize program 104! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 105 106 PetscCallA(PetscInitialize(ierr)) 107 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)) 108 PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr)) 109 110 PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only') 111 112! Initialize problem parameters 113 i5 = 5 114 lambda_max = 6.81 115 lambda_min = 0.0 116 lambda = 6.0 117 mx = 4 118 my = 4 119 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mx', mx, flg, ierr)) 120 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-my', my, flg, ierr)) 121 PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-par', lambda, flg, ierr)) 122 PetscCheckA(lambda < lambda_max .and. lambda > lambda_min, PETSC_COMM_SELF, PETSC_ERR_USER, 'Lambda out of range ') 123 N = mx*my 124 pc = PETSC_FALSE 125 PetscCallA(PetscOptionsGetBool(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-pc', pc, PETSC_NULL_BOOL, ierr)) 126 127! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 128! Create nonlinear solver context 129! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 130 131 PetscCallA(SNESCreate(PETSC_COMM_WORLD, snes, ierr)) 132 133 if (pc .eqv. PETSC_TRUE) then 134 PetscCallA(SNESSetType(snes, SNESNEWTONTR, ierr)) 135 PetscCallA(SNESNewtonTRSetPostCheck(snes, postcheck, snes, ierr)) 136 end if 137 138! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139! Create vector data structures; set function evaluation routine 140! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 141 142 PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr)) 143 PetscCallA(VecSetSizes(x, PETSC_DECIDE, N, ierr)) 144 PetscCallA(VecSetFromOptions(x, ierr)) 145 PetscCallA(VecDuplicate(x, r, ierr)) 146 147! Set function evaluation routine and vector. Whenever the nonlinear 148! solver needs to evaluate the nonlinear function, it will call this 149! routine. 150! - Note that the final routine argument is the user-defined 151! context that provides application-specific data for the 152! function evaluation routine. 153 154 PetscCallA(SNESSetFunction(snes, r, FormFunction, fdcoloring, ierr)) 155 156! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 157! Create matrix data structure; set Jacobian evaluation routine 158! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 159 160! Create matrix. Here we only approximately preallocate storage space 161! for the Jacobian. See the users manual for a discussion of better 162! techniques for preallocating matrix memory. 163 164 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-snes_mf', matrix_free, ierr)) 165 if (.not. matrix_free) then 166 PetscCallA(MatCreateSeqAIJ(PETSC_COMM_WORLD, N, N, i5, PETSC_NULL_INTEGER_ARRAY, J, ierr)) 167 end if 168 169! 170! This option will cause the Jacobian to be computed via finite differences 171! efficiently using a coloring of the columns of the matrix. 172! 173 fd_coloring = .false. 174 PetscCallA(PetscOptionsHasName(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-snes_fd_coloring', fd_coloring, ierr)) 175 if (fd_coloring) then 176 177! 178! This initializes the nonzero structure of the Jacobian. This is artificial 179! because clearly if we had a routine to compute the Jacobian we won't need 180! to use finite differences. 181! 182 PetscCallA(FormJacobian(snes, x, J, J, 0, ierr)) 183! 184! Color the matrix, i.e. determine groups of columns that share no common 185! rows. These columns in the Jacobian can all be computed simultaneously. 186! 187 PetscCallA(MatColoringCreate(J, mc, ierr)) 188 PetscCallA(MatColoringSetType(mc, MATCOLORINGNATURAL, ierr)) 189 PetscCallA(MatColoringSetFromOptions(mc, ierr)) 190 PetscCallA(MatColoringApply(mc, iscoloring, ierr)) 191 PetscCallA(MatColoringDestroy(mc, ierr)) 192! 193! Create the data structure that SNESComputeJacobianDefaultColor() uses 194! to compute the actual Jacobians via finite differences. 195! 196 PetscCallA(MatFDColoringCreate(J, iscoloring, fdcoloring, ierr)) 197 PetscCallA(MatFDColoringSetFunction(fdcoloring, FormFunction, fdcoloring, ierr)) 198 PetscCallA(MatFDColoringSetFromOptions(fdcoloring, ierr)) 199 PetscCallA(MatFDColoringSetUp(J, iscoloring, fdcoloring, ierr)) 200! 201! Tell SNES to use the routine SNESComputeJacobianDefaultColor() 202! to compute Jacobians. 203! 204 PetscCallA(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, fdcoloring, ierr)) 205 PetscCallA(SNESGetJacobian(snes, J, PETSC_NULL_MAT, PETSC_NULL_FUNCTION, PETSC_NULL_INTEGER, ierr)) 206 PetscCallA(SNESGetJacobian(snes, J, PETSC_NULL_MAT, PETSC_NULL_FUNCTION, fdcoloring, ierr)) 207 PetscCallA(ISColoringDestroy(iscoloring, ierr)) 208 209 else if (.not. matrix_free) then 210 211! Set Jacobian matrix data structure and default Jacobian evaluation 212! routine. Whenever the nonlinear solver needs to compute the 213! Jacobian matrix, it will call this routine. 214! - Note that the final routine argument is the user-defined 215! context that provides application-specific data for the 216! Jacobian evaluation routine. 217! - The user can override with: 218! -snes_fd : default finite differencing approximation of Jacobian 219! -snes_mf : matrix-free Newton-Krylov method with no preconditioning 220! (unless user explicitly sets preconditioner) 221! -snes_mf_operator : form matrix from which to construct the preconditioner as set by the user, 222! but use matrix-free approx for Jacobian-vector 223! products within Newton-Krylov method 224! 225 PetscCallA(SNESSetJacobian(snes, J, J, FormJacobian, 0, ierr)) 226 end if 227 228! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 229! Customize nonlinear solver; set runtime options 230! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 231 232! Set runtime options (e.g., -snes_monitor -snes_rtol <rtol> -ksp_type <type>) 233 234 PetscCallA(SNESSetFromOptions(snes, ierr)) 235 236! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 237! Evaluate initial guess; then solve nonlinear system. 238! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 239 240! Note: The user should initialize the vector, x, with the initial guess 241! for the nonlinear solver prior to calling SNESSolve(). In particular, 242! to employ an initial guess of zero, the user should explicitly set 243! this vector to zero by calling VecSet(). 244 245 PetscCallA(FormInitialGuess(x, ierr)) 246 PetscCallA(SNESSolve(snes, PETSC_NULL_VEC, x, ierr)) 247 PetscCallA(SNESGetIterationNumber(snes, its, ierr)) 248 write (outputString, *) its 249 PetscCallA(PetscPrintf(PETSC_COMM_WORLD, 'Number of SNES iterations = '//trim(outputString)//'\n', ierr)) 250 251! PetscDraw contour plot of solution 252 253 xl = 300 254 yl = 0 255 width = 300 256 height = 300 257 PetscCallA(PetscDrawCreate(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, 'Solution', xl, yl, width, height, draw, ierr)) 258 PetscCallA(PetscDrawSetFromOptions(draw, ierr)) 259 260 PetscCallA(VecGetArrayRead(x, lx_v, ierr)) 261 imx = int(mx, kind=kind(imx)) 262 imy = int(my, kind=kind(imy)) 263 PetscCallA(PetscDrawTensorContour(draw, imx, imy, PETSC_NULL_REAL_ARRAY, PETSC_NULL_REAL_ARRAY, lx_v, ierr)) 264 PetscCallA(VecRestoreArrayRead(x, lx_v, ierr)) 265 266! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 267! Free work space. All PETSc objects should be destroyed when they 268! are no longer needed. 269! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 270 271 if (.not. matrix_free) PetscCallA(MatDestroy(J, ierr)) 272 if (fd_coloring) PetscCallA(MatFDColoringDestroy(fdcoloring, ierr)) 273 274 PetscCallA(VecDestroy(x, ierr)) 275 PetscCallA(VecDestroy(r, ierr)) 276 PetscCallA(SNESDestroy(snes, ierr)) 277 PetscCallA(PetscDrawDestroy(draw, ierr)) 278 PetscCallA(PetscFinalize(ierr)) 279 end 280 281! --------------------------------------------------------------------- 282! 283! FormInitialGuess - Forms initial approximation. 284! 285! Input Parameter: 286! X - vector 287! 288! Output Parameters: 289! X - vector 290! ierr - error code 291! 292! Notes: 293! This routine serves as a wrapper for the lower-level routine 294! "ApplicationInitialGuess", where the actual computations are 295! done using the standard Fortran style of treating the local 296! vector data as a multidimensional array over the local mesh. 297! This routine merely accesses the local vector data via 298! VecGetArray() and VecRestoreArray(). 299! 300 subroutine FormInitialGuess(X, ierr) 301 use petscsnes 302 implicit none 303 304! Input/output variables: 305 Vec X 306 PetscErrorCode ierr 307 308! Declarations for use with local arrays: 309 PetscScalar, pointer :: lx_v(:) 310 311 ierr = 0 312 313! Get a pointer to vector data. 314! - VecGetArray() returns a pointer to the data array. 315! - You MUST call VecRestoreArray() when you no longer need access to 316! the array. 317 318 PetscCallA(VecGetArray(X, lx_v, ierr)) 319 320! Compute initial guess 321 322 PetscCallA(ApplicationInitialGuess(lx_v, ierr)) 323 324! Restore vector 325 326 PetscCallA(VecRestoreArray(X, lx_v, ierr)) 327 328 end 329 330! ApplicationInitialGuess - Computes initial approximation, called by 331! the higher level routine FormInitialGuess(). 332! 333! Input Parameter: 334! x - local vector data 335! 336! Output Parameters: 337! f - local vector data, f(x) 338! ierr - error code 339! 340! Notes: 341! This routine uses standard Fortran-style computations over a 2-dim array. 342! 343 subroutine ApplicationInitialGuess(x, ierr) 344 use petscksp 345 implicit none 346 347! Common blocks: 348 PetscReal lambda 349 PetscInt mx, my 350 PetscBool fd_coloring 351 common/params/lambda, mx, my, fd_coloring 352 353! Input/output variables: 354 PetscScalar x(mx, my) 355 PetscErrorCode ierr 356 357! Local variables: 358 PetscInt i, j 359 PetscReal temp1, temp, hx, hy, one 360 361! Set parameters 362 363 ierr = 0 364 one = 1.0 365 hx = one/(mx - 1) 366 hy = one/(my - 1) 367 temp1 = lambda/(lambda + one) 368 369 do 20 j = 1, my 370 temp = min(j - 1, my - j)*hy 371 do 10 i = 1, mx 372 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 373 x(i, j) = 0.0 374 else 375 x(i, j) = temp1*sqrt(min(min(i - 1, mx - i)*hx, temp)) 376 end if 37710 continue 37820 continue 379 380 end 381 382! --------------------------------------------------------------------- 383! 384! FormFunction - Evaluates nonlinear function, F(x). 385! 386! Input Parameters: 387! snes - the SNES context 388! X - input vector 389! dummy - optional user-defined context, as set by SNESSetFunction() 390! (not used here) 391! 392! Output Parameter: 393! F - vector with newly computed function 394! 395! Notes: 396! This routine serves as a wrapper for the lower-level routine 397! "ApplicationFunction", where the actual computations are 398! done using the standard Fortran style of treating the local 399! vector data as a multidimensional array over the local mesh. 400! This routine merely accesses the local vector data via 401! VecGetArray() and VecRestoreArray(). 402! 403 subroutine FormFunction(snes, X, F, fdcoloring, ierr) 404 use petscsnes 405 implicit none 406 407! Input/output variables: 408 SNES snes 409 Vec X, F 410 PetscErrorCode ierr 411 MatFDColoring fdcoloring 412 413! Common blocks: 414 PetscReal lambda 415 PetscInt mx, my 416 PetscBool fd_coloring 417 common/params/lambda, mx, my, fd_coloring 418 419! Declarations for use with local arrays: 420 PetscScalar, pointer :: lx_v(:), lf_v(:) 421 PetscInt, pointer :: indices(:) 422 423! Get pointers to vector data. 424! - VecGetArray() returns a pointer to the data array. 425! - You MUST call VecRestoreArray() when you no longer need access to 426! the array. 427 428 PetscCallA(VecGetArrayRead(X, lx_v, ierr)) 429 PetscCallA(VecGetArray(F, lf_v, ierr)) 430 431! Compute function 432 433 PetscCallA(ApplicationFunction(lx_v, lf_v, ierr)) 434 435! Restore vectors 436 437 PetscCallA(VecRestoreArrayRead(X, lx_v, ierr)) 438 PetscCallA(VecRestoreArray(F, lf_v, ierr)) 439 440 PetscCallA(PetscLogFlops(11.0d0*mx*my, ierr)) 441! 442! fdcoloring is in the common block and used here ONLY to test the 443! calls to MatFDColoringGetPerturbedColumns() and MatFDColoringRestorePerturbedColumns() 444! 445 if (fd_coloring) then 446 PetscCallA(MatFDColoringGetPerturbedColumns(fdcoloring, PETSC_NULL_INTEGER, indices, ierr)) 447 print *, 'Indices from GetPerturbedColumns' 448 write (*, 1000) indices 4491000 format(50i4) 450 PetscCallA(MatFDColoringRestorePerturbedColumns(fdcoloring, PETSC_NULL_INTEGER, indices, ierr)) 451 end if 452 end 453 454! --------------------------------------------------------------------- 455! 456! ApplicationFunction - Computes nonlinear function, called by 457! the higher level routine FormFunction(). 458! 459! Input Parameter: 460! x - local vector data 461! 462! Output Parameters: 463! f - local vector data, f(x) 464! ierr - error code 465! 466! Notes: 467! This routine uses standard Fortran-style computations over a 2-dim array. 468! 469 subroutine ApplicationFunction(x, f, ierr) 470 use petscsnes 471 implicit none 472 473! Common blocks: 474 PetscReal lambda 475 PetscInt mx, my 476 PetscBool fd_coloring 477 common/params/lambda, mx, my, fd_coloring 478 479! Input/output variables: 480 PetscScalar x(mx, my), f(mx, my) 481 PetscErrorCode ierr 482 483! Local variables: 484 PetscScalar two, one, hx, hy 485 PetscScalar hxdhy, hydhx, sc 486 PetscScalar u, uxx, uyy 487 PetscInt i, j 488 489 ierr = 0 490 one = 1.0 491 two = 2.0 492 hx = one/(mx - 1) 493 hy = one/(my - 1) 494 sc = hx*hy*lambda 495 hxdhy = hx/hy 496 hydhx = hy/hx 497 498! Compute function 499 500 do 20 j = 1, my 501 do 10 i = 1, mx 502 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 503 f(i, j) = x(i, j) 504 else 505 u = x(i, j) 506 uxx = hydhx*(two*u - x(i - 1, j) - x(i + 1, j)) 507 uyy = hxdhy*(two*u - x(i, j - 1) - x(i, j + 1)) 508 f(i, j) = uxx + uyy - sc*exp(u) 509 end if 51010 continue 51120 continue 512 513 end 514 515! --------------------------------------------------------------------- 516! 517! FormJacobian - Evaluates Jacobian matrix. 518! 519! Input Parameters: 520! snes - the SNES context 521! x - input vector 522! dummy - optional user-defined context, as set by SNESSetJacobian() 523! (not used here) 524! 525! Output Parameters: 526! jac - Jacobian matrix 527! jac_prec - optionally different matrix used to construct the preconditioner (not used here) 528! 529! Notes: 530! This routine serves as a wrapper for the lower-level routine 531! "ApplicationJacobian", where the actual computations are 532! done using the standard Fortran style of treating the local 533! vector data as a multidimensional array over the local mesh. 534! This routine merely accesses the local vector data via 535! VecGetArray() and VecRestoreArray(). 536! 537 subroutine FormJacobian(snes, X, jac, jac_prec, dummy, ierr) 538 use petscsnes 539 implicit none 540 541! Input/output variables: 542 SNES snes 543 Vec X 544 Mat jac, jac_prec 545 PetscErrorCode ierr 546 integer dummy 547 548! Common blocks: 549 PetscReal lambda 550 PetscInt mx, my 551 PetscBool fd_coloring 552 common/params/lambda, mx, my, fd_coloring 553 554! Declarations for use with local array: 555 PetscScalar, pointer :: lx_v(:) 556 557! Get a pointer to vector data 558 559 PetscCallA(VecGetArrayRead(X, lx_v, ierr)) 560 561! Compute Jacobian entries 562 563 PetscCallA(ApplicationJacobian(lx_v, jac, jac_prec, ierr)) 564 565! Restore vector 566 567 PetscCallA(VecRestoreArrayRead(X, lx_v, ierr)) 568 569! Assemble matrix 570 571 PetscCallA(MatAssemblyBegin(jac_prec, MAT_FINAL_ASSEMBLY, ierr)) 572 PetscCallA(MatAssemblyEnd(jac_prec, MAT_FINAL_ASSEMBLY, ierr)) 573 574 end 575 576! --------------------------------------------------------------------- 577! 578! ApplicationJacobian - Computes Jacobian matrix, called by 579! the higher level routine FormJacobian(). 580! 581! Input Parameters: 582! x - local vector data 583! 584! Output Parameters: 585! jac - Jacobian matrix 586! jac_prec - optionally different matrix used to construct the preconditioner (not used here) 587! ierr - error code 588! 589! Notes: 590! This routine uses standard Fortran-style computations over a 2-dim array. 591! 592 subroutine ApplicationJacobian(x, jac, jac_prec, ierr) 593 use petscsnes 594 implicit none 595 596! Common blocks: 597 PetscReal lambda 598 PetscInt mx, my 599 PetscBool fd_coloring 600 common/params/lambda, mx, my, fd_coloring 601 602! Input/output variables: 603 PetscScalar x(mx, my) 604 Mat jac, jac_prec 605 PetscErrorCode ierr 606 607! Local variables: 608 PetscInt i, j, row(1), col(5), i1, i5 609 PetscScalar two, one, hx, hy 610 PetscScalar hxdhy, hydhx, sc, v(5) 611 612! Set parameters 613 614 i1 = 1 615 i5 = 5 616 one = 1.0 617 two = 2.0 618 hx = one/(mx - 1) 619 hy = one/(my - 1) 620 sc = hx*hy 621 hxdhy = hx/hy 622 hydhx = hy/hx 623 624! Compute entries of the Jacobian matrix 625! - Here, we set all entries for a particular row at once. 626! - Note that MatSetValues() uses 0-based row and column numbers 627! in Fortran as well as in C. 628 629 do 20 j = 1, my 630 row(1) = (j - 1)*mx - 1 631 do 10 i = 1, mx 632 row(1) = row(1) + 1 633! boundary points 634 if (i == 1 .or. j == 1 .or. i == mx .or. j == my) then 635 PetscCallA(MatSetValues(jac_prec, i1, row, i1, row, [one], INSERT_VALUES, ierr)) 636! interior grid points 637 else 638 v(1) = -hxdhy 639 v(2) = -hydhx 640 v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(x(i, j)) 641 v(4) = -hydhx 642 v(5) = -hxdhy 643 col(1) = row(1) - mx 644 col(2) = row(1) - 1 645 col(3) = row(1) 646 col(4) = row(1) + 1 647 col(5) = row(1) + mx 648 PetscCallA(MatSetValues(jac_prec, i1, row, i5, col, v, INSERT_VALUES, ierr)) 649 end if 65010 continue 65120 continue 652 653 end 654 655! 656!/*TEST 657! 658! build: 659! requires: !single !complex 660! 661! test: 662! args: -snes_monitor_short -nox -snes_type newtontr -ksp_gmres_cgs_refinement_type refine_always 663! 664! test: 665! suffix: 2 666! args: -snes_monitor_short -nox -snes_fd -ksp_gmres_cgs_refinement_type refine_always 667! 668! test: 669! suffix: 3 670! args: -snes_monitor_short -nox -snes_fd_coloring -mat_coloring_type sl -ksp_gmres_cgs_refinement_type refine_always 671! filter: sort -b 672! filter_output: sort -b 673! 674! test: 675! suffix: 4 676! args: -pc -par 6.807 -nox 677! 678!TEST*/ 679