Re: eigs and ARPACK

[Thomas Treichl]

David Bateman schrieb:
> John W. Eaton wrote:
> > On  3-Dec-2008, Søren Hauberg wrote:
> >
> > |   It was just pointed out on the Octave-Forge list that ARPACK seems to
> > | have been released under a standard BSD license [1]. Doesn't this mean
> > | that we can move the code from the 'arpack' package into core Octave,
> > | which would give us the 'eigs' and 'svds' functions?
> > | | Søren
> > | | [1] http://www.caam.rice.edu/software/ARPACK/RiceBSD.txt
> >
> > Yes, this is great news.  Would someone like to prepare a patch?
> >
> > jwe
>
> Well attached is an untested port of the octave-forge code. I don't know 
> if I'll have the time to test this code before a few days and so I'm 
> sending it even if it is untested. It also doesn't support single 
> precision matrices and probably should..
>
> Regards
> David

Hi David,

I've had a first look at your patch especially because I wanted to prepare the 
build scripts for Octave.app on MacOSX if ARPACK becomes another dependency of 
Octave. Please understand that I just want to help by reporting my experiences 
and not to criticize your great work ;)

The first experience is that the sources can be compiled even if the 
'./configure' check for 'libarpack' fails. I thought this should be mentioned, too.

The '-larpack' check in './configure' failed for me because the linker needs the 
values of my $FLIBS (which currently are one -L... option and the option -lf2c). 
I don't know if this is the best solution but I changed your line of 
'configure.in' from

  AC_CHECK_LIB(arpack, F77_FUNC(dseupd,DSEUPD), [ARPACK_LIBS=...

into

  save_LIBS="$LIBS"; LIBS="$LIBS $FLIBS"
  AC_CHECK_LIB(arpack, F77_FUNC(dseupd,DSEUPD), [ARPACK_LIBS="...
  LIBS="$save_LIBS"

then the test worked and '-larpack' was found and used correctly.

Linking of liboctave failed for me. I have seen that you added 'eigs-base.cc' to 
'TEMPLATE_SRC' in 'liboctave/Makefile.in'. I also have seen that some F77_FUNCs 
are used inside of 'eigs-base.cc'. I added '$(ARPACK_LIBS)' to 'LINK_DEPS' in 
file 'Makefile.in' then linking was successful.

Compilation then finished but some tests failed and some crashed Octave. I also 
should notice that *all* tests of line 723ff '%% Real positive definite tests, n 
must be even' successfully pass.

I have not had a closer look right now what is failing but I can send some of 
the output that I see. At least these tests crash Octave:

  line 983:

    %!test
    %! d1 = eigs (A, k, 'sm');
    %! assert (abs(d1), abs(d0(1:k)), 1e-12);

  line 1009:

    %!test
    %! d1 = eigs(A, speye(n), k, 'lm');
    %! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
    %!test
    %! opts.cholB=true;
    %! d1 = eigs(A, speye(n), k, 'lm', opts);
    %! assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);

  line 1023:

    %!test
    %! opts.cholB=true;
    %! d1 = eigs(A, speye(n), k, 4.1, opts);
    %! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
    %! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);
    %!test
    %! opts.cholB=true;
    %! q = [2:n,1];
    %! opts.permB=q;
    %! d1 = eigs(A, speye(n)(q,q), k, 4.1, opts);
    %! assert (abs(abs(d1)), abs(eigs(A,k,4.1)), 1e-12);
    %! assert (sort(imag(abs(d1))), sort(imag(eigs(A,k,4.1))), 1e-12);

the message is always the same

  octave:1> test eigs
  complex division by zero
  Abort trap
  Me:~/

The output of some of the tests that fail for me (but do not crash Octave) are:

  ***** test
 d1 = eigs (A, k);
 assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
!!!!! test failed
eigs: error -8 in dnaupdshared variables {
  n =  20
  k =  4
  A =

Compressed Column Sparse (rows = 20, cols = 20, nnz = 56 [14%])
)

  (1, 1) ->  1
  (3, 1) ->  1
  (2, 2) ->  2
  (4, 2) ->  1
  (1, 3) -> -1
  (3, 3) ->  3
  (5, 3) ->  1
  (2, 4) -> -1
  (4, 4) ->  4
  (6, 4) ->  1
  (3, 5) -> -1
  (5, 5) ->  5
  (7, 5) ->  1
  (4, 6) -> -1
  (6, 6) ->  6
  (8, 6) ->  1
  (5, 7) -> -1
  (7, 7) ->  7
  (9, 7) ->  1
  (6, 8) -> -1
  (8, 8) ->  8
  (10, 8) ->  1
  (7, 9) -> -1
  (9, 9) ->  9
  (11, 9) ->  1
  (8, 10) -> -1
  (10, 10) ->  10
  (12, 10) ->  1
  (9, 11) -> -1
  (11, 11) ->  11
  (13, 11) ->  1
  (10, 12) -> -1
  (12, 12) ->  12
  (14, 12) ->  1
  (11, 13) -> -1
  (13, 13) ->  13
  (15, 13) ->  1
  (12, 14) -> -1
  (14, 14) ->  14
  (16, 14) ->  1
  (13, 15) -> -1
  (15, 15) ->  15
  (17, 15) ->  1
  (14, 16) -> -1
  (16, 16) ->  16
  (18, 16) ->  1
  (15, 17) -> -1
  (17, 17) ->  17
  (19, 17) ->  1
  (16, 18) -> -1
  (18, 18) ->  18
  (20, 18) ->  1
  (17, 19) -> -1
  (19, 19) ->  19
  (18, 20) -> -1
  (20, 20) ->  20

  d0 =

      1.5933
      2.5933
      2.9036
      3.9036
      5.0031
      6.0031
      6.9999
      7.9999
      9.0000
     10.0000
     11.0000
     12.0000
     13.0001
     14.0001
     14.9969
     15.9969
     17.0964
     18.0964
     18.4067
     19.4067

}
  ***** test
 d1 = eigs (A,k+1);
 assert (abs(d1), abs(d0(end:-1:(end-k))),1e-12);
!!!!! test failed
eigs: error -8 in dnaupdshared variables {
  n =  20
  k =  4
  A =

Compressed Column Sparse (rows = 20, cols = 20, nnz = 56 [14%])
)

  (1, 1) ->  1
  (3, 1) ->  1
  (2, 2) ->  2
  (4, 2) ->  1
  (1, 3) -> -1
  (3, 3) ->  3
  (5, 3) ->  1
  (2, 4) -> -1
  (4, 4) ->  4
  (6, 4) ->  1
  (3, 5) -> -1
  (5, 5) ->  5
  (7, 5) ->  1
  (4, 6) -> -1
  (6, 6) ->  6
  (8, 6) ->  1
  (5, 7) -> -1
  (7, 7) ->  7
  (9, 7) ->  1
  (6, 8) -> -1
  (8, 8) ->  8
  (10, 8) ->  1
  (7, 9) -> -1
  (9, 9) ->  9
  (11, 9) ->  1
  (8, 10) -> -1
  (10, 10) ->  10
  (12, 10) ->  1
  (9, 11) -> -1
  (11, 11) ->  11
  (13, 11) ->  1
  (10, 12) -> -1
  (12, 12) ->  12
  (14, 12) ->  1
  (11, 13) -> -1
  (13, 13) ->  13
  (15, 13) ->  1
  (12, 14) -> -1
  (14, 14) ->  14
  (16, 14) ->  1
  (13, 15) -> -1
  (15, 15) ->  15
  (17, 15) ->  1
  (14, 16) -> -1
  (16, 16) ->  16
  (18, 16) ->  1
  (15, 17) -> -1
  (17, 17) ->  17
  (19, 17) ->  1
  (16, 18) -> -1
  (18, 18) ->  18
  (20, 18) ->  1
  (17, 19) -> -1
  (19, 19) ->  19
  (18, 20) -> -1
  (20, 20) ->  20

  d0 =

      1.5933
      2.5933
      2.9036
      3.9036
      5.0031
      6.0031
      6.9999
      7.9999
      9.0000
     10.0000
     11.0000
     12.0000
     13.0001
     14.0001
     14.9969
     15.9969
     17.0964
     18.0964
     18.4067
     19.4067

}
  ***** test
 d1 = eigs (A, k, 'lm');
 assert (abs(d1), abs(d0(end:-1:(end-k+1))), 1e-12);
