Fortran Interface for Finite Element nonlinear solver

[AndreaBSItaly]

Dear All,
I've written a FE Fortran code (for academic use only) that, by now, use PANUA Paradiso. I would replace the Pardiso interface with STRUMPACK, because it seems to me that PANUA Paradiso is becoming a commercial package. If I've correctly understood, both PARDISO and STRUMPACK should use the same matrix format for sparse matrices.

Basically, I've two type of matrices:

• Real and symmetric indefinite, diagonal or Bunch-Kaufman pivoting (MTYPE=-2 with PARDISO). (In this case I've to store only the upper part of the matrix)
• Real and non-symmetric complete supernode pivoting (MTYPE=11 with PARDISO).

I've some troubles with the fortran interface. In particular, I'm getting some memory leak. Do you know if the format for the matrix is exactly the same?
With Pardiso I simply set mtype (-2 or 11), the "phase" parameter (symbolic and numeric factorization, solve, iterative refinement)
and then I call the external subroutine as
CALL pardiso (pt, maxfct, mnum, mtype, phase,NGDL,KMATR,Kia,
1 Kja, idum, nrhs, iparm, msglvl, F, u, error)

With STRUMPACK I'm doing:

 REAL*8, DIMENSION(:), POINTER, INTENT(INOUT) :: KMatr
  INTEGER, DIMENSION(:), POINTER, INTENT(INOUT) :: KJA 
  INTEGER, DIMENSION(:), POINTER, INTENT(INOUT) :: KIA 
  REAL*8, DIMENSION(:), TARGET, INTENT(INOUT) :: F, U
  LOGICAL, INTENT(IN) :: Symm
  LOGICAL, INTENT(INOUT) :: FLAGNEGEIG !Used to write output of
                                       !negative eigenvalue

  INTEGER(c_int) :: error 
  type(STRUMPACK_SparseSolver) :: S
  INTEGER, target :: NGDL
  INTEGER :: N


  NGDL=SIZE(KIA)-1 !SIZE OF THE MATRIX THAT I WOULD LIKE TO SOLVE
  N=SIZE(KMatr)  !NUMBER OF NONZERO ITEMS


  call STRUMPACK_init_mt(S, STRUMPACK_DOUBLE, STRUMPACK_MT, 
 #          0, c_null_ptr, 1)   
  call STRUMPACK_set_matching(S, STRUMPACK_MATCHING_NONE);
  call STRUMPACK_set_reordering_method(S, STRUMPACK_METIS)

c Set compression method. Other options include NONE, HSS, HODLR,
c LOSSY, LOSSLESS. HODLR is only supported in parallel, and only
c supports double precision (including complex double).
call STRUMPACK_set_compression(S, STRUMPACK_BLR);

c Set the block size and relative compression tolerances for BLR
c compression.
call STRUMPACK_set_compression_leaf_size(S, 64);
call STRUMPACK_set_compression_rel_tol(S, dble(1.e-2));

c Only sub-blocks in the sparse triangular factors corresponing to
c separators larger than this minimum separator size will be
c compressed. For performance, this value should probably be larger
c than 128. This value should be larger for HODLR/HODBF, than for
c BLR, since HODLR/HODBF have larger constants in the complexity.
c For an n x n 2D domain, the largest separator will correspond to
c an n x n sub-block in the sparse factors.
call STRUMPACK_set_compression_min_sep_size(S, 300);

  call STRUMPACK_set_csr_matrix
 #  (S, c_loc(NGDL), c_loc(KIa), c_loc(KJa), c_loc(KMatr), 1);
  PRINT *, 'ORDINAMENTO'
  error = STRUMPACK_reorder(S)

  PRINT *, 'SOLUTORE'
  error = STRUMPACK_solve(S, c_loc(F), c_loc(u), 0);

Can anyone help me? Thank you in advance!

Andrea

[AndreaBSItaly]

I think I resolved my issue. The format is not the same (pardiso use FORTRAN indices from 1, while STRUMPACK uses indices from 0), and it seems that STRUMPACK needs the all the nonzero terms (even in the case of a symmetric matrix)
Only a few simple questions. I'm not an expert of linear system solvers.

1. Do you have any suggestion to tune or speedup the solver?
   Now I do:
   call STRUMPACK_init_mt(S, STRUMPACK_DOUBLE, STRUMPACK_MT,

   0, c_null_ptr, 1)

   call STRUMPACK_set_reordering_method(S, STRUMPACK_METIS)
   call STRUMPACK_set_compression(S, STRUMPACK_BLR);
   call STRUMPACK_set_compression_leaf_size(S, 64);
   call STRUMPACK_set_compression_rel_tol(S, dble(1.e-2));
   call STRUMPACK_set_compression_min_sep_size(S, 300);
   call STRUMPACK_set_csr_matrix

   (S, c_loc(NGDL), c_loc(Ia), c_loc(Ja), c_loc(KMatr), 1);

   error = STRUMPACK_solve(S, c_loc(F), c_loc(u), 1);

2. Is it possible to avoid the output (without modifying your source code?)
   Thank you again!
   Andrea

[pghysels]

To disable the output, set the last argument for

  call STRUMPACK_init_mt(S, STRUMPACK_DOUBLE, STRUMPACK_MT, 
 #          0, c_null_ptr, 1)   

to 0.
Or call

call STRUMPACK_set_verbose(S, 0)

[pghysels]

To speed up the solver, you can either try to use the GPU (for now this only works without compression, ie, STRUMPACK_set_compression(S, STRUMPACK_NONE);)
or you need to tune these parameters:

call STRUMPACK_set_compression_leaf_size(S, 64);
call STRUMPACK_set_compression_rel_tol(S, dble(1.e-2));
call STRUMPACK_set_compression_min_sep_size(S, 300);

Setting a larger min_sep_size (say 1000) and a larger leaf_size (fi 256) will probably lead to better performance, but at the cost of higher memory usage.
The compression_rel_tol is already relatively large, but you could try
call STRUMPACK_set_compression_rel_tol(S, dble(1.e-1));, although this will lead to slower convergence in the iterative solver.

[AndreaBSItaly]

Dear Pieter, thank you very much for your help Andrea

…

Il giorno 20 lug 2023, alle ore 09:28, Pieter Ghysels ***@***.***> ha scritto: the

-- Informativa sulla Privacy: https://www.unibs.it/it/node/1452 <https://www.unibs.it/it/node/1452>