text:
*> \brief \b DDRVGBX
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA,
*                          AFB, LAFB, ASAV, B, BSAV, X, XACT, S, WORK,
*                          RWORK, IWORK, NOUT )
*
*       .. Scalar Arguments ..
*       LOGICAL            TSTERR
*       INTEGER            LA, LAFB, NN, NOUT, NRHS
*       DOUBLE PRECISION   THRESH
*       ..
*       .. Array Arguments ..
*       LOGICAL            DOTYPE( * )
*       INTEGER            IWORK( * ), NVAL( * )
*       DOUBLE PRECISION   A( * ), AFB( * ), ASAV( * ), B( * ), BSAV( * ),
*      $                   RWORK( * ), S( * ), WORK( * ), X( * ),
*      $                   XACT( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DDRVGB tests the driver routines DGBSV, -SVX, and -SVXX.
*>
*> Note that this file is used only when the XBLAS are available,
*> otherwise ddrvgb.f defines this subroutine.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] DOTYPE
*> \verbatim
*>          DOTYPE is LOGICAL array, dimension (NTYPES)
*>          The matrix types to be used for testing.  Matrices of type j
*>          (for 1 <= j <= NTYPES) are used for testing if DOTYPE(j) =
*>          .TRUE.; if DOTYPE(j) = .FALSE., then type j is not used.
*> \endverbatim
*>
*> \param[in] NN
*> \verbatim
*>          NN is INTEGER
*>          The number of values of N contained in the vector NVAL.
*> \endverbatim
*>
*> \param[in] NVAL
*> \verbatim
*>          NVAL is INTEGER array, dimension (NN)
*>          The values of the matrix column dimension N.
*> \endverbatim
*>
*> \param[in] NRHS
*> \verbatim
*>          NRHS is INTEGER
*>          The number of right hand side vectors to be generated for
*>          each linear system.
*> \endverbatim
*>
*> \param[in] THRESH
*> \verbatim
*>          THRESH is DOUBLE PRECISION
*>          The threshold value for the test ratios.  A result is
*>          included in the output file if RESULT >= THRESH.  To have
*>          every test ratio printed, use THRESH = 0.
*> \endverbatim
*>
*> \param[in] TSTERR
*> \verbatim
*>          TSTERR is LOGICAL
*>          Flag that indicates whether error exits are to be tested.
*> \endverbatim
*>
*> \param[out] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LA)
*> \endverbatim
*>
*> \param[in] LA
*> \verbatim
*>          LA is INTEGER
*>          The length of the array A.  LA >= (2*NMAX-1)*NMAX
*>          where NMAX is the largest entry in NVAL.
*> \endverbatim
*>
*> \param[out] AFB
*> \verbatim
*>          AFB is DOUBLE PRECISION array, dimension (LAFB)
*> \endverbatim
*>
*> \param[in] LAFB
*> \verbatim
*>          LAFB is INTEGER
*>          The length of the array AFB.  LAFB >= (3*NMAX-2)*NMAX
*>          where NMAX is the largest entry in NVAL.
*> \endverbatim
*>
*> \param[out] ASAV
*> \verbatim
*>          ASAV is DOUBLE PRECISION array, dimension (LA)
*> \endverbatim
*>
*> \param[out] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (NMAX*NRHS)
*> \endverbatim
*>
*> \param[out] BSAV
*> \verbatim
*>          BSAV is DOUBLE PRECISION array, dimension (NMAX*NRHS)
*> \endverbatim
*>
*> \param[out] X
*> \verbatim
*>          X is DOUBLE PRECISION array, dimension (NMAX*NRHS)
*> \endverbatim
*>
*> \param[out] XACT
*> \verbatim
*>          XACT is DOUBLE PRECISION array, dimension (NMAX*NRHS)
*> \endverbatim
*>
*> \param[out] S
*> \verbatim
*>          S is DOUBLE PRECISION array, dimension (2*NMAX)
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension
*>                      (NMAX*max(3,NRHS,NMAX))
*> \endverbatim
*>
*> \param[out] RWORK
*> \verbatim
*>          RWORK is DOUBLE PRECISION array, dimension
*>                      (max(2*NMAX,NMAX+2*NRHS))
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (2*NMAX)
*> \endverbatim
*>
*> \param[in] NOUT
*> \verbatim
*>          NOUT is INTEGER
*>          The unit number for output.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \ingroup double_lin
*
*  =====================================================================
      SUBROUTINE DDRVGB( DOTYPE, NN, NVAL, NRHS, THRESH, TSTERR, A, LA,
     $                   AFB, LAFB, ASAV, B, BSAV, X, XACT, S, WORK,
     $                   RWORK, IWORK, NOUT )
*
*  -- LAPACK test routine --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*
*     .. Scalar Arguments ..
      LOGICAL            TSTERR
      INTEGER            LA, LAFB, NN, NOUT, NRHS
      DOUBLE PRECISION   THRESH
