Actual source code: rosenbrock1f.F90

      1: !  Program usage: mpiexec -n 1 rosenbrock1f [-help] [all TAO options]
      2: !
      3: !  Description:  This example demonstrates use of the TAO package to solve an
      4: !  unconstrained minimization problem on a single processor.  We minimize the
      5: !  extended Rosenbrock function:
      6: !       sum_{i=0}^{n/2-1} (alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2)
      7: !
      8: !  The C version of this code is rosenbrock1.c
      9: !

     11: !

     13: ! ----------------------------------------------------------------------
     14: !
     15: #include "petsc/finclude/petsctao.h"
     16:       use petsctao
     17:       implicit none

     19: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     20: !                   Variable declarations
     21: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     22: !
     23: !  See additional variable declarations in the file rosenbrock1f.h

     25:       PetscErrorCode ierr    ! used to check for functions returning nonzeros
     26:       type(tVec) x       ! solution vector
     27:       type(tMat) H       ! hessian matrix
     28:       type(tTao) ta     ! TAO_SOVER context
     29:       PetscBool flg
     30:       PetscInt i2, i1
     31:       PetscMPIInt size
     32:       PetscReal zero
     33:       PetscReal alpha
     34:       PetscInt n
     35:       common/params/alpha, n

     37: !  Note: Any user-defined Fortran routines (such as FormGradient)
     38: !  MUST be declared as external.

     40:       external FormFunctionGradient, FormHessian

     42:       zero = 0.0d0
     43:       i2 = 2
     44:       i1 = 1

     46: !  Initialize TAO and PETSc
     47:       PetscCallA(PetscInitialize(ierr))

     49:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
     50:       PetscCheckA(size == 1, PETSC_COMM_SELF, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only')

     52: !  Initialize problem parameters
     53:       n = 2
     54:       alpha = 99.0d0

     56: ! Check for command line arguments to override defaults
     57:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
     58:       PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-alpha', alpha, flg, ierr))

     60: !  Allocate vectors for the solution and gradient
     61:       PetscCallA(VecCreateSeq(PETSC_COMM_SELF, n, x, ierr))

     63: !  Allocate storage space for Hessian
     64:       PetscCallA(MatCreateSeqBAIJ(PETSC_COMM_SELF, i2, n, n, i1, PETSC_NULL_INTEGER_ARRAY, H, ierr))

     66:       PetscCallA(MatSetOption(H, MAT_SYMMETRIC, PETSC_TRUE, ierr))

     68: !  The TAO code begins here

     70: !  Create TAO solver
     71:       PetscCallA(TaoCreate(PETSC_COMM_SELF, ta, ierr))
     72:       PetscCallA(TaoSetType(ta, TAOLMVM, ierr))

     74: !  Set routines for function, gradient, and hessian evaluation
     75:       PetscCallA(TaoSetObjectiveAndGradient(ta, PETSC_NULL_VEC, FormFunctionGradient, 0, ierr))
     76:       PetscCallA(TaoSetHessian(ta, H, H, FormHessian, 0, ierr))

     78: !  Optional: Set initial guess
     79:       PetscCallA(VecSet(x, zero, ierr))
     80:       PetscCallA(TaoSetSolution(ta, x, ierr))

     82: !  Check for TAO command line options
     83:       PetscCallA(TaoSetFromOptions(ta, ierr))

     85: !  SOLVE THE APPLICATION
     86:       PetscCallA(TaoSolve(ta, ierr))

     88: !  TaoView() prints ierr about the TAO solver; the option
     89: !      -tao_view
     90: !  can alternatively be used to activate this at runtime.
     91: !      PetscCallA(TaoView(ta,PETSC_VIEWER_STDOUT_SELF,ierr))

     93: !  Free TAO data structures
     94:       PetscCallA(TaoDestroy(ta, ierr))

     96: !  Free PETSc data structures
     97:       PetscCallA(VecDestroy(x, ierr))
     98:       PetscCallA(MatDestroy(H, ierr))

    100:       PetscCallA(PetscFinalize(ierr))
    101:     end

    103: ! --------------------------------------------------------------------
    104: !  FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
    105: !
    106: !  Input Parameters:
    107: !  tao - the Tao context
    108: !  X   - input vector
    109: !  dummy - not used
    110: !
    111: !  Output Parameters:
    112: !  G - vector containing the newly evaluated gradient
    113: !  f - function value

    115:     subroutine FormFunctionGradient(ta, X, f, G, dummy, ierr)
    116:       use petsctao
    117:       implicit none

    119:       type(tTao) ta
    120:       type(tVec) X, G
    121:       PetscReal f
    122:       PetscErrorCode ierr
    123:       PetscInt dummy

    125:       PetscReal ff, t1, t2
    126:       PetscInt i, nn
    127:       PetscReal, pointer :: g_v(:), x_v(:)
    128:       PetscReal alpha
    129:       PetscInt n
    130:       common/params/alpha, n

    132:       ierr = 0
    133:       nn = n/2
    134:       ff = 0

    136: !     Get pointers to vector data
    137:       PetscCall(VecGetArrayRead(X, x_v, ierr))
    138:       PetscCall(VecGetArray(G, g_v, ierr))

    140: !     Compute G(X)
    141:       do i = 0, nn - 1
    142:         t1 = x_v(1 + 2*i + 1) - x_v(1 + 2*i)*x_v(1 + 2*i)
    143:         t2 = 1.0 - x_v(1 + 2*i)
    144:         ff = ff + alpha*t1*t1 + t2*t2
    145:         g_v(1 + 2*i) = -4*alpha*t1*x_v(1 + 2*i) - 2.0*t2
    146:         g_v(1 + 2*i + 1) = 2.0*alpha*t1
    147:       end do

    149: !     Restore vectors
    150:       PetscCall(VecRestoreArrayRead(X, x_v, ierr))
    151:       PetscCall(VecRestoreArray(G, g_v, ierr))

    153:       f = ff
    154:       PetscCall(PetscLogFlops(15.0d0*nn, ierr))

    156:     end

    158: !
    159: ! ---------------------------------------------------------------------
    160: !
    161: !  FormHessian - Evaluates Hessian matrix.
    162: !
    163: !  Input Parameters:
    164: !  tao     - the Tao context
    165: !  X       - input vector
    166: !  dummy   - optional user-defined context, as set by SNESSetHessian()
    167: !            (not used here)
    168: !
    169: !  Output Parameters:
    170: !  H      - Hessian matrix
    171: !  PrecH  - optionally different matrix used to compute the preconditioner (not used here)
    172: !  ierr   - error code
    173: !
    174: !  Note: Providing the Hessian may not be necessary.  Only some solvers
    175: !  require this matrix.

    177:     subroutine FormHessian(ta, X, H, PrecH, dummy, ierr)
    178:       use petsctao
    179:       implicit none

    181: !  Input/output variables:
    182:       type(tTao) ta
    183:       type(tVec) X
    184:       type(tMat) H, PrecH
    185:       PetscErrorCode ierr
    186:       PetscInt dummy

    188:       PetscReal v(0:1, 0:1)
    189:       PetscBool assembled

    191: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    192: ! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr)
    193: ! will return an array of doubles referenced by x_array offset by x_index.
    194: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    195: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    196:       PetscReal, pointer :: x_v(:)
    197:       PetscInt i, nn, ind(0:1), i2
    198:       PetscReal alpha
    199:       PetscInt n
    200:       common/params/alpha, n

    202:       ierr = 0
    203:       nn = n/2
    204:       i2 = 2

    206: !  Zero existing matrix entries
    207:       PetscCall(MatAssembled(H, assembled, ierr))
    208:       if (assembled .eqv. PETSC_TRUE) PetscCall(MatZeroEntries(H, ierr))

    210: !  Get a pointer to vector data

    212:       PetscCall(VecGetArrayRead(X, x_v, ierr))

    214: !  Compute Hessian entries

    216:       do i = 0, nn - 1
    217:         v(1, 1) = 2.0*alpha
    218:         v(0, 0) = -4.0*alpha*(x_v(1 + 2*i + 1) - 3*x_v(1 + 2*i)*x_v(1 + 2*i)) + 2
    219:         v(1, 0) = -4.0*alpha*x_v(1 + 2*i)
    220:         v(0, 1) = v(1, 0)
    221:         ind(0) = 2*i
    222:         ind(1) = 2*i + 1
    223:         PetscCall(MatSetValues(H, i2, ind, i2, ind, reshape(v, [i2*i2]), INSERT_VALUES, ierr))
    224:       end do

    226: !  Restore vector

    228:       PetscCall(VecRestoreArrayRead(X, x_v, ierr))

    230: !  Assemble matrix

    232:       PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY, ierr))
    233:       PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY, ierr))

    235:       PetscCall(PetscLogFlops(9.0d0*nn, ierr))

    237:     end

    239: !
    240: !/*TEST
    241: !
    242: !   build:
    243: !      requires: !complex
    244: !
    245: !   test:
    246: !      args: -tao_monitor_short -tao_type ntr -tao_gatol 1.e-5
    247: !      requires: !single
    248: !
    249: !TEST*/