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