text:
! ================================================================================================ !
!  Collimation Common Variables
! ================================================================================================ !
module coll_common

  use floatPrecision
  use numerical_constants, only : zero

  implicit none

  ! Logical Flags
  logical, save :: coll_debug         = .true.
  logical, save :: dowrite_impact     = .false.
  logical, save :: dowrite_dist       = .false.
  logical, save :: dowrite_secondary  = .false.
  logical, save :: dowrite_amplitude  = .false.
  logical, save :: dowrite_tracks     = .false.
  logical, save :: dowrite_efficiency = .false.
  logical, save :: dowrite_crycoord   = .false.
  logical, save :: coll_hasCrystal    = .false.

  ! Various Variables
  integer, save :: rnd_seed   = 0

  ! Collimation Particle Arrays
  real(kind=fPrec), allocatable, save :: rcx(:)
  real(kind=fPrec), allocatable, save :: rcxp(:)
  real(kind=fPrec), allocatable, save :: rcy(:)
  real(kind=fPrec), allocatable, save :: rcyp(:)
  real(kind=fPrec), allocatable, save :: rcp(:)
  real(kind=fPrec), allocatable, save :: rcs(:)

  ! Process index for interaction with crystals
  integer, allocatable, save :: cry_proc(:)
  integer, allocatable, save :: cry_proc_prev(:)
  integer, allocatable, save :: cry_proc_tmp(:)

  ! Pencil Beam
  integer,          save :: ipencil       = 0
  integer,          save :: pencil_distr  = 0
  real(kind=fPrec), save :: pencil_offset = zero
  real(kind=fPrec), save :: pencil_rmsx   = zero
  real(kind=fPrec), save :: pencil_rmsy   = zero
  real(kind=fPrec), save :: xp_pencil0    = zero
  real(kind=fPrec), save :: yp_pencil0    = zero
  real(kind=fPrec), allocatable, save :: x_pencil(:)
  real(kind=fPrec), allocatable, save :: y_pencil(:)
  real(kind=fPrec), allocatable, save :: xp_pencil(:)
  real(kind=fPrec), allocatable, save :: yp_pencil(:)
  real(kind=fPrec), allocatable, save :: pencil_dx(:)

  ! Other Arrays
  integer,          allocatable, save :: cn_impact(:)
  integer,          allocatable, save :: cn_absorbed(:)
  real(kind=fPrec), allocatable, save :: caverage(:)
  real(kind=fPrec), allocatable, save :: csigma(:)
  real(kind=fPrec), allocatable, save :: gap_rms_error(:)
  real(kind=fPrec), allocatable, save :: csum(:)
  real(kind=fPrec), allocatable, save :: csqsum(:)

  ! Output File Names
  character(len=12), parameter :: coll_survivalFile   = "survival.dat"
  character(len=12), parameter :: coll_gapsFile       = "collgaps.dat"
  character(len=10), parameter :: coll_impactFile     = "impact.dat"
  character(len=11), parameter :: coll_tracksFile     = "tracks2.dat"
  character(len=17), parameter :: coll_positionsFile  = "CollPositions.dat"
  character(len=20), parameter :: coll_pencilFile     = "pencilbeam_distr.dat"
  character(len=16), parameter :: coll_ellipseFile    = "coll_ellipse.dat"
  character(len=15), parameter :: coll_allImpactFile  = "all_impacts.dat"
  character(len=19), parameter :: coll_allAbsorbFile  = "all_absorptions.dat"
  character(len=16), parameter :: coll_scatterFile    = "Coll_Scatter.dat"
  character(len=16), parameter :: coll_fstImpactFile  = "FirstImpacts.dat"
  character(len=17), parameter :: coll_flukImpFile    = "FLUKA_impacts.dat"
  character(len=21), parameter :: coll_flukImpAllFile = "FLUKA_impacts_all.dat"
  character(len=17), parameter :: coll_sigmaSetFile   = "sigmasettings.out"
  character(len=16), parameter :: coll_settingsFile   = "collsettings.dat"
  character(len=16), parameter :: coll_jawProfileFile = "jaw_profiles.dat"
  character(len=13), parameter :: coll_ampFile        = "amplitude.dat"
  character(len=17), parameter :: coll_orbitCheckFile = "orbitchecking.dat"
  character(len=16), parameter :: coll_summaryFile    = "coll_summary.dat"
  character(len=14), parameter :: coll_efficFile      = "efficiency.dat"
  character(len=19), parameter :: coll_efficDPFile    = "efficiency_dpop.dat"
  character(len=17), parameter :: coll_effic2DFile    = "efficiency_2d.dat"
  character(len=16), parameter :: coll_cryEntFile     = "cry_entrance.dat"
  character(len=12), parameter :: coll_cryExitFile    = "cry_exit.dat"
  character(len=19), parameter :: coll_cryInterFile   = "cry_interaction.dat"

  ! Output File Units
  integer, save :: outlun              = -1
  integer, save :: coll_survivalUnit   = -1
  integer, save :: coll_gapsUnit       = -1
  integer, save :: coll_impactUnit     = -1
  integer, save :: coll_tracksUnit     = -1
  integer, save :: coll_positionsUnit  = -1
  integer, save :: coll_pencilUnit     = -1
  integer, save :: coll_ellipseUnit    = -1
  integer, save :: coll_allImpactUnit  = -1
  integer, save :: coll_allAbsorbUnit  = -1
  integer, save :: coll_scatterUnit    = -1
  integer, save :: coll_fstImpactUnit  = -1
  integer, save :: coll_flukImpUnit    = -1
  integer, save :: coll_flukImpAllUnit = -1
  integer, save :: coll_sigmaSetUnit   = -1
  integer, save :: coll_settingsUnit   = -1
  integer, save :: coll_jawProfileUnit = -1
  integer, save :: coll_ampUnit        = -1
  integer, save :: coll_orbitCheckUnit = -1
  integer, save :: coll_summaryUnit    = -1
  integer, save :: coll_efficUnit      = -1
  integer, save :: coll_efficDPUnit    = -1
  integer, save :: coll_effic2DUnit    = -1
  integer, save :: coll_cryEntUnit     = -1
  integer, save :: coll_cryExitUnit    = -1
  integer, save :: coll_cryInterUnit   = -1

contains

subroutine coll_expandArrays(npart_new)

  use mod_alloc
  use numerical_constants

  integer, intent(in) :: npart_new

  call alloc(rcx,  npart_new, zero, "rcx")
  call alloc(rcxp, npart_new, zero, "rcxp")
  call alloc(rcy,  npart_new, zero, "rcy")
  call alloc(rcyp, npart_new, zero, "rcyp")
  call alloc(rcp,  npart_new, zero, "rcp")
  call alloc(rcs,  npart_new, zero, "rcs")

  call alloc(cry_proc, npart_new, -1, "cry_proc")
  call alloc(cry_proc_prev, npart_new, -1, "cry_proc_prev")
  call alloc(cry_proc_tmp, npart_new, -1, "cry_proc_tmp")

end subroutine coll_expandArrays

subroutine coll_expandNColl(nColl)

  use mod_alloc
  use numerical_constants

  integer, intent(in) :: nColl

  call alloc(cn_impact,     nColl, 0,    "cn_impact")
  call alloc(cn_absorbed,   nColl, 0,    "cn_absorbed")
  call alloc(caverage,      nColl, zero, "caverage")
  call alloc(csigma,        nColl, zero, "csigma")
  call alloc(gap_rms_error, nColl, zero, "gap_rms_error")
  call alloc(csum,          nColl, zero, "csum")
  call alloc(csqsum,        nColl, zero, "csqsum")
  call alloc(x_pencil,      nColl, zero, "x_pencil")
  call alloc(y_pencil,      nColl, zero, "y_pencil")
  call alloc(xp_pencil,     nColl, zero, "xp_pencil")
  call alloc(yp_pencil,     nColl, zero, "yp_pencil")
  call alloc(pencil_dx,     nColl, zero, "pencil_dx")

end subroutine coll_expandNColl

end module coll_common

! ================================================================================================ !
!  Collimator Material Variables
! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!  Cross section inputs and material property database
! ================================================================================================ !
module coll_materials

  use floatPrecision
  use numerical_constants, only : zero, one, c1e10

  implicit none

  integer, parameter :: nmat  = 16 ! Total number of materials
  integer, parameter :: nrmat = 14 ! Number of real materials

  ! pp cross-sections and parameters for energy dependence
  real(kind=fPrec), parameter :: pptref = 0.04_fPrec
  real(kind=fPrec), parameter :: freeco = 1.618_fPrec

  ! Collimator Material Arrays
  character(4),     public, save :: colmats(nmat)   ! Material Names
  real(kind=fPrec), public, save :: exenergy(nmat)  ! Mean excitation energy [GeV]
  real(kind=fPrec), public, save :: anuc(nmat)      ! Atomic mass
  real(kind=fPrec), public, save :: zatom(nmat)     ! Atomic Z
  real(kind=fPrec), public, save :: rho(nmat)       ! Density
  real(kind=fPrec), public, save :: emr(nmat)       ! Nuclear radius
  real(kind=fPrec), public, save :: hcut(nmat)      ! T cut (upper)
  real(kind=fPrec), public, save :: radl(nmat)      ! Remaining length
  real(kind=fPrec), public, save :: bnref(nmat)     ! Nuclear elastic slope from Schiz et al., PRD 21(3010)1980
  real(kind=fPrec), public, save :: csect(0:5,nmat) ! Cross section
  real(kind=fPrec), public, save :: xintl(nmat)     ! Interaction length
  real(kind=fPrec), public, save :: bn(nmat)        ! Nuclear elastic related
  real(kind=fPrec), public, save :: freep(nmat)     ! Number of nucleons involved in single scattering
  real(kind=fPrec), public, save :: cgen(200,nmat)  ! Used by FUNLUX / Rutherford routine

  ! All cross-sections are in barns. Nuclear values from RPP at 20 GeV
  ! Coulomb is integerated above t=tLcut[Gev2] (+-1% out Gauss mcs)

  ! In Cs and CsRef,1st index: Cross-sections for processes
  ! 0:Total, 1:absorption, 2:nuclear elastic, 3:pp or pn elastic
  ! 4:Single Diffractive pp or pn, 5:Coulomb for t above mcs

  ! Claudia 2013: updated cross section values. Unit: Barn. New 2013:
  real(kind=fPrec), public, parameter :: csref(0:5,nmat) = reshape([ &
    [0.271_fPrec, 0.192_fPrec, zero, zero, zero, 0.0035e-2_fPrec], & ! BE
    [0.643_fPrec, 0.418_fPrec, zero, zero, zero, 0.0340e-2_fPrec], & ! AL
    [1.253_fPrec, 0.769_fPrec, zero, zero, zero, 0.1530e-2_fPrec], & ! CU
    [2.765_fPrec, 1.591_fPrec, zero, zero, zero, 0.7680e-2_fPrec], & ! W
    [3.016_fPrec, 1.724_fPrec, zero, zero, zero, 0.9070e-2_fPrec], & ! PB
    [0.337_fPrec, 0.232_fPrec, zero, zero, zero, 0.0076e-2_fPrec], & ! C
    [0.337_fPrec, 0.232_fPrec, zero, zero, zero, 0.0076e-2_fPrec], & ! C2
    [0.664_fPrec, 0.430_fPrec, zero, zero, zero, 0.0390e-2_fPrec], & ! Si
    [1.388_fPrec, 0.844_fPrec, zero, zero, zero, 0.1860e-2_fPrec], & ! Ge
    [0.362_fPrec, 0.247_fPrec, zero, zero, zero, 0.0094e-2_fPrec], & ! MoGR
    [0.572_fPrec, 0.370_fPrec, zero, zero, zero, 0.0279e-2_fPrec], & ! CuCD
    [1.713_fPrec, 1.023_fPrec, zero, zero, zero, 0.2650e-2_fPrec], & ! Mo
    [1.246_fPrec, 0.765_fPrec, zero, zero, zero, 0.1390e-2_fPrec], & ! Glid
    [2.548_fPrec, 1.473_fPrec, zero, zero, zero, 0.5740e-2_fPrec], & ! Iner
    [       zero,        zero, zero, zero, zero,            zero], & ! VA
    [       zero,        zero, zero, zero, zero,            zero]  & ! BL
  ], shape=[6,nmat])

  ! Cprob to choose an interaction in iChoix
  real(kind=fPrec), public, save :: cprob(0:5,nmat) = reshape([ &
    [zero, zero, zero, zero, zero, one], & ! BE
    [zero, zero, zero, zero, zero, one], & ! AL
    [zero, zero, zero, zero, zero, one], & ! CU
    [zero, zero, zero, zero, zero, one], & ! W
    [zero, zero, zero, zero, zero, one], & ! PB
    [zero, zero, zero, zero, zero, one], & ! C
    [zero, zero, zero, zero, zero, one], & ! C2
    [zero, zero, zero, zero, zero, one], & ! Si
    [zero, zero, zero, zero, zero, one], & ! Ge
    [zero, zero, zero, zero, zero, one], & ! MoGR
    [zero, zero, zero, zero, zero, one], & ! CuCD
    [zero, zero, zero, zero, zero, one], & ! Mo
    [zero, zero, zero, zero, zero, one], & ! Glid
    [zero, zero, zero, zero, zero, one], & ! Iner
    [zero, zero, zero, zero, zero, one], & ! VA
    [zero, zero, zero, zero, zero, one]  & ! BL
  ], shape=[6,nmat])

  ! Electron density and plasma energy
  real(kind=fPrec), public, save :: edens(nmat) = zero
  real(kind=fPrec), public, save :: pleng(nmat) = zero

contains

