topic Thank you very much indeed. in Intel® oneAPI Math Kernel Library

INTEL MKL ERROR : PARAMETER 5 WAS INCORRECT ON ENTRY TO MKL_DDIASM

Abdolreza_A_ — Mon, 08 May 2017 18:32:01 GMT

Hello to all of my new firends.

I am writing a code for "Turbulent Negarively Bouyant Jet" in these days for my thesis in BSc.

I would like to use MKL_DDIASM for part of my code but when I RUN the code this message appeares and all of my RHS result (argument C) will be zero. The message is : INTEL MKL ERROR : PARAMETER 5 WAS INCORRECT ON ENTRY TO MKL_DDIASM

Then I tried to make a simple code to understand exactly what DDIASM do. I wrote 2 simple codes which I bring here. The first code works well but the secnod code doesn't. I notice that the second code have a general sparse matrix "A".

When I define MATDESCRA (1) = "G" then this message will apear : INTEL MKL ERROR : PARAMETER 5 WAS INCORRECT ON ENTRY TO MKL_DDIASM

my 1st answer : What's the problem with "G" for MATDESCRA (1) ????!!!

my 2nd answer : what should I do to solve the system of equations in my 2nd simple code to obtain correct answers for that with DDIASM? (How must I define MATDESCRA exactly?)

PROGRAM 2nd_simple_code
    IMPLICIT NONE
    CHARACTER(1) , DIMENSION ( 4 ) :: MATDESCRA
    REAL (8) , DIMENSION ( 3 , 3 ) :: A
    REAL (8) , DIMENSION ( 3 , 2 ) :: B , C
    REAL (8) , DIMENSION ( 3 , 3 ) :: VAL
    INTEGER I , J

    
    CHARACTER(1) TRANSA
    
    INTEGER :: M , N , LVAL , NDIAG , LDB , LDC
    REAL (8) :: ALPHA

    INTEGER , DIMENSION ( 3 ) :: IDIAG  
    
    
    A = 0
    B = 0
    C = 0
    VAL = 0
    IDIAG = 0
    
    A ( 1 , 1 ) = 1.0D0
    A ( 2 , 2 ) = 1.0D0
    A ( 3 , 3 ) = 1.0D0
    
    A ( 2 , 1 ) = 2.0D0
    A ( 1 , 3 ) = 5.0D0
    
    B ( 1 , 1 ) = 16.0D0
    B ( 1 , 2 ) = 34.0D0

    B ( 2 , 1 ) = 4.0D0
    B ( 2 , 2 ) = 13.0D0

    B ( 3 , 1 ) = 3.0D0
    B ( 3 , 2 ) = 6.0D0
    

    VAL ( 2 , 1 ) = 2.0D0
    VAL ( 1 , 2 ) = 1.0D0
    VAL ( 2 , 2 ) = 1.0D0
    VAL ( 3 , 2 ) = 1.0D0
    VAL ( 1 , 3 ) = 5.0D0
    
    
M = 3
N = 2
LVAL = 3
NDIAG = 3
LDB = 3
LDC = 3    
ALPHA = 1.0D0
TRANSA = "N"
MATDESCRA (1) = 'G'
MATDESCRA (2) = 'U'
MATDESCRA (3) = 'U'
MATDESCRA (4) = 'F'
    
IDIAG (1) = -1
IDIAG (2) = 0
IDIAG (3) = 2
    
    CALL MKL_DDIASM ( TRANSA , M , N , ALPHA , MATDESCRA , VAL , LVAL , IDIAG , NDIAG , B , LDB , C , LDC )
    

    DO I = 1 , 3
    WRITE ( * , * ) (C (I ,J), J = 1 , 2)
    END DO
READ(*,*)

END PROGRAM 2nd_simple_code

*****************************************************************************************************

PROGRAM 1st_simple_code
    IMPLICIT NONE
    CHARACTER(1) , DIMENSION ( 4 ) :: MATDESCRA
    REAL (8) , DIMENSION ( 3 , 3 ) :: A
    REAL (8) , DIMENSION ( 3 , 2 ) :: B , C
    REAL (8) , DIMENSION ( 3 , 1 ) :: VAL
    INTEGER I , J

    
    CHARACTER(1) TRANSA
    
    INTEGER :: M , N , LVAL , NDIAG , LDB , LDC
    REAL (8) :: ALPHA

    INTEGER , DIMENSION ( 1 ) :: IDIAG  
    
    
    A = 0
    B = 0
    C = 0
    
    A ( 1 , 1 ) = 1.0D0
    A ( 2 , 2 ) = 1.0D0
    A ( 3 , 3 ) = 1.0D0
    
    B ( 1 , 1 ) = 1.0D0
    B ( 1 , 2 ) = 4.0D0

    B ( 2 , 1 ) = 2.0D0
    B ( 2 , 2 ) = 5.0D0

    B ( 3 , 1 ) = 3.0D0
    B ( 3 , 2 ) = 6.0D0
    

    VAL ( 1 , 1 ) = 1.0D0
    VAL ( 2 , 1 ) = 1.0D0
    VAL ( 3 , 1 ) = 1.0D0

    
    
