topic PROGRAM DGESV_VS_DDNSCSR in Intel® oneAPI Math Kernel Library

problem with PARDISO and DSS

Abdolreza_A_ — Fri, 11 Aug 2017 09:14:42 GMT

Hello my dear friends ...

I wrote a FORTRAN code for simulating a transient 2D laminar flow in fintite differences method.

I have 2 problems with my code. One of them is with PARDISO and the other is with DSS.

for PARDISO :
My code solves 3 systems of linear equations respectively in each iteration :

The matrices of the coefficients for the first system and the second one are the same. The systems' difference is just the RHS.
The matrices of the coefficients for the third system is different. but the it's dimension is the same.
I used 2 approaches (i.e. 2 SOLVERs in one code in eatch iteration to campare the results of the same system of equations between them):

In the first part of each iteration, I make the matrix of coefficients of each of the 3 systems of equations "by myself" (for example A , B , C) and then I send the system to DGESV to solve. Fortunately it gives me the best results.
In the second part of each iteration, I send the the matrix of coefficients (A , B , C) to the DDNSCSR to make it in CSR form.Then I send it to the PARDISO and it gives me almost the same results as which DGESV had returnet to me.
MY PROBLEM IS : when I make the CSR format by myself and send it to the PARDISO, the results are not the same !!!

I should note that the CSR format which I made is EXACTLY ( !!! ) the same as the CSR format which was made by DDNSCSR ( !!! )
I mean exactly : the inlet CSR format matrices which was sent to PARIDISO, are exactly the same in 2 approach. (I had checked it completely with a separate code which I had writtenbefore)

----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

for DSS : I don't know anything about the applying the DSS solver. I will send my code (it is the same code but I substitude PARDISO with DSS). I don't know why it didn't work at all. ( I gave the solver from internet)
Thank for your help.

PROGRAM DGESV_VS_DDNSCSR

Abdolreza_A_ — Fri, 11 Aug 2017 09:16:34 GMT

PROGRAM DGESV_VS_DDNSCSR_PARDISO

IMPLICIT NONE

! ****************************************************************************************

CHARACTER (len=30) :: FILE_NAME_1 , FILE_NAME_2

CHARACTER (len=1024) :: FINAL_NAME_1 , FINAL_NAME_2

CHARACTER (len=1024) :: FORMAT_STRING

! ****************************** ORDINARY GENERAL VARIABLES ******************************

INTEGER, PARAMETER :: N = 51 ! GENERAL N MUST BE ODD

INTEGER, PARAMETER :: NN = N**2 ! GENERAL

INTEGER, PARAMETER :: N2 = (N/2)+1 ! GENERAL

DOUBLE PRECISION , PARAMETER :: DELTA_X = 1.0D0 / (N-1) ! GENERAL

DOUBLE PRECISION , PARAMETER :: DELTA_Y = 1.0D0 / (N-1) ! GENERAL

DOUBLE PRECISION , PARAMETER :: DELTA_T = 0.01D0 ! GENERAL

DOUBLE PRECISION , PARAMETER :: REYNOLDS = 160.0D0 ! GENERAL

INTEGER :: I , J , ITERATION , IPRINT ! GENERAL

DOUBLE PRECISION :: XITERATION ! GENERAL

INTEGER :: INFO_DGESV , INFO_DDNSCSR

DOUBLE PRECISION :: C_A_DGESV , C_B_DGESV , C_C_DGESV , C_D_DGESV , C_E_DGESV , SUM_DGESV , TEMP1_DGESV

DOUBLE PRECISION , DIMENSION ( NN ) :: RHS_U_DGESV , RHS_V_DGESV , RHS_P_DGESV , IPIV_DGESV , TEMPORARY_DGESV

DOUBLE PRECISION , DIMENSION ( N , N ) :: U_N_DGESV , V_N_DGESV , U_INT_DGESV , V_INT_DGESV , U_N1_DGESV , V_N1_DGESV , P_N1_DGESV

DOUBLE PRECISION , DIMENSION ( NN , NN ) :: A_DGESV , B_DGESV , C_DGESV

INTEGER :: INFO

DOUBLE PRECISION :: C_A , C_B , C_C , C_D , C_E , SUM

DOUBLE PRECISION , DIMENSION ( NN ) :: RHS_P

DOUBLE PRECISION , DIMENSION ( N , N ) :: U_N , V_N , U_INT , V_INT , U_N1 , V_N1 , P_N1

DOUBLE PRECISION , DIMENSION ( NN , 2 ) :: RHS_U_V

INTEGER , PARAMETER :: NON_ZERO_OF_A = 5 * (N**2) - (10*N) + 8 ! (12503)

INTEGER , PARAMETER :: NON_ZERO_OF_B = 5 * (N**2) - ( 4*N) - 4 ! (12797)

! ****************************** FOR "DDNSCSR" TRANSFORMER SUBROUTINE ******************************

INTEGER , DIMENSION (6) :: JOB_A_DDNSCSR , JOB_B_DDNSCSR

DOUBLE PRECISION , DIMENSION (NON_ZERO_OF_A) :: A_CSR_DDNSCSR ! == A_CSR FOR "DDNSCSR" TRANSFORMER SUBROUTINE

DOUBLE PRECISION , DIMENSION (NON_ZERO_OF_B) :: B_CSR_DDNSCSR ! == B_CSR FOR "DDNSCSR" TRANSFORMER SUBROUTINE

INTEGER , DIMENSION (NON_ZERO_OF_A) :: J_A_DDNSCSR ! == COLUMN_A FOR "DDNSCSR" TRANSFORMER SUBROUTINE

INTEGER , DIMENSION (NON_ZERO_OF_B) :: J_B_DDNSCSR ! == COLUMN_B FOR "DDNSCSR" TRANSFORMER SUBROUTINE

INTEGER , DIMENSION ( NN + 1 ) :: I_A_DDNSCSR ! == RowINDEX_A FOR "DDNSCSR" TRANSFORMER SUBROUTINE

