ralna / spral

Sparse Parallel Robust Algorithms Library
https://ralna.github.io/spral/
Other
105 stars 27 forks source link

Random matrices generated are not positive definite #215

Open chrhansk opened 3 months ago

chrhansk commented 3 months ago

I have been trying to generate random matrices according to the documentation using the following driver:

#include <stdio.h>
#include "spral.h"

int main(void) {
  int state = 0;
  int matrix_type = SPRAL_MATRIX_REAL_SYM_PSDEF;

  int m=3, n=3, nnz=3;
  int ptr[n+1], row[nnz];
  double val[nnz];

  /* Generate matrix */
  printf("Generating a %d x %d non-singular matrix with %d non-zeroes\n",
         m, n, nnz);

  if(spral_random_matrix_generate(&state, matrix_type, m, n, nnz, ptr, row, val,
                                  SPRAL_RANDOM_MATRIX_NONSINGULAR | SPRAL_RANDOM_MATRIX_SORT)) {
    printf("Error return from spral_random_matrix_generate()\n");
    return 1;
  }

  /* Print matrix using utility routine from SPRAL_MATRIX_UTILS package */
  printf("Generated matrix:\n");
  spral_print_matrix(-1, matrix_type, m, n, ptr, row, val, 0);

  /* Success */
  return 0;
}

The output of the program is the following:

Generating a 3 x 3 non-singular matrix with 3 non-zeroes
Generated matrix:
Real symmetric positive definite matrix, dimension 3x3 with 3 entries.
0:   9.9999E-01                          
1:               -3.1031E-01             
2:                             3.9037E-01

This matrix is symmetric but indefinite. It states in the documentation that In the positive-definite case, a post-processing step sums the absolute values of all the entries in each column and replaces the diagonal with this value. I don't see this happening in the code though. Is this functionality not implemented yet?

jfowkes commented 3 months ago

This matrix is symmetric but indefinite. It states in the documentation that In the positive-definite case, a post-processing step sums the absolute values of all the entries in each column and replaces the diagonal with this value. I don't see this happening in the code though. Is this functionality not implemented yet?

Most likely, the principal author of SPRAL left in a hurry and unfortunately not all functionality was implemented.