diff --git a/lib/heuristics/dichotomious_approach.rb b/lib/heuristics/dichotomious_approach.rb
index 1a5fe362d..1ec4485e1 100644
--- a/lib/heuristics/dichotomious_approach.rb
+++ b/lib/heuristics/dichotomious_approach.rb
@@ -89,14 +89,10 @@ def self.dichotomious_heuristic(service_vrp, job = nil, &block)
             if sub_service_vrp[:vrp].resolution_minimum_duration
               sub_service_vrp[:vrp].resolution_minimum_duration *= sub_service_vrp[:vrp].services.size / service_vrp[:vrp].services.size.to_f * 2
             end
-            matrix_indices = sub_service_vrp[:vrp].points.map{ |point|
-              service_vrp[:vrp].points.find{ |r_point| point.id == r_point.id }.matrix_index
-            }
-            SplitClustering.update_matrix_index(sub_service_vrp[:vrp])
-            SplitClustering.update_matrix(service_vrp[:vrp].matrices, sub_service_vrp[:vrp], matrix_indices)
+
             result = OptimizerWrapper.define_process(sub_service_vrp, job, &block)
 
-            transfer_unused_vehicles(result, sub_service_vrps) if index.zero? && result
+            transfer_unused_vehicles(service_vrp, result, sub_service_vrps) if index.zero? && result
 
             result
           }
@@ -129,27 +125,23 @@ def self.dichotomious_heuristic(service_vrp, job = nil, &block)
       result
     end
 
-    def self.transfer_unused_vehicles(result, sub_service_vrps)
+    def self.transfer_unused_vehicles(service_vrp, result, sub_service_vrps)
+      original_vrp = service_vrp[:vrp]
       sv_zero = sub_service_vrps[0][:vrp]
       sv_one = sub_service_vrps[1][:vrp]
 
-      # First transfer empty vehicles that appear in the routes
-      result[:routes].each{ |r|
-        next if !r[:activities].select{ |a| a[:service_id] }.empty?
-
-        vehicle = sv_zero.vehicles.find{ |v| v.id == r[:vehicle_id] }
-        sv_one.vehicles << vehicle
-        sv_zero.vehicles -= [vehicle]
-        sv_one.points += sv_zero.points.select{ |p| p.id == vehicle.start_point_id || p.id == vehicle.end_point_id }
-      }
-
-      # Then transfer the vehicles which do not appear in the routes.
+      original_matrix_indices = nil
       sv_zero.vehicles.each{ |vehicle|
-        next if result[:routes].any?{ |r| r[:vehicle_id] == vehicle.id }
+        route = result[:routes].find{ |r| r[:vehicle_id] == vehicle.id }
+
+        # Transfer the vehicles which do not appear in the routes or the empty vehicles that appear in the routes
+        next if route && !route[:activities].select{ |a| a[:service_id] }.empty?
 
         sv_one.vehicles << vehicle
         sv_zero.vehicles -= [vehicle]
-        sv_one.points += sv_zero.points.select{ |p| p.id == vehicle.start_point_id || p.id == vehicle.end_point_id }
+        vehicle_points = sv_zero.points.select{ |p| p.id == vehicle.start_point_id || p.id == vehicle.end_point_id }
+
+        update_sv_one_matrix(sv_one, original_vrp, original_matrix_indices, vehicle, vehicle_points)
       }
 
       # Transfer unsued vehicle limit to the other side as well
@@ -158,6 +150,45 @@ def self.transfer_unused_vehicles(result, sub_service_vrps)
       sv_one.resolution_vehicle_limit += sv_zero_unused_vehicle_limit
     end
 
