porting arapck code to octave

classic Classic list List threaded Threaded
14 messages Options
Reply | Threaded
Open this post in threaded view
|

porting arapck code to octave

David Bateman-2
I've been looking at porting the eigs function to octave and basically
have it finished.. However I have a bizarre issue in that

eigs (diag(1:100),4,'LM')

doesn't work due to an internal error in an ARPACK function.. I'd
remember something like this and that I'd patch arpack to avoid it on my
machine in the past. However, I thought I'd investigate the issue
further and modified an ARPACK example function "dssimp" to treat the
exact same case as above and no bug manifested..

I then tried reimplementing dssimp  as an oct-file, trying to keep it
close to the same functionality as the fortran code (ie arrays
initialized to zero, etc), and the result was I saw the same bug in the
DSAUPD function as previously. I've been staring at this code for hours
and don't see what the issue is. I attach my modified dssimp.f and the
code that duplicates it in octave with "dssimp(diag(1:100))"... Anyone
want to confirm that

gfortran -o dssimp dssimp.f -larpack -lblas -llapack
./dssimp

works and gives the 4 largest eigenvalues [100, 99, 98, 97] as expected.
whereas

mkoctfile dssimp.cc -larpack
octave --eval "dssimp(diag(1:100))"

doesn't... Anyone want to look at the code and see if they can see an
issue that I missed?

Regards
David




--
David Bateman                                [hidden email]
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)


#include <octave/oct.h>
#include "f77-fcn.h"

extern "C"
{

  F77_RET_T
  F77_FUNC (dsaupd, DSAUPD) (octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
                             const octave_idx_type&, F77_CONST_CHAR_ARG_DECL,
                             const octave_idx_type&, const double&,
                             double*, const octave_idx_type&, double*,
                             const octave_idx_type&, octave_idx_type*,
                             octave_idx_type*, double*, double*,
                             const octave_idx_type&, octave_idx_type&
                             F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (dseupd, DSEUPD) (const int&, F77_CONST_CHAR_ARG_DECL,
                             int*, double*, double*,
                             const octave_idx_type&, const double&,
                             F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
                             F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
                             const double&, double*, const octave_idx_type&,
                             double*, const octave_idx_type&, octave_idx_type*,
                             octave_idx_type*, double*, double*,
                             const octave_idx_type&, octave_idx_type&
                             F77_CHAR_ARG_LEN_DECL F77_CHAR_ARG_LEN_DECL
                             F77_CHAR_ARG_LEN_DECL);

  F77_RET_T
  F77_FUNC (dgemv, DGEMV) (F77_CONST_CHAR_ARG_DECL,
                           const octave_idx_type&, const octave_idx_type&, const double&,
                           const double*, const octave_idx_type&, const double*,
                           const octave_idx_type&, const double&, double*,
                           const octave_idx_type&
                           F77_CHAR_ARG_LEN_DECL);

}

#define maxn 256
#define maxnev 10
#define maxncv 25
#define ldv maxn

static bool
vector_product (const Matrix& m, const double *x, double *y)
{
  octave_idx_type nr = m.rows ();
  octave_idx_type nc = m.cols ();

  F77_XFCN (dgemv, DGEMV, (F77_CONST_CHAR_ARG2 ("N", 1),
                           nr, nc, 1.0,  m.data (), nr,
                           x, 1, 0.0, y, 1
                           F77_CHAR_ARG_LEN (1)));

  if (f77_exception_encountered)
    {
      (*current_liboctave_error_handler)
        ("dssimp: unrecoverable error in dgemv");
      return false;
    }
  else
    return true;
}

DEFUN_DLD (dssimp, args, nargout, "no help")
{
  int nargin = args.length ();
  octave_value_list retval;

  if (nargin != 1)
    print_usage ();
  else
    {
      Matrix m = args(0).matrix_value ();

      if (!error_state)
        {

#define MY_LOCAL_BUFFER(T, V, N) \
          OCTAVE_LOCAL_BUFFER (T, V, N); \
          for (octave_idx_type iii = 0; iii < N; iii++) \
            V[iii] = 0.;

          MY_LOCAL_BUFFER (double, v, ldv*maxncv);  

          MY_LOCAL_BUFFER (double, workl, maxncv * (maxncv+8));
          MY_LOCAL_BUFFER (double, workd, 3 * maxn);
          MY_LOCAL_BUFFER (double, d, maxncv*2);
          MY_LOCAL_BUFFER (double, resid, maxn);
          MY_LOCAL_BUFFER (double, ax, maxn);

          MY_LOCAL_BUFFER (int, select, maxncv);

          MY_LOCAL_BUFFER (octave_idx_type, iparam, 11);
          MY_LOCAL_BUFFER (octave_idx_type, ipntr, 11);

          for (int i = 0; i < 11; i++)
            {
              iparam[i] = 0;
              ipntr[i] = 0;
            }
          iparam[0] = 1;
          iparam[2] = 300;
          iparam[6] = 1;

          octave_idx_type nev = 4;
          octave_idx_type ncv = 20;
          octave_idx_type ido = 0;
          octave_idx_type lworkl = ncv * ( ncv + 8);
          char bmat = 'I';
          std::string which = "LM";
          double tol = 0.;
          octave_idx_type info = 0;
         
          octave_idx_type n = m.rows();
          octave_idx_type lldv = ldv;

          while (true)
            {
              F77_FUNC (dsaupd, DSAUPD)
                (ido, F77_CONST_CHAR_ARG2 (&bmat, 1), n,
                 F77_CONST_CHAR_ARG2 ((which.c_str ()), 2),
                 nev, tol, resid, ncv, v, lldv, iparam,
                 ipntr, workd, workl, lworkl, info
                 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));

              if (f77_exception_encountered)
                {
                  (*current_liboctave_error_handler)
                    ("dssimp: unrecoverable exception encountered in dsaupd");
                  return retval;
                }

              if (ido == -1 || ido == 1)
                {
#if 0
                  if (!vector_product (m, workd + ipntr[0] - 1,
                                       workd + ipntr[1] - 1))
                    break;
#else
                  // v = diag(1:100) * w
                  for (octave_idx_type i = 0; i < n; i++)
                    workd[i+ipntr[1]-1] = static_cast<double>(i) * workd[i+ipntr[0]-1];
#endif
                }
              else
                break;
            }

          if (info < 0)
            error ("dssimp: error %d in dsaupd", info);
          else
            {
              int rvec = 1;

              Matrix eig_vec (n, nev);
              double *z = eig_vec.fortran_vec ();

              ColumnVector eig_val (nev);
              double *d = eig_val.fortran_vec ();

              Array<int> s (ncv);
              int *sel = s.fortran_vec ();
              double sigma;

              F77_FUNC (dseupd, DSEUPD)
                (rvec, F77_CONST_CHAR_ARG2 ("A", 1), sel, d, z, n, sigma,
                 F77_CONST_CHAR_ARG2 (&bmat, 1), n,
                 F77_CONST_CHAR_ARG2 ((which.c_str()), 2),
                 nev, tol, resid, ncv, v, n, iparam, ipntr, workd, workl,
                 lworkl, info
                 F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(1) F77_CHAR_ARG_LEN(2));

              if (f77_exception_encountered)
                {
                  (*current_liboctave_error_handler)
                    ("dssimp: unrecoverable exception encountered in dseupd");
                  return retval;
                }

              retval(1) = eig_val;
              retval(2) = eig_vec;
            }
        }
    }  

  return retval;
}

      program dssimp
