round trip check of precision

Source/Fortran/DataTypesModule.F90

@@ -2,14 +2,19 @@
 !> A module to store specifications for basic data types.
 MODULE DataTypesModule
   USE NTMPIModule
-  USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_DOUBLE, C_DOUBLE_COMPLEX, C_LONG
+  USE, INTRINSIC :: ISO_C_BINDING, ONLY : C_DOUBLE, C_DOUBLE_COMPLEX, C_LONG, &
+    & C_FLOAT
   IMPLICIT NONE
   PRIVATE
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> The precision of floating point numbers we will use in this program.
   INTEGER, PARAMETER, PUBLIC :: NTREAL = C_DOUBLE
+  !> A low precision type for mixed-precision application.
+  INTEGER, PARAMETER, PUBLIC :: NTLOWP = C_FLOAT
   !> MPI floating point datatype with the precision we will use in this program.
-  INTEGER, PUBLIC :: MPINTREAL = MPI_DOUBLE_PRECISION
+  INTEGER, PARAMETER, PUBLIC :: MPINTREAL = MPI_DOUBLE_PRECISION
+  !~> MPI low precision floating point for mixed-precision application.
+  INTEGER, PARAMETER, PUBLIC :: MPILOWP = MPI_REAL
   !> The complex numbers we will use in this program.
   INTEGER, PARAMETER, PUBLIC :: NTCOMPLEX = C_DOUBLE_COMPLEX
   !> MPI complex datatype with the precision we will use in this program.

Source/Fortran/SMatrixAlgebraModule.F90

@@ -10,7 +10,7 @@ MODULE SMatrixAlgebraModule
        & ConstructMatrixMemoryPool
   USE SMatrixModule, ONLY: Matrix_lsr, Matrix_lsc, DestructMatrix, CopyMatrix, &
        & TransposeMatrix, ConjugateMatrix, ConstructMatrixFromTripletList, &
-       & ConstructEmptyMatrix
+       & ConstructEmptyMatrix, RoundTripLowP
   USE SVectorModule, ONLY : AddSparseVectors, PairwiseMultiplyVectors
   USE TripletListModule, ONLY: TripletList_r, TripletList_c, SortTripletList, &
        & DestructTripletList, ConstructTripletList
@@ -237,6 +237,7 @@ SUBROUTINE GemmMatrix_lsr(matA, matB, matC, IsATransposed_in, &
          & INTENT(INOUT), TARGET :: blocked_memory_pool_in
     !! Intermediate Data
     TYPE(Matrix_lsr) :: matAB
+    TYPE(Matrix_lsr) :: matAL, matBL
     LOGICAL :: IsATransposed, IsBTransposed
     REAL(NTREAL) :: alpha
     REAL(NTREAL) :: beta
@@ -272,6 +273,7 @@ SUBROUTINE GemmMatrix_lsc(matA, matB, matC, IsATransposed_in, &
          & INTENT(INOUT), TARGET :: blocked_memory_pool_in
     !! Intermediate Data
     TYPE(Matrix_lsc) :: matAB
+    TYPE(Matrix_lsc) :: matAL, matBL
     LOGICAL :: IsATransposed, IsBTransposed
     REAL(NTREAL) :: alpha
     REAL(NTREAL) :: beta

Source/Fortran/SMatrixModule.F90

@@ -1,7 +1,7 @@
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !> A module for handling locally stored CSR matrices.
 MODULE SMatrixModule
-  USE DataTypesModule, ONLY: NTREAL, NTCOMPLEX, NTLONG
+  USE DataTypesModule, ONLY: NTREAL, NTCOMPLEX, NTLONG, NTLOWP
   USE MatrixMarketModule, ONLY : ParseMMHeader, WriteMMSize, WriteMMLine, &
        & MAX_LINE_LENGTH
   USE TripletListModule, ONLY: TripletList_r, TripletList_c, SortTripletList, &
@@ -51,6 +51,7 @@ MODULE SMatrixModule
   PUBLIC :: ConjugateMatrix
   PUBLIC :: PrintMatrix
   PUBLIC :: MatrixToTripletList
+  PUBLIC :: RoundTripLowP
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   INTERFACE ConstructEmptyMatrix
      MODULE PROCEDURE ConstructEmptyMatrixSub_lsr
@@ -123,6 +124,10 @@ MODULE SMatrixModule
      MODULE PROCEDURE ConvertMatrixType_lsrtolsc
      MODULE PROCEDURE ConvertMatrixType_lsctolsr
   END INTERFACE ConvertMatrixType
+  INTERFACE RoundTripLowP
+     MODULE PROCEDURE RoundTripLowP_lsr
+     MODULE PROCEDURE RoundTripLowP_lsc
+  END INTERFACE RoundTripLowP
 CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> A subroutine type wrapper for the constructor.
   PURE SUBROUTINE ConstructEmptyMatrixSub_lsr(this, rows, columns, zero_in)
@@ -603,4 +608,21 @@ SUBROUTINE ConvertMatrixType_lsctolsr(rin, cout)
 
   END SUBROUTINE ConvertMatrixType_lsctolsr
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  PURE SUBROUTINE RoundTripLowP_lsr(matA, matB)
+    TYPE(Matrix_lsr), INTENT(IN) :: matA
+    TYPE(Matrix_lsr), INTENT(INOUT) :: matB
+    INTEGER :: II
+
+    CALL CopyMatrix(matA, matB)
+    DO II = 1, SIZE(matB%values)
+      matB%values(II) = REAL(matB%values(II), KIND=NTLOWP)
+    END DO
+  END SUBROUTINE RoundTripLowP_lsr
+  PURE SUBROUTINE RoundTripLowP_lsc(matA, matB)
+    TYPE(Matrix_lsc), INTENT(IN) :: matA
+    TYPE(Matrix_lsc), INTENT(INOUT) :: matB
+    INTEGER :: II
+
+    CALL CopyMatrix(matA, matB)
+  END SUBROUTINE RoundTripLowP_lsc
 END MODULE SMatrixModule

Source/Fortran/sparse_includes/GemmMatrix.f90

@@ -55,9 +55,12 @@
      sparsity_estimate = 1e-8
   END IF
 
+  CALL RoundTripLowP(matA, matAL)
+  CALL RoundTripLowP(matB, matBL)
+
   !! Decide whether to do dense or sparse version.
   IF (MIN(sparsity_a, sparsity_b) .GT. sparsity_threshold) THEN
-     CALL DenseBranch(matA, matB, matAB, IsATransposed, IsBTransposed, &
+     CALL DenseBranch(matAL, matBL, matAB, IsATransposed, IsBTransposed, &
           & alpha, threshold)
   ELSE
      !! Setup the memory pool
@@ -77,10 +80,10 @@
      END IF
      !! Multiply
      IF (pool_flag) THEN
-        CALL SparseBranch(matA, matB, matAB, IsATransposed, IsBTransposed, &
+        CALL SparseBranch(matAL, matBL, matAB, IsATransposed, IsBTransposed, &
              & alpha, threshold, blocked_memory_pool_in)
      ELSE
-        CALL SparseBranch(matA, matB, matAB, IsATransposed, IsBTransposed, &
+        CALL SparseBranch(matAL, matBL, matAB, IsATransposed, IsBTransposed, &
              & alpha, threshold, blocked_memory_pool)
      END IF
   END IF
@@ -98,4 +101,6 @@
   END IF
 
   CALL DestructMatrix(matAB)
+  CALL DestructMatrix(matAL)
+  CALL DestructMatrix(matBL)
   CALL DestructMatrixMemoryPool(blocked_memory_pool)