INTEGER , DIMENSION ( NN + 1 ) :: I_B_DDNSCSR ! == RowINDEX_B FOR "DDNSCSR" TRANSFORMER SUBROUTINE

! ****************************** TO DEFINE THE VARIABLES NEEDED FOR PARDISO ******************************

! ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR, IndexROW_A, COLUMN_A, PERM_A, NRHS_A, IPARM_A,

! MSGLVL_A, B_A, X_A, ERROR_A)

INTEGER*8 , DIMENSION ( 64 ) :: PT_A = 0

INTEGER :: MAXFCT_A = 1

INTEGER :: MNUM_A = 1

INTEGER :: MTYPE_A = 11

INTEGER :: PHASE_A

INTEGER , PARAMETER :: N_A = NN

INTEGER , DIMENSION ( N_A ) :: PERM_A

INTEGER , PARAMETER :: NRHS_A = 2

INTEGER , DIMENSION ( 64 ) :: IPARM_A = 0

INTEGER :: MSGLVL_A = 1

DOUBLE PRECISION , DIMENSION ( N_A * NRHS_A ) :: B_A

DOUBLE PRECISION , DIMENSION ( N_A * NRHS_A ) :: X_A

INTEGER :: ERROR_A = 0

! ****************************** TO DEFINE THE VARIABLES NEEDED FOR PARDISO ******************************

! ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR, IndexROW_B, COLUMN_B, PERM_B, NRHS_B, IPARM_B,

! MSGLVL_B, B_B, X_B, ERROR_B)

INTEGER*8 , DIMENSION ( 64 ) :: PT_B = 0

INTEGER :: MAXFCT_B = 1

INTEGER :: MNUM_B = 1

INTEGER :: MTYPE_B = 11

INTEGER :: PHASE_B

INTEGER , PARAMETER :: N_B = NN

INTEGER , DIMENSION ( N_B ) :: PERM_B

INTEGER , PARAMETER :: NRHS_B = 1

INTEGER , DIMENSION ( 64 ) :: IPARM_B = 0

INTEGER :: MSGLVL_B = 1

DOUBLE PRECISION , DIMENSION ( N_B * NRHS_B ) :: B_B

DOUBLE PRECISION , DIMENSION ( N_B * NRHS_B ) :: X_B

INTEGER :: ERROR_B = 0

! ****************************** TO MAKE THE MARICES WHICH WE NEED ******************************

JOB_A_DDNSCSR (1) = 0

JOB_A_DDNSCSR (2) = 1

JOB_A_DDNSCSR (3) = 1

JOB_A_DDNSCSR (4) = 2

JOB_A_DDNSCSR (5) = NON_ZERO_OF_A ! NON_ZERO_OF_A = 5 * (N**2) - (10*N) + 8 = 12503

JOB_A_DDNSCSR (6) = 1

JOB_B_DDNSCSR (1) = 0

JOB_B_DDNSCSR (2) = 1

JOB_B_DDNSCSR (3) = 1

JOB_B_DDNSCSR (4) = 2

JOB_B_DDNSCSR (5) = NON_ZERO_OF_B ! NON_ZERO_OF_B = 5 * (N**2) - ( 4*N) - 4 = 12797

JOB_B_DDNSCSR (6) = 1

U_N_DGESV = 0.0D0

V_N_DGESV = 0.0D0

U_INT_DGESV = 0.0D0

V_INT_DGESV = 0.0D0

U_N1_DGESV = 0.0D0

V_N1_DGESV = 0.0D0

P_N1_DGESV = 0.0D0

U_N = 0.0D0

V_N = 0.0D0

U_INT = 0.0D0

V_INT = 0.0D0

U_N1 = 0.0D0

V_N1 = 0.0D0

P_N1 = 0.0D0

IPRINT = 1

DO ITERATION = 1 , 1000

!**************************************** TO MAKE THE MATRICES ( IN EACH ITERATION )

A_DGESV = 0.0D0

B_DGESV = 0.0D0

C_DGESV = 0.0D0

RHS_U_DGESV = 0.0D0

RHS_V_DGESV = 0.0D0

RHS_P_DGESV = 0.0D0

RHS_U_V = 0.0D0

RHS_P = 0.0D0

! ****************************** TO MAKE PART 1 : DIRICHLET BOUNDARY FOR " U " AND " V " ON THE LEFT SIDE OF THE DOMAIN

DO I = 1 , N

U_N_DGESV ( I , 1 ) = 0.0D0

V_N_DGESV ( I , 1 ) = 0.0D0

END DO

U_N_DGESV ( N2 - 1 , 1 ) = 0.75D0 * 2.0D0

U_N_DGESV ( N2 , 1 ) = 2.0D0

U_N_DGESV ( N2 + 1 , 1 ) = 0.75D0 * 2.0D0

DO I = 1 , N

A_DGESV ( (I-1)*N + 1 , (I-1)*N + 1 ) = 1.0D0

RHS_U_DGESV ( (I-1)*N + 1 ) = 0.0D0

RHS_V_DGESV ( (I-1)*N + 1 ) = 0.0D0

END DO

RHS_U_DGESV ( (N2-2)*N + 1 ) = 0.75D0 * 2.0D0

RHS_U_DGESV ( (N2-1)*N + 1 ) = 2.0D0

RHS_U_DGESV ( (N2-0)*N + 1 ) = 0.75D0 * 2.0D0

! ****************************** TO MAKE PART 1 FOR "PARDISO"

DO I = 1 , N

U_N ( I , 1 ) = 0.0D0

V_N ( I , 1 ) = 0.0D0

END DO

U_N ( N2 - 1 , 1 ) = 0.75D0 * 2.0D0

U_N ( N2 , 1 ) = 2.0D0

U_N ( N2 + 1 , 1 ) = 0.75D0 * 2.0D0

RHS_U_V ( (N2-2)*N + 1 , 1 ) = 0.75D0 * 2.0D0

RHS_U_V ( (N2-1)*N + 1 , 1 ) = 2.0D0

RHS_U_V ( (N2-0)*N + 1 , 1 ) = 0.75D0 * 2.0D0

