text:
*DECK DWNNLS
      SUBROUTINE DWNNLS (W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE,
     +   IWORK, WORK)
C***BEGIN PROLOGUE  DWNNLS
C***PURPOSE  Solve a linearly constrained least squares problem with
C            equality constraints and nonnegativity constraints on
C            selected variables.
C***LIBRARY   SLATEC
C***CATEGORY  K1A2A
C***TYPE      DOUBLE PRECISION (WNNLS-S, DWNNLS-D)
C***KEYWORDS  CONSTRAINED LEAST SQUARES, CURVE FITTING, DATA FITTING,
C             EQUALITY CONSTRAINTS, INEQUALITY CONSTRAINTS,
C             NONNEGATIVITY CONSTRAINTS, QUADRATIC PROGRAMMING
C***AUTHOR  Hanson, R. J., (SNLA)
C           Haskell, K. H., (SNLA)
C***DESCRIPTION
C
C     Abstract
C
C     This subprogram solves a linearly constrained least squares
C     problem.  Suppose there are given matrices E and A of
C     respective dimensions ME by N and MA by N, and vectors F
C     and B of respective lengths ME and MA.  This subroutine
C     solves the problem
C
C               EX = F, (equations to be exactly satisfied)
C
C               AX = B, (equations to be approximately satisfied,
C                        in the least squares sense)
C
C               subject to components L+1,...,N nonnegative
C
C     Any values ME.GE.0, MA.GE.0 and 0.LE. L .LE.N are permitted.
C
C     The problem is reposed as problem DWNNLS
C
C               (WT*E)X = (WT*F)
C               (   A)    (   B), (least squares)
C               subject to components L+1,...,N nonnegative.
C
C     The subprogram chooses the heavy weight (or penalty parameter) WT.
C
C     The parameters for DWNNLS are
C
C     INPUT.. All TYPE REAL variables are DOUBLE PRECISION
C
C     W(*,*),MDW,  The array W(*,*) is double subscripted with first
C     ME,MA,N,L    dimensioning parameter equal to MDW.  For this
C                  discussion let us call M = ME + MA.  Then MDW
C                  must satisfy MDW.GE.M.  The condition MDW.LT.M
C                  is an error.
C
C                  The array W(*,*) contains the matrices and vectors
C
C                       (E  F)
C                       (A  B)
C
C                  in rows and columns 1,...,M and 1,...,N+1
C                  respectively.  Columns 1,...,L correspond to
C                  unconstrained variables X(1),...,X(L).  The
C                  remaining variables are constrained to be
C                  nonnegative. The condition L.LT.0 or L.GT.N is
C                  an error.
C
C     PRGOPT(*)    This double precision array is the option vector.
C                  If the user is satisfied with the nominal
C                  subprogram features set
C
C                  PRGOPT(1)=1 (or PRGOPT(1)=1.0)
C
C                  Otherwise PRGOPT(*) is a linked list consisting of
C                  groups of data of the following form
C
C                  LINK
C                  KEY
C                  DATA SET
C
C                  The parameters LINK and KEY are each one word.
C                  The DATA SET can be comprised of several words.
C                  The number of items depends on the value of KEY.
C                  The value of LINK points to the first
C                  entry of the next group of data within
C                  PRGOPT(*).  The exception is when there are
C                  no more options to change.  In that
C                  case LINK=1 and the values KEY and DATA SET
C                  are not referenced. The general layout of
C                  PRGOPT(*) is as follows.
C
C               ...PRGOPT(1)=LINK1 (link to first entry of next group)
C               .  PRGOPT(2)=KEY1 (key to the option change)
C               .  PRGOPT(3)=DATA VALUE (data value for this change)
C               .       .
C               .       .
C               .       .
C               ...PRGOPT(LINK1)=LINK2 (link to the first entry of
C               .                       next group)
C               .  PRGOPT(LINK1+1)=KEY2 (key to the option change)
C               .  PRGOPT(LINK1+2)=DATA VALUE
C               ...     .
C               .       .
C               .       .
C               ...PRGOPT(LINK)=1 (no more options to change)
C
C                  Values of LINK that are nonpositive are errors.
C                  A value of LINK.GT.NLINK=100000 is also an error.
C                  This helps prevent using invalid but positive
C                  values of LINK that will probably extend
C                  beyond the program limits of PRGOPT(*).
C                  Unrecognized values of KEY are ignored.  The
C                  order of the options is arbitrary and any number
C                  of options can be changed with the following
C                  restriction.  To prevent cycling in the
C                  processing of the option array a count of the
C                  number of options changed is maintained.
C                  Whenever this count exceeds NOPT=1000 an error
C                  message is printed and the subprogram returns.
C
C                  OPTIONS..
C
C                  KEY=6
C                         Scale the nonzero columns of the
C                  entire data matrix
C                  (E)
C                  (A)
C                  to have length one. The DATA SET for
C                  this option is a single value.  It must
C                  be nonzero if unit length column scaling is
C                  desired.
C
C                  KEY=7
C                         Scale columns of the entire data matrix
C                  (E)
C                  (A)
C                  with a user-provided diagonal matrix.
C                  The DATA SET for this option consists
C                  of the N diagonal scaling factors, one for
C                  each matrix column.
C
C                  KEY=8
C                         Change the rank determination tolerance from
C                  the nominal value of SQRT(SRELPR).  This quantity
C                  can be no smaller than SRELPR, The arithmetic-
C                  storage precision.  The quantity used
C                  here is internally restricted to be at
C                  least SRELPR.  The DATA SET for this option
C                  is the new tolerance.
C
C                  KEY=9
C                         Change the blow-up parameter from the
C                  nominal value of SQRT(SRELPR).  The reciprocal of
C                  this parameter is used in rejecting solution
C                  components as too large when a variable is
C                  first brought into the active set.  Too large
C                  means that the proposed component times the
C                  reciprocal of the parameter is not less than
C                  the ratio of the norms of the right-side
C                  vector and the data matrix.
C                  This parameter can be no smaller than SRELPR,
C                  the arithmetic-storage precision.
C
C                  For example, suppose we want to provide
C                  a diagonal matrix to scale the problem
C                  matrix and change the tolerance used for
C                  determining linear dependence of dropped col
C                  vectors.  For these options the dimensions of
C                  PRGOPT(*) must be at least N+6.  The FORTRAN
C                  statements defining these options would
C                  be as follows.
C
C                  PRGOPT(1)=N+3 (link to entry N+3 in PRGOPT(*))
C                  PRGOPT(2)=7 (user-provided scaling key)
C
C                  CALL DCOPY(N,D,1,PRGOPT(3),1) (copy the N
C                  scaling factors from a user array called D(*)
C                  into PRGOPT(3)-PRGOPT(N+2))
C
C                  PRGOPT(N+3)=N+6 (link to entry N+6 of PRGOPT(*))
C                  PRGOPT(N+4)=8 (linear dependence tolerance key)
C                  PRGOPT(N+5)=... (new value of the tolerance)
C
C                  PRGOPT(N+6)=1 (no more options to change)
C
C
C     IWORK(1),    The amounts of working storage actually allocated
C     IWORK(2)     for the working arrays WORK(*) and IWORK(*),
C                  respectively.  These quantities are compared with
C                  the actual amounts of storage needed for DWNNLS( ).
C                  Insufficient storage allocated for either WORK(*)
C                  or IWORK(*) is considered an error.  This feature
C                  was included in DWNNLS( ) because miscalculating
C                  the storage formulas for WORK(*) and IWORK(*)
C                  might very well lead to subtle and hard-to-find
C                  execution errors.
C
C                  The length of WORK(*) must be at least
C
C                  LW = ME+MA+5*N
C                  This test will not be made if IWORK(1).LE.0.
C
C                  The length of IWORK(*) must be at least
C
C                  LIW = ME+MA+N
C                  This test will not be made if IWORK(2).LE.0.
C
C     OUTPUT.. All TYPE REAL variables are DOUBLE PRECISION
C
C     X(*)         An array dimensioned at least N, which will
C                  contain the N components of the solution vector
C                  on output.
C
C     RNORM        The residual norm of the solution.  The value of
C                  RNORM contains the residual vector length of the
C                  equality constraints and least squares equations.
C
C     MODE         The value of MODE indicates the success or failure
C                  of the subprogram.
C
C                  MODE = 0  Subprogram completed successfully.
C
C                       = 1  Max. number of iterations (equal to
C                            3*(N-L)) exceeded. Nearly all problems
C                            should complete in fewer than this
C                            number of iterations. An approximate
C                            solution and its corresponding residual
C                            vector length are in X(*) and RNORM.
C
C                       = 2  Usage error occurred.  The offending
C                            condition is noted with the error
C                            processing subprogram, XERMSG( ).
C
C     User-designated
C     Working arrays..
C
C     WORK(*)      A double precision working array of length at least
C                  M + 5*N.
C
C     IWORK(*)     An integer-valued working array of length at least
C                  M+N.
C
C***REFERENCES  K. H. Haskell and R. J. Hanson, An algorithm for
C                 linear least squares problems with equality and
C                 nonnegativity constraints, Report SAND77-0552, Sandia
C                 Laboratories, June 1978.
C               K. H. Haskell and R. J. Hanson, Selected algorithms for
C                 the linearly constrained least squares problem - a
C                 users guide, Report SAND78-1290, Sandia Laboratories,
C                 August 1979.
C               K. H. Haskell and R. J. Hanson, An algorithm for
C                 linear least squares problems with equality and
C                 nonnegativity constraints, Mathematical Programming
C                 21 (1981), pp. 98-118.
C               R. J. Hanson and K. H. Haskell, Two algorithms for the
C                 linearly constrained least squares problem, ACM
C                 Transactions on Mathematical Software, September 1982.
C               C. L. Lawson and R. J. Hanson, Solving Least Squares
C                 Problems, Prentice-Hall, Inc., 1974.
C***ROUTINES CALLED  DWNLSM, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   790701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890618  Completely restructured and revised.  (WRB & RWC)
C   891006  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900510  Convert XERRWV calls to XERMSG calls, change Prologue
C           comments to agree with WNNLS.  (RWC)
C   920501  Reformatted the REFERENCES section.  (WRB)
C***END PROLOGUE  DWNNLS
      INTEGER IWORK(*), L, L1, L2, L3, L4, L5, LIW, LW, MA, MDW, ME,
     *     MODE, N
      DOUBLE PRECISION  PRGOPT(*), RNORM, W(MDW,*), WORK(*), X(*)
      CHARACTER*8 XERN1
C***FIRST EXECUTABLE STATEMENT  DWNNLS
      MODE = 0
      IF (MA+ME.LE.0 .OR. N.LE.0) RETURN