c
c     This example program is intended to illustrate the
c     simplest case of using ARPACK in considerable detail.  
c     This code may be used to understand basic usage of ARPACK
c     and as a template for creating an interface to ARPACK.  
c  
c     This code shows how to use ARPACK to find a few eigenvalues
c     (lambda) and corresponding eigenvectors (x) for the standard
c     eigenvalue problem:
c          
c                        A*x = lambda*x
c
c     where A is an n by n real symmetric matrix.
c
c     The main points illustrated here are
c
c        1) How to declare sufficient memory to find NEV
c           eigenvalues of largest magnitude.  Other options
c           are available.
c
c        2) Illustration of the reverse communication interface
c           needed to utilize the top level ARPACK routine DSAUPD
c           that computes the quantities needed to construct
c           the desired eigenvalues and eigenvectors(if requested).
c
c        3) How to extract the desired eigenvalues and eigenvectors
c           using the ARPACK routine DSEUPD.
c
c     The only thing that must be supplied in order to use this
c     routine on your problem is to change the array dimensions
c     appropriately, to specify WHICH eigenvalues you want to compute
c     and to supply a matrix-vector product
c
c                         w <-  Av
c
c     in place of the call to AV( ) below.
c
c     Once usage of this routine is understood, you may wish to explore
c     the other available options to improve convergence, to solve generalized
c     problems, etc.  Look at the file ex-sym.doc in DOCUMENTS directory.
c     This codes implements  
c
c\Example-1
c     ... Suppose we want to solve A*x = lambda*x in regular mode,
c         where A is derived from the central difference discretization
c         of the 2-dimensional Laplacian on the unit square with
c         zero Dirichlet boundary condition.
c     ... OP = A  and  B = I.
c     ... Assume "call av (n,x,y)" computes y = A*x
c     ... Use mode 1 of DSAUPD.
c
c\BeginLib
c
c\Routines called:
c     dsaupd  ARPACK reverse communication interface routine.
c     dseupd  ARPACK routine that returns Ritz values and (optionally)
c             Ritz vectors.
c     dnrm2   Level 1 BLAS that computes the norm of a vector.
c     daxpy   Level 1 BLAS that computes y <- alpha*x+y.
c
c\Author
c     Richard Lehoucq
c     Danny Sorensen
c     Chao Yang
c     Dept. of Computational &
c     Applied Mathematics
c     Rice University
c     Houston, Texas
c
c\SCCS Information: @(#)
c FILE: ssimp.F   SID: 2.6   DATE OF SID: 10/17/00   RELEASE: 2
c
c\Remarks
c     1. None
c
c\EndLib
c
c-----------------------------------------------------------------------
c
c     %------------------------------------------------------%
c     | Storage Declarations:                                |
c     |                                                      |
c     | The maximum dimensions for all arrays are            |
c     | set here to accommodate a problem size of            |
c     | N .le. MAXN                                          |
c     |                                                      |
c     | NEV is the number of eigenvalues requested.          |
c     |     See specifications for ARPACK usage below.       |
c     |                                                      |
c     | NCV is the largest number of basis vectors that will |
c     |     be used in the Implicitly Restarted Arnoldi      |
c     |     Process.  Work per major iteration is            |
c     |     proportional to N*NCV*NCV.                       |
c     |                                                      |
c     | You must set:                                        |
c     |                                                      |
c     | MAXN:   Maximum dimension of the A allowed.          |
c     | MAXNEV: Maximum NEV allowed.                         |
c     | MAXNCV: Maximum NCV allowed.                         |
c     %------------------------------------------------------%
c
      integer          maxn, maxnev, maxncv, ldv
      parameter       (maxn=256, maxnev=10, maxncv=25,
     $                 ldv=maxn )
c
c     %--------------%
c     | Local Arrays |
c     %--------------%
c
      Double precision
     &                 v(ldv,maxncv), workl(maxncv*(maxncv+8)),
     &                 workd(3*maxn), d(maxncv,2), resid(maxn),
     &                 ax(maxn)
      logical          select(maxncv)
      integer          iparam(11), ipntr(11)
c
c     %---------------%
c     | Local Scalars |
c     %---------------%
c
      character        bmat*1, which*2
      integer          ido, n, nev, ncv, lworkl, info, ierr,
     &                 j, nx, ishfts, maxitr, mode1, nconv
      logical          rvec
      Double precision      
     &                 tol, sigma
c
c     %------------%
c     | Parameters |
c     %------------%
c
      Double precision
     &                 zero
      parameter        (zero = 0.0D+0)
c  
c     %-----------------------------%
c     | BLAS & LAPACK routines used |
c     %-----------------------------%
c
      Double precision          
     &                 dnrm2
      external         dnrm2, daxpy
c
c     %--------------------%
c     | Intrinsic function |
c     %--------------------%
c
      intrinsic        abs
c
c     %-----------------------%
c     | Executable Statements |
c     %-----------------------%
c
c     %-------------------------------------------------%
c     | The following include statement and assignments |
c     | initiate trace output from the internal         |
c     | actions of ARPACK.  See debug.doc in the        |
c     | DOCUMENTS directory for usage.  Initially, the  |
c     | most useful information will be a breakdown of  |
c     | time spent in the various stages of computation |
c     | given by setting msaupd = 1.                    |
c     %-------------------------------------------------%
c
      include 'debug.h'
      ndigit = -3
      logfil = 6
      msgets = 0
      msaitr = 0
      msapps = 0
      msaupd = 1
      msaup2 = 0
      mseigt = 0
      mseupd = 0
c    
c     %-------------------------------------------------%
c     | The following sets dimensions for this problem. |
c     %-------------------------------------------------%
c
      nx = 10
      n = nx*nx
c
c     %-----------------------------------------------%
c     |                                               |
c     | Specifications for ARPACK usage are set       |
c     | below:                                        |
c     |                                               |
c     |    1) NEV = 4  asks for 4 eigenvalues to be   |  
c     |       computed.                               |
c     |                                               |
c     |    2) NCV = 20 sets the length of the Arnoldi |
c     |       factorization                           |
c     |                                               |
c     |    3) This is a standard problem              |
c     |         (indicated by bmat  = 'I')            |
c     |                                               |
c     |    4) Ask for the NEV eigenvalues of          |
c     |       largest magnitude                       |
c     |         (indicated by which = 'LM')           |
c     |       See documentation in DSAUPD for the     |
c     |       other options SM, LA, SA, LI, SI.       |
c     |                                               |
c     | Note: NEV and NCV must satisfy the following  |
c     | conditions:                                   |
c     |              NEV <= MAXNEV                    |
c     |          NEV + 1 <= NCV <= MAXNCV             |
c     %-----------------------------------------------%
c
      nev   = 4
      ncv   = 20
      bmat  = 'I'
      which = 'LM'
c
      if ( n .gt. maxn ) then
         print *, ' ERROR with _SSIMP: N is greater than MAXN '
         go to 9000
      else if ( nev .gt. maxnev ) then
         print *, ' ERROR with _SSIMP: NEV is greater than MAXNEV '
         go to 9000
      else if ( ncv .gt. maxncv ) then
         print *, ' ERROR with _SSIMP: NCV is greater than MAXNCV '
         go to 9000
      end if
c
c     %-----------------------------------------------------%
c     |                                                     |
c     | Specification of stopping rules and initial         |
c     | conditions before calling DSAUPD                    |
c     |                                                     |
c     | TOL  determines the stopping criterion.             |
c     |                                                     |
c     |      Expect                                         |
c     |           abs(lambdaC - lambdaT) < TOL*abs(lambdaC) |
c     |               computed   true                       |
c     |                                                     |
c     |      If TOL .le. 0,  then TOL <- macheps            |
c     |           (machine precision) is used.              |
c     |                                                     |
c     | IDO  is the REVERSE COMMUNICATION parameter         |
c     |      used to specify actions to be taken on return  |
c     |      from DSAUPD. (See usage below.)                |
c     |                                                     |
c     |      It MUST initially be set to 0 before the first |
c     |      call to DSAUPD.                                |
c     |                                                     |
c     | INFO on entry specifies starting vector information |
c     |      and on return indicates error codes            |
c     |                                                     |
c     |      Initially, setting INFO=0 indicates that a     |
c     |      random starting vector is requested to         |
c     |      start the ARNOLDI iteration.  Setting INFO to  |
c     |      a nonzero value on the initial call is used    |
c     |      if you want to specify your own starting       |
c     |      vector (This vector must be placed in RESID.)  |
c     |                                                     |
c     | The work array WORKL is used in DSAUPD as           |
c     | workspace.  Its dimension LWORKL is set as          |
c     | illustrated below.                                  |
c     |                                                     |
c     %-----------------------------------------------------%
c
      lworkl = ncv*(ncv+8)
      tol = zero
      info = 0
      ido = 0
