Dpocon and Dgecon mismatch

[btracey]

Our tests at Gonum for Dpocon compare to the result for Dgecon. For most matrices (randomly generated), the two outputs match to machine precision.

The matrix
[2.9995576045549965, -2.0898894566158663, 3.965560740124006,
-2.0898894566158663, 1.9634729526261008, -2.8681002706874104,
3.965560740124006, -2.8681002706874104, 5.502416670471008]

Gives a significantly different answer for the two routines. Dgecon gives 0.020928829230190116 while Dpocon gives 0.024054837369015203. The output is the same for both routines in OpenBLAS and gonum/lapack, which suggests a problem with the algorithm rather than the implementation.

[xianyi]

Thank you for the feedback. Yes, it looks like an algorithm issue.

[btracey]

Should a bug with netlib be filed?

[xianyi]

@btracey , I think so.

[btracey]

Filed http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4817

[martin-frbg]

As far as I can tell there has been no activity on the netlib bug since it was filed over two years ago. From the wording of the documentation of both functions "an estimate is obtained for norm(inv(A))..." it is not clear to me if there is any implied guarantee that both should always return the same result, even if it might be convenient for cross-checking between implementations. What is the current stance of Gonum on this issue - worthwile to pursue, or an oddity that was sufficiently worked around
by modifying the test ?

[vladimir-ch]

Finally I've had time to look into this more closely, motivated also by seeing similar differences recently in the condition number estimation for symmetric positive definite band matrices. The main computational part in the estimation is the DLACN2 routine which is an iterative algorithm to estimate the 1-norm of a matrix. The Fortran implementation has a reference to a paper where the algorithm, some theory and numerical experiments are presented. The first point in the summary of their experiments says:

Approximately three-fifths of the estimates were exact (i.e., the computed relative error was of the order of the unit roundoff). The proportion of exact estimates varies greatly with the type of matrix.

They also report that in all their tests the estimate was correct to within a factor of 10 which is what I use in Gonum in testlapack.DpbconTest and the test passes. Tests in the reference allow a more generous difference of factor 30.

In other words, a not-too-large difference is normal and I think you can close this issue.

[martin-frbg]

Great, thanks for researching this.

[vladimir-ch]

One more remark: my comment above of course does not answer why there can be a difference between DGECON and DPOCON, only why there can be a difference between the estimated and exact condition number. The former comes most likely from DGETRF using pivoting which leads to DLACN2 producing different sequences of vectors that may converge to different results. If the above matrix is permuted so that DGETRF doesn't do it, both DPOCON and DGECON give the same (within eps) estimate for the condition number.