 
             parallel<system>::for_each(thrust::make_counting_iterator(0u),
                                        thrust::make_counting_iterator(0u) + baq.matrix_index.size() / stride,
-                                       stride_hmm_matrix_index<system>(baq.matrix_index));
+                                       stride_hmm_index<system>(baq.matrix_index));
+
+            // compute the scaling index
+            parallel<system>::for_each(thrust::make_counting_iterator(0u),
+                                       thrust::make_counting_iterator(0u) + num_matrices,
+                                       compute_hmm_scaling_factor_size_strided<system>(context, batch.device, active_baq_read_list, temp_active_list));
+
+            parallel<system>::inclusive_scan(temp_active_list.begin(),
+                                             num_matrices,
+                                             baq.scaling_index.begin() + 1,
+                                             thrust::plus<uint32>());
+
+            parallel<system>::for_each(thrust::make_counting_iterator(0u),
+                                       thrust::make_counting_iterator(0u) + baq.scaling_index.size() / stride,
+                                       stride_hmm_index<system>(baq.scaling_index));
         }
 
         uint32 matrix_len = baq.matrix_index[baq.matrix_index.size() - 1];
+        uint32 scaling_len = baq.scaling_index[baq.scaling_index.size() - 1];
 
         baq.forward.resize(matrix_len);
         baq.backward.resize(matrix_len);
-
-        // compute the index and size of the HMM scaling factors
-        baq.scaling_index.resize(num_active + 1);
-        // first offset is zero
-        thrust::fill_n(baq.scaling_index.begin(), 1, 0);
-        parallel<system>::inclusive_scan(thrust::make_transform_iterator(active_baq_read_list.begin(),
-                                                                         compute_hmm_scaling_factor_size<system>(context, batch.device)),
-                                         num_active,
-                                         baq.scaling_index.begin() + 1,
-                                         thrust::plus<uint32>());
-
-        uint32 scaling_len = baq.scaling_index[num_active];
         baq.scaling.resize(scaling_len);
 
 //        fprintf(stderr, "reads: %u\n", batch.num_reads);

            }
        }

//        printf("matrix index %u size %u\n", index, ret);
        matrix_size_output[index] = ret;
    }
};

template <target_system system>
struct stride_hmm_matrix_index
{
    typename vector<system, uint32>::view matrix_index;

    stride_hmm_matrix_index(typename vector<system, uint32>::view matrix_index)
        : matrix_index(matrix_index)
    { }

    CUDA_HOST_DEVICE void operator() (const uint32 matrix_group_index)
    {
        constexpr uint32 stride = baq_stride<system>::stride;

        auto initial_element = matrix_index[matrix_group_index * stride];
        for(uint32 i = 1; i < stride; i++)
        {
            matrix_index[matrix_group_index * stride + i] = initial_element + i;
        }
    }
};

        const uint32 matrix_offset = matrix_index % baq_stride<system>::stride;
        const uint32 matrix_base = matrix_index - matrix_offset;

        const uint32 scaling_offset = scaling_index % baq_stride<system>::stride;
        const uint32 scaling_base = scaling_index - scaling_offset;

        // set up matrix and scaling factor pointers
        forwardMatrix = matrix_iterator(&ctx.baq.forward[matrix_base + matrix_offset]);
        backwardMatrix = matrix_iterator(&ctx.baq.backward[matrix_base + matrix_offset]);
        scalingFactors = matrix_iterator(&ctx.baq.scaling[scaling_base + scaling_offset]);

        // get the windows for the current read
        const auto& hmm_reference_window = ctx.baq.hmm_reference_windows[read_index];