+    def self.update_sv_one_matrix(sv_one, original_vrp, original_matrix_indices, vehicle, vehicle_points)
+      vehicle_points.each{ |new_point|
+        point_exists = sv_one.points.find{ |p| p.id == new_point.id }
+
+        if point_exists
+          vehicle.start_point = point_exists if vehicle.start_point_id == new_point.id
+          vehicle.end_point = point_exists if vehicle.end_point_id == new_point.id
+          next
+        end
+
+        new_point.matrix_index = sv_one.points.size
+
+        original_matrix_indices ||= sv_one.points.map{ |p| original_vrp.points.find{ |pi| pi.id == p.id }.matrix_index }
+        new_point_original_matrix_index = original_vrp.points.find{ |pi| pi.id == new_point.id }.matrix_index
+
+        # Update the matrix
+        sv_one.matrices.each_with_index{ |sv_one_matrix, m_index|
+          %i[time distance value].each{ |dimension|
+            d_matrix = sv_one_matrix.send(dimension)
+            next unless d_matrix
+
+            original_matrix = original_vrp.matrices[m_index].send(dimension)
+
+            # existing points to new_point
+            d_matrix.each_with_index{ |row, r_index|
+              row << original_matrix[original_matrix_indices[r_index]][new_point_original_matrix_index]
+            }
+            # new_point to existing points
+            d_matrix << original_matrix[new_point_original_matrix_index].values_at(*original_matrix_indices)
+            # new_point to new_point
+            d_matrix.last << 0
+          }
+        }
+
+        original_matrix_indices << new_point_original_matrix_index
+        sv_one.points << new_point
+      }
+    end
+
     def self.dicho_level_coeff(service_vrp)
       balance = 0.66666
       level_approx = Math.log(service_vrp[:vrp].resolution_dicho_division_vehicle_limit / (service_vrp[:vrp].resolution_vehicle_limit || service_vrp[:vrp].vehicles.size).to_f, balance)
diff --git a/lib/interpreters/split_clustering.rb b/lib/interpreters/split_clustering.rb
index e829888e1..5cd3d3379 100644
--- a/lib/interpreters/split_clustering.rb
+++ b/lib/interpreters/split_clustering.rb
@@ -600,10 +600,13 @@ def self.remove_poorly_populated_routes(vrp, result, limit)
       log 'Some routes are emptied due to poor workload -- time or quantity.', level: :warn if emptied_routes
     end
 