! ****************************** TO MAKE PART 2 : ALL NODES OF THE DOMAIN EXCEPT THE NODES WHICH ARE ON THE BOUNDARIES

DO I = 2 , N-1

DO J = 2 , N-1

C_A_DGESV = ( V_N_DGESV (I,J) / ( 2 * DELTA_Y ) ) - ( ( 1 / REYNOLDS ) * ( 1 / ( DELTA_Y * DELTA_Y ) ) ) ! UPPER NODE

C_B_DGESV = -1 * ( U_N_DGESV (I,J) / ( 2 * DELTA_X ) ) - ( (1 / REYNOLDS) * ( 1/(DELTA_X * DELTA_X ) ) ) ! LEFT NODE

C_C_DGESV = ( 1 / DELTA_T ) + ( ( 1 / REYNOLDS ) * ( ( 2 / ( DELTA_X * DELTA_X ) ) + ( 2 / ( DELTA_Y * DELTA_Y ) ) ) ) ! THE NODE ITSELF

C_D_DGESV = ( U_N_DGESV (I,J) / ( 2 * DELTA_X ) ) - ( ( 1 / REYNOLDS ) * ( 1 / ( DELTA_X * DELTA_X ) ) ) ! RIGHT NODE

C_E_DGESV = -1 * ( V_N_DGESV (I,J) / ( 2 * DELTA_Y ) ) - ( ( 1 / REYNOLDS ) * ( 1 / (DELTA_Y * DELTA_Y ) ) ) ! LOWER NODE

A_DGESV ( (I-1)*N + J , (I-1)*N + J - N ) = C_A_DGESV ! COEFFICIENT FOR THE UPPER NODE : ( I - 1 , J )

A_DGESV ( (I-1)*N + J , (I-1)*N + J - 1 ) = C_B_DGESV ! COEFFICIENT FOR THE LEFT NODE : ( I , J - 1 )

A_DGESV ( (I-1)*N + J , (I-1)*N + J ) = C_C_DGESV ! COEFFICIENT FOR THE NODE ITSELF : ( I , J )

A_DGESV ( (I-1)*N + J , (I-1)*N + J + 1 ) = C_D_DGESV ! COEFFICIENT FOR THE RIGHT NODE : ( I , J + 1 )

A_DGESV ( (I-1)*N + J , (I-1)*N + J + N ) = C_E_DGESV ! COEFFICIENT FOR THE LOWER NODE : ( I + 1 , J )

RHS_U_DGESV ( (I-1)*N + J ) = U_N_DGESV ( I , J ) / DELTA_T

RHS_V_DGESV ( (I-1)*N + J ) = V_N_DGESV ( I , J ) / DELTA_T

END DO

! ****************************** TO MAKE PART 2 FOR "PARDISO"

DO I = 2 , N-1

DO J = 2 , N-1

C_A = ( V_N (I,J) / ( 2 * DELTA_Y ) ) - ( ( 1 / REYNOLDS ) * ( 1 / ( DELTA_Y * DELTA_Y ) ) ) ! UPPER NODE

C_B = -1 * ( U_N (I,J) / ( 2 * DELTA_X ) ) - ( ( 1 / REYNOLDS ) * ( 1 / ( DELTA_X * DELTA_X ) ) ) ! LEFT NODE

C_C = ( 1 / DELTA_T ) + ( ( 1 / REYNOLDS ) * ( ( 2 / ( DELTA_X * DELTA_X ) ) + ( 2 / ( DELTA_Y * DELTA_Y ) ) ) ) ! THE NODE ITSELF

C_D = ( U_N (I,J) / ( 2 * DELTA_X ) ) - ( ( 1 / REYNOLDS ) * ( 1 / ( DELTA_X * DELTA_X ) ) ) ! RIGHT NODE

C_E = -1 * ( V_N (I,J) / ( 2 * DELTA_Y ) ) - ( ( 1 / REYNOLDS ) * ( 1 / ( DELTA_Y * DELTA_Y ) ) ) ! LOWER NODE

RHS_U_V ( (I-1)*N + J , 1 ) = U_N ( I , J ) / DELTA_T

RHS_U_V ( (I-1)*N + J , 2 ) = V_N ( I , J ) / DELTA_T

END DO

! ****************************** TO MAKE PART 3 : THE HOMOGENOUS BOUNDARY CONDITION FOR THE UPPER & THE LOWER SIDE OF THE DOMAIN

! (ROUND U / ROUND Y = ROUND V / ROUND Y = 0) (2nd Order Neumann)

DO J = 2 , N-1

A_DGESV ( J , J ) = -3.0D0 ! FOR POINT ( 1 , J )

A_DGESV ( J , J + N ) = 4.0D0 ! FOR POINT ( 2 , J )

A_DGESV ( J , J + 2*N ) = -1.0D0 ! FOR POINT ( 3 , J )

RHS_U_DGESV ( J ) = 0.0D0

RHS_V_DGESV ( J ) = 0.0D0

A_DGESV ( (N-1)*N + J , (N-1)*N + J ) = -3.0D0 ! FOR POINT ( N , J )

A_DGESV ( (N-1)*N + J , (N-1)*N + J - N ) = 4.0D0 ! FOR POINT ( N - 1 , J )

A_DGESV ( (N-1)*N + J , (N-1)*N + J - 2*N ) = -1.0D0 ! FOR POINT ( N - 2 , J )

RHS_U_DGESV ( (N-1)*N + J ) = 0.0D0

RHS_V_DGESV ( (N-1)*N + J ) = 0.0D0

END DO

! ****************************** TO MAKE PART 3 FOR "PARDISO"

DO J = 2 , N-1

RHS_U_V ( J , 1 ) = 0.0D0

RHS_U_V ( J , 2 ) = 0.0D0

RHS_U_V ( (N-1)*N + J , 1 ) = 0.0D0

RHS_U_V ( (N-1)*N + J , 2 ) = 0.0D0

END DO

! ****************************** TO MAKE PART 4 : THE HOMOGENOUS BOUNDARY CONDITION ON THE RIGHT SIDE OF THE DOMAIN

