Pardiso does not keep matrix factorization in memory (?) for subsequent "solve" computations

[Ioannis_K_]

Hello,

I am trying to call pardiso through a subroutine in a large program of mine.
The subroutine receives the coefficient array from the main program as a dummy argument Kffv.

This is the subroutine first line, together with an explanation of various arguments (some of them are not important for the issue that I am reporting, as they are not directly used with pardiso):

===========================================================================================

subroutine gauselim(Kffv,bvec,n,nrows,perm,rowIndex,Icol,Kdiag,Idiag,formstffness)

c     INPUT ARGUMENTS: 
c     ----------------                
c      bvec:    RHS vector
c      formstffness: logical variable (see my explanation below
c      Icol:    integer vector with dimension "n": Icol(iv) = the column number to which component "iv" of Kffv corresponds
c      Idiag:   integer vector with dimension "nrows", containing the locations (in Kffv) of the diagonal components of the coefficient array
c      Kdiag:   real vector with dimension nrows, containing the diagonal components of the coefficient array (used for preconditioning)
c      Kffv:    vector containing nonzero values of Kff        
c      n:       dimension of Kffv
c      nrows:   number of rows in system of equations   
c      perm:    integer vector, not used here (only passed for consistency for calling direct solver)
c      rowIndex: integer vector, containing the row information for coefficient array: rowIndex(i) = which component of Kffv corresponds to the 1st nonzero component of row "i", rowIndex(nrow+1) = (component of Kffv which corresponds to the last component of row nrow) - rowIndex(1)
c
c     OUTPUT ARGUMENTS:   
c     -----------------
c      bvec: solution vector      

===================================================================================

I need to call this subroutine for different values of bvec (which are calculated during program execution and are not all known simultaneously, so each time that I call the subroutine I must solve with a single RHS vector bvec. The coefficient array Kffv may or may not change. This is controlled by the logical argument formstffness: if formstffness = .true., then we have a new coefficient array and we need to factorize it before solving.
if formstffness = .false., then we have not changed the coefficient array from the previous factorization (and we do not need to factorize from scratch - we merely need to use the pre-existing factorization to directly solve).

I have the following lines of code in my program (after initializing the various arguments of pardiso):

=================================================================================================

       if(formstffness) then    ! check if coefficient array has changed, at which case we must factorize from scratch:
          ! first, clear memory:
          call pardiso (pt, maxfct, mnum, mtype, -1, nrows, Kffv,rowIndex, icol, perm, nrhs, iparm, msglvl, bvec, xvec, ierror)
           
        ! now, conduct numerical factorization:   
               phase = 12
      
      call pardiso (pt, maxfct, mnum, mtype, phase, nrows, Kffv,  rowIndex, icol, perm, nrhs, iparm, msglvl, bvec, xvec, ierror)
          endif

          ! finally, solve:
           phase = 33
      
      call pardiso (pt, maxfct, mnum, mtype, phase, nrows, Kffv,  rowIndex, icol, perm, nrhs, iparm, msglvl, bvec, xvec, ierror)

============================================================================

This code works fine the first time the routine is called (where formsffness = .true.)
However, the second time that the routine is called (with formsffness = .false.) it does not give a solution, which means that I have somehow lost the factorization.
Do you have any ideas/thoughts on what I may be doing wrong in my routine?
 

[mecej4]

Many of the arguments passed to Pardiso in the calls that you show are local variables whose declarations you did not show. If you wish to factorize once and solve many times, as you described, some of these variables (in particular, the PT array) may need to be given the SAVE attribute or have proper values set before calling Pardiso; without that, their values may be undefined.

[Kirill_V_Intel]

Hello Ioannis,

Please note (as also pointed by mecej4) that in order to re-use the factorization you need to store the PT array. From the code snippet which you give, it looks like you don't save it after the call to pardiso factorization phase. I'd suggest you add PT as a parameter to your subroutine gauselim(), save it after factorization and re-use it in the subsequent calls when you want to keep the same factorization and only solve the system with a new rhs.

Let us know if you still have problems after trying this.

Best,
Kirill

[Ioannis_K_]

Hello Kirill,

The problem was indeed due to the local nature of some arguments. I followed mecej4's suggestion, placed the various arrays/parameters of the factorization in a module with the SAVE attribute, and now everything works perfectly fine.

Many thanks for the help.