main.tex 42.3 KB
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084
%------------------------------- Script de Duiks --------------------------%

% Include package and setting
\input{packageNSetting}


\begin{document}

%========================================= TITRE =====================================

\begin{lstlisting}
MUD2

c
c     file mud2.d
c
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c     *                                                               *
c     *                  copyright (c) 2008 by UCAR                   *
c     *                                                               *
c     *       University Corporation for Atmospheric Research         *
c     *                                                               *
c     *                      all rights reserved                      *
c     *                                                               *
c     *                     MUDPACK  version 5.0.1                    *
c     *                                                               *
c     *                 A Fortran Package of Multigrid                *
c     *                                                               *
c     *                Subroutines and Example Programs               *
c     *                                                               *
c     *      for Solving Elliptic Partial Differential Equations      *
c     *                                                               *
c     *                             by                                *
c     *                                                               *
c     *                         John Adams                            *
c     *                                                               *
c     *                             of                                *
c     *                                                               *
c     *         the National Center for Atmospheric Research          *
c     *                                                               *
c     *                Boulder, Colorado  (80307)  U.S.A.             *
c     *                                                               *
c     *                   which is sponsored by                       *
c     *                                                               *
c     *              the National Science Foundation                  *
c     *                                                               *
c     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c
c ... file mud2.d
c
c     contains documentation for:
c     subroutine mud2(iparm,fparm,work,coef,bndyc,rhs,phi,mgopt,ierror)
c     A sample fortran driver is file "tmud2.f".
c
c ... required MUDPACK files
c
c     mudcom.f
c
c ... purpose
c
c     subroutine mud2 automatically discretizes and attempts to compute
c     the second-order difference approximation to the two-dimensional
c     linear nonseparable elliptic partial differential equation on a
c     rectangle.  the approximation is generated on a uniform grid covering
c     the rectangle (see mesh description below).  boundary conditions
c     may be specified (dirchlet), periodic, or mixed derivative in any
c     combination.  the form of the pde solved is:
c
c
c          cxx(x,y)*pxx + cyy(x,y)*pyy + cx(x,y)*px + cy(x,y)*py +
c
c          ce(x,y)*p(x,y) = r(x,y).
c
c
c     pxx,pyy,px,py are second and first partial derivatives of the
c     unknown real solution function p(x,y) with respect to the
c     independent variables x,y.  cxx,cyy,cx,cy,ce are the known
c     real coefficients of the elliptic pde and r(x,y) is the known
c     real right hand side of the equation.  cxx and cyy should be
c     positive for all x,y in the solution region.  Nonseparability
c     means some of the coefficients depend on both x and y.  If
c     the PDE is separable subroutine mud2sp should be used instead
c     of mud2 or muh2.
c                                                                               
c
c ... mesh description . . .
c
c     the approximation is generated on a uniform nx by ny grid.  the grid
c     is superimposed on the rectangular solution region
c
c          [xa,xb] x [yc,yd].
c
c     let
c
c          dlx = (xb-xa)/(nx-1), dly = (yd-yc)/(ny-1)
c
c     be the uniform grid increments in the x,y directions. then
c
c          xi=xa+(i-1)*dlx,  yj=yc+(j-1)*dly
c
c     for i=1,...,nx and j=1,...,ny  denote the x,y uniform mesh points
c
c
c ... language                                                                  
c
c     fortran90/fortran77
c
c
c ... portability                                                               
c
c     mudpack5.0.1 software has been compiled and tested with Fortran77
c     and Fortran90 on a variety of platforms.
c                                                                               
c ... methods
c
c     details of the methods employeed by the solvers in mudpack are given
c     in [1,9].  [1,2,9] contain performance measurements on a variety of
c     elliptic pdes (see "references" in the file "readme").  in summary:
c
c *** discretization and solution (second-order solvers) (see [1])
c
c     the pde and boundary conditions are automatically discretized at all
c     grid levels using second-order finite difference formula.  diagonal
c     dominance at coarser grid levels is maintained in the presence of
c     nonzero first-order terms by adjusting the second-order coefficient
c     when necessary.  the resulting block tri-diagonal linear system is
c     approximated using multigrid iteration [10,11,13,15,16,18].  version
c     5.0.1 of mudpack uses only fully weighted residual restriction.  defaults
c     include cubic prolongation and w(2,1) cycles.  these can be overridden
c     with  selected multigrid options (see "mgopt").  error control based on
c     maximum relative differences is available. full multigrid cycling (fmg)
c     or cycling beginning or restarting at the finest grid level can be
c     selected. a menu of relaxation methods including gauss-seidel point,
c     line relaxation(s) (in any combination of directions) and planar
c     relaxation (for three-dimensional anisotropic problems) are provided.
c     all methods use ordering based on alternating points (red/black),
c     lines, or planes for cray vectorization and improved convergence
c     rates [14].
c
c *** higher order solution (fourth-order solvers) (see [9,19,21])
c
c     if the multigrid cycling results in a second-order estimate (i.e.,
c     discretization level error is reached) then this can be improved to a
c     fourth-order estimate using the technique of "deferred corrections."
c     the values in the solution array are used to generate a fourth-order
c     approximation to the truncation error.  second-order finite difference
c     formula are used to approximate third and fourth partial derivatives
c     of the solution function [3].  the truncation error estimate is
c     transferred down to all grid levels using weighted averaging where
c     it serves as a new right hand side.  the default multigrid options
c     are used to compute the fourth-order correction term which is added
c     to the original solution array.
c
c
c ... references (partial)
c
c
c     [1] J. Adams, "MUDPACK: Multigrid Fortran Software for the Efficient
c     Solution of Linear Elliptic Partial Differential Equations,"
c     Applied Math. and Comput. vol.34, Nov 1989, pp.113-146.
c
c     [2] J. Adams,"FMG Results with the Multigrid Software Package MUDPACK,"
c     proceedings of the fourth Copper Mountain Conference on Multigrid, SIAM,
c     1989, pp.1-12.
c     .
c     .
c     .
c     [7] J. Adams, R. Garcia, B. Gross, J. Hack, D. Haidvogel, and V. Pizzo,
c     "Applications of Multigrid Software in the Atmospheric Sciences,"
c      Mon. Wea. Rev.,vol. 120 # 7, July 1992, pp. 1447-1458.
c     .
c     .
c     .
c     [9] J. Adams, "Recent Enhancements in MUDPACK, a Multigrid Software
c     package for Elliptic Partial Differential Equations," Applied Math.
c     and Comp., 1991, vol. 43, May 1991, pp. 79-94.
c
c     [10]J. Adams, "MUDPACK-2: Multigrid Software for Approximating
c     Elliptic Partial Differential Equations on Uniform Grids with
c     any Resolution," Applied Math. and Comp., 1993, vol. 53, February
c     1993, pp. 235-249
c     .
c     .
c     .
c
c ... argument description
c                                                                               
c
c **********************************************************************
c *** input arguments *************************************************
c **********************************************************************
c
c
c ... iparm                                                                     
c
c          an integer vector of length 17 used to pass integer
c          arguments.  iparm is set internally and defined as
c          follows:
c
c
c ... intl=iparm(1)
c
c          an initialization argument.  intl=0  must be input
c          on an initial call. in this case input arguments will
c          be checked for errors and the elliptic partial differential
c          equation and boundary conditions will be discretized using
c          second order finite difference formula.
c
c ***      An approximation is NOT generated after an intl=0 call!
c          mud2 should be called with intl=1 to approximate the elliptic
c          PDE discretized by the intl=0 call.  intl=1 should also
c          be input if mud2 has been called earlier and only the
c          values in in rhs (see below) or gbdy (see bndyc below)
c          or phi (see below) have changed.  This will bypass
c          redundant pde discretization and argument checking
c          and save computational time.  Some examples of when
c          intl=1 calls should be used are:
c
c          (0) after a intl=0 argument checking and discretization call
c
c          (1) mud2 is being recalled for additional accuracy.  In
c              this case iguess=iparm(12)=1 should also be used.
c
c          (2) mud2 is being called every time step in a time dependent
c              problem (see discussion below) where the elliptic operator
c              does not depend on time.
c
c          (3) mud2 is being used to solve the same elliptic equation
c              for several different right hand sides (iguess=0 should
c              probably be used for each new righthand side).
c
c          intl = 0 must be input before calling with intl = 1 when any
c          of the following conditions hold:
c
c          (a) the initial call to mud2

