G.1.15 Mathematic_Vector Class Unit Test Program

This lightly commented program performs a unit test on the Mathematic_Vector Class.

program Unit_Test

  use Caesar_Data_Structures_Module
  use Caesar_Mathematic_Vector_Class
  use Caesar_Numbers_Module, only: one, two, three, six
  implicit none

  type(Communication_type) :: Comm
  type(Base_Structure_type) :: Row_Structure
  type(Status_type) :: status
  type(character,name_length) :: Locus_Name, MV_Name
  type(integer,1) :: Row_Length_Vector, Rows
  type(integer) :: shift
  type(real,1) :: Values
  type(real) :: A_Norm2, B_Norm2, DotP, NRows_real, PNorm
  type(integer) :: I, NRows
  type(Mathematic_Vector_type) :: A_MV, B_MV

  ! Initializations.

  call Initialize (Comm)
  call Output (Comm)
  call Initialize (status)
  call Initialize (Locus_Name)
  Locus_Name = 'Rows'
  call Initialize (MV_Name)
  MV_Name = 'Variables'
  call Initialize (Row_Length_Vector, NPEs)
  Row_Length_Vector = (/ (i**2, i = 1, NPEs) /)

  ! Local temp vectors.

  call Initialize (Rows, Row_Length_Vector(this_PE))
  call Initialize (Values, Row_Length_Vector(this_PE))

  ! Initialize base structure and mathematic vectors.
 
  call Initialize (Row_Structure, Row_Length_Vector, Locus_Name, status)
  call Initialize (A_MV, Row_Structure, MV_Name, status)
  call Duplicate (B_MV, A_MV, status)
  NRows = Length_Total(Row_Structure)

  ! Testing statements.

  ! Note that there can be no procedure calls inside a loop which
  ! is executed a different number of times on each PE, because 
  ! global communication is done for verifications within a procedure.

  shift = - First_PE(Row_Structure) + 1
  do i = First_PE(Row_Structure), Last_PE(Row_Structure)
    ! Values(i) = global i
    Values(i + shift) = changetype(real, i)
    ! Rows(i) = global i, but reversed in order on each PE
    Rows(i + shift) = Last_PE(Row_Structure) + First_PE(Row_Structure) - i
  end do

  ! Set A_MV via operator overloading.

  A_MV = Values

  ! Set B_MV via explicit call (vector of values) with half of Values, 
  ! then use the Add_Values procedure to add the other half.
  ! Note that B_MV values are in reverse order on each PE.

  Values = Values/2.d0
  call Set_Values (B_MV, Values, Rows)
  call Add_Values (B_MV, Values, Rows)

  ! Output the MVs.

  PNorm = P_Norm(A_MV, 4)
  call Output (A_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
  call Output (B_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))

  ! Check some things that are known because Values(i) = i :

  if (Average(A_MV) /= (Maximum(A_MV)+one)/two) then
    if (this_is_IO_PE) then
      write (6,*) 'Error -- Average or Maximum is incorrect.'
    end if
  end if  
  if (Maximum(A_MV) /= Infinity_Norm(A_MV)) then
    if (this_is_IO_PE) then
      write (6,*) 'Error -- Maximum or Infinity_Norm is incorrect.'
    end if
  end if  
  if (Minimum(A_MV) /= one) then
    if (this_is_IO_PE) then
      write (6,*) 'Error -- Minimum is incorrect.'
    end if
  end if  
  if (Sum(A_MV) /= One_Norm(A_MV)) then
    if (this_is_IO_PE) then
      write (6,*) 'Error -- Sum or One_Norm is incorrect.'
    end if
  end if  

  ! Tweak one value of A_MV and re-output.

  call Set_Values (A_MV, 10.d0, Last_PE(Row_Structure))
  call Output (A_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))

  ! Set A_MV = B_MV (values are reversed on each PE).

  Values = B_MV
  A_MV = Values

  ! Orthogonality test (note A_MV = B_MV from above).

  if (A_MV .Orthogonal. B_MV) then
    if (this_is_IO_PE) then
      write (6,*) 'Error -- Vectors are orthogonal.'
    end if
  end if

  ! Dot product test. 
  ! This makes use of the fact that
  !
  !  --
  !  \     2     1     3      2
  !  /    k   =  - (2 n  + 3 n  + n)
  !  --          6
  ! k=1,n

  A_MV = Values             ! Make sure A_MV = B_MV. Note that value order is 
  B_MV = Values             ! reversed on each PE, but that won't affect the
                            ! dot product.
  A_Norm2 = Two_Norm(A_MV)  ! These are calculated here so that the
  B_Norm2 = Two_Norm(B_MV)  ! Cauchy-Schwartz inequality check will be
                            ! done in the dot product.
  DotP = A_MV .DotProduct. B_MV
  NRows_real = changetype(real,NRows)
  if (this_is_IO_PE) then
    write (6,100) 'Dot product of A_MV and B_MV = ', DotP
100 format (/,a,f15.1)
    if (.not. VeryClose(DotP, &
        (two * NRows_real**3 + three * NRows_real**2 + NRows_real)/six)) then
      write (6,*) 'Error -- Dot product is incorrect.'
    end if
  end if

  ! Also need 
  ! 1: call Add (MV1, MV2, MV3) or call Add_to_MV (MV1, MV2)
  ! 2: MV1 = MV2
  ! 3: call Scale (MV, scalar)
  ! 4: MV1 == MV2

  ! Check state of the mathematic vectors.

  VERIFY(Valid_State(A_MV),0)
  VERIFY(Valid_State(B_MV),0)

  ! Finalize mathematic vectors and communications.

  call Finalize (A_MV)
  call Finalize (B_MV)
  call Finalize (Comm)

end