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