Memory leak fix in the dense block branch. (#136)

* Memory leak fix in the dense block branch. * Another missed destructor in the Dense Branch * Another fix. * Some further cleanups to avoid memory leaks. This also removes the need to use the preprocessor on this module (can use built-in include instead). * Optimization. * Remove need for explicit transpose in dense branch * Bug fix and better testing.
william-dawson · Feb 3, 2020 · f668b9c · f668b9c
1 parent 8a20f0c
commit f668b9c
Show file tree

Hide file tree

Showing 6 changed files with 221 additions and 59 deletions.
diff --git a/Source/Fortran/DMatrixModule.F90 b/Source/Fortran/DMatrixModule.F90
@@ -5,9 +5,10 @@
 !! performance.
 MODULE DMatrixModule
   USE DataTypesModule, ONLY : NTREAL, NTCOMPLEX
-  USE SMatrixModule, ONLY : Matrix_lsr, Matrix_lsc
+  USE SMatrixModule, ONLY : Matrix_lsr, Matrix_lsc, &
+       & ConstructMatrixFromTripletList
   USE TripletListModule, ONLY : TripletList_r, TripletList_c, &
-       & AppendToTripletList, ConstructTripletList
+       & AppendToTripletList, ConstructTripletList, DestructTripletList
   USE TripletModule, ONLY : Triplet_r, Triplet_c
   IMPLICIT NONE
   PRIVATE
@@ -101,6 +102,7 @@ PURE SUBROUTINE ConstructEmptyMatrixSup_ldr(this, rows, columns)
     !> Columns of the matrix
     INTEGER, INTENT(IN) :: columns
 
+    CALL DestructMatrix(this)
     this = ConstructEmptyMatrix_ldr(rows, columns)
   END SUBROUTINE ConstructEmptyMatrixSup_ldr
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -126,7 +128,7 @@ PURE SUBROUTINE ConstructMatrixDFromS_ldr(sparse_matrix, dense_matrix)
     !! Helper Variables
     TYPE(Triplet_r) :: temporary
 
-#include "dense_includes/ConstructMatrixDFromS.f90"
+    INCLUDE "dense_includes/ConstructMatrixDFromS.f90"
 
   END SUBROUTINE ConstructMatrixDFromS_ldr
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -143,9 +145,7 @@ PURE SUBROUTINE ConstructMatrixSFromD_ldr(dense_matrix, sparse_matrix, &
     TYPE(Triplet_r) :: temporary
     TYPE(TripletList_r) :: temporary_list
 
-#define SMTYPE Matrix_lsr
-#include "dense_includes/ConstructMatrixSFromD.f90"
-#undef SMTYPE
+    INCLUDE "dense_includes/ConstructMatrixSFromD.f90"
 
   END SUBROUTINE ConstructMatrixSFromD_ldr
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -156,7 +156,7 @@ PURE SUBROUTINE CopyMatrix_ldr(matA, matB)
     !> matB = matA
     TYPE(Matrix_ldr), INTENT(INOUT) :: matB
 
-#include "dense_includes/CopyMatrix.f90"
+    INCLUDE "dense_includes/CopyMatrix.f90"
 
   END SUBROUTINE CopyMatrix_ldr
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -209,7 +209,7 @@ PURE SUBROUTINE TransposeMatrix_ldr(matA, matAT)
     !> matAT = matA^T.
     TYPE(Matrix_ldr), INTENT(INOUT) :: matAT
 
-#include "dense_includes/TransposeMatrix.f90"
+    INCLUDE "dense_includes/TransposeMatrix.f90"
 
   END SUBROUTINE TransposeMatrix_ldr
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -227,7 +227,7 @@ PURE SUBROUTINE ComposeMatrix_ldr(mat_array, block_rows, block_columns, &
     !> The composed matrix.
     TYPE(Matrix_ldr), INTENT(INOUT) :: out_matrix
 
-#include "dense_includes/ComposeMatrix.f90"
+    INCLUDE "dense_includes/ComposeMatrix.f90"
 
   END SUBROUTINE ComposeMatrix_ldr
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -247,21 +247,26 @@ PURE SUBROUTINE SplitMatrix_ldr(this, block_rows, block_columns, &
     !> Specifies the size of the columns.
     INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: block_size_column_in
 
-#include "dense_includes/SplitMatrix.f90"
+    INCLUDE "dense_includes/SplitMatrix.f90"
 
   END SUBROUTINE SplitMatrix_ldr
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> A wrapper for multiplying two dense matrices.
-  SUBROUTINE MultiplyMatrix_ldr(MatA,MatB,MatC)
+  SUBROUTINE MultiplyMatrix_ldr(MatA, MatB, MatC, IsATransposed_in, &
+      & IsBTransposed_in)
     !> The first matrix.
     TYPE(Matrix_ldr), INTENT(IN) :: MatA
     !> The second matrix.
     TYPE(Matrix_ldr), INTENT(IN) :: MatB
     !> MatC = MatA*MatB.
     TYPE(Matrix_ldr), INTENT(INOUT) :: MatC
+    !> True if A is already transposed.
+    LOGICAL, OPTIONAL, INTENT(IN) :: IsATransposed_in
+    !> True if B is already transposed.
+    LOGICAL, OPTIONAL, INTENT(IN) :: IsBTransposed_in
     !! Local variables
-    CHARACTER, PARAMETER :: TRANSA = 'N'
-    CHARACTER, PARAMETER :: TRANSB = 'N'
+    CHARACTER :: TRANSA
+    CHARACTER :: TRANSB
     INTEGER :: M
     INTEGER :: N
     INTEGER :: K
@@ -271,16 +276,55 @@ SUBROUTINE MultiplyMatrix_ldr(MatA,MatB,MatC)
     DOUBLE PRECISION, PARAMETER :: BETA = 0.0
     INTEGER :: LDC
 
-    MatC = Matrix_ldr(MatA%rows,MatB%columns)
+    !! Optional Parameters
+    TRANSA = 'N'
+    IF (PRESENT(IsATransposed_in)) THEN
+       IF (IsATransposed_in) THEN
+        TRANSA = 'T'
+       END IF
+    END IF
+    TRANSB = 'N'
+    IF (PRESENT(IsBTransposed_in)) THEN
+       IF (IsBTransposed_in) THEN
+        TRANSB = 'T'
+       END IF
+    END IF
 
     !! Setup Lapack
-    M = MatA%rows
-    N = MatB%columns
-    K = MatA%columns
-    LDA = M
-    LDB = K
+    IF (TRANSA .EQ. 'T') THEN
+      M = MatA%columns
+    ELSE
+      M = MatA%rows
+    END IF
+
+    IF (TRANSB .EQ. 'T') THEN
+      N = MatB%rows
+    ELSE
+      N = MatB%columns
+    END IF
+
+    IF (TRANSA .EQ. 'T') THEN
+      K = MatA%rows
+    ELSE
+      K = MatA%columns
+    END IF
+
+    IF (TRANSA .EQ. 'T') THEN
+      LDA = K
+    ELSE
+      LDA = M
+    END IF
+
+    IF (TRANSB .EQ. 'T') THEN
+      LDB = N
+    ELSE
+      LDB = K
+    END IF
+
     LDC = M
 
+    !! Multiply
+    CALL ConstructEmptyMatrix(MatC, M, N)
     CALL DGEMM(TRANSA, TRANSB, M, N, K, ALPHA, MatA%data, LDA, MatB%data, &
          & LDB, BETA, MatC%data, LDC)
 
@@ -295,6 +339,7 @@ PURE SUBROUTINE ConstructEmptyMatrixSup_ldc(this, rows, columns)
     !> The number of columns o the matrix.
     INTEGER, INTENT(IN) :: columns
 
+    CALL DestructMatrix(this)
     this = ConstructEmptyMatrix_ldc(rows, columns)
   END SUBROUTINE ConstructEmptyMatrixSup_ldc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -320,7 +365,7 @@ PURE SUBROUTINE ConstructMatrixDFromS_ldc(sparse_matrix, dense_matrix)
     !! Helper Variables
     TYPE(Triplet_c) :: temporary
 
-#include "dense_includes/ConstructMatrixDFromS.f90"
+    INCLUDE "dense_includes/ConstructMatrixDFromS.f90"
 
   END SUBROUTINE ConstructMatrixDFromS_ldc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -337,9 +382,7 @@ PURE SUBROUTINE ConstructMatrixSFromD_ldc(dense_matrix, sparse_matrix, &
     TYPE(Triplet_c) :: temporary
     TYPE(TripletList_c) :: temporary_list
 
-#define SMTYPE Matrix_lsc
-#include "dense_includes/ConstructMatrixSFromD.f90"
-#undef SMTYPE
+    INCLUDE "dense_includes/ConstructMatrixSFromD.f90"
 
   END SUBROUTINE ConstructMatrixSFromD_ldc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -350,7 +393,7 @@ PURE SUBROUTINE CopyMatrix_ldc(matA, matB)
     !> matB = matA
     TYPE(Matrix_ldc), INTENT(INOUT) :: matB
 
-#include "dense_includes/CopyMatrix.f90"
+    INCLUDE "dense_includes/CopyMatrix.f90"
 
   END SUBROUTINE CopyMatrix_ldc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -387,12 +430,14 @@ FUNCTION MatrixNorm_ldc(this) RESULT(norm)
     REAL(NTREAL) :: norm
     !! Local Variables
     INTEGER :: II, JJ
+    COMPLEX(NTCOMPLEX) :: val, conjval
 
     norm = 0
     DO II =1, this%rows
        DO JJ = 1,  this%columns
-          norm = norm + &
-               & REAL(this%data(II,JJ)*CONJG(this%data(II,JJ)),KIND=NTREAL)
+          val = this%data(II,JJ)
+          conjval = CONJG(val)
+          norm = norm + REAL(val*conjval,KIND=NTREAL)
        END DO
     END DO
   END FUNCTION MatrixNorm_ldc
@@ -404,7 +449,7 @@ PURE SUBROUTINE TransposeMatrix_ldc(matA, matAT)
     !> matAT = matA^T.
     TYPE(Matrix_ldc), INTENT(INOUT) :: matAT
 
-#include "dense_includes/TransposeMatrix.f90"
+    INCLUDE "dense_includes/TransposeMatrix.f90"
 
   END SUBROUTINE TransposeMatrix_ldc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -422,7 +467,7 @@ PURE SUBROUTINE ComposeMatrix_ldc(mat_array, block_rows, block_columns, &
     !> The composed matrix.
     TYPE(Matrix_ldc), INTENT(INOUT) :: out_matrix
 
-#include "dense_includes/ComposeMatrix.f90"
+    INCLUDE "dense_includes/ComposeMatrix.f90"
 
   END SUBROUTINE ComposeMatrix_ldc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -442,21 +487,26 @@ PURE SUBROUTINE SplitMatrix_ldc(this, block_rows, block_columns, &
     !> Specifies the size of the columns.
     INTEGER, DIMENSION(:), INTENT(IN), OPTIONAL :: block_size_column_in
 
-#include "dense_includes/SplitMatrix.f90"
+    INCLUDE "dense_includes/SplitMatrix.f90"
 
   END SUBROUTINE SplitMatrix_ldc
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> A wrapper for multiplying two dense matrices.
-  SUBROUTINE MultiplyMatrix_ldc(MatA,MatB,MatC)
+  SUBROUTINE MultiplyMatrix_ldc(MatA, MatB, MatC, IsATransposed_in, &
+      & IsBTransposed_in)
     !> The first matrix.
     TYPE(Matrix_ldc), INTENT(IN) :: MatA
     !> The second matrix.
     TYPE(Matrix_ldc), INTENT(IN) :: MatB
     !> MatC = MatA*MatB.
     TYPE(Matrix_ldc), INTENT(INOUT) :: MatC
+    !> True if A is already transposed.
+    LOGICAL, OPTIONAL, INTENT(IN) :: IsATransposed_in
+    !> True if B is already transposed.
+    LOGICAL, OPTIONAL, INTENT(IN) :: IsBTransposed_in
     !! Local variables
-    CHARACTER, PARAMETER :: TRANSA = 'N'
-    CHARACTER, PARAMETER :: TRANSB = 'N'
+    CHARACTER :: TRANSA
+    CHARACTER :: TRANSB
     INTEGER :: M
     INTEGER :: N
     INTEGER :: K
@@ -466,16 +516,55 @@ SUBROUTINE MultiplyMatrix_ldc(MatA,MatB,MatC)
     COMPLEX*16, PARAMETER :: BETA = 0.0
     INTEGER :: LDC
 
-    MatC = Matrix_ldc(MatA%rows,MatB%columns)
+    !! Optional Parameters
+    TRANSA = 'N'
+    IF (PRESENT(IsATransposed_in)) THEN
+       IF (IsATransposed_in) THEN
+        TRANSA = 'T'
+       END IF
+    END IF
+    TRANSB = 'N'
+    IF (PRESENT(IsBTransposed_in)) THEN
+       IF (IsBTransposed_in) THEN
+        TRANSB = 'T'
+       END IF
+    END IF
 
     !! Setup Lapack
-    M = MatA%rows
-    N = MatB%columns
-    K = MatA%columns
-    LDA = M
-    LDB = K
+    IF (TRANSA .EQ. 'T') THEN
+      M = MatA%columns
+    ELSE
+      M = MatA%rows
+    END IF
+
+    IF (TRANSB .EQ. 'T') THEN
+      N = MatB%rows
+    ELSE
+      N = MatB%columns
+    END IF
+
+    IF (TRANSA .EQ. 'T') THEN
+      K = MatA%rows
+    ELSE
+      K = MatA%columns
+    END IF
+
+    IF (TRANSA .EQ. 'T') THEN
+      LDA = K
+    ELSE
+      LDA = M
+    END IF
+
+    IF (TRANSB .EQ. 'T') THEN
+      LDB = N
+    ELSE
+      LDB = K
+    END IF
+
     LDC = M
 
+    !! Multiply
+    CALL ConstructEmptyMatrix(MatC, M, N)
     CALL ZGEMM(TRANSA, TRANSB, M, N, K, ALPHA, MatA%data, LDA, MatB%data, &
          & LDB, BETA, MatC%data, LDC)
 

diff --git a/Source/Fortran/SMatrixAlgebraModule.F90 b/Source/Fortran/SMatrixAlgebraModule.F90
@@ -3,7 +3,8 @@
 MODULE SMatrixAlgebraModule
   USE DataTypesModule, ONLY : NTREAL, NTCOMPLEX
   USE DMatrixModule, ONLY : Matrix_ldr, Matrix_ldc, ConstructMatrixDFromS, &
-       & ConstructMatrixSFromD, MultiplyMatrix, DestructMatrix
+       & ConstructMatrixSFromD, CopyMatrix, MultiplyMatrix, TransposeMatrix, &
+       & DestructMatrix
   USE MatrixMemoryPoolModule, ONLY : MatrixMemoryPool_lr, MatrixMemoryPool_lc, &
        & DestructMatrixMemoryPool, CheckMemoryPoolValidity, SetPoolSparsity, &
        & ConstructMatrixMemoryPool
@@ -420,8 +421,6 @@ SUBROUTINE DenseBranch_lsr(matA, matB, matC, IsATransposed, IsBTransposed, &
     !> Threshold for flushing values.
     REAL(NTREAL), INTENT(IN) :: threshold
     !! Local Data
-    TYPE(Matrix_lsr) :: untransposedMatA
-    TYPE(Matrix_lsr) :: untransposedMatB
     TYPE(Matrix_ldr) :: DenseA
     TYPE(Matrix_ldr) :: DenseB
     TYPE(Matrix_ldr) :: DenseC
@@ -447,8 +446,6 @@ SUBROUTINE DenseBranch_lsc(matA, matB, matC, IsATransposed, IsBTransposed, &
     !> Threshold for flushing values.
     REAL(NTREAL), INTENT(IN) :: threshold
     !! Local Data
-    TYPE(Matrix_lsc) :: untransposedMatA
-    TYPE(Matrix_lsc) :: untransposedMatB
     TYPE(Matrix_ldc) :: DenseA
     TYPE(Matrix_ldc) :: DenseB
     TYPE(Matrix_ldc) :: DenseC

diff --git a/Source/Fortran/dense_includes/ComposeMatrix.f90 b/Source/Fortran/dense_includes/ComposeMatrix.f90
@@ -22,6 +22,7 @@
   !! Allocate Memory
   CALL ConstructEmptyMatrix(out_matrix, out_columns, out_rows)
 
+  !! Copy
   DO JJ = 1, block_columns
      DO II = 1, block_rows
         out_matrix%data(row_offsets(II):row_offsets(II+1)-1, &

diff --git a/Source/Fortran/dense_includes/ConstructMatrixSFromD.f90 b/Source/Fortran/dense_includes/ConstructMatrixSFromD.f90
@@ -32,4 +32,6 @@
      END DO
   END IF
 
-  sparse_matrix = SMTYPE(temporary_list, rows, columns)
+  CALL ConstructMatrixFromTripletList(sparse_matrix, temporary_list, &
+       & rows, columns)
+  CALL DestructTripletList(temporary_list)
diff --git a/Source/Fortran/sparse_includes/DenseBranch.f90 b/Source/Fortran/sparse_includes/DenseBranch.f90
@@ -1,21 +1,10 @@
-  !! Handle Transposed Case
-  IF (IsATransposed) THEN
-     CALL TransposeMatrix(matA,untransposedMatA)
-  ELSE
-     CALL CopyMatrix(matA,untransposedMatA)
-  END IF
-  IF (IsBTransposed) THEN
-     CALL TransposeMatrix(matB,untransposedMatB)
-  ELSE
-     CALL CopyMatrix(matB,untransposedMatB)
-  END IF
-
   !! Convert Forward
-  CALL ConstructMatrixDFromS(untransposedMatA, DenseA)
-  CALL ConstructMatrixDFromS(untransposedMatB, DenseB)
+  CALL ConstructMatrixDFromS(matA, DenseA)
+  CALL ConstructMatrixDFromS(matB, DenseB)
 
   !! Multiply
-  CALL MultiplyMatrix(DenseA, DenseB, DenseC)
+  CALL MultiplyMatrix(DenseA, DenseB, DenseC, &
+       & IsATransposed_in = IsATransposed, IsBTransposed_in = IsBTransposed)
 
   !! Convert Back
   CALL ConstructMatrixSFromD(DenseC, matC, threshold)