M = 3
N = 2
LVAL = 3
NDIAG = 1
LDB = 3
LDC = 3    
ALPHA = 1.0D0
TRANSA = "N"
MATDESCRA (1) = 'D'
MATDESCRA (2) = 'U'
MATDESCRA (3) = 'U'
MATDESCRA (4) = 'F'
    
IDIAG (1) = 0

CALL MKL_DDIASM ( TRANSA , M , N , ALPHA , MATDESCRA , VAL , LVAL , IDIAG , NDIAG , B , LDB , C , LDC )
    

    DO I = 1 , 3
    WRITE ( * , * ) (C (I ,J), J = 1 , 2)
    END DO
READ(*,*)

END PROGRAM 1st_simple_code

**************************************************************

Thanks a lot dear friends.

The reference manual page for

mecej4 — Tue, 09 May 2017 00:41:36 GMT

The reference manual page for the routine in question says:

A is a sparse upper or lower triangular matrix with unit or non-unit main diagonal

If you specify 'G' as the first element of MATDESCRA, you are in violation of this limitation, because G signifies "general", which includes non-triangular matrices. See the table https://software.intel.com/en-us/node/468532#TBL2-7 for a list of the permitted values.

Ok thanks. But what should I

Abdolreza_A_ — Tue, 09 May 2017 03:05:35 GMT

Ok thanks. But what should I do now??? My opinion : I should split "A" into 3 parts : L + D + U Then I muat call DDIASM for 3 times and at first call I muat introduce part "L", Then at second call I must introduce part " D" and at third call I must introduce part "U". Is this true??? And for obtaining the right solve for system of equations, what should I do? My opinion : I must add 3 solves from 3 call directly. Is this right???

Quote:(name withheld) wrote:

mecej4 — Tue, 09 May 2017 03:37:52 GMT

(name withheld) wrote:

My opinion : I should split "A" into 3 parts : L + D + U

Note that the L-D-U decomposition is A = L D U (the product), not L + D + U (the sum).

Have you analysed the structure of the equations that you obtained by discretizing the PDE and applying the boundary conditions? Is diagonal storage the best representation? If not, we do not need to be talking about SPBLAS routines.

Thank you very much. Yes I

Abdolreza_A_ — Tue, 09 May 2017 04:55:05 GMT