C
      IF (IWORK(1).GT.0) THEN
         LW = ME + MA + 5*N
         IF (IWORK(1).LT.LW) THEN
            WRITE (XERN1, '(I8)') LW
            CALL XERMSG ('SLATEC', 'DWNNLS', 'INSUFFICIENT STORAGE ' //
     *         'ALLOCATED FOR WORK(*), NEED LW = ' // XERN1, 2, 1)
            MODE = 2
            RETURN
         ENDIF
      ENDIF
C
      IF (IWORK(2).GT.0) THEN
         LIW = ME + MA + N
         IF (IWORK(2).LT.LIW) THEN
            WRITE (XERN1, '(I8)') LIW
            CALL XERMSG ('SLATEC', 'DWNNLS', 'INSUFFICIENT STORAGE ' //
     *         'ALLOCATED FOR IWORK(*), NEED LIW = ' // XERN1, 2, 1)
            MODE = 2
            RETURN
         ENDIF
      ENDIF
C
      IF (MDW.LT.ME+MA) THEN
         CALL XERMSG ('SLATEC', 'DWNNLS',
     *      'THE VALUE MDW.LT.ME+MA IS AN ERROR', 1, 1)
         MODE = 2
         RETURN
      ENDIF
C
      IF (L.LT.0 .OR. L.GT.N) THEN
         CALL XERMSG ('SLATEC', 'DWNNLS',
     *      'L.GE.0 .AND. L.LE.N IS REQUIRED', 2, 1)
         MODE = 2
         RETURN
      ENDIF
C
C     THE PURPOSE OF THIS SUBROUTINE IS TO BREAK UP THE ARRAYS
C     WORK(*) AND IWORK(*) INTO SEPARATE WORK ARRAYS
C     REQUIRED BY THE MAIN SUBROUTINE DWNLSM( ).
C
      L1 = N + 1
      L2 = L1 + N
      L3 = L2 + ME + MA
      L4 = L3 + N
      L5 = L4 + N
C
      CALL DWNLSM(W, MDW, ME, MA, N, L, PRGOPT, X, RNORM, MODE, IWORK,
     *            IWORK(L1), WORK(1), WORK(L1), WORK(L2), WORK(L3),
     *            WORK(L4), WORK(L5))
      RETURN
      END
*DECK DWNLSM
      SUBROUTINE DWNLSM (W, MDW, MME, MA, N, L, PRGOPT, X, RNORM, MODE,
     +   IPIVOT, ITYPE, WD, H, SCALE, Z, TEMP, D)
C***BEGIN PROLOGUE  DWNLSM
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DWNNLS
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (WNLSM-S, DWNLSM-D)
C***AUTHOR  Hanson, R. J., (SNLA)
C           Haskell, K. H., (SNLA)
C***DESCRIPTION
C
C     This is a companion subprogram to DWNNLS.
C     The documentation for DWNNLS has complete usage instructions.
C
C     In addition to the parameters discussed in the prologue to
C     subroutine DWNNLS, the following work arrays are used in
C     subroutine DWNLSM  (they are passed through the calling
C     sequence from DWNNLS for purposes of variable dimensioning).
C     Their contents will in general be of no interest to the user.
C
C     Variables of type REAL are DOUBLE PRECISION.
C
C         IPIVOT(*)
C            An array of length N.  Upon completion it contains the
C         pivoting information for the cols of W(*,*).
C
C         ITYPE(*)
C            An array of length M which is used to keep track
C         of the classification of the equations.  ITYPE(I)=0
C         denotes equation I as an equality constraint.
C         ITYPE(I)=1 denotes equation I as a least squares
C         equation.
C
C         WD(*)
C            An array of length N.  Upon completion it contains the
C         dual solution vector.
C
C         H(*)
C            An array of length N.  Upon completion it contains the
C         pivot scalars of the Householder transformations performed
C         in the case KRANK.LT.L.
C
C         SCALE(*)
C            An array of length M which is used by the subroutine
C         to store the diagonal matrix of weights.
C         These are used to apply the modified Givens
C         transformations.
C
C         Z(*),TEMP(*)
C            Working arrays of length N.
C
C         D(*)
C            An array of length N that contains the
C         column scaling for the matrix (E).
C                                       (A)
C
C***SEE ALSO  DWNNLS
C***ROUTINES CALLED  D1MACH, DASUM, DAXPY, DCOPY, DH12, DNRM2, DROTM,
C                    DROTMG, DSCAL, DSWAP, DWNLIT, IDAMAX, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   790701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890618  Completely restructured and revised.  (WRB & RWC)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900328  Added TYPE section.  (WRB)
C   900510  Fixed an error message.  (RWC)
C   900604  DP version created from SP version.  (RWC)
C   900911  Restriction on value of ALAMDA included.  (WRB)
C***END PROLOGUE  DWNLSM
      INTEGER IPIVOT(*), ITYPE(*), L, MA, MDW, MME, MODE, N
      DOUBLE PRECISION D(*), H(*), PRGOPT(*), RNORM, SCALE(*), TEMP(*),
     *   W(MDW,*), WD(*), X(*), Z(*)
C
      EXTERNAL D1MACH, DASUM, DAXPY, DCOPY, DH12, DNRM2, DROTM, DROTMG,
     *   DSCAL, DSWAP, DWNLIT, IDAMAX, XERMSG
      DOUBLE PRECISION D1MACH, DASUM, DNRM2
      INTEGER IDAMAX
C
      DOUBLE PRECISION ALAMDA, ALPHA, ALSQ, AMAX, BLOWUP, BNORM,
     *   DOPE(3), DRELPR, EANORM, FAC, SM, SPARAM(5), T, TAU, WMAX, Z2,
     *   ZZ
      INTEGER I, IDOPE(3), IMAX, ISOL, ITEMP, ITER, ITMAX, IWMAX, J,
     *   JCON, JP, KEY, KRANK, L1, LAST, LINK, M, ME, NEXT, NIV, NLINK,
     *   NOPT, NSOLN, NTIMES
      LOGICAL DONE, FEASBL, FIRST, HITCON, POS
C
      SAVE DRELPR, FIRST
      DATA FIRST /.TRUE./
C***FIRST EXECUTABLE STATEMENT  DWNLSM
C
C     Initialize variables.
C     DRELPR is the precision for the particular machine
C     being used.  This logic avoids resetting it every entry.
C
      IF (FIRST) DRELPR = D1MACH(4)
      FIRST = .FALSE.
C
C     Set the nominal tolerance used in the code.
C
      TAU = SQRT(DRELPR)
C
      M = MA + MME
      ME = MME
      MODE = 2
C
C     To process option vector
C
      FAC = 1.D-4
C
C     Set the nominal blow up factor used in the code.
C
      BLOWUP = TAU
C
C     The nominal column scaling used in the code is
C     the identity scaling.
C
      CALL DCOPY (N, 1.D0, 0, D, 1)
C
C     Define bound for number of options to change.
C
      NOPT = 1000
C
C     Define bound for positive value of LINK.
C
      NLINK = 100000
      NTIMES = 0
      LAST = 1
      LINK = PRGOPT(1)
      IF (LINK.LE.0 .OR. LINK.GT.NLINK) THEN
         CALL XERMSG ('SLATEC', 'DWNLSM',
     +      'IN DWNNLS, THE OPTION VECTOR IS UNDEFINED', 3, 1)
         RETURN
      ENDIF
C
  100 IF (LINK.GT.1) THEN
         NTIMES = NTIMES + 1
         IF (NTIMES.GT.NOPT) THEN
         CALL XERMSG ('SLATEC', 'DWNLSM',
     +      'IN DWNNLS, THE LINKS IN THE OPTION VECTOR ARE CYCLING.',
     +      3, 1)
            RETURN
         ENDIF
C
         KEY = PRGOPT(LAST+1)
         IF (KEY.EQ.6 .AND. PRGOPT(LAST+2).NE.0.D0) THEN
            DO 110 J = 1,N
               T = DNRM2(M,W(1,J),1)
               IF (T.NE.0.D0) T = 1.D0/T
               D(J) = T
  110       CONTINUE
         ENDIF
C
         IF (KEY.EQ.7) CALL DCOPY (N, PRGOPT(LAST+2), 1, D, 1)
         IF (KEY.EQ.8) TAU = MAX(DRELPR,PRGOPT(LAST+2))
         IF (KEY.EQ.9) BLOWUP = MAX(DRELPR,PRGOPT(LAST+2))
C
         NEXT = PRGOPT(LINK)
         IF (NEXT.LE.0 .OR. NEXT.GT.NLINK) THEN
            CALL XERMSG ('SLATEC', 'DWNLSM',
     +         'IN DWNNLS, THE OPTION VECTOR IS UNDEFINED', 3, 1)
            RETURN
         ENDIF
C
         LAST = LINK
         LINK = NEXT
         GO TO 100
      ENDIF
C
      DO 120 J = 1,N
         CALL DSCAL (M, D(J), W(1,J), 1)
  120 CONTINUE
C
C     Process option vector
C
      DONE = .FALSE.
      ITER = 0
      ITMAX = 3*(N-L)
      MODE = 0
      NSOLN = L
      L1 = MIN(M,L)
C
C     Compute scale factor to apply to equality constraint equations.
C
      DO 130 J = 1,N
         WD(J) = DASUM(M,W(1,J),1)
  130 CONTINUE
C
      IMAX = IDAMAX(N,WD,1)
      EANORM = WD(IMAX)
      BNORM = DASUM(M,W(1,N+1),1)
      ALAMDA = EANORM/(DRELPR*FAC)
C
C     On machines, such as the VAXes using D floating, with a very
C     limited exponent range for double precision values, the previously
C     computed value of ALAMDA may cause an overflow condition.
C     Therefore, this code further limits the value of ALAMDA.
C
      ALAMDA = MIN(ALAMDA,SQRT(D1MACH(2)))
C
C     Define scaling diagonal matrix for modified Givens usage and
C     classify equation types.
C
      ALSQ = ALAMDA**2
      DO 140 I = 1,M
C
C        When equation I is heavily weighted ITYPE(I)=0,
C        else ITYPE(I)=1.
C
         IF (I.LE.ME) THEN
            T = ALSQ
            ITEMP = 0
         ELSE
            T = 1.D0
            ITEMP = 1
         ENDIF
         SCALE(I) = T
         ITYPE(I) = ITEMP
  140 CONTINUE
C
C     Set the solution vector X(*) to zero and the column interchange
C     matrix to the identity.
C
      CALL DCOPY (N, 0.D0, 0, X, 1)
      DO 150 I = 1,N
         IPIVOT(I) = I
  150 CONTINUE
C
C     Perform initial triangularization in the submatrix
C     corresponding to the unconstrained variables.
C     Set first L components of dual vector to zero because
C     these correspond to the unconstrained variables.
C
      CALL DCOPY (L, 0.D0, 0, WD, 1)
C
C     The arrays IDOPE(*) and DOPE(*) are used to pass
C     information to DWNLIT().  This was done to avoid
C     a long calling sequence or the use of COMMON.
C
      IDOPE(1) = ME
      IDOPE(2) = NSOLN
      IDOPE(3) = L1
C
      DOPE(1) = ALSQ
      DOPE(2) = EANORM
      DOPE(3) = TAU
      CALL DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE, RNORM,
     +            IDOPE, DOPE, DONE)
      ME    = IDOPE(1)
      KRANK = IDOPE(2)
      NIV   = IDOPE(3)
C
C     Perform WNNLS algorithm using the following steps.
C
C     Until(DONE)
C        compute search direction and feasible point
C        when (HITCON) add constraints
C        else perform multiplier test and drop a constraint
C        fin
C     Compute-Final-Solution
C
C     To compute search direction and feasible point,
C     solve the triangular system of currently non-active
C     variables and store the solution in Z(*).
C
C     To solve system
C     Copy right hand side into TEMP vector to use overwriting method.
C
  160 IF (DONE) GO TO 330
      ISOL = L + 1
      IF (NSOLN.GE.ISOL) THEN
         CALL DCOPY (NIV, W(1,N+1), 1, TEMP, 1)
         DO 170 J = NSOLN,ISOL,-1
            IF (J.GT.KRANK) THEN
               I = NIV - NSOLN + J
            ELSE
               I = J
            ENDIF
C
            IF (J.GT.KRANK .AND. J.LE.L) THEN
               Z(J) = 0.D0
            ELSE
               Z(J) = TEMP(I)/W(I,J)
               CALL DAXPY (I-1, -Z(J), W(1,J), 1, TEMP, 1)
            ENDIF
  170    CONTINUE
      ENDIF
C
C     Increment iteration counter and check against maximum number
C     of iterations.
C
      ITER = ITER + 1
      IF (ITER.GT.ITMAX) THEN
         MODE = 1
         DONE = .TRUE.
      ENDIF
C
C     Check to see if any constraints have become active.
C     If so, calculate an interpolation factor so that all
C     active constraints are removed from the basis.
C
      ALPHA = 2.D0
      HITCON = .FALSE.
      DO 180 J = L+1,NSOLN
         ZZ = Z(J)
         IF (ZZ.LE.0.D0) THEN
            T = X(J)/(X(J)-ZZ)
            IF (T.LT.ALPHA) THEN
               ALPHA = T
               JCON = J
            ENDIF
            HITCON = .TRUE.
         ENDIF
  180 CONTINUE
C
C     Compute search direction and feasible point
C
      IF (HITCON) THEN
C
C        To add constraints, use computed ALPHA to interpolate between
C        last feasible solution X(*) and current unconstrained (and
C        infeasible) solution Z(*).
C
         DO 190 J = L+1,NSOLN
            X(J) = X(J) + ALPHA*(Z(J)-X(J))
  190    CONTINUE
         FEASBL = .FALSE.
C
C        Remove column JCON and shift columns JCON+1 through N to the
C        left.  Swap column JCON into the N th position.  This achieves
C        upper Hessenberg form for the nonactive constraints and
C        leaves an upper Hessenberg matrix to retriangularize.
C
  200    DO 210 I = 1,M
            T = W(I,JCON)
            CALL DCOPY (N-JCON, W(I, JCON+1), MDW, W(I, JCON), MDW)
            W(I,N) = T
  210    CONTINUE
C
C        Update permuted index vector to reflect this shift and swap.
C
         ITEMP = IPIVOT(JCON)
         DO 220 I = JCON,N - 1
            IPIVOT(I) = IPIVOT(I+1)
  220    CONTINUE
         IPIVOT(N) = ITEMP
C
C        Similarly permute X(*) vector.
C
         CALL DCOPY (N-JCON, X(JCON+1), 1, X(JCON), 1)
         X(N) = 0.D0
         NSOLN = NSOLN - 1
         NIV = NIV - 1
C
C        Retriangularize upper Hessenberg matrix after adding
C        constraints.
C
         I = KRANK + JCON - L
         DO 230 J = JCON,NSOLN
            IF (ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.0) THEN
C
C              Zero IP1 to I in column J
C
               IF (W(I+1,J).NE.0.D0) THEN
                  CALL DROTMG (SCALE(I), SCALE(I+1), W(I,J), W(I+1,J),
     +                         SPARAM)
                  W(I+1,J) = 0.D0
                  CALL DROTM (N+1-J, W(I,J+1), MDW, W(I+1,J+1), MDW,
     +                        SPARAM)
               ENDIF
            ELSEIF (ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.1) THEN
C
C              Zero IP1 to I in column J
C
               IF (W(I+1,J).NE.0.D0) THEN
                  CALL DROTMG (SCALE(I), SCALE(I+1), W(I,J), W(I+1,J),
     +                         SPARAM)
                  W(I+1,J) = 0.D0
                  CALL DROTM (N+1-J, W(I,J+1), MDW, W(I+1,J+1), MDW,
     +                        SPARAM)
               ENDIF
            ELSEIF (ITYPE(I).EQ.1 .AND. ITYPE(I+1).EQ.0) THEN
               CALL DSWAP (N+1, W(I,1), MDW, W(I+1,1), MDW)
               CALL DSWAP (1, SCALE(I), 1, SCALE(I+1), 1)
               ITEMP = ITYPE(I+1)
               ITYPE(I+1) = ITYPE(I)
               ITYPE(I) = ITEMP
C
C              Swapped row was formerly a pivot element, so it will
C              be large enough to perform elimination.
C              Zero IP1 to I in column J.
C
               IF (W(I+1,J).NE.0.D0) THEN
                  CALL DROTMG (SCALE(I), SCALE(I+1), W(I,J), W(I+1,J),
     +                         SPARAM)
                  W(I+1,J) = 0.D0
                  CALL DROTM (N+1-J, W(I,J+1), MDW, W(I+1,J+1), MDW,
     +                        SPARAM)
               ENDIF
            ELSEIF (ITYPE(I).EQ.0 .AND. ITYPE(I+1).EQ.1) THEN
               IF (SCALE(I)*W(I,J)**2/ALSQ.GT.(TAU*EANORM)**2) THEN
C
C                 Zero IP1 to I in column J
C
                  IF (W(I+1,J).NE.0.D0) THEN
                     CALL DROTMG (SCALE(I), SCALE(I+1), W(I,J),
     +                            W(I+1,J), SPARAM)
                     W(I+1,J) = 0.D0
                     CALL DROTM (N+1-J, W(I,J+1), MDW, W(I+1,J+1), MDW,
     +                           SPARAM)
                  ENDIF
               ELSE
                  CALL DSWAP (N+1, W(I,1), MDW, W(I+1,1), MDW)
                  CALL DSWAP (1, SCALE(I), 1, SCALE(I+1), 1)
                  ITEMP = ITYPE(I+1)
                  ITYPE(I+1) = ITYPE(I)
                  ITYPE(I) = ITEMP
                  W(I+1,J) = 0.D0
               ENDIF
            ENDIF
            I = I + 1
  230    CONTINUE
C
C        See if the remaining coefficients in the solution set are
C        feasible.  They should be because of the way ALPHA was
C        determined.  If any are infeasible, it is due to roundoff
C        error.  Any that are non-positive will be set to zero and
C        removed from the solution set.
C
         DO 240 JCON = L+1,NSOLN
            IF (X(JCON).LE.0.D0) GO TO 250
  240    CONTINUE
         FEASBL = .TRUE.
  250    IF (.NOT.FEASBL) GO TO 200
      ELSE
C
C        To perform multiplier test and drop a constraint.
C
         CALL DCOPY (NSOLN, Z, 1, X, 1)
         IF (NSOLN.LT.N) CALL DCOPY (N-NSOLN, 0.D0, 0, X(NSOLN+1), 1)
C
C        Reclassify least squares equations as equalities as necessary.
C
         I = NIV + 1
  260    IF (I.LE.ME) THEN
            IF (ITYPE(I).EQ.0) THEN
               I = I + 1
            ELSE
               CALL DSWAP (N+1, W(I,1), MDW, W(ME,1), MDW)
               CALL DSWAP (1, SCALE(I), 1, SCALE(ME), 1)
               ITEMP = ITYPE(I)
               ITYPE(I) = ITYPE(ME)
               ITYPE(ME) = ITEMP
               ME = ME - 1
            ENDIF
            GO TO 260
         ENDIF
C
C        Form inner product vector WD(*) of dual coefficients.
C
         DO 280 J = NSOLN+1,N
            SM = 0.D0
            DO 270 I = NSOLN+1,M
               SM = SM + SCALE(I)*W(I,J)*W(I,N+1)
  270       CONTINUE
            WD(J) = SM
  280    CONTINUE
C
C        Find J such that WD(J)=WMAX is maximum.  This determines
C        that the incoming column J will reduce the residual vector
C        and be positive.
C
  290    WMAX = 0.D0
         IWMAX = NSOLN + 1
         DO 300 J = NSOLN+1,N
            IF (WD(J).GT.WMAX) THEN
               WMAX = WD(J)
               IWMAX = J
            ENDIF
  300    CONTINUE
         IF (WMAX.LE.0.D0) GO TO 330
C
C        Set dual coefficients to zero for incoming column.
C
         WD(IWMAX) = 0.D0
C
C        WMAX .GT. 0.D0, so okay to move column IWMAX to solution set.
C        Perform transformation to retriangularize, and test for near
C        linear dependence.
C
C        Swap column IWMAX into NSOLN-th position to maintain upper
C        Hessenberg form of adjacent columns, and add new column to
C        triangular decomposition.
C
         NSOLN = NSOLN + 1
         NIV = NIV + 1
         IF (NSOLN.NE.IWMAX) THEN
            CALL DSWAP (M, W(1,NSOLN), 1, W(1,IWMAX), 1)
            WD(IWMAX) = WD(NSOLN)
            WD(NSOLN) = 0.D0
            ITEMP = IPIVOT(NSOLN)
            IPIVOT(NSOLN) = IPIVOT(IWMAX)
            IPIVOT(IWMAX) = ITEMP
         ENDIF
C
C        Reduce column NSOLN so that the matrix of nonactive constraints
C        variables is triangular.
C
         DO 320 J = M,NIV+1,-1
            JP = J - 1
C
C           When operating near the ME line, test to see if the pivot
C           element is near zero.  If so, use the largest element above
C           it as the pivot.  This is to maintain the sharp interface
C           between weighted and non-weighted rows in all cases.
C
            IF (J.EQ.ME+1) THEN
               IMAX = ME
               AMAX = SCALE(ME)*W(ME,NSOLN)**2
               DO 310 JP = J - 1,NIV,-1
                  T = SCALE(JP)*W(JP,NSOLN)**2
                  IF (T.GT.AMAX) THEN
                     IMAX = JP
                     AMAX = T
                  ENDIF
  310          CONTINUE
               JP = IMAX
            ENDIF
C
            IF (W(J,NSOLN).NE.0.D0) THEN
               CALL DROTMG (SCALE(JP), SCALE(J), W(JP,NSOLN),
     +                      W(J,NSOLN), SPARAM)
               W(J,NSOLN) = 0.D0
               CALL DROTM (N+1-NSOLN, W(JP,NSOLN+1), MDW, W(J,NSOLN+1),
     +                     MDW, SPARAM)
            ENDIF
  320    CONTINUE
C
C        Solve for Z(NSOLN)=proposed new value for X(NSOLN).  Test if
C        this is nonpositive or too large.  If this was true or if the
C        pivot term was zero, reject the column as dependent.
C
         IF (W(NIV,NSOLN).NE.0.D0) THEN
            ISOL = NIV
            Z2 = W(ISOL,N+1)/W(ISOL,NSOLN)
            Z(NSOLN) = Z2
            POS = Z2 .GT. 0.D0
            IF (Z2*EANORM.GE.BNORM .AND. POS) THEN
               POS = .NOT. (BLOWUP*Z2*EANORM.GE.BNORM)
            ENDIF
C
C           Try to add row ME+1 as an additional equality constraint.
C           Check size of proposed new solution component.
C           Reject it if it is too large.
C
         ELSEIF (NIV.LE.ME .AND. W(ME+1,NSOLN).NE.0.D0) THEN
            ISOL = ME + 1
            IF (POS) THEN
C
C              Swap rows ME+1 and NIV, and scale factors for these rows.
C
               CALL DSWAP (N+1, W(ME+1,1), MDW, W(NIV,1), MDW)
               CALL DSWAP (1, SCALE(ME+1), 1, SCALE(NIV), 1)
               ITEMP = ITYPE(ME+1)
               ITYPE(ME+1) = ITYPE(NIV)
               ITYPE(NIV) = ITEMP
               ME = ME + 1
            ENDIF
         ELSE
            POS = .FALSE.
         ENDIF
C
         IF (.NOT.POS) THEN
            NSOLN = NSOLN - 1
            NIV = NIV - 1
         ENDIF
         IF (.NOT.(POS.OR.DONE)) GO TO 290
      ENDIF
      GO TO 160
C
C     Else perform multiplier test and drop a constraint.  To compute
C     final solution.  Solve system, store results in X(*).
C
C     Copy right hand side into TEMP vector to use overwriting method.
C
  330 ISOL = 1
      IF (NSOLN.GE.ISOL) THEN
         CALL DCOPY (NIV, W(1,N+1), 1, TEMP, 1)
         DO 340 J = NSOLN,ISOL,-1
            IF (J.GT.KRANK) THEN
               I = NIV - NSOLN + J
            ELSE
               I = J
            ENDIF
C
            IF (J.GT.KRANK .AND. J.LE.L) THEN
               Z(J) = 0.D0
            ELSE
               Z(J) = TEMP(I)/W(I,J)
               CALL DAXPY (I-1, -Z(J), W(1,J), 1, TEMP, 1)
            ENDIF
  340    CONTINUE
      ENDIF
C
C     Solve system.
C
      CALL DCOPY (NSOLN, Z, 1, X, 1)
C
C     Apply Householder transformations to X(*) if KRANK.LT.L
C
      IF (KRANK.LT.L) THEN
         DO 350 I = 1,KRANK
            CALL DH12 (2, I, KRANK+1, L, W(I,1), MDW, H(I), X, 1, 1, 1)
  350    CONTINUE
      ENDIF
C
C     Fill in trailing zeroes for constrained variables not in solution.
C
      IF (NSOLN.LT.N) CALL DCOPY (N-NSOLN, 0.D0, 0, X(NSOLN+1), 1)
C
C     Permute solution vector to natural order.
C
      DO 380 I = 1,N
         J = I
  360    IF (IPIVOT(J).EQ.I) GO TO 370
         J = J + 1
         GO TO 360
C
  370    IPIVOT(J) = IPIVOT(I)
         IPIVOT(I) = J
         CALL DSWAP (1, X(J), 1, X(I), 1)
  380 CONTINUE
C
C     Rescale the solution using the column scaling.
C
      DO 390 J = 1,N
         X(J) = X(J)*D(J)
  390 CONTINUE
C
      DO 400 I = NSOLN+1,M
         T = W(I,N+1)
         IF (I.LE.ME) T = T/ALAMDA
         T = (SCALE(I)*T)*T
         RNORM = RNORM + T
  400 CONTINUE
C
      RNORM = SQRT(RNORM)
      RETURN
      END
*DECK DWNLIT
      SUBROUTINE DWNLIT (W, MDW, M, N, L, IPIVOT, ITYPE, H, SCALE,
     +   RNORM, IDOPE, DOPE, DONE)
C***BEGIN PROLOGUE  DWNLIT
C***SUBSIDIARY
C***PURPOSE  Subsidiary to DWNNLS
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (WNLIT-S, DWNLIT-D)
C***AUTHOR  Hanson, R. J., (SNLA)
C           Haskell, K. H., (SNLA)
C***DESCRIPTION
C
C     This is a companion subprogram to DWNNLS( ).
C     The documentation for DWNNLS( ) has complete usage instructions.
C
C     Note  The M by (N+1) matrix W( , ) contains the rt. hand side
C           B as the (N+1)st col.
C
C     Triangularize L1 by L1 subsystem, where L1=MIN(M,L), with
C     col interchanges.
C
C***SEE ALSO  DWNNLS
C***ROUTINES CALLED  DCOPY, DH12, DROTM, DROTMG, DSCAL, DSWAP, DWNLT1,
C                    DWNLT2, DWNLT3, IDAMAX
C***REVISION HISTORY  (YYMMDD)
C   790701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890618  Completely restructured and revised.  (WRB & RWC)
C   890620  Revised to make WNLT1, WNLT2, and WNLT3 subroutines.  (RWC)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C   900604  DP version created from SP version. .  (RWC)
C***END PROLOGUE  DWNLIT
      INTEGER IDOPE(*), IPIVOT(*), ITYPE(*), L, M, MDW, N
      DOUBLE PRECISION DOPE(*), H(*), RNORM, SCALE(*), W(MDW,*)
      LOGICAL DONE
C
      EXTERNAL DCOPY, DH12, DROTM, DROTMG, DSCAL, DSWAP, DWNLT1,
     *   DWNLT2, DWNLT3, IDAMAX
      INTEGER IDAMAX
      LOGICAL DWNLT2
C
      DOUBLE PRECISION ALSQ, AMAX, EANORM, FACTOR, HBAR, RN, SPARAM(5),
     *   T, TAU
      INTEGER I, I1, IMAX, IR, J, J1, JJ, JP, KRANK, L1, LB, LEND, ME,
     *   MEND, NIV, NSOLN
      LOGICAL INDEP, RECALC
C
C***FIRST EXECUTABLE STATEMENT  DWNLIT
      ME    = IDOPE(1)
      NSOLN = IDOPE(2)
      L1    = IDOPE(3)
C
      ALSQ   = DOPE(1)
      EANORM = DOPE(2)
      TAU    = DOPE(3)
C
      LB     = MIN(M-1,L)
      RECALC = .TRUE.
      RNORM  = 0.D0
      KRANK  = 0
C
C     We set FACTOR=1.0 so that the heavy weight ALAMDA will be
C     included in the test for column independence.
C
      FACTOR = 1.D0
      LEND = L
      DO 180 I=1,LB
C
C        Set IR to point to the I-th row.
C
         IR = I
         MEND = M
         CALL DWNLT1 (I, LEND, M, IR, MDW, RECALC, IMAX, HBAR, H, SCALE,
     +                W)
C
C        Update column SS and find pivot column.
C
         CALL DWNLT3 (I, IMAX, M, MDW, IPIVOT, H, W)
C
C        Perform column interchange.
C        Test independence of incoming column.
C
  130    IF (DWNLT2(ME, MEND, IR, FACTOR, TAU, SCALE, W(1,I))) THEN
C
C           Eliminate I-th column below diagonal using modified Givens
C           transformations applied to (A B).
C
C           When operating near the ME line, use the largest element
C           above it as the pivot.
C
            DO 160 J=M,I+1,-1
               JP = J-1
               IF (J.EQ.ME+1) THEN
                  IMAX = ME
                  AMAX = SCALE(ME)*W(ME,I)**2
                  DO 150 JP=J-1,I,-1
                     T = SCALE(JP)*W(JP,I)**2
                     IF (T.GT.AMAX) THEN
                        IMAX = JP
                        AMAX = T
                     ENDIF
  150             CONTINUE
                  JP = IMAX
               ENDIF
C
               IF (W(J,I).NE.0.D0) THEN
                  CALL DROTMG (SCALE(JP), SCALE(J), W(JP,I), W(J,I),
     +                         SPARAM)
                  W(J,I) = 0.D0
                  CALL DROTM (N+1-I, W(JP,I+1), MDW, W(J,I+1), MDW,
     +                        SPARAM)
               ENDIF
  160       CONTINUE
         ELSE IF (LEND.GT.I) THEN
C
C           Column I is dependent.  Swap with column LEND.
C           Perform column interchange,
C           and find column in remaining set with largest SS.
C
            CALL DWNLT3 (I, LEND, M, MDW, IPIVOT, H, W)
            LEND = LEND - 1
            IMAX = IDAMAX(LEND-I+1, H(I), 1) + I - 1
            HBAR = H(IMAX)
            GO TO 130
         ELSE
            KRANK = I - 1
            GO TO 190
         ENDIF
  180 CONTINUE
      KRANK = L1
C
  190 IF (KRANK.LT.ME) THEN
         FACTOR = ALSQ
         DO 200 I=KRANK+1,ME
            CALL DCOPY (L, 0.D0, 0, W(I,1), MDW)
  200    CONTINUE
C
C        Determine the rank of the remaining equality constraint
C        equations by eliminating within the block of constrained
C        variables.  Remove any redundant constraints.
C
         RECALC = .TRUE.
         LB = MIN(L+ME-KRANK, N)
         DO 270 I=L+1,LB
            IR = KRANK + I - L
            LEND = N
            MEND = ME
            CALL DWNLT1 (I, LEND, ME, IR, MDW, RECALC, IMAX, HBAR, H,
     +                   SCALE, W)
C
C           Update col ss and find pivot col
C
            CALL DWNLT3 (I, IMAX, M, MDW, IPIVOT, H, W)
C
C           Perform column interchange
C           Eliminate elements in the I-th col.
C
            DO 240 J=ME,IR+1,-1
               IF (W(J,I).NE.0.D0) THEN
                 CALL DROTMG (SCALE(J-1), SCALE(J), W(J-1,I), W(J,I),
     +                        SPARAM)
                  W(J,I) = 0.D0
                  CALL DROTM (N+1-I, W(J-1,I+1), MDW,W(J,I+1), MDW,
     +                        SPARAM)
               ENDIF
  240       CONTINUE
C
C           I=column being eliminated.
C           Test independence of incoming column.
C           Remove any redundant or dependent equality constraints.
C
            IF (.NOT.DWNLT2(ME, MEND, IR, FACTOR,TAU,SCALE,W(1,I))) THEN
               JJ = IR
               DO 260 IR=JJ,ME
                  CALL DCOPY (N, 0.D0, 0, W(IR,1), MDW)
                  RNORM = RNORM + (SCALE(IR)*W(IR,N+1)/ALSQ)*W(IR,N+1)
                  W(IR,N+1) = 0.D0
                  SCALE(IR) = 1.D0
C
C                 Reclassify the zeroed row as a least squares equation.
C
                  ITYPE(IR) = 1
  260          CONTINUE
C
C              Reduce ME to reflect any discovered dependent equality
C              constraints.
C
               ME = JJ - 1
               GO TO 280
            ENDIF
  270    CONTINUE
      ENDIF
C
C     Try to determine the variables KRANK+1 through L1 from the
C     least squares equations.  Continue the triangularization with
C     pivot element W(ME+1,I).
C
  280 IF (KRANK.LT.L1) THEN
         RECALC = .TRUE.
C
C        Set FACTOR=ALSQ to remove effect of heavy weight from
C        test for column independence.
C
         FACTOR = ALSQ
         DO 350 I=KRANK+1,L1
C
C           Set IR to point to the ME+1-st row.
C
            IR = ME+1
            LEND = L
            MEND = M
            CALL DWNLT1 (I, L, M, IR, MDW, RECALC, IMAX, HBAR, H, SCALE,
     +                   W)
C
C           Update column SS and find pivot column.
C
            CALL DWNLT3 (I, IMAX, M, MDW, IPIVOT, H, W)
C
C           Perform column interchange.
C           Eliminate I-th column below the IR-th element.
C
            DO 320 J=M,IR+1,-1
               IF (W(J,I).NE.0.D0) THEN
                 CALL DROTMG (SCALE(J-1), SCALE(J), W(J-1,I), W(J,I),
     +                        SPARAM)
                  W(J,I) = 0.D0
                  CALL DROTM (N+1-I, W(J-1,I+1),  MDW, W(J,I+1), MDW,
     +                        SPARAM)
               ENDIF
  320       CONTINUE
C
C           Test if new pivot element is near zero.
C           If so, the column is dependent.
C           Then check row norm test to be classified as independent.
C
            T = SCALE(IR)*W(IR,I)**2
            INDEP = T .GT. (TAU*EANORM)**2
            IF (INDEP) THEN
               RN = 0.D0
               DO 340 I1=IR,M
                  DO 330 J1=I+1,N
                     RN = MAX(RN, SCALE(I1)*W(I1,J1)**2)
  330             CONTINUE
  340          CONTINUE
               INDEP = T .GT. RN*TAU**2
            ENDIF
C
C           If independent, swap the IR-th and KRANK+1-th rows to
C           maintain the triangular form.  Update the rank indicator
C           KRANK and the equality constraint pointer ME.
C
            IF (.NOT.INDEP) GO TO 360
            CALL DSWAP(N+1, W(KRANK+1,1), MDW, W(IR,1), MDW)
            CALL DSWAP(1, SCALE(KRANK+1), 1, SCALE(IR), 1)
C
C           Reclassify the least square equation as an equality
C           constraint and rescale it.
C
            ITYPE(IR) = 0
            T = SQRT(SCALE(KRANK+1))
            CALL DSCAL(N+1, T, W(KRANK+1,1), MDW)
            SCALE(KRANK+1) = ALSQ
            ME = ME+1
            KRANK = KRANK+1
  350    CONTINUE
      ENDIF
C
C     If pseudorank is less than L, apply Householder transformation.
C     from right.
C
  360 IF (KRANK.LT.L) THEN
         DO 370 J=KRANK,1,-1
            CALL DH12 (1, J, KRANK+1, L, W(J,1), MDW, H(J), W, MDW, 1,
     +                J-1)
  370    CONTINUE
      ENDIF
C
      NIV = KRANK + NSOLN - L
      IF (L.EQ.N) DONE = .TRUE.
C
C     End of initial triangularization.
C
      IDOPE(1) = ME
      IDOPE(2) = KRANK
      IDOPE(3) = NIV
      RETURN
      END
*DECK DWNLT1
      SUBROUTINE DWNLT1 (I, LEND, MEND, IR, MDW, RECALC, IMAX, HBAR, H,
     +   SCALE, W)
C***BEGIN PROLOGUE  DWNLT1
C***SUBSIDIARY
C***PURPOSE  Subsidiary to WNLIT
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (WNLT1-S, DWNLT1-D)
C***AUTHOR  Hanson, R. J., (SNLA)
C           Haskell, K. H., (SNLA)
C***DESCRIPTION
C
C     To update the column Sum Of Squares and find the pivot column.
C     The column Sum of Squares Vector will be updated at each step.
C     When numerically necessary, these values will be recomputed.
C
C***SEE ALSO  DWNLIT
C***ROUTINES CALLED  IDAMAX
C***REVISION HISTORY  (YYMMDD)
C   790701  DATE WRITTEN
C   890620  Code extracted from WNLIT and made a subroutine.  (RWC))
C   900604  DP version created from SP version.  (RWC)
C***END PROLOGUE  DWNLT1
      INTEGER I, IMAX, IR, LEND, MDW, MEND
      DOUBLE PRECISION H(*), HBAR, SCALE(*), W(MDW,*)
      LOGICAL RECALC