! ( ROUND U / ROUND Y = 0 ) & ( ROUND V / ROUND X = 0 ) (2nd Order Neumann)

DO I = 1 , N

A_DGESV ( I*N , I*N ) = -3.0D0 ! FOR POINT ( I , N )

A_DGESV ( I*N , I*N - 1 ) = 4.0D0 ! FOR POINT ( I , N-1 )

A_DGESV ( I*N , I*N - 2 ) = -1.0D0 ! FOR POINT ( I , N-2 )

RHS_U_DGESV ( I*N ) = 0.0D0

RHS_V_DGESV ( I*N ) = 0.0D0

END DO

! ****************************** TO MAKE PART 4 FOR "PARDISO"

DO I = 1 , N

RHS_U_V ( I*N , 1 ) = 0.0D0

RHS_U_V ( I*N , 2 ) = 0.0D0

END DO

!******************** INITIALIZING THE VARIABLES WHICH NEEDED FOR CALLING PARDISO (FOR VELOCITIES)

IPARM_A ( 1 ) = 1

IPARM_A ( 2 ) = 3

IPARM_A ( 3 ) = 1

IPARM_A ( 4 ) = 0

IPARM_A ( 5 ) = 0

IPARM_A ( 6 ) = 0

IPARM_A ( 7 ) = 0

IPARM_A ( 8 ) = -1000 ! 10

IPARM_A ( 9 ) = 0

IPARM_A ( 10 ) = 13

IPARM_A ( 11 ) = 0

IPARM_A ( 12 ) = 0

IPARM_A ( 13 ) = 0

IPARM_A ( 14 ) = 0

IPARM_A ( 15 ) = 0

IPARM_A ( 16 ) = 0

IPARM_A ( 17 ) = 0

IPARM_A ( 18 ) = 0

IPARM_A ( 19 ) = 0

IPARM_A ( 20 ) = -1 ! 0

IPARM_A ( 21 ) = 0

IPARM_A ( 22 ) = 0

IPARM_A ( 23 ) = 0

IPARM_A ( 24 ) = 1 ! 0

IPARM_A ( 25 ) = 0

IPARM_A ( 26 ) = 0

IPARM_A ( 27 ) = 1

IPARM_A ( 28 ) = 0

IPARM_A ( 29 ) = 0

IPARM_A ( 30 ) = 0

IPARM_A ( 31 ) = 0

IPARM_A ( 32 ) = 0

IPARM_A ( 33 ) = 0

IPARM_A ( 34 ) = 0

IPARM_A ( 35 ) = 0

IPARM_A ( 36 ) = 0

IPARM_A ( 37 ) = 0

IPARM_A ( 38 ) = 0

IPARM_A ( 39 ) = 0

IPARM_A ( 40 ) = 0

IPARM_A ( 41 ) = 0

IPARM_A ( 42 ) = 0

IPARM_A ( 43 ) = 0

IPARM_A ( 44 ) = 0

IPARM_A ( 45 ) = 0

IPARM_A ( 46 ) = 0

IPARM_A ( 47 ) = 0

IPARM_A ( 48 ) = 0

IPARM_A ( 49 ) = 0

IPARM_A ( 50 ) = 0

IPARM_A ( 51 ) = 0

IPARM_A ( 52 ) = 0

IPARM_A ( 53 ) = 0

IPARM_A ( 54 ) = 0

IPARM_A ( 55 ) = 0

IPARM_A ( 56 ) = 0

IPARM_A ( 57 ) = 0

IPARM_A ( 58 ) = 0

IPARM_A ( 59 ) = 0

IPARM_A ( 60 ) = 0

IPARM_A ( 61 ) = 0

IPARM_A ( 62 ) = 0

IPARM_A ( 63 ) = 0

IPARM_A ( 64 ) = 0

! ****************************** TO MAKE THE RIGHT HAND SIDE MATRIX "B_A" FOR PARDISO (FOR VELOCITIES)

DO I = 1 , NN

DO J = 1 , 2

B_A ( (J-1)*NN + I ) = RHS_U_V ( I , J )

END DO

! ****************************** TO SOLVE THE SYSTEM OF EQUATIONS TO OBTAIN INTERMEDIATE VELOCITIES

INFO_DDNSCSR = 0

CALL MKL_DDNSCSR ( JOB_A_DDNSCSR , NN , NN , A_DGESV , NN , A_CSR_DDNSCSR , J_A_DDNSCSR , I_A_DDNSCSR , INFO_DDNSCSR )

C_DGESV = A_DGESV ! REPLACE THE ARRAYS OF MATRICE "A" INTO MATRICE "C" BECAUSE "DGESV" WILL CHANGE "A" !!!

INFO_DGESV = 0

CALL DGESV ( NN , 1 , A_DGESV , NN , IPIV_DGESV , RHS_U_DGESV , NN , INFO_DGESV )

INFO_DGESV = 0

CALL DGESV ( NN , 1 , C_DGESV , NN , IPIV_DGESV , RHS_V_DGESV , NN , INFO_DGESV )

DO I = 1 , N

DO J = 1 , N

U_INT_DGESV ( I , J ) = RHS_U_DGESV ( (I-1)*N + J )

V_INT_DGESV ( I , J ) = RHS_V_DGESV ( (I-1)*N + J )

END DO

DO I = 1 , N

U_INT_DGESV ( I , 1 ) = 0.0D0

V_INT_DGESV ( I , 1 ) = 0.0D0

END DO

U_INT_DGESV ( N2 - 1 , 1 ) = 0.75D0 * 2.0D0

U_INT_DGESV ( N2 , 1 ) = 2.0D0

U_INT_DGESV ( N2 + 1 , 1 ) = 0.75D0 * 2.0D0

!-----------------------------------------------------------------

! Only Analysis, Reordering And Symbolic Factorization, This Step Also Allocates All Memory That Is Necessary For The Factorization :

PHASE_A = 11

CALL PARDISO ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR_DDNSCSR, I_A_DDNSCSR, J_A_DDNSCSR, PERM_A, NRHS_A, IPARM_A, MSGLVL_A, B_A, X_A, ERROR_A)