Thank you very much. Yes I analysed the structure of the equations and saw that diagonal storage format is almost most efficient approach for me. But because of my lack of mathematic knowledge (I'm sorry about that) I repeat my answer in other words : Should I now decompose "A" into 3 parts : "L" + "D" + "U" and then call DDIASM for 3 times ? ( 1st call for part "U" of "A" + 2nd call for part "D" of "A" + 3rd call for part "U" of "A" >>> so I will have 3 different answers which obtain from these 3 call respectively) then what should I do? Must I calculate sum of these 3 answers directly or must I use another way? I don't know what should I do exactly? Thanks a lot for your patience !!

If your simultaneous

mecej4 — Tue, 09 May 2017 16:36:38 GMT

If your simultaneous equations cannot be represented by a triangular matrix, you cannot use MKL_DDIAS. Whether diagonal storage is convenient or not is of little significance if you do not have a solver that matches your problem requirements.

Secondly, and I repeat what I said in #4, L, D and U are factors whose product, not sum, equals A.

Thank you very much indeed.

Abdolreza_A_ — Tue, 09 May 2017 19:04:11 GMT

Thank you very much indeed. So with thia kind of matrix which I have, which subroutine do you prefer?

Hello

Ying_H_Intel — Wed, 10 May 2017 03:38:36 GMT

Hello

Welcome to MKL forum. If you want to solve one sparse equation, you may need to use the PARDISO interface.

You can find all example code in MKL install folder. For example : C:\Program Files (x86)\IntelSWTools\compilers_and_libraries_2017.2.187\windows\mkl\examples\examples_core_f.zip => unzip it \solverf\source\ pardiso_unsym_f.f. Which will solve sparse matrix and get the X.

And you can read MKL user guide and developer guide to find more information:
https://software.intel.com/en-us/articles/intel-math-kernel-library-documentation

Best Regards,

Ying

!===============================================================================

*   Content : Intel(R) MKL PARDISO Fortran example
*
********************************************************************************
C----------------------------------------------------------------------
C Example program to show the use of the "PARDISO" routine
C for nonsymmetric linear systems
C---------------------------------------------------------------------
C This program can be downloaded from the following site:
C www.pardiso-project.org
C
C (C) Olaf Schenk, Department of Computer Science,
C University of Basel, Switzerland.
C Email: olaf.schenk@unibas.ch
C
C---------------------------------------------------------------------
      PROGRAM pardiso_unsym
        IMPLICIT NONE
        include 'mkl_pardiso.fi'
C.. Internal solver memory pointer for 64-bit architectures
C.. INTEGER*8 pt(64)
C.. Internal solver memory pointer for 32-bit architectures
C.. INTEGER*4 pt(64)
C.. This is OK in both cases
        TYPE(MKL_PARDISO_HANDLE) pt(64)
C.. All other variables
        INTEGER maxfct, mnum, mtype, phase, n, nrhs, error, msglvl
        INTEGER iparm(64)
        INTEGER ia(6)
        INTEGER ja(13)
        REAL*8 a(13)
        REAL*8 b(5)
        REAL*8 x(5)
        INTEGER i, idum(1)
        REAL*8 ddum(1)
C.. Fill all arrays containing matrix data.
        DATA n /5/, nrhs /1/, maxfct /1/, mnum /1/
        DATA ia /1,4,6,9,12,14/
        DATA ja
     1 /1,2, 4,
     2   1,2,
     3       3,4,5,
     4   1, 3,4,
     5     2,    5/
        DATA a
     1 / 1.d0,-1.d0,      -3.d0,
     2   -2.d0, 5.d0,
     3                4.d0, 6.d0, 4.d0,
     4   -4.d0,       2.d0, 7.d0,
     5          8.d0,            -5.d0/
C..
C.. Setup PARDISO control parameter
C..
        DO i = 1, 64
           iparm(i) = 0
        END DO
        iparm(1) = 1 ! no solver default
        iparm(2) = 2 ! fill-in reordering from METIS
        iparm(3) = 1 ! numbers of processors
        iparm(4) = 0 ! no iterative-direct algorithm
        iparm(5) = 0 ! no user fill-in reducing permutation
        iparm(6) = 0 ! =0 solution on the first n components of x
        iparm(7) = 0 ! not in use
        iparm(8) = 9 ! numbers of iterative refinement steps
        iparm(9) = 0 ! not in use
        iparm(10) = 13 ! perturb the pivot elements with 1E-13
        iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
        iparm(12) = 0 ! not in use
        iparm(13) = 1 ! maximum weighted matching algorithm is switched-on (default for non-symmetric)
        iparm(14) = 0 ! Output: number of perturbed pivots
        iparm(15) = 0 ! not in use
        iparm(16) = 0 ! not in use
        iparm(17) = 0 ! not in use
        iparm(18) = -1 ! Output: number of nonzeros in the factor LU
        iparm(19) = -1 ! Output: Mflops for LU factorization
        iparm(20) = 0 ! Output: Numbers of CG Iterations
        error = 0 ! initialize error flag
        msglvl = 1 ! print statistical information
        mtype = 11 ! real unsymmetric
C.. Initialize the internal solver memory pointer. This is only
C necessary for the FIRST call of the PARDISO solver.
        DO i = 1, 64
            pt(i)%DUMMY = 0
        END DO
C.. Reordering and Symbolic Factorization, This step also allocates
C   all memory that is necessary for the factorization
        phase = 11 ! only reordering and symbolic factorization
        CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     & idum, nrhs, iparm, msglvl, ddum, ddum, error)
        WRITE(*,*) 'Reordering completed ... '
        IF (error .NE. 0) THEN
            WRITE(*,*) 'The following ERROR was detected: ', error
            STOP 1
        END IF
        WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
        WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
C.. Factorization.
        phase = 22 ! only factorization
        CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     & idum, nrhs, iparm, msglvl, ddum, ddum, error)
        WRITE(*,*) 'Factorization completed ... '
        IF (error .NE. 0) THEN
            WRITE(*,*) 'The following ERROR was detected: ', error
            STOP 1
        END IF
C.. Back substitution and iterative refinement
        iparm(8) = 2 ! max numbers of iterative refinement steps
        phase = 33 ! only factorization
        DO i = 1, n
            b(i) = 1.0D0
        END DO
        CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     & idum, nrhs, iparm, msglvl, b, x, error)
        WRITE(*,*) 'Solve completed ... '
        WRITE(*,*) 'The solution of the system is '
        DO i = 1, n
            WRITE(*,*) ' x(',i,') = ', x(i)
        END DO
C.. Termination and release of memory
        phase = -1 ! release internal memory
        CALL pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum,
     & idum, idum, nrhs, iparm, msglvl, ddum, ddum, error)
      END PROGRAM pardiso_unsym

thanks you very much indeed

Abdolreza_A_ — Wed, 10 May 2017 16:02:27 GMT

thanks you very much indeed Ying H. !! So Should I use this PARDISO directly to my code without any editing on ipharms??

Hello,

Ying_H_Intel — Thu, 11 May 2017 00:56:55 GMT

Hello,

Right, you can use the code directly with your input sparse matrix ( ia, ja and a). and later, you may change the parameters of MKL function according to your needs and matrix type.

You may refer to MKL developer reference manual and MKL developer guide, which are basic docs for MKL usage and MKL functionality.

https://software.intel.com/en-us/articles/intel-math-kernel-library-documentation

Best Regards,

Ying