C
      EXTERNAL IDAMAX
      INTEGER IDAMAX
C
      INTEGER J, K
C
C***FIRST EXECUTABLE STATEMENT  DWNLT1
      IF (IR.NE.1 .AND. (.NOT.RECALC)) THEN
C
C        Update column SS=sum of squares.
C
         DO 10 J=I,LEND
            H(J) = H(J) - SCALE(IR-1)*W(IR-1,J)**2
   10    CONTINUE
C
C        Test for numerical accuracy.
C
         IMAX = IDAMAX(LEND-I+1, H(I), 1) + I - 1
         RECALC = (HBAR+1.E-3*H(IMAX)) .EQ. HBAR
      ENDIF
C
C     If required, recalculate column SS, using rows IR through MEND.
C
      IF (RECALC) THEN
         DO 30 J=I,LEND
            H(J) = 0.D0
            DO 20 K=IR,MEND
               H(J) = H(J) + SCALE(K)*W(K,J)**2
   20       CONTINUE
   30    CONTINUE
C
C        Find column with largest SS.
C
         IMAX = IDAMAX(LEND-I+1, H(I), 1) + I - 1
         HBAR = H(IMAX)
      ENDIF
      RETURN
      END
*DECK DWNLT2
      LOGICAL FUNCTION DWNLT2 (ME, MEND, IR, FACTOR, TAU, SCALE, WIC)
