Rank-1 update to LU matrices

[Sunny_S_1]

Hi,

I have a matrix Y for which I need LU factors that I can get using DSS interface Routines like dss_factor_complex. However, the diagonal entries of this matrix need to be updated after each iteration of a loop that gets me LU factors.

Since Y is very large, and LU factorization will take a lot of computational effort. Is there a way to update LU factors when Y updates without having to do LU factorization on every iteration?

I am trying to do this:

iteration1: Y (original) > reorder > factor (LU)> solve

iteration2: Y' (updated) > update LU (from iteration1) > solve

Iteration3: Y"(updated) > update LU (from iteration2) > solve

Thanks.

Sunny

[Gennady_F_Intel · Accepted solution]

sunny, that's not possible with DSS API. Please have a look at the PARDISO API - iparm(4) - Preconditioned CGS/CG.  Based on your explanation above, it should helps.  

[Sunny_S_1]

Thanks Gennady for your answer.

I tried your recommendation with my code, but ran into another problem. I am trying to follow an example in the compiler folder named "pardiso_unsym_complexf.f" very closely. However, when I call pardiso, I am getting an ERROR = -1 suggesting an inconsistent input. But I checked that my arguments were consistent with PARDISO_DC.

Note that I have a complex matrix (A) that I am trying to factorise. To help you give a context, I am providing an excerpt from my code:

Declaration:

    !DOUBLE PRECISION
    INTEGER, PARAMETER :: DP = 8
    !INTERNAL SOLVER MEMORY POINTER
    TYPE(MKL_PARDISO_HANDLE), ALLOCATABLE :: PT(:)
    !OTHER PARADISO VARIABLES
    INTEGER MAXFCT, MNUM, MTYPE, PHASE, N, NRHS, ERROR, MSGLVL, NNZ     
    INTEGER ERROR1
    INTEGER, ALLOCATABLE :: IPARM(:)
    INTEGER, ALLOCATABLE :: IA(:)
    INTEGER, ALLOCATABLE :: JA(:)
    COMPLEX(KIND=DP), ALLOCATABLE :: A( : )
    COMPLEX(KIND=DP), ALLOCATABLE :: B( : )
    COMPLEX(KIND=DP), ALLOCATABLE :: X( : )
    INTEGER II, JJ, NI, IIA, KEY, IDUM(1), IBUS, ISQBUS, KSQBUS, SIZEYDIAG, SIZEYOFF
    COMPLEX(KIND=DP) :: DDUM(1)

Parameters:

    IPARM(1) = 1 ! no solver default
    IPARM(2) = 2 ! fill-in reordering from METIS
    iparm(3) = 1 ! numbers of processors, value of OMP_NUM_THREADS
    IF (IMODE == 0) THEN
        IPARM(4) = 60 ! no iterative-direct algorithm - first pass 0 (60); second onwards 2 (62)
    ELSE
        IPARM(4) = 62
    ENDIF
    IPARM(5) = 0 ! no user fill-in reducing permutation
    IPARM(6) = 0 ! =0 solution on the first n compoments of x
    IPARM(8) = 9 ! numbers of iterative refinement steps
    IPARM(10) = 13 ! perturbe the pivot elements with 1E-13
    IPARM(11) = 1 ! use nonsymmetric permutation and scaling MPS
!    IPARM(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
    iparm(13) = 1 ! maximum weighted matching algorithm is switched-on (default for non-symmetric).
    IPARM(14) = 0 ! Output: number of perturbed pivots
    IPARM(18) = -1 ! Output: number of nonzeros in the factor LU
    IPARM(19) = -1 ! Output: Mflops for LU factorization
    IPARM(20) = -1 ! Output: Numbers of CG Iterations (-1)
    
    ERROR = 0 ! INTIALISE ERROR FLAG
    MSGLVL = 1 !PRINT STATISTICAL INFORMATION
!    MTYPE = 6 !COMPLEX, SYMMETRIC
!    MTYPE = 3 !COMPLEX, STRUCTURALLY SYMMETRIC
    mtype     = 13      ! complex unsymmetric matrix

Finally, the call:

    !INITIALIZE THE INTERNAL SOLVER MEMORY POINTER - ONLY NECESSARY FOR THE FIRST CALL OF THE PARDISO SOLVER
    ALLOCATE(PT(64))
    DO II = 1, 64
        PT(II)%DUMMY = 0
    END DO
    
    !REORDER AND SYMBOLIC FACTORISATION
    PHASE = 11
    
    CALL PARDISO(PT, MAXFCT, MNUM, MTYPE, PHASE, N, A, IA, JA, IDUM, NRHS, IPARM, MSGLVL, DDUM, DDUM, ERROR)

Thanks.

[Gennady_F_Intel]

please add matrix checker - iparm(27) = 1 and see the message you will get

[Gennady_F_Intel]

I have also check how works the original test you mentioned - pardiso_unsym_complex_f

It works fine on my side. here is the output I see:

pardiso_unsym_complex_f.exe

=== PARDISO: solving a complex nonsymetric system ===
1-based array indexing is turned ON
PARDISO double precision computation is turned ON
METIS algorithm at reorder step is turned ON
Scaling is turned ON
Matching is turned ON

Summary: ( reordering phase )
================

Times:
======
Time spent in calculations of symmetric matrix portrait (fulladj): 0.000012 s
Time spent in reordering of the initial matrix (reorder)         : 0.000144 s
Time spent in symbolic factorization (symbfct)                   : 0.000301 s
Time spent in data preparations for factorization (parlist)      : 0.000001 s
Time spent in allocation of internal data structures (malloc)    : 0.007281 s
Time spent in additional calculations                            : 0.000039 s
Total time spent                                                 : 0.007778 s

.....................

.....................

 Solve completed ...

 The solution of the system is
  x(           1 ) =  (0.174767984570877,2.117724204435874E-002)
  x(           2 ) =  (-0.176470588235294,-0.294117647058824)
  x(           3 ) =  (4.932176149148183E-002,2.959819993571200E-002)
  x(           4 ) =  (4.298112687698030E-002,-3.140928502548559E-002)
  x(           5 ) =  (-0.120858887817422,-0.170859530697525)
  x(           6 ) =  (-0.369347476695596,-8.614593378335146E-004)
  x(           7 ) =  (9.161041465766635E-002,0.125361620057859)
  x(           8 ) =  (0.223940855030537,0.139427836708454)

[Sunny_S_1]

Thanks Gennady. My ja was not ordered in ascending order. I fixed it and now my code works.

Thanks again.