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