C  5 JUL 05 - MWS - SELECT NEW ATOM,BASIS,EFP,PCM,DAF DIMENSIONS
C 12 NOV 98 - GDF - CHANGE BIT PACKING TO ISHIFT
C 28 SEP 97 - MWS - CONVERT PARALLEL CALLS FROM TCGMSG TO DDI
C 24 JUN 97 - MWS - EDAVID: COPY ALL IMPROVED VECTORS BEFORE EDTCHK
C 14 FEB 97 - MWS - EDAVID: FLUSH PRINT BUFFER EVERY ITERATION
C 17 OCT 96 - MWS - EDAVID: FORCE THREE CI ITERATIONS, NO MATTER WHAT.
C 29 SEP 96 - KRG - ADDED SOME ERROR MESSAGES BEFORE ABRT CALLS
C 26 JUL 95 - FPR - EDTCHK: USE CVGETH TOLERANCE IN LOOP
C 10 NOV 94 - MWS - REMOVE FTNCHEK WARNINGS
C 13 DEC 93 - TLW - DO CI DIAGONALIZATION IN PARALLEL
C 16 JUL 93 - MWS - USE SEPARATE STORAGE FOR I/O BUFFERS
C 11 JAN 92 - TLW - MAKE READS PARALLEL
C 10 JAN 92 - TLW - CHANGE REWINDS TO CALL SEQREW
C  6 JAN 92 - TLW - MAKE WRITES PARALLEL;ADD COMMON PAR
C  9 NOV 91 - MWS - EDNEWB,EDTRUN,FRARBR,FRTRBT: CI VECTOR FILE HEADER
C 26 JUN 91 - MWS - EDAVID: TRY TO FORCE AT LEAST 3 DAVIDSON ITERS
C 30 MAR 90 - MWS - DIVIDED GUGDG INTO GUGDGA AND GUGDGB
C
C*MODULE GUGDGB  *DECK EDAVID
      SUBROUTINE EDAVID(LUPRNT,LUMATA,LUVECS,NLARGE,IDLORT,IDHIRT,
     *                  MXBASV,MXMPRV,NEXTRA,MXITER,KPRINT,LENA,LENBAB,
     *                  MXVECT,IBASIS,CVGETH,MODE, ESHIFT,ADIAGS,EIGVAL,
     *                  XVECT,BVECT,AXB,AMATRX,BTAB,CVECT,WRK,
     *                  IWRK,IMCONV,MPRVME,XX1,XX2,IXX,NEMEMX,INDX)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL CVGART,PRINT,PRINT0,PRINT1,PRINT2,PRINT3,PRINT4,PRINTE
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION BVECT(NLARGE,*),AXB(NLARGE,*),XVECT(NLARGE,MXVECT),
     *          CVECT(MXBASV,MXVECT,2),BTAB(LENBAB,2),EIGVAL(MXBASV),
     *          WRK(MXBASV,8),IWRK(MXBASV),IMCONV(MXVECT),
     *          MPRVME(MXMPRV),AMATRX(LENA),ADIAGS(NLARGE),INDX(NLARGE),
     *          XX1(NEMEMX),XX2(NEMEMX),IXX(NEMEMX)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
  001 FORMAT(17H CONVERGENCE TOLS,1P,5(I8,E10.2)/(17X,5(I8,E10.2)))
  002 FORMAT(1X,I3,I4,F12.8,I3,5(F17.9,I3)/(23X,5(F17.9,I3)))
  003 FORMAT(1X,'AFTER',I5,' ITERATIONS THE EIGENVALUES ARE -'/
     *       (5(3H  (,I2,1H),F20.12)))
  980 FORMAT(1X,'****** DAVIDSON METHOD NOT CONVERGED AFTER',I4
     *      ,' ITERATIONS ******'/)
  981 FORMAT(1X,'****** ERROR IN EDAVID:  DIAGONALIZATION OF SUBSPACE'
     *      ,' FAILED ON ITERATION',I5,' WITH ERROR VALUE',I5)
  982 FORMAT(1X,'****** ERROR IN EDAVID:  INCONSISTENT DIMENSIONS'/
     *       '        NCVECT=',I7,' MXVECT=',I7)
  983 FORMAT(1X,'****** ERROR IN EDAVID:  FINAL VECTORS ARE LINEARLY'
     *      ,' DEPENDENT  NDROP=',I5)