! ================================================================================================ !
!  Initialise Material Database Values
!  V.K. Berglyd Olsen, BE-ABP-HSS
!  Created: 2019-09-19
!  Updated: 2019-09-19
!  Rewritten to make it easier to add new materials.
! ================================================================================================ !
subroutine collmat_init

  use crcoall

  integer iMat

  iMat = 1
  colmats(iMat)  = "BE"
  exenergy(iMat) =  63.7e-9_fPrec
  anuc(iMat)     =   9.01_fPrec
  zatom(iMat)    =   4.00_fPrec
  rho(iMat)      =   1.848_fPrec
  emr(iMat)      =   0.22_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.353_fPrec
  bnref(iMat)    =  74.7_fPrec

  iMat = iMat + 1
  colmats(iMat)  = "AL"
  exenergy(iMat) = 166.0e-9_fPrec
  anuc(iMat)     =  26.98_fPrec
  zatom(iMat)    =  13.00_fPrec
  rho(iMat)      =   2.70_fPrec
  emr(iMat)      =   0.302_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.089_fPrec
  bnref(iMat)    = 120.3_fPrec

  ! Copper
  iMat = iMat + 1
  colmats(iMat)  = "CU"
  exenergy(iMat) = 322.0e-9_fPrec
  anuc(iMat)     =  63.55_fPrec
  zatom(iMat)    =  29.00_fPrec
  rho(iMat)      =   8.96_fPrec
  emr(iMat)      =   0.366_fPrec
  hcut(iMat)     =   0.01_fPrec
  radl(iMat)     =   0.0143_fPrec
  bnref(iMat)    = 217.8_fPrec

  ! Tungsten
  iMat = iMat + 1
  colmats(iMat)  = "W"
  exenergy(iMat) = 727.0e-9_fPrec
  anuc(iMat)     = 183.85_fPrec
  zatom(iMat)    =  74.00_fPrec
  rho(iMat)      =  19.30_fPrec
  emr(iMat)      =   0.520_fPrec
  hcut(iMat)     =   0.01_fPrec
  radl(iMat)     =   0.0035_fPrec
  bnref(iMat)    = 440.3_fPrec

  ! Lead
  iMat = iMat + 1
  colmats(iMat)  = "PB"
  exenergy(iMat) = 823.0e-9_fPrec
  anuc(iMat)     = 207.19_fPrec
  zatom(iMat)    =  82.00_fPrec
  rho(iMat)      =  11.35_fPrec
  emr(iMat)      =   0.542_fPrec
  hcut(iMat)     =   0.01_fPrec
  radl(iMat)     =   0.0056_fPrec
  bnref(iMat)    = 455.3_fPrec

  ! Carbon
  iMat = iMat + 1
  colmats(iMat)  = "C"
  exenergy(iMat) =  78.0e-9_fPrec
  anuc(iMat)     =  12.01_fPrec
  zatom(iMat)    =   6.00_fPrec
  rho(iMat)      =   1.67_fPrec
  emr(iMat)      =   0.25_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.2557_fPrec
  bnref(iMat)    =  70.0_fPrec

  ! Carbon2
  iMat = iMat + 1
  colmats(iMat)  = "C2"
  exenergy(iMat) =  78.0e-9_fPrec
  anuc(iMat)     =  12.01_fPrec
  zatom(iMat)    =   6.00_fPrec
  rho(iMat)      =   4.52_fPrec
  emr(iMat)      =   0.25_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.094_fPrec
  bnref(iMat)    =  70.0_fPrec

  ! Silicon
  iMat = iMat + 1
  colmats(iMat)  = "Si"
  exenergy(iMat) = 173.0e-9_fPrec
  anuc(iMat)     =  28.08_fPrec
  zatom(iMat)    =  14.00_fPrec
  rho(iMat)      =   2.33_fPrec
  emr(iMat)      =   0.441_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     = one
  bnref(iMat)    = 120.14_fPrec

  ! Germanium
  iMat = iMat + 1
  colmats(iMat)  = "Ge"
  exenergy(iMat) = 350.0e-9_fPrec
  anuc(iMat)     =  72.63_fPrec
  zatom(iMat)    =  32.00_fPrec
  rho(iMat)      =   5.323_fPrec
  emr(iMat)      =   0.605_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     = one
  bnref(iMat)    = 226.35_fPrec

  iMat = iMat + 1
  colmats(iMat)  = "MoGR"
  exenergy(iMat) =  87.1e-9_fPrec
  anuc(iMat)     =  13.53_fPrec
  zatom(iMat)    =   6.65_fPrec
  rho(iMat)      =   2.500_fPrec
  emr(iMat)      =   0.25_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.1193_fPrec
  bnref(iMat)    =  76.7_fPrec

  iMat = iMat + 1
  colmats(iMat)  = "CuCD"
  exenergy(iMat) = 152.9e-9_fPrec
  anuc(iMat)     =  25.24_fPrec
  zatom(iMat)    =  11.90_fPrec
  rho(iMat)      =   5.40_fPrec
  emr(iMat)      =   0.308_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.0316_fPrec
  bnref(iMat)    = 115.0_fPrec

  iMat = iMat + 1
  colmats(iMat)  = "Mo"
  exenergy(iMat) = 424.0e-9_fPrec
  anuc(iMat)     =  95.96_fPrec
  zatom(iMat)    =  42.00_fPrec
  rho(iMat)      =  10.22_fPrec
  emr(iMat)      =   0.481_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.0096_fPrec
  bnref(iMat)    = 273.9_fPrec

  iMat = iMat + 1
  colmats(iMat)  = "Glid"
  exenergy(iMat) = 320.8e-9_fPrec
  anuc(iMat)     =  63.15_fPrec
  zatom(iMat)    =  28.80_fPrec
  rho(iMat)      =   8.93_fPrec
  emr(iMat)      =   0.418_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.0144_fPrec
  bnref(iMat)    = 208.7_fPrec

  iMat = iMat + 1
  colmats(iMat)  = "Iner"
  exenergy(iMat) = 682.2e-9_fPrec
  anuc(iMat)     = 166.70_fPrec
  zatom(iMat)    =  67.70_fPrec
  rho(iMat)      =  18.00_fPrec
  emr(iMat)      =   0.578_fPrec
  hcut(iMat)     =   0.02_fPrec
  radl(iMat)     =   0.00385_fPrec
  bnref(iMat)    = 392.1_fPrec

  if(iMat > nrmat) then
    write(lerr,"(a)") "COLL> ERROR Variable imat > nrmat in collmat_init. Please increase nrmat."
    !call prror
  end if

  ! The following two must always be the last two materials

  ! Vacuum
  iMat = iMat + 1
  colmats(iMat)  = "VA"
  exenergy(iMat) = zero
  anuc(iMat)     = zero
  zatom(iMat)    = zero
  rho(iMat)      = zero
  emr(iMat)      = zero
  hcut(iMat)     = zero
  radl(iMat)     = 1.0e12_fPrec
  bnref(iMat)    = zero

  ! Black Absorber
  iMat = iMat + 1
  colmats(iMat)  = "BL"
  exenergy(iMat) = c1e10
  anuc(iMat)     = zero
  zatom(iMat)    = zero
  rho(iMat)      = zero
  emr(iMat)      = zero
  hcut(iMat)     = zero
  radl(iMat)     = 1.0e12_fPrec
  bnref(iMat)    = zero

  if(iMat > nmat) then
    write(lerr,"(a)") "COLL> ERROR Variable imat > nmat in collmat_init. Please increase nmat."
    !call prror
  end if

end subroutine collmat_init

! ================================================================================================ !
!  V.K. Berglyd Olsen, BE-ABP-HSS
!  Created: 2019-09-16
!  Updated: 2019-09-16
!  Get collimator material number from name (case sensitive)
! ================================================================================================ !
integer function collmat_getCollMatID(matName)

  character(len=*), intent(in) :: matName
  integer i, matID

  matID = -1
  do i=1,nmat
    if(colmats(i) == matName) then
      matID = i
      exit
    end if
  end do
  collmat_getCollMatID = matID

end function collmat_getCollMatID

end module coll_materials

max_stars_repo_name: xsuite/xcol

max_issues_repo_name: xsuite/xcol

max_forks_repo_name: xsuite/xcol


      SUBROUTINE DGBSL(ABD, LDA, N, ML, MU, IPVT, B, JOB)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
!...Translated by Pacific-Sierra Research 77to90  4.3E  21:52:17  11/14/01
!...Switches:
!-----------------------------------------------
!   I n t e r f a c e   B l o c k s
!-----------------------------------------------
      USE ddot_I
      USE daxpy_I
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: LDA
      INTEGER , INTENT(IN) :: N
      INTEGER , INTENT(IN) :: ML
      INTEGER , INTENT(IN) :: MU
      INTEGER , INTENT(IN) :: JOB
      INTEGER , INTENT(IN) :: IPVT(1)
      REAL(DOUBLE)  :: ABD(LDA,1)
      REAL(DOUBLE)  :: B(1)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: K, KB, L, LA, LB, LM, M, NM1
      REAL(DOUBLE) :: T
!-----------------------------------------------
!
!     DGBSL SOLVES THE DOUBLE PRECISION BAND SYSTEM
!     A * X = B  OR  TRANS(A) * X = B
!     USING THE FACTORS COMPUTED BY DGBCO OR DGBFA.
!
!     ON ENTRY
!
!        ABD     DOUBLE PRECISION(LDA, N)
!                THE OUTPUT FROM DGBCO OR DGBFA.
!
!        LDA     INTEGER
!                THE LEADING DIMENSION OF THE ARRAY  ABD .
!
!        N       INTEGER
!                THE ORDER OF THE ORIGINAL MATRIX.
!
!        ML      INTEGER
!                NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
!
!        MU      INTEGER
!                NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
!
!        IPVT    INTEGER(N)
!                THE PIVOT VECTOR FROM DGBCO OR DGBFA.
!
!        B       DOUBLE PRECISION(N)
!                THE RIGHT HAND SIDE VECTOR.
!
!        JOB     INTEGER
!                = 0         TO SOLVE  A*X = B ,
!                = NONZERO   TO SOLVE  TRANS(A)*X = B , WHERE
!                            TRANS(A)  IS THE TRANSPOSE.
!
!     ON RETURN
!
!        B       THE SOLUTION VECTOR  X .
!
!     ERROR CONDITION
!
!        A DIVISION BY ZERO WILL OCCUR IF THE INPUT FACTOR CONTAINS A
!        ZERO ON THE DIAGONAL.  TECHNICALLY THIS INDICATES SINGULARITY
!        BUT IT IS OFTEN CAUSED BY IMPROPER ARGUMENTS OR IMPROPER
!        SETTING OF LDA .  IT WILL NOT OCCUR IF THE SUBROUTINES ARE
!        CALLED CORRECTLY AND IF DGBCO HAS SET RCOND .GT. 0.0
!        OR DGBFA HAS SET INFO .EQ. 0 .
!
!     TO COMPUTE  INVERSE(A) * C  WHERE  C  IS A MATRIX
!     WITH  P  COLUMNS
!           CALL DGBCO(ABD,LDA,N,ML,MU,IPVT,RCOND,Z)
!           IF (RCOND IS TOO SMALL) GO TO ...
!           DO 10 J = 1, P
!              CALL DGBSL(ABD,LDA,N,ML,MU,IPVT,C(1,J),0)
!        10 CONTINUE
!
!     LINPACK. THIS VERSION DATED 08/14/78 .
!     CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
!
!     SUBROUTINES AND FUNCTIONS
!
!     BLAS DAXPY,DDOT
!     FORTRAN MIN0
!
!     INTERNAL VARIABLES
!
!
      M = MU + ML + 1
      NM1 = N - 1
      IF (JOB == 0) THEN
!
!        JOB = 0 , SOLVE  A * X = B
!        FIRST SOLVE L*Y = B
!
         IF (ML /= 0) THEN
            IF (NM1 >= 1) THEN
               DO K = 1, NM1
                  LM = MIN0(ML,N - K)
                  L = IPVT(K)
                  T = B(L)
                  IF (L /= K) THEN
                     B(L) = B(K)
                     B(K) = T
                  ENDIF
                  CALL DAXPY (LM, T, ABD(M+1,K), 1, B(K+1), 1)
               END DO
            ENDIF
         ENDIF
         DO KB = 1, N
            K = N + 1 - KB
            B(K) = B(K)/ABD(M,K)
            LM = MIN0(K,M) - 1
            LA = M - LM
            LB = K - LM
            T = -B(K)
            CALL DAXPY (LM, T, ABD(LA,K), 1, B(LB), 1)
         END DO
      ELSE
         DO K = 1, N
            LM = MIN0(K,M) - 1
            LA = M - LM
            LB = K - LM
            T = DDOT(LM,ABD(LA,K),1,B(LB),1)
            B(K) = (B(K)-T)/ABD(M,K)
         END DO
!
!        NOW SOLVE TRANS(L)*X = Y
!
         IF (ML /= 0) THEN
            IF (NM1 >= 1) THEN
               DO KB = 1, NM1
                  K = N - KB
                  LM = MIN0(ML,N - K)
                  B(K) = B(K) + DDOT(LM,ABD(M+1,K),1,B(K+1),1)
                  L = IPVT(K)
                  IF (L == K) CYCLE
                  T = B(L)
                  B(L) = B(K)
                  B(K) = T
               END DO
            ENDIF
         ENDIF
      ENDIF
      RETURN
      END SUBROUTINE DGBSL


text:
! Simple tests symbol handling with use statement.
module test2008_57_A
   integer i_var

   interface operator (.add.)
      module procedure real_add
   end interface

   contains
      function real_add(a,b)
         integer :: real_add
         integer, intent(in) :: a,b
         real_add = 1
      end function real_add

end module

function foo()
   use test2008_57_A, only : operator(.add.)

   integer a,b,c

 ! Make a call to the function from module "A"
   a = b .add. c

end function

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!===============================================================================
module helper_module

  implicit none

  interface check_value
    module procedure check_values_int
    module procedure check_values_real
  end interface check_value