c
c     %---------------------------------------------------%
c     | Specification of Algorithm Mode:                  |
c     |                                                   |
c     | This program uses the exact shift strategy        |
c     | (indicated by setting PARAM(1) = 1).              |
c     | IPARAM(3) specifies the maximum number of Arnoldi |
c     | iterations allowed.  Mode 1 of DSAUPD is used     |
c     | (IPARAM(7) = 1). All these options can be changed |
c     | by the user. For details see the documentation in |
c     | DSAUPD.                                           |
c     %---------------------------------------------------%
c
      ishfts = 1
      maxitr = 300
      mode1 = 1
c
      iparam(1) = ishfts
c                
      iparam(3) = maxitr
c                  
      iparam(7) = mode1
c
c     %------------------------------------------------%
c     | M A I N   L O O P (Reverse communication loop) |
c     %------------------------------------------------%
c
 10   continue
c
c        %---------------------------------------------%
c        | Repeatedly call the routine DSAUPD and take |
c        | actions indicated by parameter IDO until    |
c        | either convergence is indicated or maxitr   |
c        | has been exceeded.                          |
c        %---------------------------------------------%
c
         call dsaupd ( ido, bmat, n, which, nev, tol, resid,
     &                 ncv, v, ldv, iparam, ipntr, workd, workl,
     &                 lworkl, info )
c
         if (ido .eq. -1 .or. ido .eq. 1) then
c
c           %--------------------------------------%
c           | Perform matrix vector multiplication |
c           |              y <--- OP*x             |
c           | The user should supply his/her own   |
c           | matrix vector multiplication routine |
c           | here that takes workd(ipntr(1)) as   |
c           | the input, and return the result to  |
c           | workd(ipntr(2)).                     |
c           %--------------------------------------%
c
            call av2 (nx, workd(ipntr(1)), workd(ipntr(2)))
c
c           %-----------------------------------------%
c           | L O O P   B A C K to call DSAUPD again. |
c           %-----------------------------------------%
c
            go to 10
c
         end if
c
c     %----------------------------------------%
c     | Either we have convergence or there is |
c     | an error.                              |
c     %----------------------------------------%
c
      if ( info .lt. 0 ) then
c
c        %--------------------------%
c        | Error message. Check the |
c        | documentation in DSAUPD. |
c        %--------------------------%
c
         print *, ' '
         print *, ' Error with _saupd, info = ', info
         print *, ' Check documentation in _saupd '
         print *, ' '
c
      else
c
c        %-------------------------------------------%
c        | No fatal errors occurred.                 |
c        | Post-Process using DSEUPD.                |
c        |                                           |
c        | Computed eigenvalues may be extracted.    |  
c        |                                           |
c        | Eigenvectors may be also computed now if  |
c        | desired.  (indicated by rvec = .true.)    |
c        |                                           |
c        | The routine DSEUPD now called to do this  |
c        | post processing (Other modes may require  |
c        | more complicated post processing than     |
c        | mode1.)                                   |
c        |                                           |
c        %-------------------------------------------%
c          
          rvec = .true.
c
          call dseupd ( rvec, 'All', select, d, v, ldv, sigma,
     &         bmat, n, which, nev, tol, resid, ncv, v, ldv,
     &         iparam, ipntr, workd, workl, lworkl, ierr )
c
c         %----------------------------------------------%
c         | Eigenvalues are returned in the first column |
c         | of the two dimensional array D and the       |
c         | corresponding eigenvectors are returned in   |
c         | the first NCONV (=IPARAM(5)) columns of the  |
c         | two dimensional array V if requested.        |
c         | Otherwise, an orthogonal basis for the       |
c         | invariant subspace corresponding to the      |
c         | eigenvalues in D is returned in V.           |
c         %----------------------------------------------%
c
          if ( ierr .ne. 0) then
c
c            %------------------------------------%
c            | Error condition:                   |
c            | Check the documentation of DSEUPD. |
c            %------------------------------------%
c
             print *, ' '
             print *, ' Error with _seupd, info = ', ierr
             print *, ' Check the documentation of _seupd. '
             print *, ' '
c
          else
c
             nconv =  iparam(5)
             do 20 j=1, nconv
c
c               %---------------------------%
c               | Compute the residual norm |
c               |                           |
c               |   ||  A*x - lambda*x ||   |
c               |                           |
c               | for the NCONV accurately  |
c               | computed eigenvalues and  |
c               | eigenvectors.  (iparam(5) |
c               | indicates how many are    |
c               | accurate to the requested |
c               | tolerance)                |
c               %---------------------------%
c
                call av2(nx, v(1,j), ax)
                call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
                d(j,2) = dnrm2(n, ax, 1)
                d(j,2) = d(j,2) / abs(d(j,1))
c
 20          continue
c
c            %-----------------------------%
c            | Display computed residuals. |
c            %-----------------------------%
c
             call dmout(6, nconv, 2, d, maxncv, -6,
     &            'Ritz values and relative residuals')
          end if
c
c         %-------------------------------------------%
c         | Print additional convergence information. |
c         %-------------------------------------------%
c
          if ( info .eq. 1) then
             print *, ' '
             print *, ' Maximum number of iterations reached.'
             print *, ' '
          else if ( info .eq. 3) then
             print *, ' '
             print *, ' No shifts could be applied during implicit',
     &                ' Arnoldi update, try increasing NCV.'
             print *, ' '
          end if      
c
          print *, ' '
          print *, ' _SSIMP '
          print *, ' ====== '
          print *, ' '
          print *, ' Size of the matrix is ', n
          print *, ' The number of Ritz values requested is ', nev
          print *, ' The number of Arnoldi vectors generated',
     &             ' (NCV) is ', ncv
          print *, ' What portion of the spectrum: ', which
          print *, ' The number of converged Ritz values is ',
     &               nconv
          print *, ' The number of Implicit Arnoldi update',
     &             ' iterations taken is ', iparam(3)
          print *, ' The number of OP*x is ', iparam(9)
          print *, ' The convergence criterion is ', tol
          print *, ' '
c
      end if
c
c     %---------------------------%
c     | Done with program dssimp. |
c     %---------------------------%
c
 9000 continue
c
      end
c
      subroutine av2 (nx, v, w)
      integer nx, j
      double precision  v(nx*nx), w(nx*nx)
      do 10 j = 1, nx * nx
         w(j) = dble(j) * v(j)
 10   continue

      end

c
c ------------------------------------------------------------------
c     matrix vector subroutine
c
c     The matrix used is the 2 dimensional discrete Laplacian on unit
c     square with zero Dirichlet boundary condition.
c
c     Computes w <--- OP*v, where OP is the nx*nx by nx*nx block
c     tridiagonal matrix
c
c                  | T -I          |
c                  |-I  T -I       |
c             OP = |   -I  T       |
c                  |        ...  -I|
c                  |           -I T|
c
c     The subroutine TV is called to computed y<---T*x.
c
      subroutine av (nx, v, w)
      integer           nx, j, lo, n2
      Double precision
     &                  v(nx*nx), w(nx*nx), one, h2
      parameter         ( one = 1.0D+0 )
c
      call tv(nx,v(1),w(1))
      call daxpy(nx, -one, v(nx+1), 1, w(1), 1)