C
C* 23 FEB 1984 - STE
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE) DATE: NOV 1982
C*
C*    PURPOSE -
C*       FIND FIRST FEW EIGENVALUES AND VECTORS USING DAVIDSON'S METHOD
C*                      *                             *****
C*
C*    METHOD -
C*       SOLVE  A*X = E*X  FOR ROOTS IDLORT TO IDHIRT (ASCENDING ORDER)
C*       USING THE DAVIDSON METHOD (J. COMP. PHYS. 17, PP. 87-94, 1975).
C*       THESE ROUTINES FOLLOW THE NOTATION USED BY DAVIDSON IN "MATRIX
C*       EIGENVECTOR METHODS" PRESENTED AT THE NATO ADVANCED STUDY INST-
C*       ITUTE ON METHODS IN COMPUTATIONAL MOLECULAR PHYSICS HELD IN
C*       BAD WINDSHEIM, GERMANY, IN AUGUST 1982.
C*       THE METHOD IS AS FOLLOWS:
C*       THE APPROXIMATE EIGENVECTORS  X  ARE GIVEN BY
C*                 X = B*C
C*       WHERE  B  IS A SET OF BASIS VECTORS WHICH INCREASE IN NUMBER
C*       BY NIMPRV VECTORS EACH ITERATION AND  C  ARE EIGENVECTORS
C*       OF B(TRANSPOSE)*A*B.  THE EIGENVALUES ASSOCIATED WITH  C
C*       ARE THE APPROXIMATE EIGENVALUES.  IF THE RESIDUAL VECTOR IS
C*          R = (A - E)*X  =  (A*B)*C - (B)*C*E,
C*       THEN THE NEXT NEWTON-RAPHSON BASIS VECTOR WOULD BE R/(E-A).
C*       DAVIDSON USES FIRST ORDER PERTURBATION THEORY TO OBTAIN
C*       THE CORRECTION VECTOR
C*                 Q = R / (E - A(I,I)).
C*       WHEN THE CORRECTION VECTORS ARE ORTHONORMAL TO THE PREVIOUS
C*       BASIS VECTORS AND THEMSELVES, THEY BECOME THE NEXT BASIS
C*       VECTORS.
C*       WHEN THE NUMBER OF BASIS VECTORS BECOMES TOO LARGE (OR ONE
C*       OF SEVERAL ROOTS CONVERGES), THE PROCESS IS RESTARTED USING
C*       NEW APPROXIMATE VECTORS  X  AS THE BASIS VECTORS.
C*       CONVERGENCE IS REACHED ON EACH ROOT WHEN THE SUM OF THE
C*       SQUARES OF THE LAST NIMPRV ELEMENTS OF THE CONTRACTION
C*       VECTOR  C  IS BELOW THE REQUESTED THRESHOLD.
C*
C*    LIMITATIONS -
C*       THE METHOD WORKS BEST FOR DIAGONALLY DOMINANT MATRICES WHERE
C*       A GOOD GUESS TO THE EIGENVECTOR IS A UNIT VECTOR.  IT IS
C*       POSSIBLE TO "CONVERGE" TO THE WRONG ROOT WHEN THE APPROXIMATE
C*       VECTOR IS ACCIDENTALY CLOSE TO A HIGHER ROOT, ALLOWING THE
C*       RESIDUAL VECTOR TO BECOME SMALL ENOUGH THAT A WEAK CONVERGENCE
C*       THRESHOLD WILL BE SATISFIED.
C*       A STRONG CONVERGENCE REQUIREMENT WILL USUALLY FORCE A MOVE
C*       TO THE CORRECT ROOT AS THE ORTHOGONALIZED "NOISE" OF THE SMALL
C*       RESIDUAL FINDS A LOWER ROOT.
C*       IF THE MAXIMUM SIZE OF THE EXPANSION BASIS IS "TOO SMALL,"
C*       AN INSUFFICIENT PORTION OF THE ORIGINAL SPACE MAY BE SPANNED,
C*       RESULTING IN TOLERANCES THAT ARE LESS THAN THE TRUE ACCURACY.
C*       IF A SUFFICIENTLY LARGE EXPANSION BASIS IS NOT FEASIBLE, THE
C*       CONVERGENCE THRESHOLD SHOULD BE SUFFICIENTLY SMALL TO ENSURE
C*       THE DESIRED ACCURACY.
C*
C*    SPECIAL FEATURES -
C*       THE NUMBER OF ROOTS EACH ITERATION IS UNDER USER CONTROL,
C*       WHEREAS THE ORIGINAL DAVIDSON METHOD IMPROVED JUST ONE ROOT
C*       AND THE LIU MODIFICATION IMPROVED ALL ROOTS (EVEN AFTER SOME
C*       HAD CONVERGED).
C*
C*    ON ENTRY -
C*       LUPRNT - INTEGER (LOGICAL UNIT NO.)
C*                PRINT FILE
C*                IF POSITIVE, PRINT INFORMATION REQUESTED BY KPRINT
C*                IF NEGATIVE, PRINT ONLY ERROR MESSAGES
C*                IF ZERO, SUPPRESS ALL PRINTING, INCLUDING ERROR MESS.
C*       LUMATA - INTEGER (LOGICAL UNIT NO.)
C*                INPUT MATRIX  A  FILE
C*       LUVECS - INTEGER (LOGICAL UNIT NO.)
C*                IF ZERO, SUFFICIENT MEMORY IS AVAILABLE TO HAVE
C*                         BVECT  AND  AXB  MEMORY RESIDENT.
C*                OTHERWISE  BVECT  WILL BE STORED ON UNIT LUVECS
C*                        AND  AXB  WILL BE STORED ON UNIT LUVECS+3
C*                IN THIS CASE THE ARRAYS BVECT AND AXB MAY BE
C*                EQUIVALENCED BY CALL, BUT BVECT(1,MXVECT+1) MUST BE
C*                EQUIVALENCED TO XVECT(1,1).
C*       NLARGE - INTEGER
C*                DIMENSION OF REAL SYMMETRIC INPUT MATRIX  A.
C*       IDLORT - INTEGER
C*                NO. FOR LOWEST ROOT TO BE OBTAINED (LOWEST IS 1)
C*       IDHIRT - INTEGER
C*                NO. FOR HIGHEST ROOT TO BE OBTAINED (HIGHEST = NLARGE)
C*       MXBASV - INTEGER
C*                MAXIMUM NUMBER OF BASIS VECTORS (MXBASV > IDHIRT)
C*                THE NO. OF BASIS VECTORS IS TRUNCATED TO IDHIRT WHEN
C*                THERE ARE MXBASV OF THEM.
C*       MXMPRV - INTEGER
C*                MAXIMUM NO. OF VECTORS TO SIMULTANEOUSLY IMPROVE
C*       NEXTRA - INTEGER
C*                THIS MANY EXTRA EXPANSION BASIS VECTORS WILL BE
C*                USED ON THE FIRST ITERATION.  THE NUMBER OF EXTRA
C*                VECTORS IS DECREASED BY ONE EACH ITERATION UNTIL
C*                NO EXTRA VECTORS ARE PRESENT.  THIS IS USEFUL IN
C*                "CAPTURING" VECTORS FOR HIGHER ROOTS.
C*       MXITER - INTEGER
C*                MAXIMUM NUMBER OF ITERATIONS BEFORE QUITTING
C*       KPRINT - INTEGER
C*                PRINTED OUTPUT CONTROL MASK (NO PRINTING IF ZERO)
C*                BIT  0  FINAL EIGENVALUES
C*                BIT  1  FINAL TOLERANCES
C*                BIT  2  EIGENVALUES AND TOLERANCES AT EACH TRUNCATION
C*                BIT  3  EIGENVALUES EVERY ITERATION
C*                BIT  4  TOLERANCES EVERY ITERATION
C*                E.G., KPRINT=10 PRINTS EIGENVALUES EVERY ITERATION
C*                AND FINAL TOLERANCES.
C*       LENA   - INTEGER
C*                LENGTH OF BUFFER USED TO PROCESS MATRIX  A.
C*       LENBAB - INTEGER
C*                WORDS OF MEMORY AVAILABLE TO STORE TRANSFORMED  A.
C*                = MXBASV*(MXBASV+1)/2
C*       MXVECT - INTEGER
C*                MAXIMUM NO. OF VECTORS THAT MAY SAFELY BE MEMORY
C*                RESIDENT (ALLOWS FOR NEXTRA VECTORS).
C*                = MAX(IDHIRT, MXMPRV+NEXTRA)
C*       IBASIS - INTEGER
C*                INITIAL NO. OF BASIS VECTORS IN  BVECT.
C*       CVGETH - W.P. REAL
C*                RELATIVE EIGENVALUE CONVERGENCE THRESHOLD
C*                (VECTOR ACCURACY IS SQUARE ROOT OF THIS VALUE)
C*       MODE   - INTEGER
C*                DETERMINES WHICH ROUTINE IS CALLED TO DO A * B
C*                = 1,  A  IS MEMORY RESIDENT IN  AMATRX
C*                = 3,  A  IS ON LUMATA IN GAMESS/GUGA FORMAT
C*       ESHIFT - W.P. REAL
C*                AMOUNT DIAGONAL ELEMENTS HAVE BEEN SHIFTED BY.
C*                (MUST BE ADDED TO EIGVAL TO OBTAIN TRUE EIGENVALUES)
C*       ADIAGS - W.P. REAL (NLARGE)
C*                SHIFTED DIAGONAL ELEMENTS OF MATRIX A
C*       INDX   - INTEGER (NLARGE)
C*                INDEXING ARRAY FOR GAMESS/GUGA OFF-DIAGONAL ELEMENTS
C*       AMATRX - W.P. REAL (LENA)  (WHEN MODE = 1)
C*                REAL SYMMETRIC MATRIX  A
C*       BVECT  - W.P. REAL (NLARGE,MXBASV)
C*                EXPANSION BASIS VECTORS (B)
C*                IF LUVECS GT 0, DIMENSION IS (NLARGE,2*MXVECT) AND
C*                THE SECOND SET OF VECTORS MAY BE EQUIVALENCED TO XVECT
C*       AXB    - W.P. REAL (NLARGE,MXBASV)
C*                A * B VECTORS
C*                IF LUVECS GT 0, DIMENSION IS (NLARGE,MXVECT)
C*                AND MAY BE EQUIVALENCED TO BVECT.
C*       AMATRX - W.P. REAL (LENA)
C*                INPUT BUFFER FOR REAL SYMMETRIC MATRIX  A
C*       BTAB   - W.P. REAL (LENBAB,2)
C*                EXPANSION BASIS REPRESENTATION OF  A  (TWO COPIES)
C*                (B(TRANSPOSE) * A * B)
C*       CVECT  - W.P. REAL (MXBASV,MXVECT,2)
C*                CONTRACTION VECTORS (EIGENVECTORS OF BTAB) AND
C*                CONTRACTION VECTORS FROM PREVIOUS ITERATION.
C*
C*    ON EXIT -
C*       MXITER - INTEGER
C*                NEGATIVE IF METHOD HAS NOT CONVERGED IN * ITERATIONS
C*       EIGVAL - W.P. REAL (MXBASV)
C*                EIGENVALUES + ESHIFT ARE IN POSITIONS IDLORT TO IDHIRT
C*       XVECT  - W.P. REAL (NLARGE,IDHIRT)
C*                REQUESTED EIGENVECTORS OF  A. IF IDLORT IS NOT EQUAL
C*                TO ONE, THE VECTORS 1 TO IDLORT-1 IN XVECT ARE NOT
C*                ACCURATE.
C*
C*    WORK AREAS -
C*       XVECT  - W.P. REAL (NLARGE, MXVECT)
C*                (EQUIVALENCED TO BVECT(1,MXVECT+1) WHEN LUVECS GT 0)
C*       WRK    - W.P. REAL (MXBASV,8)
C*                WRK(*,*) = WORK SPACE TO DIAGONALIZE BTAB.
C*                WRK(*,1) = EIGENVALUES + ESHIFT
C*                WRK(*,2) = CONVERGENCE TOLERANCES
C*                WRK(*,3) = LENGTH**2 OF UNNORMALIZED RESIDUAL VECTORS
C*       IWRK   - INTEGER (MXBASV)
C*                WORK SPACE TO DIAGONALIZE BTAB.
C*       IMCONV - INTEGER (MXVECT)
C*                ROOT CONVERGED = 1, NOT CONVERGED = 0
C*       MPRVME - INTEGER (MXMPRV+NEXTRA)
C*                LIST OF ROOTS TO BE SIMULTANEOUSLY IMPROVED
C*
C           SET UP PRINTING VARIABLES
C
      LUEMSG=ABS(LUPRNT)
      PRINT  = LUPRNT.GT.0 .AND. MASWRK
      PRINTE = LUEMSG.NE.0 .AND. MASWRK
      PRINT0 = PRINT .AND. LOGAND(KPRINT,2**0).NE.0
      PRINT1 = PRINT .AND. LOGAND(KPRINT,2**1).NE.0
      PRINT2 = PRINT .AND. LOGAND(KPRINT,2**2).NE.0
      PRINT3 = PRINT .AND. LOGAND(KPRINT,2**3).NE.0
      PRINT4 = PRINT .AND. LOGAND(KPRINT,2**4).NE.0
C
C           INITIALIZE CONVERGENCE AND ROOT SELECTION ARRAYS
C
      NBASIS=IBASIS
      NIMPRV=0
      MAXIMP=MXMPRV+NEXTRA
      NCVECT=MAX(IDHIRT,MAXIMP)
      IF(NCVECT.NE.MXVECT) THEN
         IF(PRINTE) WRITE(LUEMSG,982) NCVECT,MXVECT
         CALL ABRT
      END IF
      DO 110 K=1,MXVECT
         IMCONV(K)=1
         IF(K.LT.IDLORT) GO TO 110
            IMCONV(K)=0
            IF(NIMPRV.EQ.MAXIMP) GO TO 110
               NIMPRV=NIMPRV+1
               MPRVME(NIMPRV)=K
  110 CONTINUE
      MAXBAS=MIN(MXBASV,NLARGE)
      NHITRI=(IDHIRT*IDHIRT+IDHIRT)/2
      IJBTAB=(NBASIS*NBASIS+NBASIS)/2
      NXTRRT=NEXTRA
      VERMAX = CVGETH
      IDMAX = 0
C
C      ********* BEGIN LOOP 190 OVER DIAGONALIZATION ITERATIONS *******
C
      CVGESV = CVGETH
      VERMSV = VERMAX
      DO 190 ITER=1,MXITER
C
C             MAKE SURE WE ALWAYS DO AT LEAST THREE CI ITERATIONS
C
         IF(ITER.EQ.1) THEN
            CVGETH = 0.000001D+00*CVGETH
            VERMAX = 0.000001D+00*VERMAX
         END IF
         IF(ITER.EQ.4) THEN
            CVGETH = CVGESV
            VERMAX = VERMSV
         END IF