! CALL PARDISO ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR, RowINDEX_A, COLUMN_A, PERM_A, NRHS_A, IPARM_A, MSGLVL_A, B_A, X_A, ERROR_A)

IF (ERROR_A .NE. 0) THEN

WRITE ( * , * ) 'The following ERROR was detected in Analysis, Reordering and Symbolic Factorization for solving VELOCITIES at ITERATION : ', ERROR_A , ITERATION

READ ( * , * )

END IF

!-----------------------------------------------------------------

! Only Numerical Factorization :

PHASE_A = 22

CALL PARDISO ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR_DDNSCSR, I_A_DDNSCSR, J_A_DDNSCSR, PERM_A, NRHS_A, IPARM_A, MSGLVL_A, B_A, X_A, ERROR_A)

! CALL PARDISO ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR, RowINDEX_A, COLUMN_A, PERM_A, NRHS_A, IPARM_A, MSGLVL_A, B_A, X_A, ERROR_A)

IF (ERROR_A .NE. 0) THEN

WRITE ( * , * ) 'The following ERROR was detected in Numerical Factorization for solving VELOCITIES at ITERATION : ', ERROR_A , ITERATION

READ ( * , * )

END IF

!-----------------------------------------------------------------

! Solving With Back Substitution And Iterative Refinement :

PHASE_A = 33

CALL PARDISO ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR_DDNSCSR, I_A_DDNSCSR, J_A_DDNSCSR, PERM_A, NRHS_A, IPARM_A, MSGLVL_A, B_A, X_A, ERROR_A)

! CALL PARDISO ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR, RowINDEX_A, COLUMN_A, PERM_A, NRHS_A, IPARM_A, MSGLVL_A, B_A, X_A, ERROR_A)

!-----------------------------------------------------------------

! Termination And Release Internal Memory :

PHASE_A = -1

CALL PARDISO ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR_DDNSCSR, I_A_DDNSCSR, J_A_DDNSCSR, PERM_A, NRHS_A, IPARM_A, MSGLVL_A, B_A, X_A, ERROR_A)

! CALL PARDISO ( PT_A, MAXFCT_A, MNUM_A, MTYPE_A, PHASE_A, N_A, A_CSR, RowINDEX_A, COLUMN_A, PERM_A, NRHS_A, IPARM_A, MSGLVL_A, B_A, X_A, ERROR_A)

! ****************************** TO MAKE THE INTERMEDIATE VELOCITIES "U_INT" , "V_INT" FOR PARDISO

DO I = 1 , N

DO J = 1 , N

U_INT ( I , J ) = X_A ( (I-1)*N + J )

V_INT ( I , J ) = X_A ( (I-1)*N + J + NN )

END DO

DO I = 1 , N

U_INT ( I , 1 ) = 0.0D0

V_INT ( I , 1 ) = 0.0D0

END DO

U_INT ( N2 - 1 , 1 ) = 0.75D0 * 2.0D0

U_INT ( N2 , 1 ) = 2.0D0

U_INT ( N2 + 1 , 1 ) = 0.75D0 * 2.0D0

!***************************** TO CALCULATE THE EDITED PRESSURE "P_N1" FROM OBTAINED INTERMEDIATE VELOCITIES

!***************************** (AT THE END OF THE FIRST STEP IN EACH IETRATION)

SUM_DGESV = 0.0D0

DO I = 2 , N-1

DO J = 2 , N-1

RHS_P_DGESV ( (I-1)*N + J ) = (DELTA_X / (2*DELTA_T)) * ( U_INT_DGESV (I,J+1) - U_INT_DGESV (I,J-1) + V_INT_DGESV (I-1,J) - V_INT_DGESV (I+1,J) )

SUM_DGESV = SUM_DGESV + RHS_P_DGESV ( N*(I-1)+J )

END DO

SUM_DGESV = SUM_DGESV / ((N - 2)*(N - 2))

RHS_P_DGESV = RHS_P_DGESV - SUM_DGESV

!***************************** FOR PARDISO

SUM = 0.0D0

DO I = 2 , N-1

DO J = 2 , N-1

RHS_P ( (I-1)*N + J ) = (DELTA_X / (2*DELTA_T)) * ( U_INT(I,J+1) - U_INT(I,J-1) + V_INT(I-1,J) - V_INT(I+1,J) )

SUM = SUM + RHS_P ( (I-1)*N + J )

END DO

SUM = SUM / ( (N-2) * (N-2) )

RHS_P = RHS_P - SUM

!********** TO DETERMINE THE OTHER COEFFICIENTS IN THE MATRIX OF COEFFICCNETS FOR PRESSURES IN EACH ITERATION

DO I = 2 , N-1

DO J = 2 , N-1

B_DGESV ( (I-1)*N + J , (I-1)*N + J - N ) = 1 ! COEFFICIENT FOR THE UPPER NODE : ( I - 1 , J )

B_DGESV ( (I-1)*N + J , (I-1)*N + J - 1 ) = 1 ! COEFFICIENT FOR THE LEFT NODE : ( I , J - 1 )

B_DGESV ( (I-1)*N + J , (I-1)*N + J ) = -4 ! COEFFICIENT FOR THE NODE ITSELF : ( I , J )

B_DGESV ( (I-1)*N + J , (I-1)*N + J + 1 ) = 1 ! COEFFICIENT FOR THE RIGHT NODE : ( I , J + 1 )

B_DGESV ( (I-1)*N + J , (I-1)*N + J + N ) = 1 ! COEFFICIENT FOR THE LOWER NODE : ( I + 1 , J )

END DO

DO I = N2 , N2

DO J = N2 , N2

B_DGESV ( (I-1)*N + J , (I-1)*N + J - N ) = 0 ! COEFFICIENT FOR THE UPPER NODE : ( I - 1 , J )

B_DGESV ( (I-1)*N + J , (I-1)*N + J - 1 ) = 0 ! COEFFICIENT FOR THE LEFT NODE : ( I , J - 1 )

