Different RNG seeding strategies in complex vs real versions of parallel getv0()

[szhorvat]

I noticed that in PARPACK, the initial seed is chosen differently in Z/C-prefix complex functions vs D/S-prefix real ones.

In pcgetv0 (and pzgetv0), we have:

https://github.com/opencollab/arpack-ng/blob/master/PARPACK/SRC/MPI/pcgetv0.f#L218

c     %-----------------------------------%
c     | Generate a seed on each processor |
c     | using process id (myid).          |
c     | Note: the seed must be between 1  |
c     | and 4095.  iseed(4) must be odd.  |
c     %-----------------------------------%
c
      call MPI_COMM_RANK(comm, myid, ierr)
      igen = 1000 + 2*myid + 1
      if (igen .gt. 4095) then
         write(0,*) 'Error in p_getv0: seed exceeds 4095!'
      end if
c
      iseed(1) = igen/1000
      igen     = mod(igen,1000)
      iseed(2) = igen/100
      igen     = mod(igen,100)
      iseed(3) = igen/10
      iseed(4) = 7

I.e., each parallel thread gets a unique seed, making sure they don't generate the same sequence of random numbers.

---

But in psgetv0 / pdgetv0, we have the same fixed initial seed as in the serial version:

https://github.com/opencollab/arpack-ng/blob/master/PARPACK/SRC/MPI/pdgetv0.f#L215

c     %-----------------------------------%
c     | Initialize the seed of the LAPACK |
c     | random number generator           |
c     %-----------------------------------%
c
      iseed(1) = 1
      iseed(2) = 3
      iseed(3) = 5
      iseed(4) = 7

This does not look correct to me. The seed should be different for each parallel thread, just like in the Z/C versions.

---

This is not related to the earlier seeding issues I reported. I do not use PARPACK myself, this is just something I noticed while looking at the code.

This goes back to the original PARPACK, and is not something that was changed in ARPACK-NG. To be accurate, the setting of iseed(4) was changed (unnecessarily, as it was already odd), but the general setup wasn't.

Note that the original PARPACK came as two archives: the original code, and a patch. Only the patch has the special seeding based on the process ID. It seems to me that they simply forgot to include the D/S versions in the patch.

[szhorvat]

Also, this does not look right:

      if (igen .gt. 4095) then
         write(0,*) 'Error in p_getv0: seed exceeds 4095!'
      end if

The requirement on the elements of iseed are:

1. Each element must be in the integer range [0, 4096), representing a 12-bit value
2. The last element, iseed(4), must be odd

To achieve this, it is not necessary for igen to be < 4096. iseed(2), iseed(3) and iseed(4) are clearly single-digit numbers. iseed(1) will be less than 4096 if igen < 4096000.

---

I'll note that the original seeding PARPACK looked like this:

         iseed(1) = igen/1000
         igen     = mod(igen,1000)
         iseed(2) = igen/100
         igen     = mod(igen,100)
         iseed(3) = igen/10
         iseed(4) = mod(igen,10)

In APRACK-NG this is changed to iseed(4) = 7, presumably to ensure that it's odd.

This is unnecessary since igen is set up to be odd (1000 + 2*myid + 1), therefore it will stay odd modulo powers of 10 as well.

[fghoussen]

original PARPACK came as two archives: the original code, and a patch. Only the patch has the special seeding based on the process ID. It seems to me that they simply forgot to include the D/S versions in the patch.

You're likely right. Feel free to propose a PR to keep things symmetric.

[fghoussen]

@szhorvat: do you plan to PR a patch before the next release?

[szhorvat]

I won't have a computer for the next 10 days

[szhorvat]

In APRACK-NG this is changed to iseed(4) = 7, presumably to ensure that it's odd.

This is unnecessary since igen is set up to be odd (1000 + 2*myid + 1), therefore it will stay odd modulo powers of 10 as well.

I believe this is actually harmful as it fixes the least significant digits, causing the seed to be potentially same on multiple parallel threads. I will fix this as well in the upcoming patch.