diff --git a/Source/C/PSMatrix_c.h b/Source/C/PSMatrix_c.h
index 2448272c..3b0e93f8 100644
--- a/Source/C/PSMatrix_c.h
+++ b/Source/C/PSMatrix_c.h
@@ -63,5 +63,6 @@ double MatrixNorm_ps_wrp(const int *ih_this);
 double MeasureAsymmetry_ps_wrp(const int *ih_this);
 void MatrixTrace_ps_wrp(const int *ih_this, double *trace_val);
 int IsIdentity_ps_wrp(const int *ih_this);
+void SymmetrizeMatrix_ps_wrp(int *ih_this);
 
 #endif
diff --git a/Source/CPlusPlus/PSMatrix.cc b/Source/CPlusPlus/PSMatrix.cc
index 255a7499..ea4c89e4 100644
--- a/Source/CPlusPlus/PSMatrix.cc
+++ b/Source/CPlusPlus/PSMatrix.cc
@@ -232,6 +232,11 @@ double Matrix_ps::MeasureAsymmetry() const {
   return MeasureAsymmetry_ps_wrp(ih_this);
 }
 
+//////////////////////////////////////////////////////////////////////////////
+void Matrix_ps::Symmetrize() {
+  SymmetrizeMatrix_ps_wrp(ih_this);
+}
+
 //////////////////////////////////////////////////////////////////////////////
 Matrix_ps::~Matrix_ps() { DestructMatrix_ps_wrp(ih_this); }
 } // namespace NTPoly
diff --git a/Source/CPlusPlus/PSMatrix.h b/Source/CPlusPlus/PSMatrix.h
index f2a14ff9..0fbf20a5 100644
--- a/Source/CPlusPlus/PSMatrix.h
+++ b/Source/CPlusPlus/PSMatrix.h
@@ -169,6 +169,8 @@ class Matrix_ps {
   double MeasureAsymmetry() const;
   //! compute the trace of a matrix.
   double Trace() const;
+  //! symmetrize the matrix
+  void Symmetrize();
 
 public:
   //! Destructor.
diff --git a/Source/Fortran/PSMatrixAlgebraModule.F90 b/Source/Fortran/PSMatrixAlgebraModule.F90
index 8e05ecd3..72d0b4a0 100644
--- a/Source/Fortran/PSMatrixAlgebraModule.F90
+++ b/Source/Fortran/PSMatrixAlgebraModule.F90
@@ -35,6 +35,7 @@ MODULE PSMatrixAlgebraModule
   PUBLIC :: MatrixTrace
   PUBLIC :: SimilarityTransform
   PUBLIC :: MeasureAsymmetry
+  PUBLIC :: SymmetrizeMatrix
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   INTERFACE MatrixSigma
      MODULE PROCEDURE MatrixSigma_ps
@@ -543,6 +544,21 @@ FUNCTION MeasureAsymmetry(this) RESULT(norm_value)
     CALL DestructMatrix(tmat)
 
   END FUNCTION MeasureAsymmetry
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Make the matrix symmetric
+  SUBROUTINE SymmetrizeMatrix(this)
+    !> The matrix to symmetrize
+    TYPE(Matrix_ps), INTENT(INOUT) :: this
+    !! Local variables
+    TYPE(Matrix_ps) :: tmat
+
+    CALL TransposeMatrix(this, tmat)
+    CALL ConjugateMatrix(tmat)
+    CALL IncrementMatrix(tmat, this, alpha_in=1.0_NTREAL)
+    CALL ScaleMatrix(this, 0.5_NTREAL)
+    CALL DestructMatrix(tmat)
+
+  END SUBROUTINE SymmetrizeMatrix
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> Transform a matrix B = P * A * P^-1
   !! This routine will check if P is the identity matrix, and if so
diff --git a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90 b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90
index aff5b8ed..19deea30 100644
--- a/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90
+++ b/Source/Wrapper/PSMatrixAlgebraModule_wrp.F90
@@ -21,6 +21,7 @@ MODULE PSMatrixAlgebraModule_wrp
   PUBLIC :: MatrixNorm_ps_wrp
   PUBLIC :: MatrixTrace_ps_wrp
   PUBLIC :: MeasureAsymmetry_ps_wrp
+  PUBLIC :: SymmetrizeMatrix_ps_wrp
 CONTAINS!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !> Matrix B = alpha*Matrix A + Matrix B (AXPY)
   SUBROUTINE IncrementMatrix_ps_wrp(ih_matA, ih_matB, alpha_in,threshold_in) &
@@ -156,5 +157,15 @@ FUNCTION MeasureAsymmetry_ps_wrp(ih_this) &
     h_this = TRANSFER(ih_this,h_this)
     norm_value = MeasureAsymmetry(h_this%DATA)
   END FUNCTION MeasureAsymmetry_ps_wrp
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+  !> Symmetrize a matrix
+  SUBROUTINE SymmetrizeMatrix_ps_wrp(ih_matA) &
+        & BIND(c,NAME="SymmetrizeMatrix_ps_wrp")
+    INTEGER(KIND=c_int), INTENT(INOUT) :: ih_matA(SIZE_wrp)
+    TYPE(Matrix_ps_wrp) :: h_matA
+
+    h_matA = TRANSFER(ih_matA,h_matA)
+    CALL SymmetrizeMatrix(h_matA%DATA)
+  END SUBROUTINE SymmetrizeMatrix_ps_wrp
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 END MODULE PSMatrixAlgebraModule_wrp
diff --git a/UnitTests/test_psmatrixalgebra.py b/UnitTests/test_psmatrixalgebra.py
index aeb3cc3b..f9a1882c 100644
--- a/UnitTests/test_psmatrixalgebra.py
+++ b/UnitTests/test_psmatrixalgebra.py
@@ -285,7 +285,7 @@ def test_reverse(self):
             self.check_result()
 
     def test_asymmetry(self):
-        '''Test routines to permute a matrix.'''
+        '''Test routines to measure asymmetry of a matrix.'''
         from numpy import inf
         from scipy.linalg import norm
         for param in self.parameters:
@@ -305,6 +305,29 @@ def test_asymmetry(self):
 
             self.check_floats(ref, comp)
 
+    def test_symmetrize(self):
+        '''Test routines to symmetrize a matrix.'''
+        from numpy import inf
+        from scipy.linalg import norm
+        for param in self.parameters:
+            matrix1 = param.create_matrix(snum=1, complex=self.complex1)
+            self.write_matrix(matrix1, self.input_file1)
+
+            comm.barrier()
+
+            if param.sparsity > 0.0:
+                ntmatrix1 = nt.Matrix_ps(self.input_file1, False)
+            else:
+                ntmatrix1 = nt.Matrix_ps(param.rows)
+
+            self.CheckMat = 0.5 * (matrix1 + matrix1.H)
+            ntmatrix1 = nt.Matrix_ps(self.input_file1, False)
+            ntmatrix1.Symmetrize()
+            ntmatrix1.WriteToMatrixMarket(self.result_file)
+            comm.barrier()
+
+            self.check_result()
+
 
 class TestPSMatrixAlgebra_r(TestPSMatrixAlgebra, unittest.TestCase):
     '''Special routines for real algebra'''