B_DGESV ( (I-1)*N + J , (I-1)*N + J ) = 1 ! COEFFICIENT FOR THE NODE ITSELF : ( I , J )

B_DGESV ( (I-1)*N + J , (I-1)*N + J + 1 ) = 0 ! COEFFICIENT FOR THE RIGHT NODE : ( I , J + 1 )

B_DGESV ( (I-1)*N + J , (I-1)*N + J + N ) = 0 ! COEFFICIENT FOR THE LOWER NODE : ( I + 1 , J )

RHS_P_DGESV ( (I-1)*N + J ) = 0

END DO

! ****************************** TO MAKE PART 0 : THE ARTIFICIAL BOUNDARY CONDITION FOR PRESSURES

! ****************************** ON THE LEFT SIDE AND THE RIGHT SIDE OF THE DOMAIN

DO I = 2 , N-1

B_DGESV ( (I-1)*N + 1 , (I-1)*N + 1 ) = -4.0D0 ! FOR POINT ( I , 1 )

B_DGESV ( (I-1)*N + 1 , (I-1)*N + 1 - N ) = 1.0D0 ! FOR POINT ( I - 1 , 1 )

B_DGESV ( (I-1)*N + 1 , (I-1)*N + 1 + N ) = 1.0D0 ! FOR POINT ( I + 1 , 1 )

B_DGESV ( (I-1)*N + 1 , (I-1)*N + 1 + 1 ) = 2.0D0 ! FOR POINT ( I , 2 )

RHS_P_DGESV ( (I-1)*N + 1 ) = 0.0D0

END DO

DO I = 2 , N-1

B_DGESV ( I*N , I*N ) = -4.0D0 ! FOR POINT ( I , N )

B_DGESV ( I*N , I*N - 1 ) = 2.0D0 ! FOR POINT ( I , N-1 )

B_DGESV ( I*N , I*N - N ) = 1.0D0 ! FOR POINT ( I - 1 , N )

B_DGESV ( I*N , I*N + N ) = 1.0D0 ! FOR POINT ( I + 1 , N )

RHS_P_DGESV ( I*N ) = 0.0D0

END DO

! ****************************** TO MAKE PART 0 FOR PARDISO

! IT DID IN "B_CSR", "COLUMN_B" AND "RowINDEX_B"

! ****************************** TO MAKE PART 1 FOR PARDISO

DO I = 2 , N-1 ! FOR THE LEFT SIDE OF THE DOMAIN

RHS_P ( (I-1)*N + 1 ) = 0.0D0

END DO

DO I = 2 , N-1 ! FOR THE RIGHT SIDE OF THE DOMAIN

RHS_P ( I*N ) = 0.0D0

END DO

! ****************************** TO MAKE PART 2 : THE ARTIFICIAL BOUNDARY CONDITION FOR PRESSURES

! ****************************** ON THE UPPER SIDE AND THE LOWER SIDE OF THE DOMAIN

DO J = 2 , N-1

B_DGESV ( J , J ) = -4.0D0 ! FOR POINT ( 1 , J )

B_DGESV ( J , J + N ) = 2.0D0 ! FOR POINT ( 2 , J )

B_DGESV ( J , J + 1 ) = 1.0D0 ! FOR POINT ( 1 , J + 1 )

B_DGESV ( J , J - 1 ) = 1.0D0 ! FOR POINT ( 1 , J - 1 )

RHS_P_DGESV ( J ) = 0.0D0

B_DGESV ( (N-1)*N + J , (N-1)*N + J ) = -4.0D0 ! FOR POINT ( N , J )

B_DGESV ( (N-1)*N + J , (N-1)*N + J + 1 ) = 1.0D0 ! FOR POINT ( N , J + 1 )

B_DGESV ( (N-1)*N + J , (N-1)*N + J - 1 ) = 1.0D0 ! FOR POINT ( N , J - 1 )

B_DGESV ( (N-1)*N + J , (N-1)*N + J - N ) = 2.0D0 ! FOR POINT ( N - 1 , J )

RHS_P_DGESV ( (N-1)*N + J ) = 0.0D0

END DO

! ****************************** TO MAKE PART 2 FOR PARDISO

DO J = 2 , N-1

RHS_P ( J ) = 0.0D0

RHS_P ( (N-1)*N + J ) = 0.0D0

END DO

! ****************************** TO MAKE PART 3 : AVERAGING OF TWO ADJACENT NODES FOR 4 CORNERS

!*** FOR POINT ( 1 , 1 ) OF THE DOMAIN

B_DGESV ( 1 , 1 ) = 1.0D0 ! FOR POINT ( 1 , 1 )

B_DGESV ( 1 , 2 ) = -0.5D0 ! FOR POINT ( 1 , 2 )

B_DGESV ( 1 , 1 + N ) = -0.5D0 ! FOR POINT ( 2 , 1 )

RHS_P_DGESV ( 1 ) = 0.0D0

!*** FOR POINT ( 1 , N ) OF THE DOMAIN

B_DGESV ( N , N ) = 1.0D0 ! FOR POINT ( 1 , N )

B_DGESV ( N , N - 1 ) = -0.5D0 ! FOR POINT ( 1 , N - 1 )

B_DGESV ( N , N + N ) = -0.5D0 ! FOR POINT ( 2 , N )

RHS_P_DGESV ( 51 ) = 0.0D0

!*** FOR POINT ( N , 1 ) OF THE DOMAIN

B_DGESV ( (N-1)*N + 1 , (N-1)*N + 1 ) = 1.0D0 ! FOR POINT ( N , 1 )

B_DGESV ( (N-1)*N + 1 , (N-1)*N + 2 ) = -0.5D0 ! FOR POINT ( N , 2 )

B_DGESV ( (N-1)*N + 1 , (N-1)*N + 1 - N ) = -0.5D0 ! FOR POINT ( N - 1 , 1 )

RHS_P_DGESV ( (N-1)*N + 1 ) = 0.0D0

!*** FOR POINT ( N , N ) OF THE DOMAIN

B_DGESV ( N*N , N*N ) = 1.0D0 ! FOR POINT ( N , N )

