pjd.f

a software code for computing selected eigenvalues of large sparse symmetric matrices

*              INFO=0 if normal termination. 
*              INFO>0 if allowed maximum number of matrix vector 
*                        multplications performed without finding all
*                        wanted eigenvalues & eigenvectors. 
*              INFO<0 if an error occurred - see printed output for
*                        details
*              (INFO=-54321: the computation of the preconditioner
*                        failed).
*
*  GAP     (output) DOUBLE PRECISION
*          The estimated distance between the set of NEIG computed
*          eigenvalues and the remaining part of the spectrum; may be
*          inaccurate.
*
*  ===========================================================
*
       integer keep, IJOB, NDX1, NDX2
*
        keep=ICNTL(1)
        ICNTL(1)=1
        IJOB=0
*
        CALL PJDRVCOM( N,a,ja,ia, EIGS, RES, X, LX, NEIG, 
     $                       SIGMA, ISEARCH, NINIT, MADSPACE, ITER, TOL,
     $                       SHIFT, DROPTOL, MEM, ICNTL,
     $                       IJOB, NDX1, NDX2, IPRINT, INFO, GAP)
*
       ICNTL(1)=keep
*
      RETURN
      END
*
*======================================================================*
*
*======================================================================*
*
       SUBROUTINE PJDREVCOM( N,a,ja,ia, EIGS, RES, X, LX, NEIG, 
     $                       SIGMA, ISEARCH, NINIT, MADSPACE, ITER, TOL,
     $                       SHIFT, DROPTOL, MEM, ICNTL,
     $                       IJOB, NDX1, NDX2, IPRINT, INFO, GAP)
      implicit none  
*
*     .. Scalar Arguments ..
      integer            N, LX, NEIG, ISEARCH, NINIT, MADSPACE
      integer            ITER, ICNTL(4), IJOB, NDX1, NDX2, IPRINT, INFO
      DOUBLE PRECISION   SIGMA, TOL, SHIFT, DROPTOL, MEM, GAP
*     ..
*     .. Array Arguments ..
      integer            ja(*), ia(*)
      DOUBLE PRECISION   a(*), EIGS(*), RES(*), X(*)
*     ..
*
*  Arguments
*  =========
*
*    N,A,JA,IA,                       |   see comments in subroutine PJD
*
*       In addition:
*
*            A,JA,IA need not to define exactly the matrix whose
*               eigenvalues are wanted; instead it may be some 
*               approximation; this is consistent because the matrix 
*               passed to PJDREVCOM is only used to define a
*               preconditioner, whereas matrix vector multiplications
*               are performed by a user provided routine via the reverse
*               communication protocol (see below).
*
*            A,JA,IA should be compliant with the format described in PJD
*            However: 
*                  zero diagonal entries need not to be referenced in
*                        the structure (although they may);
*                  if a diagonal preconditioning is wanted, one
*                        should set ja(1) negative;
*                        if ja(1)<0, a(1),...,a(N) is supposed
*                        to carry the diagonal of the matrix;
*                        then, ja does not need to have a length greater
*                        than 1, ia is neither referenced, and a,ja,ia
*                        are unchanged on output. 
*
*
*    EIGS, RES, X, LX, NEIG,           |
*    SIGMA, ISEARCH, NINIT, MADSPACE,  |  
*    ITER, TOL, SHIFT, DROPTOL, MEM    |  see comments in subroutine PJD
*                                      |
*    IPRINT, INFO, GAP                 |
*
*
* ICNTL    (input/output) integer
*          some control parameters
*          ICNTL(1) should be set to zero (default value), except if X 
*                   overwrites the arrays in A,JA,IA, in which case one
*                   should set ICNTL(1)=2 (this tells that the matrix
*                   cannot be refactored once the eigenvalue computation
*                   started)
*          other entries in ICNTL: see comments in subroutine PJD
*
*
*  IJOB    (input/output) integer. 
*          Used to communicate job code between the levels.
*          Input: one should use IJOB=0 on the first call, 
*                 and leave IJOB unchanged on subsequent calls
*          Output:
*               IJOB=0: work done - terminate
*               IJOB=1: compute X(NDX2:NDX2+N-1)= A*X(NDX1:NDX1+N-1)
*                       (call to matrix vector multiplication routine:
*                        MATVEC) and return to PJDREVCOM leaving IJOB
*                       (and other parameters) unchanged.
*                       
*  NDX1    (output) integer.
*  NDX2    Indicate indices into X() for the needed MATVEC when IJOB=1.
*
*  ============================================================
*
       CALL PJDRVCOM( N,a,ja,ia, EIGS, RES, X, LX, NEIG, 
     $                       SIGMA, ISEARCH, NINIT, MADSPACE, ITER, TOL,
     $                       SHIFT, DROPTOL, MEM, ICNTL,
     $                       IJOB, NDX1, NDX2, IPRINT, INFO, GAP)
      RETURN
      END