contains

  subroutine construct_fakeCopy( ptr, arr, lb, ub )
    ! input parameters
    integer, intent(in) :: lb(2), ub(2)
    real, dimension(:,:), intent(inout), pointer :: ptr
    real, dimension( lb(1):ub(1), lb(2):ub(2) ), intent(in), target  :: arr

    ptr( lb(1):ub(1), lb(2):ub(2) ) => arr

  end subroutine construct_fakeCopy

  function check_values_int( val1, val2 )
    integer :: val1, val2
    logical :: check_values_int

    if( val1 == val2 ) then
      check_values_int = .true.
    else
      check_values_int = .false.
    endif
  end function check_values_int

  function check_values_real( val1, val2 )
    real :: val1, val2
    logical :: check_values_real

    if( val1 == val2 ) then
      check_values_real = .true.
    else
      check_values_real = .false.
    endif
  end function check_values_real

end module helper_module
!===============================================================================
!===============================================================================
program dummy_array_pointer

  use helper_module

  implicit none

  ! boundaries for the real array
  integer, parameter, dimension(2) :: lb = (/ -1, 2 /)
  integer, parameter, dimension(2) :: ub = (/  4, 10 /)

  ! test value - to check correctness
  real, parameter :: test_val = 5

  ! real array - will be filled with values to check correctness of pointer
  real, dimension( lb(1):ub(1), lb(2):ub(2) ) :: test_arr

  ! pointer for connecting to array
  real, pointer :: ptr_arr(:,:)

  ! loop variables
  integer :: i, j

  ! set values in array
  do j = lb(2), ub(2)
    do i = lb(1), ub(1)
      test_arr(i,j) = test_val
    enddo
  enddo

  ! create an association of test_arr
  call construct_fakeCopy( ptr_arr, test_arr, lb, ub )

  print *, "checking lb(1) = ", check_value( lbound(ptr_arr, dim=1), lb(1) )
  print *, "checking lb(2) = ", check_value( lbound(ptr_arr, dim=2), lb(2) )

  print *, "checking ub(1) = ", check_value( ubound(ptr_arr, dim=1), ub(1) )
  print *, "checking ub(2) = ", check_value( ubound(ptr_arr, dim=2), ub(2) )

  do j = lb(2), ub(2)
    do i = lb(1), ub(1)
      if ( .not. check_value( ptr_arr(i, j), test_val ) ) then
        print *, "wrong value at ", i, j
        stop
      endif
    enddo
  enddo

  nullify(ptr_arr)

end program dummy_array_pointer

max_stars_repo_name: ahelm/fortran-snippets

max_issues_repo_name: ahelm/fortran-snippets

max_forks_repo_name: ahelm/fortran-snippets


text:
      subroutine ksopt (isel,x,obj,g,df,dg,nomax,ngmax,work)
      implicit double precision (a-h,o-z)
      common /kscomm/ rdf   ,adf   ,fdl   ,fdm   ,rho   ,drho  ,rhomax,
     1                fun0  ,slope ,delx  ,alpha ,alpmax,a1    ,a2    ,
     2                a3    ,a4    ,f1    ,f2    ,f3    ,f4    ,alim  ,
     3                atest ,ftest ,ifscl ,ifoff ,isx   ,isx0  ,isxlb ,
     4                isxub ,iscl  ,ig0   ,idf   ,islp  ,iobj0 ,iy    ,
     5                ip    ,ih    ,ihess ,iside ,isact ,idobj ,idg   ,
     6                itmp1 ,itmp2 ,inext ,jnext ,jsel  ,itcnt ,icntr ,
     7                icnta ,isdflg,isdrst,ifncl ,nunit ,ndv   ,ncon  ,
     8                nobj  ,nside ,nscale,iprnt ,itmax ,igrad ,limit
      dimension x(*),obj(*),g(*),df(nomax,*),dg(ngmax,*),work(*)
c
c          KSOPT is a general purpose multiobjective optimization
c          program which performs a constrained to unconstrained
c          transformation on the given optimization problem.
c          The problem is then solved using a Davidon-Fletcher-Powell
c          unconstrained optimization method.
c
c          input parameters are
c               obj   -- current set of objective function values
c               g     -- current set of constraints
c               df    -- derivatives of objective funcs. w.r.t. x-vector
c               dg    -- derivatives of constraints w.r.t. x-vector
c               nomax -- first dimension of df matrix
c               ngmax -- first dimension of dg matrix
c               work  -- working array from routine ksinit
c                        which must not be altered during optimization.
c
c          output parameters are
c               isel  -- flag requesting user supplied information
c                        =0 -- minimization terminated
c                        =1 -- user supplies obj and g
c                        =2 -- user supplies df and dg
c               x     -- new set of design variables
c
c          author   - Gregory A. Wrenn
c          location - Lockheed Engineering and Sciences Co.
c                     144 Research Drive
c                     Hampton, Va. 23666
c
c          last modification - 19 July 1996
c
c
c          choose a path based on the value of jsel
c
      call kscomg (work(1))
      iprnt3 = iprnt / 100
      ileft  = mod(iprnt,100)
      iprnt2 = ileft / 10
      iprnt1 = mod(ileft,10)
      nodim  = nobj
      if (nodim .le. 0) nodim = 1
      ncdim  = ncon
      if (ncdim .le. 0) ncdim = 1
c
      isel = 0
      if (jsel .le. 0) return
c
      go to (100,200,300,400,500,600) jsel
c
c          begin an iteration
c
  100 continue
      isel = 1
      jsel = 2
      inext = 0
      if (itcnt .eq. 1) go to 900
c
c          compute each objective function's scale and offset
c
  200 continue
      call ksando (obj,g,work(ifscl),work(ifoff),ncon,nobj)
c
c          calculate initial unconstrained function
c
      call ksfun (fun0,obj,g,rho,work(ifscl),work(ifoff),ncon,nobj,
     1            work(itmp1))
      if (igrad .eq. 1) go to 350
c
c          get finite difference gradients of obj and constraints
c
  300 continue
      jsel = 3
      isel = 1
      call ksgrad (inext,work(isx),work(isx0),work(isxlb),work(isxub),
     1             g,work(ig0),obj,work(iobj0),work(idobj),work(idg),
     2             work(iscl),delx,ndv,ncon,nobj,nside,fdl,fdm,
     3             nodim,ncdim)
      if (inext .ne. 0) go to 900
c
      do 310 i = 1,nobj
        obj(i) = work(iobj0 + i - 1)
  310 continue
c
      if (ncon .le. 0) go to 330
      do 320 i = 1,ncon
        g(i) = work(ig0 + i - 1)
  320 continue
  330 continue
      go to 440
c
c          get analytical gradients of obj and constraints
c
  350 continue
      do 360 i = 1,nobj
        work(iobj0 + i - 1) = obj(i)
  360 continue
c
      if (ncon .eq. 0) go to 380
      do 370 i = 1,ncon
        work(ig0 + i - 1) = g(i)
  370 continue
  380 continue
      isel = 2
      jsel = 4
      go to 900
c
c          scale analytical gradients of obj and constraints
c
  400 continue
      do 430 i = 1,ndv
        ss = work(iscl + i - 1)
c
        do 410 l = 1,nobj
          j = (i - 1) * nobj + l
          work(idobj + j - 1) = df(l,i) * ss
  410   continue
c
        if (ncon .eq. 0) go to 430
        do 420 k = 1,ncon
          j = (i - 1) * ncon + k
          work(idg + j - 1) = dg(k,i) * ss
  420   continue
  430 continue
c
c          compute gradient of unconstrained function
c          to be minimizied
c
  440 continue
      call ksdfun (work(idf),work(iobj0),work(ifscl),work(ifoff),
     1             work(idobj),work(ig0),work(idg),rho,ndv,ncon,nobj,
     2             work(itmp1),work(itmp2),nodim,ncdim)
c
c          compute an initial or restart hessian matrix
c
      call kshess (work(ihess),work(iobj0),work(ifscl),work(ifoff),
     1             work(idobj),work(ig0),work(idg),rho,
     2             ndv,ncon,nobj,work(itmp1),work(itmp2),nodim,ncdim)
c
c          check for side constraint violation
c
      call ksside (work(isx),work(isxlb),work(isxub),work(iside),
     1             work(idf),ndv,nside)
c
c          compute search direction
c
      call ksdfp (work(isx),work(iside),work(isact),work(idf),ndv,
     1            work(islp),slope,work(iy),work(ip),work(ih),
     2            work(ihess),isdflg)
c
c          print beginning of iteration information
c
      ipflag = 1
      call ksprnt (ipflag,iprnt1,iprnt2,x,work(iobj0),work(ig0),
     1             work(idobj),work(idg),work(iside),work(iscl),
     2             nodim,ncdim,work(itmp1),work(1))
c
c          compute unconstrained function to be minimized
c
  500 continue
      if (inext .eq. 0) go to 510
      call ksfun (fun,obj,g,rho,work(ifscl),work(ifoff),ncon,nobj,
     1            work(itmp1))
  510 continue
c
c          perform one-dimensional search
c
      jsel = 5
      isel = 1
      call ksoned (inext,jnext,work(isx),work(isx0),work(isxlb),
     1             work(isxub),fun,fun0,work(islp),slope,alpha,alpmax,
     2             ndv,a1,a2,a3,a4,f1,f2,f3,f4,alim,atest,ftest,
     3             nside,limit,nunit,iprnt3,work(iscl),work(itmp1),
     4             isdflg)
      if (inext .ne. 0) go to 900
      jsel = 6
      isel = 1
      go to 900
c
c          test for termination criteria
c
  600 continue
      call ksfun (fun,obj,g,rho,work(ifscl),work(ifoff),ncon,nobj,
     1            work(itmp1))
      if (itcnt .ge. itmax) go to 700
      af    = abs(fun)
      adelf = abs(fun - fun0)
      rdelf = adelf
      if (af .gt. 1.0e-12) rdelf = adelf / af
      icntr = icntr + 1
      icnta = icnta + 1
      if (rdelf .ge. rdf) icntr = 0
      if (adelf .ge. adf) icnta = 0
      if (icntr .ge. 3 .or. icnta .ge. 3) go to 700
c
c          go to next iteration
c
      itcnt = itcnt + 1
c
c          increment rho when necessary
c
      if (icntr .eq. 0 .and. icnta .eq. 0) go to 650
      rho = rho + drho
      if (rho .gt. rhomax) rho = rhomax
  650 continue
c
c          set flag for hessian matrix restart
c
      isdflg = isdflg + 1
      if (isdflg .gt. isdrst) isdflg = 0
c
c          re-scale when necessary
c
      iscale = 0
      if (nscale .gt. 0) iscale = nscale
      if (iscale .ne. 0) itflag = mod(itcnt - 1,iscale)
      if (iscale .eq. 0) itflag = 1
c
c          currently re-scale only when hessian matrix is reset
c          28 March 1991  g.a.w.
c
      itflag = 1
      if (isdflg .eq. 0) itflag = 0
      if (itflag .eq. 0) call ksscal (work(isx),work(isx0),work(isxlb),
     1                                work(isxub),work(iscl),
     2                                ndv,nside,nscale)
      go to 100
c
c          terminate
c
  700 continue
      isel = 0
c
      ipflag = 2
      call ksprnt (ipflag,iprnt1,iprnt2,x,obj,g,work(idobj),work(idg),
     1             work(iside),work(iscl),nodim,ncdim,work(itmp1),
     2             work(1))
c
c          un-scale design variables and return to calling routine
c
  900 continue
      call ksxlim (work(isx),work(isxlb),work(isxub),ndv,nside)
      call ksunsc (x,work(isx),work(iscl),ndv)
      if (isel .eq. 1) ifncl = ifncl + 1
      call kscomp (work(1))
      return
c
      end

ext: f

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
    !
    ! This program may be freely redistributed under the
    ! condition that the copyright notices (including this
    ! entire header) are not removed, and no compensation
    ! is received through use of the software.  Private,
    ! research, and institutional use is free.  You may
    ! distribute modified versions of this code UNDER THE
    ! CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE
    ! TO IT IN THE SAME FILE REMAIN UNDER COPYRIGHT OF THE
    ! ORIGINAL AUTHOR, BOTH SOURCE AND OBJECT CODE ARE
    ! MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR
    ! NOTICE IS GIVEN OF THE MODIFICATIONS.  Distribution
    ! of this code as part of a commercial system is
    ! permissible ONLY BY DIRECT ARRANGEMENT WITH THE
    ! AUTHOR.  (If you are not directly supplying this
    ! code to a customer, and you are instead telling them
    ! how they can obtain it for free, then you are not
    ! required to make any arrangement with me.)
    !
    ! Disclaimer:  Neither I nor: Columbia University, the
    ! National Aeronautics and Space Administration, nor
    ! the Massachusetts Institute of Technology warrant
    ! or certify this code in any way whatsoever.  This
    ! code is provided "as-is" to be used at your own risk.
    !
    !

    !
    ! UTIL1D.f90: util. func. for 1-dim. grid manipulation.
    !
    ! Darren Engwirda
    ! 31-Mar-2019
    ! de2363 [at] columbia [dot] edu
    !
    !

    subroutine linspace(xxll,xxuu,npos,xpos)

    !
    ! XXLL  lower-bound grid position.
    ! NNEW  upper-bound grid position.
    ! NPOS  no. edges in the grid.
    ! XPOS  array of grid edges. XPOS has length NPOS .
    !

        implicit none

        real*8 , intent(in)  :: xxll,xxuu
        integer, intent(in)  :: npos
        real*8 , intent(out) :: xpos(:)

        integer   :: ipos
        real*8    :: xdel

        xpos(   1) = xxll
        xpos(npos) = xxuu

        xdel = (xxuu-xxll) / (npos - 1)

        do  ipos = +2, npos-1

            xpos(ipos) = (ipos-1) * xdel

        end do

        return

    end  subroutine

    subroutine rndspace(xxll,xxuu,npos,xpos, &
        &               frac)

    !
    ! XXLL  lower-bound grid position.
    ! NNEW  upper-bound grid position.
    ! NPOS  no. edges in the grid.
    ! XPOS  array of grid edges. XPOS has length NPOS .
    ! FRAC  fractional perturbation of cell, OPTIONAL .
    !

        implicit none

        real*8 , intent(in)  :: xxll,xxuu
        integer, intent(in)  :: npos
        real*8 , intent(out) :: xpos(:)
        real*8 , intent(in), optional :: frac

        integer   :: ipos
        real*8    :: xdel,rand,move

        if (present(frac)) then
            move = +frac
        else
            move = 0.33d0
        end if

        xpos(   1) = xxll
        xpos(npos) = xxuu

        xdel = (xxuu-xxll) / (npos - 1)

        do  ipos = +2, npos-1

            xpos(ipos) = (ipos-1) * xdel

        end do

        do ipos = +2, npos-1

            call random_number (rand)

            rand = 2.d0 * (rand-.5d0)

            move = rand *  move

            xpos(ipos) = &
        &       xpos(ipos) + move * xdel

        end do

        return

    end  subroutine

