! "$Id: eptorsion2f.F 1.52 05/05/12 13:59:03-05:00 sarich@zorak.(none) $";
!
! Program usage: mpirun -np <proc> eptorsion2f [all TAO options]
!
! Description: This example demonstrates use of the TAO package to solve
! unconstrained minimization problems in parallel. This example is based
! on the Elastic-Plastic Torsion (dept) problem from the MINPACK-2 test suite.
! The command line options are:
! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction
! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction
! -par <param>, where <param> = angle of twist per unit length
!
!/*T
! Concepts: TAO - Solving an unconstrained minimization problem
! Routines: TaoInitialize(); TaoFinalize();
! Routines: TaoCreate(); TaoDestroy();
! Routines: TaoApplicationCreate(); TaoAppDestroy();
! Routines: TaoAppSetObjectiveAndGradientRoutine();
! Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
! Routines: TaoSetApplication(); TaoSetOptions();
! Routines: TaoAppSetInitialSolutionVec(); TaoSolveApplication(); TaoDestroy();
! Routines: TaoGetSolutionStatus();
! Processors: n
!T*/
!
! ----------------------------------------------------------------------
!
! Elastic-plastic torsion problem.
!
! The elastic plastic torsion problem arises from the determination
! of the stress field on an infinitely long cylindrical bar, which is
! equivalent to the solution of the following problem:
! min{ .5 * integral(||gradient(v(x))||^2 dx) - C * integral(v(x) dx)}
! where C is the torsion angle per unit length.
!
! The C version of this code is eptorsion2.c
!
! ----------------------------------------------------------------------
implicit none
#include "eptorsion2f.h"
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Variable declarations
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! See additional variable declarations in the file eptorsion2f.h
!
integer info ! used to check for functions returning nonzeros
Vec x ! solution vector
Mat H ! hessian matrix
integer Nx, Ny ! number of processes in x- and y- directions
TAO_SOLVER tao ! TAO_SOLVER solver context
TAO_APPLICATION torsionapp ! TAO application context (PETSc)
TaoTerminateReason reason
PetscTruth flg
integer iter ! iteration information
PetscScalar ff,gnorm,cnorm,xdiff
! Note: Any user-defined Fortran routines (such as FormGradient)
! MUST be declared as external.
external FormInitialGuess,FormFunctionGradient,ComputeHessian
! Initialize TAO, PETSc contexts
call PetscInitialize(PETSC_NULL_CHARACTER,info)
call TaoInitialize(PETSC_NULL_CHARACTER,info)
! Specify default parameters
param = 5.0d0
mx = 10
my = 10
Nx = PETSC_DECIDE
Ny = PETSC_DECIDE
! Check for any command line arguments that might override defaults
call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-mx",mx,flg,info)
call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-my",my,flg,info)
call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-par", &
& param,flg,info)
! Set up distributed array and vectors
call DACreate2d(MPI_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_BOX, &
& mx,my,Nx,Ny,1,1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &
& da,info)
! Create vectors
call DACreateGlobalVector(da,x,info)
call DACreateLocalVector(da,localX,info)
! Create Hessian
call DAGetMatrix(da,MATAIJ,H,info)
call MatSetOption(H,MAT_SYMMETRIC,info)
! The TAO code begins here
! Create TAO solver
call TaoCreate(MPI_COMM_WORLD,'tao_cg_fr',tao,info)
call TaoApplicationCreate(MPI_COMM_WORLD,torsionapp,info)
! Set routines for function and gradient evaluation
! TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
call TaoAppSetObjectiveAndGradientRo(torsionapp, &
& FormFunctionGradient,TAO_NULL_OBJECT,info)
call TaoAppSetHessianMat(torsionapp,H,H,info)
call TaoAppSetHessianRoutine(torsionapp,ComputeHessian, &
& TAO_NULL_OBJECT,info)
! Set initial guess
call FormInitialGuess(x,info)
call TaoAppSetInitialSolutionVec(torsionapp,x,info)
! Check for any TAO command line options
call TaoSetOptions(torsionapp, tao,info)
! Now that the PETSc application is all set, attach to tao context
call TaoSetApplication(tao,torsionapp,info)
! SOLVE THE APPLICATION
call TaoSolveApplication(torsionapp,tao,info)
! Get information on termination
call TaoGetSolutionStatus(tao,iter,ff,gnorm,cnorm,xdiff, &
& reason,info)
if (reason .lt. 0) then
print *,'TAO did not terminate successfully'
call TaoView(tao,info)
endif
! Free TAO data structures
call TaoDestroy(tao,info)
call TaoAppDestroy(torsionapp,info);
! Free PETSc data structures
call VecDestroy(x,info)
call VecDestroy(localX,info)
call MatDestroy(H,info)
call DADestroy(da,info)
! Finalize TAO and PETSc
call PetscFinalize(info)
call TaoFinalize(info)
end
! ---------------------------------------------------------------------
!
! FormInitialGuess - Computes an initial approximation to the solution.
!
! Input Parameters:
! X - vector
!
! Output Parameters:
! X - vector
! info - error code
!
subroutine FormInitialGuess(X,info)
implicit none
! mx, my are defined in eptorsion2f.h
#include "eptorsion2f.h"
! Input/output variables:
Vec X
integer info
! Local variables:
integer i, j, k, xe, ye
PetscScalar temp, val, hx, hy
integer xs, ys, xm, ym, gxm, gym, gxs, gys
hx = 1.0d0/(mx + 1)
hy = 1.0d0/(my + 1)
! Get corner information
call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
& PETSC_NULL_INTEGER,info)
call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, &
& gxm,gym,PETSC_NULL_INTEGER,info)
! Compute initial guess over locally owned part of mesh
xe = xs+xm
ye = ys+ym
do j=ys,ye-1
temp = min(j+1,my-j)*hy
do i=xs,xe-1
k = (j-gys)*gxm + i-gxs
val = min((min(i+1,mx-i))*hx,temp)
call VecSetValuesLocal(X,1,k,val,ADD_VALUES,info)
end do
end do
return
end
! ---------------------------------------------------------------------
!
! FormFunctionGradient - Evaluates gradient G(X).
!
! Input Parameters:
! tao - the TAO_SOLVER context
! X - input vector
! dummy - optional user-defined context (not used here)
!
! Output Parameters:
! f - the function value at X
! G - vector containing the newly evaluated gradient
! info - error code
!
! Notes:
! This routine serves as a wrapper for the lower-level routine
! "ApplicationGradient", where the actual computations are
! done using the standard Fortran style of treating the local
! input vector data as an array over the local mesh.
!
subroutine FormFunctionGradient(taoapp,X,f,G,dummy,info)
implicit none
! da, mx, my, param, localX declared in eptorsion2f.h
#include "eptorsion2f.h"
! Input/output variables:
TAO_APPLICATION taoapp
Vec X, G
PetscScalar f
integer dummy, info
! Declarations for use with local array:
! PETSc's VecGetArray acts differently in Fortran than it does in C.
! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (PetscOffset) x_index, info)
! will return an array of doubles referenced by x_array offset by x_index.
! i.e., to reference the kth element of X, use x_array(k + x_index).
! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice.
PetscScalar lx_v(0:1)
PetscOffset lx_i
! Local variables:
PetscScalar zero, p5, area, cdiv3, val, flin, fquad,floc
PetscScalar v, vb, vl, vr, vt, dvdx, dvdy, hx, hy
integer xe, ye, xsm, ysm, xep, yep, i, j, k, ind
integer xs, ys, xm, ym, gxs, gys, gxm, gym
info = 0
cdiv3 = param/3.0d0
zero = 0.0d0
p5 = 0.5d0
hx = 1.0d0/(mx + 1)
hy = 1.0d0/(my + 1)
fquad = zero
flin = zero
! Initialize gradient to zero
call VecSet(G,zero,info)
! Scatter ghost points to local vector
call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)
! Get corner information
call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
& PETSC_NULL_INTEGER,info)
call DAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER, &
& gxm,gym,PETSC_NULL_INTEGER,info)
! Get pointer to vector data.
call VecGetArray(localX,lx_v,lx_i,info)
! Set local loop dimensions
xe = xs+xm
ye = ys+ym
if (xs .eq. 0) then
xsm = xs-1
else
xsm = xs
endif
if (ys .eq. 0) then
ysm = ys-1
else
ysm = ys
endif
if (xe .eq. mx) then
xep = xe+1
else
xep = xe
endif
if (ye .eq. my) then
yep = ye+1
else
yep = ye
endif
! Compute local gradient contributions over the lower triangular elements
do j = ysm, ye-1
do i = xsm, xe-1
k = (j-gys)*gxm + i-gxs
v = zero
vr = zero
vt = zero
if (i .ge. 0 .and. j .ge. 0) v = lx_v(lx_i+k)
if (i .lt. mx-1 .and. j .gt. -1) vr = lx_v(lx_i+k+1)
if (i .gt. -1 .and. j .lt. my-1) vt = lx_v(lx_i+k+gxm)
dvdx = (vr-v)/hx
dvdy = (vt-v)/hy
if (i .ne. -1 .and. j .ne. -1) then
ind = k
val = - dvdx/hx - dvdy/hy - cdiv3
call VecSetValuesLocal(G,1,k,val,ADD_VALUES,info)
endif
if (i .ne. mx-1 .and. j .ne. -1) then
ind = k+1
val = dvdx/hx - cdiv3
call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
endif
if (i .ne. -1 .and. j .ne. my-1) then
ind = k+gxm
val = dvdy/hy - cdiv3
call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
endif
fquad = fquad + dvdx*dvdx + dvdy*dvdy
flin = flin - cdiv3 * (v+vr+vt)
end do
end do
! Compute local gradient contributions over the upper triangular elements
do j = ys, yep-1
do i = xs, xep-1
k = (j-gys)*gxm + i-gxs
vb = zero
vl = zero
v = zero
if (i .lt. mx .and. j .gt. 0) vb = lx_v(lx_i+k-gxm)
if (i .gt. 0 .and. j .lt. my) vl = lx_v(lx_i+k-1)
if (i .lt. mx .and. j .lt. my) v = lx_v(lx_i+k)
dvdx = (v-vl)/hx
dvdy = (v-vb)/hy
if (i .ne. mx .and. j .ne. 0) then
ind = k-gxm
val = - dvdy/hy - cdiv3
call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
endif
if (i .ne. 0 .and. j .ne. my) then
ind = k-1
val = - dvdx/hx - cdiv3
call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
endif
if (i .ne. mx .and. j .ne. my) then
ind = k
val = dvdx/hx + dvdy/hy - cdiv3
call VecSetValuesLocal(G,1,ind,val,ADD_VALUES,info)
endif
fquad = fquad + dvdx*dvdx + dvdy*dvdy
flin = flin - cdiv3*(vb + vl + v)
end do
end do
! Restore vector
call VecRestoreArray(localX,lx_v,lx_i,info)
! Assemble gradient vector
call VecAssemblyBegin(G,info)
call VecAssemblyEnd(G,info)
! Scale the gradient
area = p5*hx*hy
floc = area *(p5*fquad+flin)
call VecScale(G,area,info)
! Sum function contributions from all processes
call MPI_Allreduce(floc,f,1,MPI_DOUBLE_PRECISION,MPI_SUM, &
& MPI_COMM_WORLD,info)
call PetscLogFlops((ye-ysm)*(xe-xsm)*20+(xep-xs)*(yep-ys)*16, &
& info)
return
end
subroutine ComputeHessian(taoapp, X, H, Hpre, flag, dummy, info)
implicit none
#include "eptorsion2f.h"
TAO_APPLICATION taoapp
Vec X
Mat H,Hpre
MatStructure flag
integer dummy,info
integer i,j,k
integer col(0:4),row
integer xs,xm,gxs,gxm,ys,ym,gys,gym
PetscScalar v(0:4)
! Get local grid boundaries
call DAGetCorners(da,xs,ys,TAO_NULL_INTEGER,xm,ym, &
& TAO_NULL_INTEGER,info)
call DAGetGhostCorners(da,gxs,gys,TAO_NULL_INTEGER,gxm,gym, &
& TAO_NULL_INTEGER,info)
do j=ys,ys+ym-1
do i=xs,xs+xm-1
row = (j-gys)*gxm + (i-gxs)
k = 0
if (j .gt. gys) then
v(k) = -1.0d0
col(k) = row-gxm
k = k + 1
endif
if (i .gt. gxs) then
v(k) = -1.0d0
col(k) = row - 1
k = k +1
endif
v(k) = 4.0d0
col(k) = row
k = k + 1
if (i+1 .lt. gxs + gxm) then
v(k) = -1.0d0
col(k) = row + 1
k = k + 1
endif
if (j+1 .lt. gys + gym) then
v(k) = -1.0d0
col(k) = row + gxm
k = k + 1
endif
call MatSetValuesLocal(H,1,row,k,col,v,INSERT_VALUES,info)
enddo
enddo
! Assemble matrix
call MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY,info)
call MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY,info)
! Tell the matrix we will never add a new nonzero location to the
! matrix. If we do it will generate an error.
call MatSetOption(H,MAT_NEW_NONZERO_LOCATION_ERR,info)
call MatSetOption(H,MAT_SYMMETRIC,info)
call PetscLogFlops(9*xm*ym + 49*xm,info)
info = 0
return
end