Added sparse matrix debug output

dekoning-lab · Mar 3, 2017 · d63f9ed · d63f9ed
1 parent 0f471c7
commit d63f9ed
Show file tree

Hide file tree

Showing 4 changed files with 44 additions and 13 deletions.
diff --git a/include/dkl_sparse_matrix.h b/include/dkl_sparse_matrix.h
@@ -28,4 +28,6 @@ typedef struct csr_sparse_matrix {
  */
 bool csr_sparse_is_correct(csr_sparse_matrix *A);
 
+double *csr_to_dense(csr_sparse_matrix *A);
+
 #endif /* DKL_SPARSE_MATRIX_H */
diff --git a/src/dkl_moran.c b/src/dkl_moran.c
@@ -135,19 +135,30 @@ csr_sparse_matrix *moran_matrix_csr(wf_parameters *wf) {
   }
 
   #ifdef DEBUG
-  for(DKL_INT i = 0; i < nnz; i++) {
-    printf("%f,", Q->data[i]);
+  // for(DKL_INT i = 0; i < nnz; i++) {
+  //   printf("%f,", Q->data[i]);
+  // }
+  // printf("\n\n");
+  //
+  // for(DKL_INT i = 0; i < nnz; i++) {
+  //   printf("%" PRId64 ",", Q->cols[i]);
+  // }
+  // printf("\n\n");
+  //
+  // for(DKL_INT i = 0; i < Q->nrows + 1; i++) {
+  //   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");
 
-  for(DKL_INT i = 0; i < nnz; i++) {
-    printf("%" PRId64 ",", Q->cols[i]);
-  }
-  printf("\n\n");
-
-  for(DKL_INT i = 0; i < Q->nrows + 1; i++) {
-    printf("%" PRId64 ",", Q->row_index[i]);
-  }
   printf("\n\n");
 
 

diff --git a/src/dkl_sparse_matrix.c b/src/dkl_sparse_matrix.c
@@ -22,3 +22,23 @@ bool csr_sparse_is_correct(csr_sparse_matrix *A) {
   }
   return true;
 }
+
+double *csr_to_dense(csr_sparse_matrix *A) {
+  // Copy the buffer in the CSR representation into Q
+  // https://software.intel.com/en-us/node/520848
+
+  double *buffer = dkl_alloc(A->nrows * A->ncols, double);
+  DKL_INT info;
+  DKL_INT *job = dkl_alloc(6, DKL_INT);
+  job[0] = 1; // the rectangular matrix A is restored from the CSR format
+  job[1] = 0; // zero-based indexing for the rectangular matrix A is used
+  job[2] = 0; // zero-based indexing for the matrix in CSR format is used
+  job[3] = 2; // `buffer` is a whole matrix A
+  job[4] = A->nnz; // maximum number of the non-zero elements allowed if job[0]=0
+
+  mkl_ddnscsr(job, &A->nrows, &A->ncols, buffer, &A->ncols, A->data, A->cols, A->row_index, &info);
+  assert(info == 0);
+  dkl_dealloc(job);
+
+  return buffer;
+}
diff --git a/src/dkl_wf.c b/src/dkl_wf.c
@@ -379,8 +379,6 @@ void wf_solve(wf_parameters *wf, wf_statistics *r, double zero_threshold, bool m
   pardiso_control[MKL_PARDISO_MATRIX_CHECK_OPTION] = MKL_PARDISO_MATRIX_CHECK_ENABLE;
   #endif
 
-  //}}}1
-
   // Symbolic Factorization
   pardiso_phase = MKL_PARDISO_SOLVER_PHASE_ANALYSIS;
   #ifdef DEBUG