Memory crash in DSS

[Max_T_]

Hi!

I modified DSS solver so that I could solve my own system using it but at the dss_solve_real step memory error occurs. I use malloc to allocate memory for the system while the example dss_sym_c.c uses static constants.

When I try to solve 1000 x 1000 system my program crashes with a segmentation fault. In case of larger systems 100 000 x 100 000 "MKL-DSS-DSS-Error, Out of memory" happens.

Here's the code:

#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include "mkl_dss.h"
#include "mkl_types.h"

MKL_INT
main ()
{
  MKL_INT i;
  _MKL_DSS_HANDLE_t handle;
  _INTEGER_t error;
  _CHARACTER_t statIn[] = "determinant";
  _DOUBLE_PRECISION_t statOut[5];
  MKL_INT opt = MKL_DSS_DEFAULTS;
  MKL_INT sym = MKL_DSS_SYMMETRIC;
  MKL_INT type = MKL_DSS_INDEFINITE;
  MKL_INT nRows;
  MKL_INT nCols;
  MKL_INT nNonZeros;
  MKL_INT nRhs;
  _INTEGER_t *rowIndex;
  _INTEGER_t *columns;
  _DOUBLE_PRECISION_t *values;
  _DOUBLE_PRECISION_t *rhs;
  _DOUBLE_PRECISION_t *solValues;

  FILE *In;
  MKL_INT n, nz;
  MKL_INT *I = NULL, *J = NULL;
 _DOUBLE_PRECISION_t *val = NULL, *sol = NULL, *RHS = NULL;

  if ( ( In = fopen ("matr_rhs_binary", "rb" ) ) == NULL ) {
    fprintf (stderr, "Error opening file\n");		
    return (-1);
  }

  fread ( &n, sizeof(MKL_INT), 1, In );
  printf ("n = %d\n", n);	

  fread ( &nz, sizeof(MKL_INT), 1, In );
  printf ("nz = %d\n", nz);

  I = (MKL_INT *) malloc ( sizeof (I[0]) * (n + 1));                             // csr row pointers for matrix A
  J = (MKL_INT *) malloc ( sizeof (J[0]) * nz);                                  // csr column indices for matrix A
  val = (_DOUBLE_PRECISION_t *) malloc ( sizeof (val[0]) * nz);                  // csr values for matrix A	
  sol = (_DOUBLE_PRECISION_t *) malloc ( sizeof (sol[0]) * n);
  RHS = (_DOUBLE_PRECISION_t *) malloc ( sizeof (RHS[0]) * n);

  for (i = 0; i < n + 1; i++)
    fread ( &(I), sizeof(I[0]), 1, In );
  for (i = 0; i < nz; i++)
    fread ( &(val), sizeof(val[0]), 1, In );
  for (i = 0; i < nz; i++)
    fread ( &(J), sizeof(J[0]), 1, In );
  for (i = 0; i < n; i++)
    fread ( &(RHS), sizeof(RHS[0]), 1, In);	
	
  fclose (In);

  nRows = n;
  nCols = n;
  nNonZeros = nz;
  nRhs = n;

  rowIndex = (_INTEGER_t *) malloc ( sizeof (rowIndex[0]) * (nRows + 1) );
  columns = (_INTEGER_t *) malloc ( sizeof (columns[0]) * nNonZeros );
  values = (_DOUBLE_PRECISION_t *) malloc ( sizeof (values[0]) * nNonZeros );
  rhs = (_DOUBLE_PRECISION_t *) malloc ( sizeof (rhs[0]) * nCols );
  solValues = (_DOUBLE_PRECISION_t *) malloc ( sizeof (solValues[0]) * nCols );

  for (i = 0; i < n + 1; i++)
    rowIndex = I;
  for (i = 0; i < n; i++)
    rhs = RHS;
  for (i = 0; i < nz; i++) {
    columns = J;
    values = val;
  }
/* --------------------- */
/* Initialize the solver */
/* --------------------- */
  error = dss_create (handle, opt);
  if (error != MKL_DSS_SUCCESS)
    goto printError;
/* ------------------------------------------- */
/* Define the non-zero structure of the matrix */
/* ------------------------------------------- */
  error = dss_define_structure (handle, sym, rowIndex, nRows, nCols,
				columns, nNonZeros);
  if (error != MKL_DSS_SUCCESS)
    goto printError;
/* ------------------ */
/* Reorder the matrix */
/* ------------------ */
  error = dss_reorder (handle, opt, 0);
  if (error != MKL_DSS_SUCCESS)
    goto printError;
/* ------------------ */
/* Factor the matrix */
/* ------------------ */
  error = dss_factor_real (handle, type, values);
  if (error != MKL_DSS_SUCCESS)
    goto printError; 
/* ------------------------ */
/* Get the solution vector */
/* ------------------------ */
  error = dss_solve_real (handle, opt, rhs, nRhs, solValues);
  if (error != MKL_DSS_SUCCESS)
    goto printError;
/* ------------------------ */
/* Get the determinant (not for a diagonal matrix) */
/*--------------------------*/
  if (nRows < nNonZeros)
    {
      error = dss_statistics (handle, opt, statIn, statOut);
      if (error != MKL_DSS_SUCCESS)
	goto printError;
/*-------------------------*/
/* print determinant */
/*-------------------------*/
      printf (" determinant power is %g \n", statOut[0]);
      printf (" determinant base is %g \n", statOut[1]);
      printf (" Determinant is %g \n", (pow (10.0, statOut[0])) * statOut[1]);
    }
/* -------------------------- */
/* Deallocate solver storage */
/* -------------------------- */
  error = dss_delete (handle, opt);
  if (error != MKL_DSS_SUCCESS)
    goto printError;
/* ---------------------- */
/* Print solution vector */
/* ---------------------- */
  printf (" Solution array: ");
  for (i = 0; i < nCols; i++)
    printf (" %g", solValues);
  printf ("\n");
  exit (0);
printError:
  printf ("Solver returned error code %d\n", error);
  exit (1);
}

Thanks!

[Chao_Y_Intel]

Hello, 

For the large matrix size with the out of memory error, you can use the out of core interfaces,  it will save some computation data into the disk. 
For the segmentation fault with the matrix size 1000x1000, is that only matrix that report the problem? For other similar or small size of Matrix,  can your code work for you? 
If this is only related to that matrix, we need the test data to investigate the problem.   If you can not attached in the user forum, you may submit some issue in the intel premier support  to attach the data. 

Thanks,
Chao 

[Max_T_]

Oops! I found a misprint: nRhs = 1. Now memory doesn't crash, but I get a zero solution. The problem is definately not in the matrix because I tested the program on several systems and was always getting zero.  What can be the reason of it?