Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

## Rank-1 update to LU matrices Beginner
296 Views

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

1 Solution Moderator
296 Views

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.

5 Replies Moderator
297 Views

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. Beginner
296 Views

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(:)
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. Moderator
296 Views

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

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) Beginner
296 Views

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

Thanks again. 