pjd.f

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

*======================================================================*
*
       SUBROUTINE PJD( N,a,ja,ia, EIGS, RES, X, LX, NEIG, 
     $                       SIGMA, ISEARCH, NINIT, MADSPACE, ITER, TOL,
     $                       SHIFT, DROPTOL, MEM, ICNTL, IPRINT, INFO, 
     $                       GAP)
      implicit none  
*
*     .. Scalar Arguments ..
      INTEGER            N, LX, NEIG, ISEARCH, NINIT, MADSPACE
      INTEGER            ITER, ICNTL(4), IPRINT, INFO
      DOUBLE PRECISION   SIGMA, TOL, SHIFT, DROPTOL, MEM, GAP
*     ..
*     .. Array Arguments ..
      INTEGER            JA(*), IA(*)
      DOUBLE PRECISION   A(*), EIGS(*), RES(*), X(*)
*     ..
*
*  Arguments
*  =========
*
*  N       (input) INTEGER.
*          The dimension of the matrix. Should be larger than 1.
*
*  A       (input/output) numerical values
*  IA      (input/output) pointers for every column
*  JA      (input/output) column indices, if JA(1)<0, then A is assumed 
*          to be diagonal.
*
*          Detailed description of A,IA,JA
* 
*              Note first that no more than the diagonal part of A and 
*              EITHER values from the strict upper triangular part OR 
*              values from the strict lower triangular part are needed 
*              to define the full matrix. 
*              PJD/PJDREVCOM ASSUMES THAT EACH OFF-DIAGONAL ENTRY IS
*              REFERENCED ONLY ONCE, in either the upper triangular part
*              or in the lower triangular part.
*              (Thus, for instance, a,ia,ja may contain either the upper
*              triangular part or the lower triangular part - in each
*              case including diagonal).
*
*              PJD ASSUMES IN ADDITION THAT ALL DIAGONAL ELEMENTS ARE
*              REFERENCED, EVEN WHEN SOME OF THEM ARE EQUAL TO ZERO.
*
*              OTHERWISE THE CODE WILL NOT RUN PROPERLY
*
*              On input, IA(I), I=1,...,N+1 refers to the physical start
*              of row I. In this case the entries of row I are located 
*              in A(K), where K=IA(I),...IA(I+1)-1. JA(K) carries the 
*              associated column indices. According what is written
*              above, PJD assumes that some of these JA(K) (for 
*              IA(I)<= K < IA(I+1) ) is equal to I with corresponding 
*              A(K) carrying the value of the diagonal element, possibly
*              equal to zero.
*
*              A,IA,JA are "output" parameters because on exit the 
*              entries of each row may occur in a different order (The
*              matrix is mathematically the same, but stored in 
*              different way).
*              
*  EIGS    (input/output) DOUBLE PRECISION array, dimension NEIG.
*          On input, eigenvalue estimates corresponding to provided
*              initial guesses (EIGS(i) corresponds to approximate 
*              eigenvector number i); used only if NINIT>(MADSPACE+1)/2
*              to make sure that initial approximate eigenvectors are
*              processed in the right order. 
*              Sorting is skipped if EIGS(1)=EIGS(2)= ... =EIGS(NINIT).
*              Then, if NINIT > (MADSPACE+1)/2 , initial approximations
*              should be in stored in increasing order of eigenvalues if
*              ISEARCH.LE.1 , or in increasing distance of eigenvalues
*              to SIGMA if ISEARCH.GE.2
*          On output, eigenvalues as they are computed
*              (EIGS(i) corresponds to approximate eigenvector number i).
*
*  RES     (output) DOUBLE PRECISION array, dimension NEIG.
*          Residual norms: RES(i)=|| A*x(i)-EIGS(i)*x(i) ||, 
*              where A is the matrix, || is the two norm, and x(i) is
*              the approximate eigenvector number i .
*
*  X       (input/output+workspace) DOUBLE PRECISION array, dimension LX.
*          On input, the initial guess(es) (not required, see NINIT). 
*          On output, the iterated approximate eigenvector(s). 
*              On output (input), approximate eigenvector number i is
*              (or should be) stored in X(1+N*(i-1):N*i), 
*              for i=1,...,NINIT
*
*  LX      (input) INTEGER
*          Dimension of X. Should be at least
*              N*(2*MADSPACE+NEIG+4)+3*MADSPACE**2+MAX(MADSPACE**2,NEIG)
*          If MADSPACE.GE.3, use LX not smaller than
*              N*(3*MADSPACE+NEIG+1)+3*MADSPACE**2+MAX(MADSPACE**2,NEIG)
*          to guarantee optimal performance.
*
*  NEIG    (input/output) INTEGER
*          On input, the number of eigenvalue(s) to be computed; 
*                    should be positive.
*          On output, the number of eigenvalues effectively computed
*                    with the required accuracy.
*
*  SIGMA   (input) DOUBLE PRECISION
*          If ISEARCH.LE.0: not used
*          If ISEARCH.EQ.1: estimation of the smallest eigenvalue
*                  (may speed up somewhat the algorithm if not too
*                   inaccurate)
*          If ISEARCH.GE.2: the "target", see ISEARCH
*
*  ISEARCH (input) INTEGER
*          If ISEARCH.LE.0: compute the smallest eigenvalue(s)
*          If ISEARCH.EQ.1: compute the smallest eigenvalue(s) and use
*                   SIGMA as initial guess. If one searches for the
*                   smallest eigenvalue(s) and has to rerun the
*                   algorithm for the same problem (or a problem with
*                   similar eigenvalues at the lower end), it is a good
*                   idea to set ISEARCH=1 and SIGMA=EIGS(1) (as obtained
*                   from the first run).
*          If ISEARCH.EQ.2: compute the eigenvalues closest to SIGMA
*
*  NINIT   (input) INTEGER
*          Number of initial guess(es) provided. May be set to 0.
*
*  MADSPACE (input) INTEGER
*           Maximal dimension of the search space (usually between 10 
*           and 20). Should be at least 2.
*
*  ITER    (input/output) INTEGER
*          On input, the maximum number of matrix vector multiplications; 
*                    should be positive.
*          On output, actual number of matrix vector multiplications.
*
*  TOL     (input) DOUBLE PRECISION
*          The tolerance on residual norm. Iterations to compute
*               eigenvector number i are stopped whenever
*               || A*x(i)-EIGS(i)*x(i) || .LE. TOL
*               Should be positive.
*
*  SHIFT   (input/output) DOUBLE PRECISION
*          used only if ISEARCH.EQ.1
*          On input, SHIFT is used to shift the input matrix by a
*              multiple of the identity matrix in order to compute the 
*              preconditioner. A good heuristic is obtained by setting
*              SHIFT equal to SIGMA (the smallest eigenvalue estimate)
*              minus the estimated gap between this smallest eigenvalue
*              and the next one (i.e., SHIFT approximates 2 lambda_1 -
*              lambda_2, where lambda_1 (lambda_2) is smallest (second 
*              smallest) eigenvalue). If one has no idea of this gap,
*              SHIFT may be set equal to SIGMA.
*          On output: suggested new value for the SHIFT parameter (not
*              necessarily equal to the current estimation of 2 lambda_1
*              - lambda_2). If one searches for the smallest eigenvalue(s) 
*              and has to rerun the algorithm  for the same problem (or
*              a problem with similar eigenvalues at the lower end), it 
*              is a good idea to set ISEARCH=1, SIGMA=EIGS(1) and SHIFT
*              equal to the value on output from PJD/PJDREVCOM.
*
* DROPTOL  (input/output) DOUBLE PRECISION
*          On input, drop tolerance for the multilevel incomplete 
*              factorization. A small drop tolerance will typically lead
*              to more fill-in, i.e. more memory will be consumed and 
*              the application of the preconditioner is more costly. On 
*              the other hand, the number of iteration steps is expected
*              to be less for a smaller drop tolerance.
*          On output: suggested new value for the DROPTOL parameter, 
*              that might be useful to if one has to rerun the algorithm
*              for a similar problem.
*
* MEM      (input) DOUBLE PRECISION
*          MEM prescribes the amount of memory the user is willing to 
*              spend for the preconditioner. MEM is relative to the
*              number of nonzero of the input matrix. If it turns out
*              the preconditioner that is computed does not fit into the
*              memory area that is offered by the user, PJDINIT will
*              terminate with an error message. In this case one can
*              either increase MEM (if there is more memory available)
*              or one has to increase DROPTOL.
*
* ICNTL    (input) INTEGER
*          some control parameters
*          ICNTL(1) not used (kept for compatibility with PJDREVCOM).
*          ICNTL(2) If equal to zero, then adaptive preconditioning is
*                   used, i.e., during the eigenvector computation the
*                   ILU may be recomputed (with different SHIFT and
*                   DROPTOL), if useful and possible.
*                   If not equal to zero, then the preconditioner
*                   computed initially is kept throughout. If, in
*                   addition, ICNTL(2) is equal to -1, the existing
*                   preconditioner is reused in a static way (this 
*                   option presumes that PJD was called previously and
*                   successful for the same problem). Finally, ICNTL(2)
*                   set to -2 will cause that a previously existing
*                   preconditioner will be reused in a adaptive fashion.
*          ICNTL(3) if ICNTL(2) is equal to zero and ISEARCH.NE.2,
*                   then ICNTL(3) states whether negative diagonal
*                   entries that show up in the ILU will be changed to
*                   positive ones. If set to zero (default), then up to
*                   1% of the negative diagonal entries are converted.
*                   If more negative diagonal are discovered then the
*                   algorithm searches for a new shift (and possibly a
*                   different DROPTOL if adaptive preconditioning is
*                   used).
*                   If ICNTL(3) is set to 1, then no negative diagonal
*                   entries are permitted forcing the algorithm to
*                   compute a different shift.  
*          ICNTL(4) if set to zero, default estimate for norm of the
*                   inverse factor is used. Otherwise use ICNTL(4) as
*                   bound.
*
* IPRINT   (input) INTEGER
*          Its absolute value indicates the unit number where 
*              information is to be printed (N.B.: 5 is converted to 6).
*              If zero, only error messages are printed on standard
*              output. Otherwise, its sign indicates the level of output:
*              if negative, extensive information (for expert) is
*              provided; most users should be satisfied with the 
*              information provided for positive IPRINT.
*          
* INFO     (output) INTEGER