I do not get to solve a system using sparse matrices

[felipin85]

I'm trying to solve a sparse system matrices using mkl_dss.But I have not given them get the right solution, in contrast with getrs getrf and if I succeed.
In the file PROGRAMA07021011PIVOTGAUSS.f90 I'm usinggetrfand getrs.In the file PROGRAMA07022011SPARSE.f90 I'm resolved by sparse matrices.The sparse file, once run
generates five files val.dat, col.dat, rowindex.dat, ld.dat and sp.txt.SP.txt file is the matrix Ax = b, A system.ld.txt is the vector b.While val.dat are non-zero values of A,
rowIndex col.dat your column and row index.I have found these three on hand, and appear well.
The solution is given in the file solutionsparse.dat and are huge numbers.The solution of the file PROGRAMA07021011PIVOTGAUSS.f90 (apparently correct) iscalled solutiongauss.dat.
I tried to solve the system with matlab, inverting the matrix A and multiplying by solutiongauss.dat by match.
Does anyone know why I get these results?

Compile files with this line:

ifort-openmp-O3-o $ MKLINCLUDE PROGRAMA07022011SPARSE-I-L $ MKLPATH PROGRAMA07022011SPARSE.f90-lmkl_lapack95-Wl, - start-group $ MKLPATH / libmkl_intel.a $ MKLPATH / libmkl_intel_thread.a $ MKLPATH / libmkl_core.a-Wl,- end-group-lpthread-liomp5

and with:

ifort-openmp-O3-o $ MKLINCLUDE PROGRAMA07022011PIVOTGAUSS-I-L $ MKLPATH PROGRAMA07022011PIVOTGAUSS.f90-lmkl_lapack95-Wl, - start-group $ MKLPATH / libmkl_intel.a $ MKLPATH / libmkl_intel_thread.a $ MKLPATH / libmkl_core.a-Wl,- end-group-lpthread-liomp5

[Konstantin_A_Intel]

Hi,

Your matrix is almost singular:

Singular values:

max= 52.7334033722656 min= 3.009615796461843E-014

Condition number = 1.752163961734249E+015

Theoretically, solving such matrices in double precision cannot produce stable results.

Results,

Konstantin

[felipin85]

Quoting Konstantin Arturov (Intel)

Hi,

Your matrix is almost singular:

Singular values:

max= 52.7334033722656 min= 3.009615796461843E-014

Condition number = 1.752163961734249E+015

Theoretically, solving such matrices in double precision cannot produce stable results.

Results,

Konstantin

Thanks for your response Konstantin! But, I studied of matrix A (file SP.txt) and I think that this matrix don't singular because its determinat isn't zero, is 4.9330e+03.
Moreover, in the archives of triangulation would have a bug, right?
Thanks for your help

[Konstantin_A_Intel]

Hi,

Of course, term "almost singular" is not exact.. It's better to call this matrix ill-conditioned. So, determinant is not important here.. it doesn't correlate with the condition number that's the main factor of solution stability.

Let's assume:

X_rel -relative error of solution

Cond(A) - matrix condition number

B_rel - relative error of right-hand side.

Following statement is true:

X_rel <= Cond(A)*B_rel, or moving to machine arithmetics:

X_rel <= Cond(A)*eps, where eps is the smallest number representing by double precision type and rougly it's equal to 1e-16 (as mantissa of double precision type has 53 bits and 2^-53 ~= 1e-16). So, having Cond(A)=1e+15 we can write:

X_rel <= 1.0e+15 * 1.0e-16, or

X_rel <= 1.0e-1

So, only 1 significant decimal digit is guaranteed to be correct in the relative error of solution of this matrix.

Regards,

Konstantin

[felipin85]

Ok, Konstantin, now I'm understand for you. Well, I'm changed double precision for simple precision (kind=4), but I have a problem with dss_solve_real, because the subroutine dss_solve_real forces using double precision in input vectors.