C
         CVGART=.FALSE.
         NEWB=NBASIS+1
         IF(LUVECS.GT.0) NEWB=MXVECT+1
C
C              FIND NIMPRV NEW VECTORS TO ADD TO EXPANSION BASIS.
C
         CALL EDNEWB(LUEMSG,LUVECS,NLARGE,IDHIRT,NBASIS,MAXBAS,CVGETH,
     *               EIGVAL,CVECT,AXB,ADIAGS,NIMPRV,MPRVME,BVECT,
     *               BVECT(1,NEWB),WRK(1,3),IMCONV,CVGART,ITER)
         IF(ITER.GT.3  .AND.  NIMPRV.EQ.0) GO TO 200
C
C              QUIT IF ALL CORRECTION VECTORS ARE LINEARLY DEPENDENT
C
         NEWIMP=NIMPRV
         NBAS0=NBASIS
         NEWVEC=NBASIS+1
         IF(LUVECS.GT.0) NEWVEC=1
         NBASIS=NBASIS+NIMPRV
C
C           MULTIPLY MATRIX  A  TIMES NEW BASIS VECTORS
C                 (THIS IS MOST OF THE WORK)
C
         IF(MODE.EQ.1) CALL EDPAB1(NLARGE,NIMPRV,LENA,AMATRX,ADIAGS,
     *                             BVECT(1,NEWB),AXB(1,NEWVEC))
         IF(MODE.EQ.3) CALL EDPAB3(LUMATA,NLARGE,NIMPRV,ADIAGS,
     *                             BVECT(1,NEWB),AXB(1,NEWVEC),
     *                             INDX,XX1,XX2,IXX,NEMEMX)
         NEWB=1
         IF(LUVECS.LE.0) GO TO 130
            NEWB=MXVECT+1
            LUAXB = LUVECS + 3
            DO 120 I=1,NIMPRV
               CALL SQWRIT(LUAXB,AXB(1,I),NLARGE)
  120       CONTINUE
  130    CONTINUE
C
C           COMPLETE TRANSFORMATION TO EXPANSION BASIS SPACE
C
         CALL FRTRBT(LUVECS,NLARGE,NBAS0+1,NBASIS,LENBAB,
     *               BVECT(1,NEWB),AXB,IJBTAB,BTAB)
C
C          FIND EIGENVALUES AND CONTRACTION VECTORS
C
         NXTRRT=NXTRRT-1
         IF(NXTRRT.LT.0) NXTRRT=0
         MAXIMP=MXMPRV+NXTRRT
         NCVECT=MAX(IDHIRT,MAXIMP)
         CALL DCOPY(NBAS0*NCVECT,CVECT,1,CVECT(1,1,2),1)
         CALL DCOPY(IJBTAB,BTAB,1,BTAB(1,2),1)
         IERR = 0
         CALL GLDIAG(NBASIS,NCVECT,NBASIS,BTAB(1,2),WRK,EIGVAL,CVECT
     *              ,IERR,IWRK)
         IF(IERR.NE.0) THEN
            IF(PRINTE) WRITE(LUEMSG,981) ITER,IERR
            CALL ABRT
         END IF
C
C           CHECK FOR CONVERGENCE
C
         CALL EDTCHK(NBASIS,NBAS0,IDLORT,NCVECT,MAXIMP,CVGETH,
     *              NIMPRV,MPRVME,CVECT,CVECT(1,1,2),IMCONV,WRK(1,2),
     *              CVGART,VERMAX,IDMAX,ITER)
         IF(PRINT3) THEN
            DO 140 K=IDLORT,IDHIRT
               WRK(K,1)=EIGVAL(K)+ESHIFT
  140       CONTINUE
            WRITE(LUPRNT,002)
     *            ITER,NEWIMP,VERMAX,IDMAX,(WRK(K,1),K,K=IDLORT,IDHIRT)
            CALL FLSHBF(LUPRNT)
         END IF
C
C           QUIT IF NO MORE VECTORS TO IMPROVE
C
         IF(ITER.GT.3  .AND.  NIMPRV.EQ.0) GO TO 200
C
         IF(PRINT4) WRITE(LUPRNT,001)
     *                       (IMCONV(K),WRK(K,2),K=IDLORT,IDHIRT)
         IF(NBASIS+NIMPRV.LE.MAXBAS .AND. .NOT.CVGART) GO TO 190
C
C              NO. OF BASIS VECTORS HAS REACHED THE MAXIMUM,
C              OR A ROOT HAS CONVERGED (BUT MORE TO DO)
C                   RECALCULATE (TRUNCATE) EXPANSION BASIS
C
            CALL EDTRUN(LUVECS,NLARGE,NBASIS,IDHIRT,NHITRI,IJBTAB,
     *                  CVECT,BVECT,AXB,BTAB,XVECT)
            NBASIS=IDHIRT
            IF(PRINT2) THEN
               IF(.NOT.PRINT3) WRITE(LUPRNT,002)
     *            ITER,NIMPRV,VERMAX,IDMAX,(WRK(K,1),K,K=IDLORT,IDHIRT)
               IF (MASWRK) WRITE(LUPRNT,001)
     *            (IMCONV(K),WRK(K,2),K=IDLORT,IDHIRT)
            END IF
  190 CONTINUE
      ITER = MXITER
C
C           NO CONVERGENCE WITHIN ITERATION LIMIT
C
      IF (PRINTE) THEN
         WRITE(LUEMSG,980) ITER
         WRITE(LUEMSG,001) (IMCONV(K),WRK(K,2),K=IDLORT,IDHIRT)
         WRITE(LUEMSG,003) ITER,(K,EIGVAL(K),K=IDLORT,IDHIRT)
      END IF
      MXITER=-MXITER
C
C           CONVERGENCE FOR ALL DESIRED ROOTS HAS BEEN ACHIEVED
C           MAKE SURE THE CONTRACTION VECTORS ARE ORTHONORMAL
C
  200 CONTINUE
      NDROP = 0
      CALL ONVMGS(NBASIS,IDHIRT,CVECT,NDROP)
      IF(NDROP.NE.0) THEN
         IF(PRINTE) WRITE(LUEMSG,983) NDROP
         CALL ABRT
      END IF
C
C           GENERATE FINAL EIGENVECTOR(S)
C
      CALL FRARBR(LUVECS,NLARGE,NBASIS,IDHIRT,BVECT,CVECT,XVECT)
C
C           GENERATE FINAL EIGENVALUES
C
      DO 210 K=IDLORT,IDHIRT
         EIGVAL(K)=EIGVAL(K)+ESHIFT
  210 CONTINUE
C
C           PRINT FINAL RESULTS
C
      IF(PRINT1) WRITE(LUPRNT,001)
     *                    (IMCONV(K),WRK(K,2),K=IDLORT,IDHIRT)
      IF(PRINT0) WRITE(LUPRNT,003) ITER,(K,EIGVAL(K),K=IDLORT,IDHIRT)
      RETURN
      END
C*MODULE GUGDGB  *DECK EDNEWB
C*  7 FEB 88 - STE * 12 FEB 84 - STE
C*
C*    AUTHOR:  S. T. ELBERT (AMES LABORATORY - USDOE)  DATE: JAN 1983
C*
C*    PURPOSE -
C*       WHEN FINDING EIGENVECTOR(S) USING THE DAVIDSON METHOD, THIS
C*                    *                        *
C*       ROUTINE FINDS THE NEW SET OF EXPANSION VECTORS  B
C*                         ***                           *
C*
C*    ON ENTRY -
C*       LEUMSG - INTEGER (LOGICAL UNIT NUMBER)
C*                LOGICAL UNIT FOR PRINTING ERROR MESSAGES (IF NON-ZERO)
C*       LUVECS - INTEGER (LOGICAL UNIT NUMBER)
C*                IF NEGATIVE, BVECT AND AXB ARE MEMORY RESIDENT, ELSE
C*                BVECT IS ON LOGICAL UNIT LUVECS AND
C*                AXB   IS ON LOGICAL UNIT LUVECS+1.
C*       NLARGE - INTEGER
C*                THIS IS THE DIMENSION OF THE MATRIX BEING SOLVED.
C*       IDHIRT - INTEGER
C*                THIS IS THE NUMBER OF THE HIGHEST VECTOR SOUGHT.
C*       NBASIS - INTEGER
C*                THIS IS THE CURRENT NUMBER OF EXPANSION VECTORS.
C*       MAXBAS - INTEGER
C*                MAXIMUM NO. OF EXPANSION VECTORS ALLOWED BEFORE THE
C*                EXPANSION SPACE MUST BE UPDATED (TRUNCATED).
C*       CVGETH - W.P. REAL
C*                SQUARE OF CONVERGENCE THRESHOLD VALUE.
C*       EIGVAL - W.P. REAL (IDHIRT)
C*                EIGENVALUES OF B(TRANSPOSE) * A * B
C*                (AND OF MATRIX  A  AT CONVERGENCE)
C*       CVECT  - W.P. REAL (NBASIS,IDHIRT)
C*                EIGENVECTORS OF B(TRANSPOSE) * A * B
C*       AXB    - W.P. REAL (NLARGE,NBASIS)
C*                PRODUCT OF MATRIX  A  AND EXPANSION VECTORS  B (BVECT)
C*       ADIAG  - W.P. REAL (NLARGE)
C*                DIAGONAL ELEMENTS OF MATRIX  A
C*
C*       IMPRV  - INTEGER
C*                NO. OF NEW EXPANSION VECTORS SOUGHT
C*       MPRVME - INTEGER (IMPRV)
C*                ROOTS FOR WHICH NEW EXPANSION VECTORS SHOULD BE FOUND
C*       BVECT  - W.P. REAL (NLARGE,MAXBAS)
C*                OLD EXPANSION VECTORS  B
C*
C*    ON EXIT -
C*       IMPRV  - INTEGER
C*                THE NUMBER OF NEW EXPANSION VECTORS FOUND
C*       MPRVME - INTEGER (IMPRV)
C*                ROOTS FOR WHICH NEW EXPANSION VECTORS WERE FOUND
C*       BVECT  - W.P. REAL (NLARGE,MAXBAS)
C*                NEW EXPANSION VECTORS  B  (IF ROOM)
C*       RQQP   - W.P. REAL (NLARGE,IMPRV)
C*                NEXT EXPANSION VECTOR(S)
C*       QNORM  - W.P. REAL (IMPRV)
C*                MAGNITUDE OF NEXT EXPANSION VECTOR(S) BEFORE
C*                EXPLICIT ORTHOGONALIZATION
C*       IMCONV - INTEGER (IDHIRT)
C*                CONVERGED ROOTS ARE MARKED WITH 1
C*       CVGART - LOGICAL
C*                SET TO TRUE IF THE NEXT EXPANSION VECTOR TURNS OUT
C*                TO BE A NULL VECTOR.
C*
C*
      SUBROUTINE EDNEWB(LUEMSG,LUVECS,NLARGE,IDHIRT,NBASIS,MAXBAS,
     *                  CVGETH,EIGVAL,CVECT,AXB,ADIAG,IMPRV,MPRVME,
     *                  BVECT,RQQP,QNORM,IMCONV,CVGART,ITER)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL DOORTH,CVGART,INMEM,GOPARR,DSKWRK,MASWRK
