- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
Link Copied
6 Replies
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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)
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page