PROGRAMA09022011.f90(654): error #6633: The type of the actual argument differs from the type of the dummy argument.
if (dss_solve_real( hDSS, MKL_DSS_FORWARD_SOLVE, P, 1, X) &
-------------------------------------------------------^
PROGRAMA09022011.f90(661): error #6633: The type of the actual argument differs from the type of the dummy argument.

if (dss_solve_real( hDSS, MKL_DSS_BACKWARD_SOLVE, P, 1, X ) &
--------------------------------------------------^
PROGRAMA09022011.f90(661): error #6633: The type of the actual argument differs from the type of the dummy argument.
if (dss_solve_real( hDSS, MKL_DSS_BACKWARD_SOLVE, P, 1, X ) &
--------------------------------------------------------^
PROGRAMA09022011.f90(668): error #6633: The type of the actual argument differs from the type of the dummy argument.

if (dss_solve_real( hDSS, MKL_DSS_DEFAULTS, P, 1, X ) &
--------------------------------------------^
PROGRAMA09022011.f90(668): error #6633: The type of the actual argument differs from the type of the dummy argument.
if (dss_solve_real( hDSS, MKL_DSS_DEFAULTS, P, 1, X ) &
--------------------------------------------------^
compilation aborted for PROGRAMA09022011.f90 (code 1)

Any idea how to calculate my system sparse matrix?
Thanks, for your help.

Regards

[Konstantin_A_Intel]

I'm afraid that with single precision you have much less chances to obtain correct solution as this type has only 23 binary digits for mantissa (comparing to double precision having 53 bits). Do you really have any reason to use single precision computations?

Regards,

Konstantin

[felipin85]

No, I haven't reason to use single precision. I'm using single precision because I I understood that it would improve the system resolution.
The program will have to solve very large sparse systems, from a finite difference scheme. For this case, I can resolution with getrfand getrs, but when I grow up I'll have to calculate the system with sparse.
In theory, Is it normal that the system of good results using getfrf (and getrs) and bad results using sparse?
On Other hand, How matrix ill-conditioned test them? For viewing the finite difference schemes should have a good ill-conditioned?

Sorry for the inconvenience

Thanks and Regards!

[Konstantin_A_Intel]

Hi Felipin85,

I'm sorry for the late response. It seems I have an answer why youdid get completely different solution with PARDISO/DSS and LAPACK. The problem is that you wrote a matrix in the file SP.txt in the following way:

WRITE(52,52)(P(j),j=1,total)

52 FORMAT (10000F10.3)

That means that you lost all significant digits to the left of floating point except 3 ones. That changed the matrix so that its condition number became ~36000 instead of 1e+15 in the initial matrix. That "modified" matrix was passed to dgetf and dgetrs that resulted in good solution (but it was a solution of another matrix). On the other hand, data passed to sparse solver was correct: array Val.dat contains numbers with many significant digits.

So, if you replace line

WRITE(52,52)(P(j),j=1,total)

in PROGRAMA07022011PIVOTGAUSS.f90 with more correct

WRITE(52,*)(P(j),j=1,total)

you will get about the same solution with big numbers as DSS produced.

Regards,

Konstantin

[felipin85]

Thanks Konstantin. I wondered why not give the same solution.But, I was born a doubt, the solutions with WRITE(52,52)(P(j),j=1,total) is well. I would impose three decimal places, for example, 10,430 or something similar. Then understand that are added to the right values, and therefore not "converge", but I know that the results is correct with the values of matrix sp.txt.
How I can get dss values, as I do with the file SP.txt (write (52.52)?

Thanks, and I'm sorry for the inconvenience

Regards

[Konstantin_A_Intel · Accepted solution]

The simplest way is tochange you SPARSE routineas follows (I marked modified/added lines with

!---------------------------):

OPEN(4,FILE='Col.dat')

OPEN(3,FILE='Val.dat')

OPEN(2,FILE='RowIndex.dat')

DO i=1,NVALUES

write(4,*)(COL(i))

write(3,52)(VAL(i)) !---------------------------

END DO

52 FORMAT (10000F10.3) !---------------------------

DO i=1,total

write(*,*)(P(i))

END DO

DO i=1,total+1

write(2,*)(rowindex(i))

END DO

close(4)

close(3)

close(2)

OPEN(3,FILE='Val.dat')!---------------------------

DO i=1,NVALUES !---------------------------

read(3,*)(VAL(i)) !---------------------------

END DO !---------------------------

close(3)

I hope it will help you to resolve your problems!

Best regards,

Konstantin

[felipin85]

As we say in my land: MUCHAS GRACIAS! (very thanks).I will try to solve the case, the only bad thing is having to read the files, because I did that to verify that the entry was good. Should come as vectors in the subroutine.

Anyway, thank you very much Konstantine