Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.

Showing results for

- Intel Community
- Software Development SDKs and Libraries
- Intel® oneAPI Math Kernel Library & Intel® Math Kernel Library
- example code of Direct Sparse Solver (DSS) Interface gives wrong result

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Mute
- Printer Friendly Page

Highlighted

Jiuzhou_T_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-15-2014
07:52 AM

40 Views

** 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

Accepted Solutions

Highlighted

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-16-2014
03:49 AM

40 Views

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.

7 Replies

Highlighted

Jiuzhou_T_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-15-2014
07:57 AM

40 Views

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

Highlighted

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-16-2014
03:49 AM

41 Views

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.

Highlighted

Jiuzhou_T_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-16-2014
04:19 AM

40 Views

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.

Highlighted

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-16-2014
05:13 AM

40 Views

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.

Highlighted

Jiuzhou_T_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-16-2014
07:44 AM

40 Views

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.

Highlighted

mecej4

Black Belt

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-16-2014
08:57 AM

40 Views

The procedure is not complicated.

- 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.
- Then, you supply the nonzero entries in the matrix, and do a numerical factorization.
- 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.
- If you now wish to use the same sparse pattern, but with different nonzero entries, go to Step-2. Otherwise, go to Step-1.

Highlighted
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.

Jiuzhou_T_

Beginner

- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Email to a Friend
- Report Inappropriate Content

12-16-2014
07:28 PM

40 Views

For more complete information about compiler optimizations, see our Optimization Notice.