Wrong results with mode=LR

[pv]

The following in Octave produces wrong results (arpack-ng 3.4.0):

octave:31> load Q.h5
octave:32> opts = struct(); opts.p=30; eigs(Q, 10, 'lr', opts)
ans =

   7.6562 - 0.1236i
   7.6562 + 0.1236i
   7.0731 - 0.9404i
   7.0731 + 0.9404i
   6.7639 + 0.3453i
   6.7639 - 0.3453i
   6.7362 + 1.3872i
   6.7362 - 1.3872i
   6.3870 - 2.8622i
   6.3870 + 2.8622i

octave:33> max(abs(eig(Q)))
ans =  0.99615

Behavior in Scipy is similar, so I assume this is an ARPACK issue. For small p/ncv it sometimes converges to the right result, but this appears to depend on the platform. Problematic matrix is attached.

Q.zip

[caliarim]

I confirm a similar bad behavior in plain fortran77 arpack 3.3.0 (see attachment). QLR.zip

[jfrench7]

I can confirm a similar problem with the double complex fortran77 interface (non-parallel). Appears limited to the generalized version (bmat="G"). Fails for both inverse and shift-invert approaches. The calculated eigenvalues are all very small (near zero) if no shift. If shift is provided, the calculated eigenvalues are all only a short distance from the shift (and incorrect). Also GENERALIZED takes more iterations to arrive at the incorrect eigenvalues.

Good solution (REGULAR): 1 ( -1.8011626537770042E-010, -3.1991767215096324E-007) 2 ( 1.8011641411523629E-010, 3.1991767239723558E-007) 3 ( -1.3337945980083915E-004, -1.0000991523726590 ) 4 ( -1.1135290553387030E-004, 1.0000354627812349 ) 5 ( 3.3712691258944215E-004, -1.9995440269971894 ) 6 ( -4.2092179853785746E-004, 2.0003063448149070 ) 7 ( 6.7003182745910305E-004, -2.9994437086254644 ) 8 ( -1.0980787053557880E-003, 2.9993211868679150 ) 9 ( 3.1727979066692451E-003, -3.8237599342140709 ) 10 ( 1.3019595693017493E-003, -3.9588662326796387 )

Bad solution (INVERSE GENERALIZED - that is, have to invert M in A x = lambda M x) (Hm, not converging...)

Bad solution (SHIFT-INVERT GENERALIZED) 1 ( -1.6013468986431114E-011, 3.4504365595530612E-012) 2 ( -1.5927276776254480E-011, 3.9370852045317459E-012) 3 ( 1.4794648492903501E-011, -5.2626154957930441E-012) 4 ( 1.6360771042661942E-011, -4.9785373216198306E-012) 5 ( 1.4122457823096880E-011, -4.0415333358261451E-012) 6 ( 1.3408523769846035E-011, -1.2747952241778692E-012) 7 ( 1.6282531649160774E-011, -3.4373813537907268E-012) 8 ( 1.4109761755846917E-011, -3.0399063083829589E-012) 9 ( 1.3875791374653936E-011, -1.8096441627036827E-012) 10 ( 1.4161930552677673E-011, -2.2744426194783752E-012)

Bad solution (SHIFT-INVERT GENERALIZED), shift = 0.5 i : 1 ( -1.0927525347437343E-003, 0.49980664795990370 ) 2 ( -1.1026173426873463E-003, 0.50012759287811959 ) 3 ( -1.1448648646392055E-003, 0.49973252766587539 ) 4 ( -1.1592188681961411E-003, 0.50019836647383842 ) 5 ( -1.1166709986802983E-003, 0.49958716227140459 ) 6 ( -1.1567457955799208E-003, 0.49967133644217882 ) 7 ( -1.0697545001081275E-003, 0.49947718302957839 ) 8 ( -1.1393703666596109E-003, 0.50034531331623433 ) 9 ( -1.1743616882286976E-003, 0.50025867171943461 ) 10 ( -1.0992109445576406E-003, 0.50045820341652714 )

[sylvestre]

Could you please attach the test case? Thanks

Le 1 février 2017 19:32:29 GMT+01:00, jfrench7 [email protected] a écrit :

I can confirm a similar problem with the double complex fortran77 interface (non-parallel). Appears limited to the generalized version (bmat="G"). Fails for both inverse and shift-invert approaches. The calculated eigenvalues are all very small (near zero) if no shift. If shift is provided, the calculated eigenvalues are all only a short distance from the shift (and incorrect). Also GENERALIZED takes more iterations to arrive at the incorrect eigenvalues.