max_stars_repo_name: mark-petersen/PPR

max_issues_repo_name: mark-petersen/PPR

max_forks_repo_name: mark-petersen/PPR


!*==dgges3.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief  DGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download DGGES3 + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE DGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
!                          SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR,
!                          LDVSR, WORK, LWORK, BWORK, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          JOBVSL, JOBVSR, SORT
!       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
!       ..
!       .. Array Arguments ..
!       LOGICAL            BWORK( * )
!       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
!      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
!      $                   VSR( LDVSR, * ), WORK( * )
!       ..
!       .. Function Arguments ..
!       LOGICAL            SELCTG
!       EXTERNAL           SELCTG
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> DGGES3 computes for a pair of N-by-N real nonsymmetric matrices (A,B),
!> the generalized eigenvalues, the generalized real Schur form (S,T),
!> optionally, the left and/or right matrices of Schur vectors (VSL and
!> VSR). This gives the generalized Schur factorization
!>
!>          (A,B) = ( (VSL)*S*(VSR)**T, (VSL)*T*(VSR)**T )
!>
!> Optionally, it also orders the eigenvalues so that a selected cluster
!> of eigenvalues appears in the leading diagonal blocks of the upper
!> quasi-triangular matrix S and the upper triangular matrix T.The
!> leading columns of VSL and VSR then form an orthonormal basis for the
!> corresponding left and right eigenspaces (deflating subspaces).
!>
!> (If only the generalized eigenvalues are needed, use the driver
!> DGGEV instead, which is faster.)
!>
!> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
!> or a ratio alpha/beta = w, such that  A - w*B is singular.  It is
!> usually represented as the pair (alpha,beta), as there is a
!> reasonable interpretation for beta=0 or both being zero.
!>
!> A pair of matrices (S,T) is in generalized real Schur form if T is
!> upper triangular with non-negative diagonal and S is block upper
!> triangular with 1-by-1 and 2-by-2 blocks.  1-by-1 blocks correspond
!> to real generalized eigenvalues, while 2-by-2 blocks of S will be
!> "standardized" by making the corresponding elements of T have the
!> form:
!>         [  a  0  ]
!>         [  0  b  ]
!>
!> and the pair of corresponding 2-by-2 blocks in S and T will have a
!> complex conjugate pair of generalized eigenvalues.
!>
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] JOBVSL
!> \verbatim
!>          JOBVSL is CHARACTER*1
!>          = 'N':  do not compute the left Schur vectors;
!>          = 'V':  compute the left Schur vectors.
!> \endverbatim
!>
!> \param[in] JOBVSR
!> \verbatim
!>          JOBVSR is CHARACTER*1
!>          = 'N':  do not compute the right Schur vectors;
!>          = 'V':  compute the right Schur vectors.
!> \endverbatim
!>
!> \param[in] SORT
!> \verbatim
!>          SORT is CHARACTER*1
!>          Specifies whether or not to order the eigenvalues on the
!>          diagonal of the generalized Schur form.
!>          = 'N':  Eigenvalues are not ordered;
!>          = 'S':  Eigenvalues are ordered (see SELCTG);
!> \endverbatim
!>
!> \param[in] SELCTG
!> \verbatim
!>          SELCTG is a LOGICAL FUNCTION of three DOUBLE PRECISION arguments
!>          SELCTG must be declared EXTERNAL in the calling subroutine.
!>          If SORT = 'N', SELCTG is not referenced.
!>          If SORT = 'S', SELCTG is used to select eigenvalues to sort
!>          to the top left of the Schur form.
!>          An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
!>          SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
!>          one of a complex conjugate pair of eigenvalues is selected,
!>          then both complex eigenvalues are selected.
!>
!>          Note that in the ill-conditioned case, a selected complex
!>          eigenvalue may no longer satisfy SELCTG(ALPHAR(j),ALPHAI(j),
!>          BETA(j)) = .TRUE. after ordering. INFO is to be set to N+2
!>          in this case.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
!> \endverbatim
!>
!> \param[in,out] A
!> \verbatim
!>          A is DOUBLE PRECISION array, dimension (LDA, N)
!>          On entry, the first of the pair of matrices.
!>          On exit, A has been overwritten by its generalized Schur
!>          form S.
!> \endverbatim
!>
!> \param[in] LDA
!> \verbatim
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!> \endverbatim
!>
!> \param[in,out] B
!> \verbatim
!>          B is DOUBLE PRECISION array, dimension (LDB, N)
!>          On entry, the second of the pair of matrices.
!>          On exit, B has been overwritten by its generalized Schur
!>          form T.
!> \endverbatim
!>
!> \param[in] LDB
!> \verbatim
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!> \endverbatim
!>
!> \param[out] SDIM
!> \verbatim
!>          SDIM is INTEGER
!>          If SORT = 'N', SDIM = 0.
!>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
!>          for which SELCTG is true.  (Complex conjugate pairs for which
!>          SELCTG is true for either eigenvalue count as 2.)
!> \endverbatim
!>
!> \param[out] ALPHAR
!> \verbatim
!>          ALPHAR is DOUBLE PRECISION array, dimension (N)
!> \endverbatim
!>
!> \param[out] ALPHAI
!> \verbatim
!>          ALPHAI is DOUBLE PRECISION array, dimension (N)
!> \endverbatim
!>
!> \param[out] BETA
!> \verbatim
!>          BETA is DOUBLE PRECISION array, dimension (N)
!>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
!>          be the generalized eigenvalues.  ALPHAR(j) + ALPHAI(j)*i,
!>          and  BETA(j),j=1,...,N are the diagonals of the complex Schur
!>          form (S,T) that would result if the 2-by-2 diagonal blocks of
!>          the real Schur form of (A,B) were further reduced to
!>          triangular form using 2-by-2 complex unitary transformations.
!>          If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
!>          positive, then the j-th and (j+1)-st eigenvalues are a
!>          complex conjugate pair, with ALPHAI(j+1) negative.
!>
!>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
!>          may easily over- or underflow, and BETA(j) may even be zero.
!>          Thus, the user should avoid naively computing the ratio.
!>          However, ALPHAR and ALPHAI will be always less than and
!>          usually comparable with norm(A) in magnitude, and BETA always
!>          less than and usually comparable with norm(B).
!> \endverbatim
!>
!> \param[out] VSL
!> \verbatim
!>          VSL is DOUBLE PRECISION array, dimension (LDVSL,N)
!>          If JOBVSL = 'V', VSL will contain the left Schur vectors.
!>          Not referenced if JOBVSL = 'N'.
!> \endverbatim
!>
!> \param[in] LDVSL
!> \verbatim
!>          LDVSL is INTEGER
!>          The leading dimension of the matrix VSL. LDVSL >=1, and
!>          if JOBVSL = 'V', LDVSL >= N.
!> \endverbatim
!>
!> \param[out] VSR
!> \verbatim
!>          VSR is DOUBLE PRECISION array, dimension (LDVSR,N)
!>          If JOBVSR = 'V', VSR will contain the right Schur vectors.
!>          Not referenced if JOBVSR = 'N'.
!> \endverbatim
!>
!> \param[in] LDVSR
!> \verbatim
!>          LDVSR is INTEGER
!>          The leading dimension of the matrix VSR. LDVSR >= 1, and
!>          if JOBVSR = 'V', LDVSR >= N.
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!> \endverbatim
!>
!> \param[in] LWORK
!> \verbatim
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!> \endverbatim
!>
!> \param[out] BWORK
!> \verbatim
!>          BWORK is LOGICAL array, dimension (N)
!>          Not referenced if SORT = 'N'.
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          = 1,...,N:
!>                The QZ iteration failed.  (A,B) are not in Schur
!>                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
!>                be correct for j=INFO+1,...,N.
!>          > N:  =N+1: other than QZ iteration failed in DHGEQZ.
!>                =N+2: after reordering, roundoff changed values of
!>                      some complex eigenvalues so that leading
!>                      eigenvalues in the Generalized Schur form no
!>                      longer satisfy SELCTG=.TRUE.  This could also
!>                      be caused due to scaling.
!>                =N+3: reordering failed in DTGSEN.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date January 2015
!
!> \ingroup doubleGEeigen
!
!  =====================================================================
      SUBROUTINE DGGES3(Jobvsl,Jobvsr,Sort,SELCTG,N,A,Lda,B,Ldb,Sdim,   &
     &                  Alphar,Alphai,Beta,Vsl,Ldvsl,Vsr,Ldvsr,Work,    &
     &                  Lwork,Bwork,Info)
      IMPLICIT NONE
