- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

Link Copied

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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.

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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)

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content

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

Thanks again.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page