diff --git a/src/main.f08 b/src/main.f08
index 6a930a20..684f4ee9 100644
--- a/src/main.f08
+++ b/src/main.f08
@@ -8,7 +8,7 @@
 !! KU Leuven, Belgium.
 program legolas
   use mod_global_variables, only: dp, str_len, initialise_globals
-  use mod_matrix_structure, only: matrix_t
+  use mod_matrix_structure, only: matrix_t, new_matrix
   use mod_equilibrium, only: set_equilibrium
   use mod_inspections, only: do_equilibrium_inspections
   use mod_matrix_manager, only: build_matrices
@@ -48,7 +48,10 @@ program legolas
   settings = new_settings()
 
   call timer%start_timer()
-  call initialisation()
+
+  call read_user_parfile()
+  call print_startup_info()
+
   grid = new_grid(settings)
   background = new_background()
   physics = new_physics(settings, background)
@@ -60,13 +63,17 @@ program legolas
   call do_equilibrium_inspections(settings, grid, background, physics)
 
   call timer%start_timer()
+  matrix_A = new_matrix(nb_rows=settings%dims%get_dim_matrix(), label="A")
+  matrix_B = new_matrix(nb_rows=settings%dims%get_dim_matrix(), label="B")
   call build_matrices(matrix_B, matrix_A, settings, grid, background, physics)
   timer%matrix_time = timer%end_timer()
 
   call logger%info("solving eigenvalue problem...")
   call timer%start_timer()
+  call do_eigenvalue_problem_allocations()
   call solve_evp(matrix_A, matrix_B, settings, omega, right_eigenvectors)
   timer%evp_time = timer%end_timer()
+  call logger%info("done.")
 
   call timer%start_timer()
   eigenfunctions = new_eigenfunctions(settings, grid, background)
@@ -98,23 +105,24 @@ program legolas
 
 contains
 
-  !> Subroutine responsible for all initialisations.
-  !! Allocates and initialises main and global variables, then the equilibrium state
-  !! and eigenfunctions are initialised and the equilibrium is set.
-  subroutine initialisation()
-    use mod_matrix_structure, only: new_matrix
-    use mod_input, only: read_parfile, get_parfile
-    use mod_console, only: print_logo
-
-    character(len=5*str_len)  :: parfile
-    integer   :: nb_evs
-
+  subroutine read_user_parfile()
+    use mod_input, only: get_parfile, read_parfile
+    character(len=5*str_len) :: parfile
     call get_parfile(parfile)
     call read_parfile(parfile, settings)
+  end subroutine read_user_parfile
 
+
+  subroutine print_startup_info()
+    use mod_console, only: print_logo
     call print_logo()
     call logger%info("the physics type is " // settings%get_physics_type())
     call logger%info("the state vector is " // str(settings%get_state_vector()))
+  end subroutine print_startup_info
+
+
+  subroutine do_eigenvalue_problem_allocations()
+    integer :: nb_evs
 
     select case(settings%solvers%get_solver())
     case ("arnoldi")
@@ -124,10 +132,8 @@ subroutine initialisation()
     case default
       nb_evs = settings%dims%get_dim_matrix()
     end select
-    call logger%debug("setting #eigenvalues to " // str(nb_evs))
+    call logger%debug("allocating eigenvalue array of size " // str(nb_evs))
     allocate(omega(nb_evs))
-    matrix_A = new_matrix(nb_rows=settings%dims%get_dim_matrix(), label="A")
-    matrix_B = new_matrix(nb_rows=settings%dims%get_dim_matrix(), label="B")
 
     ! Arnoldi solver needs this, since it always calculates an orthonormal basis
     if ( &
@@ -144,7 +150,7 @@ subroutine initialisation()
       call logger%debug("allocating eigenvector arrays as dummy")
       allocate(right_eigenvectors(2, 2))
     end if
-  end subroutine initialisation
+  end subroutine do_eigenvalue_problem_allocations
 
 
   !> Initialises and calculates the eigenfunctions if requested.
@@ -158,16 +164,16 @@ end subroutine get_eigenfunctions
   !> Deallocates all main variables, then calls the cleanup
   !! routines of all relevant subroutines to do the same thing.
   subroutine cleanup()
-    call matrix_A%delete_matrix()
-    call matrix_B%delete_matrix()
     deallocate(omega)
     if (allocated(right_eigenvectors)) deallocate(right_eigenvectors)
 
-    call settings%delete()
+    call matrix_A%delete_matrix()
+    call matrix_B%delete_matrix()
+    call eigenfunctions%delete()
+    call physics%delete()
     call grid%delete()
     call background%delete()
-    call physics%delete()
-    call eigenfunctions%delete()
+    call settings%delete()
   end subroutine cleanup
 
 
diff --git a/src/settings/mod_io_settings.f08 b/src/settings/mod_io_settings.f08
index 5de91c5c..ec10b100 100644
--- a/src/settings/mod_io_settings.f08
+++ b/src/settings/mod_io_settings.f08
@@ -93,7 +93,6 @@ end function should_compute_eigenvectors
 
   pure subroutine set_all_io_to_false(this)
     class(io_settings_t), intent(inout) :: this
-    this%write_matrices = .false.
     this%write_eigenvectors = .false.
     this%write_residuals = .false.
     this%write_eigenfunctions = .false.