Good solution (REGULAR): 1 ( -1.8011626537770042E-010, -3.1991767215096324E-007) 2 ( 1.8011641411523629E-010, 3.1991767239723558E-007) 3 ( -1.3337945980083915E-004, -1.0000991523726590 ) 4 ( -1.1135290553387030E-004, 1.0000354627812349 ) 5 ( 3.3712691258944215E-004, -1.9995440269971894 ) 6 ( -4.2092179853785746E-004, 2.0003063448149070 ) 7 ( 6.7003182745910305E-004, -2.9994437086254644 ) 8 ( -1.0980787053557880E-003, 2.9993211868679150 ) 9 ( 3.1727979066692451E-003, -3.8237599342140709 ) 10 ( 1.3019595693017493E-003, -3.9588662326796387 )

Bad solution (INVERSE GENERALIZED - that is, have to invert M in A x = lambda M x) (Hm, not converging...)

Bad solution (SHIFT-INVERT GENERALIZED) 1 ( -1.6013468986431114E-011, 3.4504365595530612E-012) 2 ( -1.5927276776254480E-011, 3.9370852045317459E-012) 3 ( 1.4794648492903501E-011, -5.2626154957930441E-012) 4 ( 1.6360771042661942E-011, -4.9785373216198306E-012) 5 ( 1.4122457823096880E-011, -4.0415333358261451E-012) 6 ( 1.3408523769846035E-011, -1.2747952241778692E-012) 7 ( 1.6282531649160774E-011, -3.4373813537907268E-012) 8 ( 1.4109761755846917E-011, -3.0399063083829589E-012) 9 ( 1.3875791374653936E-011, -1.8096441627036827E-012) 10 ( 1.4161930552677673E-011, -2.2744426194783752E-012)

Bad solution (SHIFT-INVERT GENERALIZED), shift = 0.5 i : 1 ( -1.0927525347437343E-003, 0.49980664795990370 ) 2 ( -1.1026173426873463E-003, 0.50012759287811959 ) 3 ( -1.1448648646392055E-003, 0.49973252766587539 ) 4 ( -1.1592188681961411E-003, 0.50019836647383842 ) 5 ( -1.1166709986802983E-003, 0.49958716227140459 ) 6 ( -1.1567457955799208E-003, 0.49967133644217882 ) 7 ( -1.0697545001081275E-003, 0.49947718302957839 ) 8 ( -1.1393703666596109E-003, 0.50034531331623433 ) 9 ( -1.1743616882286976E-003, 0.50025867171943461 ) 10 ( -1.0992109445576406E-003, 0.50045820341652714 )

-- You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub: https://github.com/opencollab/arpack-ng/issues/51#issuecomment-276740304

-- Envoyé de mon appareil Android avec K-9 Mail. Veuillez excuser ma brièveté.

[jfrench7]

If I print the A and M matrices in ASCII, will that suffice?

• Jonathan

On 2/1/2017 12:34 PM, Sylvestre Ledru wrote:

Could you please attach the test case? Thanks

Le 1 février 2017 19:32:29 GMT+01:00, jfrench7 [email protected] a écrit :

I can confirm a similar problem with the double complex fortran77 interface (non-parallel). Appears limited to the generalized version (bmat="G"). Fails for both inverse and shift-invert approaches. The calculated eigenvalues are all very small (near zero) if no shift. If shift is provided, the calculated eigenvalues are all only a short distance from the shift (and incorrect). Also GENERALIZED takes more iterations to arrive at the incorrect eigenvalues.

Good solution (REGULAR): 1 ( -1.8011626537770042E-010, -3.1991767215096324E-007) 2 ( 1.8011641411523629E-010, 3.1991767239723558E-007) 3 ( -1.3337945980083915E-004, -1.0000991523726590 ) 4 ( -1.1135290553387030E-004, 1.0000354627812349 ) 5 ( 3.3712691258944215E-004, -1.9995440269971894 ) 6 ( -4.2092179853785746E-004, 2.0003063448149070 ) 7 ( 6.7003182745910305E-004, -2.9994437086254644 ) 8 ( -1.0980787053557880E-003, 2.9993211868679150 ) 9 ( 3.1727979066692451E-003, -3.8237599342140709 ) 10 ( 1.3019595693017493E-003, -3.9588662326796387 )

Bad solution (INVERSE GENERALIZED - that is, have to invert M in A x = lambda M x) (Hm, not converging...)

Bad solution (SHIFT-INVERT GENERALIZED) 1 ( -1.6013468986431114E-011, 3.4504365595530612E-012) 2 ( -1.5927276776254480E-011, 3.9370852045317459E-012) 3 ( 1.4794648492903501E-011, -5.2626154957930441E-012) 4 ( 1.6360771042661942E-011, -4.9785373216198306E-012) 5 ( 1.4122457823096880E-011, -4.0415333358261451E-012) 6 ( 1.3408523769846035E-011, -1.2747952241778692E-012) 7 ( 1.6282531649160774E-011, -3.4373813537907268E-012) 8 ( 1.4109761755846917E-011, -3.0399063083829589E-012) 9 ( 1.3875791374653936E-011, -1.8096441627036827E-012) 10 ( 1.4161930552677673E-011, -2.2744426194783752E-012)

