tests do not check quality of solutions

[dimpase]

Expected behavior

If I do something silly with the test result, e.g.

--- a/EXAMPLES/SYM/ssdrv1.f
+++ b/EXAMPLES/SYM/ssdrv1.f
@@ -296,7 +296,7 @@ c               %---------------------------%
 c
                 call av(nx, v(1,j), ax)
                 call saxpy(n, -d(j,1), v(1,j), 1, ax, 1)
-                d(j,2) = snrm2(n, ax, 1)
+                d(j,2) = snrm2(n, ax, 1)+10000
                 d(j,2) = d(j,2) / abs(d(j,1))
 c
  20          continue

then ssdrv1 "check" should fail, as this is a huge error in the test result.

Actual behavior

It does not.

Where/how to reproduce the problem

modify the code as above and run make check

[fghoussen]

This is likely part of the reference code (caam.rice.edu): you may have stuffs like that all over the place. Impossible to fix them all.
Not sure this is a good move : arpack-ng is meant only to ease/allow compilation/portability on modern archs. Not a dev project.

[dimpase]

On Wed, 13 Oct 2021, 18:10 Franck HOUSSEN, ***@***.***> wrote: This is likely part of the reference code (caam.rice.edu): you may have stuffs like that all over the place. Impossible to fix them all. Not sure this is a good move : arpack-ng is meant only to ease/allow compilation/portability on modern archs. Not a dev project.

ARPACK's dependencies evolve, Fortran compilers evolve, hardware evolves. It's naive to hope that the ability to build and do not crash during the run is a sufficient test, obviously not. With human eyes not checking these outputs much, it's getting even more important to have it fixed. Surely downstream projects do some implicit or explicit testing of ARPACK, but how much is missing? We just discovered a coding error in scipy tests of ARPACK interface, which resulted in the generalised eigenproblems not tested at all, and it was there for like 10 years...

…

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub <#324 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAJXYHCPGJWBB63A5CSWG3LUGW4OTANCNFSM5FYRLXBQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[dimpase]

Impossible to fix them all

One can semi-automatically collect values for bounds in all Fortran programs in EXAMPLES and TESTS that call Xmout (X in ('d','s')) (just modify Xmout) and use them in testing, to parse outputs of tests, for correctness.

[denis-bz]

