Next: G.3 Solver Class Code Up: G.2 ELL_Matrix Class Code Previous: G.2.14 Set_Values_ELL_Matrix Procedure
This lightly commented program performs a unit test on the ELL_Matrix Class.
program Unit_Test
use Caesar_Data_Structures_Module
use Caesar_Mathematic_Vector_Class
use Caesar_ELL_Matrix_Class
use Caesar_Numbers_Module, only: zero, half, one, two, three, six, ten
implicit none
type(Communication_type) :: Comm
type(Base_Structure_type) :: Row_Structure, Column_Structure
type(Status_type) :: status
type(character,name_length) :: Locus_Name, ELLM_Name, MV_Name
type(integer,1) :: Row_Length_Vector, Rows
type(integer,2) :: Columns
type(integer) :: shift, first_plus_last_PE_RS
type(real,2) :: Values
type(real,1) :: X_BNV
type(real) :: AvgB, B_Norm2, MaxB, MinB, NRows_real, ResidualNorm1, &
ResidualNorm2
type(integer) :: i, j, NRows, MaxNonzeros
type(ELL_Matrix_type) :: A_ELLM, B_ELLM
type(Mathematic_Vector_type) :: X_MV, B_MV, Residual_MV
! Initializations.
call Initialize (Comm)
call Output (Comm)
call Initialize (status)
call Initialize (Locus_Name)
call Initialize (ELLM_Name)
call Initialize (MV_Name)
call Initialize (Row_Length_Vector, NPEs)
Row_Length_Vector = (/ (i**2, i = 1, NPEs) /)
call Initialize (MaxNonzeros)
MaxNonzeros = MIN(5, SUM(Row_Length_Vector))
! Local temp vectors.
call Initialize (Rows, Row_Length_Vector(this_PE))
call Initialize (Values, Row_Length_Vector(this_PE), MaxNonzeros)
call Initialize (Columns, Row_Length_Vector(this_PE), MaxNonzeros)
call Initialize (X_BNV, Row_Length_Vector(this_PE))
! Initialize base structure, ELL matrices and mathematic vectors.
Locus_Name = 'Rows'
call Initialize (Row_Structure, Row_Length_Vector, Locus_Name, status)
Locus_Name = 'Columns'
call Initialize (Column_Structure, Row_Length_Vector, Locus_Name, status)
ELLM_Name = 'Coefficients'
call Initialize (A_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
ELLM_Name, status)
call Initialize (B_ELLM, MaxNonzeros, Row_Structure, Column_Structure, &
ELLM_Name, status)
MV_Name = 'Variables'
call Initialize (X_MV, Column_Structure, MV_Name, status)
MV_Name = 'Equations'
call Initialize (B_MV, Row_Structure, MV_Name, 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
first_plus_last_PE_RS = Last_PE(Row_Structure) + First_PE(Row_Structure)
do i = First_PE(Row_Structure), Last_PE(Row_Structure)
! Values(i,:) = global i
Values(i + shift,:) = changetype(real, i) / MaxNonzeros
! Rows(i) = global i, but reversed in order on each PE
Rows(i + shift) = first_plus_last_PE_RS - i
! Columns(i,:) = max((global i)-4,1) + (0,1,2,3,4)
do j = 1, MaxNonzeros
Columns(i + shift,j) = MAX(i-4,1) + j-1
end do
end do
VERIFY(Max_Nonzeros(A_ELLM)==MaxNonzeros,2)
! Set A_ELLM all at once:
!
! 1 PE: A_ELLM = [ 1 ]
! >1 PEs: A_ELLM = [ 0.2 0.2 0.2 0.2 0.2 ]
! [ 0.4 0.4 0.4 0.4 0.4 ]
! [ 0.6 0.6 0.6 0.6 0.6 ]
! [ 0.8 0.8 0.8 0.8 0.8 ]
! [ 1.0 1.0 1.0 1.0 1.0 ]
! [ 1.2 1.2 1.2 1.2 1.2 ]
! [ 1.4 1.4 1.4 1.4 1.4 ]
! [ 1.6 1.6 1.6 1.6 1.6 ]
! [ 1.8 1.8 1.8 1.8 1.8 ]
! [ 2.0 2.0 2.0 2.0 2.0 ]
! and matrix N = sum(i**2, i = 1, NPEs)
call Set_Values (A_ELLM, Values, Columns)
! Set B_ELLM via 'partial replacement' call (even though we are
! replacing the whole thing) with half of Values, then use the
! Add_Values procedure to add the other half.
! Note that B_ELLM values are in reverse order on each PE.
Values = Values/2.d0
call Set_Values (B_ELLM, Values, Rows, Columns)
call Add_Values (B_ELLM, Values, Rows, Columns)
! Output the ELLMs.
call Output (A_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Output (B_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! MatVec test.
X_BNV = one ! Stuff X with 1's.
X_MV = X_BNV
call MatVec (A_ELLM, X_MV, B_MV)
call Output (B_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
MaxB = Maximum(B_MV)
MinB = Minimum(B_MV)
AvgB = Average(B_MV)
if (.not. VeryClose(MinB, one) ) then
if (this_is_IO_PE) then
write (6,*) 'Error in Minimum B => ', MinB
end if
end if
if (.not. VeryClose(MaxB, changetype(real,NRows)) ) then
if (this_is_IO_PE) then
write (6,*) 'Error in Maximum B => ', MaxB
end if
end if
if (.not. VeryClose(AvgB, (MaxB+MinB)*half) ) then
if (this_is_IO_PE) then
write (6,*) 'Error in Average B => ', AvgB
end if
end if
! This should use already-communicated values to get the same answer.
call MatVec (A_ELLM, X_MV, B_MV)
call Output (B_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Norm test for resultant of MatVec (which should be {1,2,3,4,...}).
! This makes use of the fact that
!
! --
! \ 2 1 3 2
! / k = - (2 n + 3 n + n)
! -- 6
! k=1,n
B_Norm2 = Two_Norm(B_MV)
NRows_real = changetype(real,NRows)
if (this_is_IO_PE) then
write (6,100) 'Norm test:'
write (6,100) ' ||B||_2 = ', B_Norm2
100 format (/,a,f15.1)
if (.not. VeryClose(B_Norm2**2, &
(two * NRows_real**3 + three * NRows_real**2 + NRows_real)/six)) then
write (6,*) 'Error -- ||B||_2 is incorrect.'
end if
end if
! This MatVec should require new communication.
X_BNV = ten ! Stuff X with 10's.
X_MV = X_BNV
call MatVec (A_ELLM, X_MV, B_MV)
call Output (B_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
if (.not. VeryClose(Average(B_MV), AvgB*ten) ) then
if (this_is_IO_PE) then
write (6,*) 'Error -- Two consecutive MatVecs give incorrect answers.'
end if
end if
! Check some things that are known:
if (.not. VeryClose(Average(A_ELLM), &
changetype(real,half*(NRows+1)/NRows)) ) then
if (this_is_IO_PE) then
write (6,*) 'Error -- Average is incorrect.'
end if
end if
if (.not. Sum(A_ELLM) >= One_Norm(A_ELLM)) then
if (this_is_IO_PE) then
write (6,*) 'Error -- Sum or One_Norm is incorrect.'
end if
end if
if (.not. Sum(A_ELLM) >= Infinity_Norm(A_ELLM)) then
if (this_is_IO_PE) then
write (6,*) 'Error -- Sum or Infinity_Norm is incorrect.'
end if
end if
if (Minimum(A_ELLM) /= one/changetype(real,Max_Nonzeros(A_ELLM))) then
if (this_is_IO_PE) then
write (6,*) 'Error -- Minimum is incorrect.'
end if
end if
! Tweak one value of A_ELLM on each PE and re-output.
call Set_Values (A_ELLM, 10.d0, Last_PE(Row_Structure), &
Last_PE(Row_Structure))
call Output (A_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Check state of the ELL matrices.
VERIFY(Valid_State(A_ELLM),0)
VERIFY(Valid_State(B_ELLM),0)
! Test the Harwell-Boeing reader.
call Finalize (X_MV) ! Note: must finalize X_MV before A_ELLM because
! the MatVec with A_ELLM created an OV in X_MV
! based on the Data_Index in A_ELLM.
call Finalize (B_MV)
call Finalize (A_ELLM)
call Finalize (B_ELLM)
call Finalize (Row_Structure)
call Finalize (Column_Structure)
if (this_is_IO_PE) then
open (unit=8, &
File='source/class/linear_algebra/battery/Augustus_Prob0_MH_K2_8x8x8.hb')
end if
call Read_Harwell_Boeing (A_ELLM, RHS_MV=B_MV, Solution_MV=X_MV, &
Row_Structure=Row_Structure, &
Column_Structure=Column_Structure, Unit=8, &
status=status)
call Output (A_ELLM, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Output (X_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Output (B_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
! Test Residual and MatVec.
call Initialize (Residual_MV, Row_Structure, 'Residual', status)
call Residual (Residual_MV, A_ELLM, X_MV, B_MV)
ResidualNorm1 = Norm(Residual_MV)
call MatVec (A_ELLM, X_MV, B_MV)
if (this_is_IO_PE) then
write (6,101) 'The next mathematic vector should be the same as', &
'the previous one, except for roundoff.'
101 format (/,a,/,a)
end if
call Output (B_MV, MAX(1, NRows/10), MIN(NRows, NRows/10+50))
call Residual (Residual_MV, A_ELLM, X_MV, B_MV)
ResidualNorm2 = Norm(Residual_MV)
if (this_is_IO_PE) then
! Correct value for ResidualNorm1 is:
! 8.145846809344343E-13 on Intel/Absoft
! 8.145895568345390E-13 on Intel/Lahey
if (ResidualNorm1 < 8.14d-13 .or. &
ResidualNorm1 > 8.15d-13) then
write (6,*) 'Residual Norm 1 out of bounds: ', ResidualNorm1
end if
if (.not. VeryClose(ResidualNorm2,zero)) then
write (6,*) 'Error: Residual Norm 2 is nonzero: ', ResidualNorm2
end if
end if
! Also need
! 1: call Add (ELLM1, ELLM2, ELLM3) or call Add_to_ELLM (ELLM1, ELLM2)
! 2: ELLM1 = ELLM2
! 3: call Scale (ELLM, scalar)
! 4: ELLM1 == ELLM2
! Check state of the ELL matrices.
VERIFY(Valid_State(A_ELLM),0)
! Finalize objects and communications.
call Finalize (X_MV) ! Note: must finalize X_MV before A_ELLM because
! the MatVec with A_ELLM created an OV in X_MV
! based on the Data_Index in A_ELLM.
call Finalize (A_ELLM)
call Finalize (Rows)
call Finalize (Values)
call Finalize (Columns)
call Finalize (X_BNV)
call Finalize (B_MV)
call Finalize (Row_Structure)
call Finalize (Column_Structure)
call Finalize (Comm)
end