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

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

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

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.

Link Copied

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

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

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

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.

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

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.

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

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.

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

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.

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

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.

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

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