C
      DIMENSION RQQP(NLARGE,IMPRV),BVECT(NLARGE,*),AXB(NLARGE,*),
     *          ADIAG(NLARGE),CVECT(NBASIS,IDHIRT),EIGVAL(MAXBAS),
     *          QNORM(IMPRV),MPRVME(IMPRV),IMCONV(IDHIRT)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ZERO=0.0D+00, ONE=1.0D+00, SMALL=1.0D-06)
C
C>001 FORMAT(' EDNEWB: ',A/(8(2I3,1P,E10.2)))
C>002 FORMAT(' EDNEWB: ',A/(8(I3, 1P,E10.2)))
C>003 FORMAT(' EDNEWB: ',A/(10I3))
  981 FORMAT(1H0,80(1H*)/' DISASTER HAS STRUCK - UNABLE TO GENERATE'
     *      ,' NEW SET OF ORTHONORMAL EXPANSION BASIS VECTORS'/
     *       ' THE NORMS OF THE FIRST ORDER CORRECTION VECTORS ARE -'/
     *      (1P,10E10.2))
  982 FORMAT(' EDNEWB:  PASS',I2,' THROUGH IRNDER HAS QPNORM =',1P,E10.2
     *      ,' QNORM(',I3,') =',E10.2,' FOR ROOT',I3)
C
      LUAXB = 0
      INMEM = .TRUE.
      IF(LUVECS .GT. 0) THEN
         INMEM = .FALSE.
         LUAXB = LUVECS + 3
      END IF
C
C           LOOP OVER UNCONVERGED VECTORS
C
  100 CONTINUE
      IF(.NOT.INMEM) THEN
         CALL SEQREW(LUVECS)
         READ(LUVECS)
         CALL SEQREW(LUAXB)
      END IF
C
C              GENERATE RESIDUAL VECTOR  R  IN RQQP = (AB)C - (B)CE
C
      CALL VCLR(RQQP,1,NLARGE*IMPRV)
C
      DO 150 K=1,NBASIS
         IF(.NOT.INMEM) THEN
            M=1
            CALL SQREAD(LUAXB,AXB,NLARGE)
            IF(NLARGE.EQ.0) THEN
               IF (MASWRK) WRITE(LUEMSG,*)
     *            'EDNEWB: UNEXPECTED EOF ON FILE LUAXB'
               CALL ABRT
               STOP
            END IF
         ELSE
            M=K
         END IF
         DO 130 L=1,IMPRV
            IDMPRV=MPRVME(L)
            CALL DAXPY(NLARGE,CVECT(K,IDMPRV),AXB(1,M),1,RQQP(1,L),1)
  130    CONTINUE
C
         IF(.NOT.INMEM) THEN
            CALL SQREAD(LUVECS,BVECT,NLARGE)
            IF(NLARGE.EQ.0) THEN
               IF (MASWRK) WRITE(LUEMSG,*)
     *                  'EDNEWB(1): UNEXPECTED EOF ON CI VECTOR FILE'
               CALL ABRT
               STOP
            END IF
         END IF
         DO 140 L=1,IMPRV
            IDMPRV=MPRVME(L)
            DZ=EIGVAL(IDMPRV)
            CALL DAXPY(NLARGE,-DZ*CVECT(K,IDMPRV),BVECT(1,M),1,
     *                 RQQP(1,L),1)
  140    CONTINUE
  150 CONTINUE
C
C           GENERATE FIRST ORDER CORRECTION  Q  IN  RQQP = R/(E-A(I,I))
C
      DO 170 L=1,IMPRV
         IDMPRV=MPRVME(L)
         DZ=EIGVAL(IDMPRV)
         QNORM(L)=ZERO
         DO 160 I=1,NLARGE
            D = DZ - ADIAG(I)
            IF(ABS(D).LT.SMALL)D=-SMALL
            RQQP(I,L)=RQQP(I,L)/D
            QNORM(L)=QNORM(L) + (RQQP(I,L)*RQQP(I,L))
  160    CONTINUE
  170 CONTINUE
C
C>    IF (MASWRK) THEN
C>    WRITE(LUEMSG,001)
C>   *   'CVECT',(((JR,JC,CVECT(JR,JC)),JR=1,NBASIS),JC=1,IDHIRT)
C>    WRITE(LUEMSG,001)
C>   *   'BVECT',(((JR,JC,BVECT(JR,JC)),JR=1,NLARGE),JC=1,NBASIS)
C>    WRITE(LUEMSG,001)
C>   *   'AXB',(((JR,JC,AXB(JR,JC)),JR=1,NLARGE),JC=1,NBASIS)
C>    WRITE(LUEMSG,001)
C>   *   'Q',(((JR,JC,RQQP(JR,JC)),JR=1,NLARGE),JC=1,IMPRV)
C>    WRITE(LUEMSG,002) 'ADIAG',((JJ,ADIAG(JJ)),JJ=1,NLARGE)
C>    WRITE(LUEMSG,002) 'EIGVAL',((JJ,EIGVAL(JJ)),JJ=1,MAXBAS)
C>    WRITE(LUEMSG,002) 'QNORM',((JJ,QNORM(JJ)),JJ=1,IMPRV)
C>    WRITE(LUEMSG,003) 'MPRVME',(MPRVME(J),J=1,IMPRV)
C>    ENDIF
C
C                 EXPLICIT ORTHOGONALIZATION OF  QP  IN  RQQP
C     (LOOP IS REPEATED ONLY IF SUFFICIENT ROUNDING ERROR
C           MAY HAVE OCCURED ON THE PREVIOUS PASS)
C
      DOORTH=.TRUE.
      DO 290 IRNDER=1,5
         IF(.NOT.DOORTH) GO TO 290
            DOORTH=.FALSE.
            IF(.NOT.INMEM) THEN
               CALL SEQREW(LUVECS)
               READ(LUVECS)
            END IF
            DO 230 K=1,NBASIS
               IF(.NOT.INMEM) THEN
                  M=1
                  CALL SQREAD(LUVECS,BVECT,NLARGE)
                  IF(NLARGE.EQ.0) THEN
                     IF (MASWRK) WRITE(LUEMSG,*)
     *                  'EDNEWB(2): UNEXPECTED EOF ON CI VECTOR FILE'
                     CALL ABRT
                     STOP
                  END IF
               ELSE
                  M=K
               END IF
C
C              ORTHOGONALIZE TO PREVIOUS BASIS VECTORS
C
               DO 220 L=1,IMPRV
                  S=DDOT(NLARGE,RQQP(1,L),1,BVECT(1,M),1)
                  CALL DAXPY(NLARGE,-S,BVECT(1,M),1,RQQP(1,L),1)
  220          CONTINUE
  230       CONTINUE
C
            DO 260 L=1,IMPRV
               IDMPRV=MPRVME(L)
               QPNORM=DDOT(NLARGE,RQQP(1,L),1,RQQP(1,L),1)
C
C                 CHECK CONTRIBUTION OF THE NEXT VECTOR TO THE SPACE OF
C                 ROOT IDMPRV. IF STILL SMALL AFTER REORTHOGONALIZATION,
C                 MARK THIS ROOT AS CONVERGED (IMCONV(IDMPRV)=1).
C
               IF(ITER.GT.3) THEN
               IF(QPNORM.LE.CVGETH*QNORM(L)) IMCONV(IDMPRV)=IRNDER-1
               IF(QNORM(L)+QPNORM.EQ.QNORM(L)) IMCONV(IDMPRV)=IRNDER-1
               END IF
C
               IF(IMCONV(IDMPRV).EQ.1) THEN
                  SCALE=ZERO
                  CVGART=.TRUE.
               ELSE IF (QPNORM.NE.ZERO) THEN
                  SCALE=ONE/SQRT(QPNORM)
               ELSE
                  SCALE=ONE
               END IF
C
C                 IF NORM SHRANK TOO MUCH, DO IT AGAIN
C
               IF(QPNORM.LT.0.1D+00*QNORM(L)) THEN
                  DOORTH=.TRUE.
                  IF(IRNDER .GE. 3 .AND. LUEMSG .GT. 0 .AND. MASWRK)
     *               WRITE(LUEMSG,982) IRNDER,QPNORM,L,QNORM(L),IDMPRV
               END IF
               QNORM(L) = QPNORM
