1*c4762a1bSJed Brown! Program usage: mpiexec -n <proc> plate2f [all TAO options] 2*c4762a1bSJed Brown! 3*c4762a1bSJed Brown! This example demonstrates use of the TAO package to solve a bound constrained 4*c4762a1bSJed Brown! minimization problem. This example is based on a problem from the 5*c4762a1bSJed Brown! MINPACK-2 test suite. Given a rectangular 2-D domain and boundary values 6*c4762a1bSJed Brown! along the edges of the domain, the objective is to find the surface 7*c4762a1bSJed Brown! with the minimal area that satisfies the boundary conditions. 8*c4762a1bSJed Brown! The command line options are: 9*c4762a1bSJed Brown! -mx <xg>, where <xg> = number of grid points in the 1st coordinate direction 10*c4762a1bSJed Brown! -my <yg>, where <yg> = number of grid points in the 2nd coordinate direction 11*c4762a1bSJed Brown! -bmx <bxg>, where <bxg> = number of grid points under plate in 1st direction 12*c4762a1bSJed Brown! -bmy <byg>, where <byg> = number of grid points under plate in 2nd direction 13*c4762a1bSJed Brown! -bheight <ht>, where <ht> = height of the plate 14*c4762a1bSJed Brown! 15*c4762a1bSJed Brown!!/*T 16*c4762a1bSJed Brown! Concepts: TAO^Solving a bound constrained minimization problem 17*c4762a1bSJed Brown! Routines: TaoCreate(); 18*c4762a1bSJed Brown! Routines: TaoSetType(); TaoSetObjectiveAndGradientRoutine(); 19*c4762a1bSJed Brown! Routines: TaoSetHessianRoutine(); 20*c4762a1bSJed Brown! Routines: TaoSetVariableBoundsRoutine(); 21*c4762a1bSJed Brown! Routines: TaoSetInitialVector(); 22*c4762a1bSJed Brown! Routines: TaoSetFromOptions(); 23*c4762a1bSJed Brown! Routines: TaoSolve(); 24*c4762a1bSJed Brown! Routines: TaoDestroy(); 25*c4762a1bSJed Brown! Processors: n 26*c4762a1bSJed Brown!T*/ 27*c4762a1bSJed Brown 28*c4762a1bSJed Brown module mymodule 29*c4762a1bSJed Brown#include "petsc/finclude/petscdmda.h" 30*c4762a1bSJed Brown#include "petsc/finclude/petsctao.h" 31*c4762a1bSJed Brown use petscdmda 32*c4762a1bSJed Brown use petsctao 33*c4762a1bSJed Brown 34*c4762a1bSJed Brown Vec localX, localV 35*c4762a1bSJed Brown Vec Top, Left 36*c4762a1bSJed Brown Vec Right, Bottom 37*c4762a1bSJed Brown DM dm 38*c4762a1bSJed Brown PetscReal bheight 39*c4762a1bSJed Brown PetscInt bmx, bmy 40*c4762a1bSJed Brown PetscInt mx, my, Nx, Ny, N 41*c4762a1bSJed Brown end module 42*c4762a1bSJed Brown 43*c4762a1bSJed Brown 44*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 45*c4762a1bSJed Brown! Variable declarations 46*c4762a1bSJed Brown! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 47*c4762a1bSJed Brown! 48*c4762a1bSJed Brown! Variables: 49*c4762a1bSJed Brown! (common from plate2f.h): 50*c4762a1bSJed Brown! Nx, Ny number of processors in x- and y- directions 51*c4762a1bSJed Brown! mx, my number of grid points in x,y directions 52*c4762a1bSJed Brown! N global dimension of vector 53*c4762a1bSJed Brown use mymodule 54*c4762a1bSJed Brown implicit none 55*c4762a1bSJed Brown 56*c4762a1bSJed Brown PetscErrorCode ierr ! used to check for functions returning nonzeros 57*c4762a1bSJed Brown Vec x ! solution vector 58*c4762a1bSJed Brown PetscInt m ! number of local elements in vector 59*c4762a1bSJed Brown Tao tao ! Tao solver context 60*c4762a1bSJed Brown Mat H ! Hessian matrix 61*c4762a1bSJed Brown ISLocalToGlobalMapping isltog ! local to global mapping object 62*c4762a1bSJed Brown PetscBool flg 63*c4762a1bSJed Brown PetscInt i1,i3,i7 64*c4762a1bSJed Brown 65*c4762a1bSJed Brown 66*c4762a1bSJed Brown external FormFunctionGradient 67*c4762a1bSJed Brown external FormHessian 68*c4762a1bSJed Brown external MSA_BoundaryConditions 69*c4762a1bSJed Brown external MSA_Plate 70*c4762a1bSJed Brown external MSA_InitialPoint 71*c4762a1bSJed Brown! Initialize Tao 72*c4762a1bSJed Brown 73*c4762a1bSJed Brown i1=1 74*c4762a1bSJed Brown i3=3 75*c4762a1bSJed Brown i7=7 76*c4762a1bSJed Brown 77*c4762a1bSJed Brown 78*c4762a1bSJed Brown call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 79*c4762a1bSJed Brown if (ierr .ne. 0) then 80*c4762a1bSJed Brown print*,'Unable to initialize PETSc' 81*c4762a1bSJed Brown stop 82*c4762a1bSJed Brown endif 83*c4762a1bSJed Brown 84*c4762a1bSJed Brown! Specify default dimensions of the problem 85*c4762a1bSJed Brown mx = 10 86*c4762a1bSJed Brown my = 10 87*c4762a1bSJed Brown bheight = 0.1 88*c4762a1bSJed Brown 89*c4762a1bSJed Brown! Check for any command line arguments that override defaults 90*c4762a1bSJed Brown 91*c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 92*c4762a1bSJed Brown & '-mx',mx,flg,ierr) 93*c4762a1bSJed Brown CHKERRA(ierr) 94*c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 95*c4762a1bSJed Brown & '-my',my,flg,ierr) 96*c4762a1bSJed Brown CHKERRA(ierr) 97*c4762a1bSJed Brown 98*c4762a1bSJed Brown bmx = mx/2 99*c4762a1bSJed Brown bmy = my/2 100*c4762a1bSJed Brown 101*c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 102*c4762a1bSJed Brown & '-bmx',bmx,flg,ierr) 103*c4762a1bSJed Brown CHKERRA(ierr) 104*c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 105*c4762a1bSJed Brown & '-bmy',bmy,flg,ierr) 106*c4762a1bSJed Brown CHKERRA(ierr) 107*c4762a1bSJed Brown call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 108*c4762a1bSJed Brown & '-bheight',bheight,flg,ierr) 109*c4762a1bSJed Brown CHKERRA(ierr) 110*c4762a1bSJed Brown 111*c4762a1bSJed Brown 112*c4762a1bSJed Brown! Calculate any derived values from parameters 113*c4762a1bSJed Brown N = mx*my 114*c4762a1bSJed Brown 115*c4762a1bSJed Brown! Let Petsc determine the dimensions of the local vectors 116*c4762a1bSJed Brown Nx = PETSC_DECIDE 117*c4762a1bSJed Brown NY = PETSC_DECIDE 118*c4762a1bSJed Brown 119*c4762a1bSJed Brown! A two dimensional distributed array will help define this problem, which 120*c4762a1bSJed Brown! derives from an elliptic PDE on a two-dimensional domain. From the 121*c4762a1bSJed Brown! distributed array, create the vectors 122*c4762a1bSJed Brown 123*c4762a1bSJed Brown call DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE, & 124*c4762a1bSJed Brown & DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, & 125*c4762a1bSJed Brown & mx,my,Nx,Ny,i1,i1,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & 126*c4762a1bSJed Brown & dm,ierr) 127*c4762a1bSJed Brown CHKERRA(ierr) 128*c4762a1bSJed Brown call DMSetFromOptions(dm,ierr) 129*c4762a1bSJed Brown call DMSetUp(dm,ierr) 130*c4762a1bSJed Brown 131*c4762a1bSJed Brown! Extract global and local vectors from DM; The local vectors are 132*c4762a1bSJed Brown! used solely as work space for the evaluation of the function, 133*c4762a1bSJed Brown! gradient, and Hessian. Duplicate for remaining vectors that are 134*c4762a1bSJed Brown! the same types. 135*c4762a1bSJed Brown 136*c4762a1bSJed Brown call DMCreateGlobalVector(dm,x,ierr) 137*c4762a1bSJed Brown CHKERRA(ierr) 138*c4762a1bSJed Brown call DMCreateLocalVector(dm,localX,ierr) 139*c4762a1bSJed Brown CHKERRA(ierr) 140*c4762a1bSJed Brown call VecDuplicate(localX,localV,ierr) 141*c4762a1bSJed Brown CHKERRA(ierr) 142*c4762a1bSJed Brown 143*c4762a1bSJed Brown! Create a matrix data structure to store the Hessian. 144*c4762a1bSJed Brown! Here we (optionally) also associate the local numbering scheme 145*c4762a1bSJed Brown! with the matrix so that later we can use local indices for matrix 146*c4762a1bSJed Brown! assembly 147*c4762a1bSJed Brown 148*c4762a1bSJed Brown call VecGetLocalSize(x,m,ierr) 149*c4762a1bSJed Brown CHKERRA(ierr) 150*c4762a1bSJed Brown call MatCreateAIJ(PETSC_COMM_WORLD,m,m,N,N,i7,PETSC_NULL_INTEGER, & 151*c4762a1bSJed Brown & i3,PETSC_NULL_INTEGER,H,ierr) 152*c4762a1bSJed Brown CHKERRA(ierr) 153*c4762a1bSJed Brown 154*c4762a1bSJed Brown call MatSetOption(H,MAT_SYMMETRIC,PETSC_TRUE,ierr) 155*c4762a1bSJed Brown CHKERRA(ierr) 156*c4762a1bSJed Brown call DMGetLocalToGlobalMapping(dm,isltog,ierr) 157*c4762a1bSJed Brown CHKERRA(ierr) 158*c4762a1bSJed Brown call MatSetLocalToGlobalMapping(H,isltog,isltog,ierr) 159*c4762a1bSJed Brown CHKERRA(ierr) 160*c4762a1bSJed Brown 161*c4762a1bSJed Brown 162*c4762a1bSJed Brown! The Tao code begins here 163*c4762a1bSJed Brown! Create TAO solver and set desired solution method. 164*c4762a1bSJed Brown! This problems uses bounded variables, so the 165*c4762a1bSJed Brown! method must either be 'tao_tron' or 'tao_blmvm' 166*c4762a1bSJed Brown 167*c4762a1bSJed Brown call TaoCreate(PETSC_COMM_WORLD,tao,ierr) 168*c4762a1bSJed Brown CHKERRA(ierr) 169*c4762a1bSJed Brown call TaoSetType(tao,TAOBLMVM,ierr) 170*c4762a1bSJed Brown CHKERRA(ierr) 171*c4762a1bSJed Brown 172*c4762a1bSJed Brown! Set minimization function and gradient, hessian evaluation functions 173*c4762a1bSJed Brown 174*c4762a1bSJed Brown call TaoSetObjectiveAndGradientRoutine(tao, & 175*c4762a1bSJed Brown & FormFunctionGradient,0,ierr) 176*c4762a1bSJed Brown CHKERRA(ierr) 177*c4762a1bSJed Brown 178*c4762a1bSJed Brown call TaoSetHessianRoutine(tao,H,H,FormHessian, & 179*c4762a1bSJed Brown & 0, ierr) 180*c4762a1bSJed Brown CHKERRA(ierr) 181*c4762a1bSJed Brown 182*c4762a1bSJed Brown! Set Variable bounds 183*c4762a1bSJed Brown call MSA_BoundaryConditions(ierr) 184*c4762a1bSJed Brown CHKERRA(ierr) 185*c4762a1bSJed Brown call TaoSetVariableBoundsRoutine(tao,MSA_Plate, & 186*c4762a1bSJed Brown & 0,ierr) 187*c4762a1bSJed Brown CHKERRA(ierr) 188*c4762a1bSJed Brown 189*c4762a1bSJed Brown! Set the initial solution guess 190*c4762a1bSJed Brown call MSA_InitialPoint(x, ierr) 191*c4762a1bSJed Brown CHKERRA(ierr) 192*c4762a1bSJed Brown call TaoSetInitialVector(tao,x,ierr) 193*c4762a1bSJed Brown CHKERRA(ierr) 194*c4762a1bSJed Brown 195*c4762a1bSJed Brown! Check for any tao command line options 196*c4762a1bSJed Brown call TaoSetFromOptions(tao,ierr) 197*c4762a1bSJed Brown CHKERRA(ierr) 198*c4762a1bSJed Brown 199*c4762a1bSJed Brown! Solve the application 200*c4762a1bSJed Brown call TaoSolve(tao,ierr) 201*c4762a1bSJed Brown CHKERRA(ierr) 202*c4762a1bSJed Brown 203*c4762a1bSJed Brown! Free TAO data structures 204*c4762a1bSJed Brown call TaoDestroy(tao,ierr) 205*c4762a1bSJed Brown CHKERRA(ierr) 206*c4762a1bSJed Brown 207*c4762a1bSJed Brown! Free PETSc data structures 208*c4762a1bSJed Brown call VecDestroy(x,ierr) 209*c4762a1bSJed Brown call VecDestroy(Top,ierr) 210*c4762a1bSJed Brown call VecDestroy(Bottom,ierr) 211*c4762a1bSJed Brown call VecDestroy(Left,ierr) 212*c4762a1bSJed Brown call VecDestroy(Right,ierr) 213*c4762a1bSJed Brown call MatDestroy(H,ierr) 214*c4762a1bSJed Brown call VecDestroy(localX,ierr) 215*c4762a1bSJed Brown call VecDestroy(localV,ierr) 216*c4762a1bSJed Brown call DMDestroy(dm,ierr) 217*c4762a1bSJed Brown 218*c4762a1bSJed Brown! Finalize TAO 219*c4762a1bSJed Brown 220*c4762a1bSJed Brown call PetscFinalize(ierr) 221*c4762a1bSJed Brown 222*c4762a1bSJed Brown end 223*c4762a1bSJed Brown 224*c4762a1bSJed Brown! --------------------------------------------------------------------- 225*c4762a1bSJed Brown! 226*c4762a1bSJed Brown! FormFunctionGradient - Evaluates function f(X). 227*c4762a1bSJed Brown! 228*c4762a1bSJed Brown! Input Parameters: 229*c4762a1bSJed Brown! tao - the Tao context 230*c4762a1bSJed Brown! X - the input vector 231*c4762a1bSJed Brown! dummy - optional user-defined context, as set by TaoSetFunction() 232*c4762a1bSJed Brown! (not used here) 233*c4762a1bSJed Brown! 234*c4762a1bSJed Brown! Output Parameters: 235*c4762a1bSJed Brown! fcn - the newly evaluated function 236*c4762a1bSJed Brown! G - the gradient vector 237*c4762a1bSJed Brown! info - error code 238*c4762a1bSJed Brown! 239*c4762a1bSJed Brown 240*c4762a1bSJed Brown 241*c4762a1bSJed Brown subroutine FormFunctionGradient(tao,X,fcn,G,dummy,ierr) 242*c4762a1bSJed Brown use mymodule 243*c4762a1bSJed Brown implicit none 244*c4762a1bSJed Brown 245*c4762a1bSJed Brown! Input/output variables 246*c4762a1bSJed Brown 247*c4762a1bSJed Brown Tao tao 248*c4762a1bSJed Brown PetscReal fcn 249*c4762a1bSJed Brown Vec X, G 250*c4762a1bSJed Brown PetscErrorCode ierr 251*c4762a1bSJed Brown PetscInt dummy 252*c4762a1bSJed Brown 253*c4762a1bSJed Brown PetscInt i,j,row 254*c4762a1bSJed Brown PetscInt xs, xm 255*c4762a1bSJed Brown PetscInt gxs, gxm 256*c4762a1bSJed Brown PetscInt ys, ym 257*c4762a1bSJed Brown PetscInt gys, gym 258*c4762a1bSJed Brown PetscReal ft,zero,hx,hy,hydhx,hxdhy 259*c4762a1bSJed Brown PetscReal area,rhx,rhy 260*c4762a1bSJed Brown PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3 261*c4762a1bSJed Brown PetscReal d4,d5,d6,d7,d8 262*c4762a1bSJed Brown PetscReal df1dxc,df2dxc,df3dxc,df4dxc 263*c4762a1bSJed Brown PetscReal df5dxc,df6dxc 264*c4762a1bSJed Brown PetscReal xc,xl,xr,xt,xb,xlt,xrb 265*c4762a1bSJed Brown 266*c4762a1bSJed Brown 267*c4762a1bSJed Brown! PETSc's VecGetArray acts differently in Fortran than it does in C. 268*c4762a1bSJed Brown! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 269*c4762a1bSJed Brown! will return an array of doubles referenced by x_array offset by x_index. 270*c4762a1bSJed Brown! i.e., to reference the kth element of X, use x_array(k + x_index). 271*c4762a1bSJed Brown! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 272*c4762a1bSJed Brown PetscReal g_v(0:1),x_v(0:1) 273*c4762a1bSJed Brown PetscReal top_v(0:1),left_v(0:1) 274*c4762a1bSJed Brown PetscReal right_v(0:1),bottom_v(0:1) 275*c4762a1bSJed Brown PetscOffset g_i,left_i,right_i 276*c4762a1bSJed Brown PetscOffset bottom_i,top_i,x_i 277*c4762a1bSJed Brown 278*c4762a1bSJed Brown ft = 0.0 279*c4762a1bSJed Brown zero = 0.0 280*c4762a1bSJed Brown hx = 1.0/real(mx + 1) 281*c4762a1bSJed Brown hy = 1.0/real(my + 1) 282*c4762a1bSJed Brown hydhx = hy/hx 283*c4762a1bSJed Brown hxdhy = hx/hy 284*c4762a1bSJed Brown area = 0.5 * hx * hy 285*c4762a1bSJed Brown rhx = real(mx) + 1.0 286*c4762a1bSJed Brown rhy = real(my) + 1.0 287*c4762a1bSJed Brown 288*c4762a1bSJed Brown 289*c4762a1bSJed Brown! Get local mesh boundaries 290*c4762a1bSJed Brown call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, & 291*c4762a1bSJed Brown & PETSC_NULL_INTEGER,ierr) 292*c4762a1bSJed Brown call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, & 293*c4762a1bSJed Brown & gxm,gym,PETSC_NULL_INTEGER,ierr) 294*c4762a1bSJed Brown 295*c4762a1bSJed Brown! Scatter ghost points to local vector 296*c4762a1bSJed Brown call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr) 297*c4762a1bSJed Brown call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr) 298*c4762a1bSJed Brown 299*c4762a1bSJed Brown! Initialize the vector to zero 300*c4762a1bSJed Brown call VecSet(localV,zero,ierr) 301*c4762a1bSJed Brown 302*c4762a1bSJed Brown! Get arrays to vector data (See note above about using VecGetArray in Fortran) 303*c4762a1bSJed Brown call VecGetArray(localX,x_v,x_i,ierr) 304*c4762a1bSJed Brown call VecGetArray(localV,g_v,g_i,ierr) 305*c4762a1bSJed Brown call VecGetArray(Top,top_v,top_i,ierr) 306*c4762a1bSJed Brown call VecGetArray(Bottom,bottom_v,bottom_i,ierr) 307*c4762a1bSJed Brown call VecGetArray(Left,left_v,left_i,ierr) 308*c4762a1bSJed Brown call VecGetArray(Right,right_v,right_i,ierr) 309*c4762a1bSJed Brown 310*c4762a1bSJed Brown! Compute function over the locally owned part of the mesh 311*c4762a1bSJed Brown do j = ys,ys+ym-1 312*c4762a1bSJed Brown do i = xs,xs+xm-1 313*c4762a1bSJed Brown row = (j-gys)*gxm + (i-gxs) 314*c4762a1bSJed Brown xc = x_v(row+x_i) 315*c4762a1bSJed Brown xt = xc 316*c4762a1bSJed Brown xb = xc 317*c4762a1bSJed Brown xr = xc 318*c4762a1bSJed Brown xl = xc 319*c4762a1bSJed Brown xrb = xc 320*c4762a1bSJed Brown xlt = xc 321*c4762a1bSJed Brown 322*c4762a1bSJed Brown if (i .eq. 0) then !left side 323*c4762a1bSJed Brown xl = left_v(j - ys + 1 + left_i) 324*c4762a1bSJed Brown xlt = left_v(j - ys + 2 + left_i) 325*c4762a1bSJed Brown else 326*c4762a1bSJed Brown xl = x_v(row - 1 + x_i) 327*c4762a1bSJed Brown endif 328*c4762a1bSJed Brown 329*c4762a1bSJed Brown if (j .eq. 0) then !bottom side 330*c4762a1bSJed Brown xb = bottom_v(i - xs + 1 + bottom_i) 331*c4762a1bSJed Brown xrb = bottom_v(i - xs + 2 + bottom_i) 332*c4762a1bSJed Brown else 333*c4762a1bSJed Brown xb = x_v(row - gxm + x_i) 334*c4762a1bSJed Brown endif 335*c4762a1bSJed Brown 336*c4762a1bSJed Brown if (i + 1 .eq. gxs + gxm) then !right side 337*c4762a1bSJed Brown xr = right_v(j - ys + 1 + right_i) 338*c4762a1bSJed Brown xrb = right_v(j - ys + right_i) 339*c4762a1bSJed Brown else 340*c4762a1bSJed Brown xr = x_v(row + 1 + x_i) 341*c4762a1bSJed Brown endif 342*c4762a1bSJed Brown 343*c4762a1bSJed Brown if (j + 1 .eq. gys + gym) then !top side 344*c4762a1bSJed Brown xt = top_v(i - xs + 1 + top_i) 345*c4762a1bSJed Brown xlt = top_v(i - xs + top_i) 346*c4762a1bSJed Brown else 347*c4762a1bSJed Brown xt = x_v(row + gxm + x_i) 348*c4762a1bSJed Brown endif 349*c4762a1bSJed Brown 350*c4762a1bSJed Brown if ((i .gt. gxs ) .and. (j + 1 .lt. gys + gym)) then 351*c4762a1bSJed Brown xlt = x_v(row - 1 + gxm + x_i) 352*c4762a1bSJed Brown endif 353*c4762a1bSJed Brown 354*c4762a1bSJed Brown if ((j .gt. gys) .and. (i + 1 .lt. gxs + gxm)) then 355*c4762a1bSJed Brown xrb = x_v(row + 1 - gxm + x_i) 356*c4762a1bSJed Brown endif 357*c4762a1bSJed Brown 358*c4762a1bSJed Brown d1 = xc-xl 359*c4762a1bSJed Brown d2 = xc-xr 360*c4762a1bSJed Brown d3 = xc-xt 361*c4762a1bSJed Brown d4 = xc-xb 362*c4762a1bSJed Brown d5 = xr-xrb 363*c4762a1bSJed Brown d6 = xrb-xb 364*c4762a1bSJed Brown d7 = xlt-xl 365*c4762a1bSJed Brown d8 = xt-xlt 366*c4762a1bSJed Brown 367*c4762a1bSJed Brown df1dxc = d1 * hydhx 368*c4762a1bSJed Brown df2dxc = d1 * hydhx + d4 * hxdhy 369*c4762a1bSJed Brown df3dxc = d3 * hxdhy 370*c4762a1bSJed Brown df4dxc = d2 * hydhx + d3 * hxdhy 371*c4762a1bSJed Brown df5dxc = d2 * hydhx 372*c4762a1bSJed Brown df6dxc = d4 * hxdhy 373*c4762a1bSJed Brown 374*c4762a1bSJed Brown d1 = d1 * rhx 375*c4762a1bSJed Brown d2 = d2 * rhx 376*c4762a1bSJed Brown d3 = d3 * rhy 377*c4762a1bSJed Brown d4 = d4 * rhy 378*c4762a1bSJed Brown d5 = d5 * rhy 379*c4762a1bSJed Brown d6 = d6 * rhx 380*c4762a1bSJed Brown d7 = d7 * rhy 381*c4762a1bSJed Brown d8 = d8 * rhx 382*c4762a1bSJed Brown 383*c4762a1bSJed Brown f1 = sqrt(1.0 + d1*d1 + d7*d7) 384*c4762a1bSJed Brown f2 = sqrt(1.0 + d1*d1 + d4*d4) 385*c4762a1bSJed Brown f3 = sqrt(1.0 + d3*d3 + d8*d8) 386*c4762a1bSJed Brown f4 = sqrt(1.0 + d3*d3 + d2*d2) 387*c4762a1bSJed Brown f5 = sqrt(1.0 + d2*d2 + d5*d5) 388*c4762a1bSJed Brown f6 = sqrt(1.0 + d4*d4 + d6*d6) 389*c4762a1bSJed Brown 390*c4762a1bSJed Brown ft = ft + f2 + f4 391*c4762a1bSJed Brown 392*c4762a1bSJed Brown df1dxc = df1dxc / f1 393*c4762a1bSJed Brown df2dxc = df2dxc / f2 394*c4762a1bSJed Brown df3dxc = df3dxc / f3 395*c4762a1bSJed Brown df4dxc = df4dxc / f4 396*c4762a1bSJed Brown df5dxc = df5dxc / f5 397*c4762a1bSJed Brown df6dxc = df6dxc / f6 398*c4762a1bSJed Brown 399*c4762a1bSJed Brown g_v(row + g_i) = 0.5 * (df1dxc + df2dxc + df3dxc + df4dxc + & 400*c4762a1bSJed Brown & df5dxc + df6dxc) 401*c4762a1bSJed Brown enddo 402*c4762a1bSJed Brown enddo 403*c4762a1bSJed Brown 404*c4762a1bSJed Brown! Compute triangular areas along the border of the domain. 405*c4762a1bSJed Brown if (xs .eq. 0) then ! left side 406*c4762a1bSJed Brown do j=ys,ys+ym-1 407*c4762a1bSJed Brown d3 = (left_v(j-ys+1+left_i) - left_v(j-ys+2+left_i)) & 408*c4762a1bSJed Brown & * rhy 409*c4762a1bSJed Brown d2 = (left_v(j-ys+1+left_i) - x_v((j-gys)*gxm + x_i)) & 410*c4762a1bSJed Brown & * rhx 411*c4762a1bSJed Brown ft = ft + sqrt(1.0 + d3*d3 + d2*d2) 412*c4762a1bSJed Brown enddo 413*c4762a1bSJed Brown endif 414*c4762a1bSJed Brown 415*c4762a1bSJed Brown 416*c4762a1bSJed Brown if (ys .eq. 0) then !bottom side 417*c4762a1bSJed Brown do i=xs,xs+xm-1 418*c4762a1bSJed Brown d2 = (bottom_v(i+1-xs+bottom_i)-bottom_v(i-xs+2+bottom_i)) & 419*c4762a1bSJed Brown & * rhx 420*c4762a1bSJed Brown d3 = (bottom_v(i-xs+1+bottom_i)-x_v(i-gxs+x_i))*rhy 421*c4762a1bSJed Brown ft = ft + sqrt(1.0 + d3*d3 + d2*d2) 422*c4762a1bSJed Brown enddo 423*c4762a1bSJed Brown endif 424*c4762a1bSJed Brown 425*c4762a1bSJed Brown 426*c4762a1bSJed Brown if (xs + xm .eq. mx) then ! right side 427*c4762a1bSJed Brown do j=ys,ys+ym-1 428*c4762a1bSJed Brown d1 = (x_v((j+1-gys)*gxm-1+x_i)-right_v(j-ys+1+right_i))*rhx 429*c4762a1bSJed Brown d4 = (right_v(j-ys+right_i) - right_v(j-ys+1+right_i))*rhy 430*c4762a1bSJed Brown ft = ft + sqrt(1.0 + d1*d1 + d4*d4) 431*c4762a1bSJed Brown enddo 432*c4762a1bSJed Brown endif 433*c4762a1bSJed Brown 434*c4762a1bSJed Brown 435*c4762a1bSJed Brown if (ys + ym .eq. my) then 436*c4762a1bSJed Brown do i=xs,xs+xm-1 437*c4762a1bSJed Brown d1 = (x_v((gym-1)*gxm+i-gxs+x_i) - top_v(i-xs+1+top_i))*rhy 438*c4762a1bSJed Brown d4 = (top_v(i-xs+1+top_i) - top_v(i-xs+top_i))*rhx 439*c4762a1bSJed Brown ft = ft + sqrt(1.0 + d1*d1 + d4*d4) 440*c4762a1bSJed Brown enddo 441*c4762a1bSJed Brown endif 442*c4762a1bSJed Brown 443*c4762a1bSJed Brown 444*c4762a1bSJed Brown if ((ys .eq. 0) .and. (xs .eq. 0)) then 445*c4762a1bSJed Brown d1 = (left_v(0 + left_i) - left_v(1 + left_i)) * rhy 446*c4762a1bSJed Brown d2 = (bottom_v(0+bottom_i)-bottom_v(1+bottom_i))*rhx 447*c4762a1bSJed Brown ft = ft + sqrt(1.0 + d1*d1 + d2*d2) 448*c4762a1bSJed Brown endif 449*c4762a1bSJed Brown 450*c4762a1bSJed Brown if ((ys + ym .eq. my) .and. (xs + xm .eq. mx)) then 451*c4762a1bSJed Brown d1 = (right_v(ym+1+right_i) - right_v(ym+right_i))*rhy 452*c4762a1bSJed Brown d2 = (top_v(xm+1+top_i) - top_v(xm + top_i))*rhx 453*c4762a1bSJed Brown ft = ft + sqrt(1.0 + d1*d1 + d2*d2) 454*c4762a1bSJed Brown endif 455*c4762a1bSJed Brown 456*c4762a1bSJed Brown ft = ft * area 457*c4762a1bSJed Brown call MPI_Allreduce(ft,fcn,1,MPIU_SCALAR, & 458*c4762a1bSJed Brown & MPIU_SUM,PETSC_COMM_WORLD,ierr) 459*c4762a1bSJed Brown 460*c4762a1bSJed Brown 461*c4762a1bSJed Brown! Restore vectors 462*c4762a1bSJed Brown call VecRestoreArray(localX,x_v,x_i,ierr) 463*c4762a1bSJed Brown call VecRestoreArray(localV,g_v,g_i,ierr) 464*c4762a1bSJed Brown call VecRestoreArray(Left,left_v,left_i,ierr) 465*c4762a1bSJed Brown call VecRestoreArray(Top,top_v,top_i,ierr) 466*c4762a1bSJed Brown call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr) 467*c4762a1bSJed Brown call VecRestoreArray(Right,right_v,right_i,ierr) 468*c4762a1bSJed Brown 469*c4762a1bSJed Brown! Scatter values to global vector 470*c4762a1bSJed Brown call DMLocalToGlobalBegin(dm,localV,INSERT_VALUES,G,ierr) 471*c4762a1bSJed Brown call DMLocalToGlobalEnd(dm,localV,INSERT_VALUES,G,ierr) 472*c4762a1bSJed Brown 473*c4762a1bSJed Brown call PetscLogFlops(70.0d0*xm*ym,ierr) 474*c4762a1bSJed Brown 475*c4762a1bSJed Brown return 476*c4762a1bSJed Brown end !FormFunctionGradient 477*c4762a1bSJed Brown 478*c4762a1bSJed Brown 479*c4762a1bSJed Brown 480*c4762a1bSJed Brown 481*c4762a1bSJed Brown 482*c4762a1bSJed Brown! ---------------------------------------------------------------------------- 483*c4762a1bSJed Brown! 484*c4762a1bSJed Brown! 485*c4762a1bSJed Brown! FormHessian - Evaluates Hessian matrix. 486*c4762a1bSJed Brown! 487*c4762a1bSJed Brown! Input Parameters: 488*c4762a1bSJed Brown!. tao - the Tao context 489*c4762a1bSJed Brown!. X - input vector 490*c4762a1bSJed Brown!. dummy - not used 491*c4762a1bSJed Brown! 492*c4762a1bSJed Brown! Output Parameters: 493*c4762a1bSJed Brown!. Hessian - Hessian matrix 494*c4762a1bSJed Brown!. Hpc - optionally different preconditioning matrix 495*c4762a1bSJed Brown!. flag - flag indicating matrix structure 496*c4762a1bSJed Brown! 497*c4762a1bSJed Brown! Notes: 498*c4762a1bSJed Brown! Due to mesh point reordering with DMs, we must always work 499*c4762a1bSJed Brown! with the local mesh points, and then transform them to the new 500*c4762a1bSJed Brown! global numbering with the local-to-global mapping. We cannot work 501*c4762a1bSJed Brown! directly with the global numbers for the original uniprocessor mesh! 502*c4762a1bSJed Brown! 503*c4762a1bSJed Brown! MatSetValuesLocal(), using the local ordering (including 504*c4762a1bSJed Brown! ghost points!) 505*c4762a1bSJed Brown! - Then set matrix entries using the local ordering 506*c4762a1bSJed Brown! by calling MatSetValuesLocal() 507*c4762a1bSJed Brown 508*c4762a1bSJed Brown subroutine FormHessian(tao, X, Hessian, Hpc, dummy, ierr) 509*c4762a1bSJed Brown use mymodule 510*c4762a1bSJed Brown implicit none 511*c4762a1bSJed Brown 512*c4762a1bSJed Brown Tao tao 513*c4762a1bSJed Brown Vec X 514*c4762a1bSJed Brown Mat Hessian,Hpc 515*c4762a1bSJed Brown PetscInt dummy 516*c4762a1bSJed Brown PetscErrorCode ierr 517*c4762a1bSJed Brown 518*c4762a1bSJed Brown PetscInt i,j,k,row 519*c4762a1bSJed Brown PetscInt xs,xm,gxs,gxm 520*c4762a1bSJed Brown PetscInt ys,ym,gys,gym 521*c4762a1bSJed Brown PetscInt col(0:6) 522*c4762a1bSJed Brown PetscReal hx,hy,hydhx,hxdhy,rhx,rhy 523*c4762a1bSJed Brown PetscReal f1,f2,f3,f4,f5,f6,d1,d2,d3 524*c4762a1bSJed Brown PetscReal d4,d5,d6,d7,d8 525*c4762a1bSJed Brown PetscReal xc,xl,xr,xt,xb,xlt,xrb 526*c4762a1bSJed Brown PetscReal hl,hr,ht,hb,hc,htl,hbr 527*c4762a1bSJed Brown 528*c4762a1bSJed Brown! PETSc's VecGetArray acts differently in Fortran than it does in C. 529*c4762a1bSJed Brown! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 530*c4762a1bSJed Brown! will return an array of doubles referenced by x_array offset by x_index. 531*c4762a1bSJed Brown! i.e., to reference the kth element of X, use x_array(k + x_index). 532*c4762a1bSJed Brown! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 533*c4762a1bSJed Brown PetscReal right_v(0:1),left_v(0:1) 534*c4762a1bSJed Brown PetscReal bottom_v(0:1),top_v(0:1) 535*c4762a1bSJed Brown PetscReal x_v(0:1) 536*c4762a1bSJed Brown PetscOffset x_i,right_i,left_i 537*c4762a1bSJed Brown PetscOffset bottom_i,top_i 538*c4762a1bSJed Brown PetscReal v(0:6) 539*c4762a1bSJed Brown PetscBool assembled 540*c4762a1bSJed Brown PetscInt i1 541*c4762a1bSJed Brown 542*c4762a1bSJed Brown i1=1 543*c4762a1bSJed Brown 544*c4762a1bSJed Brown! Set various matrix options 545*c4762a1bSJed Brown call MatSetOption(Hessian,MAT_IGNORE_OFF_PROC_ENTRIES, & 546*c4762a1bSJed Brown & PETSC_TRUE,ierr) 547*c4762a1bSJed Brown 548*c4762a1bSJed Brown! Get local mesh boundaries 549*c4762a1bSJed Brown call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, & 550*c4762a1bSJed Brown & PETSC_NULL_INTEGER,ierr) 551*c4762a1bSJed Brown call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, & 552*c4762a1bSJed Brown & PETSC_NULL_INTEGER,ierr) 553*c4762a1bSJed Brown 554*c4762a1bSJed Brown! Scatter ghost points to local vectors 555*c4762a1bSJed Brown call DMGlobalToLocalBegin(dm,X,INSERT_VALUES,localX,ierr) 556*c4762a1bSJed Brown call DMGlobalToLocalEnd(dm,X,INSERT_VALUES,localX,ierr) 557*c4762a1bSJed Brown 558*c4762a1bSJed Brown! Get pointers to vector data (see note on Fortran arrays above) 559*c4762a1bSJed Brown call VecGetArray(localX,x_v,x_i,ierr) 560*c4762a1bSJed Brown call VecGetArray(Top,top_v,top_i,ierr) 561*c4762a1bSJed Brown call VecGetArray(Bottom,bottom_v,bottom_i,ierr) 562*c4762a1bSJed Brown call VecGetArray(Left,left_v,left_i,ierr) 563*c4762a1bSJed Brown call VecGetArray(Right,right_v,right_i,ierr) 564*c4762a1bSJed Brown 565*c4762a1bSJed Brown! Initialize matrix entries to zero 566*c4762a1bSJed Brown call MatAssembled(Hessian,assembled,ierr) 567*c4762a1bSJed Brown if (assembled .eqv. PETSC_TRUE) call MatZeroEntries(Hessian,ierr) 568*c4762a1bSJed Brown 569*c4762a1bSJed Brown 570*c4762a1bSJed Brown rhx = real(mx) + 1.0 571*c4762a1bSJed Brown rhy = real(my) + 1.0 572*c4762a1bSJed Brown hx = 1.0/rhx 573*c4762a1bSJed Brown hy = 1.0/rhy 574*c4762a1bSJed Brown hydhx = hy/hx 575*c4762a1bSJed Brown hxdhy = hx/hy 576*c4762a1bSJed Brown! compute Hessian over the locally owned part of the mesh 577*c4762a1bSJed Brown 578*c4762a1bSJed Brown do i=xs,xs+xm-1 579*c4762a1bSJed Brown do j=ys,ys+ym-1 580*c4762a1bSJed Brown row = (j-gys)*gxm + (i-gxs) 581*c4762a1bSJed Brown 582*c4762a1bSJed Brown xc = x_v(row + x_i) 583*c4762a1bSJed Brown xt = xc 584*c4762a1bSJed Brown xb = xc 585*c4762a1bSJed Brown xr = xc 586*c4762a1bSJed Brown xl = xc 587*c4762a1bSJed Brown xrb = xc 588*c4762a1bSJed Brown xlt = xc 589*c4762a1bSJed Brown 590*c4762a1bSJed Brown if (i .eq. gxs) then ! Left side 591*c4762a1bSJed Brown xl = left_v(left_i + j - ys + 1) 592*c4762a1bSJed Brown xlt = left_v(left_i + j - ys + 2) 593*c4762a1bSJed Brown else 594*c4762a1bSJed Brown xl = x_v(x_i + row -1 ) 595*c4762a1bSJed Brown endif 596*c4762a1bSJed Brown 597*c4762a1bSJed Brown if (j .eq. gys) then ! bottom side 598*c4762a1bSJed Brown xb = bottom_v(bottom_i + i - xs + 1) 599*c4762a1bSJed Brown xrb = bottom_v(bottom_i + i - xs + 2) 600*c4762a1bSJed Brown else 601*c4762a1bSJed Brown xb = x_v(x_i + row - gxm) 602*c4762a1bSJed Brown endif 603*c4762a1bSJed Brown 604*c4762a1bSJed Brown if (i+1 .eq. gxs + gxm) then !right side 605*c4762a1bSJed Brown xr = right_v(right_i + j - ys + 1) 606*c4762a1bSJed Brown xrb = right_v(right_i + j - ys) 607*c4762a1bSJed Brown else 608*c4762a1bSJed Brown xr = x_v(x_i + row + 1) 609*c4762a1bSJed Brown endif 610*c4762a1bSJed Brown 611*c4762a1bSJed Brown if (j+1 .eq. gym+gys) then !top side 612*c4762a1bSJed Brown xt = top_v(top_i +i - xs + 1) 613*c4762a1bSJed Brown xlt = top_v(top_i + i - xs) 614*c4762a1bSJed Brown else 615*c4762a1bSJed Brown xt = x_v(x_i + row + gxm) 616*c4762a1bSJed Brown endif 617*c4762a1bSJed Brown 618*c4762a1bSJed Brown if ((i .gt. gxs) .and. (j+1 .lt. gys+gym)) then 619*c4762a1bSJed Brown xlt = x_v(x_i + row - 1 + gxm) 620*c4762a1bSJed Brown endif 621*c4762a1bSJed Brown 622*c4762a1bSJed Brown if ((i+1 .lt. gxs+gxm) .and. (j .gt. gys)) then 623*c4762a1bSJed Brown xrb = x_v(x_i + row + 1 - gxm) 624*c4762a1bSJed Brown endif 625*c4762a1bSJed Brown 626*c4762a1bSJed Brown d1 = (xc-xl)*rhx 627*c4762a1bSJed Brown d2 = (xc-xr)*rhx 628*c4762a1bSJed Brown d3 = (xc-xt)*rhy 629*c4762a1bSJed Brown d4 = (xc-xb)*rhy 630*c4762a1bSJed Brown d5 = (xrb-xr)*rhy 631*c4762a1bSJed Brown d6 = (xrb-xb)*rhx 632*c4762a1bSJed Brown d7 = (xlt-xl)*rhy 633*c4762a1bSJed Brown d8 = (xlt-xt)*rhx 634*c4762a1bSJed Brown 635*c4762a1bSJed Brown f1 = sqrt( 1.0 + d1*d1 + d7*d7) 636*c4762a1bSJed Brown f2 = sqrt( 1.0 + d1*d1 + d4*d4) 637*c4762a1bSJed Brown f3 = sqrt( 1.0 + d3*d3 + d8*d8) 638*c4762a1bSJed Brown f4 = sqrt( 1.0 + d3*d3 + d2*d2) 639*c4762a1bSJed Brown f5 = sqrt( 1.0 + d2*d2 + d5*d5) 640*c4762a1bSJed Brown f6 = sqrt( 1.0 + d4*d4 + d6*d6) 641*c4762a1bSJed Brown 642*c4762a1bSJed Brown 643*c4762a1bSJed Brown hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ & 644*c4762a1bSJed Brown & (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2) 645*c4762a1bSJed Brown 646*c4762a1bSJed Brown hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ & 647*c4762a1bSJed Brown & (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4) 648*c4762a1bSJed Brown 649*c4762a1bSJed Brown ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ & 650*c4762a1bSJed Brown & (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4) 651*c4762a1bSJed Brown 652*c4762a1bSJed Brown hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ & 653*c4762a1bSJed Brown & (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2) 654*c4762a1bSJed Brown 655*c4762a1bSJed Brown hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6) 656*c4762a1bSJed Brown htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3) 657*c4762a1bSJed Brown 658*c4762a1bSJed Brown hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + & 659*c4762a1bSJed Brown & hxdhy*(1.0+d8*d8)/(f3*f3*f3) + & 660*c4762a1bSJed Brown & hydhx*(1.0+d5*d5)/(f5*f5*f5) + & 661*c4762a1bSJed Brown & hxdhy*(1.0+d6*d6)/(f6*f6*f6) + & 662*c4762a1bSJed Brown & (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)- & 663*c4762a1bSJed Brown & 2*d1*d4)/(f2*f2*f2) + (hxdhy*(1.0+d2*d2)+ & 664*c4762a1bSJed Brown & hydhx*(1.0+d3*d3)-2*d2*d3)/(f4*f4*f4) 665*c4762a1bSJed Brown 666*c4762a1bSJed Brown hl = hl * 0.5 667*c4762a1bSJed Brown hr = hr * 0.5 668*c4762a1bSJed Brown ht = ht * 0.5 669*c4762a1bSJed Brown hb = hb * 0.5 670*c4762a1bSJed Brown hbr = hbr * 0.5 671*c4762a1bSJed Brown htl = htl * 0.5 672*c4762a1bSJed Brown hc = hc * 0.5 673*c4762a1bSJed Brown 674*c4762a1bSJed Brown k = 0 675*c4762a1bSJed Brown 676*c4762a1bSJed Brown if (j .gt. 0) then 677*c4762a1bSJed Brown v(k) = hb 678*c4762a1bSJed Brown col(k) = row - gxm 679*c4762a1bSJed Brown k=k+1 680*c4762a1bSJed Brown endif 681*c4762a1bSJed Brown 682*c4762a1bSJed Brown if ((j .gt. 0) .and. (i .lt. mx-1)) then 683*c4762a1bSJed Brown v(k) = hbr 684*c4762a1bSJed Brown col(k) = row-gxm+1 685*c4762a1bSJed Brown k=k+1 686*c4762a1bSJed Brown endif 687*c4762a1bSJed Brown 688*c4762a1bSJed Brown if (i .gt. 0) then 689*c4762a1bSJed Brown v(k) = hl 690*c4762a1bSJed Brown col(k) = row - 1 691*c4762a1bSJed Brown k = k+1 692*c4762a1bSJed Brown endif 693*c4762a1bSJed Brown 694*c4762a1bSJed Brown v(k) = hc 695*c4762a1bSJed Brown col(k) = row 696*c4762a1bSJed Brown k=k+1 697*c4762a1bSJed Brown 698*c4762a1bSJed Brown if (i .lt. mx-1) then 699*c4762a1bSJed Brown v(k) = hr 700*c4762a1bSJed Brown col(k) = row + 1 701*c4762a1bSJed Brown k=k+1 702*c4762a1bSJed Brown endif 703*c4762a1bSJed Brown 704*c4762a1bSJed Brown if ((i .gt. 0) .and. (j .lt. my-1)) then 705*c4762a1bSJed Brown v(k) = htl 706*c4762a1bSJed Brown col(k) = row + gxm - 1 707*c4762a1bSJed Brown k=k+1 708*c4762a1bSJed Brown endif 709*c4762a1bSJed Brown 710*c4762a1bSJed Brown if (j .lt. my-1) then 711*c4762a1bSJed Brown v(k) = ht 712*c4762a1bSJed Brown col(k) = row + gxm 713*c4762a1bSJed Brown k=k+1 714*c4762a1bSJed Brown endif 715*c4762a1bSJed Brown 716*c4762a1bSJed Brown! Set matrix values using local numbering, defined earlier in main routine 717*c4762a1bSJed Brown call MatSetValuesLocal(Hessian,i1,row,k,col,v,INSERT_VALUES, & 718*c4762a1bSJed Brown & ierr) 719*c4762a1bSJed Brown 720*c4762a1bSJed Brown 721*c4762a1bSJed Brown 722*c4762a1bSJed Brown enddo 723*c4762a1bSJed Brown enddo 724*c4762a1bSJed Brown 725*c4762a1bSJed Brown! restore vectors 726*c4762a1bSJed Brown call VecRestoreArray(localX,x_v,x_i,ierr) 727*c4762a1bSJed Brown call VecRestoreArray(Left,left_v,left_i,ierr) 728*c4762a1bSJed Brown call VecRestoreArray(Right,right_v,right_i,ierr) 729*c4762a1bSJed Brown call VecRestoreArray(Top,top_v,top_i,ierr) 730*c4762a1bSJed Brown call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr) 731*c4762a1bSJed Brown 732*c4762a1bSJed Brown 733*c4762a1bSJed Brown! Assemble the matrix 734*c4762a1bSJed Brown call MatAssemblyBegin(Hessian,MAT_FINAL_ASSEMBLY,ierr) 735*c4762a1bSJed Brown call MatAssemblyEnd(Hessian,MAT_FINAL_ASSEMBLY,ierr) 736*c4762a1bSJed Brown 737*c4762a1bSJed Brown call PetscLogFlops(199.0d0*xm*ym,ierr) 738*c4762a1bSJed Brown 739*c4762a1bSJed Brown return 740*c4762a1bSJed Brown end 741*c4762a1bSJed Brown 742*c4762a1bSJed Brown 743*c4762a1bSJed Brown 744*c4762a1bSJed Brown 745*c4762a1bSJed Brown 746*c4762a1bSJed Brown! Top,Left,Right,Bottom,bheight,mx,my,bmx,bmy,H, defined in plate2f.h 747*c4762a1bSJed Brown 748*c4762a1bSJed Brown! ---------------------------------------------------------------------------- 749*c4762a1bSJed Brown! 750*c4762a1bSJed Brown!/* 751*c4762a1bSJed Brown! MSA_BoundaryConditions - calculates the boundary conditions for the region 752*c4762a1bSJed Brown! 753*c4762a1bSJed Brown! 754*c4762a1bSJed Brown!*/ 755*c4762a1bSJed Brown 756*c4762a1bSJed Brown subroutine MSA_BoundaryConditions(ierr) 757*c4762a1bSJed Brown use mymodule 758*c4762a1bSJed Brown implicit none 759*c4762a1bSJed Brown 760*c4762a1bSJed Brown PetscErrorCode ierr 761*c4762a1bSJed Brown PetscInt i,j,k,limit,maxits 762*c4762a1bSJed Brown PetscInt xs, xm, gxs, gxm 763*c4762a1bSJed Brown PetscInt ys, ym, gys, gym 764*c4762a1bSJed Brown PetscInt bsize, lsize 765*c4762a1bSJed Brown PetscInt tsize, rsize 766*c4762a1bSJed Brown PetscReal one,two,three,tol 767*c4762a1bSJed Brown PetscReal scl,fnorm,det,xt 768*c4762a1bSJed Brown PetscReal yt,hx,hy,u1,u2,nf1,nf2 769*c4762a1bSJed Brown PetscReal njac11,njac12,njac21,njac22 770*c4762a1bSJed Brown PetscReal b, t, l, r 771*c4762a1bSJed Brown PetscReal boundary_v(0:1) 772*c4762a1bSJed Brown PetscOffset boundary_i 773*c4762a1bSJed Brown logical exitloop 774*c4762a1bSJed Brown PetscBool flg 775*c4762a1bSJed Brown 776*c4762a1bSJed Brown limit=0 777*c4762a1bSJed Brown maxits = 5 778*c4762a1bSJed Brown tol=1e-10 779*c4762a1bSJed Brown b=-0.5 780*c4762a1bSJed Brown t= 0.5 781*c4762a1bSJed Brown l=-0.5 782*c4762a1bSJed Brown r= 0.5 783*c4762a1bSJed Brown xt=0 784*c4762a1bSJed Brown yt=0 785*c4762a1bSJed Brown one=1.0 786*c4762a1bSJed Brown two=2.0 787*c4762a1bSJed Brown three=3.0 788*c4762a1bSJed Brown 789*c4762a1bSJed Brown 790*c4762a1bSJed Brown call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, & 791*c4762a1bSJed Brown & PETSC_NULL_INTEGER,ierr) 792*c4762a1bSJed Brown call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER, & 793*c4762a1bSJed Brown & gxm,gym,PETSC_NULL_INTEGER,ierr) 794*c4762a1bSJed Brown 795*c4762a1bSJed Brown bsize = xm + 2 796*c4762a1bSJed Brown lsize = ym + 2 797*c4762a1bSJed Brown rsize = ym + 2 798*c4762a1bSJed Brown tsize = xm + 2 799*c4762a1bSJed Brown 800*c4762a1bSJed Brown 801*c4762a1bSJed Brown call VecCreateMPI(PETSC_COMM_WORLD,bsize,PETSC_DECIDE,Bottom,ierr) 802*c4762a1bSJed Brown call VecCreateMPI(PETSC_COMM_WORLD,tsize,PETSC_DECIDE,Top,ierr) 803*c4762a1bSJed Brown call VecCreateMPI(PETSC_COMM_WORLD,lsize,PETSC_DECIDE,Left,ierr) 804*c4762a1bSJed Brown call VecCreateMPI(PETSC_COMM_WORLD,rsize,PETSC_DECIDE,Right,ierr) 805*c4762a1bSJed Brown 806*c4762a1bSJed Brown hx= (r-l)/(mx+1) 807*c4762a1bSJed Brown hy= (t-b)/(my+1) 808*c4762a1bSJed Brown 809*c4762a1bSJed Brown do j=0,3 810*c4762a1bSJed Brown 811*c4762a1bSJed Brown if (j.eq.0) then 812*c4762a1bSJed Brown yt=b 813*c4762a1bSJed Brown xt=l+hx*xs 814*c4762a1bSJed Brown limit=bsize 815*c4762a1bSJed Brown call VecGetArray(Bottom,boundary_v,boundary_i,ierr) 816*c4762a1bSJed Brown 817*c4762a1bSJed Brown 818*c4762a1bSJed Brown elseif (j.eq.1) then 819*c4762a1bSJed Brown yt=t 820*c4762a1bSJed Brown xt=l+hx*xs 821*c4762a1bSJed Brown limit=tsize 822*c4762a1bSJed Brown call VecGetArray(Top,boundary_v,boundary_i,ierr) 823*c4762a1bSJed Brown 824*c4762a1bSJed Brown elseif (j.eq.2) then 825*c4762a1bSJed Brown yt=b+hy*ys 826*c4762a1bSJed Brown xt=l 827*c4762a1bSJed Brown limit=lsize 828*c4762a1bSJed Brown call VecGetArray(Left,boundary_v,boundary_i,ierr) 829*c4762a1bSJed Brown 830*c4762a1bSJed Brown elseif (j.eq.3) then 831*c4762a1bSJed Brown yt=b+hy*ys 832*c4762a1bSJed Brown xt=r 833*c4762a1bSJed Brown limit=rsize 834*c4762a1bSJed Brown call VecGetArray(Right,boundary_v,boundary_i,ierr) 835*c4762a1bSJed Brown endif 836*c4762a1bSJed Brown 837*c4762a1bSJed Brown 838*c4762a1bSJed Brown do i=0,limit-1 839*c4762a1bSJed Brown 840*c4762a1bSJed Brown u1=xt 841*c4762a1bSJed Brown u2=-yt 842*c4762a1bSJed Brown k = 0 843*c4762a1bSJed Brown exitloop = .false. 844*c4762a1bSJed Brown do while (k .lt. maxits .and. (.not. exitloop) ) 845*c4762a1bSJed Brown 846*c4762a1bSJed Brown nf1=u1 + u1*u2*u2 - u1*u1*u1/three-xt 847*c4762a1bSJed Brown nf2=-u2 - u1*u1*u2 + u2*u2*u2/three-yt 848*c4762a1bSJed Brown fnorm=sqrt(nf1*nf1+nf2*nf2) 849*c4762a1bSJed Brown if (fnorm .gt. tol) then 850*c4762a1bSJed Brown njac11=one+u2*u2-u1*u1 851*c4762a1bSJed Brown njac12=two*u1*u2 852*c4762a1bSJed Brown njac21=-two*u1*u2 853*c4762a1bSJed Brown njac22=-one - u1*u1 + u2*u2 854*c4762a1bSJed Brown det = njac11*njac22-njac21*njac12 855*c4762a1bSJed Brown u1 = u1-(njac22*nf1-njac12*nf2)/det 856*c4762a1bSJed Brown u2 = u2-(njac11*nf2-njac21*nf1)/det 857*c4762a1bSJed Brown else 858*c4762a1bSJed Brown exitloop = .true. 859*c4762a1bSJed Brown endif 860*c4762a1bSJed Brown k=k+1 861*c4762a1bSJed Brown enddo 862*c4762a1bSJed Brown 863*c4762a1bSJed Brown boundary_v(i + boundary_i) = u1*u1-u2*u2 864*c4762a1bSJed Brown if ((j .eq. 0) .or. (j .eq. 1)) then 865*c4762a1bSJed Brown xt = xt + hx 866*c4762a1bSJed Brown else 867*c4762a1bSJed Brown yt = yt + hy 868*c4762a1bSJed Brown endif 869*c4762a1bSJed Brown 870*c4762a1bSJed Brown enddo 871*c4762a1bSJed Brown 872*c4762a1bSJed Brown 873*c4762a1bSJed Brown if (j.eq.0) then 874*c4762a1bSJed Brown call VecRestoreArray(Bottom,boundary_v,boundary_i,ierr) 875*c4762a1bSJed Brown elseif (j.eq.1) then 876*c4762a1bSJed Brown call VecRestoreArray(Top,boundary_v,boundary_i,ierr) 877*c4762a1bSJed Brown elseif (j.eq.2) then 878*c4762a1bSJed Brown call VecRestoreArray(Left,boundary_v,boundary_i,ierr) 879*c4762a1bSJed Brown elseif (j.eq.3) then 880*c4762a1bSJed Brown call VecRestoreArray(Right,boundary_v,boundary_i,ierr) 881*c4762a1bSJed Brown endif 882*c4762a1bSJed Brown 883*c4762a1bSJed Brown enddo 884*c4762a1bSJed Brown 885*c4762a1bSJed Brown 886*c4762a1bSJed Brown! Scale the boundary if desired 887*c4762a1bSJed Brown call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 888*c4762a1bSJed Brown & '-bottom',scl,flg,ierr) 889*c4762a1bSJed Brown if (flg .eqv. PETSC_TRUE) then 890*c4762a1bSJed Brown call VecScale(Bottom,scl,ierr) 891*c4762a1bSJed Brown endif 892*c4762a1bSJed Brown 893*c4762a1bSJed Brown call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 894*c4762a1bSJed Brown & '-top',scl,flg,ierr) 895*c4762a1bSJed Brown if (flg .eqv. PETSC_TRUE) then 896*c4762a1bSJed Brown call VecScale(Top,scl,ierr) 897*c4762a1bSJed Brown endif 898*c4762a1bSJed Brown 899*c4762a1bSJed Brown call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 900*c4762a1bSJed Brown & '-right',scl,flg,ierr) 901*c4762a1bSJed Brown if (flg .eqv. PETSC_TRUE) then 902*c4762a1bSJed Brown call VecScale(Right,scl,ierr) 903*c4762a1bSJed Brown endif 904*c4762a1bSJed Brown 905*c4762a1bSJed Brown call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 906*c4762a1bSJed Brown & '-left',scl,flg,ierr) 907*c4762a1bSJed Brown if (flg .eqv. PETSC_TRUE) then 908*c4762a1bSJed Brown call VecScale(Left,scl,ierr) 909*c4762a1bSJed Brown endif 910*c4762a1bSJed Brown 911*c4762a1bSJed Brown 912*c4762a1bSJed Brown return 913*c4762a1bSJed Brown end 914*c4762a1bSJed Brown 915*c4762a1bSJed Brown! ---------------------------------------------------------------------------- 916*c4762a1bSJed Brown! 917*c4762a1bSJed Brown!/* 918*c4762a1bSJed Brown! MSA_Plate - Calculates an obstacle for surface to stretch over 919*c4762a1bSJed Brown! 920*c4762a1bSJed Brown! Output Parameter: 921*c4762a1bSJed Brown!. xl - lower bound vector 922*c4762a1bSJed Brown!. xu - upper bound vector 923*c4762a1bSJed Brown! 924*c4762a1bSJed Brown!*/ 925*c4762a1bSJed Brown 926*c4762a1bSJed Brown subroutine MSA_Plate(tao,xl,xu,dummy,ierr) 927*c4762a1bSJed Brown use mymodule 928*c4762a1bSJed Brown implicit none 929*c4762a1bSJed Brown 930*c4762a1bSJed Brown Tao tao 931*c4762a1bSJed Brown Vec xl,xu 932*c4762a1bSJed Brown PetscErrorCode ierr 933*c4762a1bSJed Brown PetscInt i,j,row 934*c4762a1bSJed Brown PetscInt xs, xm, ys, ym 935*c4762a1bSJed Brown PetscReal lb,ub 936*c4762a1bSJed Brown PetscInt dummy 937*c4762a1bSJed Brown 938*c4762a1bSJed Brown! PETSc's VecGetArray acts differently in Fortran than it does in C. 939*c4762a1bSJed Brown! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (PetscOffset) x_index, ierr) 940*c4762a1bSJed Brown! will return an array of doubles referenced by x_array offset by x_index. 941*c4762a1bSJed Brown! i.e., to reference the kth element of X, use x_array(k + x_index). 942*c4762a1bSJed Brown! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 943*c4762a1bSJed Brown PetscReal xl_v(0:1) 944*c4762a1bSJed Brown PetscOffset xl_i 945*c4762a1bSJed Brown 946*c4762a1bSJed Brown lb = PETSC_NINFINITY 947*c4762a1bSJed Brown ub = PETSC_INFINITY 948*c4762a1bSJed Brown 949*c4762a1bSJed Brown if (bmy .lt. 0) bmy = 0 950*c4762a1bSJed Brown if (bmy .gt. my) bmy = my 951*c4762a1bSJed Brown if (bmx .lt. 0) bmx = 0 952*c4762a1bSJed Brown if (bmx .gt. mx) bmx = mx 953*c4762a1bSJed Brown 954*c4762a1bSJed Brown 955*c4762a1bSJed Brown call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, & 956*c4762a1bSJed Brown & PETSC_NULL_INTEGER,ierr) 957*c4762a1bSJed Brown 958*c4762a1bSJed Brown call VecSet(xl,lb,ierr) 959*c4762a1bSJed Brown call VecSet(xu,ub,ierr) 960*c4762a1bSJed Brown 961*c4762a1bSJed Brown call VecGetArray(xl,xl_v,xl_i,ierr) 962*c4762a1bSJed Brown 963*c4762a1bSJed Brown 964*c4762a1bSJed Brown do i=xs,xs+xm-1 965*c4762a1bSJed Brown 966*c4762a1bSJed Brown do j=ys,ys+ym-1 967*c4762a1bSJed Brown 968*c4762a1bSJed Brown row=(j-ys)*xm + (i-xs) 969*c4762a1bSJed Brown 970*c4762a1bSJed Brown if (i.ge.((mx-bmx)/2) .and. i.lt.(mx-(mx-bmx)/2) .and. & 971*c4762a1bSJed Brown & j.ge.((my-bmy)/2) .and. j.lt.(my-(my-bmy)/2)) then 972*c4762a1bSJed Brown xl_v(xl_i+row) = bheight 973*c4762a1bSJed Brown 974*c4762a1bSJed Brown endif 975*c4762a1bSJed Brown 976*c4762a1bSJed Brown enddo 977*c4762a1bSJed Brown enddo 978*c4762a1bSJed Brown 979*c4762a1bSJed Brown 980*c4762a1bSJed Brown call VecRestoreArray(xl,xl_v,xl_i,ierr) 981*c4762a1bSJed Brown 982*c4762a1bSJed Brown return 983*c4762a1bSJed Brown end 984*c4762a1bSJed Brown 985*c4762a1bSJed Brown 986*c4762a1bSJed Brown 987*c4762a1bSJed Brown 988*c4762a1bSJed Brown 989*c4762a1bSJed Brown! ---------------------------------------------------------------------------- 990*c4762a1bSJed Brown! 991*c4762a1bSJed Brown!/* 992*c4762a1bSJed Brown! MSA_InitialPoint - Calculates an obstacle for surface to stretch over 993*c4762a1bSJed Brown! 994*c4762a1bSJed Brown! Output Parameter: 995*c4762a1bSJed Brown!. X - vector for initial guess 996*c4762a1bSJed Brown! 997*c4762a1bSJed Brown!*/ 998*c4762a1bSJed Brown 999*c4762a1bSJed Brown subroutine MSA_InitialPoint(X, ierr) 1000*c4762a1bSJed Brown use mymodule 1001*c4762a1bSJed Brown implicit none 1002*c4762a1bSJed Brown 1003*c4762a1bSJed Brown Vec X 1004*c4762a1bSJed Brown PetscErrorCode ierr 1005*c4762a1bSJed Brown PetscInt start,i,j 1006*c4762a1bSJed Brown PetscInt row 1007*c4762a1bSJed Brown PetscInt xs,xm,gxs,gxm 1008*c4762a1bSJed Brown PetscInt ys,ym,gys,gym 1009*c4762a1bSJed Brown PetscReal zero, np5 1010*c4762a1bSJed Brown 1011*c4762a1bSJed Brown! PETSc's VecGetArray acts differently in Fortran than it does in C. 1012*c4762a1bSJed Brown! Calling VecGetArray((Vec) X, (PetscReal) x_array(0:1), (integer) x_index, ierr) 1013*c4762a1bSJed Brown! will return an array of doubles referenced by x_array offset by x_index. 1014*c4762a1bSJed Brown! i.e., to reference the kth element of X, use x_array(k + x_index). 1015*c4762a1bSJed Brown! Notice that by declaring the arrays with range (0:1), we are using the C 0-indexing practice. 1016*c4762a1bSJed Brown PetscReal left_v(0:1),right_v(0:1) 1017*c4762a1bSJed Brown PetscReal bottom_v(0:1),top_v(0:1) 1018*c4762a1bSJed Brown PetscReal x_v(0:1) 1019*c4762a1bSJed Brown PetscOffset left_i, right_i, top_i 1020*c4762a1bSJed Brown PetscOffset bottom_i,x_i 1021*c4762a1bSJed Brown PetscBool flg 1022*c4762a1bSJed Brown PetscRandom rctx 1023*c4762a1bSJed Brown 1024*c4762a1bSJed Brown zero = 0.0 1025*c4762a1bSJed Brown np5 = -0.5 1026*c4762a1bSJed Brown 1027*c4762a1bSJed Brown call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & 1028*c4762a1bSJed Brown & '-start', start,flg,ierr) 1029*c4762a1bSJed Brown 1030*c4762a1bSJed Brown if ((flg .eqv. PETSC_TRUE) .and. (start .eq. 0)) then ! the zero vector is reasonable 1031*c4762a1bSJed Brown call VecSet(X,zero,ierr) 1032*c4762a1bSJed Brown 1033*c4762a1bSJed Brown elseif ((flg .eqv. PETSC_TRUE) .and. (start .gt. 0)) then ! random start -0.5 < xi < 0.5 1034*c4762a1bSJed Brown call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr) 1035*c4762a1bSJed Brown do i=0,start-1 1036*c4762a1bSJed Brown call VecSetRandom(X,rctx,ierr) 1037*c4762a1bSJed Brown enddo 1038*c4762a1bSJed Brown 1039*c4762a1bSJed Brown call PetscRandomDestroy(rctx,ierr) 1040*c4762a1bSJed Brown call VecShift(X,np5,ierr) 1041*c4762a1bSJed Brown 1042*c4762a1bSJed Brown else ! average of boundary conditions 1043*c4762a1bSJed Brown 1044*c4762a1bSJed Brown! Get Local mesh boundaries 1045*c4762a1bSJed Brown call DMDAGetCorners(dm,xs,ys,PETSC_NULL_INTEGER,xm,ym, & 1046*c4762a1bSJed Brown & PETSC_NULL_INTEGER,ierr) 1047*c4762a1bSJed Brown call DMDAGetGhostCorners(dm,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, & 1048*c4762a1bSJed Brown & PETSC_NULL_INTEGER,ierr) 1049*c4762a1bSJed Brown 1050*c4762a1bSJed Brown 1051*c4762a1bSJed Brown 1052*c4762a1bSJed Brown! Get pointers to vector data 1053*c4762a1bSJed Brown call VecGetArray(Top,top_v,top_i,ierr) 1054*c4762a1bSJed Brown call VecGetArray(Bottom,bottom_v,bottom_i,ierr) 1055*c4762a1bSJed Brown call VecGetArray(Left,left_v,left_i,ierr) 1056*c4762a1bSJed Brown call VecGetArray(Right,right_v,right_i,ierr) 1057*c4762a1bSJed Brown 1058*c4762a1bSJed Brown call VecGetArray(localX,x_v,x_i,ierr) 1059*c4762a1bSJed Brown 1060*c4762a1bSJed Brown! Perform local computations 1061*c4762a1bSJed Brown do j=ys,ys+ym-1 1062*c4762a1bSJed Brown do i=xs,xs+xm-1 1063*c4762a1bSJed Brown row = (j-gys)*gxm + (i-gxs) 1064*c4762a1bSJed Brown x_v(x_i + row) = ((j+1)*bottom_v(bottom_i +i-xs+1)/my & 1065*c4762a1bSJed Brown & + (my-j+1)*top_v(top_i+i-xs+1)/(my+2) + & 1066*c4762a1bSJed Brown & (i+1)*left_v(left_i+j-ys+1)/mx + & 1067*c4762a1bSJed Brown & (mx-i+1)*right_v(right_i+j-ys+1)/(mx+2))*0.5 1068*c4762a1bSJed Brown enddo 1069*c4762a1bSJed Brown enddo 1070*c4762a1bSJed Brown 1071*c4762a1bSJed Brown! Restore vectors 1072*c4762a1bSJed Brown call VecRestoreArray(localX,x_v,x_i,ierr) 1073*c4762a1bSJed Brown 1074*c4762a1bSJed Brown call VecRestoreArray(Left,left_v,left_i,ierr) 1075*c4762a1bSJed Brown call VecRestoreArray(Top,top_v,top_i,ierr) 1076*c4762a1bSJed Brown call VecRestoreArray(Bottom,bottom_v,bottom_i,ierr) 1077*c4762a1bSJed Brown call VecRestoreArray(Right,right_v,right_i,ierr) 1078*c4762a1bSJed Brown 1079*c4762a1bSJed Brown call DMLocalToGlobalBegin(dm,localX,INSERT_VALUES,X,ierr) 1080*c4762a1bSJed Brown call DMLocalToGlobalEnd(dm,localX,INSERT_VALUES,X,ierr) 1081*c4762a1bSJed Brown 1082*c4762a1bSJed Brown endif 1083*c4762a1bSJed Brown 1084*c4762a1bSJed Brown return 1085*c4762a1bSJed Brown end 1086*c4762a1bSJed Brown 1087*c4762a1bSJed Brown! 1088*c4762a1bSJed Brown!/*TEST 1089*c4762a1bSJed Brown! 1090*c4762a1bSJed Brown! build: 1091*c4762a1bSJed Brown! requires: !complex 1092*c4762a1bSJed Brown! 1093*c4762a1bSJed Brown! test: 1094*c4762a1bSJed Brown! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4 1095*c4762a1bSJed Brown! filter: sort -b 1096*c4762a1bSJed Brown! filter_output: sort -b 1097*c4762a1bSJed Brown! requires: !single 1098*c4762a1bSJed Brown! 1099*c4762a1bSJed Brown! test: 1100*c4762a1bSJed Brown! suffix: 2 1101*c4762a1bSJed Brown! nsize: 2 1102*c4762a1bSJed Brown! args: -tao_smonitor -mx 8 -my 6 -bmx 3 -bmy 3 -bheight 0.2 -tao_type bqnls -tao_gatol 1.e-4 1103*c4762a1bSJed Brown! filter: sort -b 1104*c4762a1bSJed Brown! filter_output: sort -b 1105*c4762a1bSJed Brown! requires: !single 1106*c4762a1bSJed Brown! 1107*c4762a1bSJed Brown!TEST*/ 1108