Bad solution (SHIFT-INVERT GENERALIZED), shift = 0.5 i : 1 ( -1.0927525347437343E-003, 0.49980664795990370 ) 2 ( -1.1026173426873463E-003, 0.50012759287811959 ) 3 ( -1.1448648646392055E-003, 0.49973252766587539 ) 4 ( -1.1592188681961411E-003, 0.50019836647383842 ) 5 ( -1.1166709986802983E-003, 0.49958716227140459 ) 6 ( -1.1567457955799208E-003, 0.49967133644217882 ) 7 ( -1.0697545001081275E-003, 0.49947718302957839 ) 8 ( -1.1393703666596109E-003, 0.50034531331623433 ) 9 ( -1.1743616882286976E-003, 0.50025867171943461 ) 10 ( -1.0992109445576406E-003, 0.50045820341652714 )

-- You are receiving this because you are subscribed to this thread. Reply to this email directly or view it on GitHub: https://github.com/opencollab/arpack-ng/issues/51#issuecomment-276740304

-- Envoyé de mon appareil Android avec K-9 Mail. Veuillez excuser ma brièveté.

— You are receiving this because you commented. Reply to this email directly, view it on GitHub https://github.com/opencollab/arpack-ng/issues/51#issuecomment-276740730, or mute the thread https://github.com/notifications/unsubscribe-auth/AYSStUrgcKOm_CjqP38rqBtEeaP0qEeTks5rYNAagaJpZM4LDWjX.

[pv]

@jfrench7: maybe open a separate ticket for that --- this issue was only for the regular case without shift-invert, so it's probably not the same as yours. Also, there's a fortran77 interface test case for this one above.

[jfrench7]

Ok - my point was that it just doesn't work for either the generalized
inverse or the generalized shift-invert modes (vs standard regular or standard shift-invert modes).

On 2/1/2017 1:39 PM, Pauli Virtanen wrote:

@jfrench7 https://github.com/jfrench7: maybe open a separate ticket for that --- this issue was only for the regular case without shift-invert.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/opencollab/arpack-ng/issues/51#issuecomment-276759074, or mute the thread https://github.com/notifications/unsubscribe-auth/AYSStdHbSxcjvCVWzwo7-jdowl4h2I17ks5rYN94gaJpZM4LDWjX.

[caliarim]

Hi. If you check the orthogonality of the vectors v after the second call of dnaitr in dnaup2, you will see that orthogonality is lost (some vectors have scalar product of size 1e-5). Although dnaitr features double orthogonalization, it seems that classical Gram Schmidt is not appropriate. So, I implemented the modified version, and the original problem has disappeared. I attach the new dnaitr.f file. It is an important modification, so it should be extensively tested (make check passes). In particular, generalized problems should be tested. Then, the modification should be propagated to Xnaitr.f. If I have some feedback, I can try to do a pull request. dnaitr.f.zip

[caliarim]

My previous attempt does not work in the generalized case. Still, the problem is the orthogonalization. There are some possibile solutions:

1. a second orthogonalization is made in some cases, depending on a test. If you change 0.717 with 0.817 in the test in dnaitr (two-line change), the bug here described is solved.

2. it is possibile to switch from classical Gram Schmidt to modified Gram Schmidt in the standard case. It could be slower, since it requires some level-1 operations instead of 2 level-2 operations. But it is more stable, and the second orthogonalization should be less necessary on the average. It is not possible to use modified Gram Schmidt in the generalized case without introducing more storage and computational cost. In this case, you have to stay with classical Gram Schmidt and you could use 0.817 instead of 0.717.

3. always perform a second orthogonalization with classical Gram Schmidt (one-line change). It is the safer and the more expensive solution. If you always force a second orthogonalization there is no need of modified Gram Schmidt.

I would be slightly in favor of 2).

[caliarim]

any comment on this?

[svrnwnsch]

I have no clue whats best to implement (considering speed) but I think its more important to have correct results then having a fast (but wrong) calculation.

[fghoussen]

You may play/test arpack changing options on this specific problem (arpackmm #157) => it may shed some light and help understand if the problem comes from arpack, octave or scypi...

#157 is not merged yet, but you could clone at my side and test with that clone:

~> git clone https://github.com/fghoussen/arpack-ng
~> git checkout arpackmm
~> sudo apt-get install libeigen3-dev
~> ./bootstrap; ./configure --enable-icb-exmm; make all check