C
C                 NORMALIZE ORTHOGONAL CORRECTION VECTOR QP
C                 THIS BECOMES THE NEW VECTOR FOR THE NEXT ITERATION
C
               CALL DSCAL(NLARGE,SCALE,RQQP(1,L),1)
               IF(SCALE.EQ.ZERO) GO TO 260
C
C                 ORTHOGONALIZE NEXT SET OF VECTORS TO EACH OTHER
C
                  DO 250 M = L+1,IMPRV
                     S=DDOT(NLARGE,RQQP(1,L),1,RQQP(1,M),1)
                     CALL DAXPY(NLARGE,-S,RQQP(1,L),1,RQQP(1,M),1)
  250             CONTINUE
  260       CONTINUE
  290 CONTINUE
C
      IF (DOORTH) THEN
         DO 295 L = 1,IMPRV
            IF (QNORM(L) .LT. 0.9D+00  .OR.  QNORM(L) .GT. 1.1D+00) THEN
C
C       IF PROGRAM GETS HERE, THERE IS PROBABLY A MEMORY OVERLAP PROBLEM
C
               IF (LUEMSG .GT. 0 .AND. MASWRK) WRITE(LUEMSG,981) QNORM
               CALL ABRT
            END IF
  295    CONTINUE
      END IF
C
C           PUT NEW IMPROVEMENT VECTORS IN  BVECT  AND COUNT THEM
C
      NIMPRV=0
      DO 320 L=1,IMPRV
         IDMPRV=MPRVME(L)
         IF(IMCONV(IDMPRV).EQ.1) GO TO 320
            NIMPRV=NIMPRV+1
            MPRVME(NIMPRV)=IDMPRV
            M=NBASIS+NIMPRV
            IF(.NOT.INMEM) M=L
            CALL DCOPY(NLARGE,RQQP(1,L),1,BVECT(1,M),1)
            IF(.NOT.INMEM) CALL SQWRIT(LUVECS,BVECT(1,M),NLARGE)
  320 CONTINUE
      IF(NIMPRV.GT.0) THEN
         IMPRV=NIMPRV
         RETURN
      END IF
C
C           NO NEW BASIS VECTOR WAS FOUND, ANY MORE UNCONVERGED ROOTS?
C
      DO 340 I=1,IDHIRT
         IF(IMCONV(I).EQ.1) GO TO 340
            NIMPRV=NIMPRV+1
            MPRVME(NIMPRV)=I
            IF(NIMPRV.EQ.IMPRV) GO TO 100
  340 CONTINUE
      IMPRV=NIMPRV
      IF(NIMPRV.NE.0) GO TO 100
      RETURN
      END
C*MODULE GUGDGB  *DECK EDPAB1
C*  2 APR 1983 - STE
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)
C*
C*    PURPOSE -
C*       WHEN FINDING EIGENVECTOR(S) USING THE DAVIDSON METHOD, THIS
C*                    *                        *
C*       ROUTINE FORMS THE MATRIX PRODUCT OF  A  AND  B
C*       FOR CASE 1 (A IS MEMORY RESIDENT).
C*                *
C*
C*    ON ENTRY -
C*       NLARGE - INTEGER
C*                DIMENSION OF INPUT MATRIX A.
C*       NIMPRV - INTEGER
C*                NUMBER OF BASIS VECTORS TO BE MULTIPLIED.
C*       LENA   - INTEGER
C*                NO. OF VALUES OF  A  THAT CAN BE HELD IN MEMORY.
C*       AMATRX - W.P. REAL ( NLARGE*(NLARGE+1)/2 )
C*                INPUT MATRIX A
C*       ADIAG  - W.P. REAL (NLARGE)
C*                VECTOR CONTAINING THE DIAGONAL ELEMENTS OF AMATRX
C*       BVECT - W.P. REAL (NLARGE,NIMPRV)
C*                IMPROVEMENT BASIS VECTOR(S)
C*
C*    ON EXIT -
C*       AXB    - W.P. REAL (NLARGE,NIMPRV)
C*                PRODUCT VECTOR(S).
C*
      SUBROUTINE EDPAB1(NLARGE,NIMPRV,LENA,AMATRX,ADIAG,BVECT,AXB)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      DIMENSION BVECT(NLARGE,NIMPRV),AXB(NLARGE,NIMPRV),ADIAG(NLARGE),
     *          AMATRX(LENA)
C
C        INITIALIZE  AXB(1,*)
C
      A11=ADIAG(1)
      DO 100 NV=1,NIMPRV
         AXB(1,NV)=BVECT(1,NV)*A11
  100 CONTINUE
C
C           FORM AXB = A * B
C
      NIJ=1
      DO 200 NI=2,NLARGE
         ANN=ADIAG(NI)
         DO 105 NV=1,NIMPRV
            AXB(NI,NV)=BVECT(NI,NV)*ANN
  105    CONTINUE
         NIM1=NI-1
         M=MOD(NIM1,3)+1
         GO TO (130,120,110),M
C
  110       CONTINUE
            DO 115 NV=1,NIMPRV
               AXB(NI,NV)=AXB(NI,NV) + BVECT( 2,NV) * AMATRX(NIJ+2)
     *                               + BVECT( 1,NV) * AMATRX(NIJ+1)
  115       CONTINUE
            DO 117 NV=1,NIMPRV
               AXB( 1,NV)=AXB( 1,NV) + BVECT(NI,NV) * AMATRX(NIJ+1)
  117       CONTINUE
            DO 118 NV=1,NIMPRV
               AXB( 2,NV)=AXB( 2,NV) + BVECT(NI,NV) * AMATRX(NIJ+2)
  118       CONTINUE
            GO TO 130
C
  120       CONTINUE
            DO 125 NV=1,NIMPRV
               AXB(NI,NV)=AXB(NI,NV) + BVECT( 1,NV) * AMATRX(NIJ+1)
  125    CONTINUE
            DO 127 NV=1,NIMPRV
               AXB( 1,NV)=AXB( 1,NV) + BVECT(NI,NV) * AMATRX(NIJ+1)
  127    CONTINUE
C
  130    CONTINUE
         NIJ=NIJ+M-1
         IF(NI.LE.3) GO TO 190
            NIJV=NIJ
            DO 180 NV=1,NIMPRV
               BNI=BVECT(NI,NV)
               AXBIV=AXB(NI,NV)
               NIJV=NIJ
               DO 170 NJ=M,NIM1,3
                  AXBIV=AXBIV     + BVECT(NJ  ,NV) * AMATRX(NIJV+1)
     *                            + BVECT(NJ+1,NV) * AMATRX(NIJV+2)
     *                            + BVECT(NJ+2,NV) * AMATRX(NIJV+3)
                  AXB(NJ  ,NV)= AXB(NJ  ,NV) + BNI * AMATRX(NIJV+1)
                  AXB(NJ+1,NV)= AXB(NJ+1,NV) + BNI * AMATRX(NIJV+2)
                  AXB(NJ+2,NV)= AXB(NJ+2,NV) + BNI * AMATRX(NIJV+3)
                  NIJV=NIJV+3
  170          CONTINUE
               AXB(NI,NV)=AXBIV
  180       CONTINUE
            NIJ=NIJV
  190    CONTINUE
         NIJ=NIJ+1
  200 CONTINUE
      RETURN
      END
C*MODULE GUGDGB  *DECK EDPAB3
      SUBROUTINE EDPAB3(LUMATA,NLARGE,NIMPRV,ADIAGS,BVECT,AXB,
     *                  INDX,XX1,IX,IXX,NEMEMX)
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION ADIAGS(NLARGE),BVECT(*),AXB(*),INDX(NLARGE),
     *          XX1(NEMEMX),IX(NEMEMX),IXX(NEMEMX)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      PARAMETER (ONE=1.0D+00)
C
C     ----- CALCULATE AXB = A * BVECT -----
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)  DATE: MAR 1983
C*
C*    PURPOSE -
C*       WHEN FINDING EIGENVECTOR(S) USING THE DAVIDSON METHOD, THIS
C*                    *                        *
C*       ROUTINE FORMS THE MATRIX PRODUCT OF  A  AND  B
C*                                *           *       *
C*       FOR CASE 3 (GAMESS/GUGA FILE)
C*                *
C*
C*    ON ENTRY -
C*       LUMATA - INTEGER (LOGICAL UNIT NO.)
C*                ENERGY MATRIX FILE
C*       NLARGE - INTEGER
C*                DIMENSION OF INPUT MATRIX A.
C*       NIMPRV - INTEGER
C*                NUMBER OF BASIS VECTORS TO BE MULTIPLIED.
C*       NEMEMX - INTEGER
C*                LENGHT OF  A  MATRIX BUFFER.
C*       ADIAGS - W.P. REAL (NLARGE)
C*                VECTOR CONTAINING THE DIAGONAL ELEMENTS OF AMATRX
C*       BVECT  - W.P. REAL (NLARGE,NIMPRV)
C*                IMPROVEMENT BASIS VECTOR(S)
C*    ON EXIT -
C*       AXB    - W.P. REAL (NLARGE,NIMPRV)
C*                PRODUCT VECTOR(S).
C
      MX = 0
*I64  MASK32 = 2**32 - 1
      IF (NIMPRV .GT. 1) GO TO 400
C
C     ----- GET DIAGONAL CI-ELEMENTS CONTRIBUTION -----
C
      IF (GOPARR) CALL DSCAL(NLARGE,ONE/NPROC,ADIAGS,1)
      DO 110 IWKS = 1,NLARGE
         AXB(IWKS) = ADIAGS(IWKS)*BVECT(IWKS)
  110 CONTINUE
      DPROC = NPROC
      IF (GOPARR) CALL DSCAL(NLARGE,DPROC,ADIAGS,1)
