! "$Id: plate2f.F 1.47 05/05/12 13:59:03-05:00 sarich@zorak.(none) $";
!
! Program usage: mpirun -np <proc> plate2f [all TAO options]
!
! This example demonstrates use of the TAO package to solve a bound constrained
! minimization problem. This example is based on a problem from the
! MINPACK-2 test suite. Given a rectangular 2-D domain and boundary values
! along the edges of the domain, the objective is to find the surface
! with the minimal area that satisfies the boundary conditions.
! 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
! -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction
! -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction
! -bheight <ht>, where <ht> = height of the plate
!
!/*T
! Concepts: TAO - Solving a bound constrained minimization problem
! Routines: TaoInitialize(); TaoFinalize();
! Routines: TaoCreate(); TaoDestroy();
! Routines: TaoAppSetObjectiveAndGradientRoutine();
! Routines: TaoAppSetHessianMat(); TaoAppSetHessianRoutine();
! Routines: TaoAppSetInitialSolutionVec(); TaoAppSetVariableBounds();
! Routines: TaoSetApplication(); TaoSetOptions();
! Routines: TaoApplicationCreate(); TaoSolve();
! Routines: TaoView(); TaoAppDestroy();
! Processors: n
!T*/
implicit none
#include "plate2f.h"
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
! Variable declarations
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!
! Variables:
! (common from plate2f.h):
! Nx, Ny number of processors in x- and y- directions
! mx, my number of grid points in x,y directions
! N global dimension of vector
integer info ! used to check for functions returning nonzeros
Vec x ! solution vector
Vec xl, xu ! lower and upper bounds vectorsp
integer m ! number of local elements in vector
TAO_SOLVER tao ! TAO_SOLVER solver context
TAO_APPLICATION plateapp ! PETSc application
Mat H ! Hessian matrix
ISLocalToGlobalMapping isltog ! local to global mapping object
PetscTruth flg
external FormFunctionGradient
external FormHessian
external MSA_BoundaryConditions
external MSA_Plate
external MSA_InitialPoint
! Initialize Tao
call PetscInitialize(PETSC_NULL_CHARACTER,info)
call TaoInitialize(PETSC_NULL_CHARACTER,info)
! Specify default dimensions of the problem
mx = 10
my = 10
bheight = 0.1
! Check for any command line arguments that override defaults
call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-mx",mx,flg,info)
call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-my",my,flg,info)
bmx = mx/2
bmy = my/2
call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmx",bmx,flg,info)
call PetscOptionsGetInt(TAO_NULL_CHARACTER,"-bmy",bmy,flg,info)
call PetscOptionsGetReal(TAO_NULL_CHARACTER,"-bheight",bheight, &
& flg,info)
! Calculate any derived values from parameters
N = mx*my
! Let Petsc determine the dimensions of the local vectors
Nx = PETSC_DECIDE
NY = PETSC_DECIDE
! A two dimensional distributed array will help define this problem, which
! derives from an elliptic PDE on a two-dimensional domain. From the
! distributed array, create the 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)
! Extract global and local vectors from DA; The local vectors are
! used solely as work space for the evaluation of the function,
! gradient, and Hessian. Duplicate for remaining vectors that are
! the same types.
call DACreateGlobalVector(da,x,info)
call DACreateLocalVector(da,localX,info)
call VecDuplicate(localX,localV,info)
! Create a matrix data structure to store the Hessian.
! Here we (optionally) also associate the local numbering scheme
! with the matrix so that later we can use local indices for matrix
! assembly
call VecGetLocalSize(x,m,info)
call MatCreateMPIAIJ(MPI_COMM_WORLD,m,m,N,N,7,PETSC_NULL_INTEGER, &
& 3,PETSC_NULL_INTEGER,H,info)
call MatSetOption(H,MAT_SYMMETRIC,info)
call DAGetISLocalToGlobalMapping(da,isltog,info)
call MatSetLocalToGlobalMapping(H,isltog,info)
! The Tao code begins here
! Create TAO solver and set desired solution method.
! This problems uses bounded variables, so the
! method must either be 'tao_tron' or 'tao_blmvm'
call TaoCreate(MPI_COMM_WORLD,'tao_blmvm',tao,info)
call TaoApplicationCreate(MPI_COMM_WORLD,plateapp,info)
! Set minimization function and gradient, hessian evaluation functions
! TaoAppSetObjectiveAndGradientRoutine is shortened to 31 chars to comply with some compilers
call TaoAppSetObjectiveAndGradientRo(plateapp, &
& FormFunctionGradient,PETSC_NULL_OBJECT,info)
call TaoAppSetHessianMat(plateapp,H,H,info)
call TaoAppSetHessianRoutine(plateapp,FormHessian, &
& PETSC_NULL_OBJECT, info)
! Set Variable bounds
call MSA_BoundaryConditions(info)
call VecDuplicate(x,xl,info)
call VecDuplicate(x,xu,info)
call MSA_Plate(xl,xu,info)
call TaoAppSetVariableBounds(plateapp,xl,xu,info)
! Set the initial solution guess
call MSA_InitialPoint(x, info)
call TaoAppSetInitialSolutionVec(plateapp,x,info)
! Now that the PETSc application is all set, attach to tao context
call TaoSetApplication(tao,plateapp,info)
! Check for any tao command line options
call TaoSetOptions(plateapp,tao,info)
! Solve the application
call TaoSolveApplication(plateapp,tao,info)
! View TAO solver information
! call TaoView(tao,info)
! Free TAO data structures
call TaoDestroy(tao,info)
call TaoAppDestroy(plateapp,info)
! Free PETSc data structures
call VecDestroy(x,info)
call VecDestroy(xl,info)
call VecDestroy(xu,info)
call VecDestroy(Top,info)
call VecDestroy(Bottom,info)
call VecDestroy(Left,info)
call VecDestroy(Right,info)
call MatDestroy(H,info)
call VecDestroy(localX,info)
call VecDestroy(localV,info)
call DADestroy(da,info)
! Finalize TAO
call TaoFinalize(info)
call PetscFinalize(info)
end
! ---------------------------------------------------------------------
!
! FormFunctionGradient - Evaluates function f(X).
!
! Input Parameters:
! tao - the TAO_SOLVER context
! X - the input vector
! dummy - optional user-defined context, as set by TaoSetFunction()
! (not used here)
!
! Output Parameters:
! fcn - the newly evaluated function
! G - the gradient vector
! info - error code
!
subroutine FormFunctionGradient(tao,X,fcn,G,dummy,info)
implicit none
! da, localX, localG, Top, Bottom, Left, Right defined in plate2f.h
#include "plate2f.h"
! Input/output variables
TAO_SOLVER tao
PetscScalar fcn
Vec X, G
integer dummy, info
integer i,j,row
integer xs, xm, gxs, gxm, ys, ym, gys, gym
PetscScalar ft,zero,hx,hy,hydhx,hxdhy,area,rhx,rhy
PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8
PetscScalar df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc
PetscScalar xc,xl,xr,xt,xb,xlt,xrb
! 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 g_v(0:1),x_v(0:1), top_v(0:1)
PetscScalar left_v(0:1), right_v(0:1),bottom_v(0:1)
PetscOffset g_i,left_i,right_i,bottom_i,top_i
PetscOffset x_i
ft = 0.0d0
zero = 0.0d0
hx = 1.0d0/(mx + 1)
hy = 1.0d0/(my + 1)
hydhx = hy/hx
hxdhy = hx/hy
area = 0.5d0 * hx * hy
rhx = mx + 1.0d0
rhy = my + 1.0d0
! Get local mesh boundaries
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)
! Scatter ghost points to local vector
call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)
! Initialize the vector to zero
call VecSet(localV,zero,info)
! Get arrays to vector data (See note above about using VecGetArray in Fortran)
call VecGetArray(localX,x_v,x_i,info)
call VecGetArray(localV,g_v,g_i,info)
call VecGetArray(Top,top_v,top_i,info)
call VecGetArray(Bottom,bottom_v,bottom_i,info)
call VecGetArray(Left,left_v,left_i,info)
call VecGetArray(Right,right_v,right_i,info)
! Compute function over the locally owned part of the mesh
do j = ys,ys+ym-1
do i = xs,xs+xm-1
row = (j-gys)*gxm + (i-gxs)
xc = x_v(row+x_i)
xt = xc
xb = xc
xr = xc
xl = xc
xrb = xc
xlt = xc
if (i .eq. 0) then !left side
xl = left_v(j - ys + 1 + left_i)
xlt = left_v(j - ys + 2 + left_i)
else
xl = x_v(row - 1 + x_i)
endif
if (j .eq. 0) then !bottom side
xb = bottom_v(i - xs + 1 + bottom_i)
xrb = bottom_v(i - xs + 2 + bottom_i)
else
xb = x_v(row - gxm + x_i)
endif
if (i + 1 .eq. gxs + gxm) then !right side
xr = right_v(j - ys + 1 + right_i)
xrb = right_v(j - ys + right_i)
else
xr = x_v(row + 1 + x_i)
endif
if (j + 1 .eq. gys + gym) then !top side
xt = top_v(i - xs + 1 + top_i)
xlt = top_v(i - xs + top_i)
else
xt = x_v(row + gxm + x_i)
endif
if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then
xlt = x_v(row - 1 + gxm + x_i)
endif
if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then
xrb = x_v(row + 1 - gxm + x_i)
endif
d1 = xc-xl
d2 = xc-xr
d3 = xc-xt
d4 = xc-xb
d5 = xr-xrb
d6 = xrb-xb
d7 = xlt-xl
d8 = xt-xlt
df1dxc = d1 * hydhx
df2dxc = d1 * hydhx + d4 * hxdhy
df3dxc = d3 * hxdhy
df4dxc = d2 * hydhx + d3 * hxdhy
df5dxc = d2 * hydhx
df6dxc = d4 * hxdhy
d1 = d1 * rhx
d2 = d2 * rhx
d3 = d3 * rhy
d4 = d4 * rhy
d5 = d5 * rhy
d6 = d6 * rhx
d7 = d7 * rhy
d8 = d8 * rhx
f1 = sqrt(1.0d0 + d1*d1 + d7*d7)
f2 = sqrt(1.0d0 + d1*d1 + d4*d4)
f3 = sqrt(1.0d0 + d3*d3 + d8*d8)
f4 = sqrt(1.0d0 + d3*d3 + d2*d2)
f5 = sqrt(1.0d0 + d2*d2 + d5*d5)
f6 = sqrt(1.0d0 + d4*d4 + d6*d6)
ft = ft + f2 + f4
df1dxc = df1dxc / f1
df2dxc = df2dxc / f2
df3dxc = df3dxc / f3
df4dxc = df4dxc / f4
df5dxc = df5dxc / f5
df6dxc = df6dxc / f6
g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + &
& df5dxc + df6dxc)
enddo
enddo
! Compute triangular areas along the border of the domain.
if (xs .eq. 0) then ! left side
do j=ys,ys+ym-1
d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) &
& * rhy
d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) &
& * rhx
ft = ft + sqrt(1.0d0 + d3*d3 + d2*d2)
enddo
endif
if (ys .eq. 0) then !bottom side
do i=xs,xs+xm-1
d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) &
& * rhx
d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy
ft = ft + sqrt(1.0 + d3*d3 + d2*d2)
enddo
endif
if (xs + xm .eq. mx) then ! right side
do j=ys,ys+ym-1
d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx
d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy
ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
enddo
endif
if (ys + ym .eq. my) then
do i=xs,xs+xm-1
d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy
d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx
ft = ft + sqrt(1.0d0 + d1*d1 + d4*d4)
enddo
endif
if ((ys .eq. 0) .and. (xs .eq. 0)) then
d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy
d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx
ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
endif
if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then
d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy
d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx
ft = ft + sqrt(1.0d0 + d1*d1 + d2*d2)
endif
ft = ft * area
call MPI_Allreduce(ft,fcn,1,MPI_DOUBLE_PRECISION, &
& MPI_SUM,MPI_COMM_WORLD,info)
! Restore vectors
call VecRestoreArray(localX,x_v,x_i,info)
call VecRestoreArray(localV,g_v,g_i,info)
call VecRestoreArray(Left,left_v,left_i,info)
call VecRestoreArray(Top,top_v,top_i,info)
call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
call VecRestoreArray(Right,right_v,right_i,info)
! Scatter values to global vector
call DALocalToGlobal(da,localV,INSERT_VALUES,G,info)
call PetscLogFlops(70*xm*ym,info)
return
end !FormFunctionGradient
! ----------------------------------------------------------------------------
!
!/*
! FormHessian - Evaluates Hessian matrix.
!
! Input Parameters:
!. tao - the TAO_SOLVER context
!. X - input vector
!. dummy - not used
!
! Output Parameters:
!. Hessian - Hessian matrix
!. Hpc - optionally different preconditioning matrix
!. flag - flag indicating matrix structure
!
! Notes:
! Due to mesh point reordering with DAs, we must always work
! with the local mesh points, and then transform them to the new
! global numbering with the local-to-global mapping. We cannot work
! directly with the global numbers for the original uniprocessor mesh!
!
! Two methods are available for imposing this transformation
! when setting matrix entries:
! (A) MatSetValuesLocal(), using the local ordering (including
! ghost points!)
! - Do the following two steps once, before calling TaoSolve()
! - Use DAGetISLocalToGlobalMapping() to extract the
! local-to-global map from the DA
! - Associate this map with the matrix by calling
! MatSetLocalToGlobalMapping()
! - Then set matrix entries using the local ordering
! by calling MatSetValuesLocal()
! (B) MatSetValues(), using the global ordering
! - Use DAGetGlobalIndices() to extract the local-to-global map
! - Then apply this map explicitly yourself
! - Set matrix entries using the global ordering by calling
! MatSetValues()
! Option (A) seems cleaner/easier in many cases, and is the procedure
! used in this example.
*/
subroutine FormHessian(tao, X, Hessian, Hpc, flg, dummy, info)
implicit none
! da,Top,Left,Right,Bottom,mx,my,localX defined in plate2f.h
#include "plate2f.h"
TAO_SOLVER tao
Vec X
Mat Hessian,Hpc
MatStructure flg
integer dummy,info
integer i,j,k,row
integer xs,xm,gxs,gxm,ys,ym,gys,gym,col(0:6)
PetscScalar hx,hy,hydhx,hxdhy,rhx,rhy
PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8
PetscScalar xc,xl,xr,xt,xb,xlt,xrb
PetscScalar hl,hr,ht,hb,hc,htl,hbr
! 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 right_v(0:1),left_v(0:1),bottom_v(0:1),top_v(0:1)
PetscScalar x_v(0:1)
PetscOffset x_i,right_i,left_i,bottom_i,top_i
PetscScalar v(0:6)
PetscTruth assembled
! Set various matrix options
call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES,info)
call MatSetOption(Hessian,MAT_COLUMNS_SORTED,info)
call MatSetOption(Hessian,MAT_ROWS_SORTED,info)
! Get local mesh boundaries
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)
! Scatter ghost points to local vectors
call DAGlobalToLocalBegin(da,X,INSERT_VALUES,localX,info)
call DAGlobalToLocalEnd(da,X,INSERT_VALUES,localX,info)
! Get pointers to vector data (see note on Fortran arrays above)
call VecGetArray(localX,x_v,x_i,info)
call VecGetArray(Top,top_v,top_i,info)
call VecGetArray(Bottom,bottom_v,bottom_i,info)
call VecGetArray(Left,left_v,left_i,info)
call VecGetArray(Right,right_v,right_i,info)
! Initialize matrix entries to zero
call MatAssembled(Hessian,assembled,info)
if (assembled .eq. PETSC_TRUE) call MatZeroEntries(Hessian,info)
rhx = mx + 1.0
rhy = my + 1.0
hx = 1.0/rhx
hy = 1.0/rhy
hydhx = hy/hx
hxdhy = hx/hy
! compute Hessian over the locally owned part of the mesh
do i=xs,xs+xm-1
do j=ys,ys+ym-1
row = (j-gys)*gxm + (i-gxs)
xc = x_v(row + x_i)
xt = xc
xb = xc
xr = xc
xl = xc
xrb = xc
xlt = xc
if (i .eq. gxs) then ! Left side
xl = left_v(left_i + j - ys + 1)
xlt = left_v(left_i + j - ys + 2)
else
xl = x_v(x_i + row -1 )
endif
if (j .eq. gys) then ! bottom side
xb = bottom_v(bottom_i + i - xs + 1)
xrb = bottom_v(bottom_i + i - xs + 2)
else
xb = x_v(x_i + row - gxm)
endif
if (i+1 .eq. gxs + gxm) then !right side
xr = right_v(right_i + j - ys + 1)
xrb = right_v(right_i + j - ys)
else
xr = x_v(x_i + row + 1)
endif
if (j+1 .eq. gym+gys) then !top side
xt = top_v(top_i +i - xs + 1)
xlt = top_v(top_i + i - xs)
else
xt = x_v(x_i + row + gxm)
endif
if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then
xlt = x_v(x_i + row - 1 + gxm)
endif
if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then
xrb = x_v(x_i + row + 1 - gxm)
endif
d1 = (xc-xl)*rhx
d2 = (xc-xr)*rhx
d3 = (xc-xt)*rhy
d4 = (xc-xb)*rhy
d5 = (xrb-xr)*rhy
d6 = (xrb-xb)*rhx
d7 = (xlt-xl)*rhy
d8 = (xlt-xt)*rhx
f1 = sqrt( 1.0d0 + d1*d1 + d7*d7)
f2 = sqrt( 1.0d0 + d1*d1 + d4*d4)
f3 = sqrt( 1.0d0 + d3*d3 + d8*d8)
f4 = sqrt( 1.0d0 + d3*d3 + d2*d2)
f5 = sqrt( 1.0d0 + d2*d2 + d5*d5)
f6 = sqrt( 1.0d0 + d4*d4 + d6*d6)
hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ &
& (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2)
hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ &
& (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4)
ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ &
& (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4)
hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ &
& (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2)
hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6)
htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3)
hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + &
& hxdhy*(1.0+d8*d8)/(f3*f3*f3) + &
& hydhx*(1.0+d5*d5)/(f5*f5*f5) + &
& hxdhy*(1.0+d6*d6)/(f6*f6*f6) + &
& (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- &
& 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ &
& hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4)
hl = hl * 0.5
hr = hr * 0.5
ht = ht * 0.5
hb = hb * 0.5
hbr = hbr * 0.5
htl = htl * 0.5
hc = hc * 0.5
k = 0
if (j .gt. 0) then
v(k) = hb
col(k) = row - gxm
k=k+1
endif
if ((j .gt. 0) .and. (i .lt. mx-1)) then
v(k) = hbr
col(k) = row-gxm+1
k=k+1
endif
if (i .gt. 0) then
v(k) = hl
col(k) = row - 1
k = k+1
endif
v(k) = hc
col(k) = row
k=k+1
if (i .lt. mx-1) then
v(k) = hr
col(k) = row + 1
k=k+1
endif
if ((i .gt. 0) .and. (j .lt. my-1)) then
v(k) = htl
col(k) = row + gxm - 1
k=k+1
endif
if (j .lt. my-1) then
v(k) = ht
col(k) = row + gxm
k=k+1
endif
! Set matrix values using local numbering, defined earlier in main routine
call MatSetValuesLocal(Hessian,1,row,k,col,v,INSERT_VALUES, &
& info)
enddo
enddo
! restore vectors
call VecRestoreArray(localX,x_v,x_i,info)
call VecRestoreArray(Left,left_v,left_i,info)
call VecRestoreArray(Right,right_v,right_i,info)
call VecRestoreArray(Top,top_v,top_i,info)
call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
! Assemble the matrix
call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,info)
call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,info)
call PetscLogFlops(199*xm*ym,info)
return
end
! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h
! ----------------------------------------------------------------------------
!
!/*
! MSA_BoundaryConditions - calculates the boundary conditions for the region
!
!
!*/
subroutine MSA_BoundaryConditions(info)
implicit none
! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy defined in plate2f.h
#include "plate2f.h"
integer i,j,k,limit,info,maxits
integer xs, xm, gxs, gxm, ys, ym, gys, gym
integer bsize, lsize, tsize, rsize
PetscScalar one, two, three, tol
PetscScalar scl,fnorm,det,xt,yt,hx,hy
PetscScalar u1,u2,nf1,nf2,njac11,njac12,njac21,njac22
PetscScalar b, t, l, r
PetscScalar boundary_v(0:1)
PetscOffset boundary_i
logical exitloop
TaoTruth flg
limit=0
maxits = 5
tol=1e-10
b=-0.5d0
t= 0.5d0
l=-0.5d0
r= 0.5d0
xt=0
yt=0
one=1.0d0
two=2.0d0
three=3.0d0
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)
bsize = xm + 2
lsize = ym + 2
rsize = ym + 2
tsize = xm + 2
call VecCreateMPI(MPI_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,info)
call VecCreateMPI(MPI_COMM_WORLD,tsize,PETSC_DECIDE,Top,info)
call VecCreateMPI(MPI_COMM_WORLD,lsize,PETSC_DECIDE,Left,info)
call VecCreateMPI(MPI_COMM_WORLD,rsize,PETSC_DECIDE,Right,info)
hx= (r-l)/(mx+1)
hy= (t-b)/(my+1)
do j=0,3
if (j.eq.0) then
yt=b
xt=l+hx*xs
limit=bsize
call VecGetArray(Bottom,boundary_v,boundary_i,info)
elseif (j.eq.1) then
yt=t
xt=l+hx*xs
limit=tsize
call VecGetArray(Top,boundary_v,boundary_i,info)
elseif (j.eq.2) then
yt=b+hy*ys
xt=l
limit=lsize
call VecGetArray(Left,boundary_v,boundary_i,info)
elseif (j.eq.3) then
yt=b+hy*ys
xt=r
limit=rsize
call VecGetArray(Right,boundary_v,boundary_i,info)
endif
do i=0,limit-1
u1=xt
u2=-yt
k = 0
exitloop = .false.
do while (k .lt. maxits .and. (.not. exitloop) )
nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt
nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt
fnorm=sqrt(nf1*nf1+nf2*nf2)
if (fnorm .gt. tol) then
njac11=one+u2*u2-u1*u1
njac12=two*u1*u2
njac21=-two*u1*u2
njac22=-one - u1*u1 + u2*u2
det = njac11*njac22-njac21*njac12
u1 = u1-(njac22*nf1-njac12*nf2)/det
u2 = u2-(njac11*nf2-njac21*nf1)/det
else
exitloop = .true.
endif
k=k+1
enddo
boundary_v(i + boundary_i) = u1*u1-u2*u2
if ((j .eq. 0) .or. (j .eq. 1)) then
xt = xt + hx
else
yt = yt + hy
endif
enddo
if (j.eq.0) then
call VecRestoreArray(Bottom,boundary_v,boundary_i,info)
elseif (j.eq.1) then
call VecRestoreArray(Top,boundary_v,boundary_i,info)
elseif (j.eq.2) then
call VecRestoreArray(Left,boundary_v,boundary_i,info)
elseif (j.eq.3) then
call VecRestoreArray(Right,boundary_v,boundary_i,info)
endif
enddo
! Scale the boundary if desired
call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-bottom", &
& scl,flg,info)
if (flg .eq. PETSC_TRUE) then
call VecScale(scl,Bottom,info)
endif
call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-top", &
& scl,flg,info)
if (flg .eq. PETSC_TRUE) then
call VecScale(scl,Top,info)
endif
call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-right", &
& scl,flg,info)
if (flg .eq. PETSC_TRUE) then
call VecScale(scl,Right,info)
endif
call PetscOptionsGetReal(PETSC_NULL_CHARACTER,"-left", &
& scl,flg,info)
if (flg .eq. PETSC_TRUE) then
call VecScale(scl,Left,info)
endif
return
end
! ----------------------------------------------------------------------------
!
!/*
! MSA_Plate - Calculates an obstacle for surface to stretch over
!
! Output Parameter:
!. xl - lower bound vector
!. xu - upper bound vector
!
!*/
subroutine MSA_Plate(xl,xu,info)
implicit none
! mx,my,bmx,bmy,da,bheight defined in plate2f.h
#include "plate2f.h"
Vec xl,xu
integer i,j,row,info
integer xs, xm, ys, ym
PetscScalar lb,ub
! 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 xl_v(0:1)
PetscOffset xl_i
lb = -1.0d300
ub = 1.0d300
if (bmy .lt. 0) bmy = 0
if (bmy .gt. my) bmy = my
if (bmx .lt. 0) bmx = 0
if (bmx .gt. mx) bmx = mx
call DAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
& PETSC_NULL_INTEGER,info)
call VecSet(xl,lb,info)
call VecSet(xu,ub,info)
call VecGetArray(xl,xl_v,xl_i,info)
do i=xs,xs+xm-1
do j=ys,ys+ym-1
row=(j-ys)*xm + (i-xs)
if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. &
& j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then
xl_v(xl_i+row) = bheight
endif
enddo
enddo
call VecRestoreArray(xl,xl_v,xl_i,info)
return
end
! ----------------------------------------------------------------------------
!
!/*
! MSA_InitialPoint - Calculates an obstacle for surface to stretch over
!
! Output Parameter:
!. X - vector for initial guess
!
!*/
subroutine MSA_InitialPoint(X, info)
implicit none
! mx,my,localX,da,Top,Left,Bottom,Right defined in plate2f.h
#include "plate2f.h"
Vec X
integer start,i,j,info
integer row,xs,xm,gxs,gxm,ys,ym,gys,gym
PetscScalar zero, np5
! PETSc's VecGetArray acts differently in Fortran than it does in C.
! Calling VecGetArray((Vec) X, (PetscScalar) x_array(0:1), (integer) 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 left_v(0:1), right_v(0:1), top_v(0:1)
PetscScalar x_v(0:1), bottom_v(0:1)
PetscOffset left_i, right_i, top_i, bottom_i, x_i
PetscTruth flg
PetscRandom rctx
zero = 0.0d0
np5 = -0.5d0
call PetscOptionsGetInt(PETSC_NULL_CHARACTER,"-start", &
& start,flg,info)
if ((flg .eq. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable
call VecSet(X,zero,info)
elseif ((flg .eq. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5
call PetscRandomCreate(MPI_COMM_WORLD,RANDOM_DEFAULT,rctx,info)
do i=0,start-1
call VecSetRandom(X,rctx,info)
enddo
call PetscRandomDestroy(rctx,info)
call VecShift(X,np5,info)
else ! average of boundary conditions
! Get Local mesh boundaries
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 pointers to vector data
call VecGetArray(Top,top_v,top_i,info)
call VecGetArray(Bottom,bottom_v,bottom_i,info)
call VecGetArray(Left,left_v,left_i,info)
call VecGetArray(Right,right_v,right_i,info)
call VecGetArray(localX,x_v,x_i,info)
! Perform local computations
do j=ys,ys+ym-1
do i=xs,xs+xm-1
row = (j-gys)*gxm + (i-gxs)
x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my &
& + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + &
& (i+1)*left_v(left_i+j-ys+1)/mx + &
& (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5
enddo
enddo
! Restore vectors
call VecRestoreArray(localX,x_v,x_i,info)
call VecRestoreArray(Left,left_v,left_i,info)
call VecRestoreArray(Top,top_v,top_i,info)
call VecRestoreArray(Bottom,bottom_v,bottom_i,info)
call VecRestoreArray(Right,right_v,right_i,info)
call DALocalToGlobal(da,localX,INSERT_VALUES,X,info)
endif
return
end