free'ing after allocation seems to crash my code in interplay with intel's mkl_sparse library

69 Views Asked by At

I am currently updating a piece of Code that eventually should be used inside a loop for time evolution. I have encountered some problems thatseem to come down to the notorious pointer allocation and free'ing thing.

I have an algorithm that creates large matrices and stores them in binary files, and another function, load_sparse_matrix(), that reads the binary files and returns a sparse matrix of the sparse_matrix_t format. Note that in order to prevent stack overflow, I allocate 3 diffrent heap pointers that temporarily store the matrix values. In the minimal example, I've commented the fread-part out and instead created an identy matrix in lines 55 resp 67. The sparse matrix is succesfully created and returned, and obviously, before I return, I'd like to free the 3 pointer arrays to prevent a memory leak (lines 96 -97).

I now want to perform a trivial matrix-vector-multiplication ( Identity Matrix * 0-Vector, line 30), but this operation causes Segmentation Fault. When, in turn, I don't free the pointers from the load_sparse_matrix()-function, the matrix-vector-multiplication is performed correclty, at the cost of a memory leak.

What am I missing? It's been a few days and I can't wrap my head around it. Maybe someone with a more intricate understanding of pointers and the sparse matrix library can find the error at eye sight.

Please find the not-working minimal example below, and feel free to try it on the Online-IDE replit, where I have added to relevant libraries and makefiles: Code

#include <complex.h>
#include <math.h>
#include <mkl_types.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MKL_Complex16 complex
#include "mkl_spblas.h"

#define L 4
#define N 16

// void loadCRS_int(char *desc, char *rowcolumn, int num_values, MKL_INT *read);
// void loadCRS_float(char *desc, int num_values, double complex *read);
sparse_matrix_t load_sparse_matrix(char *observable_name);

int main() {
  char *observable_name = "s^z_4";

  sparse_matrix_t A = load_sparse_matrix(observable_name);
  struct matrix_descr descr_matrix = {SPARSE_MATRIX_TYPE_HERMITIAN,
                                      SPARSE_FILL_MODE_UPPER,
                                      SPARSE_DIAG_NON_UNIT};

  double complex *A_init_state = calloc(N, sizeof(double complex));
  double complex *phi = calloc(N, sizeof(double complex));
  double complex alpha = 1.0 + I * 0.0, beta = 0.0 + I * 0.0;

  mkl_sparse_z_mv(SPARSE_OPERATION_NON_TRANSPOSE, alpha, A, descr_matrix, phi,
                  beta, A_init_state);

  free(A_init_state);
  free(phi);
  mkl_sparse_destroy(A);
  return 0;
}

sparse_matrix_t load_sparse_matrix(char *observable_name) {
  MKL_INT *observable_rowPtr = calloc(N + 1, sizeof(MKL_INT));
  // char matrix_CRS_buffer[512];
  // sprintf(matrix_CRS_buffer, "CRS_files/matrix_CRS_L=%02d_%s_row.dat", L,
  // observable_name); // builds variable dependend file names for observable

  // FILE *file = fopen(matrix_CRS_buffer, "r");
  // if (file == NULL) {
  //     printf("file does not exist\n");
  //     printf("filename %s\n", matrix_CRS_buffer);
  // }
  // fread(observable_rowPtr, sizeof(MKL_INT), N+1, file);
  // fclose(file);

  // loadCRS_int(observable_name, "row", N + 1 , observable_rowPtr);

  for (int i = 0; i <= N; i++) {
    observable_rowPtr[i] = i;
  }

  int num_values = observable_rowPtr[N];

  double complex *observable_val = calloc(num_values, sizeof(double complex));
  MKL_INT *observable_colInd = calloc(num_values, sizeof(MKL_INT));

  // loadCRS_int(observable_name, "column", num_values , observable_colInd);
  // loadCRS_float(observable_name, num_values, observable_val);

  for (int i = 0; i < num_values; i++) {
    observable_colInd[i] = i;
    observable_val[i] = 1.0 * i;
  }

  // init sparse matrix
  sparse_matrix_t A;
  sparse_status_t status;
  int exit_status = 0;

  status = mkl_sparse_z_create_csr(
      &A, SPARSE_INDEX_BASE_ZERO, N, N, observable_rowPtr,
      observable_rowPtr + 1, observable_colInd,
      observable_val); // Create handle with matrix stored in CSR format
  if (status != SPARSE_STATUS_SUCCESS) {
    printf(" Error in mkl_sparse_d_create_csr: %d \n", status);
    exit_status = 1;
    return NULL;
  }

  status = mkl_sparse_optimize(A); // Analyze sparse matrix; choose proper
                                   // kernels and workload balancing strategy
  if (status != SPARSE_STATUS_SUCCESS) {
    printf(" Error in mkl_sparse_optimize: %d \n", status);
    exit_status = 1;
    return NULL;
  }

  free(observable_val);
  free(observable_colInd);
  free(observable_rowPtr);
  return A;
}

// void loadCRS_int(char *desc, char *rowcolumn, int num_values, MKL_INT *read)
// {
//   char matrix_CRS_buffer[512];
//   sprintf(matrix_CRS_buffer, "CRS_files/matrix_CRS_L=%02d_%s_%s.dat", L,
//   desc,
//           rowcolumn); // builds variable dependend file names for observable

//   FILE *file = fopen(matrix_CRS_buffer, "r");
//   if (file == NULL) {
//     printf("file does not exist\n");
//     printf("filename %s\n", matrix_CRS_buffer);
//     return;
//   }
//   fread(read, sizeof(int), num_values, file);
//   fclose(file);

//   return;
// }

// void loadCRS_float(char *desc, int num_values, double complex *read) {
//   char matrix_CRS_buffer[512];
//   sprintf(matrix_CRS_buffer, "CRS_files/matrix_CRS_L=%02d_%s_row.dat", L,
//           desc); // builds variable dependend file names for observable

//   FILE *file = fopen(matrix_CRS_buffer, "r");
//   if (file == NULL) {
//     printf("file does not exist\n");
//     printf("filename %s\n", matrix_CRS_buffer);
//     return;
//   }
//   fread(read, sizeof(double complex), num_values, file);
//   fclose(file);

//   return;
// }
0

There are 0 best solutions below