Arpack random generator cannot be seeded

[pv]

Expected behavior

It would nice to be able to set the random seed of Arpack's internal random number generator, for cases where more strict fp reproducibility is aimed at (needs of course also compiler and library support to be complete). There could be a new subroutine call in the API that sets it.

Actual behavior

Currently, the rng seeding is done on the first call to each of *getv0.f routines, https://github.com/opencollab/arpack-ng/blob/c43cb868546e8f987cf19b3c2613c56be3865ee3/SRC/dgetv0.f#L202-L208

Specifying an initial vector is not sufficient, since new random vectors are sometimes needed along the run, https://github.com/opencollab/arpack-ng/blob/c43cb868546e8f987cf19b3c2613c56be3865ee3/SRC/dsaitr.f#L409

Where/how to reproduce the problem

• arpack-ng: as of 2019-07

[fghoussen]

You may propose a PR. If so, think about also providing / updating isoc binding API.

[fghoussen]

BTW, not sure to get the pb right. Seeds are fixed: so generated random numbers should be the same if you provide the same initial vector, no ?

[pv]

The seed is initialized only the first time *getv0 is called in a process (save variables), so you'll get different results every time you invoke arpack in the same process. This means e.g. that libraries that use arpack cannot control the randomness, since they don't control the whole-program execution.

[fghoussen]

OK, I can try to PR something about that. So you need something to reset seed at every iteration of the RCI-while loop, right ?

[pv]

I have a branch here: https://github.com/opencollab/arpack-ng/compare/master...pv:rng-seed It should retain the old behavior exactly by default, but I did not yet get around to testing it in real world.

[pv]

By the way, are the iso_c_bindings as they currently are done in arpack-ng correct for the logical type? Adding the interface declaration the compiler complains that logical(kind=c_bool) and logical are not compatible types, and indeed false seems to behave incorrectly if no conversion is done and interface is omitted?

[fghoussen]

I have a branch here: master...pv:rng-seed

Looks a bit fuzzy (many impacts) and incomplete (missing all p*getv0.f). Tried to PR some fix : not tested, @pv could you test if this behaves like you need ? (debug/stat common have already been initialized. I realized that, here, it is about save... Which may not behave exactly like common). If this works, this should cover all cases.

By the way, are the iso_c_bindings as they currently are done in arpack-ng correct for the logical type?

I believe yes. ICB is mainly "typing variables". Not sure how/when mixing ICB and interface is safe. Fortran is an every day pain (personal opinion)... So I avoid to play near border lines (but you may be right, this may be possible to mix them)

[pv]

The current #223 won't do it --- save variables are local to each routine, so there are four different RNG states. Wanting to preserve backward compatibility with that is what necessitates the complexity, so I guess you'd necessarily converge to a similar solution; common blocks could simplify it, though.

PARPACK is more complicated and a similar seeding approach may not be appropriate there -- needs more thinking (and not sure I have the time), so it's not done.

For the other issue, to state it more clearly: the current icb definitions look wrong --- on my platform (gfortran), logical(kind=c_bool) appears to be logical(1) whereas logical becomes logical(4), so the calling convention goes wrong without interface declarations / variable conversions. You can e.g. follow execution of icb_test_c in debugger inside dseupd_c, and notice the value of rvec is not correctly passed. (Or check it adding print statements; most clearly seen changing rvec=false in icb_test_c.c and trying to print rvec at beginning of dseupd.f.)

[fghoussen]

Can't afford more time on this. Hope somebody can.