C
C     ----- GET OFF-DIAGONAL CI-ELEMENTS CONTRIBUTION -----
C
      NUWK = 0
      CALL SEQREW(LUMATA)
      READ (LUMATA)
C
  120 CONTINUE
         CALL PKREAD(LUMATA,XX1,IX,IXX,MX,NEMEMX)
         NX = ABS(MX)
         IF (NX .EQ. 0) GO TO 390
         DO 380 N = 1,NX
            VAL = XX1(N)
C
            NPACK = N
*I32        IPACK = IX( 2*NPACK - 1 )
*I32        JPACK = IX( 2*NPACK     )
*I64        LABEL = IX(NPACK)
*I64        IPACK = ISHFT( LABEL, -32 )
*I64        JPACK = IAND( LABEL, MASK32 )
            IUWK = IPACK
            JUWK = JPACK
C
            IF (IUWK.EQ.0) THEN
               NUWK = INT(VAL)
               GO TO 380
            END IF
C
*I32        NLWK = IXX(N)
*I64        IF (MOD(N,2) .EQ. 0) THEN
*I64          LABEL = IXX( N/2 )
*I64          NLWK  = ISHFT( LABEL, -32 )
*I64        ELSE
*I64          LABEL = IXX( N/2+1 )
*I64          NLWK  = IAND( LABEL, MASK32 )
*I64        END IF
C
C                  MOST OF THE TIME, NUWK AND NLWK ARE BOTH 1
C
            IF(NUWK.GT.1  .OR.  NLWK.GT.1) GO TO 140
            II=INDX(IUWK)
            JJ=INDX(JUWK)
            AXB(II) = AXB(II) + VAL*BVECT(JJ)
            AXB(JJ) = AXB(JJ) + VAL*BVECT(II)
            GO TO 380
C
C                  CASE FOR NLWK=1, NUWK>1
C
  140       CONTINUE
            IF(NLWK.GT.1) GO TO 200
            II=INDX(IUWK)
            JJ=INDX(JUWK)
            DO 150 J=1,NUWK
               AXB(II) = AXB(II) + VAL*BVECT(JJ)
               AXB(JJ) = AXB(JJ) + VAL*BVECT(II)
               II=II+1
               JJ=JJ+1
  150       CONTINUE
            GO TO 380
C
C                  CASE FOR NLWK>1, NUWK=1
C
  200       CONTINUE
            IF(NUWK.GT.1) GO TO 240
            DO 220 I=1,NLWK
               II=INDX(IUWK)
               JJ=INDX(JUWK)
               AXB(II) = AXB(II) + VAL*BVECT(JJ)
               AXB(JJ) = AXB(JJ) + VAL*BVECT(II)
               IUWK=IUWK+1
               JUWK=JUWK+1
  220       CONTINUE
            GO TO 380
C
C                  GENERAL CASE FOR NLWK>1, NUWK>1
C
  240       CONTINUE
            DO 260 I=1,NLWK
               II=INDX(IUWK)
               JJ=INDX(JUWK)
               DO 250 J=1,NUWK
                  AXB(II) = AXB(II) + VAL*BVECT(JJ)
                  AXB(JJ) = AXB(JJ) + VAL*BVECT(II)
                  II=II+1
                  JJ=JJ+1
  250          CONTINUE
               IUWK=IUWK+1
               JUWK=JUWK+1
  260       CONTINUE
C
  380    CONTINUE
      IF (MX .GT. 0) GO TO 120
  390 CONTINUE
      IF (GOPARR) CALL DDI_GSUMF(1865,AXB,NLARGE)
      RETURN
C     ******
C
C              IMPROVE MORE THAN ONE STATE
C
C
  400 CONTINUE
C
C     ----- GET DIAGONAL CI-ELEMENTS CONTRIBUTION -----
C
      IF (GOPARR) CALL DSCAL(NLARGE,ONE/NPROC,ADIAGS,1)
      IJ=0
      DO 420 IMPRV=1,NIMPRV
         DO 410 IWKS=1,NLARGE
            IJ=IJ+1
            AXB(IJ) = ADIAGS(IWKS)*BVECT(IJ)
  410    CONTINUE
  420 CONTINUE
      DPROC = NPROC
      IF (GOPARR) CALL DSCAL(NLARGE,DPROC,ADIAGS,1)
C
C     ----- GET OFF-DIAGONAL CI-ELEMENTS CONTRIBUTION -----
C
      NUWK = 0
      CALL SEQREW(LUMATA)
      READ (LUMATA)
  500 CONTINUE
         CALL PKREAD(LUMATA,XX1,IX,IXX,MX,NEMEMX)
         NX = ABS(MX)
         IF (NX .EQ. 0) GO TO 890
         DO 880 N=1,NX
            VAL = XX1(N)
C
            NPACK = N
*I32        IPACK = IX( 2*NPACK - 1 )
*I32        JPACK = IX( 2*NPACK     )
*I64        LABEL = IX(NPACK)
*I64        IPACK = ISHFT( LABEL, -32 )
*I64        JPACK = IAND( LABEL, MASK32 )
            IUWK = IPACK
            JUWK = JPACK
C
            IF(IUWK.EQ.0) THEN
               NUWK = INT(VAL)
               GO TO 880
            END IF
C
*I32        NLWK = IXX(N)
*I64        IF (MOD(N,2) .EQ. 0) THEN
*I64           LABEL = IXX( N/2 )
*I64           NLWK  = ISHFT( LABEL, -32 )
*I64        ELSE
*I64           LABEL = IXX( N/2+1 )
*I64           NLWK  = IAND( LABEL, MASK32 )
*I64        END IF
C
C               MOST OF THE TIME, NLWK AND NUWK ARE 1
C
            IF(NUWK.GT.1  .OR.  NLWK.GT.1) GO TO 600
            II = INDX(IUWK)
            JJ = INDX(JUWK)
            NSTAT=0
            DO 550 IMPRV=1,NIMPRV
               AXB(II+NSTAT) = AXB(II+NSTAT)+VAL*BVECT(JJ+NSTAT)
               AXB(JJ+NSTAT) = AXB(JJ+NSTAT)+VAL*BVECT(II+NSTAT)
               NSTAT = NSTAT+NLARGE
  550       CONTINUE
            GO TO 880
C
C                  CASE FOR NLWK=1, NUWK>1
C
  600       CONTINUE
            IF(NLWK.GT.1) GO TO 700
            II = INDX(IUWK)
            JJ = INDX(JUWK)
            DO 660 J=1,NUWK
               NSTAT=0
               DO 650 IMPRV=1,NIMPRV
                  AXB(II+NSTAT) = AXB(II+NSTAT)+VAL*BVECT(JJ+NSTAT)
                  AXB(JJ+NSTAT) = AXB(JJ+NSTAT)+VAL*BVECT(II+NSTAT)
                  NSTAT=NSTAT+NLARGE
  650          CONTINUE
               II=II+1
               JJ=JJ+1
  660       CONTINUE
            GO TO 880
C
C                  CASE FOR NLWK>1, NUWK=1
C
  700       CONTINUE
            IF(NUWK.GT.1) GO TO 800
            DO 770 I=1,NLWK
               II = INDX(IUWK)
               JJ = INDX(JUWK)
               NSTAT=0
               DO 750 IMPRV=1,NIMPRV
                  AXB(II+NSTAT) = AXB(II+NSTAT)+VAL*BVECT(JJ+NSTAT)
                  AXB(JJ+NSTAT) = AXB(JJ+NSTAT)+VAL*BVECT(II+NSTAT)
                  NSTAT=NSTAT+NLARGE
  750          CONTINUE
               IUWK=IUWK+1
               JUWK=JUWK+1
  770       CONTINUE
            GO TO 880
C
C                  GENERAL CASE, NUWK>1 AND NLWK>1
C
  800       CONTINUE
            DO 870 I=1,NLWK
               II = INDX(IUWK)
               JJ = INDX(JUWK)
               DO 860 J=1,NUWK
                  NSTAT=0
                  DO 850 IMPRV=1,NIMPRV
                     AXB(II+NSTAT) = AXB(II+NSTAT)+VAL*BVECT(JJ+NSTAT)
                     AXB(JJ+NSTAT) = AXB(JJ+NSTAT)+VAL*BVECT(II+NSTAT)
                     NSTAT=NSTAT+NLARGE
  850             CONTINUE
                  II=II+1
                  JJ=JJ+1
  860          CONTINUE
               IUWK=IUWK+1
               JUWK=JUWK+1
  870       CONTINUE
  880    CONTINUE
      IF (MX .GT. 0) GO TO 500
  890 CONTINUE
      IF (GOPARR) CALL DDI_GSUMF(1865,AXB,NLARGE*NIMPRV)
      RETURN
      END