c
      do 10 j = 2, nx-1
         lo = (j-1)*nx
         call tv(nx, v(lo+1), w(lo+1))
         call daxpy(nx, -one, v(lo-nx+1), 1, w(lo+1), 1)
         call daxpy(nx, -one, v(lo+nx+1), 1, w(lo+1), 1)
  10  continue
c
      lo = (nx-1)*nx
      call tv(nx, v(lo+1), w(lo+1))
      call daxpy(nx, -one, v(lo-nx+1), 1, w(lo+1), 1)
c
c     Scale the vector w by (1/h^2), where h is the mesh size
c
      n2 = nx*nx
      h2 = one / dble((nx+1)*(nx+1))
      call dscal(n2, one/h2, w, 1)
      return
      end
c
c-------------------------------------------------------------------
      subroutine tv (nx, x, y)
c
      integer           nx, j
      Double precision
     &                  x(nx), y(nx), dd, dl, du
c
      Double precision
     &                  one, four
      parameter         (one = 1.0D+0, four = 4.0D+0)
c
c     Compute the matrix vector multiplication y<---T*x
c     where T is a nx by nx tridiagonal matrix with DD on the
c     diagonal, DL on the subdiagonal, and DU on the superdiagonal.
c    
c
      dd  = four
      dl  = -one
      du  = -one
c
      y(1) =  dd*x(1) + du*x(2)
      do 10 j = 2,nx-1
         y(j) = dl*x(j-1) + dd*x(j) + du*x(j+1)
 10   continue
      y(nx) =  dl*x(nx-1) + dd*x(nx)
      return
      end

Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

Jaroslav Hajek-2
On Sat, Dec 20, 2008 at 12:27 AM, David Bateman <[hidden email]> wrote:

> I've been looking at porting the eigs function to octave and basically have
> it finished.. However I have a bizarre issue in that
>
> eigs (diag(1:100),4,'LM')
>
> doesn't work due to an internal error in an ARPACK function.. I'd remember
> something like this and that I'd patch arpack to avoid it on my machine in
> the past. However, I thought I'd investigate the issue further and modified
> an ARPACK example function "dssimp" to treat the exact same case as above
> and no bug manifested..
>
> I then tried reimplementing dssimp  as an oct-file, trying to keep it close
> to the same functionality as the fortran code (ie arrays initialized to
> zero, etc), and the result was I saw the same bug in the DSAUPD function as
> previously. I've been staring at this code for hours and don't see what the
> issue is. I attach my modified dssimp.f and the code that duplicates it in
> octave with "dssimp(diag(1:100))"... Anyone want to confirm that
>
> gfortran -o dssimp dssimp.f -larpack -lblas -llapack
> ./dssimp
>
> works and gives the 4 largest eigenvalues [100, 99, 98, 97] as expected.
> whereas
>
> mkoctfile dssimp.cc -larpack
> octave --eval "dssimp(diag(1:100))"
>
> doesn't... Anyone want to look at the code and see if they can see an issue
> that I missed?
>
> Regards
> David
>
>
>
>
> --
> David Bateman                                [hidden email]
> 35 rue Gambetta                              +33 1 46 04 02 18 (Home)
> 92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)
>
>
>      program dssimp
> c
> c     This example program is intended to illustrate the
> c     simplest case of using ARPACK in considerable detail.
> c     This code may be used to understand basic usage of ARPACK
> c     and as a template for creating an interface to ARPACK.
> c
> c     This code shows how to use ARPACK to find a few eigenvalues
> c     (lambda) and corresponding eigenvectors (x) for the standard
> c     eigenvalue problem:
> c
> c                        A*x = lambda*x
> c
> c     where A is an n by n real symmetric matrix.
> c
> c     The main points illustrated here are
> c
> c        1) How to declare sufficient memory to find NEV
> c           eigenvalues of largest magnitude.  Other options
> c           are available.
> c
> c        2) Illustration of the reverse communication interface
> c           needed to utilize the top level ARPACK routine DSAUPD
> c           that computes the quantities needed to construct
> c           the desired eigenvalues and eigenvectors(if requested).
> c
> c        3) How to extract the desired eigenvalues and eigenvectors
> c           using the ARPACK routine DSEUPD.
> c
> c     The only thing that must be supplied in order to use this
> c     routine on your problem is to change the array dimensions
> c     appropriately, to specify WHICH eigenvalues you want to compute
> c     and to supply a matrix-vector product
> c
> c                         w <-  Av
> c
> c     in place of the call to AV( ) below.
> c
> c     Once usage of this routine is understood, you may wish to explore
> c     the other available options to improve convergence, to solve
> generalized
> c     problems, etc.  Look at the file ex-sym.doc in DOCUMENTS directory.
> c     This codes implements
> c
> c\Example-1
> c     ... Suppose we want to solve A*x = lambda*x in regular mode,
> c         where A is derived from the central difference discretization
> c         of the 2-dimensional Laplacian on the unit square with
> c         zero Dirichlet boundary condition.
> c     ... OP = A  and  B = I.
> c     ... Assume "call av (n,x,y)" computes y = A*x
> c     ... Use mode 1 of DSAUPD.
> c
> c\BeginLib
> c
> c\Routines called:
> c     dsaupd  ARPACK reverse communication interface routine.
> c     dseupd  ARPACK routine that returns Ritz values and (optionally)
> c             Ritz vectors.
> c     dnrm2   Level 1 BLAS that computes the norm of a vector.
> c     daxpy   Level 1 BLAS that computes y <- alpha*x+y.
> c
> c\Author
> c     Richard Lehoucq
> c     Danny Sorensen
> c     Chao Yang
> c     Dept. of Computational &
> c     Applied Mathematics
> c     Rice University
> c     Houston, Texas
> c
> c\SCCS Information: @(#)
> c FILE: ssimp.F   SID: 2.6   DATE OF SID: 10/17/00   RELEASE: 2
> c
> c\Remarks
> c     1. None
> c
> c\EndLib
> c
> c-----------------------------------------------------------------------
> c
> c     %------------------------------------------------------%
> c     | Storage Declarations:                                |
> c     |                                                      |
> c     | The maximum dimensions for all arrays are            |
> c     | set here to accommodate a problem size of            |
> c     | N .le. MAXN                                          |
> c     |                                                      |
> c     | NEV is the number of eigenvalues requested.          |
> c     |     See specifications for ARPACK usage below.       |
> c     |                                                      |
> c     | NCV is the largest number of basis vectors that will |
> c     |     be used in the Implicitly Restarted Arnoldi      |
> c     |     Process.  Work per major iteration is            |
> c     |     proportional to N*NCV*NCV.                       |
> c     |                                                      |
> c     | You must set:                                        |
> c     |                                                      |
> c     | MAXN:   Maximum dimension of the A allowed.          |
> c     | MAXNEV: Maximum NEV allowed.                         |
> c     | MAXNCV: Maximum NCV allowed.                         |
> c     %------------------------------------------------------%
> c
>      integer          maxn, maxnev, maxncv, ldv
>      parameter       (maxn=256, maxnev=10, maxncv=25,
>     $                 ldv=maxn )
> c
> c     %--------------%
> c     | Local Arrays |
> c     %--------------%
> c
>      Double precision
>     &                 v(ldv,maxncv), workl(maxncv*(maxncv+8)),
>     &                 workd(3*maxn), d(maxncv,2), resid(maxn),
>     &                 ax(maxn)
>      logical          select(maxncv)
>      integer          iparam(11), ipntr(11)
> c
> c     %---------------%
> c     | Local Scalars |
> c     %---------------%
> c
>      character        bmat*1, which*2
>      integer          ido, n, nev, ncv, lworkl, info, ierr,
>     &                 j, nx, ishfts, maxitr, mode1, nconv
>      logical          rvec
>      Double precision
>     &                 tol, sigma
> c
> c     %------------%
> c     | Parameters |
> c     %------------%
> c
>      Double precision
>     &                 zero
>      parameter        (zero = 0.0D+0)
> c
> c     %-----------------------------%
> c     | BLAS & LAPACK routines used |
> c     %-----------------------------%
> c
>      Double precision
>     &                 dnrm2
>      external         dnrm2, daxpy
> c
> c     %--------------------%
> c     | Intrinsic function |
> c     %--------------------%
> c
>      intrinsic        abs
> c
> c     %-----------------------%
> c     | Executable Statements |
> c     %-----------------------%
> c
> c     %-------------------------------------------------%
> c     | The following include statement and assignments |
> c     | initiate trace output from the internal         |
> c     | actions of ARPACK.  See debug.doc in the        |
> c     | DOCUMENTS directory for usage.  Initially, the  |
> c     | most useful information will be a breakdown of  |
> c     | time spent in the various stages of computation |
> c     | given by setting msaupd = 1.                    |
> c     %-------------------------------------------------%
> c
>      include 'debug.h'
>      ndigit = -3
>      logfil = 6
>      msgets = 0
>      msaitr = 0
>      msapps = 0
>      msaupd = 1
>      msaup2 = 0
>      mseigt = 0
>      mseupd = 0
> c
> c     %-------------------------------------------------%
> c     | The following sets dimensions for this problem. |
> c     %-------------------------------------------------%
> c
>      nx = 10
>      n = nx*nx
> c
> c     %-----------------------------------------------%
> c     |                                               |
> c     | Specifications for ARPACK usage are set       |
> c     | below:                                        |
> c     |                                               |
> c     |    1) NEV = 4  asks for 4 eigenvalues to be   |
> c     |       computed.                               |
> c     |                                               |
> c     |    2) NCV = 20 sets the length of the Arnoldi |
> c     |       factorization                           |
> c     |                                               |
> c     |    3) This is a standard problem              |
> c     |         (indicated by bmat  = 'I')            |
> c     |                                               |
> c     |    4) Ask for the NEV eigenvalues of          |
> c     |       largest magnitude                       |
> c     |         (indicated by which = 'LM')           |
> c     |       See documentation in DSAUPD for the     |
> c     |       other options SM, LA, SA, LI, SI.       |
> c     |                                               |
> c     | Note: NEV and NCV must satisfy the following  |
> c     | conditions:                                   |
> c     |              NEV <= MAXNEV                    |
> c     |          NEV + 1 <= NCV <= MAXNCV             |
> c     %-----------------------------------------------%
> c
>      nev   = 4
>      ncv   = 20
>      bmat  = 'I'
>      which = 'LM'
> c
>      if ( n .gt. maxn ) then
>         print *, ' ERROR with _SSIMP: N is greater than MAXN '
>         go to 9000
>      else if ( nev .gt. maxnev ) then
>         print *, ' ERROR with _SSIMP: NEV is greater than MAXNEV '
>         go to 9000
>      else if ( ncv .gt. maxncv ) then
>         print *, ' ERROR with _SSIMP: NCV is greater than MAXNCV '
>         go to 9000
>      end if
> c
> c     %-----------------------------------------------------%
> c     |                                                     |
> c     | Specification of stopping rules and initial         |
> c     | conditions before calling DSAUPD                    |
> c     |                                                     |
> c     | TOL  determines the stopping criterion.             |
> c     |                                                     |
> c     |      Expect                                         |
> c     |           abs(lambdaC - lambdaT) < TOL*abs(lambdaC) |
> c     |               computed   true                       |
> c     |                                                     |
> c     |      If TOL .le. 0,  then TOL <- macheps            |
> c     |           (machine precision) is used.              |
> c     |                                                     |
> c     | IDO  is the REVERSE COMMUNICATION parameter         |
> c     |      used to specify actions to be taken on return  |
> c     |      from DSAUPD. (See usage below.)                |
> c     |                                                     |
> c     |      It MUST initially be set to 0 before the first |
> c     |      call to DSAUPD.                                |
> c     |                                                     |
> c     | INFO on entry specifies starting vector information |
> c     |      and on return indicates error codes            |
> c     |                                                     |
> c     |      Initially, setting INFO=0 indicates that a     |
> c     |      random starting vector is requested to         |
> c     |      start the ARNOLDI iteration.  Setting INFO to  |
> c     |      a nonzero value on the initial call is used    |
> c     |      if you want to specify your own starting       |
> c     |      vector (This vector must be placed in RESID.)  |
> c     |                                                     |
> c     | The work array WORKL is used in DSAUPD as           |
> c     | workspace.  Its dimension LWORKL is set as          |
> c     | illustrated below.                                  |
> c     |                                                     |
> c     %-----------------------------------------------------%
> c
>      lworkl = ncv*(ncv+8)
>      tol = zero
>      info = 0
>      ido = 0
> c
> c     %---------------------------------------------------%
> c     | Specification of Algorithm Mode:                  |
> c     |                                                   |
> c     | This program uses the exact shift strategy        |
> c     | (indicated by setting PARAM(1) = 1).              |
> c     | IPARAM(3) specifies the maximum number of Arnoldi |
> c     | iterations allowed.  Mode 1 of DSAUPD is used     |
> c     | (IPARAM(7) = 1). All these options can be changed |
> c     | by the user. For details see the documentation in |
> c     | DSAUPD.                                           |
> c     %---------------------------------------------------%
> c
>      ishfts = 1
>      maxitr = 300
>      mode1 = 1
> c
>      iparam(1) = ishfts
> c
>      iparam(3) = maxitr
> c
>      iparam(7) = mode1
> c
> c     %------------------------------------------------%
> c     | M A I N   L O O P (Reverse communication loop) |
> c     %------------------------------------------------%
> c
>  10   continue
> c
> c        %---------------------------------------------%
> c        | Repeatedly call the routine DSAUPD and take |
> c        | actions indicated by parameter IDO until    |
> c        | either convergence is indicated or maxitr   |
> c        | has been exceeded.                          |
> c        %---------------------------------------------%
> c
>         call dsaupd ( ido, bmat, n, which, nev, tol, resid,
>     &                 ncv, v, ldv, iparam, ipntr, workd, workl,
>     &                 lworkl, info )
> c
>         if (ido .eq. -1 .or. ido .eq. 1) then
> c
> c           %--------------------------------------%
> c           | Perform matrix vector multiplication |
> c           |              y <--- OP*x             |
> c           | The user should supply his/her own   |
> c           | matrix vector multiplication routine |
> c           | here that takes workd(ipntr(1)) as   |
> c           | the input, and return the result to  |
> c           | workd(ipntr(2)).                     |
> c           %--------------------------------------%
> c
>            call av2 (nx, workd(ipntr(1)), workd(ipntr(2)))
> c
> c           %-----------------------------------------%
> c           | L O O P   B A C K to call DSAUPD again. |
> c           %-----------------------------------------%
> c
>            go to 10
> c
>         end if
> c
> c     %----------------------------------------%
> c     | Either we have convergence or there is |
> c     | an error.                              |
> c     %----------------------------------------%
> c
>      if ( info .lt. 0 ) then
> c
> c        %--------------------------%
> c        | Error message. Check the |
> c        | documentation in DSAUPD. |
> c        %--------------------------%
> c
>         print *, ' '
>         print *, ' Error with _saupd, info = ', info
>         print *, ' Check documentation in _saupd '
>         print *, ' '
> c
>      else
> c
> c        %-------------------------------------------%
> c        | No fatal errors occurred.                 |
> c        | Post-Process using DSEUPD.                |
> c        |                                           |
> c        | Computed eigenvalues may be extracted.    |
> c        |                                           |
> c        | Eigenvectors may be also computed now if  |
> c        | desired.  (indicated by rvec = .true.)    |
> c        |                                           |
> c        | The routine DSEUPD now called to do this  |
> c        | post processing (Other modes may require  |
> c        | more complicated post processing than     |
> c        | mode1.)                                   |
> c        |                                           |
> c        %-------------------------------------------%
> c
>          rvec = .true.
> c
>          call dseupd ( rvec, 'All', select, d, v, ldv, sigma,
>     &         bmat, n, which, nev, tol, resid, ncv, v, ldv,
>     &         iparam, ipntr, workd, workl, lworkl, ierr )
> c
> c         %----------------------------------------------%
> c         | Eigenvalues are returned in the first column |
> c         | of the two dimensional array D and the       |
> c         | corresponding eigenvectors are returned in   |
> c         | the first NCONV (=IPARAM(5)) columns of the  |
> c         | two dimensional array V if requested.        |
> c         | Otherwise, an orthogonal basis for the       |
> c         | invariant subspace corresponding to the      |
> c         | eigenvalues in D is returned in V.           |
> c         %----------------------------------------------%
> c
>          if ( ierr .ne. 0) then
> c
> c            %------------------------------------%
> c            | Error condition:                   |
> c            | Check the documentation of DSEUPD. |
> c            %------------------------------------%
> c
>             print *, ' '
>             print *, ' Error with _seupd, info = ', ierr
>             print *, ' Check the documentation of _seupd. '
>             print *, ' '
> c
>          else
> c
>             nconv =  iparam(5)
>             do 20 j=1, nconv
> c
> c               %---------------------------%
> c               | Compute the residual norm |
> c               |                           |
> c               |   ||  A*x - lambda*x ||   |
> c               |                           |
> c               | for the NCONV accurately  |
> c               | computed eigenvalues and  |
> c               | eigenvectors.  (iparam(5) |
> c               | indicates how many are    |
> c               | accurate to the requested |
> c               | tolerance)                |
> c               %---------------------------%
> c
>                call av2(nx, v(1,j), ax)
>                call daxpy(n, -d(j,1), v(1,j), 1, ax, 1)
>                d(j,2) = dnrm2(n, ax, 1)
>                d(j,2) = d(j,2) / abs(d(j,1))
> c
>  20          continue
> c
> c            %-----------------------------%
> c            | Display computed residuals. |
> c            %-----------------------------%
> c
>             call dmout(6, nconv, 2, d, maxncv, -6,
>     &            'Ritz values and relative residuals')
>          end if
> c
> c         %-------------------------------------------%
> c         | Print additional convergence information. |
> c         %-------------------------------------------%
> c
>          if ( info .eq. 1) then
>             print *, ' '
>             print *, ' Maximum number of iterations reached.'
>             print *, ' '
>          else if ( info .eq. 3) then
>             print *, ' '
>             print *, ' No shifts could be applied during implicit',
>     &                ' Arnoldi update, try increasing NCV.'
>             print *, ' '
>          end if
> c
>          print *, ' '
>          print *, ' _SSIMP '
>          print *, ' ====== '
>          print *, ' '
>          print *, ' Size of the matrix is ', n
>          print *, ' The number of Ritz values requested is ', nev
>          print *, ' The number of Arnoldi vectors generated',
>     &             ' (NCV) is ', ncv
>          print *, ' What portion of the spectrum: ', which
>          print *, ' The number of converged Ritz values is ',
>     &               nconv
>          print *, ' The number of Implicit Arnoldi update',
>     &             ' iterations taken is ', iparam(3)
>          print *, ' The number of OP*x is ', iparam(9)
>          print *, ' The convergence criterion is ', tol
>          print *, ' '
> c
>      end if
> c
> c     %---------------------------%
> c     | Done with program dssimp. |
> c     %---------------------------%
> c
>  9000 continue
> c
>      end
> c
>      subroutine av2 (nx, v, w)
>      integer nx, j
>      double precision  v(nx*nx), w(nx*nx)
>      do 10 j = 1, nx * nx
>         w(j) = dble(j) * v(j)
>  10   continue
>
>      end
>
> c
> c ------------------------------------------------------------------
> c     matrix vector subroutine
> c
> c     The matrix used is the 2 dimensional discrete Laplacian on unit
> c     square with zero Dirichlet boundary condition.
> c
> c     Computes w <--- OP*v, where OP is the nx*nx by nx*nx block
> c     tridiagonal matrix
> c
> c                  | T -I          |
> c                  |-I  T -I       |
> c             OP = |   -I  T       |
> c                  |        ...  -I|
> c                  |           -I T|
> c
> c     The subroutine TV is called to computed y<---T*x.
> c
>      subroutine av (nx, v, w)
>      integer           nx, j, lo, n2
>      Double precision
>     &                  v(nx*nx), w(nx*nx), one, h2
>      parameter         ( one = 1.0D+0 )
> c
>      call tv(nx,v(1),w(1))
>      call daxpy(nx, -one, v(nx+1), 1, w(1), 1)
> c
>      do 10 j = 2, nx-1
>         lo = (j-1)*nx
>         call tv(nx, v(lo+1), w(lo+1))
>         call daxpy(nx, -one, v(lo-nx+1), 1, w(lo+1), 1)
>         call daxpy(nx, -one, v(lo+nx+1), 1, w(lo+1), 1)
>  10  continue
> c
>      lo = (nx-1)*nx
>      call tv(nx, v(lo+1), w(lo+1))
>      call daxpy(nx, -one, v(lo-nx+1), 1, w(lo+1), 1)
> c
> c     Scale the vector w by (1/h^2), where h is the mesh size
> c
>      n2 = nx*nx
>      h2 = one / dble((nx+1)*(nx+1))
>      call dscal(n2, one/h2, w, 1)
>      return
>      end
> c
> c-------------------------------------------------------------------
>      subroutine tv (nx, x, y)
> c
>      integer           nx, j
>      Double precision
>     &                  x(nx), y(nx), dd, dl, du
> c
>      Double precision
>     &                  one, four
>      parameter         (one = 1.0D+0, four = 4.0D+0)
> c
> c     Compute the matrix vector multiplication y<---T*x
> c     where T is a nx by nx tridiagonal matrix with DD on the
> c     diagonal, DL on the subdiagonal, and DU on the superdiagonal.
> c
> c
>      dd  = four
>      dl  = -one
>      du  = -one
> c
>      y(1) =  dd*x(1) + du*x(2)
>      do 10 j = 2,nx-1
>         y(j) = dl*x(j-1) + dd*x(j) + du*x(j+1)
>  10   continue
>      y(nx) =  dl*x(nx-1) + dd*x(nx)
>      return
>      end
>
>
>