!*--DGGES3286
!
!  -- LAPACK driver routine (version 3.6.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     January 2015
!
!     .. Scalar Arguments ..
      CHARACTER Jobvsl , Jobvsr , Sort
      INTEGER Info , Lda , Ldb , Ldvsl , Ldvsr , Lwork , N , Sdim
!     ..
!     .. Array Arguments ..
      LOGICAL Bwork(*)
      DOUBLE PRECISION A(Lda,*) , Alphai(*) , Alphar(*) , B(Ldb,*) ,    &
     &                 Beta(*) , Vsl(Ldvsl,*) , Vsr(Ldvsr,*) , Work(*)
!     ..
!     .. Function Arguments ..
      LOGICAL SELCTG
      EXTERNAL SELCTG
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ZERO , ONE
      PARAMETER (ZERO=0.0D+0,ONE=1.0D+0)
!     ..
!     .. Local Scalars ..
      LOGICAL cursl , ilascl , ilbscl , ilvsl , ilvsr , lastsl ,        &
     &        lquery , lst2sl , wantst
      INTEGER i , icols , ierr , ihi , ijobvl , ijobvr , ileft , ilo ,  &
     &        ip , iright , irows , itau , iwrk , lwkopt
      DOUBLE PRECISION anrm , anrmto , bignum , bnrm , bnrmto , eps ,   &
     &                 pvsl , pvsr , safmax , safmin , smlnum
!     ..
!     .. Local Arrays ..
      INTEGER idum(1)
      DOUBLE PRECISION dif(2)
!     ..
!     .. External Subroutines ..
      EXTERNAL DGEQRF , DGGBAK , DGGBAL , DGGHD3 , DHGEQZ , DLABAD ,    &
     &         DLACPY , DLASCL , DLASET , DORGQR , DORMQR , DTGSEN ,    &
     &         XERBLA
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      DOUBLE PRECISION DLAMCH , DLANGE
      EXTERNAL LSAME , DLAMCH , DLANGE
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , MAX , SQRT
!     ..
!     .. Executable Statements ..
!
!     Decode the input arguments
!
      IF ( LSAME(Jobvsl,'N') ) THEN
         ijobvl = 1
         ilvsl = .FALSE.
      ELSEIF ( LSAME(Jobvsl,'V') ) THEN
         ijobvl = 2
         ilvsl = .TRUE.
      ELSE
         ijobvl = -1
         ilvsl = .FALSE.
      ENDIF
!
      IF ( LSAME(Jobvsr,'N') ) THEN
         ijobvr = 1
         ilvsr = .FALSE.
      ELSEIF ( LSAME(Jobvsr,'V') ) THEN
         ijobvr = 2
         ilvsr = .TRUE.
      ELSE
         ijobvr = -1
         ilvsr = .FALSE.
      ENDIF
!
      wantst = LSAME(Sort,'S')
!
!     Test the input arguments
!
      Info = 0
      lquery = (Lwork==-1)
      IF ( ijobvl<=0 ) THEN
         Info = -1
      ELSEIF ( ijobvr<=0 ) THEN
         Info = -2
      ELSEIF ( (.NOT.wantst) .AND. (.NOT.LSAME(Sort,'N')) ) THEN
         Info = -3
      ELSEIF ( N<0 ) THEN
         Info = -5
      ELSEIF ( Lda<MAX(1,N) ) THEN
         Info = -7
      ELSEIF ( Ldb<MAX(1,N) ) THEN
         Info = -9
      ELSEIF ( Ldvsl<1 .OR. (ilvsl .AND. Ldvsl<N) ) THEN
         Info = -15
      ELSEIF ( Ldvsr<1 .OR. (ilvsr .AND. Ldvsr<N) ) THEN
         Info = -17
      ELSEIF ( Lwork<6*N+16 .AND. .NOT.lquery ) THEN
         Info = -19
      ENDIF
!
!     Compute workspace
!
      IF ( Info==0 ) THEN
         CALL DGEQRF(N,N,B,Ldb,Work,Work,-1,ierr)
         lwkopt = MAX(6*N+16,3*N+INT(Work(1)))
         CALL DORMQR('L','T',N,N,N,B,Ldb,Work,A,Lda,Work,-1,ierr)
         lwkopt = MAX(lwkopt,3*N+INT(Work(1)))
         IF ( ilvsl ) THEN
            CALL DORGQR(N,N,N,Vsl,Ldvsl,Work,Work,-1,ierr)
            lwkopt = MAX(lwkopt,3*N+INT(Work(1)))
         ENDIF
         CALL DGGHD3(Jobvsl,Jobvsr,N,1,N,A,Lda,B,Ldb,Vsl,Ldvsl,Vsr,     &
     &               Ldvsr,Work,-1,ierr)
         lwkopt = MAX(lwkopt,3*N+INT(Work(1)))
         CALL DHGEQZ('S',Jobvsl,Jobvsr,N,1,N,A,Lda,B,Ldb,Alphar,Alphai, &
     &               Beta,Vsl,Ldvsl,Vsr,Ldvsr,Work,-1,ierr)
         lwkopt = MAX(lwkopt,2*N+INT(Work(1)))
         IF ( wantst ) THEN
            CALL DTGSEN(0,ilvsl,ilvsr,Bwork,N,A,Lda,B,Ldb,Alphar,Alphai,&
     &                  Beta,Vsl,Ldvsl,Vsr,Ldvsr,Sdim,pvsl,pvsr,dif,    &
     &                  Work,-1,idum,1,ierr)
            lwkopt = MAX(lwkopt,2*N+INT(Work(1)))
         ENDIF
         Work(1) = lwkopt
      ENDIF
!
      IF ( Info/=0 ) THEN
         CALL XERBLA('DGGES3 ',-Info)
         RETURN
      ELSEIF ( lquery ) THEN
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( N==0 ) THEN
         Sdim = 0
         RETURN
      ENDIF
!
!     Get machine constants
!
      eps = DLAMCH('P')
      safmin = DLAMCH('S')
      safmax = ONE/safmin
      CALL DLABAD(safmin,safmax)
      smlnum = SQRT(safmin)/eps
      bignum = ONE/smlnum
!
!     Scale A if max element outside range [SMLNUM,BIGNUM]
!
      anrm = DLANGE('M',N,N,A,Lda,Work)
      ilascl = .FALSE.
      IF ( anrm>ZERO .AND. anrm<smlnum ) THEN
         anrmto = smlnum
         ilascl = .TRUE.
      ELSEIF ( anrm>bignum ) THEN
         anrmto = bignum
         ilascl = .TRUE.
      ENDIF
      IF ( ilascl ) CALL DLASCL('G',0,0,anrm,anrmto,N,N,A,Lda,ierr)
!
!     Scale B if max element outside range [SMLNUM,BIGNUM]
!
      bnrm = DLANGE('M',N,N,B,Ldb,Work)
      ilbscl = .FALSE.
      IF ( bnrm>ZERO .AND. bnrm<smlnum ) THEN
         bnrmto = smlnum
         ilbscl = .TRUE.
      ELSEIF ( bnrm>bignum ) THEN
         bnrmto = bignum
         ilbscl = .TRUE.
      ENDIF
      IF ( ilbscl ) CALL DLASCL('G',0,0,bnrm,bnrmto,N,N,B,Ldb,ierr)
!
!     Permute the matrix to make it more nearly triangular
!
      ileft = 1
      iright = N + 1
      iwrk = iright + N
      CALL DGGBAL('P',N,A,Lda,B,Ldb,ilo,ihi,Work(ileft),Work(iright),   &
     &            Work(iwrk),ierr)
!
!     Reduce B to triangular form (QR decomposition of B)
!
      irows = ihi + 1 - ilo
      icols = N + 1 - ilo
      itau = iwrk
      iwrk = itau + irows
      CALL DGEQRF(irows,icols,B(ilo,ilo),Ldb,Work(itau),Work(iwrk),     &
     &            Lwork+1-iwrk,ierr)
!
!     Apply the orthogonal transformation to matrix A
!
      CALL DORMQR('L','T',irows,icols,irows,B(ilo,ilo),Ldb,Work(itau),  &
     &            A(ilo,ilo),Lda,Work(iwrk),Lwork+1-iwrk,ierr)
!
!     Initialize VSL
!
      IF ( ilvsl ) THEN
         CALL DLASET('Full',N,N,ZERO,ONE,Vsl,Ldvsl)
         IF ( irows>1 ) CALL DLACPY('L',irows-1,irows-1,B(ilo+1,ilo),   &
     &                              Ldb,Vsl(ilo+1,ilo),Ldvsl)
         CALL DORGQR(irows,irows,irows,Vsl(ilo,ilo),Ldvsl,Work(itau),   &
     &               Work(iwrk),Lwork+1-iwrk,ierr)
      ENDIF
!
!     Initialize VSR
!
      IF ( ilvsr ) CALL DLASET('Full',N,N,ZERO,ONE,Vsr,Ldvsr)
!
!     Reduce to generalized Hessenberg form
!
      CALL DGGHD3(Jobvsl,Jobvsr,N,ilo,ihi,A,Lda,B,Ldb,Vsl,Ldvsl,Vsr,    &
     &            Ldvsr,Work(iwrk),Lwork+1-iwrk,ierr)
!
!     Perform QZ algorithm, computing Schur vectors if desired
!
      iwrk = itau
      CALL DHGEQZ('S',Jobvsl,Jobvsr,N,ilo,ihi,A,Lda,B,Ldb,Alphar,Alphai,&
     &            Beta,Vsl,Ldvsl,Vsr,Ldvsr,Work(iwrk),Lwork+1-iwrk,ierr)
      IF ( ierr/=0 ) THEN
         IF ( ierr>0 .AND. ierr<=N ) THEN
            Info = ierr
         ELSEIF ( ierr>N .AND. ierr<=2*N ) THEN
            Info = ierr - N
         ELSE
            Info = N + 1
         ENDIF
         GOTO 100
      ENDIF
!
!     Sort eigenvalues ALPHA/BETA if desired
!
      Sdim = 0
      IF ( wantst ) THEN
!
!        Undo scaling on eigenvalues before SELCTGing
!
         IF ( ilascl ) THEN
            CALL DLASCL('G',0,0,anrmto,anrm,N,1,Alphar,N,ierr)
            CALL DLASCL('G',0,0,anrmto,anrm,N,1,Alphai,N,ierr)
         ENDIF
         IF ( ilbscl ) CALL DLASCL('G',0,0,bnrmto,bnrm,N,1,Beta,N,ierr)
!
!        Select eigenvalues
!
         DO i = 1 , N
            Bwork(i) = SELCTG(Alphar(i),Alphai(i),Beta(i))
         ENDDO
!
         CALL DTGSEN(0,ilvsl,ilvsr,Bwork,N,A,Lda,B,Ldb,Alphar,Alphai,   &
     &               Beta,Vsl,Ldvsl,Vsr,Ldvsr,Sdim,pvsl,pvsr,dif,       &
     &               Work(iwrk),Lwork-iwrk+1,idum,1,ierr)
         IF ( ierr==1 ) Info = N + 3
!
      ENDIF
!
!     Apply back-permutation to VSL and VSR
!
      IF ( ilvsl ) CALL DGGBAK('P','L',N,ilo,ihi,Work(ileft),           &
     &                         Work(iright),N,Vsl,Ldvsl,ierr)
!
      IF ( ilvsr ) CALL DGGBAK('P','R',N,ilo,ihi,Work(ileft),           &
     &                         Work(iright),N,Vsr,Ldvsr,ierr)
!
!     Check if unscaling would cause over/underflow, if so, rescale
!     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
!     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
!
      IF ( ilascl ) THEN
         DO i = 1 , N
            IF ( Alphai(i)/=ZERO ) THEN
               IF ( (Alphar(i)/safmax)>(anrmto/anrm) .OR.               &
     &              (safmin/Alphar(i))>(anrm/anrmto) ) THEN
                  Work(1) = ABS(A(i,i)/Alphar(i))
                  Beta(i) = Beta(i)*Work(1)
                  Alphar(i) = Alphar(i)*Work(1)
                  Alphai(i) = Alphai(i)*Work(1)
               ELSEIF ( (Alphai(i)/safmax)>(anrmto/anrm) .OR.           &
     &                  (safmin/Alphai(i))>(anrm/anrmto) ) THEN
                  Work(1) = ABS(A(i,i+1)/Alphai(i))
                  Beta(i) = Beta(i)*Work(1)
                  Alphar(i) = Alphar(i)*Work(1)
                  Alphai(i) = Alphai(i)*Work(1)
               ENDIF
            ENDIF
         ENDDO
      ENDIF
!
      IF ( ilbscl ) THEN
         DO i = 1 , N
            IF ( Alphai(i)/=ZERO ) THEN
               IF ( (Beta(i)/safmax)>(bnrmto/bnrm) .OR. (safmin/Beta(i))&
     &              >(bnrm/bnrmto) ) THEN
                  Work(1) = ABS(B(i,i)/Beta(i))
                  Beta(i) = Beta(i)*Work(1)
                  Alphar(i) = Alphar(i)*Work(1)
                  Alphai(i) = Alphai(i)*Work(1)
               ENDIF
            ENDIF
         ENDDO
      ENDIF
!
!     Undo scaling
!
      IF ( ilascl ) THEN
         CALL DLASCL('H',0,0,anrmto,anrm,N,N,A,Lda,ierr)
         CALL DLASCL('G',0,0,anrmto,anrm,N,1,Alphar,N,ierr)
         CALL DLASCL('G',0,0,anrmto,anrm,N,1,Alphai,N,ierr)
      ENDIF
!
      IF ( ilbscl ) THEN
         CALL DLASCL('U',0,0,bnrmto,bnrm,N,N,B,Ldb,ierr)
         CALL DLASCL('G',0,0,bnrmto,bnrm,N,1,Beta,N,ierr)
      ENDIF
!
      IF ( wantst ) THEN
!
!        Check if reordering is correct
!
         lastsl = .TRUE.
         lst2sl = .TRUE.
         Sdim = 0
         ip = 0
         DO i = 1 , N
            cursl = SELCTG(Alphar(i),Alphai(i),Beta(i))
            IF ( Alphai(i)==ZERO ) THEN
               IF ( cursl ) Sdim = Sdim + 1
               ip = 0
               IF ( cursl .AND. .NOT.lastsl ) Info = N + 2
            ELSEIF ( ip==1 ) THEN
!
!                 Last eigenvalue of conjugate pair
!
               cursl = cursl .OR. lastsl
               lastsl = cursl
               IF ( cursl ) Sdim = Sdim + 2
               ip = -1
               IF ( cursl .AND. .NOT.lst2sl ) Info = N + 2
            ELSE
!
!                 First eigenvalue of conjugate pair
!
               ip = 1
            ENDIF
            lst2sl = lastsl
            lastsl = cursl
         ENDDO
!
      ENDIF
!
!
 100  Work(1) = lwkopt
!
!
!     End of DGGES3
!
      END SUBROUTINE DGGES3


! PR fortran/80918
! { dg-do compile }

subroutine foo (a)
  integer :: a(*)
  !$omp task depend(inout:a)
  !$omp end task
  !$omp task depend(inout:a)
  !$omp end task
end subroutine foo


text:
C  *********************************************************
C  *                                                       *
C  *    SUBROUTINE 04.02.04.01/locat2                      *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      SUBROUTINE LOCAT2 (YLOC, XVAL)

C  LOCAT2 draws an annotation text string at the real y-location
C  corresponding to the integer index given in YLOC, and puts a
C  polymarker at the text position.

      INTEGER    YLOC

      REAL       XVAL, YLOCEL, XA(1),YA(1)

      XA(1) = XVAL
      YA(1) = YLOCEL(YLOC)
      CALL PPM (1, XA,YA)
      CALL LOCATX (YLOC, XVAL)

      END

C  *********************************************************
C  *                                                       *
C  *    SUBROUTINE 04.02.04.01/locans                      *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      SUBROUTINE LOCANS (YLOC, XVAL)

C  LOCANS draws an annotation text primitive with a distinct
C  annotation point and reference point, at the real y-location
C  corresponding to the integer index given in YLOC.

      INTEGER   YLOC
      REAL      XVAL, YLOCEL, XOFF,YOFF, NPCX,NPCY,NPCPWC

      DATA XOFF /-0.15/, YOFF /0.05/

      CALL WCNPC (0.0, 0.0, NPCX,NPCY, NPCPWC)
      CALL PATR (XVAL-XOFF, YLOCEL(YLOC)-YOFF,
     1           XOFF*NPCPWC,YOFF*NPCPWC, 'Hiho')

      END

C  *********************************************************
C  *                                                       *
C  *    SUBROUTINE 04.02.04.01/tsanup                      *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      SUBROUTINE TSANUP (TXFONT, TXPATH, TXSTR, NUMOUT, NGSTR)

C  TSANUP draws an annotation text string with various values for
C  character up and for annotation text alignment.  The annotation
C  text strings are displayed along radii of a circle and are
C  enclosed by a text extent rectangle (except for one which is
C  deliberately drawn incorrectly).  TSANUP assumes that default
C  character height has been set to 1.0.

C  Input parameters:
C    TXFONT : Text font to be used
C    TXPATH : Annotation text path to be used
C    TXSTR  : Annotation text string to be used
C  Output parameters:
C    NUMOUT : Number of annotation text primitives drawn
C    NGSTR  : Identifier of incorrect annotation text primitive

C composition type
C                preconcatenate  postconcatenate  replace
      INTEGER    PCPRE,          PCPOST,          PCREPL
      PARAMETER (PCPRE = 0,      PCPOST = 1,      PCREPL = 2)

C text alignment horizontal
      INTEGER    PAHNOR,     PALEFT,     PACENT,     PARITE
      PARAMETER (PAHNOR = 0, PALEFT = 1, PACENT = 2, PARITE = 3)
C text alignment vertical
      INTEGER    PAVNOR,     PATOP,      PACAP,     PAHALF,
     1           PABASE,     PABOTT
      PARAMETER (PAVNOR = 0, PATOP  = 1, PACAP = 2, PAHALF = 3,
     1           PABASE = 4, PABOTT = 5)
C text path
      INTEGER    PRIGHT,     PLEFT,     PUP,     PDOWN
      PARAMETER (PRIGHT = 0, PLEFT = 1, PUP = 2, PDOWN = 3)
C text precision
      INTEGER    PSTRP,     PCHARP,     PSTRKP
      PARAMETER (PSTRP = 0, PCHARP = 1, PSTRKP = 2)

C linetype
      INTEGER    PLSOLI,     PLDASH,     PLDOT,     PLDASD
      PARAMETER (PLSOLI = 1, PLDASH = 2, PLDOT = 3, PLDASD = 4)

      INTEGER    IX, RNDINT, NUMSTR, NUMOUT, NGSTR
C  numstr = number of strings to display
      PARAMETER (NUMSTR = 7)
      INTEGER    TXFONT, TXPATH, ERRIND, ALGNH(NUMSTR), ALGNV(NUMSTR)

      REAL       NTXWD,NTXHT, ORAD,IRAD, PI, XF(3,3), XA(9),YA(9)
      REAL       VTAB(PATOP:PABOTT), HTAB(PALEFT:PARITE), Z
      REAL       TXANG, SINT,COST, TXPOSX,TXPOSY, SCF, RECLNG,RECSHT
      REAL       XRT,YRT, XUP,YUP, XVEC,YVEC, INTSP, XB(2),YB(2), LABRAD
      REAL       NPCX,NPCY,NPCPWC

      PARAMETER (PI = 3.14159265)

      LOGICAL    VERT

      CHARACTER  TXSTR*(*), DIG2*2

      NUMOUT = NUMSTR
      VERT = TXPATH.EQ.PUP .OR. TXPATH.EQ.PDOWN

      CALL PSTXFN (TXFONT)
C  use <inquire text extent> on txstr to determine:
C    ntxwd  = nominal width of text extent rectangle
C    ntxht  = nominal height of text extent rectangle
C    vtab,htab = normalized heights and widths
      CALL TXEXAL (TXFONT, TXPATH, TXSTR, NTXWD, NTXHT, HTAB, VTAB)
      IF (VERT) THEN
         RECLNG = NTXHT
         RECSHT = NTXWD
      ELSE
         RECLNG = NTXWD
         RECSHT = NTXHT
      ENDIF

C  ngstr = incorrect string identifier
      NGSTR = RNDINT(1,7)
C  irad = inner radius of circle for display of annotation text strings
      IRAD = 1.1 * (NUMSTR*RECSHT) / (2*PI)
C  orad = outer radius of circle for display of annotation text strings
      ORAD = IRAD + RECLNG

C  set up transform to fit annotation text circle on screen, from MC
C  centered at origin, with radius orad to WC radius=0.4, centered
C  at 0.5,0.5
      Z = 0.0
C  label radius factor
      LABRAD = 1.2
      SCF = 0.45 / (LABRAD*ORAD)
      CALL PBLTM  (Z,Z, 0.5,0.5, Z, SCF,SCF, ERRIND, XF)
      CALL CHKINQ ('pbltm', ERRIND)
      CALL PSLMT  (XF, PCPRE)

C  set height for annotation text, apply scale factor directly,
C  since annotation text unaffected by transformation.
      CALL WCNPC (0.0, 0.0, NPCX,NPCY, NPCPWC)
      CALL PSATCH (NPCPWC * SCF)

C  annotation text alignments to be used:
C     algnh    algnv
C     ------   ------
C  1: center   base
C  2: left     cap
C  3: right    top
C  4: left     bottom
C  5: center   half
C  6: right    half
C  7: left     cap
      CALL SETVAL ('2,1,3,1,2,3,1', ALGNH)
      CALL SETVAL ('4,2,1,5,3,3,2', ALGNV)

C first draw labels and expected rectangles
      CALL PSPLCI (2)
      CALL PSTXCI (1)
      CALL PSTXAL (PACENT, PAHALF)
      CALL PSTXP (PRIGHT)

      DO 100 IX = 1,NUMSTR
C  txang = angle for annotation text string
         TXANG = IX*2*PI / NUMSTR
         COST = COS(TXANG)
         SINT = SIN(TXANG)
C  draw rectangle for expected position of annotation text string,
C     based on txang, ntxwd,ntxht
C  first draw bottom line along radius
         XA(1) = IRAD * COST
         XA(5) = XA(1)
         XA(4) = (IRAD + RECLNG) * COST

         YA(1) = IRAD * SINT
         YA(5) = YA(1)
         YA(4) = (IRAD + RECLNG) * SINT
C  unit up vector for rectangle  - this keeps left-bottom corner
C  on the radius for either horizontal or vertical path
         IF (VERT) THEN
            XVEC =  SINT
            YVEC = -COST
            INTSP = 0.5
         ELSE
            XVEC = -SINT
            YVEC =  COST
            INTSP = VTAB(PABASE) / VTAB(PATOP)
         ENDIF

         XA(2) = XA(1) + RECSHT*XVEC
         XA(3) = XA(4) + RECSHT*XVEC
         YA(2) = YA(1) + RECSHT*YVEC
         YA(3) = YA(4) + RECSHT*YVEC

         CALL PSLN (PLDOT)
         CALL PPL (5, XA,YA)

C  do baseline or centerline
         XB(1) = XA(1) + INTSP*RECSHT*XVEC
         XB(2) = XA(4) + INTSP*RECSHT*XVEC
         YB(1) = YA(1) + INTSP*RECSHT*YVEC
         YB(2) = YA(4) + INTSP*RECSHT*YVEC

         CALL PSLN (PLDASH)
         CALL PPL (2, XB, YB)

C  label each spoke
         WRITE (DIG2, '(A,I1)') '#', IX
         CALL PTX ( LABRAD*XA(4), LABRAD*YA(4), DIG2)
100   CONTINUE

      CALL PSATP (TXPATH)
      DO 200 IX = 1,NUMSTR
C  txang = angle for annotation text string
         TXANG = IX*2*PI / NUMSTR
         COST = COS(TXANG)
         SINT = SIN(TXANG)
C  txpos = requested annotation text position, based on txang,irad,
C    and annotation text alignment to be requested.
C  first get unit up,right direction
         IF (VERT) THEN
            XUP =  COST
            YUP =  SINT
            XRT =  SINT
            YRT = -COST
         ELSE
            XUP = -SINT
            YUP =  COST
            XRT =  COST
            YRT =  SINT
         ENDIF

         TXPOSX = IRAD * COST + HTAB(ALGNH(IX)) * XRT
     1                        + VTAB(ALGNV(IX)) * XUP
         TXPOSY = IRAD * SINT + HTAB(ALGNH(IX)) * YRT
     1                        + VTAB(ALGNV(IX)) * YUP
         CALL PSATAL (ALGNH(IX), ALGNV(IX))

         IF (IX .EQ. NGSTR) THEN
            TXANG = TXANG + (RNDINT(0,1)-0.5) * 0.5*PI/NUMSTR
            COST = COS(TXANG)
            SINT = SIN(TXANG)
            IF (VERT) THEN
               XUP =  COST
               YUP =  SINT
            ELSE
               XUP = -SINT
               YUP =  COST
            ENDIF
         ENDIF
         CALL PSATCU (XUP, YUP)
         CALL PATR (TXPOSX, TXPOSY, Z,Z, TXSTR)
C  next ix
200   CONTINUE

      END

C  *********************************************************
C  *                                                       *
C  *    SUBROUTINE 04.02.04.01/dsanst                      *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      SUBROUTINE DSANST (LSIZ, LANST)

C  DSANST accepts a list of annotation styles and displays at most 8
C  of them.

C marker type
      INTEGER    PPOINT,   PPLUS,   PAST,   POMARK,   PXMARK
      PARAMETER (PPOINT=1, PPLUS=2, PAST=3, POMARK=4, PXMARK=5)

      INTEGER    LSIZ, LANST(LSIZ), NUMDIS, IX, PERM(10), LDIS(10)

      REAL       XA(10),YA(10),XO(10),YO(10)
      REAL       NPCX,NPCY,NPCPWC, YINCR,YLOC, ANOFFX,ANOFFY

      CHARACTER  ANSTID*11

      CALL WCNPC (0.0, 0.0, NPCX,NPCY, NPCPWC)

      NUMDIS = MIN(8, LSIZ)
      CALL RNSET (NUMDIS, LSIZ, PERM)

C  build list of annotation styles to be displayed
      DO 100 IX = 1,NUMDIS
         LDIS(IX) = LANST(PERM(IX))
100   CONTINUE
C  and sort it
      CALL SRTIAR (NUMDIS, LDIS)

C  display in order and label all annotation styles in ldis
      YINCR = MIN (0.25, 1.0 / (LSIZ+1))
      YLOC  = 1 - YINCR
C  annotation offsets in WC
      ANOFFX = 0.2
      ANOFFY = 0.3 * YINCR
      DO 220 IX = 1,NUMDIS
C  mark reference and annotation point
         XO(IX) = 0.6
         YO(IX) = YLOC
         XA(IX) = XO(IX) - ANOFFX
         YA(IX) = YO(IX) - ANOFFY

         WRITE (ANSTID, '(A,I5,A)') 'Type#', LDIS(IX), ':'
         CALL VISLAB (ANSTID, 'R', 0.05,0.35,
     1                YLOC-0.3*YINCR, YLOC+0.3*YINCR)

         CALL PSANS (LDIS(IX))
         CALL PATR (XA(IX), YA(IX), ANOFFX*NPCPWC, ANOFFY*NPCPWC, 'abc')
         YLOC = YLOC-YINCR
220   CONTINUE
      CALL PSMK (PAST)
      CALL PPM  (NUMDIS, XA,YA)
      CALL PSMK (POMARK)
      CALL PPM  (NUMDIS, XO,YO)

      END


text:
Address(M Street) runs for only two blocks between Duke Drive on the south side and 8th Street East Eighth Street on the north side, with a cute traffic circle where it crosses Colgate.

The street is exclusively residential until it meets East Eighth where the Davis Manor Center Orange Mall is on the right.

One notable house is on the outer corner of Duke and M, where at least 14 cats hang out, fed by the man who lives there. Another is the Robot Rocket Residence.

20071026 12:59:31 nbsp I just moved about 3 houses up from there, and the cats seem to be gone Users/ThomasAnastasio

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
       FUNCTION K42R(JPICK,MOVEN,ISCORE,JSCORE,RANDOM, JA)
C BY OTTO BORUFSEN
C TYPED FROM FORTRAN BY AX, 1/25/79
      DIMENSION MHIST(2,2)
      k42r=ja    ! Added 7/27/93 to report own old value
C INITIALIZE FIRST MOVE
      IF(MOVEN.NE.1)GOTO 20
         L3MOV=0
         L3ECH=0
         IDEF=0
         ICOOP=0
         IPICK=0
         DO 10 I=1,2
         DO 10 J=1,2
10       MHIST(I,J)=0
         GO TO 500
20    IF(MOVEN.EQ.2)GOTO 25
C UPDATE MOVE HISTORY
      MHIST(I2PCK+1,JPICK+1)=MHIST(I2PCK+1,JPICK+1)+1
25    IF(IDEF.EQ.0)GOTO 30
C OPPONENT HAS BEEN PROVED "RANDOM" OR
C "DEFECTIVE",I DEFECT FOR 25 MOVES
      K42R=1
      GO TO 100
30    IF(IPICK.EQ.0.OR.JPICK.EQ.0)GOTO 40
C MUTUAL DEFECTIONS ON LAST MOVE.
      L3MOV=L3MOV+1
      IF(L3MOV.LT.3)GOTO 50
C MUTUAL DEFECTION ON
C LAST THREE MOVES.
C I COOPERATE ONCE ON NEXT MOVE.
      K42R=0
      L3MOV=0
      L3ECH=0
      GOTO 100
C ONE (OR BOTH) COOPERATED ON LAST MOVE.
40    L3MOV=0
      IF(IPICK.EQ.JPICK)GOTO 45
      IF(JPICK.NE.I2PCK.OR.IPICK.NE.J2PCK)GOTO 45
C ECHO-EFFECT ON LAST MOVE.
      L3ECH=L3ECH+1
      IF(L3ECH.LT.3)GOTO 50
C ECHO-EFFECT ON LAST THREE MOVES.
C MY NEXT DEFECTION WILL BE SUBSTITUTED
C BY A COOPERATION.
      L3ECH=0
      L3MOV=0
      ICOOP=1
      GOTO 50
45    L3ECH=0
C PLAY 'TIT FOR TAT' AS MAIN RULE.
50    K42R=JPICK
100   IF(MOD(MOVEN-2,25).NE.0.OR.MOVEN.EQ.2)GOTO 650
C ON EVERY 25 MOVES:
C CHECK IF OPPONENT SEEMS TO BE
C 'RANDOM' OR 'DEFECTIVE'.
      IDEF=0
      JNCOP=MHIST(1,1)+MHIST(2,1)
C IS OPPONENT 'RANDOM'?
      IF(JNCOP.GT.17)GOTO 155
      IF(JNCOP.LT.8)GOTO 130
      IF(100*MHIST(1,1)/JNCOP.LT.70)IDEF=1
      GOTO 155
C IS OPPONENT 'DEFECTIVE'?
130   IF(JNCOP.LT.3)IDEF=1
155   DO 160 I=1,2
      DO 160 J=1,2
160   MHIST(I,J)=0
      IF(IDEF.EQ.0)GOTO 650
C OPPONENT SEEMS TO BE
C 'RANDOM' OR 'DEFECTIVE'.
C I DEFECT FOR NEXT 25 MOVES.
      ICOOP=0
      L3MOV=0
      L3ECH=0
      GOTO 600
C I COOPERATE.
500   K42R=0
      GOTO 650
C I DEFECT.
600   K42R=1
650   IF(ICOOP.EQ.0.OR.K42R.EQ.0)GOTO 660
         ICOOP=0
         K42R=0
660   I2PCK=IPICK
      J2PCK=JPICK
      IPICK=K42R
      RETURN
      END

max_stars_repo_name: Axelrod-Python/TourExec

max_issues_repo_name: Axelrod-Python/TourExec

max_forks_repo_name: Axelrod-Python/TourExec


text:
MODULE USR_SMlib_Band_Gauss_Solver
  !
  !!  This module is part of SMLIB v. 1.1.  It contains procedures needed
  !!  to interface smlib with Lapack"s band matrix Gaussian elimination routines
  !!  gbfa and gbsl.  If Lapack is not installed on your computer, these routines
  !!  are available from NETLIB.
  !
  !! ------------------------------------------------------------------------
  !
  !!  Copyright (C) 1998 Ernst A. Meese
  !!    Refer to the file copyright.doc for details and important disclaimer.
  !
  !!  Created: May 1998
  !
  !! ------------------------------------------------------------------------
  !!  Module contents
  !
  !!   Interfaces:
  !!      LU_Factor - Generic interface for LU factorisation of a matrix.
  !!      SOLVE     - Generic interface for the solution of a linear equation system.
  !
  !!   Procedures:
  !!      LU_Factor_MSR_to_Band_LU - Factor an MSR matrix, assuming it to have
  !!           some band structure, into a band LU format USEd by Lapack.
  !!      SOLVE_Band_LU - Solve a linear equation system, given a band LU
  !!           factorisation of the matrix.
  !
  !! -----------------------------------------------------------------------
  !
  !!  We use the modules defining abstract datatypes for the matrices
  !
  use USR_SMlib_Precision         !!((16-B-USR_SMlib_Precision.f90))
  use USR_SMlib_Matrix_Arithmetic !!((19-B-USR_SMlib_Matrix_Arithmetic.f90))
  use LIB_LINPACK_Rsp             !!((04-B-LIB_LINPACK_Rsp.f))
  use LIB_LINPACK_Rdp             !!((04-B-LIB_LINPACK_Rdp.f))
  IMPLICIT NONE

  PRIVATE

  INTERFACE LU_Factor
     MODULE PROCEDURE LU_Factor_MSR_to_Band_LU
  END INTERFACE

  INTERFACE Solve
     MODULE PROCEDURE SOLVE_Band_LU
  END INTERFACE

  PUBLIC :: LU_Factor, Solve

CONTAINS

  SUBROUTINE LU_Factor_MSR_to_Band_LU (A, LU, inform )
    !
    type(MSR),            intent(IN)    :: A
    type(Band_LU),        intent(INOUT) :: LU
    integer,    optional, intent(OUT)   :: inform
    !
    integer :: m, ml, mu, i, j, k, jl, ju, info

    mu = 0 ; ml = 0
    DO i = 1, A % N
       ju = MAX( i, MAXVAL( A % JA( A % JA(i):A % JA(i+1)-1 ) ) )
       jl = MIN( i, MINVAL( A % JA( A % JA(i):A % JA(i+1)-1 ) ) )
       ml = MAX( mu, ju - i )
       mu = MAX( ml, i - jl )
    END DO

    IF ( Allocated_Matrix(LU) ) THEN
       IF ( LU % Upper_band_width /= mu .or. LU % Lower_band_width /= ml .or. size(LU % abd,2) /= A % N ) &
            Call Deallocate_Matrix ( LU )
    END IF
    IF (.not.Allocated_Matrix(LU)) Call Allocate_matrix ( LU, ml, mu, A % N )

    LU % abd = 0.0_prec
    m = ml + mu + 1
    DO i = 1, A % N
       DO j = A % JA (i), A % JA (i+1) - 1
          k = A % JA ( j ) - i + m
          LU % abd(k,i) = A % A ( j )
       END DO
    END DO
    LU % abd(m,1:A%N) = A % A (1:A%N)

    !provide overloaded GBFA instead [waw]
        CALL GBFA( LU % ABD, SIZE(LU%ABD,1), LU%N, ML, MU, LU%IPVT, INFO)

        !commented out by [waw]
        !SELECT CASE ( KIND(LU % ABD(1,1)) )
    !CASE ( single_prec )
    !!   CALL SGBFA( REAL(LU % ABD,single_prec), SIZE(LU%ABD,1), &
        !!     LU%N, ML, MU, LU%IPVT, INFO)
    !CASE ( double_prec )
    !!   CALL DGBFA( REAL(LU % ABD,double_prec), SIZE(LU%ABD,1), &
        !!     LU%N, ML, MU, LU%IPVT, INFO)
    !CASE DEFAULT
    !!   STOP "Real type has to be either single or double precision for band Gaussian solver."
    !END SELECT
    IF( PRESENT( inform ) ) inform = INFO

  END SUBROUTINE LU_Factor_MSR_to_Band_LU

  SUBROUTINE SOLVE_Band_LU( LU, b, x )
    TYPE(Band_LU),                 intent(in) :: LU
    real(prec), dimension(LU % N), intent(in) :: b
    real(prec), dimension(LU % N), intent(out):: x

    x = b

        !provide overloaded GBSL instead [waw]
        CALL GBSL( LU%ABD, SIZE(LU%ABD,1), &
           LU%N, LU%Lower_band_width, LU%Upper_band_width, LU%IPVT, &
           x , 1 )

        !commented out by [waw]
        !SELECT CASE ( KIND(LU % ABD(1,1)) )
    !CASE ( single_prec )
    !!   CALL SGBSL( REAL(LU%ABD,single_prec), SIZE(LU%ABD,1), &
        !!   LU%N, LU%Lower_band_width, LU%Upper_band_width, LU%IPVT, &
        !!   REAL(x,single_prec), 1 )
    !CASE ( double_prec )
    !!   CALL DGBSL( REAL(LU%ABD,double_prec), SIZE(LU%ABD,1), &
        !!   LU%N, LU%Lower_band_width, LU%Upper_band_width, LU%IPVT, &
        !!   REAL(x,double_prec), 1 )
    !CASE DEFAULT
    !!   STOP "Real type has to be either single or double precision for band Gaussian solver."
    !END SELECT

  END SUBROUTINE SOLVE_Band_LU

END MODULE USR_SMlib_Band_Gauss_Solver

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack


text:
Presented by UC Davis’ amazing allfemale a cappella group, The Spokes The Spokes, HellaCappella is the most highly anticipated a cappella concert of the year!

For almost a decade, HellaCappella’s audience has continued to grow, attracting music lovers across generations from the greater Davis and Sacramento community. In 2012, The Spokes raised the stakes and hosted HellaCappella at UC Davis’ ARC Pavilion ARC Pavilion, selling over 2000 seats. Then, in 2013 for the first time ever and as the first studentrun group to do so, The Spokes proudly presented the ninth annual HellaCappella at UC Davis’ Mondavi Center Mondavi Center in the beautiful, premier arts facility, Jackson Hall Jackson Hall, where they had a sold out show with 1,801 audience members.

>>> MARK YOUR CALENDARS! HellaCappella 2014 will be Friday, April 18th, 2014 in Jackson Hall in the Mondavi Center! <<<

CLICK https://www.facebook.com/events/791341117562075/ HERE FOR THE FACEBOOK EVENT! CLICK http://www.mondaviarts.org/events/event.cfm?event_id1427&season2013 HERE TO BUY TICKETS! This years event will feature ASSIGNED SEATING, so get your tickets as soon as possible!

HellaCappella 2014 will be emceed by Davis own allmale group The Afterglow The Afterglow! Our line up features new West Coast a cappella groups:

· Mens Octet, UC Berkeley
· On the Rocks, University of Oregon
· The Stanford Harmonics, Stanford University
· Dil Se, UC Berkeley
· Cloud 9, UC Santa Cruz
· The The Afterglow The Afterglow, UC Davis
· The Spokes The Spokes, UC Davis

Emceed by UC Davis’ own Birdstrike Theater Birdstrike Comedy Theatre, HellaCappella 2013s full lineup was as exciting and diverse as ever, featuring:

· Mens Octet, UC Berkeley
· On the Rocks, University of Oregon
· The Golden Overtones, UC Berkeley
· DeCadence, UC Berkeley
· Dil Se, UC Berkeley
· Cloud 9, UC Santa Cruz
· The Liquid Hotplates Liquid Hotplates, UC Davis
· The Spokes The Spokes, UC Davis

The Spokes are the only allfemale group at UC Davis. Founded in 2004, the group currently consists of 14 talented and dedicated girls whose shared passion for music and performing shines when they take the stage. Their repertoire includes a wide range of music genres, from top 40 (Never Ever Getting Back Together by Taylor Swift) to alternative (Seven Nation Army by The White Stripes) to R&B (I’ll Be Waiting by Adele) to classics (Sweet Caroline by Neil Diamond). The Spokes’ love for singing can be heard in their latest album, Shades of Pink: Volume 1, which is available on iTunes, Spotify, Pandora, and Loudr. http://www.youtube.com/watch?vLuJouRBtzB8 Here is a video of The Spokes singing their rendition of Us Against the World by Play at HellaCappella 2013!

In a continuous effort to give back to the community, The Spokes committed to raising funds for Teach For America, the leading national teaching corps dedicated to working at every level of education, policy and other professions, to ensure that all children can receive an excellent education. HellaCappella 2012 raised over $3000 for The American Cancer Society, and HellaCappella 2013 raised $2,000 to support TFA’s educational efforts!

HellaCappella 2014 will benefit http://www.womensempowerment.org/ The Womens Empowerment Center in Sacramento!

We hope to see you there next year!
For more information on The Spokes:

https://www.thedavisspokes.com Spokes website!
https://www.facebook.com/pages/TheSpokes/136852192894?frefts Spokes Facebook!
email: the.spokes@gmail.com

For more information on The Mondavi Center:
www.mondaviarts.org

 Pictures


text:
Describe Users/DanielPotter here.

20100728 11:53:46 nbsp Hi Daniel, Welcome to the Wiki!  Just to make sure you catch the replies, a couple of people were wondering why you found the message on the Igor Rimmer Bench offensive. Users/TomGarberson

20140412 14:24:54 nbsp Hi Daniel,

You are correct.  The message on the bench could offend a lot of people, but the reasons are subtle.
Im not sure who wrote the message so I cant be sure about the cause.
It could be poor wording by a nonnative English speaker, lack of poetic license by a
native English speaker, or it could have been deliberate.  In any case, you may want to have it changed.

The four sources of the problem are:

1) Igor Rimmers name should not be in all capital letters. It would look insensitive
on any memorial, Jewish or not.  The only way it may have been appropriate
is if nothing else were written afterword, just In loving memory of IGOR RIMMER
But, in this case it should read Igor Rimmer.