B_DGESV ( N*N , N*N - 1 ) = -0.5D0 ! FOR POINT ( N , N - 1 )

B_DGESV ( N*N , N*N - N ) = -0.5D0 ! FOR POINT ( N - 1 , N )

RHS_P_DGESV ( N*N ) = 0.0D0

! ****************************** TO MAKE PART 3 FOR PARDISO

! ********** FOR POINT ( 1 , 1 ) OF DOMAIN

RHS_P ( 1 ) = 0.0D0

! ********** FOR POINT ( 1 , N ) OF DOMAIN

RHS_P ( 51 ) = 0.0D0

! ********** FOR POINT ( N , 1 ) OF DOMAIN

RHS_P ( (N-1)*N + 1 ) = 0.0D0

! ********** FOR POINT ( N , N ) OF DOMAIN

RHS_P ( N*N ) = 0.0D0

!******************** INITIALIZING THE VARIABLES WHICH NEEDED FOR CALLING PARDISO (FOR PRESSURES)

IPARM_B ( 1 ) = 1

IPARM_B ( 2 ) = 3

IPARM_B ( 3 ) = 1

IPARM_B ( 4 ) = 0

IPARM_B ( 5 ) = 0

IPARM_B ( 6 ) = 0

IPARM_B ( 7 ) = 0

IPARM_B ( 8 ) = -1000 ! 10

IPARM_B ( 9 ) = 0

IPARM_B ( 10 ) = 13

IPARM_B ( 11 ) = 0

IPARM_B ( 12 ) = 0

IPARM_B ( 13 ) = 0

IPARM_B ( 14 ) = 0

IPARM_B ( 15 ) = 0

IPARM_B ( 16 ) = 0

IPARM_B ( 17 ) = 0

IPARM_B ( 18 ) = 0

IPARM_B ( 19 ) = 0

IPARM_B ( 20 ) = -1 ! 0

IPARM_B ( 21 ) = 0

IPARM_B ( 22 ) = 0

IPARM_B ( 23 ) = 0

IPARM_B ( 24 ) = 1 ! 0

IPARM_B ( 25 ) = 0

IPARM_B ( 26 ) = 0

IPARM_B ( 27 ) = 1

IPARM_B ( 28 ) = 0

IPARM_B ( 29 ) = 0

IPARM_B ( 30 ) = 0

IPARM_B ( 31 ) = 0

IPARM_B ( 32 ) = 0

IPARM_B ( 33 ) = 0

IPARM_B ( 34 ) = 0

IPARM_B ( 35 ) = 0

IPARM_B ( 36 ) = 0

IPARM_B ( 37 ) = 0

IPARM_B ( 38 ) = 0

IPARM_B ( 39 ) = 0

IPARM_B ( 40 ) = 0

IPARM_B ( 41 ) = 0

IPARM_B ( 42 ) = 0

IPARM_B ( 43 ) = 0

IPARM_B ( 44 ) = 0

IPARM_B ( 45 ) = 0

IPARM_B ( 46 ) = 0

IPARM_B ( 47 ) = 0

IPARM_B ( 48 ) = 0

IPARM_B ( 49 ) = 0

IPARM_B ( 50 ) = 0

IPARM_B ( 51 ) = 0

IPARM_B ( 52 ) = 0

IPARM_B ( 53 ) = 0

IPARM_B ( 54 ) = 0

IPARM_B ( 55 ) = 0

IPARM_B ( 56 ) = 0

IPARM_B ( 57 ) = 0

IPARM_B ( 58 ) = 0

IPARM_B ( 59 ) = 0

IPARM_B ( 60 ) = 0

IPARM_B ( 61 ) = 0

IPARM_B ( 62 ) = 0

IPARM_B ( 63 ) = 0

IPARM_B ( 64 ) = 0

! ****************************** TO MAKE THE RIGHT HAND SIDE MATRIX "B_A" FOR PARDISO (FOR PRESSURES)

DO I = 1 , NN

B_B ( I ) = RHS_P ( I )

END DO

! ****************************** TO SOLVE THE SYSTEM OF EQUATIONS TO OBTAIN FINAL PRESSURES

INFO_DDNSCSR = 0

CALL MKL_DDNSCSR ( JOB_B_DDNSCSR , NN , NN , B_DGESV , NN , B_CSR_DDNSCSR , J_B_DDNSCSR , I_B_DDNSCSR , INFO_DDNSCSR )

INFO_DGESV = 0

CALL DGESV ( NN , 1 , B_DGESV , NN , IPIV_DGESV , RHS_P_DGESV , NN , INFO_DGESV )

DO I = 1 , N

DO J = 1 , N

P_N1_DGESV ( I , J ) = RHS_P_DGESV ( N*(I-1) + J )

END DO

!-----------------------------------------------------------------

! Only Analysis, Reordering And Symbolic Factorization, This Step Also Allocates All Memory That Is Necessary For The Factorization :

PHASE_B = 11

CALL PARDISO ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR_DDNSCSR, I_B_DDNSCSR, J_B_DDNSCSR, PERM_B, NRHS_B, IPARM_B, MSGLVL_B, B_B, X_B, ERROR_B)

! CALL PARDISO ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR, RowINDEX_B, COLUMN_B, PERM_B, NRHS_B, IPARM_B, MSGLVL_B, B_B, X_B, ERROR_B)

IF (ERROR_B .NE. 0) THEN

WRITE ( * , * ) 'The following ERROR was detected in Analysis, Reordering and Symbolic Factorization for solving VELOCITIES at ITERATION : ', ERROR_B , ITERATION

READ ( * , * )

END IF

!-----------------------------------------------------------------

! Only Numerical Factorization :

PHASE_B = 22

CALL PARDISO ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR_DDNSCSR, I_B_DDNSCSR, J_B_DDNSCSR, PERM_B, NRHS_B, IPARM_B, MSGLVL_B, B_B, X_B, ERROR_B)

! CALL PARDISO ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR, RowINDEX_B, COLUMN_B, PERM_B, NRHS_B, IPARM_B, MSGLVL_B, B_B, X_B, ERROR_B)

