Community
cancel
Showing results for
Did you mean: Beginner
243 Views

## DSS symmetric solve

We have been using the following code compiled with Intell Fortran on Linux clusters linked with MKL 10.2.2.025 for a number of years.  We are trying to update to Intel 15 compilers and MKL 11.2 update 2.  The symmetric solve produces wrong answers using the new compiler and libraries.  The code has not changed.  When we include the symmetric terms and solve it as a symmetric_structure matrix we get the right answers.

Is any one else having problems with symmetric solves with the current compiler and library?

Load the data into indx, vall, and y_rt ...

lcount = 0
do row=1,n_solver

pntr(row) = lcount +1

col = row
lcount = lcount+1
indx(lcount) = col

col = row + 1
!         if( (col <= n*m).and.(au1(row) /= 0.)) then
if( (col <= n*m) ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au1(row)
endif

col = row - 1 + m
!         if( (col <= n*m).and.(au2(row) /= 0.)) then
if( (col <= n*m) ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au2(row)
endif

col = row + m
!         if( (col <= n*m).and.(au3(row) /= 0.)) then
if( (col <= n*m) ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au3(row)
endif

col = row + 1 + m
!         if( (col <= n*m).and.(au4(row) /= 0.)) then
if( (col <= n*m) ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au4(row)
endif

end do

! Take care of last pntr
pntr(n_solver+1) = lcount + 1
non_zeros=lcount

error = DSS_DEFINE_STRUCTURE( handle, MKL_DSS_SYMMETRIC, pntr, n_solver, &
n_solver, indx, non_zeros )

error = DSS_REORDER( handle, MKL_DSS_DEFAULTS, perm )

error = DSS_FACTOR_REAL( handle, MKL_DSS_INDEFINITE, vall )

nRhs=1 !number of right hand sides
error = DSS_SOLVE_REAL(handle, MKL_DSS_DEFAULTS, y_rt, nRhs, x_rt )

x_rt is correct under Intel 11/MKL 10.2.2 but wrong with intel 15/MKL 11.4 update 2

15 Replies Beginner
243 Views

Update:  It also returns a correct solution for Intel 13.1.117.  So the problem is occurring somewhere between the Intel 13 and Intel 15 releases.

I also forgot to include the compilation and link options.  We do link in FFTW and which does have some MKL interfaces in it which may overlap the MKL library.

mpiifort -fpe0 -D_UNIX -r8 -O0 -I /theory/kellerd/shu_test/Draco_i2_5_iSNB_XBT_FPEOS_Blizzard/finclude -fpp -convert big_endian -I/opt/LLE/local/fftw-3.3.1/include -I/opt/LLE/local/fftw-3.3.1/api -I/opt/LLE/local/hdf5-1.8.2/include -I/opt/LLE/local/hdf5-1.8.2/lib -c

mpiifort -o draco_hp15 hp15/*.o -fpe0 -D_UNIX  -Wl,--start-group -L/opt/LLE/local/fftw-3.3.1/lib -lfftw3  -liomp5 -lpthread -lz -lmpiif -lmpi -lmkl_core -lmkl_sequential -lmkl_intel_thread -lmkl_intel_lp64 -L/opt/LLE/local/hdf5-1.8.2/lib -lhdf5 -L/opt/LLE/local/hdf5-1.8.2/lib -lhdf5_fortran -Wl,--end-group Black Belt
243 Views

The DSS example code provided with MKL works correctly (checked with Matlab) with IFort 15.0.2 and MKL 11.2. I suggest that you try it on your system, and if it runs correctly, provide a reproducer based on your code that enables the problem to be seen. Moderator
243 Views

yes, please check the problem with the latest version of MKL 11.2 update 2 and then in the case if you will see the problem is still exists, give us the example of this code which we can compile and debug this problem. You can take the eval version of MKL for this evaluation. Beginner
243 Views

I developed a small test that demonstrates the error.  It creates the correct solution vector under Intel 13 and the incorrect solution under Intel 15 with MKL 11.2 update 2.  The input file is much smaller than our number array size but still is  623862720 bytes and attached.

>cat run
#!/bin/sh
# can be more optimized on newer architectures
ldd a.out
./a.out

>cat out.13

linux-vdso.so.1 =>  (0x00007fff5accc000)
libdl.so.2 => /lib64/libdl.so.2 (0x0000003f08000000)
libmkl_intel_lp64.so => /opt/lic/intel13/composer_xe_2013.1.117/mkl/lib/intel64/libmkl_intel_lp64.so (0x00002b60b33f4000)
libmkl_core.so => /opt/lic/intel13/composer_xe_2013.1.117/mkl/lib/intel64/libmkl_core.so (0x00002b60b4a99000)
libiomp5.so => /opt/lic/intel13/composer_xe_2013.1.117/compiler/lib/intel64/libiomp5.so (0x00002b60b5ca8000)
libmpi.so.4 => /opt/lic/intel13/impi/4.1.0.024/intel64/lib/libmpi.so.4 (0x00002b60b5faa000)
libmpigf.so.4 => /opt/lic/intel13/impi/4.1.0.024/intel64/lib/libmpigf.so.4 (0x00002b60b65d0000)
librt.so.1 => /lib64/librt.so.1 (0x0000003f08800000)
libm.so.6 => /lib64/libm.so.6 (0x0000003f07800000)
libc.so.6 => /lib64/libc.so.6 (0x0000003f07c00000)
libgcc_s.so.1 => /lib64/libgcc_s.so.1 (0x0000003f11800000)
/lib64/ld-linux-x86-64.so.2 (0x0000003f07400000)
start solving the system
TEST MKL DSS; Execution time=  0.2031250
50156093319.1081       -89969204211.7639        89514538131.7439     <<<<<<<< correct numbers

>cat out.15
linux-vdso.so.1 =>  (0x00007fffbf1ff000)
libmkl_intel_lp64.so => /opt/lic/intel15/composer_xe_2015.2.164/mkl/lib/intel64/libmkl_intel_lp64.so (0x00002b2d13f89000)
libmkl_core.so => /opt/lic/intel15/composer_xe_2015.2.164/mkl/lib/intel64/libmkl_core.so (0x00002b2d15cbd000)
libiomp5.so => /opt/lic/intel15/composer_xe_2015.2.164/compiler/lib/intel64/libiomp5.so (0x00002b2d1781c000)
libmpifort.so.12 => /opt/lic/intel15/impi/5.0.3.048/intel64/lib/libmpifort.so.12 (0x00002b2d17b6d000)
libmpi.so.12 => /opt/lic/intel15/impi/5.0.3.048/intel64/lib/libmpi.so.12 (0x00002b2d17df7000)
libdl.so.2 => /lib64/libdl.so.2 (0x0000003f08000000)
librt.so.1 => /lib64/librt.so.1 (0x0000003f08800000)
libm.so.6 => /lib64/libm.so.6 (0x0000003f07800000)
libc.so.6 => /lib64/libc.so.6 (0x0000003f07c00000)
libgcc_s.so.1 => /lib64/libgcc_s.so.1 (0x0000003f11800000)
/lib64/ld-linux-x86-64.so.2 (0x0000003f07400000)
start solving the system
TEST MKL DSS; Execution time=  0.1875000
73375.2577664697       -359933.504809541        280913.663845399    <<<<<< incorrect numbers

```INCLUDE 'mkl_dss.f90' ! Include the standard DSS "header file."
use mkl_dss
use ifport
implicit none
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
integer, parameter :: DP = KIND(1.0D0)
integer, parameter        :: MAX_NZ_ROW = 9
logical, save:: initialized=.FALSE.
integer :: m_solver,n_solver
integer :: n, m, row, col, lcount, non_zeros, error, perm(1),nRhs
type(MKL_DSS_HANDLE),save :: handle ! Allocate storage for the solver handle.

integer, allocatable, dimension(:), save:: indx, pntr
real(KIND=DP),   allocatable, dimension(:), save:: vall,x_rt,y_rt
!real(8) :: t0, t1, MPI_WTIME
real(4) :: t0, t1, ta(2)
character*(10) :: suffix
integer :: ierr
real(8) :: x0
integer(4) :: ic0,ic1,crate
!-----------------------------------------------

!      call MPI_INIT(ierr)
!DSS implementation
!  val   -  values
!  pntr  -  rowIndex
!  indx  -  colums
!  ncols =  n_solver
!  nrows =  ncols

m_solver=60   !400  !400      !120 (for test)   !=je-js+1
n_solver=2280 !200  !2000     !498 (for test)   !=ie-is+1
n_solver=n_solver*m_solver

m = m_solver
n = n_solver/m_solver

if(.NOT. initialized)then
if( .NOT. allocated(indx) ) allocate(indx(MAX_NZ_ROW*n_solver))
if( .NOT. allocated(pntr) ) allocate(pntr(n_solver+1))
if( .NOT. allocated(vall) ) allocate(vall(MAX_NZ_ROW*n_solver))
if( .NOT. allocated(x_rt) ) allocate(x_rt(n_solver))
if( .NOT. allocated(y_rt) ) allocate(y_rt(n_solver))
endif
indx(:)=0
pntr(:)=0
vall(:)=0.0d0
x_rt(:)=0.0d0
y_rt(:)=0.0d0

lcount = 0
do row=1,n_solver

pntr(row) = lcount +1

col = row
lcount = lcount+1
indx(lcount) = col

col = row + 1
if( col <= n*m ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au1
endif

col = row - 1 + m
if( col <= n*m ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au2
endif

col = row + m
if( col <= n*m ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au3
endif

col = row + 1 + m
if( col <= n*m ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au4
endif

end do
!
!     Take care of last pntr
!
pntr(n_solver+1) = lcount + 1
non_zeros=lcount

if (.NOT. initialized) then

! 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_SYMMETRIC, pntr, n_solver, &
& n_solver, indx, non_zeros )
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

initialized = .TRUE.
endif

print*,"start solving the system"
call SYSTEM_CLOCK(ic0,crate)
t0 = dble(ic0)/dble(crate)
!
! Factor the matrix.
!
error = DSS_FACTOR_REAL( handle, MKL_DSS_INDEFINITE, vall )
if (error /= MKL_DSS_SUCCESS) goto 999

!
! Solve for a vector x_rt, given an input vector y_rt as b
!
nRhs=1 !number of right hand sides
error = DSS_SOLVE_REAL(handle, MKL_DSS_DEFAULTS, y_rt, nRhs, x_rt )
if (error /= MKL_DSS_SUCCESS) goto 999

call SYSTEM_CLOCK(ic1,crate)
t1 = dble(ic1)/dble(crate)
print*,"TEST MKL DSS; Execution time=",t1-t0
print*,sum(x_rt),minval(x_rt),maxval(x_rt)

open(566,FILE="x_rt_based_on_input_"//TRIM(suffix))
do row=1,n_solver
write(566,*) x_rt(row)
enddo
stop

999 WRITE(*,*) "Solver returned error code ", error Beginner
243 Views

Noticed my typo after I posted the example.  My note said "the input file is smaller than out number array size".  I meant to say it is smaller than our typical array sizes.  The attached input DOES demonstrate the problem Black Belt
243 Views

DSS lets you compute some statistics, among them the determinant. For your "small" matrix, I obtained this:

`Determinant =  2.407 X 10**( -1057550)`

Applying Pardiso to the same matrix, I obtained:

```non-zeros in A:           683819
number of nonzeros in L  5914460
number of nonzeros in U        1```

Something is seriously wrong if there is only one nonzero element in the U factor (I need to check if I have made an error in obtaining these results). In these circumstances, your summary digest of the solution (min, max and sum of elements of X) is not useful, probably not reproducible across platforms, and may even change with slightly different compiler options. There may be similar problems with your "large" matrices. Before accepting a solution, you should probably obtain the residuals after solving the equations. For the "small" matrix, I obtained residuals as large as 0.001.

Note that the data file containing the "small" matrix is large (> 600 MB) in formatted text form. By placing the data into an unformatted file and cutting out the portions that are not used in the calculation, you could reduce it to about 10 MB. Beginner
243 Views

Thanks for looking at this for us.

I admit to really not knowing much about solvers but the issue is that the solution to an input set returns a solution that my physcis folks say is 'right' on Intel 13 is 'wrong' on Intel 15.  It would seem that if the sum, min and max of the solution differ significantly from an Intel 13 run to an Intel 15 run something significant has changed and that our 'right' answers have significantly changed with Intel 15.  It is not clear to me from your note if you think there is something 'seriously wrong' with the input or coding or from the result?

I modified the source to include additional statistics:

print*,sum(x_rt),minval(x_rt),maxval(x_rt)

call MKL_CVT_TO_NULL_TERMINATED_STR(statreq,11,"determinant")
error = DSS_STATISTICS(handle,MKL_DSS_DEFAULTS,statreq,statistic)
print*,"Determinant=",statistic(2)," X 10 ",statistic(1)

call MKL_CVT_TO_NULL_TERMINATED_STR(statreq,7,"inertia")
error = DSS_STATISTICS(handle,MKL_DSS_DEFAULTS,statreq,statistic)
print*,"pos =",statistic(1)
print*,"neg =",statistic(2)
print*,"zer =",statistic(3)

call MKL_CVT_TO_NULL_TERMINATED_STR(statreq,8,"solvemem")
error = DSS_STATISTICS(handle,MKL_DSS_DEFAULTS,statreq,statistic)
print*,"Solvemem=",statistic(1)

Intel 13:

50204591221.3627       -89969204206.9475        89514538126.9583
Determinant=   1.75565552825196       X 10   -1336812.00000000
pos =   67894.0000000000
neg =   68906.0000000000
zer =  0.000000000000000E+000
Solvemem=   69469.0000000000

Intel 15

73375.2577673518       -359933.504809540        280913.663845398
Determinant=  -2.32535455034809       X 10   -1057541.00000000
pos =   65547.0000000000
neg =   71253.0000000000
zer =  0.000000000000000E+000
Solvemem=   69867.0000000000 Black Belt
243 Views

You can further the testing if you can put together a smaller example, with as few unknowns as possible, but with the same structure as the matrices that you have shown, and approved by the "physics folks" as correct with IFort 13/MKL 10.2.2.  The problem with large size matrix problems is that you are restricted to fewer checks on correctness. On the other hand, when you whittle a problem down you may be left with an absurd remnant.

I shall be happy to share my verification code with you once I have done some checks.

Note this from #8:

Determinant=   1.76  X 10   -1336812

It is possible for a matrix to be well-conditioned and yet have a determinant nearly equal to 0. However, if the problem comes from a physics model and the variables are reasonably scaled this should not happen. Ask the physics folks if they are comfortable with this attoscopic value for the determinant. If this value is wrong because of a coding error, we can track that down and fix it. However, if the problem as formulated is nearly impossible to obtain an accurate solution for, the physics people may have a much tougher problem.

A more basic question: is it expected/reasonable that the matrix is singular? If not, there can be no zeroes on the principal diagonals of the L and U factors. Therefore, the minimum number of non-zeroes in L and U is equal to the number of equations. Beginner
243 Views

The physics guys are not around today so I tried to do what I can to further help find this problem.  I cut the grid size down to a 14x30 grid so

m =8 n=112.  It still demostrates the problem.  The solution grid when plotted follows.  Note they show that the garray1 field on the left axis is "wrong" for 15, "right" for 13.  I also attached attached the input file (still in text format) for this smaller grid size. 451535.

cat out.15_30
linux-vdso.so.1 =>  (0x00007fff86871000)
libmkl_intel_lp64.so => /opt/lic/intel15/composer_xe_2015.2.164/mkl/lib/intel64/libmkl_intel_lp64.so (0x00002b41f3e18000)
libmkl_core.so => /opt/lic/intel15/composer_xe_2015.2.164/mkl/lib/intel64/libmkl_core.so (0x00002b41f5b4c000)
libiomp5.so => /opt/lic/intel15/composer_xe_2015.2.164/compiler/lib/intel64/libiomp5.so (0x00002b41f76ab000)
libmpifort.so.12 => /opt/lic/intel15/impi/5.0.3.048/intel64/lib/libmpifort.so.12 (0x00002b41f79fc000)
libmpi.so.12 => /opt/lic/intel15/impi/5.0.3.048/intel64/lib/libmpi.so.12 (0x00002b41f7c86000)
libdl.so.2 => /lib64/libdl.so.2 (0x0000003f08000000)
librt.so.1 => /lib64/librt.so.1 (0x0000003f08800000)
libm.so.6 => /lib64/libm.so.6 (0x0000003f07800000)
libc.so.6 => /lib64/libc.so.6 (0x0000003f07c00000)
libgcc_s.so.1 => /lib64/libgcc_s.so.1 (0x0000003f11800000)
/lib64/ld-linux-x86-64.so.2 (0x0000003f07400000)
start solving the system
TEST MKL DSS; Execution time=  7.8125000E-02
295618.772828955       -312701.448947056        220058.681258707
Determinant=  -1.85272122511410       X 10   -98200.0000000000
pos =   6727.00000000000
neg =   5873.00000000000
zer =  0.000000000000000E+000
Solvemem=   5647.00000000000

cat out.13_30
linux-vdso.so.1 =>  (0x00007fff7e4b4000)
libdl.so.2 => /lib64/libdl.so.2 (0x0000003f08000000)
libmkl_intel_lp64.so => /opt/lic/intel13/composer_xe_2013.1.117/mkl/lib/intel64/libmkl_intel_lp64.so (0x00002ac989820000)
libmkl_core.so => /opt/lic/intel13/composer_xe_2013.1.117/mkl/lib/intel64/libmkl_core.so (0x00002ac98aec5000)
libiomp5.so => /opt/lic/intel13/composer_xe_2013.1.117/compiler/lib/intel64/libiomp5.so (0x00002ac98c0d4000)
libmpi.so.4 => /opt/lic/intel13/impi/4.1.0.024/intel64/lib/libmpi.so.4 (0x00002ac98c3d6000)
libmpigf.so.4 => /opt/lic/intel13/impi/4.1.0.024/intel64/lib/libmpigf.so.4 (0x00002ac98c9fc000)
librt.so.1 => /lib64/librt.so.1 (0x0000003f08800000)
libm.so.6 => /lib64/libm.so.6 (0x0000003f07800000)
libc.so.6 => /lib64/libc.so.6 (0x0000003f07c00000)
libgcc_s.so.1 => /lib64/libgcc_s.so.1 (0x0000003f11800000)
/lib64/ld-linux-x86-64.so.2 (0x0000003f07400000)
start solving the system
TEST MKL DSS; Execution time=  1.5625000E-02
-19306486435.1262       -5772347418.44442        3505444963.56421
Determinant=  -3.76569356065273       X 10   -117523.000000000
pos =   6477.00000000000
neg =   6123.00000000000
zer =  0.000000000000000E+000
Solvemem=   5603.00000000000   Black Belt
243 Views

Do you have a writeup with the P.D.E.,boundary and initial conditions specified? Beginner
243 Views

I just recognized today that I had sent in incorrect data input file.  I apologize for any wasted effort. The net is still the same however.  Intel 15 returns wrong answers.  I added a double check after the solve.  The solution returned when used to recompute the y_rt supplied has 9 errors:

Intel 13 output:

start solving the system
TEST MKL DSS; Execution time=  0.1109009
0.257042531294272       9.162025370382090E-014  4.008343651777284E-003
Determinant=   2.12161563270491       X 10   -2990.00000000000
pos =   420.000000000000
neg =  0.000000000000000E+000
zer =  0.000000000000000E+000
Solvemem=   159.000000000000

Intel 15 output:

start solving the system
TEST MKL DSS; Execution time=  7.7514648E-02
0.251832820371750       9.162025183722943E-014  4.008343651777284E-003
Determinant=   1.62949566754441       X 10   -2987.00000000000
pos =   420.000000000000
neg =  0.000000000000000E+000
zer =  0.000000000000000E+000
Solvemem=   262.000000000000
not equal         121  2.405238875359285E-019  2.065808031766965E-019                         <<<<<< double check should at least be close
not equal         151  2.405238875359036E-019  2.068078684773001E-019
not equal         181  2.405238875359422E-019  2.038278901649021E-019
not equal         211  2.405238875359251E-019  1.844924519720349E-019
not equal         241  2.405238875359413E-019 -6.431842696536678E-020
not equal         331  3.299385792934408E-020 -7.117734061904493E-014
not equal         332  9.889113991683395E-020 -6.104485185954297E-014
not equal         361  3.303422135645092E-020 -1.172086994261061E-013
not equal         391  8.509130524910173E-013  6.084941155961990E-013

A slightly modified source file and a new input file are provided.  They clearly demonstrate a problem and should be easily reproduced

```INCLUDE 'mkl_dss.f90' ! Include the standard DSS "header file."
use mkl_dss
use ifport
implicit none
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
integer, parameter :: DP = 8 !KIND(1.0D0)
integer, parameter        :: MAX_NZ_ROW = 9
logical, save:: initialized=.FALSE.
integer :: m_solver,n_solver
integer :: n, m, row, col, lcount, non_zeros, error, nRhs
type(MKL_DSS_HANDLE),save :: handle ! Allocate storage for the solver handle.

integer, allocatable, dimension(:), save:: indx, pntr, perm
real(KIND=DP),   allocatable, dimension(:), save:: vall,x_rt,y_rt,y_rt_check
!real(8) :: t0, t1, MPI_WTIME
real(4) :: t0, t1, ta(2)
character*(10) :: suffix
integer :: ierr,statreq(10)
real(8) :: x0, statistic(10)
integer(4) :: ic0,ic1,crate

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

!      call MPI_INIT(ierr)
!DSS implementation
!  val   -  values
!  pntr  -  rowIndex
!  indx  -  colums
!  ncols =  n_solver
!  nrows =  ncols

m= 30  !60   !400  !400      !120 (for test)   !=je-js+1
n= 14  !2280 !200  !2000     !498 (for test)   !=ie-is+1

n_solver = n * m

if(.NOT. initialized)then
if( .NOT. allocated(indx) ) allocate(indx(MAX_NZ_ROW*n_solver))
if( .NOT. allocated(pntr) ) allocate(pntr(n_solver+1))
if( .NOT. allocated(vall) ) allocate(vall(MAX_NZ_ROW*n_solver))
if( .NOT. allocated(x_rt) ) allocate(x_rt(n_solver))
if( .NOT. allocated(y_rt) ) allocate(y_rt(n_solver))
if( .NOT. allocated(y_rt_check) ) allocate(y_rt_check(n_solver))
if( .NOT. allocated(perm)) allocate(perm(n_solver))
endif
indx(:)=0
pntr(:)=0
vall(:)=0.0d0
x_rt(:)=0.0d0
y_rt(:)=0.0d0

open(555,FILE="input_30x14")

lcount = 0
do row=1,n_solver

pntr(row) = lcount +1

col = row
lcount = lcount+1
indx(lcount) = col

col = row + 1
if( col <= n*m ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au1
endif

col = row - 1 + m
if( col <= n*m ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au2
endif

col = row + m
if( col <= n*m ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au3
endif

col = row + 1 + m
if( col <= n*m ) then
lcount = lcount+1
indx(lcount) = col
vall(lcount) = au4
endif

end do
!
!     Take care of last pntr
!
pntr(n_solver+1) = lcount + 1
non_zeros=lcount

if (.NOT. initialized) then

! 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_SYMMETRIC, pntr, n_solver, &
& n_solver, indx, non_zeros )
if (error /= MKL_DSS_SUCCESS) goto 999

! Reorder the matrix.
!  error = DSS_REORDER( handle, MKL_DSS_DEFAULTS, perm )
error = DSS_REORDER( handle, MKL_DSS_AUTO_ORDER, perm )
if (error /= MKL_DSS_SUCCESS) goto 999

initialized = .TRUE.
endif

print*,"start solving the system"
call SYSTEM_CLOCK(ic0,crate)
t0 = dble(ic0)/dble(crate)
!
! Factor the matrix.
!
error = DSS_FACTOR_REAL( handle, MKL_DSS_INDEFINITE, vall )
if (error /= MKL_DSS_SUCCESS) goto 999

!
! Solve for a vector x_rt, given an input vector y_rt as b
!
nRhs=1 !number of right hand sides
error = DSS_SOLVE_REAL(handle, MKL_DSS_DEFAULTS, y_rt, nRhs, x_rt )
if (error /= MKL_DSS_SUCCESS) goto 999

call SYSTEM_CLOCK(ic1,crate)
t1 = dble(ic1)/dble(crate)
print*,"TEST MKL DSS; Execution time=",t1-t0
print*,sum(x_rt),minval(x_rt),maxval(x_rt)

call MKL_CVT_TO_NULL_TERMINATED_STR(statreq,11,"determinant")
error = DSS_STATISTICS(handle,MKL_DSS_DEFAULTS,statreq,statistic)
print*,"Determinant=",statistic(2)," X 10 ",statistic(1)

call MKL_CVT_TO_NULL_TERMINATED_STR(statreq,7,"inertia")
error = DSS_STATISTICS(handle,MKL_DSS_DEFAULTS,statreq,statistic)
print*,"pos =",statistic(1)
print*,"neg =",statistic(2)
print*,"zer =",statistic(3)

call MKL_CVT_TO_NULL_TERMINATED_STR(statreq,8,"solvemem")
error = DSS_STATISTICS(handle,MKL_DSS_DEFAULTS,statreq,statistic)
print*,"Solvemem=",statistic(1)

!then calculate the corresponding right-hand-side y_rt by multipying x_rt by the matrix storred in vall() array
y_rt_check(:)=0.0d0
do row=1,n_solver
do lcount=pntr(row),pntr(row+1)-1
col = indx(lcount)
y_rt_check(row) = y_rt_check(row) + vall(lcount)*x_rt(col)
if(col/=row) y_rt_check(col) = y_rt_check(col) + vall(lcount)*x_rt(row)
enddo
enddo

open(566,FILE="x_rt_based_on_input_"//TRIM(suffix))
do row=1,n_solver
write(566,*) x_rt(row)
if (abs(y_rt(row)-y_rt_check(row)) > abs(1.0e-6*y_rt(row))) print*,'not equal',row,y_rt(row), y_rt_check(row)
enddo
stop

999 WRITE(*,*) "Solver returned error code ", error
```

.

451632 Black Belt
243 Views

I do not have the specific older version, but I think that you will reproduce the old results with the new compiler if, on line 133 of #5, you replace  MKL_DSS_INDEFINITE by MKL_DSS_DEFAULTS which, according to the current documentation, sets MKL_DLL_POSITIVE_DEFINITE.

Of course, you have to check whether this is a correct setting. You can also check your documentation for the old and new versions to see if the default value for the call to DSS_FACTOR changed with the version.

This may also indicate that, as you suspected, there is an error in the newer version of MKL. This is likely since the DSS calls are mapped into equivalent PARDISO calls, and an error in doing so may have been introduced recently. You may wish to consider rewriting the solver part of your code to make direct PARDISO calls instead of going through the DSS interface. Beginner
243 Views

Looks like an error has been introduced.

Replacing MKL_DSS_INDEFINITE by MKL_DSS_DEFAULTS or MKL_DLL_POSITIVE_DEFINITE with Intel 15 results in:

start solving the system
TEST MKL DSS; Execution time=  6.2500000E-02
0.257042531294272       9.162025370385757E-014  4.008343651777285E-003
Determinant=   2.12161563270509       X 10   -2990.00000000000
pos =  -2990.00000000000
neg =   2.12161563270509
zer =  0.000000000000000E+000
Solvemem=   262.000000000000 Beginner
243 Views

Would you please share the equivalent PARDISO version of this code that works? Black Belt
243 Views

The zip file is attached. It contains the input file (from you, but included for consistency), a utility to read the input file and write an unformatted Fortran file in a format that I use when testing problems with MKL/Pardiso, and the program to read that unformated file and solve the problem with Pardiso.

The utility is self-contained, and does not need to be linked with MKL. The solver program should be linked with MKL, and if you have not previously compiled the file mkl_pardiso.f90 that is provided with MKL, you should do so first to produce the module file that is needed by the solver program.

I have tested only on Windows 8.1 X64, but there should be no changes needed for Linux. Please let me know if you run into problems. The last line of the output should be

0.257042531294272       9.162025370381688E-014  4.008343651777284E-003 