-    def self.update_matrix(original_matrices, sub_vrp, matrix_indices)
-      sub_vrp.matrices.each_with_index{ |matrix, index|
-        [:time, :distance].each{ |dimension|
-          matrix[dimension] = sub_vrp.vehicles.first.matrix_blend(original_matrices[index], matrix_indices, [dimension], cost_time_multiplier: 1, cost_distance_multiplier: 1)
+    def self.update_matrix(sub_vrp, matrix_indices)
+      sub_vrp.matrices.each{ |matrices|
+        [:time, :distance, :value].each{ |dimension|
+          matrix = matrices.send(dimension)
+          next unless matrix
+
+          matrices.send("#{dimension}=", matrix_indices.map{ |r_index| matrix[r_index].values_at(*matrix_indices) })
         }
       }
     end
@@ -643,26 +646,32 @@ def self.build_partial_service_vrp(service_vrp, partial_service_ids, available_v
             vehicle.id == r.vehicle_id && (route_week_day.nil? || vehicle_week_day_availability.include?(route_week_day))
           }
         }
-        sub_vrp.matrices.delete_if{ |m|
-          sub_vrp.vehicles.none?{ |vehicle| vehicle.matrix_id == m.id }
-        }
       end
-      sub_vrp.services = sub_vrp.services.select{ |service| partial_service_ids.include?(service.id) }.compact
-      points_ids = sub_vrp.services.map{ |s| s.activity.point.id }.uniq.compact
-      sub_vrp.rests = sub_vrp.rests.select{ |r| sub_vrp.vehicles.flat_map{ |v| v.rests.map(&:id) }.include? r.id }
-      sub_vrp.relations = sub_vrp.relations.select{ |r| r.linked_ids.all? { |id| sub_vrp.services.any? { |s| s.id == id } } }
+      sub_vrp.services.select!{ |service| partial_service_ids.include?(service.id) }
+      rest_ids = sub_vrp.vehicles.flat_map{ |v| v.rests.map(&:id) }.uniq
+      sub_vrp.rests.select!{ |r| rest_ids.include?(r.id) }
       if entity == :vehicle
         sub_vrp.relations.delete_if{ |r| r.type == :same_vehicle }
         sub_vrp.services.each{ |s| s.relations.delete_if{ |r| r.type == :same_vehicle } }
       end
-      sub_vrp.points = sub_vrp.points.select{ |p| points_ids.include? p.id }.compact
-      sub_vrp.points += sub_vrp.vehicles.flat_map{ |vehicle| [vehicle.start_point, vehicle.end_point] }.compact.uniq
+      sub_vrp.relations.select!{ |r|
+        # Split respects relations, it is enough to check only the first linked id --  [0..0].any? is to handle empties
+        r.linked_ids[0..0].any? { |sid| sub_vrp.services.any? { |s| s.id == sid } } &&
+          r.linked_vehicle_ids[0..0].any? { |vid| sub_vrp.vehicles.any? { |v| v.id == vid } }
+      }
+      points_ids = sub_vrp.services.map{ |s| s.activity.point.id }.compact |
+                   sub_vrp.vehicles.flat_map{ |vehicle| [vehicle.start_point_id, vehicle.end_point_id] }.compact
+      sub_vrp.points.select!{ |p| points_ids.include?(p.id) }
+      sub_vrp.vehicles.each{ |vehicle|
+        vehicle.start_point = sub_vrp.points.find{ |p| p.id == vehicle.start_point_id } if vehicle.start_point
+        vehicle.end_point = sub_vrp.points.find{ |p| p.id == vehicle.end_point_id } if vehicle.end_point
+      }
       sub_vrp = add_corresponding_entity_skills(entity, sub_vrp)
 
       if !sub_vrp.matrices&.empty?
         matrix_indices = sub_vrp.points.map{ |point| point.matrix_index }
         update_matrix_index(sub_vrp)
-        update_matrix(sub_vrp.matrices, sub_vrp, matrix_indices)
+        update_matrix(sub_vrp, matrix_indices)
       end
 
       log "<--- build_partial_service_vrp takes #{Time.now - tic}", level: :debug
diff --git a/test/api/v01/vrp_test.rb b/test/api/v01/vrp_test.rb
index d4a105bb4..6cae9df81 100644
--- a/test/api/v01/vrp_test.rb
+++ b/test/api/v01/vrp_test.rb
@@ -265,6 +265,7 @@ def test_block_call_under_clustering
     vrp1[:configuration][:preprocessing][:partitions] = TestHelper.vehicle_and_days_partitions
 
     vrp2 = VRP.independent_skills
+    vrp2[:matrices][0][:distance] = Oj.load(Oj.dump(vrp2[:matrices][0][:time]))
     vrp2[:points] = VRP.lat_lon_periodic[:points]
     vrp2[:services].first[:skills] = ['D']
     vrp2[:configuration][:preprocessing] = {
diff --git a/test/lib/interpreters/split_clustering_test.rb b/test/lib/interpreters/split_clustering_test.rb
index 7ccbe7f03..abd0acc8c 100644
--- a/test/lib/interpreters/split_clustering_test.rb
+++ b/test/lib/interpreters/split_clustering_test.rb
@@ -459,11 +459,10 @@ def test_correct_number_of_visits_when_concurrent_split_independent_and_max_spli
       # and then max_split again during solution process
       # should not raise "Wrong number of visits returned in result" error
       vrp = VRP.independent_skills
+      vrp[:matrices][0][:distance] = Oj.load(Oj.dump(vrp[:matrices][0][:time]))
       vrp[:points] = VRP.lat_lon_periodic[:points]
       vrp[:services].first[:skills] = ['D']
-      vrp[:configuration][:preprocessing] = {
-        max_split_size: 4
-      }
+      vrp[:configuration][:preprocessing] = { max_split_size: 4 }
 
       OptimizerWrapper.wrapper_vrp('demo', { services: { vrp: [:ortools] }}, TestHelper.create(vrp), nil)
     end
@@ -758,6 +757,7 @@ def test_results_regularity
       vrp = Marshal.dump(TestHelper.load_vrp(self)) # call load_vrp only once to not to dump for each restart
       (1..@regularity_restarts).each{ |trial|
         puts "Regularity trial: #{trial}/#{@regularity_restarts}"
+        Models.delete_all
         result = OptimizerWrapper.wrapper_vrp('ortools', { services: { vrp: [:ortools] }}, Marshal.load(vrp), nil) # rubocop: disable Security/MarshalLoad
         visits_unassigned << result[:unassigned].size
         unassigned_service_ids = result[:unassigned].collect{ |unassigned| unassigned[:original_service_id] }
@@ -801,6 +801,7 @@ def test_results_regularity_2
       vrp = Marshal.dump(TestHelper.load_vrp(self)) # call load_vrp only once to not to dump for each restart
       (1..@regularity_restarts).each{ |trial|
         OptimizerLogger.log "Regularity trial: #{trial}/#{@regularity_restarts}"
+        Models.delete_all
         result = OptimizerWrapper.wrapper_vrp('ortools', { services: { vrp: [:ortools] }}, Marshal.load(vrp), nil) # rubocop: disable Security/MarshalLoad
         visits_unassigned << result[:unassigned].size
         unassigned_service_ids = result[:unassigned].collect{ |unassigned| unassigned[:original_service_id] }
@@ -960,6 +961,17 @@ def test_list_vehicles
       assert_equal 5, Interpreters::SplitClustering.list_vehicles({ start: 0, end: 4 }, vrp.vehicles, :work_day).size
     end
 
+    def test_split_keeps_matrices_in_case_vehicles_are_moved_between_subproblems
+      problem = VRP.lat_lon_two_vehicles
+      problem[:vehicles].last[:matrix_id] = 'm2'
+      problem[:matrices] << problem[:matrices][0].merge({ id: 'm2' })
+      vrp = TestHelper.create(problem)
+
+      sub_vrp = Interpreters::SplitClustering.build_partial_service_vrp({ vrp: vrp }, vrp.services.map(&:id), [0])[:vrp]
+
+      assert_equal %w[m1 m2], sub_vrp.matrices.map(&:id), 'Split should not eliminate matrices in case vehicles are moved between subproblems'
+    end
+
     def test_split_with_vehicle_alternative_skills
       problem = VRP.lat_lon_two_vehicles
 
diff --git a/test/real_cases_periodic_test.rb b/test/real_cases_periodic_test.rb
index e536a982c..fe2f93e1d 100644
--- a/test/real_cases_periodic_test.rb
+++ b/test/real_cases_periodic_test.rb
@@ -127,6 +127,7 @@ def test_minimum_stop_in_route
       assert result[:routes].any?{ |r| r[:activities].size - 2 < 5 },
              'We expect at least one route with less than 5 services, this test is useless otherwise'
       should_remain_assigned = result[:routes].sum{ |r| r[:activities].size - 2 >= 5 ? r[:activities].size - 2 : 0 }
+      should_remain_unassigned = result[:unassigned].size
 
       # one vehicle should have at least 5 stops :
       vrp.vehicles.each{ |v| v.cost_fixed = 5 }
@@ -135,7 +136,7 @@ def test_minimum_stop_in_route
       assert result[:routes].all?{ |r| (r[:activities].size - 2).zero? || r[:activities].size - 2 >= 5 },
              'Expecting no route with less than 5 stops unless it is an empty route'
       assert_operator should_remain_assigned, :<=, (result[:routes].sum{ |r| r[:activities].size - 2 })
-      assert_equal 25, result[:unassigned].size
+      assert_operator result[:unassigned].size, :>=, should_remain_unassigned
 
       all_ids = (result[:routes].flat_map{ |route| route[:activities].collect{ |stop| stop[:service_id] } }.compact +
                 result[:unassigned].collect{ |un| un[:service_id] }).uniq
diff --git a/test/wrapper_test.rb b/test/wrapper_test.rb
index 6a57ae4af..647290b73 100644
--- a/test/wrapper_test.rb
+++ b/test/wrapper_test.rb
@@ -1043,9 +1043,9 @@ def test_input_zones
         distance: [[0, 376184], [379177, 0]]
       }],
       points: [{
-        id: 'point_0', location: { lat: 48, lon: 5 }
+        id: 'point_0', location: { lat: 48, lon: 5 }, matrix_index: 0
       }, {
-        id: 'point_1', location: { lat: 49, lon: 1 }
+        id: 'point_1', location: { lat: 49, lon: 1 }, matrix_index: 1
       }],
       zones: [{
         id: 'zone_0',