topic Rank-1 update to LU matrices in Intel® oneAPI Math Kernel Library

Rank-1 update to LU matrices

Sunny_S_1 — Thu, 19 Nov 2015 01:41:47 GMT

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

sunny, that's not possible

Gennady_F_Intel — Thu, 19 Nov 2015 12:10:35 GMT

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.

Thanks Gennady for your

Sunny_S_1 — Tue, 01 Dec 2015 01:16:41 GMT

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.

please add matrix checker -

Gennady_F_Intel — Tue, 01 Dec 2015 04:07:59 GMT

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

I have also check how works

Gennady_F_Intel — Tue, 01 Dec 2015 04:16:33 GMT

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)

Thanks Gennady. My ja was not

Sunny_S_1 — Tue, 08 Dec 2015 18:11:14 GMT

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

Thanks again.