1. rather than your MY_LOCAL_BUFFER, you can use the pre-defined
OCTAVE_LOCAL_BUFFER_INIT.

2. the following snippet:
                  for (octave_idx_type i = 0; i < n; i++)
                    workd[i+ipntr[1]-1] = static_cast<double>(i) * workd[i+ipntr[0]-1];
is obviously not multiplying by diag(1:100), but by diag(0:99).

3.
              retval(1) = eig_val;
              retval(2) = eig_vec;

should read

              retval(0) = eig_val;
              retval(1) = eig_vec;

after correcting these problems, I get the correct answer:
ans =

    97.000
    98.000
    99.000
   100.000

regards


--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

David Bateman-2
Jaroslav Hajek wrote:

> 1. rather than your MY_LOCAL_BUFFER, you can use the pre-defined
> OCTAVE_LOCAL_BUFFER_INIT.
>
> 2. the following snippet:
>  for (octave_idx_type i = 0; i < n; i++)
>    workd[i+ipntr[1]-1] = static_cast<double>(i) * workd[i+ipntr[0]-1];
> is obviously not multiplying by diag(1:100), but by diag(0:99).
>
> 3.
>      retval(1) = eig_val;
>      retval(2) = eig_vec;
>
> should read
>
>      retval(0) = eig_val;
>      retval(1) = eig_vec;
>
> after correcting these problems, I get the correct answer:
> ans =
>
>     97.000
>     98.000
>     99.000
>    100.000
>
> regards
>
>
>  