C***BEGIN PROLOGUE  DWNLT2
C***SUBSIDIARY
C***PURPOSE  Subsidiary to WNLIT
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (WNLT2-S, DWNLT2-D)
C***AUTHOR  Hanson, R. J., (SNLA)
C           Haskell, K. H., (SNLA)
C***DESCRIPTION
C
C     To test independence of incoming column.
C
C     Test the column IC to determine if it is linearly independent
C     of the columns already in the basis.  In the initial tri. step,
C     we usually want the heavy weight ALAMDA to be included in the
C     test for independence.  In this case, the value of FACTOR will
C     have been set to 1.E0 before this procedure is invoked.
C     In the potentially rank deficient problem, the value of FACTOR
C     will have been set to ALSQ=ALAMDA**2 to remove the effect of the
C     heavy weight from the test for independence.
C
C     Write new column as partitioned vector
C           (A1)  number of components in solution so far = NIV
C           (A2)  M-NIV components
C     And compute  SN = inverse weighted length of A1
C                  RN = inverse weighted length of A2
C     Call the column independent when RN .GT. TAU*SN
C
C***SEE ALSO  DWNLIT
C***ROUTINES CALLED  (NONE)
C***REVISION HISTORY  (YYMMDD)
C   790701  DATE WRITTEN
C   890620  Code extracted from WNLIT and made a subroutine.  (RWC))
C   900604  DP version created from SP version.  (RWC)
C***END PROLOGUE  DWNLT2
      DOUBLE PRECISION FACTOR, SCALE(*), TAU, WIC(*)
      INTEGER IR, ME, MEND
C
      DOUBLE PRECISION RN, SN, T
      INTEGER J
C
C***FIRST EXECUTABLE STATEMENT  DWNLT2
      SN = 0.E0
      RN = 0.E0
      DO 10 J=1,MEND
         T = SCALE(J)
         IF (J.LE.ME) T = T/FACTOR
         T = T*WIC(J)**2
C
         IF (J.LT.IR) THEN
            SN = SN + T
         ELSE
            RN = RN + T
         ENDIF
   10 CONTINUE
      DWNLT2 = RN .GT. SN*TAU**2
      RETURN
      END
*DECK DWNLT3
      SUBROUTINE DWNLT3 (I, IMAX, M, MDW, IPIVOT, H, W)
C***BEGIN PROLOGUE  DWNLT3
C***SUBSIDIARY
C***PURPOSE  Subsidiary to WNLIT
C***LIBRARY   SLATEC
C***TYPE      DOUBLE PRECISION (WNLT3-S, DWNLT3-D)
C***AUTHOR  Hanson, R. J., (SNLA)
C           Haskell, K. H., (SNLA)
C***DESCRIPTION
C
C     Perform column interchange.
C     Exchange elements of permuted index vector and perform column
C     interchanges.
C
C***SEE ALSO  DWNLIT
C***ROUTINES CALLED  DSWAP
C***REVISION HISTORY  (YYMMDD)
C   790701  DATE WRITTEN
C   890620  Code extracted from WNLIT and made a subroutine.  (RWC))
C   900604  DP version created from SP version.  (RWC)
C***END PROLOGUE  DWNLT3
      INTEGER I, IMAX, IPIVOT(*), M, MDW
      DOUBLE PRECISION H(*), W(MDW,*)
C
      EXTERNAL DSWAP
C
      DOUBLE PRECISION T
      INTEGER ITEMP
C
C***FIRST EXECUTABLE STATEMENT  DWNLT3
      IF (IMAX.NE.I) THEN
         ITEMP        = IPIVOT(I)
         IPIVOT(I)    = IPIVOT(IMAX)
         IPIVOT(IMAX) = ITEMP
C
         CALL DSWAP(M, W(1,IMAX), 1, W(1,I), 1)
C
         T       = H(IMAX)
         H(IMAX) = H(I)
         H(I)    = T
      ENDIF
      RETURN
      END

max_issues_repo_name: bslakman/RMG-Java

max_forks_repo_name: bslakman/RMG-Java


text:
!----------------------------------------------------------------
! (c) Copyright, 2018 by the Regents of the University of California.
! FieldQuantclass: 3D field quantity class in Field module of APPLICATION
!                 layer.
! Version: 2.0
! Author: Ji Qiang, LBNL
! Description: This class defines a 3-D field quantity in the accelerator.
!              The field quantity can be updated at each step.
! Comments:
!----------------------------------------------------------------
      module FieldQuantclass
        use Timerclass
        use CompDomclass
        use Pgrid2dclass
        use FFTclass
        use Transposeclass
        use PhysConstclass
        use Dataclass
        type FieldQuant
!          private
          !# of mesh points in x and y directions.
          integer :: Nx,Ny,Nz,Nxlocal,Nylocal,Nzlocal
          !# field quantity array.
          double precision, pointer, dimension(:,:,:) :: FieldQ
        end type FieldQuant
      contains
        !Initialize field class.
        subroutine construct_FieldQuant(this,innx,inny,innz,geom,grid)
        implicit none
        include 'mpif.h'
        type (FieldQuant), intent(out) :: this
        type (CompDom), intent(in) :: geom
        type (Pgrid2d), intent(in) :: grid
        integer, intent(in) :: innx, inny, innz
        integer :: myid, myidx, myidy, nptot,nproccol,nprocrow
        integer, allocatable, dimension(:,:,:) :: LocalTable

        call getsize_Pgrid2d(grid,nptot,nproccol,nprocrow)
        call getpost_Pgrid2d(grid,myid,myidy,myidx)

        allocate(LocalTable(2,0:nprocrow-1,0:nproccol-1))
        call getlctabnm_CompDom(geom,LocalTable)

        this%Nx = innx
        this%Ny = inny
        this%Nz = innz
        this%Nxlocal = innx
        if(nproccol.gt.1) then
          this%Nylocal = LocalTable(2,myidx,myidy)+2
        else
          this%Nylocal = LocalTable(2,myidx,myidy)
        endif
        if(nprocrow.gt.1) then
          this%Nzlocal = LocalTable(1,myidx,myidy)+2
        else
          this%Nzlocal = LocalTable(1,myidx,myidy)
        endif

        allocate(this%FieldQ(this%Nxlocal,this%Nylocal,this%Nzlocal))
        this%FieldQ = 0.0

        deallocate(LocalTable)

        end subroutine construct_FieldQuant

        ! set field quantity.
        subroutine set_FieldQuant(this,innx,inny,innz,geom,grid,&
                                  nprx,npry)
        implicit none
        include 'mpif.h'
        type (FieldQuant), intent(inout) :: this
        type (CompDom), intent(in) :: geom
        type (Pgrid2d), intent(in) :: grid
        integer, intent(in) :: innx, inny, innz, nprx, npry
        integer :: myid, myidx, myidy
        integer, dimension(2,0:nprx-1,0:npry-1)::LocalTable

        call getpost_Pgrid2d(grid,myid,myidy,myidx)

        call getlctabnm_CompDom(geom,LocalTable)

        this%Nx = innx
        this%Ny = inny
        this%Nz = innz
        this%Nxlocal = innx
        if(npry.gt.1) then
          this%Nylocal = LocalTable(2,myidx,myidy)+2
        else
          this%Nylocal = LocalTable(2,myidx,myidy)
        endif
        if(nprx.gt.1) then
          this%Nzlocal = LocalTable(1,myidx,myidy)+2
        else
          this%Nzlocal = LocalTable(1,myidx,myidy)
        endif
        this%Nxlocal = innx

        deallocate(this%FieldQ)
        allocate(this%FieldQ(this%Nxlocal,this%Nylocal,this%Nzlocal))
        this%FieldQ = 0.0

        end subroutine set_FieldQuant