IF (ERROR_B .NE. 0) THEN

WRITE ( * , * ) 'The following ERROR was detected in Numerical Factorization for solving VELOCITIES at ITERATION : ', ERROR_B , ITERATION

READ ( * , * )

END IF

!-----------------------------------------------------------------

! Solving With Back Substitution And Iterative Refinement :

PHASE_B = 33

CALL PARDISO ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR_DDNSCSR, I_B_DDNSCSR, J_B_DDNSCSR, PERM_B, NRHS_B, IPARM_B, MSGLVL_B, B_B, X_B, ERROR_B)

! CALL PARDISO ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR, RowINDEX_B, COLUMN_B, PERM_B, NRHS_B, IPARM_B, MSGLVL_B, B_B, X_B, ERROR_B)

!-----------------------------------------------------------------

! Termination And Release Internal Memory :

PHASE_B = -1

CALL PARDISO ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR_DDNSCSR, I_B_DDNSCSR, J_B_DDNSCSR, PERM_B, NRHS_B, IPARM_B, MSGLVL_B, B_B, X_B, ERROR_B)

! CALL PARDISO ( PT_B, MAXFCT_B, MNUM_B, MTYPE_B, PHASE_B, N_B, B_CSR, RowINDEX_B, COLUMN_B, PERM_B, NRHS_B, IPARM_B, MSGLVL_B, B_B, X_B, ERROR_B)

!***********************************************************************************************************************************************

DO I = 1 , N

DO J = 1 , N

P_N1 ( I , J ) = X_B ( (I-1)*N + J )

END DO

! ******************** TO CALCULATE "GRADIENT P" AND THEN CALCULATE "U_N1" AND "V_N1" AT THE END OF EACH ITERATION

DO I = 2 , N-1

DO J = 2 , N-1

U_N1_DGESV (I,J) = U_INT_DGESV (I,J) - ( DELTA_T * ( (P_N1_DGESV (I,J+1) - P_N1_DGESV (I,J-1)) / (2*DELTA_X) ) )

V_N1_DGESV (I,J) = V_INT_DGESV (I,J) - ( DELTA_T * ( (P_N1_DGESV (I-1,J) - P_N1_DGESV (I+1,J)) / (2*DELTA_Y) ) )

END DO

! ******************** FOR PARDISO

DO I = 2 , N-1

DO J = 2 , N-1

U_N1 (I,J) = U_INT (I,J) - ( DELTA_T * ( (P_N1 (I,J+1) - P_N1 (I,J-1)) / (2*DELTA_X) ) )

V_N1 (I,J) = V_INT (I,J) - ( DELTA_T * ( (P_N1 (I-1,J) - P_N1 (I+1,J)) / (2*DELTA_Y) ) )

END DO

! ******************** ENFORCE THE DIRICHLET BOUNDARY CONDITION OF ON THE LEFT SIDE OF THE DOMAIN

DO I = 1 , N

U_N1_DGESV ( I , 1 ) = 0

V_N1_DGESV ( I , 1 ) = 0

END DO

U_N1_DGESV ( N2-1 , 1 ) = 0.75D0 * 2.0D0

U_N1_DGESV ( N2 , 1 ) = 2.0D0

U_N1_DGESV ( N2+1 , 1 ) = 0.75D0 * 2.0D0

DO I = 1 , N

DO J = 1 , N

U_N_DGESV (I,J) = U_N1_DGESV (I,J)

V_N_DGESV (I,J) = V_N1_DGESV (I,J)

END DO

! ******************** FOR PARDISO

DO I = 1 , N

U_N1 ( I , 1 ) = 0

V_N1 ( I , 1 ) = 0

END DO

U_N1 ( N2 - 1 , 1 ) = 0.75D0 * 2.0D0

U_N1 ( N2 , 1 ) = 2.0D0

U_N1 ( N2 + 1 , 1 ) = 0.75D0 * 2.0D0

DO I = 1 , N

DO J = 1 , N

U_N (I,J) = U_N1 (I,J)

V_N (I,J) = V_N1 (I,J)

END DO

! ******************** EXPORTING RESULTS FOR TECHPLOT

XITERATION = ITERATION

IF ( ITERATION / IPRINT == XITERATION / IPRINT ) THEN

FORMAT_STRING = "(A14,I6)"

WRITE ( FILE_NAME_1 , FORMAT_STRING ) "PARDISO" , ITERATION

WRITE ( FILE_NAME_2 , FORMAT_STRING ) "DGESV" , ITERATION

FINAL_NAME_1 = TRIM(FILE_NAME_1)//".TXT"

FINAL_NAME_2 = TRIM(FILE_NAME_2)//".TXT"

OPEN ( 1 , FILE = FINAL_NAME_1 )

101 FORMAT (7F30.16)

DO I = N , 1 , -1

DO J = 1 , N

WRITE ( 1 , 101 ) DELTA_X * (J-1) , DELTA_Y * (N-I) , U_INT (I,J) , V_INT (I,J) , U_N1 (I,J) , V_N1 (I,J) , P_N1 (I,J)

END DO

CLOSE ( 1 )

OPEN ( 2 , FILE = FINAL_NAME_2 )

DO I = N , 1 , -1

DO J = 1 , N

WRITE ( 2 , 101 ) DELTA_X * (J-1) , DELTA_Y * (N-I) , U_INT_DGESV (I,J) , V_INT_DGESV (I,J) , U_N1_DGESV (I,J) , V_N1_DGESV (I,J) , P_N1_DGESV (I,J)

END DO

CLOSE ( 2 )

END IF

! ********************************************************************************

END DO

! ********************************************************************************

END PROGRAM DGESV_VS_DDNSCSR_PARDISO

Abdolreza_A_ — Fri, 11 Aug 2017 09:17:11 GMT

PROGRAM DGESV_VS_DDNSCSR_PARDISO

IMPLICIT NONE

! INCLUDE 'MKL_DSS.f90' ! Include the standard DSS "header file."

USE MKL_DSS