                                              num_active,
                                              baq.matrix_index.begin() + 1,
                                              thrust::plus<uint32>());
+
+            // 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>());
         } else {
             // round up the number of matrices to a multiple of baq_stride
             const uint32 num_matrices = ((num_active + stride - 1) / stride) * stride;
             temp_active_list.resize(num_matrices);
 
+            // allocate and initialize the matrix and scaling indices
             baq.matrix_index.resize(num_matrices + 1);
+            baq.scaling_index.resize(num_matrices + 1);
             // first offset is zero
             thrust::fill_n(baq.matrix_index.begin(), 1, 0);
+            thrust::fill_n(baq.scaling_index.begin(), 1, 0);
 
             // compute matrix sizes as the maximum size of each baq_stride groups of reads
             parallel<system>::for_each(thrust::make_counting_iterator(0u),


            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));
        }

        uint32 matrix_len = baq.matrix_index[baq.matrix_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);

        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];
