diff --git a/src/dkl_moran.c b/src/dkl_moran.c index 33b8b95..db6e7e3 100644 --- a/src/dkl_moran.c +++ b/src/dkl_moran.c @@ -76,15 +76,15 @@ csr_sparse_matrix *moran_matrix_csr(wf_parameters *wf) { moran_row(row, 1, Ne, u, v, w_AA, w_Aa, w_aa); row[2] = 1 - row[1] - row[3] - row[4]; - Q->data[0] = row[2]; - Q->data[1] = row[3]; - Q->data[2] = row[4]; + Q->data[0] = 1.0 - row[2]; + Q->data[1] = -row[3]; + Q->data[2] = -row[4]; moran_row(row, 2, Ne, u, v, w_AA, w_Aa, w_aa); - Q->data[3] = row[1]; - Q->data[4] = row[2]; - Q->data[5] = row[3]; - Q->data[6] = row[4]; + Q->data[3] = -row[1]; + Q->data[4] = 1.0 - row[2]; + Q->data[5] = -row[3]; + Q->data[6] = -row[4]; // Tail // [. . * * * *] @@ -103,15 +103,15 @@ csr_sparse_matrix *moran_matrix_csr(wf_parameters *wf) { moran_row(row, Ne - 1, Ne, u, v, w_AA, w_Aa, w_aa); row[2] = 1 - row[0] - row[1] - row[3]; - Q->data[nnz - 3] = row[0]; - Q->data[nnz - 2] = row[1]; - Q->data[nnz - 1] = row[2]; + Q->data[nnz - 3] = -row[0]; + Q->data[nnz - 2] = -row[1]; + Q->data[nnz - 1] = 1.0 - row[2]; moran_row(row, Ne - 2, Ne, u, v, w_AA, w_Aa, w_aa); - Q->data[nnz - 7] = row[0]; - Q->data[nnz - 6] = row[1]; - Q->data[nnz - 5] = row[2]; - Q->data[nnz - 4] = row[3]; + Q->data[nnz - 7] = -row[0]; + Q->data[nnz - 6] = -row[1]; + Q->data[nnz - 5] = 1.0 - row[2]; + Q->data[nnz - 4] = -row[3]; #pragma omp parallel for for(DKL_INT i = 3; i < size - 1; i++) { @@ -126,11 +126,11 @@ csr_sparse_matrix *moran_matrix_csr(wf_parameters *wf) { moran_row(row, i, Ne, u, v, w_AA, w_Aa, w_aa); - Q->data[l ] = row[0]; - Q->data[l + 1] = row[1]; - Q->data[l + 2] = row[2]; - Q->data[l + 3] = row[3]; - Q->data[l + 4] = row[4]; + Q->data[l ] = -row[0]; + Q->data[l + 1] = -row[1]; + Q->data[l + 2] = 1.0 - row[2]; + Q->data[l + 3] = -row[3]; + Q->data[l + 4] = -row[4]; } @@ -149,16 +149,6 @@ csr_sparse_matrix *moran_matrix_csr(wf_parameters *wf) { // printf("%" PRId64 ",", Q->row_index[i]); // } - double *A = csr_to_dense(Q); - - for (DKL_INT i = 0; i < Q->nrows; i++) { - for (DKL_INT j = 0; j < Q->ncols; j++) { - DKL_INT idx = i * Q->ncols + j; - printf("%f, ", A[idx]); - } - printf("\n"); - } - printf("\n\n"); diff --git a/src/dkl_wf.c b/src/dkl_wf.c index 180b2c0..86604b1 100644 --- a/src/dkl_wf.c +++ b/src/dkl_wf.c @@ -322,6 +322,18 @@ void wf_solve(wf_parameters *wf, wf_statistics *r, double zero_threshold, bool m assert(csr_sparse_is_correct(A) == true); end_time = get_current_time(); printf("Building matrix: %gs\n", end_time - start_time); + + double *T = csr_to_dense(A); + + for (DKL_INT i = 0; i < A->nrows; i++) { + for (DKL_INT j = 0; j < A->ncols; j++) { + DKL_INT idx = i * A->ncols + j; + printf("%f, ", T[idx]); + } + printf("\n"); + } + dkl_dealloc(T); + #endif if (!moran) {