*     ..
*     .. Array Arguments ..
      LOGICAL            DOTYPE( * )
      INTEGER            IWORK( * ), NVAL( * )
      DOUBLE PRECISION   A( * ), AFB( * ), ASAV( * ), B( * ), BSAV( * ),
     $                   RWORK( * ), S( * ), WORK( * ), X( * ),
     $                   XACT( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
      INTEGER            NTYPES
      PARAMETER          ( NTYPES = 8 )
      INTEGER            NTESTS
      PARAMETER          ( NTESTS = 7 )
      INTEGER            NTRAN
      PARAMETER          ( NTRAN = 3 )
*     ..
*     .. Local Scalars ..
      LOGICAL            EQUIL, NOFACT, PREFAC, TRFCON, ZEROT
      CHARACTER          DIST, EQUED, FACT, TRANS, TYPE, XTYPE
      CHARACTER*3        PATH
      INTEGER            I, I1, I2, IEQUED, IFACT, IKL, IKU, IMAT, IN,
     $                   INFO, IOFF, ITRAN, IZERO, J, K, K1, KL, KU,
     $                   LDA, LDAFB, LDB, MODE, N, NB, NBMIN, NERRS,
     $                   NFACT, NFAIL, NIMAT, NKL, NKU, NRUN, NT,
     $                   N_ERR_BNDS
      DOUBLE PRECISION   AINVNM, AMAX, ANORM, ANORMI, ANORMO, ANRMPV,
     $                   CNDNUM, COLCND, RCOND, RCONDC, RCONDI, RCONDO,
     $                   ROLDC, ROLDI, ROLDO, ROWCND, RPVGRW,
     $                   RPVGRW_SVXX
*     ..
*     .. Local Arrays ..
      CHARACTER          EQUEDS( 4 ), FACTS( 3 ), TRANSS( NTRAN )
      INTEGER            ISEED( 4 ), ISEEDY( 4 )
      DOUBLE PRECISION   RESULT( NTESTS ), BERR( NRHS ),
     $                   ERRBNDS_N( NRHS, 3 ), ERRBNDS_C( NRHS, 3 )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DGET06, DLAMCH, DLANGB, DLANGE, DLANTB,
     $                   DLA_GBRPVGRW
      EXTERNAL           LSAME, DGET06, DLAMCH, DLANGB, DLANGE, DLANTB,
     $                   DLA_GBRPVGRW
*     ..
*     .. External Subroutines ..
      EXTERNAL           ALADHD, ALAERH, ALASVM, DERRVX, DGBEQU, DGBSV,
     $                   DGBSVX, DGBT01, DGBT02, DGBT05, DGBTRF, DGBTRS,
     $                   DGET04, DLACPY, DLAQGB, DLARHS, DLASET, DLATB4,
     $                   DLATMS, XLAENV, DGBSVXX, DGBEQUB
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, MAX, MIN
*     ..
*     .. Scalars in Common ..
      LOGICAL            LERR, OK
      CHARACTER*32       SRNAMT
      INTEGER            INFOT, NUNIT
*     ..
*     .. Common blocks ..
      COMMON             / INFOC / INFOT, NUNIT, OK, LERR
      COMMON             / SRNAMC / SRNAMT
*     ..
*     .. Data statements ..
      DATA               ISEEDY / 1988, 1989, 1990, 1991 /
      DATA               TRANSS / 'N', 'T', 'C' /
      DATA               FACTS / 'F', 'N', 'E' /
      DATA               EQUEDS / 'N', 'R', 'C', 'B' /
*     ..
*     .. Executable Statements ..
*
*     Initialize constants and the random number seed.
*
      PATH( 1: 1 ) = 'Double precision'
      PATH( 2: 3 ) = 'GB'
      NRUN = 0
      NFAIL = 0
      NERRS = 0
      DO 10 I = 1, 4
         ISEED( I ) = ISEEDY( I )
   10 CONTINUE
*
*     Test the error exits
*
      IF( TSTERR )
     $   CALL DERRVX( PATH, NOUT )
      INFOT = 0
*
*     Set the block size and minimum block size for testing.
*
      NB = 1
      NBMIN = 2
      CALL XLAENV( 1, NB )
      CALL XLAENV( 2, NBMIN )
*
*     Do for each value of N in NVAL
*
      DO 150 IN = 1, NN
         N = NVAL( IN )
         LDB = MAX( N, 1 )
         XTYPE = 'N'
*
*        Set limits on the number of loop iterations.
*
         NKL = MAX( 1, MIN( N, 4 ) )
         IF( N.EQ.0 )
     $      NKL = 1
         NKU = NKL
         NIMAT = NTYPES
         IF( N.LE.0 )
     $      NIMAT = 1
*
         DO 140 IKL = 1, NKL
*
*           Do for KL = 0, N-1, (3N-1)/4, and (N+1)/4. This order makes
*           it easier to skip redundant values for small values of N.
*
            IF( IKL.EQ.1 ) THEN
               KL = 0
            ELSE IF( IKL.EQ.2 ) THEN
               KL = MAX( N-1, 0 )
            ELSE IF( IKL.EQ.3 ) THEN
               KL = ( 3*N-1 ) / 4
            ELSE IF( IKL.EQ.4 ) THEN
               KL = ( N+1 ) / 4
            END IF
            DO 130 IKU = 1, NKU
*
*              Do for KU = 0, N-1, (3N-1)/4, and (N+1)/4. This order
*              makes it easier to skip redundant values for small
*              values of N.
*
               IF( IKU.EQ.1 ) THEN
                  KU = 0
               ELSE IF( IKU.EQ.2 ) THEN
                  KU = MAX( N-1, 0 )
               ELSE IF( IKU.EQ.3 ) THEN
                  KU = ( 3*N-1 ) / 4
               ELSE IF( IKU.EQ.4 ) THEN
                  KU = ( N+1 ) / 4
               END IF
*
*              Check that A and AFB are big enough to generate this
*              matrix.
*
               LDA = KL + KU + 1
               LDAFB = 2*KL + KU + 1
               IF( LDA*N.GT.LA .OR. LDAFB*N.GT.LAFB ) THEN
                  IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $               CALL ALADHD( NOUT, PATH )
                  IF( LDA*N.GT.LA ) THEN
                     WRITE( NOUT, FMT = 9999 )LA, N, KL, KU,
     $                  N*( KL+KU+1 )
                     NERRS = NERRS + 1
                  END IF
                  IF( LDAFB*N.GT.LAFB ) THEN
                     WRITE( NOUT, FMT = 9998 )LAFB, N, KL, KU,
     $                  N*( 2*KL+KU+1 )
                     NERRS = NERRS + 1
                  END IF
                  GO TO 130
               END IF
*
               DO 120 IMAT = 1, NIMAT
*
*                 Do the tests only if DOTYPE( IMAT ) is true.
*
                  IF( .NOT.DOTYPE( IMAT ) )
     $               GO TO 120
*
*                 Skip types 2, 3, or 4 if the matrix is too small.
*
                  ZEROT = IMAT.GE.2 .AND. IMAT.LE.4
                  IF( ZEROT .AND. N.LT.IMAT-1 )
     $               GO TO 120
*
*                 Set up parameters with DLATB4 and generate a
*                 test matrix with DLATMS.
*
                  CALL DLATB4( PATH, IMAT, N, N, TYPE, KL, KU, ANORM,
     $                         MODE, CNDNUM, DIST )
                  RCONDC = ONE / CNDNUM
*
                  SRNAMT = 'DLATMS'
                  CALL DLATMS( N, N, DIST, ISEED, TYPE, RWORK, MODE,
     $                         CNDNUM, ANORM, KL, KU, 'Z', A, LDA, WORK,
     $                         INFO )
*
*                 Check the error code from DLATMS.
*
                  IF( INFO.NE.0 ) THEN
                     CALL ALAERH( PATH, 'DLATMS', INFO, 0, ' ', N, N,
     $                            KL, KU, -1, IMAT, NFAIL, NERRS, NOUT )
                     GO TO 120
                  END IF
*
*                 For types 2, 3, and 4, zero one or more columns of
*                 the matrix to test that INFO is returned correctly.
*
                  IZERO = 0
                  IF( ZEROT ) THEN
                     IF( IMAT.EQ.2 ) THEN
                        IZERO = 1
                     ELSE IF( IMAT.EQ.3 ) THEN
                        IZERO = N
                     ELSE
                        IZERO = N / 2 + 1
                     END IF
                     IOFF = ( IZERO-1 )*LDA
                     IF( IMAT.LT.4 ) THEN
                        I1 = MAX( 1, KU+2-IZERO )
                        I2 = MIN( KL+KU+1, KU+1+( N-IZERO ) )
                        DO 20 I = I1, I2
                           A( IOFF+I ) = ZERO
   20                   CONTINUE
                     ELSE
                        DO 40 J = IZERO, N
                           DO 30 I = MAX( 1, KU+2-J ),
     $                             MIN( KL+KU+1, KU+1+( N-J ) )
                              A( IOFF+I ) = ZERO
   30                      CONTINUE
                           IOFF = IOFF + LDA
   40                   CONTINUE
                     END IF
                  END IF
*
*                 Save a copy of the matrix A in ASAV.
*
                  CALL DLACPY( 'Full', KL+KU+1, N, A, LDA, ASAV, LDA )
*
                  DO 110 IEQUED = 1, 4
                     EQUED = EQUEDS( IEQUED )
                     IF( IEQUED.EQ.1 ) THEN
                        NFACT = 3
                     ELSE
                        NFACT = 1
                     END IF
*
                     DO 100 IFACT = 1, NFACT
                        FACT = FACTS( IFACT )
                        PREFAC = LSAME( FACT, 'F' )
                        NOFACT = LSAME( FACT, 'N' )
                        EQUIL = LSAME( FACT, 'E' )
*
                        IF( ZEROT ) THEN
                           IF( PREFAC )
     $                        GO TO 100
                           RCONDO = ZERO
                           RCONDI = ZERO
*
                        ELSE IF( .NOT.NOFACT ) THEN
*
*                          Compute the condition number for comparison
*                          with the value returned by DGESVX (FACT =
*                          'N' reuses the condition number from the
*                          previous iteration with FACT = 'F').
*
                           CALL DLACPY( 'Full', KL+KU+1, N, ASAV, LDA,
     $                                  AFB( KL+1 ), LDAFB )
                           IF( EQUIL .OR. IEQUED.GT.1 ) THEN
*
*                             Compute row and column scale factors to
*                             equilibrate the matrix A.
*
                              CALL DGBEQU( N, N, KL, KU, AFB( KL+1 ),
     $                                     LDAFB, S, S( N+1 ), ROWCND,
     $                                     COLCND, AMAX, INFO )
                              IF( INFO.EQ.0 .AND. N.GT.0 ) THEN
                                 IF( LSAME( EQUED, 'R' ) ) THEN
                                    ROWCND = ZERO
                                    COLCND = ONE
                                 ELSE IF( LSAME( EQUED, 'C' ) ) THEN
                                    ROWCND = ONE
                                    COLCND = ZERO
                                 ELSE IF( LSAME( EQUED, 'B' ) ) THEN
                                    ROWCND = ZERO
                                    COLCND = ZERO
                                 END IF
*
*                                Equilibrate the matrix.
*
                                 CALL DLAQGB( N, N, KL, KU, AFB( KL+1 ),
     $                                        LDAFB, S, S( N+1 ),
     $                                        ROWCND, COLCND, AMAX,
     $                                        EQUED )
                              END IF
                           END IF
*
*                          Save the condition number of the
*                          non-equilibrated system for use in DGET04.
*
                           IF( EQUIL ) THEN
                              ROLDO = RCONDO
                              ROLDI = RCONDI
                           END IF
*
*                          Compute the 1-norm and infinity-norm of A.
*
                           ANORMO = DLANGB( '1', N, KL, KU, AFB( KL+1 ),
     $                              LDAFB, RWORK )
                           ANORMI = DLANGB( 'I', N, KL, KU, AFB( KL+1 ),
     $                              LDAFB, RWORK )
*
*                          Factor the matrix A.
*
                           CALL DGBTRF( N, N, KL, KU, AFB, LDAFB, IWORK,
     $                                  INFO )
*
*                          Form the inverse of A.
*
                           CALL DLASET( 'Full', N, N, ZERO, ONE, WORK,
     $                                  LDB )
                           SRNAMT = 'DGBTRS'
                           CALL DGBTRS( 'No transpose', N, KL, KU, N,
     $                                  AFB, LDAFB, IWORK, WORK, LDB,
     $                                  INFO )
*
*                          Compute the 1-norm condition number of A.
*
                           AINVNM = DLANGE( '1', N, N, WORK, LDB,
     $                              RWORK )
                           IF( ANORMO.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
                              RCONDO = ONE
                           ELSE
                              RCONDO = ( ONE / ANORMO ) / AINVNM
                           END IF
*
*                          Compute the infinity-norm condition number
*                          of A.
*
                           AINVNM = DLANGE( 'I', N, N, WORK, LDB,
     $                              RWORK )
                           IF( ANORMI.LE.ZERO .OR. AINVNM.LE.ZERO ) THEN
                              RCONDI = ONE
                           ELSE
                              RCONDI = ( ONE / ANORMI ) / AINVNM
                           END IF
                        END IF
*
                        DO 90 ITRAN = 1, NTRAN
*
*                          Do for each value of TRANS.
*
                           TRANS = TRANSS( ITRAN )
                           IF( ITRAN.EQ.1 ) THEN
                              RCONDC = RCONDO
                           ELSE
                              RCONDC = RCONDI
                           END IF
*
*                          Restore the matrix A.
*
                           CALL DLACPY( 'Full', KL+KU+1, N, ASAV, LDA,
     $                                  A, LDA )
*
*                          Form an exact solution and set the right hand
*                          side.
*
                           SRNAMT = 'DLARHS'
                           CALL DLARHS( PATH, XTYPE, 'Full', TRANS, N,
     $                                  N, KL, KU, NRHS, A, LDA, XACT,
     $                                  LDB, B, LDB, ISEED, INFO )
                           XTYPE = 'C'
                           CALL DLACPY( 'Full', N, NRHS, B, LDB, BSAV,
     $                                  LDB )
*
                           IF( NOFACT .AND. ITRAN.EQ.1 ) THEN
*
*                             --- Test DGBSV  ---
*
*                             Compute the LU factorization of the matrix
*                             and solve the system.
*
                              CALL DLACPY( 'Full', KL+KU+1, N, A, LDA,
     $                                     AFB( KL+1 ), LDAFB )
                              CALL DLACPY( 'Full', N, NRHS, B, LDB, X,
     $                                     LDB )
*
                              SRNAMT = 'DGBSV '
                              CALL DGBSV( N, KL, KU, NRHS, AFB, LDAFB,
     $                                    IWORK, X, LDB, INFO )
*
*                             Check error code from DGBSV .
*
                              IF( INFO.NE.IZERO )
     $                           CALL ALAERH( PATH, 'DGBSV ', INFO,
     $                                        IZERO, ' ', N, N, KL, KU,
     $                                        NRHS, IMAT, NFAIL, NERRS,
     $                                        NOUT )
*
*                             Reconstruct matrix from factors and
*                             compute residual.
*
                              CALL DGBT01( N, N, KL, KU, A, LDA, AFB,
     $                                     LDAFB, IWORK, WORK,
     $                                     RESULT( 1 ) )
                              NT = 1
                              IF( IZERO.EQ.0 ) THEN
*
*                                Compute residual of the computed
*                                solution.
*
                                 CALL DLACPY( 'Full', N, NRHS, B, LDB,
     $                                        WORK, LDB )
                                 CALL DGBT02( 'No transpose', N, N, KL,
     $                                        KU, NRHS, A, LDA, X, LDB,
     $                                        WORK, LDB, RWORK,
     $                                        RESULT( 2 ) )
*
*                                Check solution from generated exact
*                                solution.
*
                                 CALL DGET04( N, NRHS, X, LDB, XACT,
     $                                        LDB, RCONDC, RESULT( 3 ) )
                                 NT = 3
                              END IF
*
*                             Print information about the tests that did
*                             not pass the threshold.
*
                              DO 50 K = 1, NT
                                 IF( RESULT( K ).GE.THRESH ) THEN
                                    IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                                 CALL ALADHD( NOUT, PATH )
                                    WRITE( NOUT, FMT = 9997 )'DGBSV ',
     $                                 N, KL, KU, IMAT, K, RESULT( K )
                                    NFAIL = NFAIL + 1
                                 END IF
   50                         CONTINUE
                              NRUN = NRUN + NT
                           END IF
*
*                          --- Test DGBSVX ---
*
                           IF( .NOT.PREFAC )
     $                        CALL DLASET( 'Full', 2*KL+KU+1, N, ZERO,
     $                                     ZERO, AFB, LDAFB )
                           CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X,
     $                                  LDB )
                           IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN
*
*                             Equilibrate the matrix if FACT = 'F' and
*                             EQUED = 'R', 'C', or 'B'.
*
                              CALL DLAQGB( N, N, KL, KU, A, LDA, S,
     $                                     S( N+1 ), ROWCND, COLCND,
     $                                     AMAX, EQUED )
                           END IF
*
*                          Solve the system and compute the condition
*                          number and error bounds using DGBSVX.
*
                           SRNAMT = 'DGBSVX'
                           CALL DGBSVX( FACT, TRANS, N, KL, KU, NRHS, A,
     $                                  LDA, AFB, LDAFB, IWORK, EQUED,
     $                                  S, S( N+1 ), B, LDB, X, LDB,
     $                                  RCOND, RWORK, RWORK( NRHS+1 ),
     $                                  WORK, IWORK( N+1 ), INFO )
*
*                          Check the error code from DGBSVX.
*
                           IF( INFO.NE.IZERO )
     $                        CALL ALAERH( PATH, 'DGBSVX', INFO, IZERO,
     $                                     FACT // TRANS, N, N, KL, KU,
     $                                     NRHS, IMAT, NFAIL, NERRS,
     $                                     NOUT )
*
*                          Compare WORK(1) from DGBSVX with the computed
*                          reciprocal pivot growth factor RPVGRW
*
                           IF( INFO.NE.0 ) THEN
                              ANRMPV = ZERO
                              DO 70 J = 1, INFO
                                 DO 60 I = MAX( KU+2-J, 1 ),
     $                                   MIN( N+KU+1-J, KL+KU+1 )
                                    ANRMPV = MAX( ANRMPV,
     $                                       ABS( A( I+( J-1 )*LDA ) ) )
   60                            CONTINUE
   70                         CONTINUE
                              RPVGRW = DLANTB( 'M', 'U', 'N', INFO,
     $                                 MIN( INFO-1, KL+KU ),
     $                                 AFB( MAX( 1, KL+KU+2-INFO ) ),
     $                                 LDAFB, WORK )
                              IF( RPVGRW.EQ.ZERO ) THEN
                                 RPVGRW = ONE
                              ELSE
                                 RPVGRW = ANRMPV / RPVGRW
                              END IF
                           ELSE
                              RPVGRW = DLANTB( 'M', 'U', 'N', N, KL+KU,
     $                                 AFB, LDAFB, WORK )
                              IF( RPVGRW.EQ.ZERO ) THEN
                                 RPVGRW = ONE
                              ELSE
                                 RPVGRW = DLANGB( 'M', N, KL, KU, A,
     $                                    LDA, WORK ) / RPVGRW
                              END IF
                           END IF
                           RESULT( 7 ) = ABS( RPVGRW-WORK( 1 ) ) /
     $                                   MAX( WORK( 1 ), RPVGRW ) /
     $                                   DLAMCH( 'E' )
*
                           IF( .NOT.PREFAC ) THEN
*
*                             Reconstruct matrix from factors and
*                             compute residual.
*
                              CALL DGBT01( N, N, KL, KU, A, LDA, AFB,
     $                                     LDAFB, IWORK, WORK,
     $                                     RESULT( 1 ) )
                              K1 = 1
                           ELSE
                              K1 = 2
                           END IF
*
                           IF( INFO.EQ.0 ) THEN
                              TRFCON = .FALSE.
*
*                             Compute residual of the computed solution.
*
                              CALL DLACPY( 'Full', N, NRHS, BSAV, LDB,
     $                                     WORK, LDB )
                              CALL DGBT02( TRANS, N, N, KL, KU, NRHS,
     $                                     ASAV, LDA, X, LDB, WORK, LDB,
     $                                     RWORK( 2*NRHS+1 ),
     $                                     RESULT( 2 ) )
*
*                             Check solution from generated exact
*                             solution.
*
                              IF( NOFACT .OR. ( PREFAC .AND.
     $                            LSAME( EQUED, 'N' ) ) ) THEN
                                 CALL DGET04( N, NRHS, X, LDB, XACT,
     $                                        LDB, RCONDC, RESULT( 3 ) )
                              ELSE
                                 IF( ITRAN.EQ.1 ) THEN
                                    ROLDC = ROLDO
                                 ELSE
                                    ROLDC = ROLDI
                                 END IF
                                 CALL DGET04( N, NRHS, X, LDB, XACT,
     $                                        LDB, ROLDC, RESULT( 3 ) )
                              END IF
*
*                             Check the error bounds from iterative
*                             refinement.
*
                              CALL DGBT05( TRANS, N, KL, KU, NRHS, ASAV,
     $                                     LDA, B, LDB, X, LDB, XACT,
     $                                     LDB, RWORK, RWORK( NRHS+1 ),
     $                                     RESULT( 4 ) )
                           ELSE
                              TRFCON = .TRUE.
                           END IF
*
*                          Compare RCOND from DGBSVX with the computed
*                          value in RCONDC.
*
                           RESULT( 6 ) = DGET06( RCOND, RCONDC )
*
*                          Print information about the tests that did
*                          not pass the threshold.
*
                           IF( .NOT.TRFCON ) THEN
                              DO 80 K = K1, NTESTS
                                 IF( RESULT( K ).GE.THRESH ) THEN
                                    IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                                 CALL ALADHD( NOUT, PATH )
                                    IF( PREFAC ) THEN
                                       WRITE( NOUT, FMT = 9995 )
     $                                    'DGBSVX', FACT, TRANS, N, KL,
     $                                    KU, EQUED, IMAT, K,
     $                                    RESULT( K )
                                    ELSE
                                       WRITE( NOUT, FMT = 9996 )
     $                                    'DGBSVX', FACT, TRANS, N, KL,
     $                                    KU, IMAT, K, RESULT( K )
                                    END IF
                                    NFAIL = NFAIL + 1
                                 END IF
   80                         CONTINUE
                              NRUN = NRUN + 7 - K1
                           ELSE
                              IF( RESULT( 1 ).GE.THRESH .AND. .NOT.
     $                            PREFAC ) THEN
                                 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                              CALL ALADHD( NOUT, PATH )
                                 IF( PREFAC ) THEN
                                    WRITE( NOUT, FMT = 9995 )'DGBSVX',
     $                                 FACT, TRANS, N, KL, KU, EQUED,
     $                                 IMAT, 1, RESULT( 1 )
                                 ELSE
                                    WRITE( NOUT, FMT = 9996 )'DGBSVX',
     $                                 FACT, TRANS, N, KL, KU, IMAT, 1,
     $                                 RESULT( 1 )
                                 END IF
                                 NFAIL = NFAIL + 1
                                 NRUN = NRUN + 1
                              END IF
                              IF( RESULT( 6 ).GE.THRESH ) THEN
                                 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                              CALL ALADHD( NOUT, PATH )
                                 IF( PREFAC ) THEN
                                    WRITE( NOUT, FMT = 9995 )'DGBSVX',
     $                                 FACT, TRANS, N, KL, KU, EQUED,
     $                                 IMAT, 6, RESULT( 6 )
                                 ELSE
                                    WRITE( NOUT, FMT = 9996 )'DGBSVX',
     $                                 FACT, TRANS, N, KL, KU, IMAT, 6,
     $                                 RESULT( 6 )
                                 END IF
                                 NFAIL = NFAIL + 1
                                 NRUN = NRUN + 1
                              END IF
                              IF( RESULT( 7 ).GE.THRESH ) THEN
                                 IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                              CALL ALADHD( NOUT, PATH )
                                 IF( PREFAC ) THEN
                                    WRITE( NOUT, FMT = 9995 )'DGBSVX',
     $                                 FACT, TRANS, N, KL, KU, EQUED,
     $                                 IMAT, 7, RESULT( 7 )
                                 ELSE
                                    WRITE( NOUT, FMT = 9996 )'DGBSVX',
     $                                 FACT, TRANS, N, KL, KU, IMAT, 7,
     $                                 RESULT( 7 )
                                 END IF
                                 NFAIL = NFAIL + 1
                                 NRUN = NRUN + 1
                              END IF
*
                           END IF
*
*                    --- Test DGBSVXX ---
*
*                    Restore the matrices A and B.
*
                     CALL DLACPY( 'Full', KL+KU+1, N, ASAV, LDA, A,
     $                          LDA )
                     CALL DLACPY( 'Full', N, NRHS, BSAV, LDB, B, LDB )

                     IF( .NOT.PREFAC )
     $                  CALL DLASET( 'Full', 2*KL+KU+1, N, ZERO, ZERO,
     $                    AFB, LDAFB )
                     CALL DLASET( 'Full', N, NRHS, ZERO, ZERO, X, LDB )
                     IF( IEQUED.GT.1 .AND. N.GT.0 ) THEN
*
*                       Equilibrate the matrix if FACT = 'F' and
*                       EQUED = 'R', 'C', or 'B'.
*
                        CALL DLAQGB( N, N, KL, KU, A, LDA, S, S( N+1 ),
     $                       ROWCND, COLCND, AMAX, EQUED )
                     END IF
*
*                    Solve the system and compute the condition number
*                    and error bounds using DGBSVXX.
*
                     SRNAMT = 'DGBSVXX'
                     N_ERR_BNDS = 3
                     CALL DGBSVXX( FACT, TRANS, N, KL, KU, NRHS, A, LDA,
     $                    AFB, LDAFB, IWORK, EQUED, S, S( N+1 ), B, LDB,
     $                    X, LDB, RCOND, RPVGRW_SVXX, BERR, N_ERR_BNDS,
     $                    ERRBNDS_N, ERRBNDS_C, 0, ZERO, WORK,
     $                    IWORK( N+1 ), INFO )
*
*                    Check the error code from DGBSVXX.
*
                     IF( INFO.EQ.N+1 ) GOTO 90
                     IF( INFO.NE.IZERO ) THEN
                        CALL ALAERH( PATH, 'DGBSVXX', INFO, IZERO,
     $                               FACT // TRANS, N, N, -1, -1, NRHS,
     $                               IMAT, NFAIL, NERRS, NOUT )
                        GOTO 90
                     END IF
*
*                    Compare rpvgrw_svxx from DGBSVXX with the computed
*                    reciprocal pivot growth factor RPVGRW
*

                     IF ( INFO .GT. 0 .AND. INFO .LT. N+1 ) THEN
                        RPVGRW = DLA_GBRPVGRW(N, KL, KU, INFO, A, LDA,
     $                       AFB, LDAFB)
                     ELSE
                        RPVGRW = DLA_GBRPVGRW(N, KL, KU, N, A, LDA,
     $                       AFB, LDAFB)
                     ENDIF

                     RESULT( 7 ) = ABS( RPVGRW-rpvgrw_svxx ) /
     $                             MAX( rpvgrw_svxx, RPVGRW ) /
     $                             DLAMCH( 'E' )
*
                     IF( .NOT.PREFAC ) THEN
*
*                       Reconstruct matrix from factors and compute
*                       residual.
*
                        CALL DGBT01( N, N, KL, KU, A, LDA, AFB, LDAFB,
     $                       IWORK, WORK, RESULT( 1 ) )
                        K1 = 1
                     ELSE
                        K1 = 2
                     END IF
*
                     IF( INFO.EQ.0 ) THEN
                        TRFCON = .FALSE.
*
*                       Compute residual of the computed solution.
*
                        CALL DLACPY( 'Full', N, NRHS, BSAV, LDB, WORK,
     $                               LDB )
                        CALL DGBT02( TRANS, N, N, KL, KU, NRHS, ASAV,
     $                               LDA, X, LDB, WORK, LDB, RWORK,
     $                               RESULT( 2 ) )
*
*                       Check solution from generated exact solution.
*
                        IF( NOFACT .OR. ( PREFAC .AND. LSAME( EQUED,
     $                      'N' ) ) ) THEN
                           CALL DGET04( N, NRHS, X, LDB, XACT, LDB,
     $                                  RCONDC, RESULT( 3 ) )
                        ELSE
                           IF( ITRAN.EQ.1 ) THEN
                              ROLDC = ROLDO
                           ELSE
                              ROLDC = ROLDI
                           END IF
                           CALL DGET04( N, NRHS, X, LDB, XACT, LDB,
     $                                  ROLDC, RESULT( 3 ) )
                        END IF
                     ELSE
                        TRFCON = .TRUE.
                     END IF
*
*                    Compare RCOND from DGBSVXX with the computed value
*                    in RCONDC.
*
                     RESULT( 6 ) = DGET06( RCOND, RCONDC )
*
*                    Print information about the tests that did not pass
*                    the threshold.
*
                     IF( .NOT.TRFCON ) THEN
                        DO 45 K = K1, NTESTS
                           IF( RESULT( K ).GE.THRESH ) THEN
                              IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                           CALL ALADHD( NOUT, PATH )
                              IF( PREFAC ) THEN
                                 WRITE( NOUT, FMT = 9995 )'DGBSVXX',
     $                                FACT, TRANS, N, KL, KU, EQUED,
     $                                IMAT, K, RESULT( K )
                              ELSE
                                 WRITE( NOUT, FMT = 9996 )'DGBSVXX',
     $                                FACT, TRANS, N, KL, KU, IMAT, K,
     $                                RESULT( K )
                              END IF
                              NFAIL = NFAIL + 1
                           END IF
 45                     CONTINUE
                        NRUN = NRUN + 7 - K1
                     ELSE
                        IF( RESULT( 1 ).GE.THRESH .AND. .NOT.PREFAC )
     $                       THEN
                           IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                        CALL ALADHD( NOUT, PATH )
                           IF( PREFAC ) THEN
                              WRITE( NOUT, FMT = 9995 )'DGBSVXX', FACT,
     $                             TRANS, N, KL, KU, EQUED, IMAT, 1,
     $                             RESULT( 1 )
                           ELSE
                              WRITE( NOUT, FMT = 9996 )'DGBSVXX', FACT,
     $                             TRANS, N, KL, KU, IMAT, 1,
     $                             RESULT( 1 )
                           END IF
                           NFAIL = NFAIL + 1
                           NRUN = NRUN + 1
                        END IF
                        IF( RESULT( 6 ).GE.THRESH ) THEN
                           IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                        CALL ALADHD( NOUT, PATH )
                           IF( PREFAC ) THEN
                              WRITE( NOUT, FMT = 9995 )'DGBSVXX', FACT,
     $                             TRANS, N, KL, KU, EQUED, IMAT, 6,
     $                             RESULT( 6 )
                           ELSE
                              WRITE( NOUT, FMT = 9996 )'DGBSVXX', FACT,
     $                             TRANS, N, KL, KU, IMAT, 6,
     $                             RESULT( 6 )
                           END IF
                           NFAIL = NFAIL + 1
                           NRUN = NRUN + 1
                        END IF
                        IF( RESULT( 7 ).GE.THRESH ) THEN
                           IF( NFAIL.EQ.0 .AND. NERRS.EQ.0 )
     $                        CALL ALADHD( NOUT, PATH )
                           IF( PREFAC ) THEN
                              WRITE( NOUT, FMT = 9995 )'DGBSVXX', FACT,
     $                             TRANS, N, KL, KU, EQUED, IMAT, 7,
     $                             RESULT( 7 )
                           ELSE
                              WRITE( NOUT, FMT = 9996 )'DGBSVXX', FACT,
     $                             TRANS, N, KL, KU, IMAT, 7,
     $                             RESULT( 7 )
                           END IF
                           NFAIL = NFAIL + 1
                           NRUN = NRUN + 1
                        END IF
*
                     END IF
   90                   CONTINUE
  100                CONTINUE
  110             CONTINUE
  120          CONTINUE
  130       CONTINUE
  140    CONTINUE
  150 CONTINUE
*
*     Print a summary of the results.
*
      CALL ALASVM( PATH, NOUT, NFAIL, NRUN, NERRS )

*     Test Error Bounds from DGBSVXX

      CALL DEBCHVXX(THRESH, PATH)

 9999 FORMAT( ' *** In DDRVGB, LA=', I5, ' is too small for N=', I5,
     $      ', KU=', I5, ', KL=', I5, / ' ==> Increase LA to at least ',
     $      I5 )
 9998 FORMAT( ' *** In DDRVGB, LAFB=', I5, ' is too small for N=', I5,
     $      ', KU=', I5, ', KL=', I5, /
     $      ' ==> Increase LAFB to at least ', I5 )
 9997 FORMAT( 1X, A, ', N=', I5, ', KL=', I5, ', KU=', I5, ', type ',
     $      I1, ', test(', I1, ')=', G12.5 )
 9996 FORMAT( 1X, A, '( ''', A1, ''',''', A1, ''',', I5, ',', I5, ',',
     $      I5, ',...), type ', I1, ', test(', I1, ')=', G12.5 )
 9995 FORMAT( 1X, A, '( ''', A1, ''',''', A1, ''',', I5, ',', I5, ',',
     $      I5, ',...), EQUED=''', A1, ''', type ', I1, ', test(', I1,
     $      ')=', G12.5 )
*
      RETURN
*
*     End of DDRVGBX
*
      END

max_stars_repo_name: quellyn/lapack

max_issues_repo_name: quellyn/lapack

max_forks_repo_name: quellyn/lapack


text:
module qutraj_linalg
  !
  ! Dummy module
  !

  use qutraj_precision
  use qutraj_general
  use qutraj_hilbert

  implicit none

  contains

  subroutine eigenvalues(rho,eig,n)
    ! Eigenvalues of dense hermitian matrix rho
    complex(wp), intent(in) :: rho(:,:)
    integer, intent(in) :: n
    real(wp), intent(out) :: eig(n)
    call error('eigenvalues: subroutine not available.')
  end subroutine

  subroutine entropy(rho,S)
    ! Calculate entropy for dense density matrix
    complex(wp), intent(in) :: rho(:,:)
    real(wp), intent(out) :: S
    call error('eigenvalues: subroutine not available.')
  end subroutine

end module

max_stars_repo_name: anubhavvardhan/qutip

max_issues_repo_name: trxw/qutip

max_forks_repo_name: trxw/qutip


text:
! Copyright (c) 1998-2004, University of California, Los Angeles, Emily A. Carter
!                       2004-2016, Princeton University, Emily A. Carter
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this list of
!     conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice, this list
!     of conditions and the following disclaimer in the documentation and/or other
!     materials provided with the distribution.
!
! 3. Neither the name of the copyright holder nor the names of its contributors may be
!     used to endorse or promote products derived from this software without specific
!     prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
! CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
! INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
! NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
! STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
! ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
! OF THE POSSIBILITY OF SUCH DAMAGE.
!**************************************************************
  !//  CI_UTILITIES_MOD - THIS MODULE CONTAINS A NUMBER OF UTILITY
  !//  FUNCTIONS WHICH ARE USED ONLY FOR CI CALCULATIONS
  !//
  !//  WRITTEN BY DEREK WALTER, 1999
  !//  WARNING: THIS CODE DOES NOT CONFORM TO Y2K STANDARDS!!
!**************************************************************

module ci_utilities_mod

  use graph_var_mod
  use molecule_var_mod
  use utilities_mod
  use global_var_mod

  contains

function vertex(i,j)

      !// THIS FUNCTIONS RETURNS THE VERTEX NUMBER OF THE VERTEX WITH J ELECTRONS
      !// IN THE ITH ROW

      integer::vertex
      integer, intent(in)::i,j

      vertex = auxiliary(i) + j

end function vertex

function belong(i,j)

      !// THIS FUNCTION CHECKS TO SEE IF THE VERTEX HAVING J ELECTRONS IN
      !// ROW I IS IN THE GRAPH

      logical::belong
      integer, intent(in)::i,j

      if ((i<0).or.(j<0).or.(i > num_orbitals)) then
          belong =.false.
          return
      endif

      belong = ((j >= left_border(i)).and.(j <= right_border(i)).and.&
                (i <= num_orbitals).and.(i >= 0))

end function belong

function add0(i,j)

      !// THIS FUNCTION CHECKS TO SEE IF VERTEX(I,J) HAS A 0 ARC HANGING OFF IT

      logical::add0
      integer, intent(in)::i,j

      if ((i<0).or.(j<0)) then
      add0 =.false.
      return
      endif

      add0 = (belong(i,j).and.(y0(vertex(i,j)) > 0))

end function add0

function add1(i,j)

      !// THIS FUNCTION CHECKS TO SEE IF VERTEX(I,J) HAS A 1 ARC HANGING OFF IT

      logical::add1
      integer, intent(in)::i,j

      if ((i<0).or.(j<0)) then
      add1 =.false.
      return
      endif
      add1 = (belong(i,j).and.(y1(vertex(i,j)) >= 0))

end function add1

function add2(i,j)

      !// THIS FUNCTION CHECKS TO SEE IF VERTEX(I,J) HAS A 2 ARC HANGING OFF IT

      logical::add2
      integer, intent(in)::i,j

      if ((i<0).or.(j<0)) then
      add2 =.false.
      return
      endif
      add2 = (belong(i,j).and.(y2(vertex(i,j)) >= 0))

end function add2

function num_paths(vertex)

      implicit none

      integer::allocatestatus      !// FOR DYNAMIC MEMORY
      integer::deallocatestatus    !// ALSO FOR DYNAMIC MEMORY
      integer::num_paths           !// THE PATHS WE ARE FINDING
      integer::level               !// THE CURRENT LEVEL
      integer::step_type           !// THE TYPE OF STEP WE TAKE
      integer::path_elecs          !// NUMBER OF ELECTRONS IN PATH
      integer, intent(in)::vertex  !// THE CURRENT TOP VERTEX
      integer, dimension(:), allocatable::step_array   !// KEEPS TRACK OF PATH WE MAKE

      allocate(step_array(opt_level:num_orbitals),stat=allocatestatus)
      call allocatecheck(allocatestatus,"step_arra")

      level = opt_level
      step_type = 2
      num_paths = 0
      path_elecs = vertex - auxiliary(opt_level)
      step_array = 0

      repeater: do

      !// FIRST TEST TO EXIT
      if ((level==opt_level).and.(step_type < 0)) exit repeater

      !// TEST IF WE HAVE A NEW PATH
      if (level == num_orbitals) then
          num_paths = num_paths + 1
          level = level - 1
          path_elecs = path_elecs - step_array(level)
          step_type = step_array(level) - 1
          cycle repeater

      !// NOW TEST TO SEE IF WE HAVE RUN OUT OF STEPS
      elseif (step_type < 0) then
          level = level - 1
          path_elecs = path_elecs - step_array(level)
          step_type = step_array(level) - 1
          cycle repeater

      !// NOW, GO THROUGH THE STEP TYPES.  FIRST 2, THEN 1, THEN 0
      elseif (step_type == 2) then
          if (add2(level,path_elecs)) then  !// THE STEP IS ALLOWED
          step_array(level) = 2
          level = level + 1
          path_elecs = path_elecs + 2
          step_type = 2
          cycle repeater
          else   !// THE STEP IS NOT ALLOWED, TRY A 1 STEP FROM THE SAME LEVEL
          step_type = 1
          cycle repeater
          endif

      !// NOW A 1 STEP
      elseif (step_type == 1) then
          if (add1(level,path_elecs)) then  !// THE STEP IS ALLOWED
          step_array(level) = 1
          level = level + 1
          path_elecs = path_elecs + 1
          step_type = 2
          cycle repeater
          else   !// THE STEP IS NOT ALLOWED, TRY A 0 STEP FROM THE SAME LEVEL
          step_type = 0
          cycle repeater

          endif

      !// NOW A 0 STEP
      elseif (step_type == 0) then
          if (add0(level,path_elecs)) then  !// THE STEP IS ALLOWED
          step_array(level) =  0
          level = level + 1
          step_type = 2
          path_elecs = path_elecs
          cycle repeater
          else   !// THE STEP IS NOT ALLOWED; STEP BY A LEVEL
          step_type = - 1
          cycle repeater
          endif

      endif

      enddo repeater

      deallocate(step_array, stat=deallocatestatus)
      call deallocatecheck(deallocatestatus,"step_arr")

end function num_paths

function num_paths2(begin_lev,begin_elec,end_lev,end_elec)

      !// IN THIS FUNCTION WE COUNT THE NUMBER OF PATHS BETWEEN
      !// TWO VERTICES

      implicit none

      integer::allocatestatus         !// FOR DYNAMIC MEMORY
      integer::deallocatestatus       !// ALSO FOR DYNAMIC MEMORY
      integer::num_paths2             !// THE PATHS WE ARE FINDING
      integer::level                  !// THE CURRENT LEVEL
      integer::step_type              !// THE TYPE OF STEP WE TAKE
      integer::path_elecs             !// NUMBER OF ELECTRONS IN PATH
      integer::num_needed             !// NUMBER OF ELECTRONS WE NEED
      integer::num_we_can_add         !// NUMBER OF ELECTRONS WE HAVE
      integer, intent(in)::begin_lev  !// BEGINNING LEVEL
      integer, intent(in)::end_lev    !// ENDING LEVEL
      integer, intent(in)::begin_elec !// BEGINNING NUMBER OF VERTICES
      integer, intent(in)::end_elec   !// ENDING NUMBER OF VERTICES
      integer, dimension(:), allocatable::step_array   !// KEEPS TRACK OF PATH WE MAKE

      allocate(step_array(begin_lev:end_lev),stat=allocatestatus)
      call allocatecheck(allocatestatus,"step_arra")

      level = begin_lev
      step_type = 2
      num_paths2 = 0
      path_elecs = begin_elec
      step_array = 0

      if (end_elec < begin_elec) return
      if (level == end_lev.and.end_elec /= begin_elec) return

      repeater: do

      !// FIRST TEST TO EXIT
      if ((level==begin_lev).and.(step_type < 0)) exit repeater

      !// TEST IF WE HAVE A NEW PATH
      if (level == end_lev) then

          if (path_elecs==end_elec) num_paths2 = num_paths2 + 1
          level = level - 1
          path_elecs = path_elecs - step_array(level)
          step_type = step_array(level) - 1
          cycle repeater

      !// NOW TEST TO SEE IF WE HAVE RUN OUT OF STEPS
      elseif (step_type < 0) then
          level = level - 1
          path_elecs = path_elecs - step_array(level)
          step_type = step_array(level) - 1
          cycle repeater

      !// NOW, GO THROUGH THE STEP TYPES.  FIRST 2, THEN 1, THEN 0
      elseif (step_type == 2) then
          if (add2(level,path_elecs).and.&
              ((path_elecs+2) <= end_elec)) then  !// THE STEP IS ALLOWED
          step_array(level) = 2
          level = level + 1
          path_elecs = path_elecs + 2
          step_type = 2
          cycle repeater
          else   !// THE STEP IS NOT ALLOWED, TRY A 1 STEP FROM THE SAME LEVEL
          step_type = 1
          cycle repeater
          endif

      !// NOW A 1 STEP
      elseif (step_type == 1) then
          num_needed = end_elec-path_elecs-1
          num_we_can_add = 2*(end_lev - level - 1)
          if (add1(level,path_elecs).and.&
              ((path_elecs+1) <= end_elec).and.&
              (num_needed <= num_we_can_add)) then  !// THE STEP IS ALLOWED
          step_array(level) = 1
          level = level + 1
          path_elecs = path_elecs + 1
          step_type = 2
          cycle repeater
          else   !// THE STEP IS NOT ALLOWED, TRY A 0 STEP FROM THE SAME LEVEL
          step_type = 0
          cycle repeater

          endif

      !// NOW A 0 STEP
      elseif (step_type == 0) then
          num_needed = end_elec-path_elecs
          num_we_can_add = 2*(end_lev - level - 1)
          if (add0(level,path_elecs).and.&
              (num_needed <= num_we_can_add)) then  !// THE STEP IS ALLOWED
          step_array(level) =  0
          level = level + 1
          step_type = 2
          path_elecs = path_elecs
          cycle repeater
          else   !// THE STEP IS NOT ALLOWED; STEP BY A LEVEL
          step_type = - 1
          cycle repeater
          endif

      endif

      enddo repeater

      deallocate(step_array, stat=deallocatestatus)
      call deallocatecheck(deallocatestatus,"step_arr")

end function num_paths2

function maximum_singles(levels, number_electrons)

      !// THIS FUNCTION TELLS YOU THE MAXIMUM NUMBER A SINGLES
      !// THAT A PATH WITH A CERTAIN NUMBER OF LEVELS AND A
      !// A CERTAIN NUMBER OF ELECTRONS CAN HAVE
      integer::maximum_singles
      integer, intent(in)::levels, number_electrons

      if (levels >= number_electrons) then
          maximum_singles = number_electrons
      elseif (number_electrons <= 2*levels) then
          maximum_singles = levels - (number_electrons - levels)
      else
          write(ioError,450) number_electrons,levels
          450 format (1x,"You have too many electrons and they won't fit!!",/,&
              1x,"Number electrons =",i4,/,&
              1x,"Number of levels =",i4)
                  stop
      endif

end function maximum_singles

function E_matrix(size,zopen_shells)

      integer,intent(in)::size
      integer,intent(in)::zopen_shells
      real(real8), dimension(size,size)::E_matrix

      integer::j  !// LOOP CONTROL VARIABLE
      integer::f2 !//  FSN(OPEN_SHELLS) - FSN(OPEN_SHELLS-2)

      f2 = fsn(zopen_shells-2)

      E_matrix = real(0.0,real8)

      do j = 1, f2
          E_matrix(j,j) = real(1.0,real8)
      enddo

      do j = f2+1,fsn(zopen_shells)
          E_matrix(j,j) = -real(1.0,real8)
      enddo

end function E_matrix

function E_matrix_as_vector(zopen_shells)

      integer,intent(in)::zopen_shells
      real(real8), dimension(fsn(zopen_shells))::E_matrix_as_vector

      integer::j
      integer::f2

      f2 = fsn(zopen_shells-2)

      do j = 1, f2
          E_matrix_as_vector(j) = real(1.0,real8)
      enddo

      do j = f2+1,fsn(zopen_shells)
          E_matrix_as_vector(j) = -real(1.0,real8)
      enddo

end function E_matrix_as_vector

!> \brief Given a set of orbital occupations, test if they match a reference
!> \param occupations   The internal orbital occupations for this configuration
!>
!> NOTE! If you pass part of a path (path%occupations) make sure you make the argument path%occupations(1:num_internal). A 0 element
!> will screw this up.
logical function is_reference(occupations)
implicit none
integer, dimension(:) :: occupations
integer :: iRef
is_reference = .false.
do iRef = 1, num_ref
    if ( arraysAreEqual(occupations(1:num_internal), references(1:num_internal,iRef))) then
        is_reference = .true.
        return
    end if
end do
end function is_reference

end module ci_utilities_mod

max_stars_repo_name: EACcodes/TigerCI

max_issues_repo_name: EACcodes/TigerCI

max_forks_repo_name: EACcodes/TigerCI


text:
!==============================================================================!
  subroutine Control_Mod_Point_For_Monitoring_Planes(b_xp, b_yp, b_zp,  &
                                                     verbose)
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  real              :: b_xp, b_yp, b_zp
  logical, optional :: verbose
!-----------------------------------[Locals]-----------------------------------!
  real :: def(3)
  real :: val(3)
!==============================================================================!

  data def / 0.0, 0.0, 0.0 /

  call Control_Mod_Read_Real_Array('POINT_FOR_MONITORING_PLANES', 3, def,  &
                                    val, verbose)
  b_xp = val(1)
  b_yp = val(2)
  b_zp = val(3)

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
c
c   Dint – version 2.0  is licensed under the Apache License, Version 2.0 (the "License");
c   you may not use Dint – version 2.0 except in compliance with the License.
c   You may obtain a copy of the License at
c       http://www.apache.org/licenses/LICENSE-2.0
c   The license is also given in the LICENSE file.
c   Unless required by applicable law or agreed to in writing, software
c   distributed under the License is distributed on an "AS IS" BASIS,
c   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
c   See the License for the specific language governing permissions and limitations under the License.
c
c -------------------------------------------------------------------------------------------
c  Dint : Direct Nonadiabatic Trajectories A code for non-Born–Oppenheimer molecular dynamics
c
c  version 2.0
c
c  A. W. Jasper
c  Argonne National Laboratory
c
c  Rui Ming Zhang
c  Tsinghua University
c
c  and
c
c  D. G. Truhlar
c  University of Minnesota
c
c  copyright  2020
c  Donald G. Truhlar and Regents of the University of Minnesota
c----------------------------------------------------------------------------------------------

      subroutine diamin2(xmin,val,xj,mm,nsurf,symb,rguess)

c Borrowed with modifications from NAT8.1.

c MNBRAK and BRENT are taken from Numerical Recipes.
c This routine is a combination of MNBRAK and BRENT.
c It minimizes the diatomic energy (obtained by calling
c DIAPOT) with respect to bond distance.  The value
c of the miminum-energy bond distance that is returned
c is XMIN, and the potential at the value is VAL.

c NOTE:  The code will automatically check for minima
c for all surfaces and for all three molecular
c arrangements.  If one of these combinations is
c unbound, the code may return an error message
c and non-physical values.  This is OK because
c the code does not use these values for molecular
c arrangements that are not energetically accessible.

      implicit none
      include 'param.f'
c mnbrak variables
      integer im,nsurf
      double precision ax,bx,cx,fa,fb,fc,func,gold,glimit,tiny,guess
      parameter (gold=1.618034d0, glimit=100.0d0,tiny=1.d-20)
      double precision dum,fu,q,r,u,ulim,xj,rguess

c brent variables
      integer itmax
      double precision brent,tol,CGOLD,ZEPS
      parameter (itmax=100,CGOLD=.3819660,ZEPS=1.0e-10)
      integer iter
      double precision a,b,d,e,etemp,fv,fw,fx,p,mm(mnat),
     &                 tol1,tol2,v,w,x,xm,val,xmin
      character*2 symb(mnat)

c hard-coded
c      guess = 2.d0
      guess = rguess
      tol = 1.d-6

c      write(6,*)"aj in diamin2"

c ***********MNBRAK***********
c initial two points are guess, and 80% of guess
      ax = guess*0.8d0
      bx = guess

c func returns y(x), where y is the function to be minimized
      call diapot2(ax,fa,xj,mm,nsurf,symb)
      call diapot2(bx,fb,xj,mm,nsurf,symb)

c make sure that fa > fb
      if (fb.gt.fa) then
        dum=ax
        ax=bx
        bx=dum
        dum=fb
        fb=fa
        fa=dum
      endif

c get third point by golden ration
      cx=bx+gold*(bx-ax)
      call diapot2(cx,fc,xj,mm,nsurf,symb)

c bracket minimum so that ax < bx < cx and fa > fb and fc > fb
 11   if (fb.ge.fc) then
        r=(bx-ax)*(fb-fc)
        q=(bx-cx)*(fb-fa)
        u=bx-((bx-cx)*q-(bx-ax)*r)/
     &     (2.0d0*dsign(max(dabs(q-r),tiny),q-r))
        ulim=bx+glimit*(cx-bx)
        if ((bx-u)*(u-cx).gt.0.0d0) then
      call diapot2(u,fu,xj,mm,nsurf,symb)
          if (fu.lt.fc)then
            ax=bx
            fa=fb
            bx=u
            fb=fu
            return
          elseif (fu.gt.fb) then
            cx=u
            fc=fu
            return
          endif
          u=cx+gold*(cx-bx)
      call diapot2(u,fu,xj,mm,nsurf,symb)
        elseif ((cx-u)*(u-ulim).gt.0.0d0) then
      call diapot2(u,fu,xj,mm,nsurf,symb)
          if (fu.lt.fc) then
            bx=cx
            cx=u
            u=cx+gold*(cx-bx)
            fb=fc
            fc=fu
      call diapot2(u,fu,xj,mm,nsurf,symb)
          endif
        elseif ((u-ulim)*(ulim-cx).ge.0.0d0) then
          u=ulim
      call diapot2(u,fu,xj,mm,nsurf,symb)
        else
          u=cx+gold*(cx-bx)
      call diapot2(u,fu,xj,mm,nsurf,symb)
        endif
        ax=bx
        bx=cx
        cx=u
        fa=fb
        fb=fc
        fc=fu
        goto 11
      endif

c uses ax,bx,cx, and the tolerance tol
c ********** BRENT ***********

      a=min(ax,cx)
      b=max(ax,cx)
      if (cx.lt.ax) a=cx
      if (cx.lt.ax) b=ax

      v=bx
      w=v
      x=v
      e=0.d0
      call diapot2(x,fx,xj,mm,nsurf,symb)
      fv=fx
      fw=fx
      do iter=1,itmax
         xm=0.5d0*(a+b)
         tol1=tol*dabs(x)+ZEPS
         tol2=2.d0*tol1
         if(dabs(x-xm).le.(tol2-.5d0*(b-a))) go to 3
         if(dabs(e).gt.tol1) then
            r=(x-w)*(fx-fv)
            q=(x-v)*(fx-fw)
            p=(x-v)*q-(x-w)*r
            q=2.d0*(q-r)
            if(q.gt.0.d0) p=-p
            q=dabs(q)
            etemp=e
            e=d
            if(dabs(p).ge.dabs(.5d0*q*etemp).or.p.le.q*(a-x).or.
     *           p.ge.q*(b-x)) go to 1
            d=p/q
            u=x+d
            if(u-a.lt.tol2 .or. b-u.lt.tol2) d=dsign(tol1,xm-x)
            go to 2
         endif
    1    if(x.ge.xm) then
            e=a-x
         else
            e=b-x
         endif
         d=CGOLD*e
    2    if(dabs(d).ge.tol1) then
            u=x+d
         else
            u=x+dsign(tol1,d)
         endif
      call diapot2(u,fu,xj,mm,nsurf,symb)
         if (fu.le.fx) then
            if (u.ge.x) then
               a=x
            else
               b=x
            endif
            v=w
            fv=fw
            w=x
            fw=fx
            x=u
            fx=fu
         else
            if(u.lt.x) then
               a=u
            else
               b=u
            endif
            if(fu.le.fw .or. w.eq.x) then
               v=w
               fv=fw
               w=u
               fw=fu
            else if(fu.le.fv .or. v.eq.x .or. v.eq.w) then
               v=u
               fv=fu
            endif
         endif
      enddo
      write(6,*)"brent exceed maximum iterations"
      x = 1.d10
      fx = 1.d10

    3 xmin=x
      val=fx

      return
      END

ext: f

max_stars_repo_name: kevinmooreiii/DiNT

max_issues_repo_name: kevinmooreiii/DiNT

max_forks_repo_name: kevinmooreiii/DiNT


text:
!***********************************************************************
!                                                                      *
      SUBROUTINE ORTHOR(J, INV)
!
!                                                                      *
!   This routine Schmidt orthogonalizes orbital  J  to  all orbitals   *
!   which  have  better  self-consistency.  Note that fixed orbitals   *
!   have the best self-consistency.                                    *
!                                                                      *
!   Call(s) to: [LIB92]: COUNT, RINT.                                  *
!                                                                      *
!   Written by Farid A Parpia, at Oxford    Last update: 22 Dec 1992   *
!                                                                      *
!***********************************************************************
! INV not initialized here but used as
! INV = INV + 1
! A bug regarding IORDER fixed and new schemes introduced. See comments
! Anyway this routine was not used anywhere.
!XHH 1997.02.21
!************************************************************************
!...Translated by Pacific-Sierra Research 77to90  4.3E  16:54:29   1/ 6/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param,  ONLY: DOUBLE
      USE parameter_def,    ONLY: NNNP
      USE DEF_C
      USE GRID_C
      USE OVL_C
      USE ORB_C
      USE ORBA_C
      USE INVT_C
      USE SCF_C
      USE IOUNIT_C
      USE FIXD_C
      USE wave_C
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE rint_I
!      USE count_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER  :: J
      INTEGER , INTENT(INOUT) :: INV
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: KRAW, K, MTP, I, NNCFF, JGG
      REAL(DOUBLE) :: EPS, OVRLAP, DNORM, FACTOR, SGN
      LOGICAL :: CHECK, CHANGED
!-----------------------------------------------
!
!XHH Added common /orba/ and /iounit/

!
!ww      EPS = ACCY*0.1D 00
      EPS = ACCY*0.01D00
!
      CHECK = .NOT.NOINVT(J)
!
!XHH
! Bug fixed. Here J is actually IORDER(J_raw), thus K should be
! treated the same way.
      CHANGED = .FALSE.
!      DO 4 K = 1,NW
      DO KRAW = 1, NW
         K = IORDER(KRAW)
!         write(istde,*) '***',kraw,k,np(k),nh(k),scnsty(k),'***'
!
!XHH orbitals with higher self-consistency are considered
!         IF ( (K .NE. J) .AND.
!     :        (NAK(K) .EQ. NAK(J)) .AND.
!     :        (SCNSTY(K) .LT. SCNSTY(J)) ) THEN
!XHH All orbitals are considered
!         IF ( (K .NE. J) .AND.
!     :        (NAK(K) .EQ. NAK(J)) ) THEN
!XHH orbitals before the current and unchanged ones are considered
         IF (.NOT.(NAK(K)==NAK(J) .AND. (K<J .OR. LFIX(K)))) CYCLE

         CHANGED = .TRUE.
!
!   Compute overlap
!
         OVRLAP = RINT(J,K,0)
!
!   Schmidt orthogonalise
!
         PZ(J) = PZ(J) - OVRLAP*PZ(K)
         MTP = MAX(MF(J),MF(K))
         PF(:MTP,J) = PF(:MTP,J) - OVRLAP*PF(:MTP,K)
         QF(:MTP,J) = QF(:MTP,J) - OVRLAP*QF(:MTP,K)
      END DO
!XHH
      IF (CHANGED) THEN
!
!   Normalise
!
         MF(J) = MTP
         DNORM = RINT(J,J,0)
         FACTOR = 1.0D00/SQRT(DNORM)
!
!   Determine if inversion is necessary
!
         IF (CHECK) THEN
            CALL COUNT (PF(:NNNP,J), MTP, NNCFF, SGN)
            IF (SGN < 0.0D00) THEN
               INV = INV + 1
               FACTOR = -FACTOR
            ENDIF
         ENDIF
!
!   Perform normalization and/or inversion
!
         PZ(J) = FACTOR*PZ(J)
         PF(2:MTP,J) = FACTOR*PF(2:MTP,J)
         QF(2:MTP,J) = FACTOR*QF(2:MTP,J)
!
!   Find new MF(J)
!
         MTP = MTP + 1
    3    CONTINUE
         MTP = MTP - 1
         IF (ABS(PF(MTP,J)) < EPS) THEN
            PF(MTP,J) = 0.0D00
            QF(MTP,J) = 0.0D00
            GO TO 3
         ELSE
            MF(J) = MTP
         ENDIF
!
!XHH two lines moved forward and added an ENDIF
!         ENDIF
!    4 CONTINUE
      ENDIF
!
      RETURN
      END SUBROUTINE ORTHOR

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
lambda x:Bool. x;
(lambda x:Bool->Bool. if x false then true else false)
(lambda x:Bool. if x then false else true);

max_stars_repo_name: brunoflores/fulluntyped

max_issues_repo_name: brunoflores/fulluntyped

max_forks_repo_name: brunoflores/fulluntyped


text:
// Generated Files
${PROJ_DIR}/llink/script/premade_examples/axi_mm_a48_d512_packet/axi_mm_a48_d512_packet_master_top.sv
${PROJ_DIR}/llink/script/premade_examples/axi_mm_a48_d512_packet/axi_mm_a48_d512_packet_master_concat.sv
${PROJ_DIR}/llink/script/premade_examples/axi_mm_a48_d512_packet/axi_mm_a48_d512_packet_master_name.sv

// Logic Link files
-f ${PROJ_DIR}/llink/rtl/llink.f

// Common Files
-f ${PROJ_DIR}/common/rtl/common.f

max_stars_repo_name: chipsalliance/aib-protocols

max_issues_repo_name: chipsalliance/aib-protocols

max_forks_repo_name: chipsalliance/aib-protocols


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2020 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine extrapolate_u1(yi,yn,ipkon,inum,kon,lakon,nfield,nk,
     &  ne,mi,ndim,orab,ielorien,co,iorienloc,cflag,
     &  vold,force,ielmat,thicke,ielprop,prop,i)
!
!     extrapolates field values at the integration points to the
!     nodes for user element i of type u1
!
!
!     INPUT:
!
!     yi(ndim,mi(1),*)   value of the variables at the integration
!                        points
!     ipkon(i)           points to the location in field kon preceding
!                        the topology of element i
!     inum(i)            < 0: node i is a network node
!                        > 0: node i is a structural node
!                        =0: node i is not used
!     kon(*)             contains the topology of all elements. The
!                        topology of element i starts at kon(ipkon(i)+1)
!                        and continues until all nodes are covered. The
!                        number of nodes depends on the element label
!     lakon(i)           contains the label of element i
!     nfield             number of variables to be extrapolated
!     nk                 maximum node number in the mesh
!     ne                 maximum element number in the mesh
!     ndim               number of variables in the integration point
!                        field to be extrapolated
!     orab(7,*)          description of all local coordinate systems.
!                        (cf. List of variables and their meaning in the
!                        User's manual)
!     ielorien(i)        orientation in element i
!     co(1..3,i)         global coordinates of node i
!     iorienloc          0: extrapolated variables requested in global
!                           coordinates
!                        1: extrapolated variables requested in local
!                           coordinates
!     cflag (char*1)     I: interpolate 3D results onto 1D/2D
!                        E: store extrapolated 1D/2D results
!                        M: store 1D section forces
!                        blank: any other case
!     vold(j,i)          value of variable j in node i at the end
!                        of the previous iteration
!     force              logical variable; if true the values to
!                        be extrapolated are force values; important for
!                        interpolation from 3D expanded structures on the
!                        original 1D/2D structure: forces across the
!                        expansion have to be summed, not interpolated
!     ielmat(i)          material of element i
!     thicke(j,i)        thickness of layer j in node i
!     ielprop(i)         properties for element i are stored in
!                        prop(ielprop(i)+1),prop(ielprop(i)+2),....
!                        (number of properties depends on the type of
!                        element)
!     prop               property field
!     i                  number of the element for which the extrapolation
!                        is to be performed
!
!     OUTPUT:
!
!     yn(nfield,*)       value of the variables at the nodes
!
      implicit none
!
      logical force
!
      character*1 cflag
      character*8 lakon(*)
!
      integer ipkon(*),inum(*),kon(*),mi(*),ne,nfield,nk,i,ndim,
     &  iorienloc,ielorien(mi(3),*),ielmat(mi(3),*),ielprop(*)
!
      real*8 yi(ndim,mi(1),*),yn(nfield,*),orab(7,*),co(3,*),prop(*),
     &  vold(0:mi(2),*),thicke(mi(3),*)
!
!
!
!     START OF THIS SUBROUTINE
!
      integer indexe,j,k,node
!
      if(iorienloc.ne.0) then
         write(*,*) '*ERROR in extrapolate_u1'
         write(*,*) '       no local orientation for variables'
         write(*,*) '       belonging to this type of element'
         write(*,*) '       allowed'
         call exit(201)
      endif
!
      if(nfield.eq.6) then
         indexe=ipkon(i)
         do j=1,2
            node=kon(indexe+j)
            do k=1,nfield
               yn(k,node)=yn(k,node)+yi(k,1,i)
            enddo
         enddo
      else
         write(*,*) '*ERROR in extrapolate_u1'
         write(*,*) '       extropolation for element of type u1'
         write(*,*) '       is only coded for fields with 6'
         write(*,*) '       entries'
         call exit(201)
      endif
!
!     END OF THIS SUBROUTINE
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
C    @(#)chkbrn.f	20.4 11/11/97
      subroutine chkbrn

      include 'ipfinc/parametr.inc'

      include 'ipfinc/brchk.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/outpt2.inc'
      include 'ipfinc/prt.inc'
C
C     BLKDTA-CHKBRN
C
      common /itot/ itot,ioerr
      dimension ibchk(3,MAXBRN)
      equivalence (ibchk,bchk)
      save

C     STORE BRANCH QUANTITIES

      if (k1 .lt. k2) then
        kpack = 200000*k1 + 10*k2 + mod(ichar(id),10)
        itot=itot+1
        ibchk(1,itot)=kpack
        bchk(2,itot)=ploss
        bchk(3,itot)=qloss
      else
        il1 = 1
        il2 = itot
        kpack = 200000*k2 + 10*k1 + mod(ichar(id),10)

        do while (il1 .le. il2)
          il = (il1+il2) / 2
          komp = kpack - ibchk(1,il)
          if (komp .lt. 0) then
            il2 = il - 1
          else if (kpack .gt. 0) then
            il1 = il + 1
          else
            go to 160
          endif
        enddo

        write (errbuf(1),150) bus(k1), base(k1), bus(k2), base(k2),
     &    id
  150   format (' Output error -- non-symmetrical ',
     1        'branch ', a8, f7.1, ' TO ', a8, f7.1, ' PAR ', a1)
        call prterx ('W',1)
        ioerr=ioerr+1
        go to 190

  160   aloss=bchk(2,il)
        bloss=bchk(3,il)
        if (sqrt((ploss-aloss)**2+(qloss-bloss)**2) .gt. 0.05) then
          aloss=ploss-aloss
          bloss=qloss-bloss
          write (errbuf(1),170) bus(k1),base(k1),bus(k2),base(k2),
     &      id,aloss,bloss
  170     format (' Output error --- branch loss quantities',
     1        'are non-symmetrical ',a8,f7.1,' TO ',a8,f7.1,
     2        ' par ',a1,' DP =',f7.2,' DQ =',f7.2)
          call prterx ('W',1)
          ioerr=ioerr+1
        endif
      endif
  190 continue
      return

C     CHECK LOSS QUANTITIES

      entry fibrck

      if (itot.eq.0) go to 220
      pinsum=0.0
      qinsum=0.0

      do i=1,itot
        pinsum=pinsum+bchk(2,i)
        qinsum=qinsum+bchk(3,i)
      enddo

      write (outbuf,210) ioerr,pinsum,qinsum
  210 format('0 LOSS CHECK: TOTAL ERRORS ',i4,' PINSUM ',f10.3,
     1  ' QINSUN ',f10.3)
      call prtout(1)

  220 continue

      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      SUBROUTINE CNMN06 (X,VLB,VUB,G,SCAL,DF,S,G1,G2,CTAM,CTBM,SLOPE,ALP
     1,A2,A3,A4,F1,F2,F3,CV1,CV2,CV3,CV4,ALPCA,ALPFES,ALPLN,ALPMIN,ALPNC
     2,ALPSAV,ALPSID,ALPTOT,ISC,N1,N2,NCAL,NVC,ICOUNT,IGOOD1,IGOOD2,IGOO
     3D3,IGOOD4,IBEST,III,NLNC,JGOTO)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      COMMON /CNMN1/ DELFUN,DABFUN,FDCH,FDCHM,CT,CTMIN,CTL,CTLMIN,ALPHAX
     1,ABOBJ1,THETA,OBJ,NDV,NCON,NSIDE,IPRINT,NFDG,NSCAL,LINOBJ,ITMAX,IT
     2RM,ICNDIR,IGOTO,NAC,INFO,INFOG,ITER,NFEASCT
      COMMON /OUTPUT/ IOUT
      DIMENSION X(N1), VLB(N1), VUB(N1), G(N2), SCAL(N1), DF(N1), S(N1),
     1 G1(N2), G2(N2), ISC(N2), NCAL(2)
C     ROUTINE TO SOLVE ONE-DIMENSIONAL SEARCH PROBLEM FOR CONSTRAINED
C     FUNCTION MINIMIZATION.
C     BY G. N. VANDERPLAATS                           AUG., 1974.
C     NASA-AMES RESEARCH CENTER, MOFFETT FIELD, CALIF.
C     OBJ = INITIAL AND FINAL FUNCTION VALUE.
C     ALP = MOVE PARAMETER.
C     SLOPE = INITIAL SLOPE.
C
C     ALPSID = MOVE TO SIDE CONSTRAINT.
C     ALPFES = MOVE TO FEASIBLE REGION.
C     ALPNC = MOVE TO NEW NON-LINEAR CONSTRAINT.
C     ALPLN = MOVE TO LINEAR CONSTRAINT.
C     ALPCA = MOVE TO RE-ENCOUNTER CURRENTLY ACTIVE CONSTRAINT.
C     ALPMIN = MOVE TO MINIMIZE FUNCTION.
C     ALPTOT = TOTAL MOVE PARAMETER.
      ZRO=0.
      IF (JGOTO.EQ.0) GO TO 10
      GO TO (140,310,520),JGOTO
10    IF (IPRINT.GE.5) WRITE (IOUT,730)
      ALPSAV=ALP
      ICOUNT=0
      ALPTOT=0.
C     TOLERANCES.
      CTAM=ABS(CTMIN)
      CTBM=ABS(CTLMIN)
C     PROPOSED MOVE.
20    CONTINUE
C     ------------------------------------------------------------------
C     *****  BEGIN SEARCH OR IMPOSE SIDE CONSTRAINT MODIFICATION  *****
C     ------------------------------------------------------------------
      A2=ALPSAV
      ICOUNT=ICOUNT+1
      ALPSID=1.0E+20
C     INITIAL ALPHA AND OBJ.
      ALP=0.
      F1=OBJ
      KSID=0
      IF (NSIDE.EQ.0) GO TO 70
C     ------------------------------------------------------------------
C     FIND MOVE TO SIDE CONSTRAINT AND INSURE AGAINST VIOLATION OF
C     SIDE CONSTRAINTS
C     ------------------------------------------------------------------
      DO 60 I=1,NDV
      SI=S(I)
      IF (ABS(SI).GT.1.0E-20) GO TO 30
C     ITH COMPONENT OF S IS SMALL.  SET TO ZERO.
      S(I)=0.
      SLOPE=SLOPE-SI*DF(I)
      GO TO 60
30    CONTINUE
      XI=X(I)
      SI=1./SI
      IF (SI.GT.0.) GO TO 40
C     LOWER BOUND.
      XI2=VLB(I)
      XI1=ABS(XI2)
      IF (XI1.LT.1.) XI1=1.
C     CONSTRAINT VALUE.
      GI=(XI2-XI)/XI1
      IF (GI.GT.-1.0E-6) GO TO 50
C     PROPOSED MOVE TO LOWER BOUND.
      ALPA=(XI2-XI)*SI
      IF (ALPA.LT.ALPSID) ALPSID=ALPA
      GO TO 60
40    CONTINUE
C     UPPER BOUND.
      XI2=VUB(I)
      XI1=ABS(XI2)
      IF (XI1.LT.1.) XI1=1.
C     CONSTRAINT VALUE.
      GI=(XI-XI2)/XI1
      IF (GI.GT.-1.0E-6) GO TO 50
C     PROPOSED MOVE TO UPPER BOUND.
      ALPA=(XI2-XI)*SI
      IF (ALPA.LT.ALPSID) ALPSID=ALPA
      GO TO 60
50    CONTINUE
C     MOVE WILL VIOLATE SIDE CONSTRAINT.  SET S(I)=0.
      SLOPE=SLOPE-S(I)*DF(I)
      S(I)=0.
      KSID=KSID+1
60    CONTINUE
C     ALPSID IS UPPER BOUND ON ALPHA.
      IF (A2.GT.ALPSID) A2=ALPSID
70    CONTINUE
C     ------------------------------------------------------------------
C               CHECK ILL-CONDITIONING
C     ------------------------------------------------------------------
      IF (KSID.EQ.NDV.OR.ICOUNT.GT.10) GO TO 710
      IF (NVC.EQ.0.AND.SLOPE.GT.0.) GO TO 710
      ALPFES=-1.
      ALPMIN=-1.
      ALPLN=1.1*ALPSID
      ALPNC=ALPSID
      ALPCA=ALPSID
      IF (NCON.EQ.0) GO TO 90
C     STORE CONSTRAINT VALUES IN G1.
      DO 80 I=1,NCON
      G1(I)=G(I)
80    CONTINUE
90    CONTINUE
C     ------------------------------------------------------------------
C                  MOVE A DISTANCE A2*S
C     ------------------------------------------------------------------
      ALPTOT=ALPTOT+A2
      DO 100 I=1,NDV
      X(I)=X(I)+A2*S(I)
100   CONTINUE
      IF (IPRINT.LT.5) GO TO 130
      WRITE (IOUT,740) A2
      IF (NSCAL.EQ.0) GO TO 120
      DO 110 I=1,NDV
110   G(I)=SCAL(I)*X(I)
      WRITE (IOUT,750) (G(I),I=1,NDV)
      GO TO 130
120   WRITE (IOUT,750) (X(I),I=1,NDV)
C     ------------------------------------------------------------------
C                   UPDATE FUNCTION AND CONSTRAINT VALUES
C     ------------------------------------------------------------------
130   NCAL(1)=NCAL(1)+1
      JGOTO=1
      RETURN
140   CONTINUE
      F2=OBJ
      IF (IPRINT.GE.5) WRITE (IOUT,760) F2
      IF (IPRINT.LT.5.OR.NCON.EQ.0) GO TO 150
      WRITE (IOUT,770)
      WRITE (IOUT,750) (G(I),I=1,NCON)
150   CONTINUE
C     ------------------------------------------------------------------
C               IDENTIFY ACCAPTABILITY OF DESIGNS F1 AND F2
C     ------------------------------------------------------------------
C     IGOOD = 0 IS ACCAPTABLE.
C     CV = MAXIMUM CONSTRAINT VIOLATION.
      IGOOD1=0
      IGOOD2=0
      CV1=0.
      CV2=0.
      NVC1=0
      IF (NCON.EQ.0) GO TO 170
      DO 160 I=1,NCON
      CC=CTAM
      IF (ISC(I).GT.0) CC=CTBM
      C1=G1(I)-CC
      C2=G(I)-CC
      IF (C2.GT.0.) NVC1=NVC1+1
      IF (C1.GT.CV1) CV1=C1
      IF (C2.GT.CV2) CV2=C2
160   CONTINUE
      IF (CV1.GT.0.) IGOOD1=1
      IF (CV2.GT.0.) IGOOD2=1
170   CONTINUE
      ALP=A2
      OBJ=F2
C     ------------------------------------------------------------------
C     IF F2 VIOLATES FEWER CONSTRAINTS THAN F1 BUT STILL HAS CONSTRAINT
C     VIOLATIONS RETURN
C     ------------------------------------------------------------------
      IF (NVC1.LT.NVC.AND.NVC1.GT.0) GO TO 710
C     ------------------------------------------------------------------
C             IDENTIFY BEST OF DESIGNS F1 ANF F2
C     ------------------------------------------------------------------
C     IBEST CORRESPONDS TO MINIMUM VALUE DESIGN.
C     IF CONSTRAINTS ARE VIOLATED, IBEST CORRESPONDS TO MINIMUM
C     CONSTRAINT VIOLATION.
      IF (IGOOD1.EQ.0.AND.IGOOD2.EQ.0) GO TO 180
C     VIOLATED CONSTRAINTS.  PICK MINIMUM VIOLATION.
      IBEST=1
      IF (CV1.GE.CV2) IBEST=2
      GO TO 190
180   CONTINUE
C     NO CONSTRAINT VIOLATION.  PICK MINIMUM F.
      IBEST=1
      IF (F2.LE.F1) IBEST=2
190   CONTINUE
      II=1
C     ------------------------------------------------------------------
C     IF CV2 IS GREATER THAN CV1, SET MOVE LIMITS TO A2.
C     PROGRAM MOD-FEB, 1981, GV.
C     ------------------------------------------------------------------
      IF(CV2.LE.CV1) GO TO 195
      ALPLN=A2
      ALPNC=A2
      ALPCA=A2
195   CONTINUE
      IF (NCON.EQ.0) GO TO 230
C     ------------------------------------------------------------------
C     *****                 2 - POINT INTERPOLATION                *****
C     ------------------------------------------------------------------
      III=0
200   III=III+1
      C1=G1(III)
      C2=G(III)
      IF (ISC(III).EQ.0) GO TO 210
C     ------------------------------------------------------------------
C                        LINEAR CONSTRAINT
C     ------------------------------------------------------------------
      IF (C1.GE.1.0E-5.AND.C1.LE.CTBM) GO TO 220
      CALL CNMN07 (II,ALP,ZRO,ZRO,C1,A2,C2,ZRO,ZRO)
      IF (ALP.LE.0.) GO TO 220
      IF (C1.GT.CTBM.AND.ALP.GT.ALPFES) ALPFES=ALP
      IF (C1.LT.CTL.AND.ALP.LT.ALPLN) ALPLN=ALP
      GO TO 220
210   CONTINUE
C     ------------------------------------------------------------------
C                     NON-LINEAR CONSTRAINT
C     ------------------------------------------------------------------
      IF (C1.GE.1.0E-5.AND.C1.LE.CTAM) GO TO 220
      CALL CNMN07 (II,ALP,ZRO,ZRO,C1,A2,C2,ZRO,ZRO)
      IF (ALP.LE.0.) GO TO 220
      IF (C1.GT.CTAM.AND.ALP.GT.ALPFES) ALPFES=ALP
      IF (C1.LT.CT.AND.ALP.LT.ALPNC) ALPNC=ALP
220   CONTINUE
      IF (III.LT.NCON) GO TO 200
230   CONTINUE
      IF (LINOBJ.GT.0.OR.SLOPE.GE.0.) GO TO 240
C     CALCULATE ALPHA TO MINIMIZE FUNCTION.
      CALL CNMN04 (II,ALPMIN,ZRO,ZRO,F1,SLOPE,A2,F2,ZRO,ZRO,ZRO,ZRO)
240   CONTINUE
C     ------------------------------------------------------------------
C                         PROPOSED MOVE
C     ------------------------------------------------------------------
C     MOVE AT LEAST FAR ENOUGH TO OVERCOME CONSTRAINT VIOLATIONS.
      A3=ALPFES
C     MOVE TO MINIMIZE FUNCTION.
      IF (ALPMIN.GT.A3) A3=ALPMIN
C     IF A3.LE.0, SET A3 = ALPSID.
      IF (A3.LE.0.) A3=ALPSID
C     LIMIT MOVE TO NEW CONSTRAINT ENCOUNTER.
      IF (A3.GT.ALPNC) A3=ALPNC
      IF (A3.GT.ALPLN) A3=ALPLN
C     MAKE A3 NON-ZERO.
      IF (A3.LE.1.0E-20) A3=1.0E-20
C     IF A3=A2=ALPSID AND F2 IS BEST, GO INVOKE SIDE CONSTRAINT
C     MODIFICATION.
      ALPB=1.-A2/A3
      ALPA=1.-ALPSID/A3
      JBEST=0
      IF (ABS(ALPB).LT.1.0E-10.AND.ABS(ALPA).LT.1.0E-10) JBEST=1
      IF (JBEST.EQ.1.AND.IBEST.EQ.2) GO TO 20
C     SIDE CONSTRAINT CHECK NOT SATISFIED.
      IF (NCON.EQ.0) GO TO 260
C     STORE CONSTRAINT VALUES IN G2.
      DO 250 I=1,NCON
      G2(I)=G(I)
250   CONTINUE
260   CONTINUE
C     IF A3=A2, SET A3=.9*A2.
      IF (ABS(ALPB).LT.1.0E-10) A3=.9*A2
C     MOVE AT LEAST .01*A2.
      IF (A3.LT.(.01*A2)) A3=.01*A2
C     LIMIT MOVE TO 5.*A2.
      IF (A3.GT.(5.*A2)) A3=5.*A2
C     LIMIT MOVE TO ALPSID.
      IF (A3.GT.ALPSID) A3=ALPSID
C     MOVE A DISTANCE A3*S.
      ALP=A3-A2
      ALPTOT=ALPTOT+ALP
      DO 270 I=1,NDV
      X(I)=X(I)+ALP*S(I)
270   CONTINUE
      IF (IPRINT.LT.5) GO TO 300
      WRITE (IOUT,780)
      WRITE (IOUT,740) A3
      IF (NSCAL.EQ.0) GO TO 290
      DO 280 I=1,NDV
280   G(I)=SCAL(I)*X(I)
      WRITE (IOUT,750) (G(I),I=1,NDV)
      GO TO 300
290   WRITE (IOUT,750) (X(I),I=1,NDV)
300   CONTINUE
C     ------------------------------------------------------------------
C              UPDATE FUNCTION AND CONSTRAINT VALUES
C     ------------------------------------------------------------------
      NCAL(1)=NCAL(1)+1
      JGOTO=2
      RETURN
310   CONTINUE
      F3=OBJ
      IF (IPRINT.GE.5) WRITE (IOUT,760) F3
      IF (IPRINT.LT.5.OR.NCON.EQ.0) GO TO 320
      WRITE (IOUT,770)
      WRITE (IOUT,750) (G(I),I=1,NCON)
320   CONTINUE
C     ------------------------------------------------------------------
C       CALCULATE MAXIMUM CONSTRAINT VIOLATION AND PICK BEST DESIGN
C     ------------------------------------------------------------------
      CV3=0.
      IGOOD3=0
      NVC1=0
      IF (NCON.EQ.0) GO TO 340
      DO 330 I=1,NCON
      CC=CTAM
      IF (ISC(I).GT.0) CC=CTBM
      C1=G(I)-CC
      IF (C1.GT.CV3) CV3=C1
      IF (C1.GT.0.) NVC1=NVC1+1
330   CONTINUE
      IF (CV3.GT.0.) IGOOD3=1
340   CONTINUE
C     DETERMINE BEST DESIGN.
      IF (IBEST.EQ.2) GO TO 360
C     CHOOSE BETWEEN F1 AND F3.
      IF (IGOOD1.EQ.0.AND.IGOOD3.EQ.0) GO TO 350
      IF (CV1.GE.CV3) IBEST=3
      GO TO 380
350   IF (F3.LE.F1) IBEST=3
      GO TO 380
360   CONTINUE
C     CHOOSE BETWEEN F2 AND F3.
      IF (IGOOD2.EQ.0.AND.IGOOD3.EQ.0) GO TO 370
      IF (CV2.GE.CV3) IBEST=3
      GO TO 380
370   IF (F3.LE.F2) IBEST=3
380   CONTINUE
      ALP=A3
      OBJ=F3
C     IF F3 VIOLATES FEWER CONSTRAINTS THAN F1 RETURN.
      IF (NVC1.LT.NVC) GO TO 710
C     IF OBJECTIVE AND ALL CONSTRAINTS ARE LINEAR, RETURN.
      IF (LINOBJ.NE.0.AND.NLNC.EQ.NCON) GO TO 710
C     IF A3 = ALPLN AND F3 IS BOTH GOOD AND BEST RETURN.
      ALPB=1.-ALPLN/A3
      IF ((ABS(ALPB).LT.1.0E-20.AND.IBEST.EQ.3).AND.(IGOOD3.EQ.0)) GO TO
     1 710
C     IF A3 = ALPSID AND F3 IS BEST, GO INVOKE SIDE CONSTRAINT
C     MODIFICATION.
      ALPA=1.-ALPSID/A3
      IF (ABS(ALPA).LT.1.0E-20.AND.IBEST.EQ.3) GO TO 20
C     ------------------------------------------------------------------
C     **********            3 - POINT INTERPOLATION            *********
C     ------------------------------------------------------------------
      ALPNC=ALPSID
      ALPCA=ALPSID
      ALPFES=-1.
      ALPMIN=-1.
C     ------------------------------------------------------------------
C     IF A3 IS GREATER THAN A2 AND CV3 IS GREATER THAN CV2, SET
C     MOVE LIMITS TO A3.  PROGRAM MOD-FEB, 1981, GV.
C     ------------------------------------------------------------------
      IF(A3.LE.A2.OR.CV3.LE.CV2) GO TO 285
      ALPLN=A3
      ALPNC=A3
      ALPCA=A3
285   CONTINUE
      IF (NCON.EQ.0) GO TO 440
      III=0
390   III=III+1
      C1=G1(III)
      C2=G2(III)
      C3=G(III)
      IF (ISC(III).EQ.0) GO TO 400
C     ------------------------------------------------------------------
C     LINEAR CONSTRAINT.  FIND ALPFES ONLY.  ALPLN SAME AS BEFORE.
C     ------------------------------------------------------------------
      IF (C1.LE.CTBM) GO TO 430
      II=1
      CALL CNMN07 (II,ALP,ZRO,ZRO,C1,A3,C3,ZRO,ZRO)
      IF (ALP.GT.ALPFES) ALPFES=ALP
      GO TO 430
400   CONTINUE
C     ------------------------------------------------------------------
C                     NON-LINEAR CONSTRAINT
C     ------------------------------------------------------------------
      II=2
      CALL CNMN07 (II,ALP,ZRO,ZRO,C1,A2,C2,A3,C3)
      IF (ALP.LE.ZRO) GO TO 430
      IF (C1.GE.CT.AND.C1.LE.0.) GO TO 410
      IF (C1.GT.CTAM.OR.C1.LT.0.) GO TO 420
C     ALP IS MINIMUM MOVE.  UPDATE FOR NEXT CONSTRAINT ENCOUNTER.
410   ALPA=ALP
      CALL CNMN07 (II,ALP,ALPA,ZRO,C1,A2,C2,A3,C3)
      IF (ALP.LT.ALPCA.AND.ALP.GE.ALPA) ALPCA=ALP
      GO TO 430
420   CONTINUE
      IF (ALP.GT.ALPFES.AND.C1.GT.CTAM) ALPFES=ALP
      IF (ALP.LT.ALPNC.AND.C1.LT.0.) ALPNC=ALP
430   CONTINUE
      IF (III.LT.NCON) GO TO 390
440   CONTINUE
      IF (LINOBJ.GT.0.OR.SLOPE.GT.0.) GO TO 450
C     ------------------------------------------------------------------
C              CALCULATE ALPHA TO MINIMIZE FUNCTION
C     ------------------------------------------------------------------
      II=3
      IF (A2.GT.A3.AND.(IGOOD2.EQ.0.AND.IBEST.EQ.2)) II=2
      CALL CNMN04 (II,ALPMIN,ZRO,ZRO,F1,SLOPE,A2,F2,A3,F3,ZRO,ZRO)
450   CONTINUE
C     ------------------------------------------------------------------
C                       PROPOSED MOVE
C     ------------------------------------------------------------------
C     MOVE AT LEAST ENOUGH TO OVERCOME CONSTRAINT VIOLATIONS.
      A4=ALPFES
C     MOVE TO MINIMIZE FUNCTION.
      IF (ALPMIN.GT.A4) A4=ALPMIN
C     IF A4.LE.0, SET A4 = ALPSID.
      IF (A4.LE.0.) A4=ALPSID
C     LIMIT MOVE TO NEW CONSTRAINT ENCOUNTER.
      IF (A4.GT.ALPLN) A4=ALPLN
      IF (A4.GT.ALPNC) A4=ALPNC
C     LIMIT MOVE TO RE-ENCOUNTER CURRENTLY ACTIVE CONSTRAINT.
      IF (A4.GT.ALPCA) A4=ALPCA
C     LIMIT A4 TO 5.*A3.
      IF (A4.GT.(5.*A3)) A4=5.*A3
C     UPDATE DESIGN.
      IF (IBEST.NE.3.OR.NCON.EQ.0) GO TO 470
C     STORE CONSTRAINT VALUES IN G2.  F3 IS BEST.  F2 IS NOT.
      DO 460 I=1,NCON
      G2(I)=G(I)
460   CONTINUE
470   CONTINUE
C     IF A4=A3 AND IGOOD1=0 AND IGOOD3=1, SET A4=.9*A3.
      ALP=A4-A3
      IF ((IGOOD1.EQ.0.AND.IGOOD3.EQ.1).AND.ABS(ALP).LT.1.0E-20) A4=.9*A
     13
C     ------------------------------------------------------------------
C                   MOVE A DISTANCE A4*S
C     ------------------------------------------------------------------
      ALP=A4-A3
      ALPTOT=ALPTOT+ALP
      DO 480 I=1,NDV
      X(I)=X(I)+ALP*S(I)
480   CONTINUE
      IF (IPRINT.LT.5) GO TO 510
      WRITE (IOUT,720)
      WRITE (IOUT,740) A4
      IF (NSCAL.EQ.0) GO TO 500
      DO 490 I=1,NDV
490   G(I)=SCAL(I)*X(I)
      WRITE (IOUT,750) (G(I),I=1,NDV)
      GO TO 510
500   WRITE (IOUT,750) (X(I),I=1,NDV)
510   CONTINUE
C     ------------------------------------------------------------------
C              UPDATE FUNCTION AND CONSTRAINT VALUES
C     ------------------------------------------------------------------
      NCAL(1)=NCAL(1)+1
      JGOTO=3
      RETURN
520   CONTINUE
      F4=OBJ
      IF (IPRINT.GE.5) WRITE (IOUT,760) F4
      IF (IPRINT.LT.5.OR.NCON.EQ.0) GO TO 530
      WRITE (IOUT,770)
      WRITE (IOUT,750) (G(I),I=1,NCON)
530   CONTINUE
C     DETERMINE ACCAPTABILITY OF F4.
      IGOOD4=0
      CV4=0.
      IF (NCON.EQ.0) GO TO 550
      DO 540 I=1,NCON
      CC=CTAM
      IF (ISC(I).GT.0) CC=CTBM
      C1=G(I)-CC
      IF (C1.GT.CV4) CV4=C1
540   CONTINUE
      IF (CV4.GT.0.) IGOOD4=1
550   CONTINUE
      ALP=A4
      OBJ=F4
C     ------------------------------------------------------------------
C                     DETERMINE BEST DESIGN
C     ------------------------------------------------------------------
      GO TO (560,610,660),IBEST
560   CONTINUE
C     CHOOSE BETWEEN F1 AND F4.
      IF (IGOOD1.EQ.0.AND.IGOOD4.EQ.0) GO TO 570
      IF (CV1.GT.CV4) GO TO 710
      GO TO 580
570   CONTINUE
      IF (F4.LE.F1) GO TO 710
580   CONTINUE
C     F1 IS BEST.
      ALPTOT=ALPTOT-A4
      OBJ=F1
      DO 590 I=1,NDV
      X(I)=X(I)-A4*S(I)
590   CONTINUE
      IF (NCON.EQ.0) GO TO 710
      DO 600 I=1,NCON
      G(I)=G1(I)
600   CONTINUE
      GO TO 710
610   CONTINUE
C     CHOOSE BETWEEN F2 AND F4.
      IF (IGOOD2.EQ.0.AND.IGOOD4.EQ.0) GO TO 620
      IF (CV2.GT.CV4) GO TO 710
      GO TO 630
620   CONTINUE
      IF (F4.LE.F2) GO TO 710
630   CONTINUE
C     F2 IS BEST.
      OBJ=F2
      A2=A4-A2
      ALPTOT=ALPTOT-A2
      DO 640 I=1,NDV
      X(I)=X(I)-A2*S(I)
640   CONTINUE
      IF (NCON.EQ.0) GO TO 710
      DO 650 I=1,NCON
      G(I)=G2(I)
650   CONTINUE
      GO TO 710
660   CONTINUE
C     CHOOSE BETWEEN F3 AND F4.
      IF (IGOOD3.EQ.0.AND.IGOOD4.EQ.0) GO TO 670
      IF (CV3.GT.CV4) GO TO 710
      GO TO 680
670   CONTINUE
      IF (F4.LE.F3) GO TO 710
680   CONTINUE
C     F3 IS BEST.
      OBJ=F3
      A3=A4-A3
      ALPTOT=ALPTOT-A3
      DO 690 I=1,NDV
      X(I)=X(I)-A3*S(I)
690   CONTINUE
      IF (NCON.EQ.0) GO TO 710
      DO 700 I=1,NCON
      G(I)=G2(I)
700   CONTINUE
710   CONTINUE
      ALP=ALPTOT
      IF (IPRINT.GE.5) WRITE (IOUT,790)
      JGOTO=0
      RETURN
C     ------------------------------------------------------------------
C                                  FORMATS
C     ------------------------------------------------------------------
C
C
720   FORMAT (/5X,25HTHREE-POINT INTERPOLATION)
730   FORMAT (/////58H* * * CONSTRAINED ONE-DIMENSIONAL SEARCH INFORMATI
     1ON * * *)
740   FORMAT (//5X,15HPROPOSED DESIGN/5X,7HALPHA =,E12.5/5X,8HX-VECTOR)
750   FORMAT (1X,8E12.4)
760   FORMAT (/5X,5HOBJ =,E13.5)
770   FORMAT (/5X,17HCONSTRAINT VALUES)
780   FORMAT (/5X,23HTWO-POINT INTERPOLATION)
790   FORMAT (/5X,35H* * * END OF ONE-DIMENSIONAL SEARCH)
      END

max_stars_repo_name: robfalck/pyoptsparse

max_issues_repo_name: robfalck/pyoptsparse

max_forks_repo_name: robfalck/pyoptsparse