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