!!!!! test failed
eigs: error -8 in dnaupdshared variables {
  n =  20
  k =  4
  A =

Compressed Column Sparse (rows = 20, cols = 20, nnz = 56 [14%])
)

  (1, 1) ->  1
  (3, 1) ->  1
  (2, 2) ->  2
  (4, 2) ->  1
  (1, 3) -> -1
  (3, 3) ->  3
  (5, 3) ->  1
  (2, 4) -> -1
  (4, 4) ->  4
  (6, 4) ->  1
  (3, 5) -> -1
  (5, 5) ->  5
  (7, 5) ->  1
  (4, 6) -> -1
  (6, 6) ->  6
  (8, 6) ->  1
  (5, 7) -> -1
  (7, 7) ->  7
  (9, 7) ->  1
  (6, 8) -> -1
  (8, 8) ->  8
  (10, 8) ->  1
  (7, 9) -> -1
  (9, 9) ->  9
  (11, 9) ->  1
  (8, 10) -> -1
  (10, 10) ->  10
  (12, 10) ->  1
  (9, 11) -> -1
  (11, 11) ->  11
  (13, 11) ->  1
  (10, 12) -> -1
  (12, 12) ->  12
  (14, 12) ->  1
  (11, 13) -> -1
  (13, 13) ->  13
  (15, 13) ->  1
  (12, 14) -> -1
  (14, 14) ->  14
  (16, 14) ->  1
  (13, 15) -> -1
  (15, 15) ->  15
  (17, 15) ->  1
  (14, 16) -> -1
  (16, 16) ->  16
  (18, 16) ->  1
  (15, 17) -> -1
  (17, 17) ->  17
  (19, 17) ->  1
  (16, 18) -> -1
  (18, 18) ->  18
  (20, 18) ->  1
  (17, 19) -> -1
  (19, 19) ->  19
  (18, 20) -> -1
  (20, 20) ->  20

  d0 =

      1.5933
      2.5933
      2.9036
      3.9036
      5.0031
      6.0031
      6.9999
      7.9999
      9.0000
     10.0000
     11.0000
     12.0000
     13.0001
     14.0001
     14.9969
     15.9969
     17.0964
     18.0964
     18.4067
     19.4067

}

That's it right now, I hope I could give some feedback...

Best regards,

  Thomas