Community
cancel
Showing results for 
Search instead for 
Did you mean: 
Abdolreza_A_
Beginner
64 Views

problem with PARDISO and DSS

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.

0 Kudos
6 Replies
Abdolreza_A_
Beginner
64 Views

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
                 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
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
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
                 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
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
                 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
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
                 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
                 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
                 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
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
                 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
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
                 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
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
    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
    END DO
    CLOSE ( 2 )
 
    
END IF
 
! ********************************************************************************
END DO
! ********************************************************************************
    END PROGRAM DGESV_VS_DDNSCSR_PARDISO
Abdolreza_A_
Beginner
64 Views

 
PROGRAM DGESV_VS_DDNSCSR_PARDISO
IMPLICIT NONE
! INCLUDE 'MKL_DSS.f90' ! Include the standard DSS "header file."
USE MKL_DSS
 
! ****************************************************************************************
 
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
 
! ***********************************************************************************************
INTEGER :: ERROR
TYPE ( MKL_DSS_HANDLE ) :: HANDLE ! Allocate storage for the solver handle.
INTEGER OERM (1)
! ****************************** 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
                 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
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
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
                 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
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
                 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
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
                 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
                 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
                 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)
 
!***********************************************************************************************************************************************
PERM (1) = 0
 
! INITIALIZE THE SOLVER :
ERROR = DSS_CREATE ( HANDLE , MKL_DSS_DEFAULTS )
IF ( ERROR /= MKL_DSS_SUCCESS ) GOTO 999
 
! DEFINE THE NON-ZERO STRUCTURE OF THE MATRIX :
ERROR = DSS_DEFINE_STRUCTURE ( HANDLE , MKL_DSS_NON_SYMMETRIC , I_B_DDNSCSR , NN , NN , J_B_DDNSCSR , NON_ZERO_OF_B )
IF ( ERROR /= MKL_DSS_SUCCESS ) GOTO 999
 
 
! REORDER THE MATRIX ;
ERROR = DSS_REORDER ( HANDLE , MKL_DSS_DEFAULTS , PERM )
IF ( ERROR /= MKL_DSS_SUCCESS ) GOTO 999
 
! FACTOR THE MATRIX :
ERROR = DSS_FACTOR_REAL ( HANDLE , MKL_DSS_DEFAULTS , B_CSR_DDNSCSR )
IF ( ERROR /= MKL_DSS_SUCCESS ) GOTO 999
 
 
! SOLVE THE PROBLEM :
ERROR = DSS_SOLVE_REAL ( HANDLE , MKL_DSS_DEFAULTS , B_B , 1 , X_B )
IF ( ERROR /= MKL_DSS_SUCCESS ) GOTO 999
 
! DEALLOCATE SOLVER STORAGE AND VARIOUS LOCAL ARRAYS :
ERROR = DSS_DELETE ( HANDLE , MKL_DSS_DEFAULTS )
 
! PRINT THE ERROR MESSAGE :
999 WRITE ( * , * ) "SOLVER RETURNED THE ERROR CODE :", ERROR
    READ  ( * , * )
 
!***********************************************************************************************************************************************
 
DO I = 1 , N
    DO J = 1 , N
        P_N1 ( I , J ) = X_B ( (I-1)*N + J )
    END DO
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
                 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
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
                 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
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
    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
    END DO
    CLOSE ( 2 )
 
    
END IF
 
! ********************************************************************************
END DO
! ********************************************************************************
    END PROGRAM DGESV_VS_DDNSCSR_PARDISO
mecej4
Black Belt
64 Views

There are several deficiencies in your post, which make it next to impossible to give you a useful response. Here are some suggestions to improve the help request.

1. Simplify the program to just what is needed to display one problem at a time. For example, make one pair of programs that solve just for the u-velocity with a small value of N, say as 5. In the first program, use DGESV and, in the second, use MKL_DDNSCSR + Pardiso. Do the two results agree to, say, seven digits? If so, repeat with larger values of N until you see a pair of solutions that are different more than you find tolerable, and report just that case.

2. Similarly, make a second pair of programs in which you solve just for the u-velocity using (i) your own assembly of the CSR representation and (ii) MKL_DDNSCSR. Do the two index arrays IA,JA, the values array V and the RHS array agree between the two versions? If they do not, report just that set of results and the programs that we can use to reproduce the problem.

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 returned 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 do not find that believable. A problem description is passed to Pardiso solely through the arguments in the CALL PARDISO statement. If two identical sequences of calls with the same values of the arguments were passed and the results differed, that would be a serious bug in Pardiso and many other users would have noticed and reported such a bug. I have used Pardiso extensively for over a decade, and have never encountered such a bug in it.

The DSS calls in MKL are provided only to help users of other compilers such as Compaq/DEC Visual Fortran and Sun Fortran which are provided with a library (CXML for Compaq VF; Sun Perflib for Sun Fortran) that contained DSS. The DSS calls in MKL are just wrappers that convert DSS arguments to the corresponding Pardiso arguments and call Pardiso. For this reason, the DSS interface is a portability feature that should be of no interest to you.

Abdolreza_A_
Beginner
64 Views

Hi dear mecej4 I do not find that believable too ! That is exactly why I have posted this problem ! I compared several times the matrices which I made myself to the matrices which were made by DNSCSR. all of them are exactly the same ! So I send my code to believe my words ! It just need copy and paste into the fortran and run (it doesn't need any import file). accept my appolofy for wasting your time. MY MATRICES >> DGESV returns the excellent results. MY MATRICES >> DNSCSR >> PARDISO returns almost the same results (just little differencies) MY OWN CSR MATRICES (WHICH DOESN'T ANY DIFFERENCES WITH WHICH HAVE MADE BY "DNSCSR") >> PARDISO : doesn't return the same results.
Abdolreza_A_
Beginner
64 Views

Morover the code will print 2 files at each iteration : THE FIRST AT EACH ITERATION : the result of DGESV THE SECOND AT EACH ITERATION : the result of PARDISO (in the first code which I have posted, it is the result of MY GENERAL MATRICES >> DDNSCSR >> PARDISO 's results **** in the second code which I have posted : the result of MY OWN CSR >> directly to PARDSIO 's results)
mecej4
Black Belt
64 Views

The second code that you posted (in #3) has syntax errors and cannot be compiled:

     The statement USE MKL_DSS is out of order. The PERM argument that you use with DSS_REORDER is not declared. 

If you actually ran the code, as you said, the version that you posted is different from the version that you were able to compile. That puts a stop to any discussion of results and their differences.

Please do not post code longer than, say, 100 lines in the body of a message. Zip the source file(s) and attach the zip to your post.

As I said in #4, it would help if your code did just enough to display a problem. I see no part of your codes where A_CSR is computed, or even declared. I do not wish to plow through long codes that do lots of extraneous calculations that are of no relevance to the MKL problem being reported.

As a first step, please construct an example with a small value of N that does just one iteration and shows different solutions for the u-velocity with 1) your own code for computing the CSR representation and 2) MKL_DDNSCSR.

Reply