Actual source code: ex52f.F90

      1: !
      2: !   Modified from ex15f.F for testing MUMPS
      3: !   Solves a linear system in parallel with KSP.
      4: !  -------------------------------------------------------------------------

      6: program main
      7: #include <petsc/finclude/petscksp.h>
      8:   use petscksp
      9:   implicit none

     11: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     12: !                   Variable declarations
     13: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     14:   Vec x, b, u
     15:   Mat A
     16:   KSP ksp
     17:   PetscScalar v, one, neg_one
     18:   PetscReal norm, tol
     19:   PetscErrorCode ierr
     20:   PetscInt i, j, II, JJ, Istart
     21:   PetscInt Iend, m, n, i1, its, five
     22:   PetscMPIInt rank
     23:   PetscBool flg
     24: #if defined(PETSC_HAVE_MUMPS)
     25:   PC pc
     26:   Mat F
     27:   PetscInt ival, icntl, infog34
     28:   PetscReal cntl, rinfo12, rinfo13, val
     29: #endif

     31: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     32: !                 Beginning of program
     33: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     34:   PetscCallA(PetscInitialize(ierr))
     35:   one = 1.0
     36:   neg_one = -1.0
     37:   i1 = 1
     38:   m = 8
     39:   n = 7
     40:   five = 5
     41:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-m', m, flg, ierr))
     42:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
     43:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))

     45: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     46: !      Compute the matrix and right-hand-side vector that define
     47: !      the linear system, Ax = b.
     48: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     49:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
     50:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*n, m*n, ierr))
     51:   PetscCallA(MatSetType(A, MATAIJ, ierr))
     52:   PetscCallA(MatSetFromOptions(A, ierr))
     53:   PetscCallA(MatMPIAIJSetPreallocation(A, five, PETSC_NULL_INTEGER_ARRAY, five, PETSC_NULL_INTEGER_ARRAY, ierr))
     54:   PetscCallA(MatSeqAIJSetPreallocation(A, five, PETSC_NULL_INTEGER_ARRAY, ierr))

     56:   PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))

     58: !  Set matrix elements for the 2-D, five-point stencil in parallel.
     59: !   - Each processor needs to insert only elements that it owns
     60: !     locally (but any non-local elements will be sent to the
     61: !     appropriate processor during matrix assembly).
     62: !   - Always specify global row and columns of matrix entries.
     63: !   - Note that MatSetValues() uses 0-based row and column numbers
     64: !     in Fortran as well as in C.
     65:   do 10, II = Istart, Iend - 1
     66:     v = -1.0
     67:     i = II/n
     68:     j = II - i*n
     69:     if (i > 0) then
     70:       JJ = II - n
     71:       PetscCallA(MatSetValues(A, i1, [II], i1, [JJ], [v], ADD_VALUES, ierr))
     72:     end if
     73:     if (i < m - 1) then
     74:       JJ = II + n
     75:       PetscCallA(MatSetValues(A, i1, [II], i1, [JJ], [v], ADD_VALUES, ierr))
     76:     end if
     77:     if (j > 0) then
     78:       JJ = II - 1
     79:       PetscCallA(MatSetValues(A, i1, [II], i1, [JJ], [v], ADD_VALUES, ierr))
     80:     end if
     81:     if (j < n - 1) then
     82:       JJ = II + 1
     83:       PetscCallA(MatSetValues(A, i1, [II], i1, [JJ], [v], ADD_VALUES, ierr))
     84:     end if
     85:     v = 4.0
     86:     PetscCallA(MatSetValues(A, i1, [II], i1, [II], [v], ADD_VALUES, ierr))
     87: 10  continue

     89: !  Assemble matrix, using the 2-step process:
     90:     PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
     91:     PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

     93: !  Create parallel vectors.
     94:     PetscCallA(VecCreateFromOptions(PETSC_COMM_WORLD, PETSC_NULL_CHARACTER, i1, PETSC_DECIDE, m*n, u, ierr))
     95:     PetscCallA(VecDuplicate(u, b, ierr))
     96:     PetscCallA(VecDuplicate(b, x, ierr))

     98: !  Set exact solution; then compute right-hand-side vector.
     99:     PetscCallA(VecSet(u, one, ierr))
    100:     PetscCallA(MatMult(A, u, b, ierr))

    102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    103: !         Create the linear solver and set various options
    104: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    105:     PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
    106:     PetscCallA(KSPSetOperators(ksp, A, A, ierr))
    107:     tol = 1.e-7
    108:     PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr))

    110: !  Test MUMPS
    111: #if defined(PETSC_HAVE_MUMPS)
    112:     PetscCallA(KSPSetType(ksp, KSPPREONLY, ierr))
    113:     PetscCallA(KSPGetPC(ksp, pc, ierr))
    114:     PetscCallA(PCSetType(pc, PCLU, ierr))
    115:     PetscCallA(PCFactorSetMatSolverType(pc, MATSOLVERMUMPS, ierr))
    116:     PetscCallA(PCFactorSetUpMatSolverType(pc, ierr))
    117:     PetscCallA(PCFactorGetMatrix(pc, F, ierr))

    119: !     sequential ordering
    120:     icntl = 7
    121:     ival = 2
    122:     PetscCallA(MatMumpsSetIcntl(F, icntl, ival, ierr))

    124: !     threshold for row pivot detection
    125:     icntl = 24
    126:     ival = 1
    127:     PetscCallA(MatMumpsSetIcntl(F, icntl, ival, ierr))
    128:     icntl = 3
    129:     val = 1.e-6
    130:     PetscCallA(MatMumpsSetCntl(F, icntl, val, ierr))

    132: !     compute determinant of A
    133:     icntl = 33
    134:     ival = 1
    135:     PetscCallA(MatMumpsSetIcntl(F, icntl, ival, ierr))
    136: #endif

    138:     PetscCallA(KSPSetFromOptions(ksp, ierr))
    139:     PetscCallA(KSPSetUp(ksp, ierr))
    140: #if defined(PETSC_HAVE_MUMPS)
    141:     icntl = 3
    142:     PetscCallA(MatMumpsGetCntl(F, icntl, cntl, ierr))
    143:     icntl = 34
    144:     PetscCallA(MatMumpsGetInfog(F, icntl, infog34, ierr))
    145:     icntl = 12
    146:     PetscCallA(MatMumpsGetRinfog(F, icntl, rinfo12, ierr))
    147:     icntl = 13
    148:     PetscCallA(MatMumpsGetRinfog(F, icntl, rinfo13, ierr))
    149:     if (rank == 0) then
    150:       write (6, 98) cntl
    151:       write (6, 99) rinfo12, rinfo13, infog34
    152:     end if
    153: 98  format('Mumps row pivot threshold = ', 1pe11.2)
    154: 99  format('Mumps determinant=(', 1pe11.2, 1pe11.2, ')*2^', i3)
    155: #endif

    157: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    158: !                      Solve the linear system
    159: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    160:     PetscCallA(KSPSolve(ksp, b, x, ierr))

    162: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    163: !                     Check solution and clean up
    164: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    165:     PetscCallA(VecAXPY(x, neg_one, u, ierr))
    166:     PetscCallA(VecNorm(x, NORM_2, norm, ierr))
    167:     PetscCallA(KSPGetIterationNumber(ksp, its, ierr))

    169:     if (rank == 0) then
    170:       if (norm > 1.e-12) then
    171:         write (6, 100) norm, its
    172:       else
    173:         write (6, 110) its
    174:       end if
    175:     end if
    176: 100 format('Norm of error ', 1pe11.4, ' iterations ', i5)
    177: 110 format('Norm of error < 1.e-12,iterations ', i5)

    179: !  Free work space.  All PETSc objects should be destroyed when they
    180: !  are no longer needed.

    182:     PetscCallA(KSPDestroy(ksp, ierr))
    183:     PetscCallA(VecDestroy(u, ierr))
    184:     PetscCallA(VecDestroy(x, ierr))
    185:     PetscCallA(VecDestroy(b, ierr))
    186:     PetscCallA(MatDestroy(A, ierr))

    188: !  Always call PetscFinalize() before exiting a program.
    189:     PetscCallA(PetscFinalize(ierr))
    190:   end

    192: !
    193: !/*TEST
    194: !
    195: !   test:
    196: !      suffix: mumps
    197: !      nsize: 3
    198: !      requires: mumps double
    199: !      output_file: output/ex52f_1.out
    200: !
    201: !TEST*/