C*MODULE GUGDGB  *DECK EDTCHK
C*  8 FEB 1984 - STE
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)  DATE: NOV 1983
C*
C*    PURPOSE -
C*       WHEN FINDING EIGENVECTOR(S) USING THE DAVIDSON METHOD, THIS
C*                    *                        *
C*       ROUTINE CALCULATES EIGENVALUE TOLERANCES AND CHECKS
C*                                     *              **  *
C*       FOR CONVERGENCE IF THE CURRENT CONTRACTION VECTOR(S)
C*       AGREES WITH THE PREVIOUS CONTRACTION VECTOR(S) AND IMPROVEMENT
C*       HAS BEEN MADE.
C*
C*    ON ENTRY -
C*       NNEW   - INTEGER
C*                DIMENSION OF CURRENT CONTRACTION VECTOR (CVECT).
C*                (NUMBER OF CURRENT BASIS VECTORS).
C*       NNOLD  - INTEGER
C*                DIMENSION OF PREVIOUS CONTRACTION VECTOR (C0VECT).
C*                NNOLD = NNEW - NIMPRV
C*       IDLORT - INTEGER
C*                NO. OF LOWEST ROOT SOUGHT.
C*       IDHIRT - INTEGER
C*                NO. OF HIGHEST ROOT SOUGHT.
C*       MXMPRV - INTEGER
C*                MAXIMUM NO. OF ROOTS THAT MAY BE IMPROVED.
C*       CVGETH - W.P. REAL
C*                EIGENVALUE CONVERGENCE THRESHOLD
C*       NIMPRV - INTEGER
C*                NO. OF ROOTS THAT HAVE BEEN IMPROVED IN CURRENT ITER.
C*       MPRVME - INTEGER (NIMPRV)
C*                LIST OF ROOTS BEING IMPROVED THIS ITERATION.
C*       CVECT  - W.P. REAL (NNEW,IDHIRT)
C*                CURRENT CONTRACTION VECTOR(S).
C*                (EIGENVECTORS OF B(T) * A * B)
C*       C0VECT - W.P. REAL (NOLD,IDHIRT)
C*                PREVIOUS CONTRACTION VECTOR(S).
C*                (CVECT FROM LAST ITERATION)
C*       IMCONV - INTEGER (IDHIRT)
C*                IS 1 IF ROOT HAS ALREADY CONVERGED.
C*
C*    ON EXIT -
C*       NIMPRV - INTEGER
C*                NO. OF ROOTS TO BE IMPROVED NEXT ITERATION.
C*       MPRVME - INTEGER (NIMPRV<=MXMPRV)
C*                LIST OF ROOTS TO BE IMPROVED NEXT ITERATION.
C*       TOLERS - W.P. REAL (IDHIRT)
C*                EIGENVALUE ACCURACY (SUM OF THE SQUARES OF THE LAST
C*                NIMPRV ELEMENTS OF EACH CONTRACTION VECTOR).
C*                HAS A NEGATIVE VALUE IF ROOT SELECTED FOR IMPROVEMENT
C*       IMCONV - INTEGER (IDHIRT)
C*                SET TO ONE WHEN A ROOT IS CONVERGED.
C*       CVGART - LOGICAL
C*                SET TO .TRUE. IF A ROOT HAS CONVERGED THIS ITERATION,
C*                OTHERWISE, SET TO .FALSE.
C*       VERMAX - W.P. REAL
C*                LARGEST ERROR FOR A CONVERGED VECTOR
C*       IDMAX  - INTEGER
C*                NUMBER OF THE CONVERGED VECTOR WITH LARGEST ERROR
C*
      SUBROUTINE EDTCHK(NNEW,NOLD,IDLORT,IDHIRT,MXMPRV,CVGETH,
     *                  NIMPRV,MPRVME,CVECT,C0VECT,IMCONV,TOLERS,CVGART,
     *                  VERMAX,IDMAX,ITER)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL CVGART
C
      DIMENSION CVECT(NNEW,IDHIRT), C0VECT(NOLD,IDHIRT)
     *         ,IMCONV(IDHIRT), TOLERS(IDHIRT), MPRVME(MXMPRV)
C
      PARAMETER (ZERO=0.0D+00, POINT9=0.9D+00)
C
C           COMPUTE EIGENVALUE TOLERANCES
C
      DO 110 I=IDLORT,IDHIRT
         TOLERS(I)=DDOT(NIMPRV,CVECT(NOLD+1,I),1,CVECT(NOLD+1,I),1)
  110 CONTINUE
C
C           ROOTS THAT CHANGE POSITION ARE NOT CONSIDERED CONVERGED
C
C           CHECK STATUS OF PREVIOUSLY "CONVERGED" ROOTS
C
      IF(ITER.LE.3) GO TO 135
      DO 120 I=IDLORT,IDHIRT
         IF(IMCONV(I).EQ.0) GO TO 120
            IMCONV(I)=0
            SUM1=DDOT(NOLD,CVECT(1,I),1,C0VECT(1,I),1)
            IF(ABS(SUM1).LE. POINT9) GO TO 120
               IF(TOLERS(I).LE.CVGETH) IMCONV(I)=1
  120 CONTINUE
C
C           CHECK STATUS OF ROOTS JUST IMPROVED
C
      DO 130 N=1,NIMPRV
         NT=MPRVME(N)
         IMCONV(NT)=0
         SUM2=DDOT(NOLD,CVECT(1,NT),1,C0VECT(1,NT),1)
         IF(ABS(SUM2).LE.POINT9) GO TO 130
            IF(TOLERS(NT).GT.CVGETH) GO TO 130
               IMCONV(NT)=1
               CVGART=.TRUE.
  130    CONTINUE
C
C        FIND UNCONVERGED ROOTS WITH LARGEST TOLERANCES
C          ORIGINALLY, CMAX WAS SET TO ZERO
C
  135 CONTINUE
      NIMPRV=0
      DO 150 N=1,MXMPRV
         NT=0
         CMAX=0.1D+00*CVGETH
         IF(ITER.LE.3) CMAX=ZERO
         DO 140 I=IDLORT,IDHIRT
            IF(IMCONV(I).EQ.1) GO TO 140
               IF(TOLERS(I).LT.CMAX) GO TO 140
                  CMAX=TOLERS(I)
                  NT=I
  140    CONTINUE
C
         IF(NT.EQ.0) GO TO 160
C
C              MARK THE CHOSEN ROOT
C
            NIMPRV=NIMPRV+1
            TOLERS(NT)=-TOLERS(NT)
            MPRVME(NIMPRV)=NT
            IF(N.GT.1) GO TO 150
               VERMAX=SQRT(CMAX)
               IDMAX =NT
  150 CONTINUE
      RETURN
C
  160 CONTINUE
      IF(NIMPRV.GT.0) RETURN
C
C           GET WORST TOLERANCE OF CONVERGED VECTORS
C
      VERMAX=ZERO
      DO 180 I=IDLORT,IDHIRT
         IF(TOLERS(I).LT.VERMAX) GO TO 180
            VERMAX=TOLERS(I)
            IDMAX=I
  180 CONTINUE
      VERMAX=SQRT(VERMAX)
      RETURN
      END
C*MODULE GUGDGB  *DECK EDTRUN
C* 14 FEB 1984 - STE
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)  DATE: NOV 1982
C*
C*    PURPOSE -
C*       WHEN FINDING EIGENVECTOR(S) USING THE DAVIDSON METHOD, THIS
C*                    *                        *
C*       ROUTINE TRUNCATES THE EXPANSION OF APPROXIMATE EIGENVECTORS,
C*               ****
C*       WHICH, IN EFFECT, RESTARTS THE PROBLEM USING THE FIRST IDHIRT
C*       APPROXIMATE VECTORS AS THE BASIS VECTORS.  IF CONVERGENCE HAS
C*       BEEN ACHIEVED, THE NEW BASIS VECTORS ARE THE FINAL VECTORS.
C*
C*    ON ENTRY -
C*       NLARGE - INTEGER
C*                DIMENSION OF THE MATRIX (A) BEING SOLVED.
C*       NBASIS - INTEGER
C*                CURRENT NUMBER OF BASIS VECTORS.
C*       IDHIRT - INTEGER
C*                NO. OF THE HIGHEST ROOT BEING SOUGHT.
C*       NHITRI - INTEGER
C*                WORDS NEEDED TO STORE BTAB (=IDHIRT*(IDHIRT+1)/2)
C*       CVECT  - W.P. REAL (NBASIS,IDHIRT)
C*                EIGENVECTORS OF B(TRANSPOSE) * A * B.
C*                (NOTE DIMENSION CHANGE ON EXIT)
C*       BVECT  - W.P. REAL (NLARGE,NBASIS)
C*                CURRENT SET OF BASIS VECTORS.
C*       AXB    - W.P. REAL (NLARGE,NBASIS)
C*                PRODUCT OF INPUT MATRIX (A) AND CURRENT BASIS VECTORS.
C*
C*    ON EXIT -
C*       IJBTAB - INTEGER
C*                = NHITRI
C*       CVECT  - W.P. REAL (IDHIRT,IDHIRT) <-- NOT THE SAME AS ENTRY
C*                UNIT MATRIX.
C*       BVECT  - W.P. REAL (NLARGE,IDHIRT)
C*                APPROX. EIGENVECTORS OF MATRIX A (NEW BASIS VECTORS)
C*       AXB    - W.P. REAL (NLARGE,IDHIRT)
C*                PRODUCT OF INPUT MATRIX (A) AND NEW BASIS VECTORS.
C*       BTAB   - W.P. REAL (NHITRI)
C*                INPUT MATRIX (A) TRANSFORMED BY NEW BASIS VECTORS.
C*       WRK    - W.P. REAL (NLARGE,IDHIRT)
C*                WORK SPACE TO HOLD PRODUCT OF (BVECT OR AXB) * CVECT
C*
      SUBROUTINE EDTRUN(LUVECS,NLARGE,NBASIS,IDHIRT,NHITRI,IJBTAB,
     *                  CVECT,BVECT,AXB,BTAB,WRK)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      COMMON /IOFILE/ IR,IW,IP,IS,IPK,IDAF,NAV,IODA(950)
      LOGICAL GOPARR,DSKWRK,MASWRK
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
      DIMENSION BVECT(NLARGE,NBASIS),AXB(NLARGE,NBASIS)
     *         ,BTAB(NHITRI), CVECT(NBASIS,IDHIRT)
     *         ,WRK(NLARGE,IDHIRT)
