G.3.7 Solve Procedure

The main documentation of the Solve Procedure contains additional explanation of this code listing.

  define([RESIDUAL_DEBUG_LEVEL],[6])

  subroutine Solve (Solver, A_ELLM, X_MV, B_MV)

    ! Use association information.

    ifdef([USE_LAMG],[
      use LAMG_Module
    ])
    ifdef([USE_MPI],[
      include(mpif.h)
    ])

    ! Input variables.

    type(Solver_type), intent(in) :: Solver             ! Solver to be set.
    type(ELL_Matrix_type), intent(inout) :: A_ELLM      ! Matrix in ELL format.
    type(Mathematic_Vector_type), intent(inout) :: B_MV ! RHS vector.

    ! Input/Output variables.


    ! Output variables.

    type(Mathematic_Vector_type), intent(inout) :: X_MV ! Solution vector.

    ! Internal variables.

    type(Status_type) :: status

    ! LAMG variables:
    ifdef([USE_LAMG],[
      type(LAMG_comm) :: LAMG_Communicator  ! LAMG communicator.
      type(LAMG_LS_Options) :: LAMG_Options ! LAMG options.
      type(integer) :: LAMG_Status          ! LAMG status.
      type(LAMG_Matrix_dcsr_r) :: A_LAMG    ! Matrix in LAMG format.
      type(real,1) :: X_LAMG                ! Solution vector in LAMG format.
      type(real,1) :: B_LAMG                ! RHS vector in LAMG format.
    ])
    ifelse(m4_eval(DEBUG_LEVEL >= RESIDUAL_DEBUG_LEVEL), 1, [
      type(Mathematic_Vector_type) :: Residual_MV  ! Residual vector.
    ])

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    ! Verify requirements.

    !VERIFY(Valid_State(input),0)  ! Quick, important test.
    !VERIFY(Valid_State(input),9)  ! Slow, unimportant test.

    ! LAMG solve.

    if (Solver%Package == 'LAMG') then

      ifdef([USE_LAMG],[

        ! Initialize LAMG communication.
  
        ifdef([USE_MPI],[
          call LAMG_Comm_Init (MPI_COMM_WORLD, LAMG_Communicator, LAMG_Status)
        ],[
          call LAMG_Comm_Init (             0, LAMG_Communicator, LAMG_Status)
        ])
  
        ! Initialize LAMG options.
  
        call LAMG_Options_Init (LAMG_Communicator, LAMG_Options, LAMG_Status)
  
        ! Extract RHS and Solution (Initial Guess) vector.
  
        call Initialize (B_LAMG, NRows_PE(A_ELLM), status)
        B_LAMG = B_MV
        call Initialize (X_LAMG, NRows_PE(A_ELLM), status)
        X_LAMG = X_MV
  
        ! Convert matrix to LAMG format.
  
        call Convert (A_LAMG, A_ELLM, LAMG_Communicator, LAMG_Options, status)
  
        ! Set LAMG options.
  
        ! Maximum allowed number of iterations.
        call LAMG_Set_Option_i (LAMG_itsmax, Solver%Maximum_Iterations, &
                                LAMG_Communicator, LAMG_Options, LAMG_Status)
        ! Convergence criterion.
        call LAMG_Set_Option_r (LAMG_tol, Solver%Epsilon, LAMG_Communicator, &
                                LAMG_Options, LAMG_Status)
        ! Package output: 0=silent, 4=verbose.
        call LAMG_Set_Option_i (LAMG_levout, Solver%LAMG_levout, &
                                LAMG_Communicator, LAMG_Options, LAMG_Status)
        ! LAMG does not currently allow any stopping tests except this one.
        VERIFY(Solver%Stopping_Test=='||r||/||b|| < Eps',5)
  
        ! Solve the linear system.
  
        call LAMG_Solve (A_LAMG, B_LAMG, X_LAMG, LAMG_Communicator, &
                         LAMG_Options, LAMG_Status)
  
        ! Set Solution vector in MV format.
  
        X_MV = X_LAMG
  
        ! Deallocate LAMG vectors and matrix.
  
        call Finalize (B_LAMG)
        call Finalize (X_LAMG)
        call LAMG_Free_dcsr_r (A_LAMG, LAMG_Communicator, LAMG_Options, &
                               LAMG_Status)
  
        ! Terminate LAMG options.
  
        call LAMG_Options_Term (LAMG_Communicator, LAMG_Options, LAMG_Status)
  
        ! Terminate LAMG communication.
  
        call LAMG_Comm_Term (LAMG_Communicator, LAMG_Status)

      ],[

        ! No LAMG available.

        write (6,*) '*********************************************'
        write (6,*) 'Error: The LAMG solver has been requested but'
        write (6,*) 'this executable was not compiled with LAMG.'
        write (6,*) '*********************************************'
        Call Abort ()

      ])

    ! Error check.

    else
      ! No variable match found.
      VERIFY(.false.,0)
    end if

    ! Verify guarantees.

    ifelse(m4_eval(DEBUG_LEVEL >= RESIDUAL_DEBUG_LEVEL), 1, [
      call Duplicate (Residual_MV, X_MV)
      call Residual (Residual_MV, A_ELLM, X_MV, B_MV)
      if (Solver%Stopping_Test=='||r||/||b|| < Eps') then
        VERIFY(Norm(Residual_MV)/Norm(B_MV)<=Solver%Epsilon,dnl
               RESIDUAL_DEBUG_LEVEL)
      end if
      call Finalize (Residual_MV)
    ])

    return
  end subroutine Solve