c          (b) any of the integer arguments other than iguess=iparm(12)
c              or maxcy=iparm(13) or mgopt have changed since the previous
c              call.
c
c          (c) any of the floating point arguments other than tolmax=
c              fparm(5) have changed since the previous call
c
c          (d) any of the coefficients input by coef (see below) have
c              changed since the previous call
c
c          (e) any of the "alfa" coefficients input by bndyc (see below)
c              have changed since the previous call.
c
c          If any of (a) through (e) are true then the elliptic PDE
c          must be discretized or rediscretized.  If none of (a)
c          through (e) holds, calls can be made with intl=1.
c          Incorrect calls with intl=1 will produce erroneous results.
c  ***     The values set in the saved work space "work" (see below) with
c          an intl=0 call must be preserved with subsequent intl=1 calls.
c
c          MUDPACK software performance should be monitored for intl=1
c          calls.  The intl=0 discretization call performance depends
c          primarily on the efficiency or lack of efficiency of the
c          user provided subroutines for pde coefficients and
c          boundary conditions.
c
c ... nxa=iparm(2)
c
c          flags boundary conditions on the edge x=xa
c
c        = 0 if p(x,y) is periodic in x on [xa,xb]
c            (i.e., p(x+xb-xa,y) = p(x,y) for all x,y
c            (if nxa=0 then nxb=0 is required, see ierror = 2)
c
c        = 1 if p(xa,y) is specified (this must be input thru phi(1,j))
c
c        = 2 if there are mixed derivative boundary conditions at x=xa
c            (see bndyc)
c
c
c ... nxb=iparm(3)
c
c          flags boundary conditions on the edge x=xb
c
c        = 0 if p(x,y) is periodic in x on [xa,xb]
c            (i.e., p(x+xb-xa,y) = p(x,y) for all x,y)
c            (if nxb=0 then nxa=0 is required, see ierror = 2)
c
c        = 1 if p(xb,y) is specified (this must be input thru phi(nx,j))
c
c        = 2 if there are mixed derivative boundary conditions at x=xb
c            (see bndyc)
c
c
c ... nyc=iparm(4)
c
c          flags boundary conditions on the edge y=yc
c
c        = 0 if p(x,y) is periodic in y on [yc,yd]
c            (i.e., p(x,y+yd-yc) = p(x,y) for all x,y
c            (if nyc=0 then nyd=0 is required, see ierror = 2)
c
c        = 1 if p(x,yc) is specified (this must be input thru phi(i,1))
c
c        = 2 if there are mixed derivative boundary conditions at y=yc
c            (see bndyc)
c
c
c ... nyd=iparm(5)
c
c          flags boundary conditions on the edge y=yd
c
c        = 0 if p(x,y) is periodic in y on [yc,yd]
c            (i.e., p(x,y+yd-yc) = p(x,y) for all x,y
c             (if nyd=0 then nyc=0 is required, see ierror = 2)
c
c        = 1 if p(x,yd) is specified (this must be input thru phi(i,ny))
c
c        = 2 if there are mixed derivative boundary conditions at y=yd
c            (see bndyc)
c
c
c *** grid size arguments
c
c
c ... ixp = iparm(6)
c
c          an integer greater than one which is used in defining the number
c          of grid points in the x direction (see nx = iparm(10)).  "ixp+1"
c          is the number of points on the coarsest x grid visited during
c          multigrid cycling.  ixp should be chosen as small as possible.
c          recommended values are the small primes 2 or 3.
c          larger values can reduce multigrid convergence rates considerably,
c          especially if line relaxation in the x direction is not used.
c          if ixp > 2 then it should be 2 or a small odd value since a power
c          of 2 factor of ixp can be removed by increasing iex = iparm(8)
c          without changing nx = iparm(10).
c
c
c ... jyq = iparm(7)
c
c          an integer greater than one which is used in defining the number
c          of grid points in the y direction (see ny = iparm(11)).  "jyq+1"
c          is the number of points on the coarsest y grid visited during
c          multigrid cycling.  jyq should be chosen as small as possible.
c          recommended values are the small primes 2 or 3.
c          larger values can reduce multigrid convergence rates considerably,
c          especially if line relaxation in the y direction is not used.
c          if jyq > 2 then it should be 2 or a small odd value since a power
c          of 2 factor of jyq can be removed by increasing jey = iparm(9)
c          without changing ny = iparm(11).
c
c
c ... iex = iparm(8)
c
c          a positive integer exponent of 2 used in defining the number
c          of grid points in the x direction (see nx = iparm(10)).
c          iex .le. 50 is required.  for efficient multigrid cycling,
c          iex should be chosen as large as possible and ixp=iparm(8)
c          as small as possible within grid size constraints when
c          defining nx.
c
c
c ... jey = iparm(9)
c
c          a positive integer exponent of 2 used in defining the number
c          of grid points in the y direction (see ny = iparm(11)).
c          jey .le. 50 is required.  for efficient multigrid cycling,
c          jey should be chosen as large as possible and jyq=iparm(7)
c          as small as possible within grid size constraints when
c          defining ny.
c
c
c
c ... nx = iparm(10)
c
c          the number of equally spaced grid points in the interval [xa,xb]
c          (including the boundaries).  nx must have the form
c
c               nx = ixp*(2**(iex-1)) + 1
c
c          where ixp = iparm(6), iex = iparm(8).
c
c
c ... ny = iparm(11)
c
c          the number of equally spaced grid points in the interval [yc,yd]
c          (including the boundaries).  ny must have the form:
c
c               ny = jyq*(2**(jey-1)) + 1
c
c          where jyq = iparm(7), jey = iparm(9).
c
c
c *** example
c
c         suppose a solution is wanted on a 33 by 97 grid.  then
c         ixp=2, jyq=6 and iex=jey=5 could be used.  a better
c         choice would be ixp=2, jyq=3, and iex=5, jey=6.
c
c *** grid size flexibility considerations:
c
c     the hybrid multigrid/direct method code muh2 provides more grid size
c     flexibility than mud2 by removing the constraint that ixp and jyq are
c     2 or 3.  This is accomplished by using a direct method whenever the
c     coarsest (ixp+1) X (jyq+1) grid is encountered in multigrid cycling.
c     If nx = ixp+1 and ny = jyq+1 then muh2 becomes a full direct method.
c     muh2 is roughly equivalent to mud2 in efficiency as long as ixp and
c     jyq remain "small" (see muh2.d).  If the problem to be approximated
c     requires a grid neither mud2 por muh2 can exactly fit then another option
c     is to generate an approximation on a "close grid" using mud2 or muh2.
c     Then transfer the result to the required grid using cubic interpolation
c     via the package "regridpack"(contact John Adams about this software)
c
c *** note
c
c     let G be the nx by ny fine grid on which the approximation is
c     generated and let n = max0(iex,jey).  in mudpack, multigrid
c     cycling is implemented on the ascending chain of grids
c
c         G(1) < ... < G(k) < ... < G(n) = G.
c
c     each G(k) (k=1,...,n) has mx(k) by my(k) grid points
c     given by:
c
c         mx(k) = ixp*[2**(max0(iex+k-n,1)-1)] + 1
c
c         my(k) = jyq*[2**(max0(jey+k-n,1)-1)] + 1
c
c
c
c ... iguess=iparm(12)
c
c          = 0 if no initial guess to the pde is provided
c
c          = 1 if an initial guess to the pde is at the finest grid
c              level is provided in phi (see below)
c
c     comments on iguess = 0 or 1 . . .
c
c     even if iguess = 0, phi must be initialized at all grid points (this
c     is not checked).  phi can be set to 0.0 at non-dirchlet grid points
c     if nothing better is available.  the values set in phi when iguess = 0
c     are passed down and serve as an initial guess to the pde at the coarsest
c     grid level where cycling commences.  in this sense, values input in
c     phi always serve as an initial guess.  setting iguess = 0 forces full
c     multigrid cycling beginning at the coarsest and finishing at the finest
c     grid level.
c
c     if iguess = 1 then the values input in phi are an initial guess to the
c     pde at the finest grid level where cycling begins.  this option should
c     be used only if a "very good" initial guess is available (as, for
c     example, when restarting from a previous iguess=0 call).
c
c     time dependent problems . . .
c
c *** assume we are solving an elliptic pde every time step in a
c     marching problem of the form:
c
c          l(p(t)) = r(t)
c
c     where the differential operator "l" has no time dependence,
c     "p(t)" is the solution and "r(t)" is the right hand side at
c     current time "t". let "dt" be the increment between time steps.
c     then p(t) can be used as an initial guess to p(t+dt) with
c     intl = 1 when solving
c
c          l(p(t+dt)) = r(t+dt).
c
c     after the first two time steps, rather than continue, it would
c     be better to define the "correction" term:
c
c          e(t,dt) = p(t+dt) - p(t)
c
c     this clearly satisfies the equation
c
c          l(e(t,dt)) = r(t+dt) - r(t).
c
c     this should be solved with iguess = 0 and intl = 1. boundary
c     conditions for e(t,dt) are obtained from the boundary conditions
c     for p(t) by subtracting given values at t from given values at
c     t+dt. for example if
c
c          d(p(t))/dx = f(t), d(p(t+dt))/dx = f(t+dt)
c
c     at some x boundary then e(t,dt) satisfies the derivative
c     boundary condition
c
c          d(e(t,dt))/dx = f(t+dt) - f(t).
c
c     e(t,dt) can be preset to 0.0 (at nondirchlet points) or (if p(t-dt)
c     is saved) to p(t)-p(t-dt).  with iguess = 0, these values will serve
c     as an initial guess to e(t,dt) at the coarsest grid level.  this
c     approach has the advantage that a full sequence of multigrid cycles,
c     beginning at the coarsest grid level, is invoked every time step in
c     solving for e(t,dt).  a few digits of accuracy in e(t,dt), which is
c     ordinarily much smaller than p(t), will yield several more digits of
c     accuracy in the final approximation:
c
c          p(t+dt) = p(t) + e(t,dt).
c
c     using this approach to integrate in time will give more accuracy
c     then using p(t) as an initial guess to p(t+dt) for all time steps.
c     it does require additional storage.
c
c     if the differential operator "l" has time dependence (either thru
c     the coefficients in the pde or the coefficients in the derivative
c     boundary conditions) then use p(t) as an initial guess to p(t+dt)
c     when solving
c
c          l(t+dt)(p(t+dt)) = r(t+dt)
c
c     with intl = 0 for all time steps (the discretization must be repeated
c     for each new "t"). either iguess = 0 (p(t) will then be an initial
c     guess at the coarsest grid level where cycles will commence) or
c     iguess = 1 (p(t) will then be an initial guess at the finest grid
c     level where cycles will remain fixed) can be tried.
c
c
c ... maxcy = iparm(13)
c
c          the exact number of cycles executed between the finest (nx by
c          ny) and the coarsest ((ixp+1) by (jyq+1)) grid levels when
c          tolmax=fparm(5)=0.0 (no error control).  when tolmax > 0.0
c          is input (error control) then maxcy is a limit on the number
c          of cycles between the finest and coarsest grid levels.  in
c          any case, at most maxcy*(iprer+ipost) relaxation sweeps are
c          are performed at the finest grid level (see iprer=mgopt(2),
c          ipost=mgopt(3) below).  when multigrid iteration is working
c          "correctly" only a few are required for convergence.  large
c          values for maxcy should not be necessary.
c
c
c ... method = iparm(14) determines the method of relaxation
c              (gauss-seidel based on alternating points or lines)
c
c          = 0 for  point relaxation
c
c          = 1 for line relaxation in the x direction
c
c          = 2 for line relaxation in the y direction
c
c          = 3 for line relaxation in both the x and y direction
c
c
c *** choice of method. . .
c
c     let fx represent the quantity cxx(x,y)/dlx**2 over the solution region.
c
c     let fy represent the quantity cyy(x,y)/dly**2 over the solution region
c
c     if fx,fy are roughly the same size and do not vary too much over
c     the solution region choose method = 0.  if this fails try method=3.
c
c     if fx is much greater than fy choose method = 1.
c
c     if fy is much greater than fx choose method = 2
c
c     if neither fx or fy dominates over the solution region and they
c     both vary considerably choose method = 3.
c
c
c ... length = iparm(15)
c
c          the length of the work space provided in vector work (see below).
c          let isx = 0 if method = 0 or method = 2
c          let isx = 3 if method = 1 or method = 3 and nxa.ne.0
c          let isx = 5 if method = 1 or method = 3 and nxa.eq.0
c          let jsy = 0 if method = 0 or method = 1
c          let jsy = 3 if method = 2 or method = 3 and nyc.ne.0
c          let jsy = 5 if method = 2 or method = 3 and nyc.eq.0
c          then . . .
c
c               length = 4*[nx*ny*(10+isx+jsy)+8*(nx+ny+2)]/3
c
c          will suffice in most cases.  the exact minimal work space
c          length required for the current nx,ny and method is output
c          in iparm(16) (even if iparm(15) is too small).  this will be
c          less then the value given by the simplified formula above
c          in most cases.
c
c
c ... fparm                                                                     
c
c          a floating point vector of length 6 used to efficiently
c          pass floating point arguments.  fparm is set internally
c          in mud2 and defined as follows . . .
c
c
c ... xa=fparm(1), xb=fparm(2)
c
c          the range of the x independent variable. xa must                     
c          be less than xb                                                      
c
c
c ... yc=fparm(3), yd=fparm(4)
c
c          the range of the y independent variable.  yc must                    
c          be less than yd.                                                     
c
c
c ... tolmax = fparm(5)
c
c          when input positive, tolmax is a maximum relative error tolerance
c          used to terminate the relaxation iterations. assume phi1(i,j)
c          and phi2(i,j) are the last two computed approximations at all
c          grid points of the finest grid level. if we define
c
c              phdif = max(abs(phi2(i,j)-phi1(i,j))) for all i,j
c
c          and
c
c              phmax = max(abs(phi2(i,j))) for all i,j
c
c          then "convergence" is considered to have occurred if and only if
c
c              phdif/phmax < tolmax.
c
c
c          if tolmax=fparm(5)=0.0 is input then there is no error control
c          and maxcy cycles from the finest grid level are executed. maxcy
c          is a limit which cannot be exceeded even with error control.
c     ***  calls with tolmax=0.0, when appropriate because of known
c          convergence behavior, are more efficient than calls with tolmax
c          positive (i.e., if possible DO NOT use error control!).
c
c ... work                                                                      
c
c          a one dimensional real saved work space (see iparm(15) for
c          length) which must be preserved from the previous call when
c          calling with intl=iparm(1)=1.
c
c ... bndyc                                                                     
c
c          a subroutine with  arguments (kbdy,xory,alfa,gbdy) which
c          are used to input mixed boundary conditions to mud2. bndyc
c          must be declared "external" in the program calling mud2.
c          the boundaries are numbered one thru four and the mixed
c          derivative boundary conditions are described below (see the
c          sample driver code "tmud2.f" for an example of how bndyc is
c          can beset up).
c                                                                               
c          * * * * * * * * * * * *  y=yd
c          *       kbdy=4        *
c          *                     *
c          *                     *
c          *                     *
c          * kbdy=1       kbdy=2 *
c          *                     *
c          *                     *
c          *                     *
c          *       kbdy=3        *
c          * * * * * * * * * * * *  y=yc
c
c          x=xa                  x=xb
c
c
c          (1) the kbdy=1 boundary
c
c          this is the edge x=xa where nxa=iparm(2)=2 flags
c          a mixed boundary condition of the form
c
c             dp/dx + alfxa(y)*p(xa,y) = gbdxa(y)
c
c          in this case kbdy=1,xory=y will be input to bndyc and
c          alfa,gbdy corresponding to alfxa(y),gbdxa(y) must be returned.
c
c
c          (2) the kbdy=2 boundary
c
c          this is the edge x=xb where nxb=iparm(3)=2 flags
c          a mixed boundary condition of the form
c
c             dp/dx + alfxb(y)*p(xb,y) = gbdxb(y)
c
c          in this case kbdy=2,xory=y, will be input to bndyc and
c          alfa,gbdy corresponding to alfxb(y),gbdxb(y) must be returned.
c
c
c          (3) the kbdy=3 boundary
c
c          this is the edge y=yc where nyc=iparm(4)=2 flags
c          a mixed boundary condition of the form
c
c             dp/dy + alfyc(x)*p(x,yc) = gbdyc(x)
c
c          in this case kbdy=3,xory=x will be input to bndyc and
c          alfa,gbdy corresponding to alfyc(x),gbdyc(x) must be returned.
c
c
c          (4) the kbdy=4 boundary
c
c          this is the edge y=yd where nyd=iparm(5)=2 flags
c          a mixed boundary condition of the form
c
c             dp/dy + alfyd(x)*p(x,yd) = gbdyd(x)
c
c          in this case kbdy=4,xory=x will be input to bndyc and
c          alfa,gbdy corresponding to alfyd(x),gbdyd(x) must be returned.
c
c                                                                               
c ***      bndyc must provide the mixed boundary condition values
c          in correspondence with those flagged in iparm(2) thru
c          iparm(5).  if all boundaries are specified or periodic
c          mud2 will never call bndyc.  even then it must be entered
c          as a dummy subroutine. bndyc must be declared "external"
c          in the routine calling mud2.  the actual name chosen may
c          be different.
c
c
c ... coef
c
c         a subroutine with arguments (x,y,cxx,cyy,cx,cy,ce) which
c         provides the known real coefficients for the elliptic pde at
c         any grid point (x,y).  the name chosen in the calling routine
c         may be different where the coefficient routine must be declared
c         "external."
c
c ... rhs                                                                       
c
c          an array dimensioned nx by ny which contains the given
c          right hand side values on the uniform 2-d mesh.
c
c              rhs(i,j) = r(xi,yj) for i=1,...,nx and j=1,...,ny
c
c ... phi                                                                       
c
c          an array dimensioned nx by ny.  on input phi must contain
c          specified boundary values. for example, if nyd=iparm(5)=1
c          then phi(i,ny) must be set equal to p(xi,yd) for i=1,...nx
c          prior to calling mud2.  these values are preserved by mud2.
c          if an initial guess is provided (iguess=iparm(11)=1) it must
c          be input thru phi.
c
c
c ***      if no initial guess is given (iguess=0) then phi must still
c          be initialized at all grid points (this is not checked).  these
c          values will serve as an initial guess to the pde at the coarsest
c          grid level after a transfer from the fine solution grid.  set phi
c          equal to to 0.0 at all internal and non-specified boundaries
c          grid points if nothing better is available.
c
c
c ... mgopt
c
c           an integer vector of length 4 which allows the user to select
c           among various multigrid options.  if mgopt(1)=0 is input then
c           a default set of multigrid arguments (chosen for robustness)
c           will be internally selected and the remaining values in mgopt
c           will be ignored.  if mgopt(1) is nonzero then the arguments
c           in mgopt are set internally and defined as follows:  (see the
c           basic coarse grid correction algorithm below)
c
c
c     kcycle = mgopt(1)
c
c            = 0 if default multigrid options are to be used
c
c            = 1 if v cycling is to be used (the least expensive per cycle)
c
c            = 2 if w cycling is to be used (the default)
c
c            > 2 if more general k cycling is to be used
c             *** warning--values larger than 2 increase
c                 the execution time per cycle considerably and
c                 result in the nonfatal error ierror = -5
c                 which indicates inefficient multigrid cycling.
c
c    iprer = mgopt(2)
c
c           the number of "pre-relaxation" sweeps executed before the
c           residual is restricted and cycling is invoked at the next
c           coarser grid level (default value is 2 whenever mgopt(1)=0)
c
c    ipost = mgopt(3)
c
c           the number of "post relaxation" sweeps executed after cycling
c           has been invoked at the next coarser grid level and the residual
c           correction has been transferred back (default value is 1
c           whenever mgopt(1)=0).
c
c *** if iprer, ipost, or (especially) kcycle is greater than 2
c     than inefficient multigrid cycling has probably been chosen and
c     the nonfatal error (see below) ierror = -5 will be set.  note
c     this warning may be overridden by any other nonzero value
c     for ierror.
c
c   intpol = mgopt(4)
c
c          = 1 if multilinear prolongation (interpolation) is used to
c              transfer residual corrections and the pde approximation
c              from coarse to fine grids within full multigrid cycling.
c
c          = 3 if multicubic prolongation (interpolation) is used to
c              transfer residual corrections and the pde approximation
c              from coarse to fine grids within full multigrid cycling.
c              (this is the default value whenever mgopt(1)=0).
c
c *** the default values (2,2,1,3) in the vector mgopt were chosen for
c     robustness.  in some cases v(2,1) cycles with linear prolongation will
c     give good results with less computation (especially in two-dimensions).
c     this  was the default and only choice in an earlier version of mudpack
c     (see [1]) and can be set with the integer vector (1,2,1,1) in mgopt.
c
c *** the schedules for one full multigrid cycle (iguess=0) using v(2,1)
c     cycles and w(2,1) cycles are depicted for a four level grid below.
c     the number of relaxation sweeps when each grid is visited are indicated.
c     the "*" stands for prolongation of the full approximation and the "."
c     stands for transfer of residuals and residual corrections within the
c     coarse grid correction algorithm (see below).  all version 5.0.1
c     mudpack solvers use only fully weighted residual restriction
c
c     one fmg with v(2,1) cycles:
c
c
c     ------------------------------2-----------------1------     level 4
c                                  * .               .
c                                 *   .             .
c     ---------------2-----------1-----2-----------1---------     level 3
c                   * .         .       .         .
c                  *   .       .         .       .
c     ------2-----1-----2-----1-----------2-----1------------     level 2
c          * .   .       .   .             .   .
c         *   . .         . .               . .
c     ---3-----3-----------3-----------------3---------------     level 1
c
c
c     one fmg with w(2,1) cycles:
c
c     ------------------------2---------------------------1--     level 4
c                            * .                         .
c     ----------2-----------1---2-----------3-----------1----     level 3
c              * .         .     .         . .         .
c     ----2---1---2---3---1-------2---3---1---2---3---1------     level 2
c        * . .     . . . .         . . . .     . . . .
c     --6---6-------6---6-----------6---6-------6---6--------     level 1
c
c
c     the form of the "recursive" coarse grid correction cycling used
c     when kcycle.ge.0 is input is described below in pseudo-algorithmic
c     language.  it is implemented non-recursively in fortran in mudpack.
c
c     algorithm cgc(k,l(k),u(k),r(k),kcycle,iprer,ipost,iresw,intpol)
c
c *** approximately solve l(k)*u(k) = r(k) using multigrid iteration
c *** k is the current grid level
c *** l(k) is the discretized pde operator at level k
c *** u(k) is the initial guess at level k
c *** r(k) is the right hand side at level k
c *** i(k,k-1) is the restriction operator from level k to level k-1
c *** (the form of i(k,k-1) depends on iresw)
c *** i(k-1,k) is the prolongation operator from level k-1 to level k
c *** (the form of i(k-1,k) depends on intpol)
c
c     begin algorithm cgc
c
c ***   pre-relax at level k
c
c     . do (i=1,iprer)
c
c     . .  relax(l(k),u(k),r(k))
c
c     . end do
c
c     . if (k > 1) then
c
c ***     restrict the residual from level k to level k-1
c
c     . . r(k-1) = i(k,k-1)(r(k)-l(k)*u(k))
c
c     . . kount = 0
c
c     . . repeat
c
c ***     solve for the residual correction at level k-1 in u(k-1)
c ***     using algorithm cgc "kcycle" times (this is the recursion)
c
c     . . . kount = kount+1
c
c     . . . invoke cgc(k-1,l(k-1),u(k-1),r(k-1),kcycle,iprer,ipost,iresw)
c
c
c     . . until (kount.eq.kcycle)
c
c ***     transfer residual correction in u(k-1) to level k
c ***     with the prolongation operator and add to u(k)
c
c     . . u(k) = u(k) + i(k-1,k)(u(k-1))
c
c     . end if
c
c ***   post relax at level k
c
c     . do (i=1,ipost)
c
c     . . relax(l(k),u(k),r(k))
c
c     . end do
c
c     . return
c
c     end algorithm cgc
c
c
c **********************************************************************
c *** output arguments ************************************************
c **********************************************************************
c
c
c ... iparm(16)  *** set for intl=0 calls only
c
c          on output iparm(16) contains the actual work space length
c          required.  this will usually be less than that given by the
c          simplified formula for length=iparm(15) (see as input argument)
c
c
c ... iparm(17)  *** set for intl=1 calls only
c
c          on output iparm(17) contains the actual number of multigrid cycles
c          between the finest and coarsest grid levels used to obtain the
c          approximation when error control (tolmax > 0.0) is set.
c
c
c ... fparm(6)   *** set for intl=1 calls with fparm(5) > 0. only
c
c          on output fparm(6) contains the final computed maximum relative
c          difference between the last two iterates at the finest grid level.
c          fparm(6) is computed only if there is error control (tolmax > 0.0)
c          assume phi1(i,j,k) and phi2(i,j,k) are the last two computed
c          values for phi(i,j,k) at all points of the finest grid level.
c          if we define
c
c               phdif = max(abs(phi2(i,j)-phi1(i,j))) over all i,j
c
c          and
c
c               phmax = max(abs(phi2(i,j)) over all i,j
c
c          then
c
c               fparm(6) = phdif/phmax
c
c          is returned whenever phmax > 0.0. in the degenerate case
c          phmax = 0.0, fparm(6) = phdif is returned.
c
c
c ... work
c
c          on output work contains intermediate values that must not
c          be destroyed if mud2 is to be called again with intl=1
c
c
c ... phi   *** for intl=1 calls only
c
c          on output phi(i,j) contains the approximation to p(xi,yj)
c          for all mesh points i = 1,...,nx and j=1,...,ny.  the last
c          computed iterate in phi is returned even if convergence is
c          not obtained
c
c
c ... ierror
c
c          For intl=iparm(1)=0 initialization calls, ierror is an
c          error flag that indicates invalid input arguments when
c          returned positive and nonfatal warnings when returned
c          negative.  Argument checking and discretization
c          is bypassed for intl=1 calls which can only return
c          ierror = -1 or 0 or 1.
c
c
c     non-fatal warnings * * *
c
c
c     =-5 if kcycle=mgopt(1) is greater than 2. values larger than 2 results
c         in an algorithm which probably does far more computation than
c         necessary.  kcycle = 1 (v cycles)  or kcycle=2 (w cycles) should
c         suffice for most problems.  ierror = -5 is also set if either
c         iprer = mgopt(2) or ipost=mgopt(3) is greater than 2.  the
c         ierror=-5 flag is overridden by any other fatal or non-fatal
c         error.
c
c     =-4 if there are dominant nonzero first order terms in the pde which
c         make it "hyperbolic" at the finest grid level. numerically, this
c         happens if:
c
c              abs(cx)*dlx > 2.*abs(cxx)   (dlx = (xb-xa)/(nx-1))
c
c                         (or)
c
c              abs(cy)*dly > 2.*abs(cyy)   (dly = (yd-yc)/(ny-1))
c
c
c         at some fine grid point (xi,yj).  if an adjustment is not made the
c         condition can lead to a matrix coming from the discretization
c         which is not diagonally dominant and divergence is possible. since
c         the condition is "likely" at coarser grid levels for pde's with
c         nonzero first order terms, the adjustments (actually first order
c         approximations to the pde)
c
c
c             cxx = amax1(cxx,0.5*abs(cx)*dx)
c
c                          (and)
c
c             cyy = amax1(cyy,0.5*abs(cy)*dy)
c
c
c         (here dx,dy are the x,y mesh sizes of the subgrid)
c
c         are made to preserve convergence of multigrid iteration. if made
c         at the finest grid level, it can lead to convergence to an
c         erroneous solution (flagged by ierror = -4).  a possible remedy
c         is to increase resolution. the ierror = -4 flag overrides the
c         nonfatal ierror = -5 flag.
c
c
c     =-3  if the continuous elliptic pde is singular.  this means the
c          boundary conditions are periodic or pure derivative at all
c          boundaries and ce(x,y) = 0.0 for all x,y.  a solution is still
c          attempted but convergence may not occur due to ill-conditioning
c          of the linear system coming from the discretization.  the
c          ierror = -3 flag overrides the ierror=-4 and ierror=-5 nonfatal
c          flags.
c
c
c     =-2  if the pde is not elliptic (i.e., cxx*cyy.le.0.0 for some (xi,yj))
c          in this case a solution is still attempted although convergence
c          may not occur due to ill-conditioning of the linear system.
c          the ierror = -2 flag overrides the ierror=-5,-4,-3 nonfatal
c          flags.
c
c
c     =-1  if convergence to the tolerance specified in tolmax=fparm(5)>0.
c          is not obtained in maxcy=iparm(13) multigrid cycles between the
c          coarsest (ixp+1,jyq+1) and finest (nx,ny) grid levels.
c          in this case the last computed iterate is still returned.
c          the ierror = -1 flag overrides all other nonfatal flags
c
c
c     no errors * * *
c
c     = 0
c
c     fatal argument errors * * *
c
c     = 1 if intl=iparm(1) is not 0 on initial call or not 0 or 1
c         on subsequent calls
c
c     = 2 if any of the boundary condition flags nxa,nxb,nyc,nyd
c         in iparm(2),iparm(3),iparm(4),iparm(5) are not 0,1 or 2
c         or if nxa,nxb or nyc,nyd are not pairwise zero.
c
c     = 3 if mino(ixp,jyq) < 2 (ixp = iparm(6), jyq = iparm(7))
c
c     = 4 if min0(iex,jey) < 1 (iex = iparm(8), jey = iparm(9)) or
c         if max0(iex,jey) > 50
c
c     = 5 if nx.ne.ixp*2**(iex-1)+1 or ny.ne.jyq*2**(jey-1)+1
c         (nx = iparm(10), ny = iparm(11))
c
c     = 6 if iguess = iparm(12) is not equal to 0 or 1
c
c     = 7 if maxcy = iparm(13) < 1
c
c     = 8 if method = iparm(14) is not 0,1,2, or 3
c
c     = 9 if length = iparm(15) is too small (see iparm(16) on output
c         for minimum required work space length)
c
c     =10 if xa >= xb or yc >= yd
c         (xa=fparm(1),xb=fparm(2),yc=fparm(3),yd=fparm(4))
c
c     =11 if tolmax = fparm(5) < 0.0
c
c     errors in setting multigrid options * * * (see also ierror=-5)
c
c     =12 if kcycle = mgopt(1) < 0 or
c         if iprer = mgopt(2) < 1 or
c         if ipost = mgopt(3) < 1 or
c         if intpol = mgopt(4) is not 1 or 3
c
c *********************************************************
\end{lstlisting}


%==================================== Bibliographie ==================================

\end{document}