*
*======================================================================*
*
*======================================================================*
*
      SUBROUTINE JDREVCOM( N, EIGS, RES, X, LX, NEIG, 
     $                       SIGMA, ISEARCH, NINIT, MADSPACE, ITER, TOL,
     $                       IJOB, NDX1, NDX2, IPRINT, INFO, GAP)
      implicit none  
*
*     .. Scalar Arguments ..
      integer            N, LX, NEIG, ISEARCH, NINIT, MADSPACE
      integer            ITER, IJOB, NDX1, NDX2, IPRINT, INFO
      DOUBLE PRECISION   SIGMA, TOL, GAP
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   EIGS( *), RES(*), X( * )
*     ..
*
*  Arguments
*  =========
*
*    N,EIGS, RES, X, LX, NEIG,        |
*    SIGMA, ISEARCH, NINIT, MADSPACE, |  
*    ITER, TOL                        |   see comments in subroutine PJD
*                                     |
*    IPRINT, INFO, GAP                |
*
*
*  IJOB    (input/output) integer. 
*          Used to communicate job code between the levels.
*          Input: one should use IJOB=0 on the first call, and leave
*                 IJOB unchanged on subsequent calls
*          Output:
*               IJOB=0: work done - terminate
*               IJOB=1: compute X(NDX2:NDX2+N-1)= A*X(NDX1:NDX1+N-1) 
*                       (call to matrix vector multiplication routine:
*                        MATVEC) and return to PJDREVCOM leaving IJOB
*                       (and other parameters) unchanged.
*               IJOB=2: solve B*X(NDX1:NDX1+N-1) = X(NDX2:NDX2+N-1) and
*                       return (call to preconditioner solve routine: 
*                       PSOLVE) and return to PJDREVCOM leaving IJOB 
*                       (and other parameters) unchanged.
*
*  NDX1    (output) integer.
*  NDX2    Indicate indices into X() for the needed MATVEC when IJOB=1,
*          or for the needed PSOLVE when IJOB=2.
*
*  ============================================================
*
      CALL       JDRVCOM ( N, EIGS, RES, X, LX, NEIG, 
     $                       SIGMA, ISEARCH, NINIT, MADSPACE, ITER, TOL,
     $                       IJOB, NDX1, NDX2, IPRINT, INFO, GAP)
      RETURN
      END
*======================================================================*
*
*======================================================================*
*
      SUBROUTINE PJDCLEANUP
      implicit none
      double precision  timefact,shift0,droptol0,diagmin,shiftmax,
     +                  memrequested,memused,slightlyless,toldiv,
     +                  condest0
      INTEGER*8         IPparam,IPPREC,IPdiag
      integer           IPnlev,factdgl,factspd,IUNIT,PR,
     +                  prvdr
      logical           reenterfirsttime
      common/pjdinitpjd/timefact,shift0,droptol0,diagmin,shiftmax,
     +                  memrequested,memused,slightlyless,toldiv,
     +                  IPparam,IPPREC,IPdiag,IPnlev,factdgl,factspd,
     +                  prvdr,IUNIT,PR,condest0,reenterfirsttime
      integer ilumem  
      common/ilupackmem/ilumem

c     make sure that PJDCLEANUP is not really executed before PJD/
c     PJDREVCOM have been called for the first time
      if (ilumem.ne.-1) then
         if (factdgl.gt.0) then
            call DGLPRECDELETE(IPdiag)
         else if (factdgl.lt.0) then
            call dsymamgdelete(IPparam,IPPREC,IPnlev,0)
         end if
c        ensure that twice calling PJDCLEANUP is captured
         ilumem=-1
      end if
      factdgl=0
c
      return
      end
*
*
*=============================================================================*