!----------------------------------------------------------------------
! update potential (solving Possion's equation) with 3D isolated
! boundary conditions.
!
        subroutine update3O_FieldQuant(this,source,fldgeom,grid,nxlc,&
          nylc,nzlc,nprocrow,nproccol,nylcr,nzlcr)
        implicit none
        include 'mpif.h'
        type (CompDom), intent(in) :: fldgeom
        double precision, dimension(nxlc,nylc,nzlc), intent(in) :: source
        type (Pgrid2d), intent(in) :: grid
        integer, intent(in) :: nxlc,nylc,nzlc,&
                               nprocrow,nproccol,nylcr,nzlcr
        type (FieldQuant), intent(inout) :: this
        double precision, dimension(3) :: msize
        double precision :: hx, hy, hz, temp
        integer, dimension(0:nprocrow-1) :: ypzstable,pztable
        integer, dimension(0:nproccol-1) :: xpystable,pytable
        double precision , dimension(nxlc,nylcr,nzlcr) :: rho
        integer :: myid,myidz,myidy,&
                   comm2d,commcol,commrow
        integer :: nxpylc2,nypzlc2
        integer :: nsxy1,nsxy2,nsyz1,nsyz2
        integer :: i,j,k,inxglb,inyglb,inzglb,innx,inny,innz
        integer, dimension(2,0:nprocrow-1,0:nproccol-1)::LocalTable
        integer, dimension(3) :: glmshnm
        integer :: jadd,kadd,ierr
        double precision :: t0
! Sherry
        double precision :: t1, t_openbc3d
        integer :: iam

        call getpost_Pgrid2d(grid,myid,myidy,myidz)
        call getcomm_Pgrid2d(grid,comm2d,commcol,commrow)

        call MPI_BARRIER(comm2d,ierr)
        call starttime_Timer(t0)

        call getlctabnm_CompDom(fldgeom,LocalTable)

        do i = 0, nprocrow-1
          pztable(i) = LocalTable(1,i,0)
        enddo
        do i = 0, nproccol-1
          pytable(i) = LocalTable(2,0,i)
        enddo

        call getmnum_CompDom(fldgeom,glmshnm)
        inxglb = glmshnm(1)
        inyglb = glmshnm(2)
        inzglb = glmshnm(3)

        innz = nzlcr
        inny = nylcr
        innx = nxlc

        if(nprocrow.gt.1) then
          kadd = 1
        else
          kadd = 0
        endif
        if(nproccol.gt.1) then
          jadd = 1
        else
          jadd = 0
        endif
        do k = 1, innz
          do j = 1, inny
            do i = 1, innx
              rho(i,j,k) = source(i,j+jadd,k+kadd)
            enddo
          enddo
        enddo

        ! +1 is from the real to complex fft.
        nsxy1 = (inxglb+1)/nproccol
        nsxy2 = (inxglb+1) - nproccol*nsxy1
        do i = 0, nproccol-1
          if(i.le.(nsxy2-1)) then
            xpystable(i) = nsxy1+1
          else
            xpystable(i) = nsxy1
          endif
        enddo

        nsyz1 = 2*inyglb/nprocrow
        nsyz2 = 2*inyglb - nprocrow*nsyz1
        do i = 0, nprocrow-1
          if(i.le.(nsyz2-1)) then
            ypzstable(i) = nsyz1+1
          else
            ypzstable(i) = nsyz1
          endif
        enddo

        nxpylc2 = xpystable(myidy)
        nypzlc2 = ypzstable(myidz)

        call getmsize_CompDom(fldgeom,msize)
        hx = msize(1)
        hy = msize(2)
        hz = msize(3)

! Sherry
!        call starttime_Timer(t1)

        ! Open boundary conditions!
        call openBC3D(innx,inny,innz,rho,hx,hy,hz, &
        nxpylc2,nypzlc2,myidz,myidy,nprocrow,nproccol,commrow,commcol,&
        comm2d,pztable,pytable,ypzstable,xpystable,&
        inxglb,inyglb,inzglb)

        do k = 1, innz
          do j = 1, inny
            do i = 1, innx
              this%FieldQ(i,j+jadd,k+kadd) = rho(i,j,k)
            enddo
          enddo
        enddo

        call MPI_BARRIER(comm2d,ierr)
        t_field = t_field + elapsedtime_Timer(t0)
! Sherry
!        call MPI_COMM_RANK(comm2d,iam,ierr)
!        t_openbc3d = elapsedtime_Timer(t1)
!        if ( iam == 0 ) then
!           print*, 'openBC3D time = ', t_openbc3d
!        endif

        end subroutine update3O_FieldQuant

        ! Solving Poisson's equation with open BCs.
        ! Sherry: now, innx, inny, innz are local dimensions
        subroutine openBC3D(innx,inny,innz,rho,hx,hy,hz,&
           nxpylc2,nypzlc2,myidz,myidy,npz,npy,commrow,commcol,comm2d, &
           pztable,pytable,ypzstable,xpystable,inxglb,&
           inyglb,inzglb)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: innx,inny,innz,inxglb,inyglb,inzglb
        integer, intent(in) :: nxpylc2,nypzlc2
        integer, intent(in) :: myidz,myidy,npz,npy,commrow,commcol,&
                               comm2d
        integer, dimension(0:npz-1), intent(in) :: pztable,ypzstable
        integer, dimension(0:npy-1), intent(in) :: pytable,xpystable
        double precision, dimension(innx,inny,innz), intent(inout) :: rho
        double precision, intent(in) :: hx, hy, hz
        double precision :: scalex,scaley,scalez
        integer :: i,j,k,n1,n2,n3,nylc22,nzlc22
        double complex, allocatable, dimension(:,:,:) :: rho2out
        double complex, allocatable, dimension(:,:,:) :: grn
        integer :: ginny,ginnz,ierr
! Sherry
        double precision :: t1, t_fft
        integer :: iam

        n1 = 2*inxglb
        n2 = 2*inyglb
        n3 = 2*inzglb

        nylc22 = nxpylc2
        nzlc22 = nypzlc2

        scalex = 1.0
        scaley = 1.0
        scalez = 1.0

        allocate(rho2out(n3,nylc22,nzlc22))

! Sherry
!        call starttime_Timer(t1)

        call fft3d1_FFT(n1,n2,n3,innz,inny,nylc22,nzlc22,&
                1,scalex,scaley,scalez,rho,ypzstable,pztable,&
        xpystable,pytable,npz,commrow,npy,commcol,comm2d,myidz,myidy, &
        rho2out)

! Sherry
!        call MPI_COMM_RANK(comm2d,iam,ierr)
!        t_fft = elapsedtime_Timer(t1)
!        if ( iam == 0 ) then
!           print*, 'fft3d1_FFT time = ', t_fft
!        endif

        !c compute FFT of the Green function on the grid:
        ! here the +1 is from the unsymmetry of green function
        ! on double-sized grid.
        if(myidz.eq.(npz-1)) then
           ginnz = innz + 1
        else
           ginnz = innz
        endif
        if(myidy.eq.(npy-1)) then
           ginny = inny + 1
        else
           ginny = inny
        endif
        allocate(grn(n3,nylc22,nzlc22))
        call greenf1tIntnew(inxglb,inyglb,inzglb,ginnz,ginny,nylc22,nzlc22, &
               hx,hy,hz,myidz,npz,commrow,myidy,npy,commcol,comm2d,&
                  ypzstable,pztable,xpystable,pytable,grn)

        ! multiply transformed charge density and transformed Green
        ! function:
        do k = 1, nzlc22
          do j = 1, nylc22
            do i = 1, n3
              rho2out(i,j,k) = rho2out(i,j,k)*grn(i,j,k)
            enddo
          enddo
        enddo

        deallocate(grn)

        ! inverse FFT:
        scalex = 1.0d0/dble(n1)
        scaley = 1.0d0/dble(n2)
        scalez = 1.0d0/dble(n3)
        call invfft3d1_FFT(n3,n2,n1,nylc22,nzlc22,inny,innz,&
               -1,scalex,scaley,scalez,rho2out,pztable,ypzstable,&
        pytable,xpystable,npz,commrow,npy,commcol,comm2d,myidz,myidy,&
        rho)

        deallocate(rho2out)

        end subroutine openBC3D

!----------------------------------------------------------------------
        ! green function for extended array.
        subroutine greenf1tIntnew(nx,ny,nz,nsizez,nsizey,nsizexy,nsizeyz,&
                  hx,hy,hz,myidx,npx,commrow,myidy,npy,commcol,comm2d,&
                   xstable,xrtable,ystable,yrtable,grnout)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: nx,ny,nz,nsizez,nsizey,nsizexy,nsizeyz
        integer, intent(in) :: myidx,myidy,npx,npy,commrow,commcol,&
                               comm2d
        integer, dimension(0:npx-1),intent(in) :: xstable,xrtable
        integer, dimension(0:npy-1),intent(in) :: ystable,yrtable
        double precision, intent(in) :: hx, hy, hz
        double complex, intent (out), &
                       dimension (2*nz,nsizexy,nsizeyz) :: grnout
        integer, dimension(0:npx-1) :: gxrtable
        integer, dimension(0:npy-1) :: gyrtable
        integer :: i,j,k,kk,ii,jj,iii,jjj,kkk,n1,n2,n3,ns1,ns2
        integer :: nblockx,nblocky,ksign,nxx,ierr
        double precision, dimension (nx+1,nsizey,nsizez) :: grn
        double precision, dimension (nx+2,nsizey+1,nsizez+1) :: grntmp
        double precision :: scalex,scaley,scalez
        double precision :: t0
        double precision, dimension(2*nx,nsizey) :: tmp1
        double complex, dimension(nx+1,nsizey) :: tmp10
        double complex, dimension(2*ny,nsizexy) :: tmp2
        double complex, dimension(2*nz,nsizexy) :: tmp3
        double complex, allocatable, dimension(:,:,:) :: x1
        double complex, allocatable, dimension(:,:,:) :: x0
        double precision :: rr,aa,bb,cc,dd,ee,ff,ss
        double complex :: gg,gg2
        double complex :: ggrr
        double precision, dimension(2) :: xx,yy,zz
        double precision, dimension(3) :: vv
        integer :: n,i0,j0,k0
        double precision :: recfourpi

        recfourpi = 1.0d0/(8.0d0*asin(1.0d0))

        call starttime_Timer(t0)

!        if(myidx.eq.0 .and. myidy.eq.0) then
!          print*,"into integrated Green function....."
!        endif
        gxrtable = xrtable
        gxrtable(npx-1) = xrtable(npx-1) + 1
        nblockx = 0
        do i = 0, myidx-1
          nblockx = nblockx + gxrtable(i)
        enddo

        gyrtable = yrtable
        gyrtable(npy-1) = yrtable(npy-1) + 1
        nblocky = 0
        do i = 0, myidy-1
          nblocky = nblocky + gyrtable(i)
        enddo

        do k0 = 1, nsizez+1
          do j0 = 1, nsizey+1
            do i0 = 1, nx+2
              jj = j0 + nblocky
              kk = k0 + nblockx
              kkk = kk - 1
              jjj = jj - 1
              iii = i0 - 1

              vv(1) = iii*hx-hx/2
              vv(2) = jjj*hy-hy/2
              vv(3) = kkk*hz-hz/2

                rr = sqrt(vv(1)**2+vv(2)**2+vv(3)**2)
                aa = vv(1)**2*vv(3)+(vv(2)**2+vv(3)**2)*vv(3) + &
                            vv(1)*vv(3)*rr
                bb = vv(1)**2*vv(2) + vv(1)*vv(2)*rr
                cc = vv(1)*(vv(1)**2+vv(2)**2)+vv(3)*vv(1)*(vv(3)+rr)
                dd = vv(3)*vv(2)*(vv(3) + rr)
                ee = vv(1)**2*vv(2)+vv(2)*(vv(2)**2+vv(3)*(vv(3)+rr))
                ff = vv(1)*vv(3)*(vv(3)+rr)
                ss = 4*vv(2)*vv(3)*log(vv(1)+rr) + &
                        4*vv(1)*vv(3)*log(vv(2)+rr) + &
                        4*vv(1)*vv(2)*log(vv(3)+rr)

                gg2 = cmplx(0.0,vv(3)**2)*log(cmplx(aa**2-bb**2,2*aa*bb)/ &
                        (aa**2+bb**2) ) + cmplx(0.0,vv(1)**2)*&
                        log(cmplx(cc**2-dd**2,2*cc*dd )/(cc**2+dd**2))+&
                        cmplx(0.0,vv(2)**2)*log(cmplx(ee**2-ff**2,2*ee*ff)/ &
                        (ee**2+ff**2) ) + ss

              !grntmp(i0,j0,k0) = recfourpi*real(gg2)/(4*hx*hy*hz) !wrong in Z code
              grntmp(i0,j0,k0) = real(gg2)/(4*hx*hy*hz)

            enddo
          enddo
        enddo

        do k0 = 1, nsizez
          do j0 = 1, nsizey
            do i0 = 1, nx+1
              grn(i0,j0,k0) = grntmp(i0+1,j0+1,k0+1)-grntmp(i0,j0+1,k0+1)-grntmp(i0+1,j0,k0+1)+&
                              grntmp(i0,j0,k0+1)-grntmp(i0+1,j0+1,k0)+grntmp(i0,j0+1,k0)+&
                              grntmp(i0+1,j0,k0)-grntmp(i0,j0,k0)
            enddo
          enddo
        enddo

!              xx(1) = iii*hx-hx/2
!              xx(2) = iii*hx+hx/2
!              yy(1) = jjj*hy-hy/2
!              yy(2) = jjj*hy+hy/2
!              zz(1) = kkk*hz-hz/2
!              zz(2) = kkk*hz+hz/2
!
!              !find the integrated Green function.
!              n = 0
!              do k = 1, 2
!                do j = 1, 2
!                  do i = 1, 2
!                    n = n+1
!                    rr(n) = sqrt(xx(i)**2+yy(j)**2+zz(k)**2)
!                    aa(n) = xx(i)**2*zz(k)+(yy(j)**2+zz(k)**2)*zz(k) + &
!                            xx(i)*zz(k)*rr(n)
!                    bb(n) = xx(i)**2*yy(j) + xx(i)*yy(j)*rr(n)
!                    cc(n) = xx(i)*(xx(i)**2+yy(j)**2)+zz(k)*xx(i)*(zz(k)+rr(n))
!                    dd(n) = zz(k)*yy(j)*(zz(k) + rr(n))
!                    ee(n) = xx(i)**2*yy(j)+yy(j)*(yy(j)**2+zz(k)*(zz(k)+rr(n)))
!                    ff(n) = xx(i)*zz(k)*(zz(k)+rr(n))
!                    ss(n) = 4*yy(j)*zz(k)*log(xx(i)+rr(n)) + &
!                            4*xx(i)*zz(k)*log(yy(j)+rr(n)) + &
!                            4*xx(i)*yy(j)*log(zz(k)+rr(n))
!                    gg(n) = cmplx(0.0,zz(k)**2)*log(cmplx(aa(n)**2-bb(n)**2,2*aa(n)*bb(n))/ &
!                            (aa(n)**2+bb(n)**2) ) + cmplx(0.0,xx(i)**2)*&
!                            log(cmplx(cc(n)**2-dd(n)**2,2*cc(n)*dd(n) )/(cc(n)**2+dd(n)**2))+&
!                            cmplx(0.0,yy(j)**2)*log(cmplx(ee(n)**2-ff(n)**2,2*ee(n)*ff(n))/ &
!                            (ee(n)**2+ff(n)**2) )
!                    gg2(n) = ss(n) +  gg(n)
!                  enddo
!                enddo
!              enddo
!              ggrr = (-gg2(1)+gg2(2)+gg2(3)-gg2(4)+gg2(5d0)-gg2(6)-gg2(7)+gg2(8))/4
!              grn(i0,j0,k0) = real(ggrr)/(hx*hy*hz)
!            enddo
!          enddo
!        enddo

        if((myidx.eq.0).and.(myidy.eq.0)) then
          if(nsizez.gt.1) then
            grn(1,1,1) = grn(1,1,2)
          else
            grn(1,1,1) = 1.0
          endif
        endif

        scalex = 1.0d0
        scaley = 1.0d0
        scalez = 1.0d0
        n1 = 2*nx
        n2 = 2*ny
        n3 = 2*nz
        ksign = 1

        nxx = n1/2 + 1
        !FFTs along y and z dimensions.
        allocate(x0(n1/2+1,nsizey,nsizez))
        do k = 1, nsizez
          do j = 1, nsizey
            do i = 1, n1/2+1
              tmp1(i,j) = grn(i,j,k)
            enddo
            do i = n1/2+2, n1
              tmp1(i,j) = grn(n1-i+2,j,k)
            enddo
          enddo

          ! FFTs along z dimensions:
          call fftrclocal_FFT(ksign,scalex,tmp1,n1,nsizey,tmp10)

          do j = 1, nsizey
            do i = 1, n1/2+1
              x0(i,j,k) = tmp10(i,j)
            enddo
          enddo
        enddo

        allocate(x1(n2/2+1,nsizexy,nsizez))

        ! FFTs along y dimensions:
!        call MPI_BARRIER(commcol,ierr)
        call trans3d_TRANSP(nxx,n2/2+1,nsizexy,nsizey,x0,x1,npy,&
                     ystable,gyrtable,commcol,nsizez)
        deallocate(x0)
        allocate(x0(n2,nsizexy,nsizez))

        do k = 1, nsizez
          do j = 1, nsizexy
            do i = 1, n2/2+1
              tmp2(i,j) = x1(i,j,k)
            enddo
            do i = n2/2+2,n2
              tmp2(i,j) = x1(n2-i+2,j,k)
            enddo
          enddo

          call fftlocal_FFT(ksign,scaley,tmp2,n2,nsizexy)

          do j = 1, nsizexy
            do i = 1, n2
              x0(i,j,k) = tmp2(i,j)
            enddo
          enddo
        enddo

        deallocate(x1)
        allocate(x1(n3/2+1,nsizexy,nsizeyz))
!        call MPI_BARRIER(commcol,ierr)
        call trans3d3_TRANSP(n2,nsizexy,nsizez,nsizeyz,x0,x1,npx,&
                      xstable,gxrtable,commrow,myidx,n3/2+1)
        deallocate(x0)

        do k = 1, nsizeyz
          do j = 1, nsizexy
            do i = 1, n3/2+1
              tmp3(i,j) = x1(i,j,k)
            enddo
            do i = n3/2+2,n3
              tmp3(i,j) = x1(n3-i+2,j,k)
            enddo
          enddo

          call fftlocal_FFT(ksign,scalez,tmp3,n3,nsizexy)

          do j = 1, nsizexy
            do i = 1, n3
              grnout(i,j,k) = tmp3(i,j)
            enddo
          enddo
        enddo

        deallocate(x1)

        t_greenf = t_greenf + elapsedtime_Timer(t0)

        end subroutine greenf1tIntnew
!--------------------------------------------------------------------

!-------------------------------------------------------------------------
        subroutine setval_FieldQuant(this,i,j,k,value)
        implicit none
        include 'mpif.h'
        type (FieldQuant), intent(out) :: this
        integer, intent(in) :: i, j, k
        double precision, intent(in) :: value

        this%FieldQ(i,j,k) = value

        end subroutine setval_FieldQuant

        double precision function get_FieldQuant(this,i,j,k)
        implicit none
        include 'mpif.h'
        type (FieldQuant), intent(in) :: this
        integer, intent(in) :: i, j, k

        get_FieldQuant = this%FieldQ(i,j,k)

        end function get_FieldQuant

        subroutine getglb_FieldQuant(this,temp)
        implicit none
        include 'mpif.h'
        type (FieldQuant), intent(in) :: this
        type (FieldQuant), intent(out) :: temp
        integer :: i, j, k, lcnz,lcny,lcnx
        double precision :: value

        lcnx = this%Nxlocal
        lcny = this%Nylocal
        lcnz = this%Nzlocal

        do k = 1, lcnz
          do j = 1, lcny
            do i = 1, lcnx
              value = get_FieldQuant(this,i,j,k)
              call setval_FieldQuant(temp,i,j,k,value)
            enddo
          enddo
        enddo

        end subroutine getglb_FieldQuant

        subroutine getlcgrid_FieldQuant(this,nxlc,nylc,nzlc)
        implicit none
        include 'mpif.h'
        type (FieldQuant), intent(in) :: this
        integer, intent(out) :: nxlc,nylc,nzlc

        nxlc = this%Nxlocal
        nylc = this%Nylocal
        nzlc = this%Nzlocal

        end subroutine

        subroutine destruct_FieldQuant(this)
        implicit none
        include 'mpif.h'
        type (FieldQuant), intent(out) :: this

        deallocate(this%FieldQ)

        end subroutine destruct_FieldQuant

!-------------------------------------------------------------------------
       !longitudinal and transverse wakefield
       subroutine wakefield_FieldQuant(Nz,xwakez,ywakez,recvdensz,exwake,eywake,ezwake,&
                                       hz,aa,gg,leng)
       implicit none
       include 'mpif.h'
       integer, intent(in) :: Nz
       double precision, intent(in) :: hz, aa, gg, leng
       double precision, dimension(Nz,2), intent(in) :: recvdensz
       double precision, dimension(Nz), intent(in) :: xwakez,ywakez
       double precision, dimension(Nz), intent(out) :: exwake,ezwake,eywake
       double precision, dimension(2*Nz,1) :: densz2n,densz2nout,&
                         greenwake,greenwakeout
       integer :: kz,twonz,one,ksign,kkzz,i
       double precision :: scale,zz00,pilc,Z0,alpha1,alpha,zz
       double precision :: coef1,coef2,densconst,offset,zzmax
       real*8 :: tmptmp

       pilc = 2*asin(1.0)

       do kz = 1, Nz
          densz2n(kz,1) = recvdensz(kz,1)
       enddo
       do kz = Nz + 1, 2*Nz
          densz2n(kz,1) = 0.0d0
       enddo

       twonz = 2*Nz
       one = 1
       ksign = 1
       scale = 1.d0
       call fftrclocal2_FFT(ksign,scale,densz2n,twonz,one,densz2nout)

       !longitudinal wakefield function
       !the following longitudinal wakefield function is from
       !"Short-range dipole wakefields in accelerating structures for the NLC",
       !K. L. F. Bane, SLAC-PUB-9663, 2003.
       pilc = 2*asin(1.0d0)
       alpha1 = 0.4648
       alpha = 1.-alpha1*sqrt(gg/leng)-(1.-2*alpha1)*(gg/leng)
       !slac wake
       !zz00 = gg*(aa/(alpha*leng))**2/8
       !fermi wake
       zz00 = 0.41*(aa/leng)**0.8*(gg/leng)**1.6*aa
       Z0 = 120*pilc
       !greenwake(1,1) = 0.0
       zz = 0.0d0
       !The 0.5 is from S+
       greenwake(1,1) = 0.5d0*Z0*Clight*exp(-sqrt(zz/zz00))/(pilc*aa*aa)
       !do kz = 1, Nz+1
       do kz = 2, Nz+1
         zz = (kz-1)*hz
         greenwake(kz,1) = Z0*Clight*exp(-sqrt(zz/zz00))/(pilc*aa*aa)
       enddo
       do kz = Nz+2, twonz
         greenwake(kz,1) = greenwake(twonz-kz+2,1)
       enddo
       do kz = 2, Nz
         greenwake(kz,1) = 0.0
       enddo

       call fftrclocal2_FFT(ksign,scale,greenwake,twonz,one,&
             greenwakeout)

       !do kz = 1, twonz
       !  greenwake(kz,1) = densz2nout(kz,1)*greenwakeout(kz,1)
       !enddo
       do kz = 1, 2
          greenwake(kz,1) = densz2nout(kz,1)*greenwakeout(kz,1)
       enddo
       do kz = 2, twonz/2
          greenwake(2*kz-1,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz-1,1)-&
                            densz2nout(2*kz,1)*greenwakeout(2*kz,1)
          greenwake(2*kz,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz,1)+&
                            densz2nout(2*kz,1)*greenwakeout(2*kz-1,1)
       enddo

       scale = 1.0/twonz
       ksign = -1
       call fftcrlocal2_FFT(ksign,scale,greenwake,twonz,one,&
            greenwakeout)

       !the "-" sign is from definition of longitudinal wake function, see Wangler's book
       do kz = 1, Nz
         ezwake(kz) = -greenwakeout(kz,1)*hz
       enddo

!------------------------------------------------------------------------------
       !calculate the transverse wakefield effects
       do kz = 1, Nz
          densz2n(kz,1) = recvdensz(kz,1)*xwakez(kz)
       enddo
       do kz = Nz + 1, 2*Nz
          densz2n(kz,1) = 0.0d0
       enddo

       twonz = 2*Nz
       one = 1
       ksign = 1
       !scale = 1./(twonz)
       scale = 1.
       call fftrclocal2_FFT(ksign,scale,densz2n,twonz,one,densz2nout)

       !tranverse wakefield
       !the following transverse wakefield function is from
       !"Short-range dipole wakefields in accelerating structures for the NLC",
       !K. L. F. Bane, SLAC-PUB-9663, 2003. here, 1.1 is an average as suggested       !on page 11 for cell 45.
       !for LCLS slac
       !zz00 = 1.1*0.169*aa*(aa/leng)**1.17*(gg/aa)**0.38
       !for Fermi Elettra
       zz00 = 1.0*0.169*aa*(aa/leng)**1.17*(gg/aa)**0.38
       coef1 = 4*Z0*Clight*zz00/(pilc*aa*aa*aa*aa)
       !densconst = 0.5e-6
       !offset = 0.001
       !zzmax = 0.002
       !coef2 = coef1/zz00*densconst*offset
       !print*,"aa wake: ",aa,leng,gg,zz00,Z0,coef1,coef2,offset,densconst
       greenwake(1,1) = 0.0
       do kz = 1, Nz+1
         zz = (kz-1)*hz
         greenwake(kz,1) = coef1*(1.0-(1.+sqrt(zz/zz00))*exp(-sqrt(zz/zz00)))
       enddo

       do kz = Nz+2, 2*Nz
         greenwake(kz,1) = greenwake(twonz-kz+2,1)
       enddo
       do kz = 1, Nz
         greenwake(kz,1) = 0.0
       enddo

       ksign = 1
       !scale = 1./(twonz)
       scale = 1.
       call fftrclocal2_FFT(ksign,scale,greenwake,twonz,one,&
                greenwakeout)

       do kz = 1, 2
          greenwake(kz,1) = densz2nout(kz,1)*greenwakeout(kz,1)
       enddo
       do kz = 2, twonz/2
          greenwake(2*kz-1,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz-1,1)-&
                            densz2nout(2*kz,1)*greenwakeout(2*kz,1)
          greenwake(2*kz,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz,1)+&
                            densz2nout(2*kz,1)*greenwakeout(2*kz-1,1)
       enddo

       !scale = 1.0
       scale = 1.0/twonz
       ksign = -1
       call fftcrlocal2_FFT(ksign,scale,greenwake,twonz,one,&
                   densz2nout)

       do kz = 1, Nz
         exwake(kz) = densz2nout(kz,1)*hz
       enddo

       !vertical "Y" wakefields
       do kz = 1, Nz
          densz2n(kz,1) = recvdensz(kz,1)*ywakez(kz)
       enddo
       do kz = Nz + 1, 2*Nz
          densz2n(kz,1) = 0.0d0
       enddo

       twonz = 2*Nz
       one = 1
       ksign = 1
       !scale = 1./(twonz)
       scale = 1.
       call fftrclocal2_FFT(ksign,scale,densz2n,twonz,one,densz2nout)

       do kz = 1, 2
          greenwake(kz,1) = densz2nout(kz,1)*greenwakeout(kz,1)
       enddo
       do kz = 2, twonz/2
          greenwake(2*kz-1,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz-1,1)-&
                            densz2nout(2*kz,1)*greenwakeout(2*kz,1)
          greenwake(2*kz,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz,1)+&
                            densz2nout(2*kz,1)*greenwakeout(2*kz-1,1)
       enddo

       !scale = 1.0
       scale = 1.0/twonz
       ksign = -1
       call fftcrlocal2_FFT(ksign,scale,greenwake,twonz,one,&
                   greenwakeout)

       do kz = 1, Nz
         eywake(kz) = greenwakeout(kz,1)*hz
       enddo

       end subroutine wakefield_FieldQuant

!-------------------------------------------------------------------------
       !longitudinal and transverse wakefield using read in discrete data
       subroutine wakefieldread_FieldQuant(Nz,xwakez,ywakez,recvdensz,exwake,eywake,ezwake,&
                                       hz,aa,gg,leng,flagbtw)
       implicit none
       include 'mpif.h'
       integer, intent(in) :: Nz,flagbtw
       double precision, intent(in) :: hz, aa, gg, leng
       double precision, dimension(Nz,2), intent(in) :: recvdensz
       double precision, dimension(Nz), intent(in) :: xwakez,ywakez
       double precision, dimension(Nz), intent(out) :: exwake,ezwake,eywake
       double precision, dimension(2*Nz,1) :: densz2n,densz2nout,&
                         greenwake,greenwakeout
       integer :: kz,twonz,one,ksign,kkzz,i,iwk,iwk1
       double precision :: scale,zz00,pilc,Z0,alpha1,alpha,zz
       double precision :: coef1,coef2,densconst,offset,zzmax
       real*8 :: tmptmp,pbtw1,pbtw2,pbtw3,pbtw4
       real*8 :: leng1,leng2,hwk,zwk,ab

       pilc = 2*dasin(1.0d0)

       if(flagbtw.eq.4) then !skip longitudinal wake
         ezwake = 0.0d0
         goto 50
       endif

       do kz = 1, Nz
          densz2n(kz,1) = recvdensz(kz,1)
       enddo
       do kz = Nz + 1, 2*Nz
          densz2n(kz,1) = 0.0d0
       enddo

       twonz = 2*Nz
       one = 1
       ksign = 1
       scale = 1.0d0
       call fftrclocal2_FFT(ksign,scale,densz2n,twonz,one,densz2nout)

       hwk = zdatwk(2)-zdatwk(1) !assumed uniform grid
       do kz = 1, Nz+1
         zz = (kz-1)*hz
         zwk = zz - zdatwk(1)
         iwk = zwk/hwk + 1
         iwk1 = iwk+1
         ab = iwk*hwk-zwk
         if(kz.eq.1) then
           greenwake(kz,1) = 0.5d0*(edatwk(iwk)*ab + edatwk(iwk1)*(1.0d0-ab))
         else
           greenwake(kz,1) = edatwk(iwk)*ab + edatwk(iwk1)*(1.0d0-ab)
         endif
       enddo

       do kz = Nz+2, twonz
         greenwake(kz,1) = greenwake(twonz-kz+2,1)
       enddo
       do kz = 2, Nz
         greenwake(kz,1) = 0.0
       enddo

       call fftrclocal2_FFT(ksign,scale,greenwake,twonz,one,&
             greenwakeout)

       !do kz = 1, twonz
       !  greenwake(kz,1) = densz2nout(kz,1)*greenwakeout(kz,1)
       !enddo
       do kz = 1, 2
          greenwake(kz,1) = densz2nout(kz,1)*greenwakeout(kz,1)
       enddo
       do kz = 2, twonz/2
          greenwake(2*kz-1,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz-1,1)-&
                            densz2nout(2*kz,1)*greenwakeout(2*kz,1)
          greenwake(2*kz,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz,1)+&
                            densz2nout(2*kz,1)*greenwakeout(2*kz-1,1)
       enddo

       scale = 1.0d0/twonz
       ksign = -1
       call fftcrlocal2_FFT(ksign,scale,greenwake,twonz,one,&
            greenwakeout)

       !the "-" sign is from definition of longitudinal wake function, see Wangler's book
       do kz = 1, Nz
         ezwake(kz) = -greenwakeout(kz,1)*hz
       enddo

50     continue

       if((flagbtw.eq.2) .or. (flagbtw.eq.3)) then !no transverse wake for standing wave cavity so far
         exwake = 0.0d0
         eywake = 0.0d0
         goto 100
       endif

!------------------------------------------------------------------------------
       !calculate the transverse wakefield effects
       do kz = 1, Nz
          densz2n(kz,1) = recvdensz(kz,1)*xwakez(kz)
       enddo
       do kz = Nz + 1, 2*Nz
          densz2n(kz,1) = 0.0d0
       enddo

       twonz = 2*Nz
       one = 1
       ksign = 1
       !scale = 1./(twonz)
       scale = 1.0d0
       call fftrclocal2_FFT(ksign,scale,densz2n,twonz,one,densz2nout)

       hwk = zdatwk(2)-zdatwk(1) !assumed uniform grid
       do kz = 1, Nz+1
         zz = (kz-1)*hz
         zwk = zz - zdatwk(1)
         iwk = zwk/hwk + 1
         iwk1 = iwk+1
         ab = iwk*hwk-zwk
         if(kz.eq.1) then
           greenwake(kz,1) = 0.5d0*(epdatwk(iwk)*ab + epdatwk(iwk1)*(1.0d0-ab))
         else
           greenwake(kz,1) = epdatwk(iwk)*ab + epdatwk(iwk1)*(1.0d0-ab)
         endif
       enddo
       do kz = Nz+2, 2*Nz
         greenwake(kz,1) = greenwake(twonz-kz+2,1)
       enddo
       do kz = 1, Nz
         greenwake(kz,1) = 0.0
       enddo

       ksign = 1
       !scale = 1./(twonz)
       scale = 1.0d0
       call fftrclocal2_FFT(ksign,scale,greenwake,twonz,one,&
                greenwakeout)

       do kz = 1, 2
          greenwake(kz,1) = densz2nout(kz,1)*greenwakeout(kz,1)
       enddo
       do kz = 2, twonz/2
          greenwake(2*kz-1,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz-1,1)-&
                            densz2nout(2*kz,1)*greenwakeout(2*kz,1)
          greenwake(2*kz,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz,1)+&
                            densz2nout(2*kz,1)*greenwakeout(2*kz-1,1)
       enddo

       !scale = 1.0
       scale = 1.0d0/twonz
       ksign = -1
       call fftcrlocal2_FFT(ksign,scale,greenwake,twonz,one,&
                   densz2nout)

       do kz = 1, Nz
         exwake(kz) = densz2nout(kz,1)*hz
       enddo

       !vertical "Y" wakefields
       do kz = 1, Nz
          densz2n(kz,1) = recvdensz(kz,1)*ywakez(kz)
       enddo
       do kz = Nz + 1, 2*Nz
          densz2n(kz,1) = 0.0d0
       enddo

       twonz = 2*Nz
       one = 1
       ksign = 1
       !scale = 1./(twonz)
       scale = 1.0d0
       call fftrclocal2_FFT(ksign,scale,densz2n,twonz,one,densz2nout)

       hwk = zdatwk(2)-zdatwk(1) !assumed uniform grid
       do kz = 1, Nz+1
         zz = (kz-1)*hz
         zwk = zz - zdatwk(1)
         iwk = zwk/hwk + 1
         iwk1 = iwk+1
         ab = iwk*hwk-zwk
         if(kz.eq.1) then
           greenwake(kz,1) = 0.5d0*(eppdatwk(iwk)*ab + eppdatwk(iwk1)*(1.0d0-ab))
         else
           greenwake(kz,1) = eppdatwk(iwk)*ab + eppdatwk(iwk1)*(1.0d0-ab)
         endif
       enddo
       do kz = Nz+2, 2*Nz
         greenwake(kz,1) = greenwake(twonz-kz+2,1)
       enddo
       do kz = 1, Nz
         greenwake(kz,1) = 0.0
       enddo

       ksign = 1
       !scale = 1./(twonz)
       scale = 1.0d0
       call fftrclocal2_FFT(ksign,scale,greenwake,twonz,one,&
                greenwakeout)

       do kz = 1, 2
          greenwake(kz,1) = densz2nout(kz,1)*greenwakeout(kz,1)
       enddo
       do kz = 2, twonz/2
          greenwake(2*kz-1,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz-1,1)-&
                            densz2nout(2*kz,1)*greenwakeout(2*kz,1)
          greenwake(2*kz,1) = densz2nout(2*kz-1,1)*greenwakeout(2*kz,1)+&
                            densz2nout(2*kz,1)*greenwakeout(2*kz-1,1)
       enddo

       !scale = 1.0
       scale = 1.0d0/twonz
       ksign = -1
       call fftcrlocal2_FFT(ksign,scale,greenwake,twonz,one,&
                   greenwakeout)

       do kz = 1, Nz
         eywake(kz) = greenwakeout(kz,1)*hz
       enddo

100    continue

       end subroutine wakefieldread_FieldQuant

        !This subroutine calculates the 1d csr wakefield including
        !entrance, stead-state, and transitions effects.
        !Here, hx, rho,...are in real units. The return ezwake is in V/m.
        !The current version uses IGF corresponding to the four cases of Saldin et al.
        subroutine csrwakeTrIGF_FieldQuant(Nx,r0,ptmin,hx,blength,rhonew,rhonewp,&
                             rhonewpp,gam,ezwake)
        implicit none
        include 'mpif.h'
        integer, intent(in) :: Nx
        real*8 :: r0,ptmin,hx,blength,gam
        real*8, dimension(Nx) :: rhonew,rhonewp,rhonewpp,ezwake
        real*8 :: xx,xxl,xx2,dx,tmprho,epstol,deltas,deltasmax,xxbar,psi,&
                  phim,pilc,yy,hx2,xconst,xk,psitmp,xslpN
        integer :: Ni,il,i,Nmax,j,islpN,islpN1
        integer :: myidlc,ierr,islp,islp1,islp0,jstart
        real*8 :: aa,uuh,uuh24,xslp,csrss,xxp,ssh,xk2,tcoef
        real*8 :: bb,cc,yh,phih,phpy2,csrtr,xblg,xxbarh,phimh,psimax
        real*8 :: psip2x2,psipx2,csrdrm1,csrdr1,phpxy2,phpy,csrss1,csrss2,&
                  csrtr1,csrtr2,csrdrm2,csrdr2
!        real*8 :: IcsrCaseA,IcsrCaseB,IcsrCaseC,IcsrCaseD

!        call MPI_COMM_RANK(MPI_COMM_WORLD, myidlc, ierr)

!        epstol = 1.0d-10 !tolerance for root finding
        epstol = 2.0d-9
        Nmax = 100 !maximum # of iteration in root finding
        pilc = 2*asin(1.0d0)
        phim = blength/r0
        xconst = 1./(4*pilc*8.854187817d-12)
        xk2 = gam/r0

        psitmp = 0.0
        ezwake(1) = 0.0
        do i = 2, Nx
           xx = ptmin + (i-1)*hx
           xxl = (xx/r0)**3*r0/24

           Ni = i
           ezwake(i) = 0.0

           xblg = (i-1)*hx

           !if(myidlc.eq.0) print*,"ii:",i,xxl,xx

           if(xx.le.0) then
           else if(xx.le.blength) then

             xslp = (xx/r0)*r0/2/gam/gam + xxl
             islp0 = (xx-xslp-ptmin)/hx
             islp1 = islp0 + 1

             !Case A
             phih = xx/r0*gam
             ! IGF integral over the interior sample points for Case A:
             do j = 2, islp1-1
                 !write(*,*) 'Inside Case A!'
                 xxp = ptmin + (j-1)*hx + hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrtr2 = IcsrCaseA(phih,ssh,xk2)
                 xxp = ptmin + (j-1)*hx - hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrtr1 = IcsrCaseA(phih,ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(j)*(csrtr1-csrtr2)
             enddo
             if(islp1.ge.2) then
             !Add the upper end IGF integral for case A
                 xxp = xx - xslp
                 ssh = (xx-xxp)*gam**3/r0
                 csrtr2 = IcsrCaseA(phih,ssh,xk2)
                 xxp = ptmin + (islp1-1)*hx - hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrtr1 = IcsrCaseA(phih,ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(islp1)*(csrtr1-csrtr2)
             !Add the lower end IGF integral for case A
                 xxp = ptmin + hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrtr2 = IcsrCaseA(phih,ssh,xk2)
                 xxp = ptmin
                 ssh = (xx-xxp)*gam**3/r0
                 csrtr1 = IcsrCaseA(phih,ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(1)*(csrtr1-csrtr2)
             ! Special case
             else if(islp1.eq.1) then
                 xxp = xx - xslp
                 ssh = (xx-xxp)*gam**3/r0
                 csrtr2 = IcsrCaseA(phih,ssh,xk2)
                 xxp = ptmin + (islp1-1)*hx
                 ssh = (xx-xxp)*gam**3/r0
                 csrtr1 = IcsrCaseA(phih,ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(islp1)*(csrtr1-csrtr2)
             endif

             ! Case B (steady-state regime)
             jstart = max(islp1,1)
             ! IGF integral over the interior sample points for Case B:
             do j = jstart+2,Ni-1
                 xxp = ptmin + (j-1)*hx + hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrss2 = IcsrCaseB(ssh,xk2)
                 xxp = ptmin + (j-1)*hx - hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrss1 = IcsrCaseB(ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(j)*(csrss1-csrss2)
             enddo
             !add end integrals
             if(jstart.le.Ni-2) then
             !Add the upper end IGF integral for case B
                 xxp = ptmin + (Ni-1)*hx
                 ssh = (xx-xxp)*gam**3/r0
                 csrss2 = IcsrCaseB(ssh,xk2)
                 xxp = ptmin + (Ni-1)*hx - hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrss1 = IcsrCaseB(ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(Ni)*(csrss1-csrss2)
             !Add the lower end IGF integral for case B
                 if(islp1.lt.1) then
                   j = 1
                   xxp = ptmin + (j-1)*hx + hx/2
                   ssh = (xx-xxp)*gam**3/r0
                   csrss2 = IcsrCaseB(ssh,xk2)
                   xxp = ptmin + (j-1)*hx
                   ssh = (xx-xxp)*gam**3/r0
                   csrss1 = IcsrCaseB(ssh,xk2)
                   ezwake(i) = ezwake(i) + rhonew(j)*(csrss1-csrss2)
                 else
                   j = islp1+1
                   xxp = ptmin + (j-1)*hx + hx/2
                   ssh = (xx-xxp)*gam**3/r0
                   csrss2 = IcsrCaseB(ssh,xk2)
                   xxp = xx-xslp
                   ssh = (xx-xxp)*gam**3/r0
                   csrss1 = IcsrCaseB(ssh,xk2)
                   ezwake(i) = ezwake(i) + rhonew(j)*(csrss1-csrss2)
                 endif
             else if(jstart.eq.Ni-1) then
                 j = Ni
                 xxp = ptmin + (j-1)*hx
                 ssh = (xx-xxp)*gam**3/r0
                 csrss2 = IcsrCaseB(ssh,xk2)
                 xxp = xx-xslp
                 ssh = (xx-xxp)*gam**3/r0
                 csrss1 = IcsrCaseB(ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(j)*(csrss1-csrss2)
             endif

           !  ezwake(i) = ezwake(i)*xconst

             !if(myidlc.eq.0) print*,"xxl: ",xxl,xx,ezwake(i),r0,ptmin,phim
             !if(myidlc.eq.0) print*,"xxl: ",xx,ezwake(i),xxl,xx-ptmin,xslp,hx

           else

             !if(myidlc.eq.0) print*,"CD:",i,islp,xslp,hx,xx,ezwake(i)

             xxbar = (xx - blength)
             xslp= (r0*phim+xxbar)/2/gam/gam + r0*phim**3/24*&
                   (r0*phim+4*xxbar)/(r0*phim+xxbar)
             islp0 = (xx-xslp-ptmin)/hx
             islp1 = islp0 + 1
             xslpN = xxbar/2.d0/gam/gam
             islpN = (xx-xslpN-ptmin)/hx
             islpN1 = islpN + 1

             ! Case C
             xxbarh = xxbar*gam/r0
             phimh = phim*gam
             ! IGF integral over the interior sample points for Case C:
             do j = 2, islp1-1
                 xxp = ptmin + (j-1)*hx + hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrdrm2 = IcsrCaseC(phimh,xxbarh,ssh,xk2)
                 xxp = ptmin + (j-1)*hx - hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrdrm1 = IcsrCaseC(phimh,xxbarh,ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(j)*(csrdrm1-csrdrm2)
             enddo
             if(islp1.ge.2) then
             !Add the upper end IGF integral for case C
                 xxp = xx - xslp
                 ssh = (xx-xxp)*gam**3/r0
                 csrdrm2 = IcsrCaseC(phimh,xxbarh,ssh,xk2)
                 xxp = ptmin + (islp1-1)*hx - hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrdrm1 = IcsrCaseC(phimh,xxbarh,ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(islp1)*(csrdrm1-csrdrm2)
             !Add the lower end IGF integral for case C
                 xxp = ptmin + hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrdrm2 = IcsrCaseC(phimh,xxbarh,ssh,xk2)
                 xxp = ptmin
                 ssh = (xx-xxp)*gam**3/r0
                 csrdrm1 = IcsrCaseC(phimh,xxbarh,ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(1)*(csrdrm1-csrdrm2)
             ! Special case
             else if(islp1.eq.1) then
                 xxp = xx - xslp
                 ssh = (xx-xxp)*gam**3/r0
                 csrdrm2 = IcsrCaseC(phimh,xxbarh,ssh,xk2)
                 xxp = ptmin
                 ssh = (xx-xxp)*gam**3/r0
                 csrdrm1 = IcsrCaseC(phimh,xxbarh,ssh,xk2)
                 ezwake(i) = ezwake(i) + rhonew(islp1)*(csrdrm1-csrdrm2)
             endif

             !Case D
             psimax = phim*gam
             xxbarh = xxbar*gam/r0
             jstart = max(islp1,0)
             !write(13,*) 'xslp and xslpN',xslp,xslpN
             ! IGF integral over the interior sample points for Case D:
             do j = jstart+2,islpN1-1
                 xxp = ptmin + (j-1)*hx + hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrdr2 = IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
                 xxp = ptmin + (j-1)*hx - hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrdr1 = IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
                 ezwake(i) = ezwake(i) + rhonew(j)*(csrdr1-csrdr2)
             enddo
             if(islpN1.ge.(jstart+2)) then
             !Add the upper end IGF integral for case D
                 ssh = xslpN*gam**3/r0
                 csrdr2 = IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
                 xxp = ptmin + (islpN1-1)*hx - hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrdr1 = IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
                 ezwake(i) = ezwake(i) + rhonew(islpN1)*(csrdr1-csrdr2)
             !Add the lower end IGF integral for case D
                 xxp = ptmin + (jstart)*hx + hx/2
                 ssh = (xx-xxp)*gam**3/r0
                 csrdr2 = IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
                 ssh = xslp*gam**3/r0
                 csrdr1 = IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
                 ezwake(i) = ezwake(i) + rhonew(jstart+1)*(csrdr1-csrdr2)
             else if(islpN1.eq.(jstart+1)) then
                 ssh = xslpN*gam**3/r0
                 csrdr2 = IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
                 ssh = xslp*gam**3/r0
                 csrdr1 = IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
                 ezwake(i) = ezwake(i) + rhonew(jstart+1)*(csrdr1-csrdr2)
             endif

           endif
           ezwake(i) = ezwake(i)*xconst
        enddo

        end subroutine csrwakeTrIGF_FieldQuant

           function IcsrCaseA(phih,ssh,xk2)
           implicit none
           double precision:: phih,ssh,xk2,IcsrCaseA
           double precision:: bb,cc,yh,phpy
           bb = 2*phih+phih**3/3-2*ssh
           cc = phih**2+phih**4/12-2*ssh*phih
           yh = (-bb+sqrt(bb*bb-4*cc))/2
           phpy = (phih+yh)
           IcsrCaseA = xk2*(-(2*phpy+phih**3)/&
                    (phpy**2+phih**4/4)+1.0d0/ssh)
           end function IcsrCaseA

           function IcsrCaseB(ssh,xk2)
           implicit none
           double precision:: ssh,xk2,IcsrCaseB
           double precision:: aa,uuh
           aa = sqrt(64.0+144.0d0*ssh**2)
           uuh = (aa+12*ssh)**(1.0d0/3.0d0)-(aa-12*ssh)**(1.0d0/3.0d0)
           IcsrCaseB = xk2*(-4*uuh*(uuh**2+8)/ &
                    ((uuh**2+4.0d0)*(uuh**2+12.0d0)))
           end function IcsrCaseB

           function IcsrCaseC(phimh,xxbarh,ssh,xk2)
           implicit none
           double precision:: phimh,xxbarh,ssh,xk2,IcsrCaseC
           double precision:: bb,cc,yh,phpxya,phpxyb
           bb = 2*(phimh+xxbarh)-2*ssh+phimh**3/3+phimh**2*xxbarh
           cc = (phimh+xxbarh)**2+phimh**2*(phimh**2+4*phimh*xxbarh)/12-&
                      2*ssh*(phimh+xxbarh)
           yh = (-bb+sqrt(bb*bb-4*cc))/2
           phpxya = (phimh+xxbarh+yh)
           phpxyb = phimh*xxbarh+phimh*phimh/2.d0
           IcsrCaseC = xk2*(-2.d0*(phpxya+phimh*phpxyb)/ &
                          (phpxya**2+phpxyb**2)+1.0d0/ssh)
           end function IcsrCaseC

           function IcsrCaseD(xxbarh,psimax,ssh,xk2,epstol,Nmax)
           implicit none
           double precision:: xxbarh,ssh,xk2,epstol,IcsrCaseD
           double precision:: psi,psipxa,psipxb,psimax,x1,x2
           integer:: Nmax
           x1 = -epstol
           x2 = psimax*(1.d0+epstol)
           call root(x1,x2,epstol,Nmax,ssh,xxbarh,psi)
           psipxa = (psi+xxbarh)
           psipxb = psi*(xxbarh+psi/2.d0)
           IcsrCaseD = xk2*(-2.d0*(psipxa+psi*psipxb)/&
                    (psipxa**2+psipxb**2)+1.0d0/ssh)
           end function IcsrCaseD

     subroutine root(x1,x2,xacc,maxit,zeta,xxh,rtsafe)
!************************************************************************
!  This routine computes the root of the function that is evaluated in
!  the subroutine 'funcd'. It is based on the subroutine 'root' of
!  Numerical Recipes 9.4, which makes use of a Newton-Raphson method
!  with root bracketing.  It has been modified to handle the two bracket
!  endpoints carefully. The routine searches for a root in the interval
!  [x1,x2] with a tolerance given by 'xacc', and returns this value
!  as 'rtsafe'.  The maximum number of iterations allowed is 'maxit'.
!  C.E.M.
!***********************************************************************
     implicit none
     double precision:: rtsafe,x1,x2,xacc
     double precision:: xxh,zeta
     integer:: j,maxit
     double precision:: df,dx,dxold,f,fh,fl,temp,xh,xl
     call funcd(x1,xxh,zeta,fl,df)
     call funcd(x2,xxh,zeta,fh,df)
     if((fl>0.d0.and.fh>0.d0).or.(fl<0.d0.and.fh<0.d0)) then
           pause 'root must be bracketed in rtsafe'
           write(*,*) 'psimax,fl,fh = ',x2,fl,fh
     endif
     if(dabs(fl)< xacc) then
       rtsafe=x1
       return
     else if(dabs(fh)< xacc) then
       rtsafe=x2
       return
     else if(fl<0.d0) then
       xl=x1
       xh=x2
     else
       xh=x1
       xl=x2
     endif
     rtsafe=0.5d0*(x1+x2)
     dxold=dabs(x2-x1)
     dx=dxold
     call funcd(rtsafe,xxh,zeta,f,df)
     do j=1,maxit
        if(((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f)>0.d0.or. &
&         dabs(2.d0*f)>dabs(dxold*df)) then
          dxold=dx
          dx=0.5d0*(xh-xl)
          rtsafe=xl+dx
          if(xl==rtsafe) return
        else
          dxold=dx
          dx=f/df
          temp=rtsafe
          rtsafe=rtsafe-dx
          if(temp==rtsafe) return
        endif
        if(abs(dx)<xacc) return
        call funcd(rtsafe,xxh,zeta,f,df)
        if(f<0.d0) then
           xl=rtsafe
        else
           xh=rtsafe
        endif
     enddo
     pause 'root finding exceeding maximum iterations'
     return
     end subroutine

     subroutine funcd(psi,xxh,deltas,f,derivf)
!**********************************************************
!  This routine evaluates the function whose root produces
!  the retarded angle psi that is required for evaluating
!  the CSR kernel in Case D of Saldin et al.  The value
!  of the function is output as 'f', and its derivative
!  is output as 'derivf'.  C.E.M.
!*********************************************************
     implicit none
     double precision:: deltas,psi,term1,xxh,gamma,f,derivf
     double precision:: alpha,kappa,tau,theta
     f = psi**4/12+xxh*psi**3/3+psi**2+(2*xxh-2*deltas)*psi-&
               2*deltas*xxh+xxh*xxh
     derivf = psi**3/3+xxh*psi**2+2*psi+2*xxh-2*deltas
     end subroutine

      end module FieldQuantclass

max_stars_repo_name: hhslepicka/IMPACT-Z

max_issues_repo_name: hhslepicka/IMPACT-Z

max_forks_repo_name: hhslepicka/IMPACT-Z


text:
module timing_module
    implicit none

    integer, dimension(8), private :: dt
    real, private :: h, m, s, ms, tt
    real, private :: last_tt

contains
    subroutine start_timing
        implicit none

        call date_and_time(VALUES=dt)
        print 100, dt(1:3), dt(5:8)
        h = real(dt(5))
        m = real(dt(6))
        s = real(dt(7))
        ms = real(dt(8))
        last_tt = 60*(60*h+m)+s+ms/1000.

    100 format (1x, i4, '/', i2, '/', i2, 1x, i2, &
        ':', i2, ':', i2, 1x, i3)
    end subroutine start_timing

    subroutine end_timing
        implicit none

        call date_and_time(VALUES=dt)
        print 100, dt(1:3), dt(5:8)

    100 format (1x, i4, '/', i2, '/', i2, 1x, i2, &
        ':', i2, ':', i2, 1x, i3)
    end subroutine end_timing

    real function time_difference()
        implicit none

        tt = 0.
        call date_and_time(VALUES=dt)
        h = real(dt(5))
        m = real(dt(6))
        s = real(dt(7))
        ms = real(dt(8))
        tt = 60*(60*h+m)+s+ms/1000.
        time_difference = tt-last_tt
        last_tt = tt
    end function time_difference
end module

max_stars_repo_name: hechengda/fortran

max_issues_repo_name: hechengda/fortran

max_forks_repo_name: hechengda/fortran