C
C           MAKE SURE THE CONTRACTION VECTORS ARE ORTHONORMAL
C
      NDROP = 0
      CALL ONVMGS(NBASIS,IDHIRT,CVECT,NDROP)
      IF(NDROP.NE.0)THEN
          IF(MASWRK)WRITE(IW,*)
     *        ' CONTRACTION VECTORS ARE NOT ORTHONORMAL IN EDTRUN'
          CALL ABRT
      END IF
C
C           UPDATE MATRIX * VECTOR
C
      LUAXB=LUVECS
      IF(LUAXB.GT.0) LUAXB = LUVECS + 3
C
      CALL FRARBR(LUAXB,NLARGE,NBASIS,IDHIRT,AXB,CVECT,WRK)
      IF(LUAXB.GT.0) THEN
         CALL SEQREW(LUAXB)
         DO 120 I=1,IDHIRT
            CALL SQWRIT(LUAXB,WRK(1,I),NLARGE)
  120    CONTINUE
      ELSE
         CALL DCOPY(NLARGE*IDHIRT,WRK,1,AXB,1)
      END IF
C
C           UPDATE VECTOR
C
      CALL FRARBR(LUVECS,NLARGE,NBASIS,IDHIRT,BVECT,CVECT,WRK)
      IF(LUVECS.GT.0) THEN
         CALL SEQREW(LUVECS)
         WRITE(LUVECS) IDHIRT,NLARGE
         DO 220 I=1,IDHIRT
            CALL SQWRIT(LUVECS,WRK(1,I),NLARGE)
  220    CONTINUE
      ELSE
         CALL DCOPY(NLARGE*IDHIRT,WRK,1,BVECT,1)
      END IF
C
C        RECALCULATE  A  IN THE EXPANSION BASIS SPACE
C              (BVECT IS STILL IN WRK)
C
      IF(LUVECS.GT.0) THEN
         CALL SEQREW(LUAXB)
         DO 240 I=1,IDHIRT
            CALL SQREAD(LUAXB,AXB(1,I),NLARGE)
  240    CONTINUE
      END IF
      IJBTAB=0
      CALL FRTRBT(LUVECS,NLARGE,1,IDHIRT,NHITRI,WRK,AXB,IJBTAB,BTAB)
C
C           SET CONTRACTION VECTORS TO UNIT VECTORS
C
      CALL RUNITV(IDHIRT,IDHIRT,CVECT)
      RETURN
      END
C*MODULE GUGDGB  *DECK FRARBR
C* 23 MAY 1983 - STE
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)  DATE: OCT 1982
C*
C*    PURPOSE -
C*       FORM PRODUCT OF RECTANGULAR MATRIX  A  AND RECTANGULAR MATRIX B
C*            *          *                   *      *                  *
C*       TO GET RECTANGULAR MATRIX  C.
C*              *
C*       (A MAY NOT BE MEMORY RESIDENT)
C*
C*    ON ENTRY -
C*       LUNIT - INTEGER (LOGICAL UNIT NUMBER)
C*                IF ZERO,  A  IS MEMORY RESIDENT, OTHERWISE  A  IS
C*                ON LOGICAL UNIT LUNIT, ONE COLUMN PER RECORD.
C*       NA     - INTEGER
C*                ROW DIMENSION OF MATRICES  A  AND  C
C*       NB     - INTEGER
C*                COLUMN DIMENSION OF MATRIX  A
C*                ROW DIMENSION OF MATRIX  B
C*       NC     - INTEGER
C*                COLUMN DIMENSION OF MATRICES  B  AND  C
C*       A      - W.P. REAL (NA,NB)
C*                LEFT HAND INPUT MATRIX
C*       B      - W.P. REAL (NB,NC)
C*                RIGHT HAND INPUT MATRIX
C*
C*    ON EXIT -
C*       C      - W.P. REAL (NA,NC)
C*                PRODUCT OF  A * B
C*       --- FOR THE CASE WHERE  A  IS NOT MEMORY RESIDENT ---
C*       A      - W.P. REAL (NA,1)
C*                LAST VECTOR OF LEFT HAND INPUT MATRIX.
C*
C*    COMPLEXITY -
C*       NA * NB * NC
C*
C*    NOTE -
C*       THIS ROUTINE IS MOST EFFICIENT IF  NA > NB,NC
C*
      SUBROUTINE FRARBR(LUNIT,NA,NB,NC,A,B,C)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION A(NA,*), B(NB,NC), C(NA,NC)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      IF(LUNIT.GT.0) THEN
         CALL SEQREW(LUNIT)
         IF(LUNIT.EQ.12) READ(LUNIT)
      END IF
C
      CALL VCLR(C,1,NA*NC)
C
      DO 120 K=1,NB
         IF(LUNIT.GT.0) THEN
            I=1
            CALL SQREAD(LUNIT,A,NA)
            IF(NA.EQ.0) THEN
               IF (MASWRK) WRITE(6,*)
     *            'FRARBR: UNEXPECTED EOF ON UNI',LUNIT
               CALL ABRT
               STOP
            END IF
         ELSE
            I=K
         END IF
         DO 110 J=1,NC
            CALL DAXPY(NA,B(K,J),A(1,I),1,C(1,J),1)
  110    CONTINUE
  120 CONTINUE
      RETURN
      END
C*MODULE GUGDGB  *DECK FRTRBT
C* 14 FEB 1984 - STE
C*
C*    AUTHOR: S. T. ELBERT (AMES LABORATORY - USDOE)  DATE:  6 APR 1983
C*
C*    PURPOSE -
C*       FORM PRODUCT OF RECTANGULAR MATRIX  A(TRANSPOSE)  AND
C*       *               *                     *
C*       RECTANGULAR MATRIX  B  TO GET SYMMETRIC (TRIANGULAR) MATRIX  C.
C*       *                   *                    *
C*       (A MAY NOT BE MEMORY RESIDENT)
C*
C*    ON ENTRY -
C*       LUVECS - INTEGER (LOGICAL UNIT NUMBER)
C*                IF NEGATIVE,  A  IS MEMORY RESIDENT, ELSE, LOGICAL
C*                UNIT WHERE A IS STORED ONE VECTOR PER RECORD.
C*       NROWS  - INTEGER
C*                ROW DIMENSION OF MATRICES  A,  B  AND  C.
C*       NCOLLO - INTEGER
C*                FIRST VECTOR OF MATRIX  B  USED IN PRODUCT.
C*       NCOLHI - INTEGER
C*                COLUMN DIMENSION OF MATRICES  A  AND  B.
C*       LENC   - INTEGER
C*                W.P. WORDS OF STORAGE AVAILABLE TO  C.
C*       A      - W.P. REAL (NROWS,NCOLHI)
C*                LEFT HAND INPUT MATRIX.
C*       B      - W.P. REAL (NROWS,NCOLHI)
C*                RIGHT HAND INPUT MATRIX.
C*       IJ     - INTEGER
C*                LAST PREVIOUSLY CALCULATED ELEMENT C(I,J)
C*                IJ = I*(I-1)/2  +  J
C*
C*    ON EXIT -
C*       IJ     - INTEGER
C*                LAST ELEMENT OF  C  JUST CALCULATED.
C*       C      - W.P. REAL (LENC)
C*                OUTPUT MATRIX IN SYMMETRIC STORAGE MODE.
C*       --- FOR THE CASE WHERE  A  IS NOT MEMORY RESIDENT ---
C*       A      - W.P. REAL (NROWS,1)
C*                LAST VECTOR OF LEFT HAND INPUT MATRIX.
C*
C*    COMPLEXITY -
C*       NROWS * (NCOLHI*(NCOLHI+1)/2 - NCOLLO*(NCOLLO-1)/2)
C*
      SUBROUTINE FRTRBT(LUVECS,NROWS,NCOLLO,NCOLHI,LENC,A,B,IJ,C)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C
      LOGICAL GOPARR,DSKWRK,MASWRK
C
      DIMENSION A(NROWS,*), B(NROWS,NCOLHI), C(LENC)
C
      COMMON /PAR   / ME,MASTER,NPROC,IBTYP,IPTIM,GOPARR,DSKWRK,MASWRK
C
      IJ0=IJ
      NLOM1=NCOLLO-1
      N0=0
      IF(LUVECS.GT.0) N0=NLOM1
C
C              DO TRIANGULAR PORTION
C
      DO 120 J=NCOLLO,NCOLHI
         IJ=IJ+NLOM1
         DO 110 I=NCOLLO,J
            IJ=IJ+1
            C(IJ)=DDOT(NROWS,A(1,I-N0),1,B(1,J-N0),1)
  110    CONTINUE
  120 CONTINUE
      IF(NCOLLO.EQ.1) RETURN
C
C              DO RECTANGULAR PORTION
C
      IF(LUVECS.GT.0) THEN
         CALL SEQREW(LUVECS)
         READ(LUVECS)
      END IF
      DO 160 I=1,NLOM1
         IJ0=IJ0+1
         IJR=IJ0
         IF(LUVECS.GT.0) THEN
            M=1
            CALL SQREAD(LUVECS,A,NROWS)
            IF(NROWS.EQ.0) THEN
               IF (MASWRK) WRITE(6,*)
     *            'ENEXPECTED EOF ON CI VECTOR FILE'
               CALL ABRT
               STOP
            END IF
         ELSE
            M=I
         END IF
         DO 150 J=NCOLLO,NCOLHI
            C(IJR)=DDOT(NROWS,A(1,M),1,B(1,J-N0),1)
            IJR=IJR+J
  150    CONTINUE
  160 CONTINUE
      RETURN
      END