2) The ampersand symbol &.  They should actually have spelled out the word
and.  It would look insensitive on any memorial, Jewish or not.

2) The word other would look insensitive on any memorial, Jewish or not.  It should read all Jews

3)The exclamation point after Never again! would also look insensitive on any memorial, Jewish
or not, in this time period.  If the message had been placed there shortly after World War II,
it may have been due to the time period.  They should have used Never again. with a period
instead of an exclamation point as in the issue is settled.

Take care. Users/Maryilyn


text:
!=======================================================================
!  Aloha_ROOTGR, Subroutine
!
!  Determines root growth
!-----------------------------------------------------------------------
!  Revision history
!
!  Written
!  09/1988 Modified by E. Alocilja & B. Baer
!  04/1989 Modified by T. Jou
!  02/08/1993 PWW Added switch block,  Header revision and minor changes
!  06/20/1994 JTR, BDB Simplified the RLNEW calculation and slowed
!                  the growth of roots in deeper soils.
!  06/25/2017 CHP Adpated for CSM v4.6
!  09/05/2020 JVJ Stage changed for inclusion in Overview
!-----------------------------------------------------------------------
!                         DEFINITIONS
!
!  L,L1   : Loop counter
!  RLDF() : A root length density factor for soil layer L used to calculate
!           new root growth distribution - unitless
!  RLNEW  : New root length to be added to the total root system length -
!           cm.  root per sq. cm. ground
!  RNFAC  : Zero to unity factor describing mineral N availability effect on
!           root growth in Layer L
!  RNLF   : Intermediate factor used to calculate distribution of new root
!           growth in the soil - unitless value between 0 and 1
!  SWDF   : Soil water deficit factor for Layer L used to calculate root
!           growth and water uptake - unitless value between 0 and 1
!  TRLDF  : An intermediate calculation used to calculate distribution of
!           new root growth in soil
!=======================================================================

      SUBROUTINE Aloha_ROOTGR (CONTROL,
     &     CUMDTT, DTT, GRORT, ISTAGE, ISWITCH, NO3, NH4,     !Input
     &     SOILPROP, SW, SWFAC,                               !Input
     &     RLV, RTDEP, RTWT)                                  !Output

      USE Aloha_mod
      IMPLICIT  NONE
      SAVE

      INTENT(IN) :: CONTROL,
     &     CUMDTT, DTT, GRORT, ISTAGE, ISWITCH, NO3, NH4,     !Input
     &     SOILPROP, SW, SWFAC                                !Input
      INTENT(OUT):: RLV, RTDEP, RTWT

      CHARACTER ISWWAT*1
      REAL      RNFAC,RLNEW,SWDF,TRLDF,RNLF,GRORT, RTWT, RTDEP, CUMDEP
      REAL  PLTPOP, DTT, CUMDTT, DEPFAC, SWFAC
      REAL, DIMENSION(NL) :: RLDF, RLV, DLAYR, ESW, SW, NO3, NH4
      INTEGER   L,L1, ISTAGE, DYNAMIC, NLAYR
      LOGICAL FIRST

      TYPE (ControlType) CONTROL
      TYPE (SoilType) SOILPROP
      TYPE (SwitchType) ISWITCH

      DLAYR = SOILPROP % DLAYR
      DYNAMIC = CONTROL % DYNAMIC