Unfortunately, that doesn't work for me as it never gets to point 3) and
in any case point 2) would just mean the eigenvalues that are found
should be [99, 98, 97, 96]... The code has an internal error in the
arpack function dsaupd (error number -9999).. However, as it seems to
work for you the question is why it works for you and not me...

Basically, I can only see compiler/linker issues as a reason at the
moment.. I compiled with gfortran version 4.3.2 and the flags

FFLAGS  = -O2 -fPIC

and linked dssimp.oct against the shared version of the library. What
did you do differently?

D.


--
David Bateman                                [hidden email]
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)

Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

Jaroslav Hajek-2
On Sat, Dec 20, 2008 at 5:23 PM, David Bateman <[hidden email]> wrote:

> Jaroslav Hajek wrote:
>>
>> 1. rather than your MY_LOCAL_BUFFER, you can use the pre-defined
>> OCTAVE_LOCAL_BUFFER_INIT.
>>
>> 2. the following snippet:
>>                  for (octave_idx_type i = 0; i < n; i++)
>>                    workd[i+ipntr[1]-1] = static_cast<double>(i) *
>> workd[i+ipntr[0]-1];
>> is obviously not multiplying by diag(1:100), but by diag(0:99).
>>
>> 3.
>>              retval(1) = eig_val;
>>              retval(2) = eig_vec;
>>
>> should read
>>
>>              retval(0) = eig_val;
>>              retval(1) = eig_vec;
>>
>> after correcting these problems, I get the correct answer:
>> ans =
>>
>>    97.000
>>    98.000
>>    99.000
>>   100.000
>>
>> regards
>>
>>
>>
>
> Unfortunately, that doesn't work for me as it never gets to point 3) and in
> any case point 2) would just mean the eigenvalues that are found should be
> [99, 98, 97, 96]... The code has an internal error in the arpack function
> dsaupd (error number -9999).. However, as it seems to work for you the
> question is why it works for you and not me...
>
> Basically, I can only see compiler/linker issues as a reason at the moment..
> I compiled with gfortran version 4.3.2 and the flags
>
> FFLAGS  = -O2 -fPIC
>
> and linked dssimp.oct against the shared version of the library. What did
> you do differently?
>

I used Intel Fortran with -O -fPIC, and I only built a static library
(where did you get export symbols table? or did you export
everything?), though I was not sure whether the linking will work
correctly (note that ARPACK should be linked to LAPACK v2 routines).
I can try with gfortran, but no sooner than after new year's eve.
In any case, if we include ARPACK, my suggestion is to work on things
a bit to investigate and remove the dependency on LAPACK v2, because
it is, IMHO, almost bound to cause trouble if ARPACK is included
within libcruft. Of course I'm up for the job, otherwise such a
suggestion would be pretty moot.
I'm not sure whether ARPACK is maintained anymore (yes, I know there
has been a license change recently, but maybe it was just because the
authors have no plans for further development).

Why it doesn't work for you is a good question. Maybe there's a buffer
overrun that in my case doesn't do any damage by chance? (Just a wild
guess).

regards

--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

John W. Eaton
Administrator
On 20-Dec-2008, Jaroslav Hajek wrote:

| I can try with gfortran, but no sooner than after new year's eve.
| In any case, if we include ARPACK, my suggestion is to work on things
| a bit to investigate and remove the dependency on LAPACK v2, because
| it is, IMHO, almost bound to cause trouble if ARPACK is included
| within libcruft. Of course I'm up for the job, otherwise such a
| suggestion would be pretty moot.
| I'm not sure whether ARPACK is maintained anymore (yes, I know there
| has been a license change recently, but maybe it was just because the
| authors have no plans for further development).

Are we going to include arpack sources in Octave?  I'd prefer to leave
it as an external dependency unless there is some good reason to do
otherwise.

| Why it doesn't work for you is a good question. Maybe there's a buffer
| overrun that in my case doesn't do any damage by chance? (Just a wild
| guess).

Both programs seem to work correctly for me.  I also thought it might
be a buffer overrun that was only showing up on some systems (32- vs
64-bit pointers, for example) so I tried it on 32- and 64-bit sytems
and it worked on both.  I also tried running it with valgrind and
there were no complaints for the call to dssimp.  I'm using Octave
3.0.1 (64-bit system) and Octave 3.0.0 (32-bit system) on a Debian
system with arpack from the Debian package.

jwe
Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

Jaroslav Hajek-2
On Sat, Dec 20, 2008 at 6:15 PM, John W. Eaton <[hidden email]> wrote:

> On 20-Dec-2008, Jaroslav Hajek wrote:
>
> | I can try with gfortran, but no sooner than after new year's eve.
> | In any case, if we include ARPACK, my suggestion is to work on things
> | a bit to investigate and remove the dependency on LAPACK v2, because
> | it is, IMHO, almost bound to cause trouble if ARPACK is included
> | within libcruft. Of course I'm up for the job, otherwise such a
> | suggestion would be pretty moot.
> | I'm not sure whether ARPACK is maintained anymore (yes, I know there
> | has been a license change recently, but maybe it was just because the
> | authors have no plans for further development).
>
> Are we going to include arpack sources in Octave?  I'd prefer to leave
> it as an external dependency unless there is some good reason to do
> otherwise.
>

The only reason is probably availability - if a GNU/Linux distro
doesn't provide ARPACK as a package, then it will probably provide
Octave without eigs unless someone volunteers to maintain ARPACK for
that distro. It's not a strong reason, probably.

> | Why it doesn't work for you is a good question. Maybe there's a buffer
> | overrun that in my case doesn't do any damage by chance? (Just a wild
> | guess).
>
> Both programs seem to work correctly for me.  I also thought it might
> be a buffer overrun that was only showing up on some systems (32- vs
> 64-bit pointers, for example) so I tried it on 32- and 64-bit sytems
> and it worked on both.  I also tried running it with valgrind and
> there were no complaints for the call to dssimp.  I'm using Octave
> 3.0.1 (64-bit system) and Octave 3.0.0 (32-bit system) on a Debian
> system with arpack from the Debian package.
>
> jwe
>



--
RNDr. Jaroslav Hajek
computing expert
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz
Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

David Bateman-2
Jaroslav Hajek wrote:

> On Sat, Dec 20, 2008 at 6:15 PM, John W. Eaton <[hidden email]> wrote:
>  
>> Are we going to include arpack sources in Octave?  I'd prefer to leave
>> it as an external dependency unless there is some good reason to do
>> otherwise.
>>
>>    
>
> The only reason is probably availability - if a GNU/Linux distro
> doesn't provide ARPACK as a package, then it will probably provide
> Octave without eigs unless someone volunteers to maintain ARPACK for
> that distro. It's not a strong reason, probably.
>
>  

Yes the ARPACK website no longer seems to be maintained. Viral Shah sent
me a whole set of patches for ARPACK in 2006 that he and Chao Yang (a
co-author of ARPACK) wrote. The link to the website these lived on
however seems to be broken..

I cc'ed both Chao and and Viral Shah at the addresses I had with the
hope of getting through to them and see what their thoughts are about
future support of ARPACK now that the license has changed to a two
clause BSD license. This gives them the possibility of publicly
releasing their patches.  If they are unwilling to support a download
site for these patches, maybe a distribution such as debian might be
willing to host the patches directly..

I agree  with John that it is best that Octave doesn't host its own copy
of arpack however.

Regards
David

--
David Bateman                                [hidden email]
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)

Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

David Bateman-2
In reply to this post by John W. Eaton
John W. Eaton wrote:

> | Why it doesn't work for you is a good question. Maybe there's a buffer
> | overrun that in my case doesn't do any damage by chance? (Just a wild
> | guess).
>
> Both programs seem to work correctly for me.  I also thought it might
> be a buffer overrun that was only showing up on some systems (32- vs
> 64-bit pointers, for example) so I tried it on 32- and 64-bit sytems
> and it worked on both.  I also tried running it with valgrind and
> there were no complaints for the call to dssimp.  I'm using Octave
> 3.0.1 (64-bit system) and Octave 3.0.0 (32-bit system) on a Debian
> system with arpack from the Debian package.
>  
ok, rinse and repeat...

I retried with various compilation and linking options and firstly it
appears that if I link to the static library everything works fine.
However if I have FFLAGS of "-fPIC -O2" and link arpack with

g++ -shared -Wl,-soname -Wl,libarpack.so -o libarpack.so  SRC/*.o
UTIL/*.o -llapack -lblas

as I would have thought was necessary, I have problems.. In fact I can't
seem to find any shared library version that I build that does work....

D.






> jwe
>
>  


--
David Bateman                                [hidden email]
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)

Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

David Bateman-2
David Bateman wrote:

> John W. Eaton wrote:
>> | Why it doesn't work for you is a good question. Maybe there's a buffer
>> | overrun that in my case doesn't do any damage by chance? (Just a wild
>> | guess).
>>
>> Both programs seem to work correctly for me.  I also thought it might
>> be a buffer overrun that was only showing up on some systems (32- vs
>> 64-bit pointers, for example) so I tried it on 32- and 64-bit sytems
>> and it worked on both.  I also tried running it with valgrind and
>> there were no complaints for the call to dssimp.  I'm using Octave
>> 3.0.1 (64-bit system) and Octave 3.0.0 (32-bit system) on a Debian
>> system with arpack from the Debian package.
>>  
> ok, rinse and repeat...
>
> I retried with various compilation and linking options and firstly it
> appears that if I link to the static library everything works fine.
> However if I have FFLAGS of "-fPIC -O2" and link arpack with
>
> g++ -shared -Wl,-soname -Wl,libarpack.so -o libarpack.so  SRC/*.o
> UTIL/*.o -llapack -lblas
>
> as I would have thought was necessary, I have problems.. In fact I
> can't seem to find any shared library version that I build that does
> work....
>
> D.
>

Ok, it seems that if I use the version of arpack from

http://mathema.tician.de/dl/software/arpack-autotools

that contains a few bug fixes and a build with libtool then I can get
eigs to work correctly... The reason I was building ARPACK myself in the
first place was that the version that seems to bee distributed with
debian isn't built with gfortran  with the version in testing.. Will
this still be the case in future versions of arpack on debian?

Regards
David

--
David Bateman                                [hidden email]
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)

Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

David Bateman-2
David Bateman wrote:

> David Bateman wrote:
>> John W. Eaton wrote:
>>> | Why it doesn't work for you is a good question. Maybe there's a
>>> buffer
>>> | overrun that in my case doesn't do any damage by chance? (Just a wild
>>> | guess).
>>>
>>> Both programs seem to work correctly for me.  I also thought it might
>>> be a buffer overrun that was only showing up on some systems (32- vs
>>> 64-bit pointers, for example) so I tried it on 32- and 64-bit sytems
>>> and it worked on both.  I also tried running it with valgrind and
>>> there were no complaints for the call to dssimp.  I'm using Octave
>>> 3.0.1 (64-bit system) and Octave 3.0.0 (32-bit system) on a Debian
>>> system with arpack from the Debian package.
>>>  
>> ok, rinse and repeat...
>>
>> I retried with various compilation and linking options and firstly it
>> appears that if I link to the static library everything works fine.
>> However if I have FFLAGS of "-fPIC -O2" and link arpack with
>>
>> g++ -shared -Wl,-soname -Wl,libarpack.so -o libarpack.so  SRC/*.o
>> UTIL/*.o -llapack -lblas
>>
>> as I would have thought was necessary, I have problems.. In fact I
>> can't seem to find any shared library version that I build that does
>> work....
>>
>> D.
>>
>
> Ok, it seems that if I use the version of arpack from
>
> http://mathema.tician.de/dl/software/arpack-autotools
>
> that contains a few bug fixes and a build with libtool then I can get
> eigs to work correctly... The reason I was building ARPACK myself in
> the first place was that the version that seems to bee distributed
> with debian isn't built with gfortran  with the version in testing..
> Will this still be the case in future versions of arpack on debian?
>
> Regards
> David
>
As this builds correctly and passes make check with the libtool version
of the shared library, I pushed my patch to savannah..

D.


--
David Bateman                                [hidden email]
35 rue Gambetta                              +33 1 46 04 02 18 (Home)
92100 Boulogne-Billancourt FRANCE            +33 6 72 01 06 33 (Mob)

Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

Thomas Weber-8
In reply to this post by David Bateman-2
On Tue, Dec 23, 2008 at 12:44:36AM +0100, David Bateman wrote:
> David Bateman wrote:
> Ok, it seems that if I use the version of arpack from
>
> http://mathema.tician.de/dl/software/arpack-autotools
>
> that contains a few bug fixes and a build with libtool then I can get  
> eigs to work correctly... The reason I was building ARPACK myself in the  
> first place was that the version that seems to bee distributed with  
> debian isn't built with gfortran  with the version in testing..

I'm not sure I'm parsing your sentence correctly, but arpack in Debian
testing is compiled with gfortran-4.3.

        Thomas
Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

David Bateman

Thomas Weber-8 wrote
On Tue, Dec 23, 2008 at 12:44:36AM +0100, David Bateman wrote:
> David Bateman wrote:
> Ok, it seems that if I use the version of arpack from
>
> http://mathema.tician.de/dl/software/arpack-autotools
>
> that contains a few bug fixes and a build with libtool then I can get  
> eigs to work correctly... The reason I was building ARPACK myself in the  
> first place was that the version that seems to bee distributed with  
> debian isn't built with gfortran  with the version in testing..

I'm not sure I'm parsing your sentence correctly, but arpack in Debian
testing is compiled with gfortran-4.3.
Is it? For my machine, that is using "debain testing", when I installed "arpack2-dev" it seemed it was compiled with g95 so I had problems linking to it..

Ok, I'll try again and see if I can get the debian distributed version of arpack to work for me...

D.
Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

Thomas Weber-8
On Tue, Dec 23, 2008 at 05:00:00AM -0800, dbateman wrote:

>
>
>
> Thomas Weber-8 wrote:
> >
> > On Tue, Dec 23, 2008 at 12:44:36AM +0100, David Bateman wrote:
> >> David Bateman wrote:
> >> Ok, it seems that if I use the version of arpack from
> >>
> >> http://mathema.tician.de/dl/software/arpack-autotools
> >>
> >> that contains a few bug fixes and a build with libtool then I can get  
> >> eigs to work correctly... The reason I was building ARPACK myself in the  
> >> first place was that the version that seems to bee distributed with  
> >> debian isn't built with gfortran  with the version in testing..
> >
> > I'm not sure I'm parsing your sentence correctly, but arpack in Debian
> > testing is compiled with gfortran-4.3.
> >
> >
>
> Is it? For my machine, that is using "debain testing", when I installed
> "arpack2-dev" it seemed it was compiled with g95 so I had problems linking
> to it..

According to
http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=422297
g95 is not in Debian. There is however a slight chance that the
maintainer has it installed on his machine. What architecture is your
machine (AMD64 or i386)?

        Thomas

Reply | Threaded
Open this post in threaded view
|

Re: porting arapck code to octave

David Bateman

Thomas Weber-8 wrote
According to
http://bugs.debian.org/cgi-bin/bugreport.cgi?bug=422297
g95 is not in Debian. There is however a slight chance that the
maintainer has it installed on his machine. What architecture is your
machine (AMD64 or i386)?

        Thomas
In fact at the moment I'm using an Intel Atom, though its particularly slow for compiles.. I an i386.. As I said I'll try again and see what results I get with the debian shared libraries of arpack..

D.