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.gt.0) then
     70:           JJ = II - n
     71:           PetscCallA(MatSetValues(A,i1,[II],i1,[JJ],[v],ADD_VALUES,ierr))
     72:         endif
     73:         if (i.lt.m-1) then
     74:           JJ = II + n
     75:           PetscCallA(MatSetValues(A,i1,[II],i1,[JJ],[v],ADD_VALUES,ierr))
     76:         endif
     77:         if (j.gt.0) then
     78:           JJ = II - 1
     79:           PetscCallA(MatSetValues(A,i1,[II],i1,[JJ],[v],ADD_VALUES,ierr))
     80:         endif
     81:         if (j.lt.n-1) then
     82:           JJ = II + 1
     83:           PetscCallA(MatSetValues(A,i1,[II],i1,[JJ],[v],ADD_VALUES,ierr))
     84:         endif
     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 .eq. 0) then
    150:          write(6,98) cntl
    151:          write(6,99) rinfo12,rinfo13,infog34
    152:       endif
    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 .eq. 0) then
    170:         if (norm .gt. 1.e-12) then
    171:            write(6,100) norm,its
    172:         else
    173:            write(6,110) its
    174:         endif
    175:       endif
    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*/