diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1ce7b08..f90adf7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -10,6 +10,11 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
  - Rectangular matrices can now be assembled by the new method
    `HYPRE.assemble!(::HYPREMatrixAssembler, i::Vector, j::Vector, a::Matrix)` where `i` are
    the rows and `j` the columns. ([#7][github-7])
+### Fixed
+ - All created HYPRE objects (`HYPREMatrix`, `HYPREVector`, and `HYPRESolver`s) are now kept
+   track of internally and explicitly `finalize`d (if they haven't been GC'd) before
+   finalizing HYPRE. This fixes a "race condition" where MPI and/or HYPRE would finalize
+   before these Julia objects are garbage collected and finalized. ([#8][github-8])
 ### Deprecated
  - The method `HYPRE.assemble!(A::HYPREMatrixAssembler, ij::Vector, a::Matrix)` have been
    deprecated in favor of `HYPRE.assemble!(A::HYPREMatrixAssembler, i::Vector, j::Vector,
@@ -37,6 +42,7 @@ Initial release of HYPRE.jl.
 [github-5]: https://github.com/fredrikekre/HYPRE.jl/pull/5
 [github-6]: https://github.com/fredrikekre/HYPRE.jl/pull/6
 [github-7]: https://github.com/fredrikekre/HYPRE.jl/pull/7
+[github-8]: https://github.com/fredrikekre/HYPRE.jl/pull/8
 
 [1.0.0]: https://github.com/fredrikekre/HYPRE.jl/releases/tag/v1.0.0
 [1.1.0]: https://github.com/fredrikekre/HYPRE.jl/compare/v1.0.0...v1.1.0
diff --git a/src/HYPRE.jl b/src/HYPRE.jl
index e54ab0f..3397076 100644
--- a/src/HYPRE.jl
+++ b/src/HYPRE.jl
@@ -38,7 +38,12 @@ function Init(; finalize_atexit=true)
     if finalize_atexit
         # TODO: MPI only calls the finalizer if not exiting due to a Julia exeption. Does
         #       the same reasoning apply here?
-        atexit(HYPRE_Finalize)
+        atexit() do
+            # Finalize any HYPRE objects that are still alive
+            foreach(finalize, keys(Internals.HYPRE_OBJECTS))
+            # Finalize the library
+            HYPRE_Finalize()
+        end
     end
     return nothing
 end
@@ -67,6 +72,7 @@ function HYPREMatrix(comm::MPI.Comm, ilower::Integer,        iupper::Integer,
     A.ijmatrix = ijmatrix_ref[]
     # Attach a finalizer
     finalizer(x -> HYPRE_IJMatrixDestroy(x.ijmatrix), A)
+    push!(Internals.HYPRE_OBJECTS, A => nothing)
     # Set storage type
     @check HYPRE_IJMatrixSetObjectType(A.ijmatrix, HYPRE_PARCSR)
     # Initialize to make ready for setting values
@@ -106,6 +112,7 @@ function HYPREVector(comm::MPI.Comm, ilower::Integer, iupper::Integer)
     b.ijvector = ijvector_ref[]
     # Attach a finalizer
     finalizer(x -> HYPRE_IJVectorDestroy(x.ijvector), b)
+    push!(Internals.HYPRE_OBJECTS, b => nothing)
     # Set storage type
     @check HYPRE_IJVectorSetObjectType(b.ijvector, HYPRE_PARCSR)
     # Initialize to make ready for setting values
diff --git a/src/Internals.jl b/src/Internals.jl
index 61b1738..2cf22c8 100644
--- a/src/Internals.jl
+++ b/src/Internals.jl
@@ -16,4 +16,6 @@ function setup_func end
 function solve_func end
 function to_hypre_data end
 
+const HYPRE_OBJECTS = WeakKeyDict{Any, Nothing}()
+
 end # module Internals
diff --git a/src/solvers.jl b/src/solvers.jl
index e54a26e..4922199 100644
--- a/src/solvers.jl
+++ b/src/solvers.jl
@@ -7,12 +7,14 @@ Abstract super type of all the wrapped HYPRE solvers.
 """
 abstract type HYPRESolver end
 
-function Internals.safe_finalizer(Destroy)
-    # Only calls the Destroy if pointer not C_NULL
-    return function(solver)
-        if solver.solver != C_NULL
-            Destroy(solver.solver)
-            solver.solver = C_NULL
+function Internals.safe_finalizer(Destroy, solver)
+    # Add the solver to object tracker for possible atexit finalizing
+    push!(Internals.HYPRE_OBJECTS, solver => nothing)
+    # Add a finalizer that only calls Destroy if pointer not C_NULL
+    finalizer(solver) do s
+        if s.solver != C_NULL
+            Destroy(s.solver)
+            s.solver = C_NULL
         end
     end
 end
@@ -109,7 +111,7 @@ mutable struct BiCGSTAB <: HYPRESolver
         @check HYPRE_ParCSRBiCGSTABCreate(comm, solver_ref)
         solver.solver = solver_ref[]
         # Attach a finalizer
-        finalizer(Internals.safe_finalizer(HYPRE_ParCSRBiCGSTABDestroy), solver)
+        Internals.safe_finalizer(HYPRE_ParCSRBiCGSTABDestroy, solver)
         # Set the options
         Internals.set_options(solver, kwargs)
         return solver
@@ -157,7 +159,7 @@ mutable struct BoomerAMG <: HYPRESolver
         @check HYPRE_BoomerAMGCreate(solver_ref)
         solver.solver = solver_ref[]
         # Attach a finalizer
-        finalizer(Internals.safe_finalizer(HYPRE_BoomerAMGDestroy), solver)
+        Internals.safe_finalizer(HYPRE_BoomerAMGDestroy, solver)
         # Set the options
         Internals.set_options(solver, kwargs)
         return solver
@@ -202,7 +204,7 @@ mutable struct FlexGMRES <: HYPRESolver
         @check HYPRE_ParCSRFlexGMRESCreate(comm, solver_ref)
         solver.solver = solver_ref[]
         # Attach a finalizer
-        finalizer(Internals.safe_finalizer(HYPRE_ParCSRFlexGMRESDestroy), solver)
+        Internals.safe_finalizer(HYPRE_ParCSRFlexGMRESDestroy, solver)
         # Set the options
         Internals.set_options(solver, kwargs)
         return solver
@@ -285,7 +287,7 @@ mutable struct GMRES <: HYPRESolver
         @check HYPRE_ParCSRGMRESCreate(comm, solver_ref)
         solver.solver = solver_ref[]
         # Attach a finalizer
-        finalizer(Internals.safe_finalizer(HYPRE_ParCSRGMRESDestroy), solver)
+        Internals.safe_finalizer(HYPRE_ParCSRGMRESDestroy, solver)
         # Set the options
         Internals.set_options(solver, kwargs)
         return solver
@@ -330,7 +332,7 @@ mutable struct Hybrid <: HYPRESolver
         @check HYPRE_ParCSRHybridCreate(solver_ref)
         solver.solver = solver_ref[]
         # Attach a finalizer
-        finalizer(Internals.safe_finalizer(HYPRE_ParCSRHybridDestroy), solver)
+        Internals.safe_finalizer(HYPRE_ParCSRHybridDestroy, solver)
         # Set the options
         Internals.set_options(solver, kwargs)
         return solver
@@ -379,7 +381,7 @@ mutable struct ILU <: HYPRESolver
         @check HYPRE_ILUCreate(solver_ref)
         solver.solver = solver_ref[]
         # Attach a finalizer
-        finalizer(Internals.safe_finalizer(HYPRE_ILUDestroy), solver)
+        Internals.safe_finalizer(HYPRE_ILUDestroy, solver)
         # Set the options
         Internals.set_options(solver, kwargs)
         return solver
@@ -426,7 +428,7 @@ mutable struct ParaSails <: HYPRESolver
         @check HYPRE_ParCSRParaSailsCreate(comm, solver_ref)
         solver.solver = solver_ref[]
         # Attach a finalizer
-        finalizer(Internals.safe_finalizer(HYPRE_ParCSRParaSailsDestroy), solver)
+        Internals.safe_finalizer(HYPRE_ParCSRParaSailsDestroy, solver)
         # Set the options
         Internals.set_options(solver, kwargs)
         return solver
@@ -461,7 +463,7 @@ mutable struct PCG <: HYPRESolver
         @check HYPRE_ParCSRPCGCreate(comm, solver_ref)
         solver.solver = solver_ref[]
         # Attach a finalizer
-        finalizer(Internals.safe_finalizer(HYPRE_ParCSRPCGDestroy), solver)
+        Internals.safe_finalizer(HYPRE_ParCSRPCGDestroy, solver)
         # Set the options
         Internals.set_options(solver, kwargs)
         return solver