Actual source code: rosenbrock1f.F

      1: ! "$Id: rosenbrock1f.F 1.52 05/05/12 13:59:03-05:00 sarich@zorak.(none) $";
      2: !
      3: !  Program usage: mpirun -np 1 rosenbrock1f [-help] [all TAO options]
      4: !
      5: !  Description:  This example demonstrates use of the TAO package to solve an
      6: !  unconstrained minimization problem on a single processor.  We minimize the
      7: !  extended Rosenbrock function:
      8: !       sum_{i=0}^{n/2-1} ( alpha*(x_{2i+1}-x_{2i}^2)^2 + (1-x_{2i})^2 )
      9: !
     10: !  The C version of this code is rosenbrock1.c
     11: !
     12: !/*T
     13: !  Concepts: TAO - Solving an unconstrained minimization problem
     14: !  Routines: TaoInitialize(); TaoFinalize(); TaoSetOptions();
     15: !  Routines: TaoApplicationCreate(); TaoSetApplication();
     16: !  Routines: TaoCreate(); TaoAppSetObjectiveAndGradientRoutine();
     17: !  Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
     18: !  Routines: TaoAppSetInitialSolutionVec();
     19: !  Routines: TaoSolveApplication(); TaoDestroy(); TaoAppDestroy();
     20: !  Routines: TaoView(); TaoGetTerminationReason();
     21: !  Processors: 1
     22: !T*/
     23: !
     24: ! ----------------------------------------------------------------------
     25: !
     26:       implicit none

     28: #include "rosenbrock1f.h"

     30: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     31: !                   Variable declarations
     32: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     33: !
     34: !  See additional variable declarations in the file rosenbrock1f.h

     36:       integer          info    ! used to check for functions returning nonzeros
     37:       Vec              x       ! solution vector
     38:       Mat              H       ! hessian matrix
     39:       TAO_SOLVER       tao     ! TAO_SOVER context
     40:       TAO_APPLICATION  taoapp   ! TAO application context
     41:       PetscTruth       flg
     42:       integer          size,rank    ! number of processes running
     43:       PetscScalar      zero
     44:       TaoTerminateReason reason

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

     51:       external FormFunctionGradient,FormHessian

     54:       zero = 0.0d0

     56: !  Initialize TAO and PETSc
     57:       call PetscInitialize(PETSC_NULL_CHARACTER,info)
     58:       call TaoInitialize(PETSC_NULL_CHARACTER,info)

     60:       call MPI_Comm_size(PETSC_COMM_WORLD,size,info)
     61:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,info)
     62:       if (size .ne. 1) then
     63:          if (rank .eq. 0) then
     64:             write(6,*) 'This is a uniprocessor example only!'
     65:          endif
     66:          SETERRQ(1,' ',info)
     67:       endif

     69: !  Initialize problem parameters
     70:       n     = 2
     71:       alpha = 99.0d0

     75: ! Check for command line arguments to override defaults
     76:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-n',n,flg,           &
     77:      &                        info)
     78:       call PetscOptionsGetReal(PETSC_NULL_CHARACTER,'-alpha',            &
     79:      &                           alpha,flg,info)

     81: !  Allocate vectors for the solution and gradient
     82:       call VecCreateSeq(PETSC_COMM_SELF,n,x,info)

     84: !  Allocate storage space for Hessian;
     85:       call MatCreateSeqBDiag(PETSC_COMM_SELF,n,n,0,2,0,0,H,info)
     86:       call MatSetOption(H,MAT_SYMMETRIC,info)

     89: !  The TAO code begins here

     91: !  Create TAO solver
     92:       call TaoCreate(PETSC_COMM_SELF,TAO_NULL_CHARACTER,tao,info)
     93:       call TaoApplicationCreate(PETSC_COMM_SELF,taoapp,info)

     95: !  Set routines for function, gradient, and hessian evaluation

     97: !     TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
     98:       call TaoAppSetObjectiveAndGradientRo(taoapp,                      &
     99:      &     FormFunctionGradient,TAO_NULL_OBJECT,info)
    100:       call TaoAppSetHessianMat(taoapp,H,H,info)
    101:       call TaoAppSetHessianRoutine(taoapp,FormHessian,TAO_NULL_OBJECT,   &
    102:      &     info)

    104: !  Optional: Set initial guess
    105:       call VecSet(x, zero, info)
    106:       call TaoAppSetInitialSolutionVec(taoapp, x, info)

    108: !  Now that the PETSc application is all set, attach to tao context
    109:       call TaoSetApplication(tao,taoapp,info)

    111: !  Check for TAO command line options
    112:       call TaoSetOptions(taoapp,tao,info)

    114: !  SOLVE THE APPLICATION
    115:       call TaoSolveApplication(taoapp,tao,info)

    117:       call TaoGetTerminationReason(tao, reason, info)
    118:       if (reason .le. 0) then
    119:          print *,'Try a different TAO method, adjust some parameters,'
    120:          print *,'or check the function evaluation routines.'
    121:       endif

    123: !  TaoView() prints info about the TAO solver; the option
    124: !      -tao_view
    125: !  can alternatively be used to activate this at runtime.
    126: !     call TaoView(tao,info)
    127:

    129: !  Free TAO data structures
    130:       call TaoDestroy(tao,info)
    131:       call TaoAppDestroy(taoapp,info)

    133: !  Free PETSc data structures
    134:       call VecDestroy(x,info)
    135:       call MatDestroy(H,info)

    137: !  Finalize TAO
    138:       call TaoFinalize(info)
    139:       call PetscFinalize(info)

    141:       end

    144: ! --------------------------------------------------------------------
    145: !  FormFunctionGradient - Evaluates the function f(X) and gradient G(X)
    146: !
    147: !  Input Parameters:
    148: !  tao - the TAO_SOLVER context
    149: !  X   - input vector
    150: !  dummy - not used
    151: !
    152: !  Output Parameters:
    153: !  G - vector containing the newly evaluated gradient
    154: !  f - function value
    155:
    156:       subroutine FormFunctionGradient(taoapp, X, f, G, dummy, info)
    157:       implicit none

    159: ! n,alpha defined in rosenbrock1f.h
    160: #include "rosenbrock1f.h"

    162:       TAO_APPLICATION  taoapp
    163:       Vec              X,G
    164:       PetscScalar      f
    165:       integer          dummy, info

    168:       PetscScalar      ff,t1,t2
    169:       integer          i,nn

    171: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    172: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
    173: ! will return an array of doubles referenced by x_array offset by x_index.
    174: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    175: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    176:       PetscScalar      g_v(0:1),x_v(0:1)
    177:       PetscOffset      g_i,x_i

    179:       info = 0
    180:       nn = n/2
    181:       ff = 0

    183: !     Get pointers to vector data
    184:       call VecGetArray(X,x_v,x_i,info)
    185:       call VecGetArray(G,g_v,g_i,info)

    188: !     Compute G(X)
    189:       do i=0,nn-1
    190:          t1 = x_v(x_i+2*i+1) - x_v(x_i+2*i)*x_v(x_i+2*i)
    191:          t2 = 1.0 - x_v(x_i + 2*i)
    192:          ff = ff + alpha*t1*t1 + t2*t2
    193:          g_v(g_i + 2*i) = -4*alpha*t1*x_v(x_i + 2*i) - 2.0*t2
    194:          g_v(g_i + 2*i + 1) = 2.0*alpha*t1
    195:       enddo

    197: !     Restore vectors
    198:       call VecRestoreArray(X,x_v,x_i,info)
    199:       call VecRestoreArray(G,g_v,g_i,info)

    201:       f = ff
    202:       call PetscLogFlops(nn*15,info)

    204:       return
    205:       end

    207: !
    208: ! ---------------------------------------------------------------------
    209: !
    210: !  FormHessian - Evaluates Hessian matrix.
    211: !
    212: !  Input Parameters:
    213: !  tao     - the TAO_SOLVER context
    214: !  X       - input vector
    215: !  dummy   - optional user-defined context, as set by SNESSetHessian()
    216: !            (not used here)
    217: !
    218: !  Output Parameters:
    219: !  H      - Hessian matrix
    220: !  PrecH  - optionally different preconditioning matrix (not used here)
    221: !  flag   - flag indicating matrix structure
    222: !  info   - error code
    223: !
    224: !  Note: Providing the Hessian may not be necessary.  Only some solvers
    225: !  require this matrix.

    227:       subroutine FormHessian(taoapp,X,H,PrecH,flag,dummy,info)
    228:       implicit none

    230: #include "rosenbrock1f.h"

    232: !  Input/output variables:
    233:       TAO_APPLICATION  taoapp
    234:       Vec              X
    235:       Mat              H, PrecH
    236:       MatStructure     flag
    237:       integer          dummy,info
    238:
    239:       PetscScalar      v(0:1,0:1)
    240:       PetscTruth assembled

    242: ! PETSc's VecGetArray acts differently in Fortran than it does in C.
    243: ! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
    244: ! will return an array of doubles referenced by x_array offset by x_index.
    245: !  i.e.,  to reference the kth element of X, use x_array(k + x_index).
    246: ! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
    247:       PetscScalar      x_v(0:1)
    248:       PetscOffset      x_i
    249:       integer          i,nn,ind(0:1)

    252:       info = 0
    253:       nn= n/2

    255: !  Zero existing matrix entries
    256:       call MatAssembled(H,assembled,info)
    257:       if (assembled .eq. PETSC_TRUE) call MatZeroEntries(H,info)

    259: !  Get a pointer to vector data

    261:       call VecGetArray(X,x_v,x_i,info)

    263: !  Compute Hessian entries

    265:       do i=0,nn-1
    266:          v(1,1) = 2.0*alpha
    267:          v(0,0) = -4.0*alpha*(x_v(x_i+2*i+1) -                          &
    268:      &                3*x_v(x_i+2*i)*x_v(x_i+2*i))+2
    269:          v(1,0) = -4.0*alpha*x_v(x_i+2*i)
    270:          v(0,1) = v(1,0)
    271:          ind(0) = 2*i
    272:          ind(1) = 2*i + 1
    273:          call MatSetValues(H,2,ind,2,ind,v,INSERT_VALUES,info)
    274:       enddo

    276: !  Restore vector

    278:       call VecRestoreArray(X,x_v,x_i,info)

    280: !  Assemble matrix

    282:       call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
    283:       call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)

    285:       call PetscLogFlops(nn*9,info)

    287: !  Set flag to indicate that the Hessian matrix retains an identical
    288: !  nonzero structure throughout all nonlinear iterations (although the
    289: !  values of the entries change). Thus, we can save some work in setting
    290: !  up the preconditioner (e.g., no need to redo symbolic factorization for
    291: !  ICC preconditioners).
    292: !   - If the nonzero structure of the matrix is different during
    293: !     successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
    294: !     must be used instead.  If you are unsure whether the matrix
    295: !     structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
    296: !   - Caution:  If you specify SAME_NONZERO_PATTERN, the software
    297: !     believes your assertion and does not check the structure
    298: !     of the matrix.  If you erroneously claim that the structure
    299: !     is the same when it actually is not, the new preconditioner
    300: !     will not function correctly.  Thus, use this optimization
    301: !     feature with caution!

    303:       flag = SAME_NONZERO_PATTERN

    305:       return
    306:       end