G.2.9 MatVec_ELL_Matrix Procedure

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

  subroutine MatVec_ELL_Matrix (ELLM, MV_in, MV_out)
  
    ! Input variables.

    ! ELL Matrix to be multiplied.
    type(ELL_Matrix_type), intent(inout) :: ELLM
    ! Mathematic Vector to be multiplied.
    type(Mathematic_Vector_type), intent(inout) :: MV_in

    ! Output variable.

    ! Result of dot product.
    type(Mathematic_Vector_type), intent(inout) :: MV_out 

    ! Internal variables.

    type(integer) :: i                     ! Loop counter.
    type(integer) :: OV_match              ! The OV match for ELLM.
    type(real,2) :: MV_in_Values_Array     ! Temp for collected MV_in Values.
    type(real) :: random_real              ! Random temporary.

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  
    ! Verify requirements.
  
    VERIFY(Valid_State(ELLM),5)            ! ELLM is valid.
    VERIFY(Valid_State(MV_in),5)           ! MV_in is valid.
    VERIFY(Valid_State(MV_out),5)          ! MV_out is valid.
    ! MV_in has same Structure as ELLM Column_Structure.
    VERIFY(ASSOCIATED(ELLM%Column_Structure,MV_in%Structure),5) 
    ! MV_out has same Structure as ELLM Row_Structure.
    VERIFY(ASSOCIATED(ELLM%Row_Structure,MV_out%Structure),5)
    VERIFY(Number_of_OVs_in_an_MV==4,5)    ! 4 is assumed in this procedure.

    ! Check to see if Index needs to be updated.

    if (.not. ELLM%Index_is_Updated) then
      call Finalize (ELLM%Index)
      call Initialize (ELLM%Index, ELLM%Column_Structure, &
                       ELLM%Row_Structure, Many_of_One_Array=ELLM%Columns)
      ELLM%Index_is_Updated = .true.
      call RANDOM_NUMBER (random_real)
      ELLM%Index_Match_Number = changetype(integer, random_real*1.e9)
    end if

    ! ### SPRONG -- need to step through the logic in this routine, 
    ! ### especially to make sure that MV_in%DV is updated correctly 
    ! ### (and not updated when it shouldn't be).

    ! Check to see if MV_in already has an Overlapped_Vector for this Matrix.

    OV_match = 0
    do i = 1, Number_of_OVs_in_an_MV
      if (Initialized(MV_in%OV(i))) then
       if (ELLM%Index_Match_Number == MV_in%Index_Match_Number(i)) then
          OV_match = i
        end if
      end if
    end do

    ! Create MV_in%OV, MV_in%DV if needed.

    if (OV_match == 0) then

      ! Determine which slot to put new OV in.

      do i = 1, Number_of_OVs_in_an_MV
        if (Initialized(MV_in%OV(i))) then
          OV_match = i + 1
        end if
      end do
      if (OV_match == 0 .or. OV_match == 5) then
        OV_match = 1
      end if

      ! Set up MV_in%DV if necessary.

      if (.not.Initialized(MV_in%DV)) then
        call Initialize (MV_in%DV, MV_in%Structure, MV_in%Dimensionality, &
                         MV_in%Name)
      end if

      ! Set up MV_in%OV(OV_match).

      if (Initialized(MV_in%OV(OV_match))) call Finalize (MV_in%OV(OV_match))
      MV_in%DV = MV_in%Values
      call Initialize (MV_in%OV(OV_match), MV_in%DV, ELLM%Index, MV_in%Name)
      MV_in%Index_Match_Number(OV_match) = ELLM%Index_Match_Number

    end if

    ! Update MV_in%DV and MV_in%OV if needed.

    call Update_DV (MV_in)
    MV_in%OV(OV_match) = MV_in%DV

    ! Extract the values from MV_in%OV.

    call Initialize (MV_in_Values_Array, Length_PE(MV_in%Structure), &
                     ELLM%Max_Nonzeros)
    MV_in_Values_Array = MV_in%OV(OV_match)

    ! Calculate the global matrix-vector product.

    MV_out%Values(:) = SUM(ELLM%Values * MV_in_Values_Array, 2) 

    ! Unset the updated? variables.

    call Set_Not_Up_to_Date (MV_out)

    ! Clean up.

    call Finalize (MV_in_Values_Array)

    ! Verify guarantees - none.

    return
  end subroutine MatVec_ELL_Matrix