text:
      SUBROUTINE TG01MD( JOB, N, M, P, A, LDA, E, LDE, B, LDB, C, LDC,
     $                   ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, NF, ND,
     $                   NIBLCK, IBLCK, TOL, IWORK, DWORK, LDWORK,
     $                   INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute orthogonal transformation matrices Q and Z which
C     reduce the regular pole pencil A-lambda*E of the descriptor system
C     (A-lambda*E,B,C) to the form (if JOB = 'F')
C
C                ( Af  *  )             ( Ef  *  )
C       Q'*A*Z = (        ) ,  Q'*E*Z = (        ) ,                 (1)
C                ( 0   Ai )             ( 0   Ei )
C
C     or to the form (if JOB = 'I')
C
C                ( Ai  *  )             ( Ei  *  )
C       Q'*A*Z = (        ) ,  Q'*E*Z = (        ) ,                 (2)
C                ( 0   Af )             ( 0   Ef )
C
C     where the pair (Af,Ef) is in a generalized real Schur form, with
C     Ef nonsingular and upper triangular and Af in real Schur form.
C     The subpencil Af-lambda*Ef contains the finite eigenvalues.
C     The pair (Ai,Ei) is in a generalized real Schur form with
C     both Ai and Ei upper triangular. The subpencil Ai-lambda*Ei,
C     with Ai nonsingular and Ei nilpotent contains the infinite
C     eigenvalues and is in a block staircase form (see METHOD).
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOB     CHARACTER*1
C             = 'F':  perform the finite-infinite separation;
C             = 'I':  perform the infinite-finite separation.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The number of rows of the matrix B, the number of columns
C             of the matrix C and the order of the square matrices A
C             and E.  N >= 0.
C
C     M       (input) INTEGER
C             The number of columns of the matrix B.  M >= 0.
C
C     P       (input) INTEGER
C             The number of rows of the matrix C.  P >= 0.
C
C     A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
C             On entry, the leading N-by-N part of this array must
C             contain the N-by-N state matrix A.
C             On exit, the leading N-by-N part of this array contains
C             the transformed state matrix Q'*A*Z in the form
C
C             ( Af  *  )                    ( Ai  *  )
C             (        ) for JOB = 'F', or  (        )  for JOB = 'I',
C             ( 0   Ai )                    ( 0   Af )
C
C             where Af is an NF-by-NF matrix in real Schur form, and Ai
C             is an (N-NF)-by-(N-NF) nonsingular and upper triangular
C             matrix. Ai has a block structure as in (3) or (4), where
C             A0,0 is ND-by-ND and Ai,i , for i = 1, ..., NIBLCK, is
C             IBLCK(i)-by-IBLCK(i). (See METHOD.)
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1,N).
C
C     E       (input/output) DOUBLE PRECISION array, dimension (LDE,N)
C             On entry, the leading N-by-N part of this array must
C             contain the N-by-N descriptor matrix E.
C             On exit, the leading N-by-N part of this array contains
C             the transformed descriptor matrix Q'*E*Z in the form
C
C             ( Ef  *  )                    ( Ei  *  )
C             (        ) for JOB = 'F', or  (        )  for JOB = 'I',
C             ( 0   Ei )                    ( 0   Ef )
C
C             where Ef is an NF-by-NF nonsingular and upper triangular
C             matrix, and Ei is an (N-NF)-by-(N-NF) nilpotent matrix in
C             an upper triangular block form as in (3) or (4).
C
C     LDE     INTEGER
C             The leading dimension of the array E.  LDE >= MAX(1,N).
C
C     B       (input/output) DOUBLE PRECISION array, dimension (LDB,M)
C             On entry, the leading N-by-M part of this array must
C             contain the N-by-M input matrix B.
C             On exit, the leading N-by-M part of this array contains
C             the transformed input matrix Q'*B.
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= MAX(1,N).
C
C     C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
C             On entry, the leading P-by-N part of this array must
C             contain the state/output matrix C.
C             On exit, the leading P-by-N part of this array contains
C             the transformed matrix C*Z.
C
C     LDC     INTEGER
C             The leading dimension of the array C.  LDC >= MAX(1,P).
C
C     ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
C             ALPHAR(1:NF) will be set to the real parts of the diagonal
C             elements of Af that would result from reducing A and E to
C             the Schur form, and then further reducing both of them to
C             triangular form using unitary transformations, subject to
C             having the diagonal of E positive real.  Thus, if Af(j,j)
C             is in a 1-by-1 block (i.e., Af(j+1,j) = Af(j,j+1) = 0),
C             then ALPHAR(j) = Af(j,j). Note that the (real or complex)
C             values (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,NF, are
C             the finite generalized eigenvalues of the matrix pencil
C             A - lambda*E.
C
C     ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
C             ALPHAI(1:NF) will be set to the imaginary parts of the
C             diagonal elements of Af that would result from reducing A
C             and E to Schur form, and then further reducing both of
C             them to triangular form using unitary transformations,
C             subject to having the diagonal of E positive real. Thus,
C             if Af(j,j) is in a 1-by-1 block (see above), then
C             ALPHAI(j) = 0. Note that the (real or complex) values
C             (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,NF, are the
C             finite generalized eigenvalues of the matrix pencil
C             A - lambda*E.
C
C     BETA    (output) DOUBLE PRECISION array, dimension (N)
C             BETA(1:NF) will be set to the (real) diagonal elements of
C             Ef that would result from reducing A and E to Schur form,
C             and then further reducing both of them to triangular form
C             using unitary transformations, subject to having the
C             diagonal of E positive real. Thus, if Af(j,j) is in a
C             1-by-1 block (see above), then BETA(j) = Ef(j,j).
C             Note that the (real or complex) values
C             (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,NF, are the
C             finite generalized eigenvalues of the matrix pencil
C             A - lambda*E.
C
C     Q       (output) DOUBLE PRECISION array, dimension (LDQ,N)
C             The leading N-by-N part of this array contains the
C             orthogonal matrix Q, which is the accumulated product of
C             the transformations applied to A, E, and B on the left.
C
C     LDQ     INTEGER
C             The leading dimension of the array Q.  LDQ >= MAX(1,N).
C
C     Z       (output) DOUBLE PRECISION array, dimension (LDZ,N)
C             The leading N-by-N part of this array contains the
C             orthogonal matrix Z, which is the accumulated product of
C             the transformations applied to A, E, and C on the right.
C
C     LDZ     INTEGER
C             The leading dimension of the array Z.  LDZ >= MAX(1,N).
C
C     NF      (output) INTEGER
C             The order of the reduced matrices Af and Ef; also, the
C             number of finite generalized eigenvalues of the pencil
C             A-lambda*E.
C
C     ND      (output) INTEGER
C             The number of non-dynamic infinite eigenvalues of the
C             matrix pair (A,E). Note: N-ND is the rank of the matrix E.
C
C     NIBLCK  (output) INTEGER
C             If ND > 0, the number of infinite blocks minus one.
C             If ND = 0, then NIBLCK = 0.
C
C     IBLCK   (output) INTEGER array, dimension (N)
C             IBLCK(i) contains the dimension of the i-th block in the
C             staircase form (3), where i = 1,2,...,NIBLCK.
C
C     Tolerances
C
C     TOL     DOUBLE PRECISION
C             A tolerance used in rank decisions to determine the
C             effective rank, which is defined as the order of the
C             largest leading (or trailing) triangular submatrix in the
C             QR factorization with column pivoting whose estimated
C             condition number is less than 1/TOL. If the user sets
C             TOL <= 0, then an implicitly computed, default tolerance
C             TOLDEF = N**2*EPS,  is used instead, where EPS is the
C             machine precision (see LAPACK Library routine DLAMCH).
C             TOL < 1.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (N)
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0, DWORK(1) returns the optimal value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The length of the array DWORK.  LDWORK >= 1, and if N > 0,
C             LDWORK >= 4*N.
C
C             If LDWORK = -1, then a workspace query is assumed; the
C             routine only calculates the optimal size of the DWORK
C             array, returns this value as the first entry of the DWORK
C             array, and no error message related to LDWORK is issued by
C             XERBLA.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             = 1:  the pencil A-lambda*E is not regular;
C             = 2:  the QZ iteration did not converge.
C
C     METHOD
C
C     For the separation of infinite structure, the reduction algorithm
C     of [1] is employed.
C     If JOB = 'F', the matrices Ai and Ei have the form
C
C           ( A0,0  A0,k ... A0,1 )         ( 0  E0,k ... E0,1 )
C      Ai = (  0    Ak,k ... Ak,1 ) ,  Ei = ( 0   0   ... Ek,1 ) ;   (3)
C           (  :     :    .    :  )         ( :   :    .    :  )
C           (  0     0   ... A1,1 )         ( 0   0   ...   0  )
C
C     if JOB = 'I' the matrices Ai and Ei have the form
C
C           ( A1,1 ... A1,k  A1,0 )         ( 0 ... E1,k  E1,0 )
C      Ai = (  :    .    :    :   ) ,  Ei = ( :  .    :    :   ) ,   (4)
C           (  :   ... Ak,k  Ak,0 )         ( : ...   0   Ek,0 )
C           (  0   ...   0   A0,0 )         ( 0 ...   0     0  )
C
C     where Ai,i, for i = 0, 1, ..., k, are nonsingular upper triangular
C     matrices. A0,0 corresponds to the non-dynamic infinite modes of
C     the system.
C
C     REFERENCES
C
C     [1] Misra, P., Van Dooren, P., and Varga, A.
C         Computation of structural invariants of generalized
C         state-space systems.
C         Automatica, 30, pp. 1921-1936, 1994.
C
C     NUMERICAL ASPECTS
C
C     The algorithm is numerically backward stable and requires
C     0( N**3 )  floating point operations.
C
C     FURTHER COMMENTS
C
C     The number of infinite poles is computed as
C
C                   NIBLCK
C        NINFP =     Sum  IBLCK(i) = N - ND - NF.
C                    i=1
C     The multiplicities of infinite poles can be computed as follows:
C     there are IBLCK(k)-IBLCK(k+1) infinite poles of multiplicity
C     k, for k = 1, ..., NIBLCK, where IBLCK(NIBLCK+1) = 0.
C     Note that each infinite pole of multiplicity k corresponds to
C     an infinite eigenvalue of multiplicity k+1.
C
C     CONTRIBUTOR
C
C     A. Varga, German Aerospace Center, DLR Oberpfaffenhofen.
C     July 1999. Based on the RASP routines SRISEP and RPDSGH.
C
C     REVISIONS
C
C     V. Sima, Dec. 2016.
C
C     KEYWORDS
C
C     Generalized eigenvalue problem, system poles, multivariable
C     system, orthogonal transformation, structural invariant.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D0, ZERO = 0.0D0 )
C     .. Scalar Arguments ..
      CHARACTER          JOB
      INTEGER            INFO, LDA, LDB, LDC, LDE, LDQ, LDWORK, LDZ, M,
     $                   N, ND, NF, NIBLCK, P
      DOUBLE PRECISION   TOL
C     .. Array Arguments ..
      INTEGER            IBLCK( * ), IWORK(*)
      DOUBLE PRECISION   A(LDA,*), ALPHAI(*), ALPHAR(*), B(LDB,*),
     $                   BETA(*),  C(LDC,*),  DWORK(*),  E(LDE,*),
     $                   Q(LDQ,*), Z(LDZ,*)
C     .. Local Scalars ..
      LOGICAL            LQUERY, TRINF
      INTEGER            I, IHI, ILO, MINWRK, NBC, NC, NR, WRKOPT
C     .. Local Arrays ..
      DOUBLE PRECISION   DUM(1)
C     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
C     .. External Subroutines ..
      EXTERNAL           DGEMM, DHGEQZ, DLACPY, TG01LD, XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC          INT, MAX, MIN
C     .. Executable Statements ..
C
C     Test the input parameters.
C
      INFO  =  0
      TRINF =  LSAME( JOB, 'I' )
      IF( .NOT.LSAME( JOB, 'F' ) .AND. .NOT.TRINF ) THEN
         INFO = -1
      ELSE IF( N.LT.0 ) THEN
         INFO = -2
      ELSE IF( M.LT.0 ) THEN
         INFO = -3
      ELSE IF( P.LT.0 ) THEN
         INFO = -4
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -6
      ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
         INFO = -8
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -10
      ELSE IF( LDC.LT.MAX( 1, P ) ) THEN
         INFO = -12
      ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
         INFO = -17
      ELSE IF( LDZ.LT.MAX( 1, N ) ) THEN
         INFO = -19
      ELSE IF( TOL.GE.ONE ) THEN
         INFO = -24
      ELSE
         LQUERY = LDWORK.EQ.-1
         IF( N.EQ.0 ) THEN
            MINWRK = 1
         ELSE
            MINWRK = 4*N
         END IF
         IF( LQUERY ) THEN
            ILO = 1
            IHI = N
C
            CALL TG01LD( JOB, 'Hessenberg', 'Identity', 'Identity', N,
     $                   0, 0, A, LDA, E, LDE, DUM, LDB, DUM, LDC, Q,
     $                   LDQ, Z, LDZ, NF, ND, NIBLCK, IBLCK, TOL, IWORK,
     $                   DWORK, -1, INFO )
            WRKOPT = MAX( MINWRK, INT( DWORK(1) ) )
            CALL DHGEQZ( 'Schur', 'Vector', 'Vector', N, ILO, IHI, A,
     $                   LDA, E, LDE, ALPHAR, ALPHAI, BETA, Q, LDQ, Z,
     $                   LDZ, DWORK, -1, INFO )
            WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
         END IF
         IF( LDWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
            INFO = -27
         END IF
      END IF
C
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'TG01MD', -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         DWORK(1) = WRKOPT
         RETURN
      END IF
C
C     Quick return if possible
C
      IF( N.EQ.0 ) THEN
         NF = 0
         ND = 0
         NIBLCK   = 0
         DWORK(1) = ONE
         RETURN
      END IF
C
C     Compute the finite-infinite separation with A in Hessenberg form
C     and E upper triangular.
C     Workspace: need  4*N.
C
      CALL TG01LD( JOB, 'Hessenberg', 'Identity', 'Identity', N, 0, 0,
     $             A, LDA, E, LDE, DUM, LDB, DUM, LDC, Q, LDQ, Z, LDZ,
     $             NF, ND, NIBLCK, IBLCK, TOL, IWORK, DWORK, LDWORK,
     $             INFO )
C
      IF( INFO.NE.0 )
     $   RETURN
      WRKOPT = MAX( MINWRK, INT( DWORK(1) ) )
C
C     Reduce the finite part to generalized real Schur form.
C     Workspace: need  N;
C                prefer larger.
C
      IF( TRINF ) THEN
         ILO = N - NF + 1
         IHI = N
      ELSE
         ILO = 1
         IHI = NF
      END IF
      CALL DHGEQZ( 'Schur', 'Vector', 'Vector', N, ILO, IHI, A, LDA,
     $             E, LDE, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, DWORK,
     $             LDWORK, INFO )
      IF( INFO.NE.0 ) THEN
         INFO = 2
         RETURN
      END IF
      WRKOPT = MAX( WRKOPT, INT( DWORK(1) ) )
C
C     B <-- Q' * B.
C     Workspace: need   N;
C                prefer N*M.
C
      NBC = MAX( 1, MIN( LDWORK/N, M ) )
      DO 10 I = 1, M, NBC
         NC = MIN( NBC, M-I+1 )
         CALL DGEMM( 'Transpose', 'No transpose', N, NC, N, ONE, Q, LDQ,
     $               B(1,I), LDB, ZERO, DWORK, N )
         CALL DLACPY( 'All', N, NC, DWORK, N, B(1,I), LDB )
   10 CONTINUE
C
C     C <-- C * Z
C     Workspace: need   N;
C                prefer P*N.
C
      NBC = MAX( 1, MIN( LDWORK/N, P ) )
      DO 20 I = 1, P, NBC
         NR = MIN( NBC, P-I+1 )
         CALL DGEMM( 'No Transpose', 'No transpose', NR, N, N, ONE,
     $              C(I,1), LDC, Z, LDZ, ZERO, DWORK, NR )
         CALL DLACPY( 'All', NR, N, DWORK, NR, C(I,1), LDC )
   20 CONTINUE
C
      DWORK(1) = WRKOPT
C
      RETURN
C *** Last line of TG01MD ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


text:
      MODULE MOD_GEOMESH

      CONTAINS

      SUBROUTINE GEOMESH(RADIUS, ENTOT, T, FAL, P, Z, RSTAR, AMU, ATMEAN, ND, NP)

!     THIS SUBROUTINE GENERATES THE GEOMETRICAL POINT MESH IN RADIUS, P AND Z
!     P and Z mesh is needed for the ray-by-ray solution of the radiative transfer equation in spherical symmetry

      INTEGER, INTENT(OUT) ::  ND, NP
      REAL*8,  INTENT(IN)  ::  RSTAR

      real*8, intent(in) :: AMU, ATMEAN

      logical, intent(in) :: FAL

      real*8, allocatable, dimension(:), intent(out) :: Z, P
      real*8, allocatable, dimension(:), intent(out) :: T, entot, radius

      CALL RGRIDM(RADIUS, ENTOT, T, FAL, RSTAR, AMU, ATMEAN, ND)

      CALL PGRID(NP, ND, RADIUS, P)

      CALL ZGRID(RADIUS, P, Z, ND, NP)

      return

      end subroutine

      SUBROUTINE RGRIDM(radius, entot, T, FAL, rstar, AMU, ATMEAN, ND)

      USE COMMON_BLOCK
      use FILE_OPERATIONS

      integer, intent(out) :: ND

      real*8,  intent(in)  :: rstar

      logical, intent(in) :: FAL

      real*8, intent(in) :: AMU, ATMEAN

      real*8, allocatable, dimension(:), intent(out) :: T, entot, radius

!     height: height in km
!     T:      Temperature in K
!     entot:  HEAVY PARTICLE DENSITY
!     radius: HEIGHT IN UNITS OF SOLAR RADII

      if (FAL)       call read_fal_mod(rstar, height, radius, T, entot, ND)

      if (.not. FAL) call read_kur_mod(amu, atmean, rstar, height, radius, T, entot, ND)

      dpn = ND

      end subroutine

      SUBROUTINE PGRID(NP, ND, R, P)
!     GRID OF IMPACT-PARAMETER POINTS

      integer,                intent(in) ::  ND
      real*8,  dimension(ND), intent(in) ::  R

      integer, intent(out) :: NP

      real*8, allocatable, dimension(:), intent(out) :: P

      real*8, dimension(13) :: cp

      data cp /0.0d0, 0.1d0, 0.2d0, 0.3d0, 0.4d0, 0.5d0, 0.6d0, 0.7d0, 0.8d0,
     $         0.91652d0, 0.97980d0, 0.99499d0, 0.99875d0/

!     NC = NUMBER OF CORE-INTERSECTING RAYS

      NC = 13; NP = ND + NC

      allocate(P(NP))

      P(1 : NC) = cp(1 : NC)

!     CORE RAYS EQUELLY SPACED

      do L = 1, ND; J = NP + 1 - L; P(J) = R(L); enddo

      return

      end subroutine

      SUBROUTINE ZGRID(RADIUS, P, Z, ND, NP)

      implicit none
!     THIS SUBROUTINE GENERATES THE GEOMETRICAL POINT MESH IN Z

      integer, intent(in)  :: ND, NP
      real*8,  intent(in)  :: RADIUS(ND), P(NP)

      real*8,  allocatable, dimension(:), intent(out) :: Z

      integer :: i, j, L, JMAX

      allocate(z(nd * np))

      z(1 : nd * np) = 0.0d0

      do L = 1, ND

         JMAX = NP + 1 - L

         do j = 1, JMAX

            i = (j - 1) * ND + L

            Z(i) = SQRT((radius(l) - P(j)) * (radius(l) + P(j)))

        enddo

      enddo

      RETURN

      END subroutine

      subroutine read_fal_mod(rstar, height, radius, T, entot, ND)

      use file_operations

      real*8,  intent(in)  :: rstar

      integer, intent(out) :: ND

      real*8, allocatable, dimension(:), intent(out) :: T, entot, radius, height

      ND = num_of_lines(fal_mod_file)

!          if (allocated(T))      deallocate(T)
!          if (allocated(entot))  deallocate(entot)
!          if (allocated(radius)) deallocate(radius)
          if (allocated(height)) deallocate(height)

      allocate(T(ND))
      allocate(entot(ND))
      allocate(radius(ND))
      allocate(height(ND))

      height = read_atm_mod(fal_mod_file, '1')
      T =      read_atm_mod(fal_mod_file, '2')
      entot =  read_atm_mod(fal_mod_file, '4')

      radius = 1.0d0 + height * 1.0d5 / rstar ! height in km, rstar and radius in cm

      end subroutine

      subroutine read_kur_mod(amu, atmean, rstar, height, radius, T, entot, ND)

      use file_operations

      real*8,  intent(in)  :: rstar

      real*8, intent(in) :: amu, atmean

      integer, intent(out) :: ND

      real*8, allocatable, dimension(:), intent(out) :: T, entot, radius, height

      real*8, allocatable, dimension(:) :: entotn, delr

      real*8, allocatable, dimension(:) :: elec_conc, rho, vturb, pressure

!     ak - Boltzmann constant
      real*8, parameter :: ak =  1.38062259d-16
      real*8, parameter :: MUN = 1.66054d-24

      ND = num_of_lines(kur_mod_file) - 1

!          if(allocated(T))         deallocate(T)
!          if(allocated(entot))     deallocate(entot)
!          if(allocated(radius))    deallocate(radius)
          if(allocated(height))    deallocate(height)
!          if(allocated(rho))       deallocate(rho)
!          if(allocated(pressure))  deallocate(pressure)
!          if(allocated(elec_conc)) deallocate(elec_conc)
!          if(allocated(vturb))     deallocate(vturb)
!          if(allocated(entotn))    deallocate(entotn)
!          if(allocated(delr))      deallocate(delr)

      allocate(T(ND))
      allocate(entot(ND))
      allocate(radius(ND))
      allocate(height(ND))
      allocate(rho(ND))
      allocate(pressure(ND))
      allocate(elec_conc(ND))
      allocate(vturb(ND))
      allocate(entotn(ND))
      allocate(delr(ND))

      rho =       read_atm_mod(kur_mod_file, '1')
      T =         read_atm_mod(kur_mod_file, '2')
      pressure =  read_atm_mod(kur_mod_file, '3')
      elec_conc = read_atm_mod(kur_mod_file, '4')
      vturb =     read_atm_mod(kur_mod_file, '7')

      pressure = rho * 10.0**4.44 ! this is from the old version and I do not understand what's it doing there

!         TAKING INTO ACCOUNT TURBULEN PRESSURE

      ENTOTN = pressure / (AK * T + 0.5 * ATMEAN * MUN * vturb**2.)

      ENTOT  = ENTOTN - elec_conc

      DO L = 1, ND - 1

         DELR(L) = (2 / (AMU * ATMEAN * RSTAR)) * (RHO(L + 1) - RHO(L)) / (ENTOT(L + 1) + ENTOT(L))

      ENDDO

      RADIUS(ND) = 1.0

      DO K = 1, ND - 1; RADIUS(ND - K) = RADIUS(ND - K + 1) + DELR(ND - K); ENDDO

      height = (radius - radius(ND)) * rstar / 1D+5 ! height in km, rstar and radius in cm

      deallocate(rho)
      deallocate(vturb)
      deallocate(elec_conc)
      deallocate(pressure)
      deallocate(delr)
      deallocate(entotn)

      end subroutine

      end module

ext: for

max_stars_repo_name: rtagirov/nessy

max_issues_repo_name: rtagirov/nessy

max_forks_repo_name: rtagirov/nessy-src


text:
!
! FNote: Command-line notebook application
! Module: Miscellaneous "helper" procedures for character string
! operations.
!

module str_utils
    implicit none

contains

    !
    !  name: replace_chars
    !  desc: Replaces characters in a character string. The string is
    !        modified in-place.
    !  @param str: allocatable character string. Will be modified!
    !  @param to_replace: character that will be replaced.
    !  @param replace_with: replacement character.
    !
    subroutine replace_chars(str, to_replace, replace_with)
        character(len=:), allocatable, intent(inout) :: str
        character, intent(in) :: to_replace
        character, intent(in) :: replace_with
        integer :: c_pos

        ! TODO: allocation check
        do
            c_pos = index(str, to_replace)

            if (c_pos == 0) then
                exit
            end if

            str(c_pos:c_pos) = replace_with
        end do
    end subroutine

end module

max_stars_repo_name: art-rasa/fnote

max_issues_repo_name: art-rasa/fnote

max_forks_repo_name: art-rasa/fnote


text:
prime_flag/reset.cpp
prime_flag/display.cpp
prime_flag/prime_numgen.cpp
prime_flag/main.cpp

ext: f


text:
Nations Giant Hamburgers sells hamburgers and pies. They have an http://www.nationsrestaurants.com/menu/index.html online menu that also has breakfast.

I didnt really think this was a big deal, Vacaville residents dont even go to Nations all that often. But Ive known of some Users/TarZxf stranger Davis residents to travel all the way out Vacaville here to go to Nations. So there you go Users/MichelleAccurso

20051218 21:00:16 nbsp There used to be a Nations in the University Mall. Users/JasonAller

20060312 10:04:10 nbsp I <3 Nations and the pies!  Im one of the people who will travel to Vacaville to get some D Users/JoAnnaRich

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C******************************************************************************
C FILE: omp_reduction.f
C DESCRIPTION:
C   OpenMP Example - Combined Parallel Loop Reduction - Fortran Version
C   This example demonstrates a sum reduction within a combined parallel loop
C   construct.  Notice that default data element scoping is assumed - there
C   are no clauses specifying shared or private variables.  OpenMP will
C   automatically make loop index variables private within team threads, and
C   global variables shared.
C AUTHOR: Blaise Barney  5/99
C LAST REVISED:
C******************************************************************************

      PROGRAM REDUCTION

      INTEGER I, N
      REAL A(100), B(100), SUM

!     Some initializations
      N = 100
      DO I = 1, N
        A(I) = I *1.0
        B(I) = A(I)
      ENDDO
      SUM = 0.0

!$OMP PARALLEL DO REDUCTION(+:SUM)
      DO I = 1, N
        SUM = SUM + (A(I) * B(I))
      ENDDO

      PRINT *, '   Sum = ', SUM
      END


text:
      subroutine sproot(t,n,c,zero,mest,m,ier)
c  subroutine sproot finds the zeros of a cubic spline s(x),which is
c  given in its normalized b-spline representation.
c
c  calling sequence:
c     call sproot(t,n,c,zero,mest,m,ier)
c
c  input parameters:
c    t    : real array,length n, containing the knots of s(x).
c    n    : integer, containing the number of knots.  n>=8
c    c    : real array,length n, containing the b-spline coefficients.
c    mest : integer, specifying the dimension of array zero.
c
c  output parameters:
c    zero : real array,lenth mest, containing the zeros of s(x).
c    m    : integer,giving the number of zeros.
c    ier  : error flag:
c      ier = 0: normal return.
c      ier = 1: the number of zeros exceeds mest.
c      ier =10: invalid input data (see restrictions).
c
c  other subroutines required: fpcuro
c
c  restrictions:
c    1) n>= 8.
c    2) t(4) < t(5) < ... < t(n-4) < t(n-3).
c       t(1) <= t(2) <= t(3) <= t(4)
c       t(n-3) <= t(n-2) <= t(n-1) <= t(n)
c
c  author :
c    p.dierckx
c    dept. computer science, k.u.leuven
c    celestijnenlaan 200a, b-3001 heverlee, belgium.
c    e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
c  latest update : march 1987
c
c ..
c ..scalar arguments..
      integer n,mest,m,ier
c  ..array arguments..
      real t(n),c(n),zero(mest)
c  ..local scalars..
      integer i,j,j1,l,n4
      real ah,a0,a1,a2,a3,bh,b0,b1,c1,c2,c3,c4,c5,d4,d5,h1,h2,
     * three,two,t1,t2,t3,t4,t5,zz
      logical z0,z1,z2,z3,z4,nz0,nz1,nz2,nz3,nz4
c  ..local array..
      real y(3)
c  ..
c  set some constants
      two = 0.2e+01
      three = 0.3e+01
c  before starting computations a data check is made. if the input data
c  are invalid, control is immediately repassed to the calling program.
      n4 = n-4
      ier = 10
      if(n.lt.8) go to 800
      j = n
      do 10 i=1,3
        if(t(i).gt.t(i+1)) go to 800
        if(t(j).lt.t(j-1)) go to 800
        j = j-1
  10  continue
      do 20 i=4,n4
        if(t(i).ge.t(i+1)) go to 800
  20  continue
c  the problem considered reduces to finding the zeros of the cubic
c  polynomials pl(x) which define the cubic spline in each knot
c  interval t(l)<=x<=t(l+1). a zero of pl(x) is also a zero of s(x) on
c  the condition that it belongs to the knot interval.
c  the cubic polynomial pl(x) is determined by computing s(t(l)),
c  s'(t(l)),s(t(l+1)) and s'(t(l+1)). in fact we only have to compute
c  s(t(l+1)) and s'(t(l+1)); because of the continuity conditions of
c  splines and their derivatives, the value of s(t(l)) and s'(t(l))
c  is already known from the foregoing knot interval.
      ier = 0
c  evaluate some constants for the first knot interval
      h1 = t(4)-t(3)
      h2 = t(5)-t(4)
      t1 = t(4)-t(2)
      t2 = t(5)-t(3)
      t3 = t(6)-t(4)
      t4 = t(5)-t(2)
      t5 = t(6)-t(3)
c  calculate a0 = s(t(4)) and ah = s'(t(4)).
      c1 = c(1)
      c2 = c(2)
      c3 = c(3)
      c4 = (c2-c1)/t4
      c5 = (c3-c2)/t5
      d4 = (h2*c1+t1*c2)/t4
      d5 = (t3*c2+h1*c3)/t5
      a0 = (h2*d4+h1*d5)/t2
      ah = three*(h2*c4+h1*c5)/t2
      z1 = .true.
      if(ah.lt.0.) z1 = .false.
      nz1 = .not.z1
      m = 0
c  main loop for the different knot intervals.
      do 300 l=4,n4
c  evaluate some constants for the knot interval t(l) <= x <= t(l+1).
        h1 = h2
        h2 = t(l+2)-t(l+1)
        t1 = t2
        t2 = t3
        t3 = t(l+3)-t(l+1)
        t4 = t5
        t5 = t(l+3)-t(l)
c  find a0 = s(t(l)), ah = s'(t(l)), b0 = s(t(l+1)) and bh = s'(t(l+1)).
        c1 = c2
        c2 = c3
        c3 = c(l)
        c4 = c5
        c5 = (c3-c2)/t5
        d4 = (h2*c1+t1*c2)/t4
        d5 = (h1*c3+t3*c2)/t5
        b0 = (h2*d4+h1*d5)/t2
        bh = three*(h2*c4+h1*c5)/t2
c  calculate the coefficients a0,a1,a2 and a3 of the cubic polynomial
c  pl(x) = ql(y) = a0+a1*y+a2*y**2+a3*y**3 ; y = (x-t(l))/(t(l+1)-t(l)).
        a1 = ah*h1
        b1 = bh*h1
        a2 = three*(b0-a0)-b1-two*a1
        a3 = two*(a0-b0)+b1+a1
c  test whether or not pl(x) could have a zero in the range
c  t(l) <= x <= t(l+1).
        z3 = .true.
        if(b1.lt.0.) z3 = .false.
        nz3 = .not.z3
        if(a0*b0.le.0.) go to 100
        z0 = .true.
        if(a0.lt.0.) z0 = .false.
        nz0 = .not.z0
        z2 = .true.
        if(a2.lt.0.) z2 = .false.
        nz2 = .not.z2
        z4 = .true.
        if(3.0*a3+a2.lt.0.) z4 = .false.
        nz4 = .not.z4
        if(.not.((z0.and.(nz1.and.(z3.or.z2.and.nz4).or.nz2.and.
     * z3.and.z4).or.nz0.and.(z1.and.(nz3.or.nz2.and.z4).or.z2.and.
     * nz3.and.nz4))))go to 200
c  find the zeros of ql(y).
 100    call fpcuro(a3,a2,a1,a0,y,j)
        if(j.eq.0) go to 200
c  find which zeros of pl(x) are zeros of s(x).
        do 150 i=1,j
          if(y(i).lt.0. .or. y(i).gt.1.0) go to 150
c  test whether the number of zeros of s(x) exceeds mest.
          if(m.ge.mest) go to 700
          m = m+1
          zero(m) = t(l)+h1*y(i)
 150    continue
 200    a0 = b0
        ah = bh
        z1 = z3
        nz1 = nz3
 300  continue
c  the zeros of s(x) are arranged in increasing order.
      if(m.lt.2) go to 800
      do 400 i=2,m
        j = i
 350    j1 = j-1
        if(j1.eq.0) go to 400
        if(zero(j).ge.zero(j1)) go to 400
        zz = zero(j)
        zero(j) = zero(j1)
        zero(j1) = zz
        j = j1
        go to 350
 400  continue
      j = m
      m = 1
      do 500 i=2,j
c   GH211017
c       if(zero(i).eq.zero(m)) go to 500
        if(abs(zero(i)-zero(m))<epsilon(zero(i))) go to 500
        m = m+1
        zero(m) = zero(i)
 500  continue
      go to 800
 700  ier = 1
 800  return
      end

max_stars_repo_name: harbik/dierckx-sys

max_issues_repo_name: harbik/dierckx-sys

max_forks_repo_name: harbik/dierckx-sys


text:
Davis Municipal Code

 Chapter 22 Motor Vehicles and Traffic

http://www.city.davis.ca.us/cmo/citycode/chapter.cfm?chapter22 Chapter 22 of the Davis Municipal Code.

Wiki Pages linked from here are created to discuss that section of the Code. See Davis Municipal Code for more information on link format if you want to create a page for discussion.

 22.01.0  IN GENERAL
  22.01.010 Definitions.
  22.01.020 Compliance with chapter.
  22.01.030 Enforcement generally.
  22.01.040 Directing of traffic.
  22.01.050 Applicability of chapter to persons riding or operating pushcarts, animals, etc.
  22.01.060 Applicability of chapter to operator of government vehicle.
  22.01.070 Exemptions as to certain vehicles.
  22.01.080 Report of accidents.
  22.01.090 Shrubbery obstructing visibility at intersection.
  22.01.100 Oneway streets and alleysSigns required.
  22.01.110 SameDesignated.
  22.01.120 Person under influence of alcoholic liquor or narcotics in or about vehicles.
  22.01.130 Clinging to vehicles prohibited.
  22.01.140 Trains not to block crossings.

 22.02.0  CITY TRAFFIC ENGINEER
  22.02.010 Office established; director of public works designated traffic engineer.
  22.02.020 Duties.

 22.03.0  DAVIS SAFETY ADVISORY COMMISSION
  22.03.010 Created.
  22.03.020 To perform all functions of former traffic committee.
  22.03.030 Purpose.
  22.03.040 Membership; appointment of members.
  22.03.050 Filling of vacancies.
  22.03.060 Automatic termination of appointment.
  22.03.070 Chairman.
  22.03.080 Meetings.
  22.03.090 Organization; quorum; procedure.
  22.03.100 Secretary.
  22.03.110 Functions.

 22.04.0  OPERATION OF VEHICLES
  22.04.010 Decrease of prima facie speed limits.
  22.04.020 Driving through, etc., railroad crossing gates, etc., prohibited.
  22.04.030 Driving through funeral processions.
  22.04.040 Driving within sidewalk area.
  22.04.050 Driving over new pavement or paint.
  22.04.060 Entering or leaving limited access roadway.
  22.04.070 Driving over fire hose.
  22.04.080 Certain vehicles prohibited within central traffic district; exception.
  22.04.090 Advertising vehicles.
  22.04.100 Horsedrawn vehicles.
  22.04.110 Truck routesRegulations as to use.
  22.04.120 SameDesignated.
  22.04.130 Commercial vehicles prohibited on certain streetsSigns required.
  22.04.140 SameStreets designated.

 22.05.0  TURNING AND BACKING MOVEMENTS
  22.05.010 Placement of turn markers; obedience to turn markers.
  22.05.020 Restricted turn signs.
  22.05.030 Backing regulations generally.

 22.06.0  TRAFFICCONTROL DEVICES
  22.06.010 Authority to install.
  22.06.020 Signs required for enforcement.
  22.06.030 Obedience required.
  22.06.040 Removal of parking enforcement markings.
  22.06.050 Lane markings.
  22.06.060 Distinctive roadway markings.
  22.06.070 Removal, relocation, etc.
  22.06.080 Traffic signals.

 22.07.0  SPECIAL STOPS
  22.07.010 Erection of stop signs.
  22.07.020 Through streets designated.
  22.07.030 Stop intersections designated.
  22.07.040 Yield intersectionsConduct of drivers.
  22.07.050 SameIntersections designated.
  22.07.060 Stop lanes designated.
  22.07.070 Vehicles emerging from alleys, driveways, etc.

 22.08.0  STOPPING, STANDING AND PARKING
  22.08.010 Applicability of article.
  22.08.020 Enforcement and prosecution.
  22.08.030 Restrictions generally.
  22.08.040 Stopping, standing, etc., in parkways prohibited.
  22.08.050 Use of streets for storage of vehicle prohibited.
  22.08.060 Use of streets for selling, repairing, washing, etc., of vehicles.
  22.08.070 Parallel parking on oneway streets.
  22.08.080 Angle parkingRegulations.
  22.08.090 SameStreets where permitted designated.
  22.08.100 Parking adjacent to schools.
  22.08.110 Commercial truck parking permitted.
  22.08.120 Parking on narrow streets.
  22.08.130 Parking on city propertyGenerally.
  22.08.140 Parking on private property.
  22.08.150 Municipal parking structureRestrictions on use.
  22.08.160 Enforcement within county parking lot.
  22.08.170 Parking on city streetsGenerally.
  22.08.180 Enforcement of parking restrictions within the school district parking lot at 526 B Street.
  22.08.190 Parking space markings.
  22.08.200 Emergency parking restrictions.
  22.08.210 Peddlers, etc., vehicles, carts, etc.
  22.08.220 Bicycle parking zones.
  22.08.230 Curb markings.
  22.08.240 Reduction of no parking restriction at fire hydrants.
  22.08.245 45Minute parking between 7a.m. and 4p.m. Monday through Friday (Saturdays, Sundays and Holidays excluded)
  22.08.250 Parking time zonesOne hour, 9:00 A.M. to 6:00 P.M.
  22.08.260 Parking time zonesOne hour, 8:00 A.M. to 10:00 P.M.
  22.08.270 Parking time zonesOne hour, 5:00 P.M. to 10:00 P.M.
  22.08.275 90Minute parking between 8 a.m. and 5 p.m. Monday through Friday, (Saturdays, Sundays and Holidays excluded)
  22.08.280 SameTwo hours, 9:00 A.M. to 6:00 P.M.
  22.08.285 90Minute parking between 8 a.m. and 5 p.m. Monday through Saturdays (Sundays and Holidays excluded)
  22.08.290 Parking time zonesTwo hours, 9:00 A.M. to noon, except Sundays and holidays.
  22.08.295 2Hour Parking Between 7 a.m. and 4 p.m. Monday Through Friday (Saturdays, Sundays and Holidays Excluded)
  22.08.300 SameParking timed zonesThree hours 8:00 A.M. to 6:00 P.M.
  22.08.305 Parking timed zones  Four hours 8:00 A.M. to 6:00 P.M.
  22.08.310 Parking prohibited at all times.
  22.08.320 Stopping, standing and parking prohibited7:00 A.M. to 6:00 P.M., except Saturdays, Sundays and holidays.
  22.08.330 Same2:00 A.M. to 6:00 A.M.
  22.08.340 Same2:00 A.M. to 6:00 A.M., except Saturdays, Sundays and holidays.
  22.08.350 Same2:00 A.M. to 6:00 A.M. on Tuesdays.
  22.08.360 Sec. 1675. Same2:00 A.M. to 6:00 A.M. on Thursdays.
  22.08.365 Same2:00 A.M. to 6:00 A.M. on Fridays.
  22.08.370 Same2:00 A.M. to 9:00 A.M. on Fridays.
  22.08.380 Same5:00 A.M. to 8:00 A.M. on Saturdays.
  22.08.390 Same2:00 A.M. to 9:00 A.M., except Saturdays, Sundays and holidays.
  22.08.400 Same9:00 A.M. to 5:00 P.M., except Saturdays, Sundays and holidays.
  22.08.410 Same7:45 A.M. to 8:45 A.M., school days only.
  22.08.420 Same8:00 A.M. to 5:00 P.M.
  22.08.430 Same1:30 P.M. to 3:30 P.M., school days only.
  22.08.440 Display of warning device by disabled commercial vehicle.

 22.09.0  LOADING ZONES
  22.09.010 Establishment and marking.
  22.09.020 Use generally.
  22.09.030 Unauthorized use prohibited.
  22.09.040 Standing in passenger loading zone.
  22.09.050 Bus stops.
  22.09.060 Taxicab stands.

 22.10.0  PEDESTRIANS
  22.10.010 Establishment of crosswalks.
  22.10.020 Use of crosswalk required in certain areas.
  22.10.030 Method of crossing roadway.
  22.10.040 Standing in roadways.

 22.11.0  PARKING PERMITS FOR PERSONS WITH IMPAIRED MOBILITY
  22.11.010 Statement of purpose.
  22.11.020 Definitions.
  22.11.030 Reserved parking space permitsGenerally.
  22.11.035 SameCrosshatched area for loading/unloading for disabled persons.
  22.11.037 SamePrivate propertySignage and accessibility of parking spaces for disabled persons.
  22.11.040 SameEligibility.
  22.11.050 SameDisplay.
  22.11.060 Enforcement.

 22.12.0  PREFERENTIAL PARKING
  22.12.010 Purpose.
  22.12.020 Findings.
  22.12.030 Definitions.
  22.12.040 Designation of preferential permit parking areas.
  22.12.050 Criteria for designation of preferential permit parking areas.
  22.12.060 Designation process.
  22.12.070 Issuance of permits.
  22.12.080 Visitor permits.
  22.12.085 Consumer permits.
  22.12.090 Posting of preferential permit parking area.
  22.12.100 Display of permits.
  22.12.110 Permit parking exemption.
  22.12.120 Application for and duration of permit.
  22.12.130 Permit fees.
  22.12.140 Penalty provisions.
  22.12.150 Revocation of permits.

 22.13.0  MOTORCYCLES AND MOTORDRIVEN CYCLES
  22.13.010 Findings.
  22.13.020 Scope.
  22.13.030 Definitions.
  22.13.040 Prohibitions.
  22.13.050 Exceptions.
  22.13.060 Conditions of permit.
  22.13.070 Regulations not exclusive.
  22.13.080 ViolationsPenalties.

 22.14.0  INTERSTATE TRUCK TERMINALS AND ACCESS
  22.14.010 FindingsIntent and purpose.
  22.14.020 Application process.
  22.14.030 Definitions.
  22.14.040 Cost recovery.
  22.14.050 Appeals process.
  22.14.060 Retrofitting ramps, intersections and roadways.

 22.15.0  TRANSPORTATION SYSTEMS MANAGEMENT REQUIREMENTS
  22.15.010 Purpose.
  22.15.020 Findings.
  22.15.030 Definitions.
  22.15.040 Requirements.
  22.15.050 Transportation management certificate (TMC).
  22.15.060 Transportation management plan (TMP).
  22.15.070 Enforcement.
  22.15.080 Appeals.

 22.16.0  ELECTRIC VEHICLES
  22.16.010 Findings
  22.16.020 Scope
  22.16.030 Definitions
  22.16.040 Electric Vehicle Parking Spaces  Generally
  22.16.050 Prohibitions
  22.16.060 Exception
  22.16.070 Noticing of Electric Vehicle Parking Space
  22.16.080 ViolationsPenalties
  22.16.090 Regulations not exclusive


text:
Im 16 years old and live here in Davis, CA. Currently attending DHS and cant wait to get out of davis.

my email is smithersc@ymail.com if you would like to contact me


text:
!	Set intercepts by a regression.

!	catobsrange maps observed to underlying categories.
!	dat is the array of responses.
!	numcat is the number of underlying categories per item.
!	numcatobs is the number of observed categories per item.
!	design is the design matrix.
!	obsweight is the vector of observation weights.
!	offset is the offset for beta.
!	tolsing is a tolerance for modified Cholesky decomposition.
!	gamma is the vector of starting value.
subroutine interceptstart(catobsrange,dat,numcat,numcatobs,design,obsweight,offset,tolsing,gamma)
implicit none
interface
	function chol(sym,tolsing)
!	chol is used to compute a modified Cholesky decomposition of the n by n matrix sym.
!	The singularity tolerance is tolsing.
		implicit none
		real(kind=8),intent(in)::sym(:,:),tolsing
	 	real(kind=8)::chol(size(sym,1),size(sym,1))
	end function chol
!		This program uses the modified Cholesky decomposition
!		in lu to solve the equation lusolve = b.
	function solve(lu,b)
		implicit none

		real(kind=8),intent(in)::lu(:,:)
		real(kind=8),intent(in)::b(:)
		real(kind=8)::solve(size(b))
		end function solve
end interface
integer,intent(in)::catobsrange(:,:),dat(:,:),numcat(:),numcatobs(:)
real(kind=8),intent(in)::design(:,:),obsweight(:),offset(:),tolsing
real(kind=8),intent(out)::gamma(:)
!	al is the allocation error flag.
!	cat looks at observed categories.
!	colsub is the number of columns in subdesign.
!	obs counts observations.
!	item counts items.
!	position and position1 count observed category groups.
!	sumnumcat is sum(numcat).
integer::al,cat,colsub,item,obs,position,position1,sumnumcat
!	gammas to use.
logical::usegamma(size(gamma))
!	ave is a average.
!	counts are marginal frequencies.
!	logs are adjusted log frequencies.
!	total is a total.
real(kind=8)::ave,counts(sum(numcat)),logs(sum(numcat)),total
!	solution is used for a solution associated with the subdesign.
!	subchol is used for a modified Cholesky decomposition associated with subdesign'subdesign.
!	subdesign is the portion of design used in computations of the gamma starting values for intercepts.
real(kind=8),allocatable::solution(:),subchol(:,:),subdesign(:,:)
counts=0.0_8
sumnumcat=sum(numcat)

!	Count marginal totals.
do obs=1,size(dat,2)
	position=1
	do item=1,size(numcat)
		if(dat(item,obs)>=0.and.dat(item,obs)<numcatobs(item))then
			cat=position+dat(item,obs)
			counts(catobsrange(1,cat):catobsrange(2,cat))=&
				counts(catobsrange(1,cat):catobsrange(2,cat))&
				+obsweight(obs)/(catobsrange(2,cat)-catobsrange(1,cat)+1.0_8)
		end if
		position=position+numcatobs(item)
	end do
end do

!	Obtain adjusted log frequencies.
position=1
do  item=1,size(dat,1)

	do cat=position,position+numcat(item)-1

		if(counts(cat)>0.0_8) then
			logs(cat)=log(counts(cat))-offset(cat)
		else
			logs(cat)=-log(2.0_8)-offset(cat)
		end if

	end do
	total=sum(counts(position:position+numcat(item)-1))
	if(total>0.0_8)then
		ave=dot_product(logs(position:position+numcat(item)-1), &
			counts(position:position+numcat(item)-1))/total
	else
		ave=0.0_8
	end if
	logs(position:position+numcat(item)-1)=logs(position:position+numcat(item)-1)-ave
	position=position+numcat(item)
end do
!	Find a matrix to use if possible.
usegamma=.false.

do position=1,size(gamma)
	if(maxval(abs(design(1:sumnumcat,position)))>0.0_8&
		.and.maxval(abs(design(sumnumcat+1:size(offset),position)))==0.0_8) &
		usegamma(position) = .true.
end do
!	Check for failure.
colsub=count(usegamma)

if(colsub==0)return
!	Set up submatrix.
allocate(subdesign(sumnumcat,colsub),subchol(colsub,colsub),solution(colsub),stat=al)
if(al/=0) stop "Allocation failed for calculation of starting values of intercepts."
position1=1
do position=1,size(gamma)
	if(usegamma(position))then
		subdesign(:,position1)=design(1:sumnumcat,position)
		position1=position1+1
	end if
end do

do position=1,colsub
	position1=1
	do  item=1,size(dat,1)
		total=sum(counts(position1:position1+numcat(item)-1))
		if(total>0.0_8)then
			ave=dot_product(subdesign(position1:position1+numcat(item)-1,position), &
				counts(position1:position1+numcat(item)-1))/total
		else
			ave=0.0_8
		end if
		subdesign(position1:position1+numcat(item)-1,position)=&
			subdesign(position1:position1+numcat(item)-1,position)-ave
		position1=position1+numcat(item)
	end do
	do position1=1,position
		subchol(position,position1)=sum(counts*subdesign(:,position)*subdesign(:,position1))
	end do
end do

subchol=chol(subchol,tolsing)
do position=1,colsub
	solution(position)=sum(subdesign(:,position)*counts*logs)
end do
solution=solve(subchol,solution)

position1=1
do position=1,size(gamma)
	if(usegamma(position))then
		gamma(position)=solution(position1)
		position1=position1+1
	end if
end do
return
end subroutine interceptstart

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


text:
Havanese is a husband and wife Graduate Students gradstudent pair.  We are long time voyeurs on daviswiki but only recently created a profile.

Howdy Havanese, and welcome to the wiki!  Its good to have voyeurs peek their faces in now and then.  http://cheshirehall.net/link/ws/smile.png   Users/JabberWokky Evan JabberWokky Edwards

20071218 09:58:34 nbsp Saw your note on the Planet Beach Tanning entry... you can change the name of any entry by clicking Edit and then (down at the bottom, next to the Save and Preview buttons), click the Rename button.  In general, most entries are named after the full name of the business or group (not the full legal name with Inc, or whatever... just the full formal name).  In cases where its almost always referred to by a shorter name (The Aggie is what everybody calls The California Aggie newspaper), a redirect is often placed in the shorter form. Users/JabberWokky

20071230 19:50:18 nbsp Hello...I am so new to this stuff.  It is really confusing!!  I made one small change but I love what you put together! Users/sarahd

20080110 18:51:42 nbsp HAVANESE!!!!!!  We are now open!  Come in!!!!!  Would you like to bring your brother in law also?  Our Sirius wont be installed until next week, but he should still like it!!! Users/sarahd


text:
Adam Thongsavat is a history student at UC Davis who is served as ASUCD President for the 20112012 term. He ran for ASUCD President on the BOLD slate on an exec ticket with Bree Rombi. He was an ASUCD Senate ASUCD Senator, who won as an independent in the Winter 2010 Election. He previously served as Chair of the ASUCD Elections Committee, was a Picnic Day Board of Directors Picnic Day Assistant Director, and a driver for Unitrans

 Campaign Statement

Hey Aggies! We are Adam Thongsavat and Bree Rombi and we are running to be your ASUCD President and Vice President. Adam is a 4th year History major from Merced, CA. He is a current Senator, former Elections Committee Chair, Picnic Day Assistant Director, and Unitrans driver. Bree is a 3rd year Communications and Spanish double major from Huntington Beach, CA. She was a former Senator, former President of Cuarto Leadership Counc...il, member of Outreach Assembly, Senator Pro Tempore, and a member of the UC Davis Varsity Bowling Team.

The message of BOLD is really simple: do things that will make a difference and have a lasting impact. UC Davis has been some of the most fun and positive years of our lives, despite the difficulties and imperfections, there is so much potential in our campus. Even though the challenges we face are great from fee hikes to budget cuts, our goal this year is to relentlessly pursue the betterment of student life and we have worked tirelessly since the day we started to make sure we live up to that purpose.

Together we have worked on things such as: Noise Ordinance Reform, funding for Campus Center for the Environment, funding for Aggie Pack, established a Renters Rights Committee, Picnic Day expansion, establishing composting in the MU, and help establish a new unitThe Pantry. We have had the best office hours record, collectively only missed one meeting, and ran the budget hearings.

We would like to dedicate our term to accomplish the following:

AN ACTION PLAN TO REACH THE UCWIDE 2020 SUSTAINABILITY MANDATE
We will collaborate with environmental leaders to create a master plan for Davis to meet this mandate and follow with an educational campaign about its significance.

PATHWAY TO CREATE STUDENT JOBS
We will work with our units to expand services and reform outdated units.

RESTORE AGGIE PRIDE AND SAVE PICNIC DAY
We will work with downtown Davis businesses to get students pumped up about game days and create a unified campus. We will also ensure Picnic Day continues to be a annual UC Davis tradition.

WORK TO MAKE AGGIETV INTO ONE OF THE COUNTRYS PREMIER STUDENTRUN COLLEGE BROADCAST STATIONS.

LOBBY THE STATE AND UNIVERSTIY AGAINST FEE INCREASES
This will include a lobby day at the capital. We will also collaborate with fellow students to educate individuals about the current budget crisis and possible solutions.

CREATE A SAFE CAMPUS ENVIRONMENT
We understand student health and wellness is one of the most important things in college and we plan to support initiatives to increase the well being of students, physical and mental.

If you have any questions please feel free to contact us!

Vote ADAM THONGSAVAT and BREE ROMBI for President and Vice President and Vote BOLD #16 for ASUCD Senate!

 Senate Campaign Statement

Dear Friends,

My name is Adam Thongsavat and I am a thirdyear history major from Merced, California.  I served as Elections Committee Chair and an Assistant Director for Picnic Day for the past 3 years.  I am currently a driver for Unitrans.

I am running for senate simply for one reason: I want to see results.  I am tired of repeated and failed platform proposals; I am tired of candidates who go to senate to preside instead of providing real leadership.  And I am tired of seeing a great campus which I love not realize its full potential.  It is time we do things that matter and will have lasting impacts long after we graduate.

The Smart Party Initiative I want to promote awareness and influence city and university policy that is more favorable towards students who chose to involve alcohol in their outings.  By working with our Student Police Relations Committee, UCDPD, DPD, and the city of Davis, I want to create a one stopshop of information including everything from individual rights in regards to hosting a party, DUI, BUI, Minor in possession, and drunk in public.  I will also strongly advocate for a more favorable and reasonable noise ordinance policy directed at student renters.

Picnic Day Expansion We should work with our programming unit Picnic Day to provide more fundraising opportunities at Picnic Day, which includes food booths and commercial tents.  I would like to see more funds injected into preweek activities which would include a community concert the night before Picnic Day.

Downtown Davis Business Partnership with Cafés and Restaurants I want to work vigorously towards a strong partnership with select downtown eateries to stay open late during high capacity days such as finals week and Picnic Day to provide late night student patrons with a place to grab a bite to eat and expand study space alternatives.  My goal is to launch a pilot program by spring quarter.

Establishment of a Renter’s Rights Committee I want to coinitiate a bipartisan effort with Raj Kumar to establish an ASUCD committee devoted to research and advocacy on rights of renters on and off campus.  Raj and I think it is imperative to have an effective body to represent nearly a third of City of Davis residents who are students.  Immediate issues would include branding the Davis Model Lease as well as establishing a framework for constant updates and revisions so it can actually be used.

I have had some great experiences at Davis the friendships I have made here are invaluable and I am so grateful for them.  That is why I am so inclined that now more than ever is the time for me to contribute to the betterment of student life it is the least I can do for a place and people which have fostered some of the most fun years of my life.

Be Bold.  Join our campaign because UC Davis can do better.  Thank you for your time and I love you UC Davis!

Please vote Adam Thongsavat #1 and Raj Kumar #2

Join my Facebook Page: http://www.facebook.com/home.php?#!/group.php?gid276811811343&refts

20100206 14:37:24 nbsp How do you feel about the proposed plan to annex the MUGA? Users/StevenDaubert


text:
      SUBROUTINE MB04BZ( JOB, COMPQ, N, A, LDA, DE, LDDE, B, LDB, FG,
     $                   LDFG, Q, LDQ, ALPHAR, ALPHAI, BETA, IWORK,
     $                   DWORK, LDWORK, ZWORK, LZWORK, BWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute the eigenvalues of a complex N-by-N skew-Hamiltonian/
C     Hamiltonian pencil aS - bH, with
C
C           (  A  D  )         (  B  F  )
C       S = (      H ) and H = (      H ).                           (1)
C           (  E  A  )         (  G -B  )
C
C     This routine computes the eigenvalues using an embedding to a real
C     skew-Hamiltonian/skew-Hamiltonian pencil aB_S - bB_T, defined as
C
C             (  Re(A)  -Im(A)  |  Re(D)  -Im(D)  )
C             (                 |                 )
C             (  Im(A)   Re(A)  |  Im(D)   Re(D)  )
C             (                 |                 )
C       B_S = (-----------------+-----------------) , and
C             (                 |      T       T  )
C             (  Re(E)  -Im(E)  |  Re(A )  Im(A ) )
C             (                 |      T       T  )
C             (  Im(E)   Re(E)  | -Im(A )  Re(A ) )
C                                                                    (2)
C             ( -Im(B)  -Re(B)  | -Im(F)  -Re(F)  )
C             (                 |                 )
C             (  Re(B)  -Im(B)  |  Re(F)  -Im(F)  )
C             (                 |                 )
C       B_T = (-----------------+-----------------) ,  T = i*H.
C             (                 |      T       T  )
C             ( -Im(G)  -Re(G)  | -Im(B )  Re(B ) )
C             (                 |      T       T  )
C             (  Re(G)  -Im(G)  | -Re(B ) -Im(B ) )
C
C     Optionally, if JOB = 'T', the pencil aB_S - bB_H (B_H = -i*B_T) is
C     transformed by a unitary matrix Q to the structured Schur form
C
C                ( BA  BD  )              ( BB  BF  )
C       B_Sout = (       H ) and B_Hout = (       H ),               (3)
C                (  0  BA  )              (  0 -BB  )
C
C     where BA and BB are upper triangular, BD is skew-Hermitian, and
C     BF is Hermitian. The embedding doubles the multiplicities of the
C     eigenvalues of the pencil aS - bH. Optionally, if COMPQ = 'C', the
C     unitary matrix Q is computed.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     JOB     CHARACTER*1
C             Specifies the computation to be performed, as follows:
C             = 'E': compute the eigenvalues only; S and H will not
C                    necessarily be transformed as in (3).
C             = 'T': put S and H into the forms in (3) and return the
C                    eigenvalues in ALPHAR, ALPHAI and BETA.
C
C     COMPQ   CHARACTER*1
C             Specifies whether to compute the unitary transformation
C             matrix Q, as follows:
C             = 'N':  Q is not computed;
C             = 'C':  compute the unitary transformation matrix Q.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the pencil aS - bH.  N >= 0, even.
C
C     A       (input/output) COMPLEX*16 array, dimension (LDA, K)
C             where K = N/2, if JOB = 'E', and K = N, if JOB = 'T'.
C             On entry, the leading N/2-by-N/2 part of this array must
C             contain the matrix A.
C             On exit, if JOB = 'T', the leading N-by-N part of this
C             array contains the upper triangular matrix BA in (3) (see
C             also METHOD). The strictly lower triangular part is not
C             zeroed, but it is preserved.
C             If JOB = 'E', this array is unchanged on exit.
C
C     LDA     INTEGER
C             The leading dimension of the array A.  LDA >= MAX(1, K).
C
C     DE      (input/output) COMPLEX*16 array, dimension
C                            (LDDE, MIN(K+1,N))
C             On entry, the leading N/2-by-N/2 lower triangular part of
C             this array must contain the lower triangular part of the
C             skew-Hermitian matrix E, and the N/2-by-N/2 upper
C             triangular part of the submatrix in the columns 2 to N/2+1
C             of this array must contain the upper triangular part of
C             the skew-Hermitian matrix D.
C             On exit, if JOB = 'T', the leading N-by-N part of this
C             array contains the skew-Hermitian matrix BD in (3) (see
C             also METHOD). The strictly lower triangular part of the
C             input matrix is preserved.
C             If JOB = 'E', this array is unchanged on exit.
C
C     LDDE    INTEGER
C             The leading dimension of the array DE.  LDDE >= MAX(1, K).
C
C     B       (input/output) COMPLEX*16 array, dimension (LDB, K)
C             On entry, the leading N/2-by-N/2 part of this array must
C             contain the matrix B.
C             On exit, if JOB = 'T', the leading N-by-N part of this
C             array contains the upper triangular matrix BB in (3) (see
C             also METHOD). The strictly lower triangular part is not
C             zeroed; the elements below the first subdiagonal of the
C             input matrix are preserved.
C             If JOB = 'E', this array is unchanged on exit.
C
C     LDB     INTEGER
C             The leading dimension of the array B.  LDB >= MAX(1, K).
C
C     FG      (input/output) COMPLEX*16 array, dimension
C                            (LDFG, MIN(K+1,N))
C             On entry, the leading N/2-by-N/2 lower triangular part of
C             this array must contain the lower triangular part of the
C             Hermitian matrix G, and the N/2-by-N/2 upper triangular
C             part of the submatrix in the columns 2 to N/2+1 of this
C             array must contain the upper triangular part of the
C             Hermitian matrix F.
C             On exit, if JOB = 'T', the leading N-by-N part of this
C             array contains the Hermitian matrix BF in (3) (see also
C             METHOD). The strictly lower triangular part of the input
C             matrix is preserved. The diagonal elements might have tiny
C             imaginary parts.
C             If JOB = 'E', this array is unchanged on exit.
C
C     LDFG    INTEGER
C             The leading dimension of the array FG.  LDFG >= MAX(1, K).
C
C     Q       (output) COMPLEX*16 array, dimension (LDQ, 2*N)
C             On exit, if COMPQ = 'C', the leading 2*N-by-2*N part of
C             this array contains the unitary transformation matrix Q
C             that reduced the matrices B_S and B_H to the form in (3).
C             However, if JOB = 'E', the reduction was possibly not
C             completed: the matrix B_H may have 2-by-2 diagonal blocks,
C             and the array Q returns the orthogonal matrix that
C             performed the partial reduction.
C             If COMPQ = 'N', this array is not referenced.
C
C     LDQ     INTEGER
C             The leading dimension of the array Q.
C             LDQ >= 1,           if COMPQ = 'N';
C             LDQ >= MAX(1, 2*N), if COMPQ = 'C'.
C
C     ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
C             The real parts of each scalar alpha defining an eigenvalue
C             of the pencil aS - bH.
C
C     ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
C             The imaginary parts of each scalar alpha defining an
C             eigenvalue of the pencil aS - bH.
C             If ALPHAI(j) is zero, then the j-th eigenvalue is real.
C
C     BETA    (output) DOUBLE PRECISION array, dimension (N)
C             The scalars beta that define the eigenvalues of the pencil
C             aS - bH.
C             Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
C             beta = BETA(j) represent the j-th eigenvalue of the pencil
C             aS - bH, in the form lambda = alpha/beta. Since lambda may
C             overflow, the ratios should not, in general, be computed.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (2*N+4)
C             On exit, IWORK(1) contains the number, q, of unreliable,
C             possibly inaccurate (pairs of) eigenvalues, and the
C             absolute values in IWORK(2), ..., IWORK(q+1) are their
C             indices, as well as of the corresponding 1-by-1 and 2-by-2
C             diagonal blocks of the arrays B and A on exit, if
C             JOB = 'T'. Specifically, a positive value is an index of
C             a real or purely imaginary eigenvalue, corresponding to a
C             1-by-1 block, while the absolute value of a negative entry
C             in IWORK is an index to the first eigenvalue in a pair of
C             consecutively stored eigenvalues, corresponding to a
C             2-by-2 block. Moreover, IWORK(q+2),..., IWORK(2*q+1)
C             contain pointers to the starting elements in DWORK where
C             each block pair is stored. Specifically, if IWORK(i+1) > 0
C             then DWORK(r) and DWORK(r+1) store corresponding diagonal
C             elements of T11 and S11, respectively, and if
C             IWORK(i+1) < 0, then DWORK(r:r+3) and DWORK(r+4:r+7) store
C             the elements of the block in T11 and S11, respectively
C             (see Section METHOD), where r = IWORK(q+1+i). Moreover,
C             IWORK(2*q+2) contains the number of the 1-by-1 blocks, and
C             IWORK(2*q+3) contains the number of the 2-by-2 blocks,
C             corresponding to unreliable eigenvalues. IWORK(2*q+4)
C             contains the total number t of the 2-by-2 blocks.
C             If INFO = 0, then q = 0, therefore IWORK(1) = 0.
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C             On exit, if INFO = 0 or INFO = 3, DWORK(1) returns the
C             optimal LDWORK, and DWORK(2) and DWORK(3) contain the
C             Frobenius norms of the matrices B_S and B_T. These norms
C             are used in the tests to decide that some eigenvalues are
C             considered as numerically unreliable. Moreover, DWORK(4),
C             ..., DWORK(3+2*s) contain the s pairs of values of the
C             1-by-1 diagonal elements of T11 and S11. The eigenvalue of
C             such a block pair is obtained from -i*T11(i,i)/S11(i,i).
C             Similarly, DWORK(4+2*s), ..., DWORK(3+2*s+8*t) contain the
C             t groups of pairs of 2-by-2 diagonal submatrices of T11
C             and S11, stored column-wise. The spectrum of such a block
C             pair is obtained from -i*ev, where ev are the eigenvalues
C             of (T11(i:i+1,i:i+1),S11(i:i+1,i:i+1)).
C             On exit, if INFO = -19, DWORK(1) returns the minimum value
C             of LDWORK.
C
C     LDWORK  INTEGER
C             The dimension of the array DWORK. If COMPQ = 'N',
C             LDWORK >= MAX( 3,  4*N*N + 3*N ), if   JOB = 'E';
C             LDWORK >= MAX( 3,  5*N*N + 3*N ), if   JOB = 'T';
C             LDWORK >= MAX( 3, 11*N*N + 2*N ), if COMPQ = 'C'.
C             For good performance LDWORK should be generally larger.
C
C             If LDWORK = -1, then a workspace query is assumed;
C             the routine only calculates the optimal size of the
C             DWORK array, returns this value as the first entry of
C             the DWORK array, and no error message related to LDWORK
C             is issued by XERBLA.
C
C     ZWORK   COMPLEX*16 array, dimension (LZWORK)
C             On exit, if INFO = 0, ZWORK(1) returns the optimal LZWORK.
C             On exit, if INFO = -21, ZWORK(1) returns the minimum value
C             of LZWORK.
C
C     LZWORK  INTEGER
C             The dimension of the array ZWORK.
C             LZWORK >= 1,       if JOB = 'E'; otherwise,
C             LZWORK >= 6*N + 4, if COMPQ = 'N';
C             LZWORK >= 8*N + 4, if COMPQ = 'C'.
C
C             If LZWORK = -1, then a workspace query is assumed;
C             the routine only calculates the optimal size of the
C             ZWORK array, returns this value as the first entry of
C             the ZWORK array, and no error message related to LZWORK
C             is issued by XERBLA.
C
C     BWORK   LOGICAL array, dimension (LBWORK)
C             LBWORK >= 0, if JOB = 'E';
C             LBWORK >= N, if JOB = 'T'.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0: succesful exit;
C             < 0: if INFO = -i, the i-th argument had an illegal value;
C             = 1: QZ iteration failed in the SLICOT Library routine
C                  MB04FD (QZ iteration did not converge or computation
C                  of the shifts failed);
C             = 2: QZ iteration failed in the LAPACK routine ZHGEQZ when
C                  trying to triangularize the 2-by-2 blocks;
C             = 3: warning: the pencil is numerically singular.
C
C     METHOD
C
C     First, T = i*H is set. Then, the embeddings, B_S and B_T, of the
C     matrices S and T, are determined and, subsequently, the SLICOT
C     Library routine MB04FD is applied to compute the structured Schur
C     form, i.e., the factorizations
C
C     ~        T  T         (  S11  S12  )
C     B_S = J Q  J  B_S Q = (          T ) and
C                           (   0   S11  )
C
C     ~        T  T         (  T11  T12  )
C     B_T = J Q  J  B_T Q = (          T ),
C                           (   0   T11  )
C
C     where Q is real orthogonal, S11 is upper triangular, and T11 is
C     upper quasi-triangular. If JOB = 'T', then the matrices above are
C                                                        ~
C     further transformed so that the 2-by-2 blocks in i*B_T are split
C     into 1-by-1 blocks.  If COMPQ = 'C', the transformations are
C     accumulated in the unitary matrix Q.
C     See also page 22 in [1] for more details.
C
C     REFERENCES
C
C     [1] Benner, P., Byers, R., Mehrmann, V. and Xu, H.
C         Numerical Computation of Deflating Subspaces of Embedded
C         Hamiltonian Pencils.
C         Tech. Rep. SFB393/99-15, Technical University Chemnitz,
C         Germany, June 1999.
C
C     NUMERICAL ASPECTS
C                                                               3
C     The algorithm is numerically backward stable and needs O(N )
C     complex floating point operations.
C
C     FURTHER COMMENTS
C
C     The returned eigenvalues are those of the pencil (-i*T11,S11),
C     where i is the purely imaginary unit.
C
C     If JOB = 'E', the returned matrix T11 is quasi-triangular. Note
C     that the off-diagonal elements of the 2-by-2 blocks of S11 are
C     zero by construction.
C
C     If JOB = 'T', the returned eigenvalues correspond to the diagonal
C     elements of BB and BA.
C
C     This routine does not perform any scaling of the matrices. Scaling
C     might sometimes be useful, and it should be done externally.
C
C     CONTRIBUTOR
C
C     V. Sima, Research Institute for Informatics, Bucharest, Romania,
C     Nov. 2011.
C
C     REVISIONS
C
C     V. Sima, July 2012, Sep. 2016, Nov. 2016, Apr. 2020.
C
C     KEYWORDS
C
C     Deflating subspace, embedded pencil, skew-Hamiltonian/Hamiltonian
C     pencil, structured Schur form.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ONE, THREE, ZERO
      PARAMETER          ( ONE = 1.0D+0, THREE = 3.0D+0, ZERO = 0.0D+0 )
      COMPLEX*16         CZERO, CONE, CIMAG
      PARAMETER          ( CZERO = ( 0.0D+0, 0.0D+0 ),
     $                      CONE = ( 1.0D+0, 0.0D+0 ),
     $                     CIMAG = ( 0.0D+0, 1.0D+0 ) )
C
C     .. Scalar Arguments ..
      CHARACTER          COMPQ, JOB
      INTEGER            INFO, LDA, LDB, LDDE, LDFG, LDQ, LDWORK,
     $                   LZWORK, N
C
C     .. Array Arguments ..
      LOGICAL            BWORK( * )
      INTEGER            IWORK( * )
      DOUBLE PRECISION   ALPHAI( * ), ALPHAR( * ), BETA( * ), DWORK( * )
      COMPLEX*16         A( LDA, * ), B( LDB, * ), DE( LDDE, * ),
     $                   FG( LDFG, * ), Q( LDQ, * ), ZWORK( * )
C
C     .. Local Scalars ..
      LOGICAL            LCMPQ, LQUERY, LTRI, UNREL
      CHARACTER*14       CMPQ, JOBF
      INTEGER            I, I1, IA, IB, IDE, IEV, IFG, IQ, IQ2, IQB, IS,
     $                   IW, IW1, IWRK, J, J1, J2, JM1, JP2, K, L, M,
     $                   MINDB, MINDW, MINZW, N2, NB, NC, NN, OPTDW,
     $                   OPTZW
      DOUBLE PRECISION   NRMBS, NRMBT
      COMPLEX*16         TMP
C
C     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
C
C     .. External Subroutines ..
      EXTERNAL           DCOPY, DLACPY, DSCAL, MB04FD, XERBLA, ZGEMM,
     $                   ZGEQRF, ZHGEQZ, ZLACPY, ZSCAL
C
C     .. Intrinsic Functions ..
      INTRINSIC          DBLE, DCMPLX, DIMAG, INT, MAX, MIN, MOD
C
C     .. Executable Statements ..
C
C     Decode the input arguments.
C
      M     = N/2
      NN    = N*N
      N2    = 2*N
      LTRI  = LSAME( JOB,   'T' )
      LCMPQ = LSAME( COMPQ, 'C' )
C
      IF( LTRI ) THEN
         K = N
      ELSE
         K = M
      END IF
C
      IF( N.EQ.0 ) THEN
         MINDW = 3
         MINZW = 1
      ELSE IF( LCMPQ ) THEN
         MINDB =  8*NN + N2
         MINDW = 11*NN + N2
         MINZW =  8*N  + 4
      ELSE
         MINDB =  4*NN + N2
         IF( LTRI ) THEN
            MINDW = 5*NN + 3*N
         ELSE
            MINDW = 4*NN + 3*N
         END IF
         MINZW = 6*N + 4
      END IF
C
      LQUERY = LDWORK.EQ.-1 .OR. LZWORK.EQ.-1
C
C     Test the input arguments.
C
      INFO = 0
      IF( .NOT.( LSAME( JOB, 'E' ) .OR. LTRI ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( LSAME( COMPQ, 'N' ) .OR. LCMPQ ) ) THEN
         INFO = -2
      ELSE IF( N.LT.0 .OR. MOD( N, 2 ).NE.0 ) THEN
         INFO = -3
      ELSE IF(  LDA.LT.MAX( 1, K ) ) THEN
         INFO = -5
      ELSE IF( LDDE.LT.MAX( 1, K ) ) THEN
         INFO = -7
      ELSE IF(  LDB.LT.MAX( 1, K ) ) THEN
         INFO = -9
      ELSE IF( LDFG.LT.MAX( 1, K ) ) THEN
         INFO = -11
      ELSE IF( LDQ.LT.1 .OR. ( LCMPQ .AND. LDQ.LT.N2 ) ) THEN
         INFO = -13
      ELSE IF( .NOT. LQUERY ) THEN
         IF( LDWORK.LT.MINDW ) THEN
            DWORK( 1 ) = MINDW
            INFO = -19
         ELSE IF( LZWORK.LT.MINZW ) THEN
            ZWORK( 1 ) = MINZW
            INFO = -21
         END IF
      END IF
C
      IF( INFO.NE.0) THEN
         CALL XERBLA( 'MB04BZ', -INFO )
         RETURN
      ELSE IF( N.GT.0 ) THEN
C
C        Compute optimal workspace.
C
         OPTZW = MINZW
         IF( LCMPQ ) THEN
            CMPQ = 'Initialize'
         ELSE
            CMPQ = 'No Computation'
         END IF
C
         IF( LTRI ) THEN
            JOBF = 'Triangularize'
            CALL ZGEQRF( N, N, ZWORK, N, ZWORK, ZWORK, -1, INFO )
            I  = INT( ZWORK( 1 ) )
            NB = MAX( I/N, 2 )
         ELSE
            JOBF = 'Eigenvalues'
         END IF
C
         IF( LQUERY ) THEN
            CALL MB04FD( JOBF, CMPQ, N2, DWORK, N, DWORK, N, DWORK, N,
     $                   DWORK, N, DWORK, N2, ALPHAI, ALPHAR, BETA,
     $                   IWORK, DWORK, -1, INFO )
            OPTDW = MAX( MINDW, MINDB + INT( DWORK( 1 ) ) )
            DWORK( 1 ) = OPTDW
            ZWORK( 1 ) = OPTZW
            RETURN
         END IF
      ELSE IF( LQUERY ) THEN
         DWORK( 1 ) = MINDW
         ZWORK( 1 ) = MINZW
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF( N.EQ.0 ) THEN
         IWORK( 1 ) = 0
         DWORK( 1 ) = THREE
         DWORK( 2 ) = ZERO
         DWORK( 3 ) = ZERO
         ZWORK( 1 ) = CONE
         RETURN
      END IF
C
C     Set up the embeddings of the matrices S and H.
C
C     Set the pointers for the inputs and outputs of MB04FD.
C     Real workspace:    need   4*N**2 + 2*N, if COMPQ = 'N';
C                               8*N**2 + 2*N, if COMPQ = 'C'.
C
      IQ = 1
      IF( LCMPQ ) THEN
         IA = IQ + N2*N2
      ELSE
         IA = 1
      END IF
C
      IDE  = IA  + NN
      IB   = IDE + NN + N
      IFG  = IB  + NN
      IWRK = IFG + NN + N
C
C     Build the embedding of A.
C
      IW = IA
      IS = IW + N*M
      DO 30 J = 1, M
         IW1 = IW
         DO 10 I = 1, M
            DWORK( IW ) = DBLE( A( I, J ) )
            IW = IW + 1
   10    CONTINUE
C
         DO 20 I = 1, M
            DWORK( IW ) =  DIMAG( A( I, J ) )
            DWORK( IS ) = -DWORK( IW )
            IW = IW + 1
            IS = IS + 1
   20    CONTINUE
         CALL DCOPY( M, DWORK( IW1 ), 1, DWORK( IS ), 1 )
         IS = IS + M
   30 CONTINUE
C
C     Build the embedding of D and E.
C
      IW = IDE
      DO 60 J = 1, M + 1
         DO 40 I = 1, M
            DWORK( IW ) = DBLE( DE( I, J ) )
            IW = IW + 1
   40    CONTINUE
C
         IW = IW + J - 1
         IS = IW
         DO 50 I = J, M
            DWORK( IW ) = DIMAG( DE( I, J ) )
            DWORK( IS ) = DWORK( IW )
            IW = IW + 1
            IS = IS + N
   50    CONTINUE
   60 CONTINUE
C
      IW1 = IW
      I1  = IW
      DO 80 J = 2, M + 1
         IS = I1
         I1 = I1 + 1
         DO 70 I = 1, J - 1
            DWORK( IW ) = -DIMAG( DE( I, J ) )
            DWORK( IS ) =  DWORK( IW )
            IW = IW + 1
            IS = IS + N
   70    CONTINUE
         IW = IW + N - J + 1
   80 CONTINUE
      CALL DLACPY( 'Full', M, M+1, DWORK( IDE ), N, DWORK( IW1-M ), N )
C
C     Build the embedding of B.
C
      IW = IB
      IS = IW + N*M
      DO 110 J = 1, M
         IW1 = IW
         DO 90 I = 1, M
            DWORK( IW ) = -DIMAG( B( I, J ) )
            IW = IW + 1
   90    CONTINUE
C
         DO 100 I = 1, M
            DWORK( IW ) =   DBLE( B( I, J ) )
            DWORK( IS ) = -DWORK( IW )
            IW = IW + 1
            IS = IS + 1
  100    CONTINUE
         CALL DCOPY( M, DWORK( IW1 ), 1, DWORK( IS ), 1 )
         IS = IS + M
  110 CONTINUE
C
C     Build the embedding of F and G.
C
      IW = IFG
      DO 140 J = 1, M + 1
         DO 120 I = 1, M
            DWORK( IW ) = -DIMAG( FG( I, J ) )
            IW = IW + 1
  120    CONTINUE
C
         IW = IW + J - 1
         IS = IW
         DO 130 I = J, M
            DWORK( IW ) =  DBLE( FG( I, J ) )
            DWORK( IS ) = DWORK( IW )
            IW = IW + 1
            IS = IS + N
  130    CONTINUE
  140 CONTINUE
C
      IW1 = IW
      I1  = IW
      DO 160 J = 2, M + 1
         IS = I1
         I1 = I1 + 1
         DO 150 I = 1, J - 1
            DWORK( IW ) = -DBLE( FG( I, J ) )
            DWORK( IS ) = DWORK( IW )
            IW = IW + 1
            IS = IS + N
  150    CONTINUE
         IW = IW + N - J + 1
  160 CONTINUE
      CALL DLACPY( 'Full', M, M+1, DWORK( IFG ), N, DWORK( IW1-M ), N )
C
C     STEP 1: Apply MB04FD to transform the extended pencil to real
C             skew-Hamiltonian/skew-Hamiltonian Schur form.
C
C     Real workspace:
C     need    4*N*N + 2*N + max(3,N),     if JOB = 'E' and COMPQ = 'N';
C             4*N*N + 2*N + max(3,N*N+N), if JOB = 'T' and COMPQ = 'N';
C            11*N*N + 2*N,                if               COMPQ = 'C'.
C     prefer larger.
C
      CALL MB04FD( JOBF, CMPQ, N2, DWORK( IA ), N, DWORK( IDE ), N,
     $             DWORK( IB ), N, DWORK( IFG ), N, DWORK( IQ ), N2,
     $             ALPHAI, ALPHAR, BETA, IWORK, DWORK( IWRK ),
     $             LDWORK-IWRK+1, INFO )
      IF( INFO.EQ.1 ) THEN
         RETURN
      ELSE IF( INFO.EQ.2 ) THEN
         INFO = 3
      END IF
      OPTDW = MAX( MINDW, MINDB + INT( DWORK( IWRK ) ) )
      NRMBS = DWORK( IWRK+1 )
      NRMBT = DWORK( IWRK+2 )
C
C     Scale the eigenvalues.
C
      CALL DSCAL( N, -ONE, ALPHAI, 1 )
C
      IF( LCMPQ ) THEN
C
C        Set the transformation matrix.
C
         IW = IQ
         DO 180 J = 1, N2
            DO 170 I = 1, N2
               Q( I, J ) = DCMPLX( DWORK( IW ) )
               IW = IW + 1
  170       CONTINUE
  180    CONTINUE
      END IF
C
C     Count the numbers of pairs of diagonal 1-by-1 and 2-by-2 blocks of
C     B and A with possibly unreliable eigenvalues.
C
      I  = IWORK( 1 )
      IS = 0
      IW = 0
      DO 190 J = 1, I
         IF( IWORK( J+1 ).GT.0 ) THEN
            IS = IS + 1
         ELSE IF( IWORK( J+1 ).LT.0 ) THEN
            IW = IW + 1
         END IF
  190 CONTINUE
C
      I = 2*I + 2
      IWORK( I )   = IS
      IWORK( I+1 ) = IW
C
      IF( LTRI ) THEN
C
C        Convert the results to complex datatype. D and F start in the
C        first column of DE and FG, respectively.
C
         IW = IA
         DO 210 J = 1, N
            DO 200 I = 1, J
               A( I, J ) = DCMPLX( DWORK( IW ) )
               IW = IW + 1
  200       CONTINUE
            IW = IW + N - J
  210    CONTINUE
C
         IW = IDE + N
         DO 230 J = 1, N
            DO 220 I = 1, J - 1
               DE( I, J ) = DCMPLX( DWORK( IW ) )
               IW = IW + 1
  220       CONTINUE
            DE( J, J ) = CZERO
            IW = IW + N - J + 1
  230    CONTINUE
C
         IW = IB
         DO 250 J = 1, N
            DO 240 I = 1, MIN( J + 1, N )
               B( I, J ) = DCMPLX( DWORK( IW ) )
               IW = IW + 1
  240       CONTINUE
            IW = IW + N - J - 1
  250    CONTINUE
C
         IW = IFG + N
         DO 270 J = 1, N
            DO 260 I = 1, J - 1
               FG( I, J ) = DCMPLX( DWORK( IW ) )
               IW = IW + 1
  260       CONTINUE
            FG( J, J ) = CZERO
            IW = IW + N - J + 1
  270    CONTINUE
C
      END IF
C
C     Count the number of diagonal 2-by-2 blocks of B and A.
C
      I1 = 0
      I  = 1
C     WHILE( I.LT.N ) DO
  280 CONTINUE
      IF ( I.LT.N ) THEN
         IF( ALPHAR( I ).NE.ZERO .AND. BETA( I ).NE.ZERO .AND.
     $       ALPHAR( I ).NE.ZERO ) THEN
            I1 = I1 + 1
            I  = I  + 2
         ELSE
            I  = I  + 1
         END IF
         GO TO 280
      END IF
C     END WHILE 280
C
      I = 2*IWORK( 1 ) + 4
      IWORK( I ) = I1
C
C     Save in DWORK the pairs of diagonal 1-by-1 and 2-by-2 blocks of
C     B and A. Also, save in IWORK the pointers to the starting element
C     of each pair corresponding to unreliable eigenvalues.
C
      CALL DCOPY( N, DWORK( IA ), N+1, DWORK( IFG ), 1 )
      DWORK( 1 ) = OPTDW
      DWORK( 2 ) = NRMBS
      DWORK( 3 ) = NRMBT
C
      K  = 4
      I  = 1
      IW = IWORK( 1 )
      J  = 1
      L  = ( N - 2*I1 )*2 + K
C     WHILE( I.LE.N ) DO
  290 CONTINUE
      IF ( I.LE.N ) THEN
         IF ( J.LE.IW )
     $      UNREL = I.EQ.ABS( IWORK( J+1 ) )
         IF( ALPHAR( I ).NE.ZERO .AND. BETA( I ).NE.ZERO .AND.
     $       ALPHAI( I ).NE.ZERO ) THEN
            IF ( UNREL ) THEN
               J = J + 1
               IWORK( J+IW ) = L
               UNREL = .FALSE.
            END IF
            CALL DLACPY( 'Full', 2, 2, DWORK( IB+(I-1)*(N+1) ), N,
     $                   DWORK( L ), 2 )
            L = L + 4
            CALL DCOPY( 2, DWORK( IFG+I-1 ), 1, DWORK( L ), 3 )
            DWORK( L+1 ) = ZERO
            DWORK( L+2 ) = ZERO
            L = L + 4
            I = I + 2
         ELSE
            IF ( UNREL ) THEN
               J = J + 1
               IWORK( J+IW ) = K
               UNREL = .FALSE.
            END IF
            DWORK( K   ) = DWORK( IB+(I-1)*(N+1) )
            DWORK( K+1 ) = DWORK( IFG+I-1 )
            K = K + 2
            I = I + 1
         END IF
         GO TO 290
      END IF
C     END WHILE 290
C
      IF( .NOT.LTRI ) THEN
C
C        Return.
C
         ZWORK( 1 ) = OPTZW
         RETURN
      END IF
C
C     Triangularize the 2-by-2 diagonal blocks in B using the complex
C     version of the QZ algorithm.
C
C     Set up pointers on the outputs of ZHGEQZ.
C     A block algorithm is used for large N.
C
      IQ2  = 1
      IEV  = 5
      IQ   = 9
      IWRK = IQ + 4*( N - 1 )
C
      J  = 1
      J1 = 1
      J2 = MIN( N, NB )
C     WHILE( J.LT.N ) DO
  300 CONTINUE
      IF( J.LT.N ) THEN
         IF( B( J+1, J ).NE.CZERO ) THEN
C
C           Triangularization step.
C           Complex workspace:    need   4*N + 6.
C           Real    workspace:    need   4*N + 5.
C
            NC  = MAX( J2-J-1, 0 )
            JM1 = MAX( J-1, 1 )
            JP2 = MIN( J+2, N )
            TMP         =  A( J+1, J )
            A( J+1, J ) = CZERO
            CALL ZHGEQZ( 'Schur Form', 'Initialize', 'Initialize', 2, 1,
     $                   2, B( J, J ), LDB, A( J, J ), LDA,
     $                   ZWORK( IEV ), ZWORK( IEV+2 ), ZWORK( IQ ), 2,
     $                   ZWORK( IQ2 ), 2, ZWORK( IWRK ), LZWORK-IWRK+1,
     $                   DWORK( 4*N+4 ), IW )
            A( J+1, J ) = TMP
            IF( IW.GT.0 ) THEN
               INFO = 2
               RETURN
            END IF
C
C           Update A, DE, B, and FG.
C           Complex workspace:    need   6*N + 4.
C
C           Update A.
C
            CALL ZGEMM(  'No Transpose', 'No Transpose', J-1, 2, 2,
     $                   CONE, A( 1, J ), LDA, ZWORK( IQ2 ), 2, CZERO,
     $                   ZWORK( IWRK ), JM1 )
            CALL ZLACPY( 'Full', J-1, 2, ZWORK( IWRK ), JM1, A( 1, J ),
     $                   LDA )
            CALL ZGEMM(  'Conjugate Transpose', 'No Transpose', 2, NC,
     $                   2, CONE, ZWORK( IQ ), 2, A( J, JP2 ), LDA,
     $                   CZERO, ZWORK( IWRK ), 2 )
            CALL ZLACPY( 'Full', 2, NC, ZWORK( IWRK ), 2, A( J, JP2 ),
     $                   LDA )
C
C           Update DE.
C
            TMP          =  DE( J+1, J )
            DE( J+1, J ) = -DE( J, J+1 )
            CALL ZGEMM(  'No Transpose', 'No Transpose', J+1, 2, 2,
     $                   CONE, DE( 1, J ), LDDE, ZWORK( IQ ), 2, CZERO,
     $                   ZWORK( IWRK ), J+1 )
            CALL ZLACPY( 'Full', J+1, 2, ZWORK( IWRK ), J+1, DE( 1, J ),
     $                   LDDE )
            CALL ZGEMM(  'Conjugate Transpose', 'No Transpose', 2,
     $                   J2-J+1, 2, CONE, ZWORK( IQ ), 2, DE( J, J ),
     $                   LDDE, CZERO, ZWORK( IWRK ), 2 )
            CALL ZLACPY( 'Full', 2, J2-J+1, ZWORK( IWRK ), 2,
     $                   DE( J, J ), LDDE )
            DE( J+1, J ) = TMP
C
C           Update B.
C
            CALL ZGEMM(  'No Transpose', 'No Transpose', J-1, 2, 2,
     $                   CONE, B( 1, J ), LDB, ZWORK( IQ2 ), 2, CZERO,
     $                   ZWORK( IWRK ), JM1 )
            CALL ZLACPY( 'Full', J-1, 2, ZWORK( IWRK ), JM1, B( 1, J ),
     $                   LDB )
            CALL ZGEMM(  'Conjugate Transpose', 'No Transpose', 2, NC,
     $                   2, CONE, ZWORK( IQ ), 2, B( J, JP2 ), LDB,
     $                   CZERO, ZWORK( IWRK ), 2 )
            CALL ZLACPY( 'Full', 2, NC, ZWORK( IWRK ), 2, B( J, JP2 ),
     $                   LDB )
C
C           Update FG.
C
            TMP          =  FG( J+1, J )
            FG( J+1, J ) = -FG( J, J+1 )
            CALL ZGEMM(  'No Transpose', 'No Transpose', J+1, 2, 2,
     $                   CONE, FG( 1, J ), LDFG, ZWORK( IQ ), 2, CZERO,
     $                   ZWORK( IWRK ), J+1 )
            CALL ZLACPY( 'Full', J+1, 2, ZWORK( IWRK ), J+1, FG( 1, J ),
     $                   LDFG )
            CALL ZGEMM(  'Conjugate Transpose', 'No Transpose', 2,
     $                   J2-J+1, 2, CONE, ZWORK( IQ ), 2, FG( J, J ),
     $                   LDFG, CZERO, ZWORK( IWRK ), 2 )
            CALL ZLACPY( 'Full', 2, J2-J+1, ZWORK( IWRK ), 2,
     $                   FG( J, J ), LDFG )
            FG( J+1, J ) = TMP
C
            IF( LCMPQ ) THEN
C
C              Update Q.
C              Complex workspace:    need   8*N + 4.
C
               CALL ZGEMM(  'No Transpose', 'No Transpose', N2, 2, 2,
     $                      CONE, Q( 1, J ), LDQ, ZWORK( IQ2 ), 2,
     $                      CZERO, ZWORK( IWRK ), N2 )
               CALL ZLACPY( 'Full', N2, 2, ZWORK( IWRK ), N2, Q( 1, J ),
     $                      LDQ )
               CALL ZGEMM(  'No Transpose', 'No Transpose', N2, 2, 2,
     $                      CONE, Q( 1, N+J ), LDQ, ZWORK( IQ ), 2,
     $                      CZERO, ZWORK( IWRK ), N2 )
               CALL ZLACPY( 'Full', N2, 2, ZWORK( IWRK ), N2,
     $                      Q( 1, N+J ), LDQ )
            END IF
C
            BWORK( J ) = .TRUE.
            J  = J  + 2
            IQ = IQ + 4
         ELSE
            BWORK( J )  = .FALSE.
            B( J+1, J ) = CZERO
            J = J + 1
         END IF
C
         IF( J.GE.J2 ) THEN
            J1 = J2 + 1
            J2 = MIN( N, J1 + NB - 1 )
            NC = J2 - J1 + 1
C
C           Update the columns J1 to J2 of A, DE, B, and FG for previous
C           transformations.
C
            I   = 1
            IQB = 9
C           WHILE( I.LT.J ) DO
  310       CONTINUE
            IF( I.LT.J ) THEN
               IF( BWORK( I ) ) THEN
                  CALL ZGEMM(  'Conjugate Transpose', 'No Transpose', 2,
     $                         NC, 2, CONE, ZWORK( IQB ), 2, A( I, J1 ),
     $                         LDA, CZERO, ZWORK( IWRK ), 2 )
                  CALL ZLACPY( 'Full', 2, NC, ZWORK( IWRK ), 2,
     $                         A( I, J1 ), LDA )
C
                  CALL ZGEMM(  'Conjugate Transpose', 'No Transpose', 2,
     $                         NC, 2, CONE, ZWORK( IQB ), 2,
     $                         DE( I, J1 ), LDDE, CZERO, ZWORK( IWRK ),
     $                         2 )
                  CALL ZLACPY( 'Full', 2, NC, ZWORK( IWRK ), 2,
     $                         DE( I, J1 ), LDDE )
C
                  CALL ZGEMM(  'Conjugate Transpose', 'No Transpose', 2,
     $                         NC, 2, CONE, ZWORK( IQB ), 2, B( I, J1 ),
     $                         LDB, CZERO, ZWORK( IWRK ), 2 )
                  CALL ZLACPY( 'Full', 2, NC, ZWORK( IWRK ), 2,
     $                         B( I, J1 ), LDB )
C
                  CALL ZGEMM(  'Conjugate Transpose', 'No Transpose', 2,
     $                         NC, 2, CONE, ZWORK( IQB ), 2,
     $                         FG( I, J1 ), LDFG, CZERO, ZWORK( IWRK ),
     $                         2 )
                  CALL ZLACPY( 'Full', 2, NC, ZWORK( IWRK ), 2,
     $                         FG( I, J1 ), LDFG )
                  IQB = IQB + 4
C
                  I = I + 2
               ELSE
                  I = I + 1
               END IF
               GO TO 310
            END IF
C           END WHILE 310
         END IF
         GO TO 300
      END IF
C     END WHILE 300
C
C     Scale B and FG by -i.
C
      DO 320 I = 1, N
         CALL ZSCAL( I, -CIMAG, B( 1, I ), 1 )
  320 CONTINUE
C
      DO 330 I = 1, N
         CALL ZSCAL( I, -CIMAG, FG( 1, I ), 1 )
  330 CONTINUE
C
      ZWORK( 1 ) = OPTZW
      RETURN
C *** Last line of MB04BZ ***
      END

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference


text:
C    @(#)bldbsekv.f	20.3 2/13/96
	integer function bldbsekv (basekv)
	real basekv
c
c	Function bldbsekv installs unique basekv into basekvs(*)
c
      	include 'ipfinc/bsekvhsh.inc'

      	integer        h, bsekvhsh, p
      	logical        status

      	status = (numbases .le. MAXBASEKVS)

      	if (status) then
           h = bsekvhsh (basekv)
           p = htable_k(h)
           do while (p .gt. 0)         !search for existing entities
              if (basekv .ne. basekvs(p)) then
                 p = nextptr_k(p)
              else
                 p = -p
              endif
           enddo
           if (p .eq. 0) then
              numbases = numbases + 1
              nextptr_k(numbases) = htable_k(h)
              htable_k(h) = numbases
              basekvs(numbases) = basekv
              bldbsekv = numbases
           else
              bldbsekv = -p         !duplicate entity flagged!
           endif
        else
           bldbsekv = 0             !overflow flagged!
        endif
        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