!=======================================================================
      SELECT CASE (DYNAMIC)
!=======================================================================
      CASE (RUNINIT, SEASINIT)
!=======================================================================
      RTWT   =  0.0
      RTDEP  = Planting % SDEPTH               ! Rooting depth = seeding depth (cm)

      PLTPOP = Planting % PLTPOP
      NLAYR  = SOILPROP % NLAYR

      FIRST = .TRUE.

      ISWWAT = ISWITCH % ISWWAT
      IF (ISWWAT .NE. 'Y') RETURN
      DO L = 1, NLAYR
        ESW(L) = SOILPROP % DUL(L) - SOILPROP % LL(L)
      ENDDO

!=======================================================================
      CASE (RATE)
!=======================================================================
!     from phenology
      SELECT CASE (ISTAGE)
        CASE (1)
          IF (ISWWAT .NE. 'Y') RETURN
          RTDEP  = RTDEP + 0.01*DTT     ! Depth of root (f) DTT
        CASE (12)                        !CASE (9) JVJ Value changed because 2 stages in vegetative phase and one stage in reproductive phase were included
          RTDEP  = RTDEP + 0.01*DTT     ! Depth of root (f) DTT
      END SELECT

      IF (ISWWAT .NE. 'Y') RETURN
      IF (GRORT < 0.0001) RETURN
