example code of Direct Sparse Solver (DSS) Interface gives wrong result

[Jiuzhou_T_]

 I am trying to use DSS routine to solver linear equation with a sparse matrix. I found the example code under the intel compiler directory named

dss_sym_f90.f90 and compiled it as ifort dss_sym_f90.f90 -o test-dss -mkl. The code solves a 5*5 linear equation and produces a wrong solution as

"Solution Array:   -326.333   983.000   163.417   398.000    61.500"  

while in fact it should be  (-1.35972222222222,        4.00000000000000       0.250000000000000 ,   6.40000000000000   ,    0.312500000000000).

Source file is pasted as below.  Could someone tell me what might be the reason of the wrong result?  

Another thing, in my application, the matrix of the linear equation is constant while the right side vector changes, is there any way that I can store the LU factorization so I don't have to do it every time ? Thanks a lot!

INCLUDE 'mkl_dss.f90' ! Include the standard DSS "header file."
PROGRAM solver_f90_test
use mkl_dss
IMPLICIT NONE
INTEGER, PARAMETER :: dp = KIND(1.0D0)
INTEGER :: error
INTEGER :: i
INTEGER, PARAMETER :: bufLen = 20
! Define the data arrays and the solution and rhs vectors.
INTEGER, ALLOCATABLE :: columns( : )
INTEGER :: nCols
INTEGER :: nNonZeros
INTEGER :: nRhs
INTEGER :: nRows
REAL(KIND=DP), ALLOCATABLE :: rhs( : )
INTEGER, ALLOCATABLE :: rowIndex( : )
REAL(KIND=DP), ALLOCATABLE :: solution( : )
REAL(KIND=DP), ALLOCATABLE :: values( : )
TYPE(MKL_DSS_HANDLE) :: handle ! Allocate storage for the solver handle.
REAL(KIND=DP),ALLOCATABLE::statOUt( : )
CHARACTER*15 statIn
INTEGER perm(1)
INTEGER buff(bufLen)
! Set the problem to be solved.
nRows = 5
nCols = 5
nNonZeros = 9
nRhs = 1
perm(1) = 0
ALLOCATE( rowIndex( nRows + 1 ) )
rowIndex = (/ 1, 6, 7, 8, 9, 10 /)
ALLOCATE( columns( nNonZeros ) )
columns = (/ 1, 2, 3, 4, 5, 2, 3, 4, 5 /)
ALLOCATE( values( nNonZeros ) )
values = (/ 9.0_DP, 1.5_DP, 6.0_DP, 0.75_DP, 3.0_DP, 0.5_DP, 12.0_DP, &
& 0.625_DP, 16.0_DP /)
ALLOCATE( rhs( nRows ) )
rhs = (/ 1.0_DP, 2.0_DP, 3.0_DP, 4.0_DP, 5.0_DP /)
! 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, rowIndex, nRows, &
& nCols, columns, nNonZeros )
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, values )
IF (error /= MKL_DSS_SUCCESS) GOTO 999
! Allocate the solution vector and solve the problem.
ALLOCATE( solution( nRows ) )
error = DSS_SOLVE_REAL(handle, MKL_DSS_DEFAULTS, rhs, nRhs, solution )
IF (error /= MKL_DSS_SUCCESS) GOTO 999
! Print Out the determinant of the matrix (no statistics for a diagonal matrix)
IF( nRows .LT. nNonZeros ) THEN
ALLOCATE(statOut( 5 ) )
statIn = 'determinant'
call mkl_cvt_to_null_terminated_str(buff,bufLen,statIn)
error = DSS_STATISTICS(handle, MKL_DSS_DEFAULTS, buff, statOut )
IF (error /= MKL_DSS_SUCCESS) GOTO 999
WRITE(*,"('pow of determinant is '(5F10.3))") ( statOut(1) )
WRITE(*,"('base of determinant is '(5F10.3))") ( statOut(2) )
WRITE(*,"('Determinant is '(5F10.3))") ( (10**statOut(1))*statOut(2) )
END IF
! Deallocate solver storage and various local arrays.
error = DSS_DELETE( handle, MKL_DSS_DEFAULTS )
IF (error /= MKL_DSS_SUCCESS ) GOTO 999
IF ( ALLOCATED( rowIndex) ) DEALLOCATE( rowIndex )
IF ( ALLOCATED( columns ) ) DEALLOCATE( columns )
IF ( ALLOCATED( values ) ) DEALLOCATE( values )
IF ( ALLOCATED( rhs ) ) DEALLOCATE( rhs )
IF ( ALLOCATED( statOut ) ) DEALLOCATE( statOut )
! Print the solution vector, deallocate it and exit
WRITE(*,"('Solution Array: '(5F10.3))") ( solution(i), i = 1, nCols )
IF ( ALLOCATED( solution ) ) DEALLOCATE( solution )
GOTO 1000
! Print an error message and exit
999 WRITE(*,*) "Solver returned error code ", error
1000 CONTINUE
END PROGRAM solver_f90_test