Fwiw, here's a simple eigcheck.py to check AV - λ MV.
How small do you want it to be, relative / absolute ?
(Does arpack-ng have both ? scipy arpack doesn't.)

#!/usr/bin/env python3
import numpy as np 

#...............................................................................
def eigcheck( A, evals, V, M=None, verbose=1, tag="" ) -> "diffs, maxdiffs":
    """ check eigenvalues and vectors
        -> diffs = A V - M V * evals  n x k,
            maxdiffs = |diffs|_maxnorm
        verbose=1 prints
    """
        # arpack rtol: |AV - λ MV| < tol |λv|, tol and |λv| both tiny ??
        # should atol too, see np.isclose 
    if evals is None or len(evals) == 0:
        print( "eigcheck: no eigvals ?" )  # arpacknoconvergence ?
        return None, None
    MV = V if M is None \
        else M.dot( V )

    diffs = A.dot( V ) - MV * evals  # residuals, n x k
    maxdiffs = np.asarray( np.linalg.norm( diffs, axis=0, ord=np.inf ))  # inf, rms ?
    if verbose:
        jmax = maxdiffs.argmax()
        emax = np.abs( evals[jmax] )
        with np.printoptions( edgeitems=5, threshold=10, formatter=dict(
                float = lambda x: "%.2g" % x )):
            print( "\neigcheck %s: max |AV - VΛ| %.2g at |λ| %.3g " % (
                    tag, maxdiffs[jmax], emax )) 
            print( "  |AV - VΛ|:", maxdiffs )  # may not be sorted

    return diffs, maxdiffs

---

Sample output, for 10k x 10k poisson2 - 4.0001 I, float32, scipy arpack with Lgmres:

eigcheck : max |AV - VΛ| 0.0024 at |λ| 0.00263
  |AV - VΛ|: [0.00018 5.7e-06 4e-05 0.00065 0.00013 4.8e-05 9.8e-05 0.00029 0.0024 7.7e-05]

(evals:  [-9.8e-05 -9.9e-05 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001 -0.0001 0.0026 0.0028]

_{Theory and practice are closer in theory than they are in practice.}

[dimpase]

one way or another, even if arpack-ng is not a "development" project, merely making sure no bitrot is happening requires making tests more robust.

Hoping that downstream does testing for you is wishful thinking -- and pushback I got proposing in scipy/scipy#14846 improving robustness of arpack tests in scipy shows this quite well.

[denis-bz]

Hi @dimpase
yeah, testing ageold spaghetti code is no fun, gets no points. What to do ?
One possibility would be a wiki for arpack etc. with
-- overview
-- basics, for beginners: don't use tol=0, v0=random; what are abstol, reltol for eigenpairs
-- brief descriptions of the labyrinth of solvers for sparse matrices / linops
-- shift-invert: a picture, the posdef trick
-- plots: everybody likes pictures

Do you know of such a wiki in any area of numerical analysis
or is that a chimera ?

Bytheway have you used Pyarpack, could one make it on macos -> openblas in say an hour ?

[dimpase]

I just tried pyarpack on macOS, with Homebrew's openblas and boost-python3, and it builds (with cmake, I didn't try autotools yet)

But I am getting errors in make test

89/108 Test  #89: icb_arpack_c_tst .................   Passed    0.19 sec
        Start  90: icb_arpack_cpp_tst
 90/108 Test  #90: icb_arpack_cpp_tst ...............Subprocess aborted***Exception:   0.14 sec
        Start  91: arpackmm_tst
 91/108 Test  #91: arpackmm_tst .....................   Passed   29.27 sec
        Start  92: pyarpackSparseBiCGDiag_tst
 92/108 Test  #92: pyarpackSparseBiCGDiag_tst .......***Exception: Illegal  1.25 sec
        Start  93: pyarpackSparseBiCGILU_tst
 93/108 Test  #93: pyarpackSparseBiCGILU_tst ........***Exception: Illegal  0.16 sec
        Start  94: pyarpackSparseCGDiag_tst
 94/108 Test  #94: pyarpackSparseCGDiag_tst .........***Exception: Illegal  0.17 sec
        Start  95: pyarpackSparseCGILU_tst
 95/108 Test  #95: pyarpackSparseCGILU_tst ..........***Exception: Illegal  0.22 sec
        Start  96: pyarpackSparseLLT_tst
 96/108 Test  #96: pyarpackSparseLLT_tst ............***Exception: Illegal  0.16 sec
        Start  97: pyarpackSparseLDLT_tst
 97/108 Test  #97: pyarpackSparseLDLT_tst ...........***Exception: Illegal  0.17 sec
        Start  98: pyarpackSparseLU_tst
 98/108 Test  #98: pyarpackSparseLU_tst .............***Exception: Illegal  0.16 sec
        Start  99: pyarpackSparseQR_tst
 99/108 Test  #99: pyarpackSparseQR_tst .............***Exception: Illegal  0.16 sec
        Start 100: pyarpackDenseLLT_tst
100/108 Test #100: pyarpackDenseLLT_tst .............***Exception: Illegal  0.17 sec
        Start 101: pyarpackDenseLDLT_tst
101/108 Test #101: pyarpackDenseLDLT_tst ............***Exception: Illegal  0.16 sec
        Start 102: pyarpackDenseLURR_tst
102/108 Test #102: pyarpackDenseLURR_tst ............***Exception: Illegal  0.16 sec
        Start 103: pyarpackDenseQRRR_tst
103/108 Test #103: pyarpackDenseQRRR_tst ............***Exception: Illegal  0.15 sec
        Start 104: pyarpackDenseLUPP_tst
104/108 Test #104: pyarpackDenseLUPP_tst ............***Exception: Illegal  0.16 sec
        Start 105: pyarpackDenseQRPP_tst
105/108 Test #105: pyarpackDenseQRPP_tst ............***Exception: Illegal  0.16 sec
        Start 106: pyarpackRestart_tst
106/108 Test #106: pyarpackRestart_tst ..............***Exception: Illegal  0.15 sec
        Start 107: icb_parpack_c_tst
107/108 Test #107: icb_parpack_c_tst ................   Passed    0.66 sec
...

[sylvestre]

Tests to verify the quality of the results are welcome.
arpack-ng was mostly started to fix some builds and packaging issues.
We introduced some tests over the years but it is far far from perfect.

There isn't any structure or sponsor behind arpack-ng.
It is mostly @fghoussen now (and I am still reading the various issues and PR) fully volunteering. Therefore, if you care about arpack-ng, please help!

A quick step, we could integrate the coverage in the CI to make sure we don't regress.

$ mkdir build
$ cd build
$ cmake -DCOVERALLS=ON -DCMAKE_BUILD_TYPE=Debug ..
$ make all check test coveralls

[denis-bz]

@sylvestre, folks,
you might like the Poisson test-matrix generator described in this gist
cheers
-- denis