C
C     The small differences between root length/weight ratios used in earlier
C     models were insignificant considering the uncertainty of the value
C     and the uncertainty of loss of assimilate by exudation and respiration.
C     A compromise value of 0.98 was choosen for all crops.
C
      RLNEW  = GRORT * Species % RLWR * PLTPOP              ! Read in from species value

      CUMDEP = 0.0
      RNFAC  = 1.0
      SWDF = 1.0

      L = 0
      DO WHILE ((CUMDEP .LT. RTDEP) .AND. (L .LT. NLAYR))
        L = L + 1
        CUMDEP = CUMDEP + DLAYR(L)
        IF (ISWITCH % ISWWAT .NE. 'N') THEN
          IF (SW(L) - SOILPROP % LL(L) .LT. 0.25*ESW(L)) THEN
            SWDF = 4.0 * (SW(L) - SOILPROP % LL(L)) / ESW(L)
            IF (SWDF .LT. 0.0) THEN
              SWDF = 0.0
            ENDIF
          ELSE
            SWDF = 1.0
          ENDIF
C
C         Made all crops so that RNFAC is constrained between 0.01 and 1.0;
C         on page 94 of Jones & Kiniry book the minimum is 0.01. - WTB
C
          IF (ISWITCH % ISWNIT .NE. 'N') THEN
            RNFAC = 1.0 - (1.17*EXP(-0.15*(NO3(L)+NH4(L))))
            RNFAC = AMAX1 (RNFAC,0.01)
          ENDIF
        ENDIF
        RLDF(L) = AMIN1(SWDF,RNFAC) * SOILPROP % WR(L) * DLAYR(L)
      END DO

      L1 = L
C
C     The following changes were made to simplify the code and make the model
C     more generic. It also takes into account some newer data provided by
C     Julio Dardenelli of Argentina.  For the first time the ceres model
C     restricts the rate of downward movement of roots with the soil property
C     -- root weighting factor -- to account for greater difficulty in growing
C     downward in hard soil. Changes made by JTR 6/16/94.
C
      DEPFAC = SQRT(SOILPROP % WR(L) * AMIN1(SWFAC * 2.0, SWDF))
      IF (CUMDTT .LT. 275.0) THEN                                ! JTR 6/17/94
         RTDEP = RTDEP + DTT * 0.1 * DEPFAC
       ELSE
         RTDEP = RTDEP + DTT * 0.2 * DEPFAC
      ENDIF

      RTDEP    = AMIN1 (RTDEP, SOILPROP % DS(NLAYR))
      RLDF(L1) = RLDF(L1)*(1.0-(CUMDEP-RTDEP)/DLAYR(L1))
      TRLDF    = 0.0

      DO  L = 1, L1
         TRLDF = TRLDF+RLDF(L)
      END DO

      IF (TRLDF .GE. RLNEW*0.00001) THEN
         RNLF = RLNEW/TRLDF
         DO L = 1, L1
            RLV(L) = RLV(L) + RLDF(L)*RNLF/DLAYR(L)-0.005*RLV(L)
            RLV(L) = AMAX1 (RLV(L),0.0)
            RLV(L) = AMIN1 (RLV(L),4.0)
         END DO
      ENDIF

!=======================================================================
!     Integration
!-----------------------------------------------------------------------
      CASE (INTEGR)
!=======================================================================
      IF (ISWWAT .NE. 'Y') RETURN

!     Initialize RLV at stage 1
      IF (ISTAGE == 1 .AND. FIRST) THEN
        FIRST = .FALSE.
        CUMDEP = 0.0
        DO L = 1, NLAYR
           CUMDEP = CUMDEP + SOILPROP % DLAYR(L)
           RLV(L) = 0.20*PLTPOP/DLAYR (L)
           IF (CUMDEP .GT. RTDEP) EXIT
        END DO

        RLV(L) = RLV(L)*(1.0-(CUMDEP-RTDEP)/DLAYR(L))
        L1 = L + 1
        IF (L1 .LT. NLAYR) THEN
           DO L = L1, NLAYR
              RLV(L) = 0.0
           END DO
        ENDIF
      ENDIF
!=======================================================================
      END SELECT
!=======================================================================
      RETURN
      END SUBROUTINE Aloha_ROOTGR
!=======================================================================
!=======================================================================

ext: for


text:
C     I is indirectly modified by many procedures via a hidden common

      program hidden

      call init

      call incr(1)

      call double

      call printi

      end

      subroutine init
      common /foo/i

      i = 0

      end

      subroutine incr(j)
      common /foo/i

      i = i + j

      end

      subroutine double
      common /foo/i

      i = 2 * i

      end

      subroutine printi
      common /foo/i

      print *, i

      end


text:
      SUBROUTINE SLARFTMOD( DIRECT, STOREV, N, K, V, LDV, TAU,
     $     T, LDT, DEPTH )
*
*  -- LAPACK auxiliary routine (version 3.4.2) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     September 2012
*
*     .. Scalar Arguments ..
      CHARACTER          DIRECT, STOREV
      INTEGER            K, LDT, LDV, N
*     ..
*     .. Array Arguments ..
      REAL               T( LDT, * ), TAU( * ), V( LDV, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, J, PREVLASTV, LASTV
*     ..
*     .. External Subroutines ..
      EXTERNAL           SGEMV_OMP_TASK, STRMV_OMP_TASK
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      EXTERNAL           LSAME
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( N.EQ.0 )
     $   RETURN
*
      IF( LSAME( DIRECT, 'F' ) ) THEN
         PREVLASTV = N
         DO I = 1, K
            PREVLASTV = MAX( I, PREVLASTV )
            IF( TAU( I ).EQ.ZERO ) THEN
*
*              H(i)  =  I
*
               DO J = 1, I
                  T( J, I ) = ZERO
               END DO
            ELSE
*
*              general case
*
               IF( LSAME( STOREV, 'C' ) ) THEN
*                 Skip any trailing zeros.
                  DO LASTV = N, I+1, -1
                     IF( V( LASTV, I ).NE.ZERO ) EXIT
                  END DO
                  DO J = 1, I-1
                     T( J, I ) = -TAU( I ) * V( I , J )
                  END DO
                  J = MIN( LASTV, PREVLASTV )
*
*                 T(1:i-1,i) := - tau(i) * V(i:j,1:i-1)**T * V(i:j,i)
*
                  CALL SGEMV_OMP_TASK( 'Transpose', J-I, I-1, -TAU( I ),
     $                        V( I+1, 1 ), LDV, V( I+1, I ), 1, ONE,
     $                        T( 1, I ), 1, DEPTH )
               ELSE
*                 Skip any trailing zeros.
                  DO LASTV = N, I+1, -1
                     IF( V( I, LASTV ).NE.ZERO ) EXIT
                  END DO
                  DO J = 1, I-1
                     T( J, I ) = -TAU( I ) * V( J , I )
                  END DO
                  J = MIN( LASTV, PREVLASTV )
*
*                 T(1:i-1,i) := - tau(i) * V(1:i-1,i:j) * V(i,i:j)**T
*
                  CALL SGEMV_OMP_TASK( 'No transpose', I-1, J-I,
     $                 -TAU( I ), V( 1, I+1 ), LDV, V( I, I+1 ), LDV,
     $                        ONE, T( 1, I ), 1, DEPTH )
               END IF
*
*              T(1:i-1,i) := T(1:i-1,1:i-1) * T(1:i-1,i)
*
               CALL STRMV_OMP_TASK( 'Upper', 'No transpose',
     $              'Non-unit', I-1, T,
     $              LDT, T( 1, I ), 1, DEPTH )
               T( I, I ) = TAU( I )
               IF( I.GT.1 ) THEN
                  PREVLASTV = MAX( PREVLASTV, LASTV )
               ELSE
                  PREVLASTV = LASTV
               END IF
            END IF
         END DO
      ELSE
         PREVLASTV = 1
         DO I = K, 1, -1
            IF( TAU( I ).EQ.ZERO ) THEN
*
*              H(i)  =  I
*
               DO J = I, K
                  T( J, I ) = ZERO
               END DO
            ELSE
*
*              general case
*
               IF( I.LT.K ) THEN
                  IF( LSAME( STOREV, 'C' ) ) THEN
*                    Skip any leading zeros.
                     DO LASTV = 1, I-1
                        IF( V( LASTV, I ).NE.ZERO ) EXIT
                     END DO
                     DO J = I+1, K
                        T( J, I ) = -TAU( I ) * V( N-K+I , J )
                     END DO
                     J = MAX( LASTV, PREVLASTV )
*
*                    T(i+1:k,i) = -tau(i) * V(j:n-k+i,i+1:k)**T * V(j:n-k+i,i)
*
                     CALL SGEMV_OMP_TASK( 'Transpose', N-K+I-J,
     $                    K-I, -TAU( I ),
     $                    V( J, I+1 ), LDV, V( J, I ), 1, ONE,
     $                    T( I+1, I ), 1, DEPTH )
                  ELSE
*                    Skip any leading zeros.
                     DO LASTV = 1, I-1
                        IF( V( I, LASTV ).NE.ZERO ) EXIT
                     END DO
                     DO J = I+1, K
                        T( J, I ) = -TAU( I ) * V( J, N-K+I )
                     END DO
                     J = MAX( LASTV, PREVLASTV )
*
*                    T(i+1:k,i) = -tau(i) * V(i+1:k,j:n-k+i) * V(i,j:n-k+i)**T
*
                     CALL SGEMV_OMP_TASK( 'No transpose', K-I, N-K+I-J,
     $                    -TAU( I ), V( I+1, J ), LDV, V( I, J ), LDV,
     $                    ONE, T( I+1, I ), 1, DEPTH )
                  END IF
*
*                 T(i+1:k,i) := T(i+1:k,i+1:k) * T(i+1:k,i)
*
                  CALL STRMV_OMP_TASK( 'Lower', 'No transpose',
     $                 'Non-unit', K-I,
     $                 T( I+1, I+1 ), LDT, T( I+1, I ), 1, DEPTH )
                  IF( I.GT.1 ) THEN
                     PREVLASTV = MIN( PREVLASTV, LASTV )
                  ELSE
                     PREVLASTV = LASTV
                  END IF
               END IF
               T( I, I ) = TAU( I )
            END IF
         END DO
      END IF
      RETURN
*
*     End of SLARFT
*
      END


text:
c
c
c     ==================================================================
      subroutine rpt3(ixyz,icoor,imp,maxm,meqn,mwaves,maux,mbc,mx,
     &                  ql,qr,aux1,aux2,aux3,asdq,
     &                  bmasdq,bpasdq)
c     ==================================================================
c
c     # Riemann solver in the transverse direction.
c     # This is a dummy routine that returns zeros and is only intended
c     # to illustrate the format of this routine.  See various example
c     # directories for better examples.

c     # This dummy routine can be used if transverse solves are not being
c     # used, i.e. if method(3) = 0.
c     #
c     # On input,
c
c     #    ql,qr is the data along some one-dimensional slice, as in rpn3
c     #         This slice is
c     #             in the x-direction if ixyz=1,
c     #             in the y-direction if ixyz=2, or
c     #             in the z-direction if ixyz=3.
c     #    asdq is an array of flux differences (A^*\Dq).
c     #         asdq(:,i) is the flux difference propagating away from
c     #         the interface between cells i-1 and i.
c     #    Note that asdq represents B^*\Dq if ixyz=2 or C^*\Dq if ixyz=3.
c
c     #    ixyz indicates the direction of the original Riemann solve,
c     #         called the x-like direction in the table below:
c
c     #               x-like direction   y-like direction   z-like direction
c     #      ixyz=1:        x                  y                  z
c     #      ixyz=2:        y                  z                  x
c     #      ixyz=3:        z                  x                  y
c
c     #    icoor indicates direction in which the transverse solve should
c     #         be performed.
c     #      icoor=2: split in the y-like direction.
c     #      icoor=3: split in the z-like direction.
c
c     #    For example,
c     #      ixyz=1, icoor=2 means asdq=A^*\Dq, and should be split in y into
c     #                        bmasdq = B^-A^*\Dq,
c     #                        bpasdq = B^+A^*\Dq.
c     #
c     #      ixyz=2, icoor=2 means asdq=B^*\Dq, and should be split in z into
c     #                        bmasdq = C^-B^*\Dq,
c     #                        bpasdq = C^+B^*\Dq.
c
c     #    The parameter imp is generally needed only if aux
c     #    arrays are being used, in order to access the appropriate
c     #    variable coefficients.
c
c
c
      implicit real*8(a-h,o-z)
      dimension     ql(meqn, 1-mbc:maxm+mbc)
      dimension     qr(meqn, 1-mbc:maxm+mbc)
      dimension   asdq(meqn, 1-mbc:maxm+mbc)
      dimension bmasdq(meqn, 1-mbc:maxm+mbc)
      dimension bpasdq(meqn, 1-mbc:maxm+mbc)
      dimension   aux1(maux, 1-mbc:maxm+mbc, 3)
      dimension   aux2(maux, 1-mbc:maxm+mbc, 3)
      dimension   aux3(maux, 1-mbc:maxm+mbc, 3)
c
      do i = 2-mbc, mx+mbc
         do m=1,meqn
            bmasdq(m,i) = 0.d0
            bpasdq(m,i) = 0.d0
         enddo
      enddo
c
c
      return
      end

max_stars_repo_name: ClimateImpactLab/classic

max_issues_repo_name: ClimateImpactLab/classic

max_forks_repo_name: cjvogl/classic