[Jiuzhou_T_]

By the way, the compiler I use is of 13.1.1 version, and my system is linux centOS 5.5 .

[mecej4 · Accepted solution]

Apparently, you used only the upper triangle of the matrix with a tool such as Matlab to compute the second "solution" in #1. Fill out the lower triangle of the symmetric matrix before giving it to Matlab, and you will get the same solution as that given by MKL/Lapack.

Regarding your second question: you do not need to store the factor(s), since that is done by DSS internally. Make a single call with DSS_FACTOR_REAL with the symmetric matrix. After that succeeds, make as many calls to DSS_SOLVE_REAL as you wish with different right hand sides from call to call. See the MKL DSS documentation for details.

[Jiuzhou_T_]

  You are right, I was really stupid that I forgot to check it is an sample code for symmetric sparse matrix. 

   In my application, I need to do LU factorization for N different sparse matrixs of the same size ( 81*81) and then solve the equation

one by one. Should I define an array of handle like  

TYPE(MKL_DSS_HANDLE) :: handle(N)   and then solve the different equations within a loop? SInce I don't know how does MKL deal with

the internal work array, do I need to concern the memeory

limit size if N gets very large, say, N>10000?

 And, last queation, for a small size sparse matrix (81*81), is it a good idea to do LU factorization and then

solve it, or maybe I should invert the sparse matrix instead?

Thank you for your reply, I really appreciate it.

[mecej4]

It is never a good idea to form the inverse and use it to solve equations, whether or not the matrix is sparse.

The statements

     "the matrix of the linear equation is constant while the right side vector changes" (in #1)

     "I need to do LU factorization for N different sparse matrixs of the same size" (in #4)

are not consistent. Clarify as to which of these applies to your project.

Since factorization consumes much more time than the subsequent solution phase(s) for a single r.h.s., if the matrix is constant the factorization should be done only once and the factors used as often as needed.

If the sequences of matrices that you mentioned in #4 are structurally the same but numerically different, MKL contains a provision for doing symbolic factorization using Pardiso.

I see no reason to allocate an array of handles. To avoid tying up memory, you should release handles and deallocate arrays that are no longer needed. Likewise, there are provisions to signal to MKL when you are done with the solution and to release memory allocated in previous calls. Please see the MKL manual for details.

[Jiuzhou_T_]

 I am sorry to confuse you. What I need to do is to solve multiple sparse linear equations, for each linear equation, the corresponding matrix is constant during my computation, and it is of size 81*81, so it is not a large matrix.  Besides, I need to solve all these linear equations again and again for many times , thus I can not release the handles that are solved because they are needed for the next time.

About the memory issue, although each matrix is only 81*81, but I got like 64*64*64 numerically different matrixs, so I don't know if I need to take special care of  this problem.

I am sorry I am really poor with English, I hope I made my problem clear this time.

[mecej4]

The procedure is not complicated.

1. Every time that the matrix changes size or the fill pattern changes, you release all old allocations and restart. You can then define the sparse matrix pattern, and call Pardiso to do a "symbolic" factorization, in which only the sparsity pattern is used.
2. Then, you supply the nonzero entries in the matrix, and do a numerical factorization.
3. Using the results of the previous step, you can solve a number of times, each time with a different r.h.s. vector. Note that we repeat Step-3 without going back to Step-2 unless the matrix changes.
4. If you now wish to use the same sparse pattern, but with different nonzero entries, go to Step-2. Otherwise, go to Step-1.

[Jiuzhou_T_]

Yes, all my constant matrixs are of the same sparse pattern, and my pesudo code looks like this do i=1, Many_times do j=1,N_equations call solve_sparse_equations(matrix_j, rhs_j) call change_rhs_vector(rhs_j) enddo enddo The thing is, I have a two lever loop for my problems. Once all the matrixs are factorized, they will be used for the outer loop many times, if I release them within the inner loop while keep their sparsity pattern, then I still have to do a numerical factorization when the outer loop index changes. However, if I do the factorization for all the matrixs once and store all their handles and related work arrays, then I can use them to solve equations without any extra work, while at the cost of some memory. Of course, if doing the numerical factorization doesn't cost much time once I did the "symbolic" factorization already, then I guess I don't have to store all these handles and work arrays at all. So,do you think I need to store the factorization handles and work arrays, or just like you suggested, I can simply do the numerical factorization.