! Copyright (c) 2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

!OPTIONS: -fopenmp

! 2.15.2 threadprivate Directive
! The threadprivate directive specifies that variables are replicated,
! with each thread having its own copy. When threadprivate variables are
! referenced in the OpenMP region, we know they are already private to
! their threads, so no new symbol needs to be created.

!DEF: /mm Module
module mm
  !$omp threadprivate (i)
contains
  !DEF: /mm/foo PUBLIC (Subroutine) Subprogram
  subroutine foo
    !DEF: /mm/foo/a ObjectEntity INTEGER(4)
    integer :: a = 3
    !$omp parallel
    !REF: /mm/foo/a
    a = 1
    !DEF: /mm/i PUBLIC (Implicit, OmpThreadprivate) ObjectEntity INTEGER(4)
    !REF: /mm/foo/a
    i = a
    !$omp end parallel
    !REF: /mm/foo/a
    print *, a
    block
      !DEF: /mm/foo/Block2/i ObjectEntity REAL(4)
      real i
      !REF: /mm/foo/Block2/i
      i = 3.14
    end block
  end subroutine foo
end module mm
!DEF: /tt MainProgram
program tt
  !REF: /mm
  use :: mm
  !DEF: /tt/foo (Subroutine) Use
  call foo
end program tt


text: ${PRJ_DIR}/hw/ip/rv_core_ibex/rtl/rv_core_ibex.sv

ext: f


text:
      real*8 function Iapqterm1(r)
      implicit none
      real*8 pa,pb,rb,t,z,g,r
      real*8 x,a,p,q,beta,ra
      common /xapqbeta/ x,a,p,q,beta,ra
      real*8 isr_tilde
      external isr_tilde
      pa=-2d0/3d0
      pb=beta-1d0
      t=ra*r**(1d0/(pa+1))
      z=t*(1d0-x)+x
      g=a*dexp(p*(1d0-x/z))*dexp(q*dsqrt(1d0-x/z))*isr_tilde(z,beta)
      g=g*z**(-1d0/3d0)
      Iapqterm1=g*(1d0-t)**pb*ra**(pa+1d0)/(pa+1d0)
      end

      real*8 function Iapqterm2(r)
      implicit none
      real*8 pa,pb,rb,t,z,g,r
      real*8 x,a,p,q,beta,ra
      common /xapqbeta/ x,a,p,q,beta,ra
      real*8 isr_tilde
      external isr_tilde
      pa=-2d0/3d0
      pb=beta-1d0
      rb=1d0-ra
      t=1d0-rb*r**(1d0/(pb+1))
      z=t*(1d0-x)+x
      g=a*dexp(p*(1d0-x/z))*dexp(q*dsqrt(1d0-x/z))*isr_tilde(z,beta)
      g=g*z**(-1d0/3d0)
      Iapqterm2=g*t**pa*rb**(pb+1d0)/(pb+1d0)
      end

      real*8 function Iapq(xpar,apar,ppar,qpar,betapar)
      implicit none
      real*8 Iapqterm1,Iapqterm2
      external Iapqterm1,Iapqterm2
      real*8 dgauss
      external dgauss
      real*8 xpar,apar,ppar,qpar,betapar
      real*8 res1,res2
      real*8 x,a,p,q,beta,ra
      common /xapqbeta/ x,a,p,q,beta,ra
      real*8 intlow,intupp,inteps
      x = xpar
      a = apar
      p = ppar
      q = qpar
      beta = betapar
      ra=0.5d0
      intlow = 0d0
      intupp = 1d0
      inteps = 1d-5
      res1 = dgauss(Iapqterm1,intlow,intupp,inteps)
      res2 = dgauss(Iapqterm2,intlow,intupp,inteps)
      Iapq = res1+res2
      end

c     This function calculate the reduced structure function, with energy
c     fraction given by "x", at scale "Qsquare"
      real*8 function eepdf_tilde_calc(x,Qsquare,n,partonid,beamid)
      implicit none
      real*8 x,Qsquare
      real*8 me
      data me /0.511d-3/
      real*8 PI
      real*8 alphaem
c     In Gmu scheme
      data alphaem/0.007562397d0/
      real*8 beta
      integer n,partonid,beamid
      real*8 a11,a01,a10,a00
      real*8 p1,q1,p2,q2,p3,q3,p4,q4
      data a11/0.8698d0/,a01/0.2863d0/,a10/0.2853d0/,a00/0.1092d0/
      data p1/-901.2d0/,q1/-19.30d0/,p2/-916.5d0/,q2/-18.70d0/
      data p3/-899.2d0/,q3/-18.03d0/,p4/-918.7d0/,q4/-17.48d0/
      real*8 Iapq,isr_tilde
      real*8 res
      data res/0d0/

      PI=4.D0*DATAN(1.D0)
      beta = alphaem/PI * (dlog(Qsquare/me/me)-1d0)

c     electron beam
      if (beamid .eq. 11) then
c     other partons are zero
          if (partonid .ne. 11) then
             res = 0d0
          else
             if (n .eq. 1) then
                res = dsqrt(a11)*isr_tilde(x,beta)
             else if (n .eq. 2) then
                res = Iapq(x,a01/dsqrt(a11),p1,q1,beta)
             else if (n .eq. 3) then
                res = dsqrt(a11)*isr_tilde(x,beta)
             else if (n .eq. 4) then
                res = Iapq(x,dsqrt(a00),p3,q3,beta)
             else
                res = 0d0
             endif
          endif
      else if (beamid .eq. -11) then
          if (partonid .ne. -11) then
              res = 0d0
          else
             if (n .eq. 1) then
                res = dsqrt(a11)*isr_tilde(x,beta)
             else if (n .eq. 2) then
                res = dsqrt(a11)*isr_tilde(x,beta)
             else if (n .eq. 3) then
                res = Iapq(x,a10/dsqrt(a11),p2,q2,beta)
             else if (n .eq. 4) then
                res = Iapq(x,dsqrt(a00),p4,q4,beta)
             else
                res = 0d0
             endif
          endif
      endif
      eepdf_tilde_calc = res
      end


program main

  use mo_kind,       only: int8=>i1, int16=>i2, int32=>i4, int64=>i8, real32=>sp, real64=>dp
#if defined (__pgiFortran__) || defined (__NAGf90Fortran__) || defined (__NAG__)
  use mo_kind,       only: real128=>dp
#else
  use mo_kind,       only: real128=>qp
#endif
  use mo_testing,    only: assert, initialize_tests, report_tests
  use mo_ansi_colors, only: color, c_red, c_green

  use mo_functional

  use mo_filter_functions
  use mo_fold_functions
  use mo_iterfold_functions
  use mo_map_functions
  use mo_unfold_functions

  implicit none

  logical, dimension(:), allocatable :: tests
  logical :: test_failed
  integer :: n, ntests
  integer, parameter :: stdout = 6

  real(kind=real32),dimension(1000) :: x

  complex(kind=real32),dimension(:),allocatable :: c4,c4_res
  complex(kind=real64),dimension(:),allocatable :: c8,c8_res
  complex(kind=real128),dimension(:),allocatable :: c16,c16_res
  complex(kind=real64) :: c8_start
  complex(kind=real128) :: c16_start
  complex(kind=real32),dimension(:),allocatable :: c_r4
  complex(kind=real64),dimension(:),allocatable :: c_r8
  complex(kind=real128),dimension(:),allocatable :: c_r16

  logical  :: isgood

  isgood = .true.

  n = 1
  ntests = 19
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(arange(1_int8,3_int8) == [1_int8,2_int8,3_int8]),&
       'arange, int8')
  n = n + 1

  tests(n) = assert(all(arange(1_int16,3_int16) == [1_int16,2_int16,3_int16]),&
       'arange, int16')
  n = n + 1

  tests(n) = assert(all(arange(1_int32,3_int32) == [1_int32,2_int32,3_int32]),&
       'arange, int32')
  n = n + 1

  tests(n) = assert(all(arange(1_int64,3_int64) == [1_int64,2_int64,3_int64]),&
       'arange, int64')
  n = n + 1

  tests(n) = assert(all(arange(1._real32,3._real32) == [1._real32,2._real32,3._real32]),&
       'arange, real32')
  n = n + 1

  tests(n) = assert(all(arange(1._real64,3._real64) == [1._real64,2._real64,3._real64]),&
       'arange, real32')
  n = n + 1

  tests(n) = assert(all(arange(1._real128,3._real128) == [1._real128,2._real128,3._real128]),&
       'arange, real128')
  n = n + 1

  tests(n) = assert(all(arange(1._real128,3._real128) == [1._real128,2._real128,3._real128]),&
       'arange, real128')
  n = n + 1

  tests(n) = assert(all(arange(cmplx(1._real32,0._real32, kind=real32),&
       cmplx(3._real32,0._real32, kind=real32),&
       cmplx(1._real32,0._real32, kind=real32))&
       == [cmplx(1._real32,0._real32, kind=real32),&
       cmplx(2._real32,0._real32, kind=real32),&
       cmplx(3._real32,0._real32, kind=real32)]),&
       'arange, complex real32')
  n = n + 1

  tests(n) = assert(all(arange(cmplx(1._real64,0._real64, kind=real64),&
       cmplx(3._real64,0._real64, kind=real64),&
       cmplx(1._real64,0._real64, kind=real64))&
       == [cmplx(1._real64,0._real64, kind=real64),&
       cmplx(2._real64,0._real64, kind=real64),&
       cmplx(3._real64,0._real64, kind=real64)]),&
       'arange, complex real64')
  n = n + 1

  tests(n) = assert(all(arange(cmplx(1._real128,0._real128, kind=real128),&
       cmplx(3._real128,0._real128, kind=real128),&
       cmplx(1._real128,0._real128, kind=real128))&
       == [cmplx(1._real128,0._real128, kind=real128),&
       cmplx(2._real128,0._real128, kind=real128),&
       cmplx(3._real128,0._real128, kind=real128)]),&
       'arange, complex real128')
  n = n + 1

  tests(n) = assert(all(arange(cmplx(1,1, kind=real32),cmplx(3,5, kind=real32),cmplx(1,2, kind=real32))&
       == [cmplx(1,1, kind=real32),cmplx(2,3, kind=real32),cmplx(3,5, kind=real32)]),&
       'arange, incrementing both parts of complex numbers')
  n = n + 1

  tests(n) = assert(all(arange(1,10) == arange(1,10,1)),&
       'arange increment equals 1 when ommited')
  n = n + 1

  tests(n) = assert(all(arange(1.,10.) == real(arange(1,10,1))),&
       'integer and real arange variants produce same values')
  n = n + 1

  tests(n) = assert(all(arange(0.,2.4,0.8) == [0.,0.8,1.6,2.4]),&
       'custom increment value')
  n = n + 1

  tests(n) = assert(all(arange(3,-1,-1) == [3,2,1,0,-1]),&
       'negative increment value')
  n = n + 1

  tests(n) = assert(size(arange(1.0,1.4,0.1)) == 5,&
       'real32-typed arange returns array of expected size')
  n = n + 1

  tests(n) = assert(size(arange(1.0_real64,1.4_real64,0.1_real64)) == 5,&
       'real64-typed arange returns array of expected size')
  n = n + 1

  tests(n) = assert(size(arange(1.0_real128,1.4_real128,0.1_real128)) == 5,&
       'real128-typed arange returns array of expected size')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 11
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(complement([1_int8,2_int8],[2_int8]) == [1]),&
       'complement, int8')
  n = n + 1

  tests(n) = assert(all(complement([1_int16,2_int16],[2_int16]) == [1]),&
       'complement, int16')
  n = n + 1

  tests(n) = assert(all(complement([1_int32,2_int32],[2_int32]) == [1]),&
       'complement, int32')
  n = n + 1

  tests(n) = assert(all(complement([1_int64,2_int64],[2_int64]) == [1]),&
       'complement, int64')
  n = n + 1

  tests(n) = assert(all(complement([1._real32,2._real32],[2._real32]) == [1]),&
       'complement, real32')
  n = n + 1

  tests(n) = assert(all(complement([1._real64,2._real64],[2._real64]) == [1]),&
       'complement, real64')
  n = n + 1

  tests(n) = assert(all(complement([1._real128,2._real128],[2._real128]) == [1]),&
       'complement, real128')
  n = n + 1

  tests(n) = assert(all(complement([cmplx(1._real32,0._real32, kind=real32),&
       cmplx(2._real32,0._real32, kind=real32)],&
       [cmplx(2._real32,0._real32, kind=real32)])&
       == [cmplx(1._real32,0._real32, kind=real32)]),'complement, complex real32')
  n = n + 1

  tests(n) = assert(all(complement([cmplx(1._real64,0._real64, kind=real64),&
       cmplx(2._real64,0._real64, kind=real64)],&
       [cmplx(2._real64,0._real64, kind=real64)])&
       == [cmplx(1._real64,0._real64, kind=real64)]),'complement, complex real64')
  n = n + 1

  tests(n) = assert(all(complement([cmplx(1._real64,0._real64, kind=real64),&
       cmplx(2._real64,0._real64, kind=real64)],&
       [cmplx(2._real64,0._real64, kind=real64)])&
       == [cmplx(1._real64,0._real64, kind=real64)]),'complement, complex real64')
  n = n + 1

  tests(n) = assert(all(complement([1,2],[2]) == ([1,2].complement.[2])),&
       'complement operator, x.complement.y')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 11
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(filter(gt3lt5_i1,[3_int8,4_int8,5_int8]) == [4]),&
       'filter, int8')
  n = n + 1

  tests(n) = assert(all(filter(gt3lt5_i2,[3_int16,4_int16,5_int16]) == [4]),&
       'filter, int16')
  n = n + 1

  tests(n) = assert(all(filter(gt3lt5_i4,[3,4,5]) == [4]),&
       'filter, int32')
  n = n + 1

  tests(n) = assert(all(filter(gt3lt5_i8,[3_int64,4_int64,5_int64]) == [4]),&
       'filter, int64')
  n = n + 1

  tests(n) = assert(all(filter(gt3lt5_r4,[3.,4.,5.]) == [4]),&
       'filter, real32')
  n = n + 1

  tests(n) = assert(all(filter(gt3lt5_r8,[3._real64,4._real64,5._real64]) == [4]),&
       'filter, real64')
  n = n + 1

  tests(n) = assert(all(filter(gt3lt5_r16,[3._real128,4._real128,5._real128]) == [4]),&
       'filter, real128')
  n = n + 1

  tests(n) = assert(all(filter(gt3lt5_c4,&
       [cmplx(3.,0., kind=real32),cmplx(4.,0., kind=real32),cmplx(5.,0., kind=real32)]) == [cmplx(4.,0., kind=real32)]),&
       'filter, complex real32')
  n = n + 1

  ! Need to assign to a variable first because cmplx() by default
  ! returns single-precision complex number which breaks the generic
  ! interface
  c8 = [cmplx(3._real64,0._real64, kind=real64),&
       cmplx(4._real64,0._real64, kind=real64),&
       cmplx(5._real64,0._real64, kind=real64)]
  c16 = [cmplx(3._real128,0._real128, kind=real128),&
       cmplx(4._real128,0._real128, kind=real128),&
       cmplx(5._real128,0._real128, kind=real128)]

  tests(n) = assert(all(filter(gt3lt5_c8,c8) == [cmplx(4.,0., kind=real64)]),&
       'filter, complex real64')
  n = n + 1

  tests(n) = assert(all(filter(gt3lt5_c16,c16) == [cmplx(4.,0., kind=real128)]),&
       'filter, complex real128')
  n = n + 1

  tests(n) = assert(size(filter(gt3lt5_i4,[1,2,3,5,6])) == 0,&
       'filter returns empty array')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 10
  call initialize_tests(tests,ntests)

  tests(n) = assert(foldl(sum_i1,0_int8,[1_int8,2_int8,3_int8,4_int8,5_int8]) == 15,&
       'foldl, int8')
  n = n + 1

  tests(n) = assert(foldl(sum_i2,0_int16,[1_int16,2_int16,3_int16,4_int16,5_int16]) == 15,&
       'foldl, int16')
  n = n + 1

  tests(n) = assert(foldl(sum_i4,0_int32,[1_int32,2_int32,3_int32,4_int32,5_int32]) == 15,&
       'foldl, int32')
  n = n + 1

  tests(n) = assert(foldl(sum_i8,0_int64,[1_int64,2_int64,3_int64,4_int64,5_int64]) == 15,&
       'foldl, int64')
  n = n + 1

  tests(n) = assert(foldl(sum_r4,0._real32,[1._real32,2._real32,3._real32,4._real32,5._real32]) == 15,&
       'foldl, real32')
  n = n + 1

  tests(n) = assert(foldl(sum_r8,0._real64,[1._real64,2._real64,3._real64,4._real64,5._real64]) == 15,&
       'foldl, real64')
  n = n + 1

  tests(n) = assert(foldl(sum_r16,0._real128,[1._real128,2._real128,3._real128,4._real128,5._real128]) == 15,&
       'foldl, real128')
  n = n + 1

  c4 = arange(cmplx(1,0, kind=real32),cmplx(5,0, kind=real32))
  c8 = c4
  c16 = c4

  c8_start = cmplx(0,0)
  c16_start = c8_start

  tests(n) = assert(foldl(sum_c4,cmplx(0.,0., kind=real32),c4) == cmplx(15,0, kind=real32),&
       'foldl, complex real32')
  n = n + 1

  tests(n) = assert(foldl(sum_c8,c8_start,c8) == cmplx(15._real64,0._real64, kind=real64),&
       'foldl, complex real64')
  n = n + 1

  tests(n) = assert(foldl(sum_c16,c16_start,c16) == cmplx(15._real128,0._real128, kind=real128),&
       'foldl, complex real128')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 10
  call initialize_tests(tests,ntests)

  tests(n) = assert(foldr(sum_i1,0_int8,[1_int8,2_int8,3_int8,4_int8,5_int8]) == 15,&
       'foldr, int8')
  n = n + 1

  tests(n) = assert(foldr(sum_i2,0_int16,[1_int16,2_int16,3_int16,4_int16,5_int16]) == 15,&
       'foldr, int16')
  n = n + 1

  tests(n) = assert(foldr(sum_i4,0_int32,[1_int32,2_int32,3_int32,4_int32,5_int32]) == 15,&
       'foldr, int32')
  n = n + 1

  tests(n) = assert(foldr(sum_i8,0_int64,[1_int64,2_int64,3_int64,4_int64,5_int64]) == 15,&
       'foldr, int64')
  n = n + 1

  tests(n) = assert(foldr(sum_r4,0._real32,[1._real32,2._real32,3._real32,4._real32,5._real32]) == 15,&
       'foldr, real32')
  n = n + 1

  tests(n) = assert(foldr(sum_r8,0._real64,[1._real64,2._real64,3._real64,4._real64,5._real64]) == 15,&
       'foldr, real64')
  n = n + 1

  tests(n) = assert(foldr(sum_r16,0._real128,[1._real128,2._real128,3._real128,4._real128,5._real128]) == 15,&
       'foldr, real128')
  n = n + 1

  c4 = arange(cmplx(1,0, kind=real32),cmplx(5,0, kind=real32))
  c8 = c4
  c16 = c4

  c8_start = cmplx(0,0, kind=real32)
  c16_start = c8_start

  tests(n) = assert(foldr(sum_c4,cmplx(0.,0., kind=real32),c4) == cmplx(15,0, kind=real32),&
       'foldr, complex real32')
  n = n + 1

  tests(n) = assert(foldr(sum_c8,c8_start,c8) == cmplx(15._real64,0._real64, kind=real64),&
       'foldr, complex real64')
  n = n + 1

  tests(n) = assert(foldr(sum_c16,c16_start,c16) == cmplx(15._real128,0._real128, kind=real128),&
       'foldr, complex real128')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 10
  call initialize_tests(tests,ntests)

  tests(n) = assert(foldt(sum_i1,0_int8,[1_int8,2_int8,3_int8,4_int8,5_int8]) == 15,&
       'foldt, int8')
  n = n + 1

  tests(n) = assert(foldt(sum_i2,0_int16,[1_int16,2_int16,3_int16,4_int16,5_int16]) == 15,&
       'foldt, int16')
  n = n + 1

  tests(n) = assert(foldt(sum_i4,0_int32,[1_int32,2_int32,3_int32,4_int32,5_int32]) == 15,&
       'foldt, int32')
  n = n + 1

  tests(n) = assert(foldt(sum_i8,0_int64,[1_int64,2_int64,3_int64,4_int64,5_int64]) == 15,&
       'foldt, int64')
  n = n + 1

  tests(n) = assert(foldt(sum_r4,0._real32,[1._real32,2._real32,3._real32,4._real32,5._real32]) == 15,&
       'foldt, real32')
  n = n + 1

  tests(n) = assert(foldt(sum_r8,0._real64,[1._real64,2._real64,3._real64,4._real64,5._real64]) == 15,&
       'foldt, real64')
  n = n + 1

  tests(n) = assert(foldt(sum_r16,0._real128,[1._real128,2._real128,3._real128,4._real128,5._real128]) == 15,&
       'foldt, real128')
  n = n + 1

  c4 = arange(cmplx(1,0, kind=real32),cmplx(5,0, kind=real32))
  c8 = c4
  c16 = c4

  c8_start = cmplx(0,0, kind=real32)
  c16_start = c8_start

  tests(n) = assert(foldt(sum_c4,cmplx(0.,0., kind=real32),c4) == cmplx(15,0, kind=real32),&
       'foldt, complex real32')
  n = n + 1

  tests(n) = assert(foldt(sum_c8,c8_start,c8) == cmplx(15._real64,0._real64, kind=real64),&
       'foldt, complex real64')
  n = n + 1

  tests(n) = assert(foldt(sum_c16,c16_start,c16) == cmplx(15._real128,0._real128, kind=real128),&
       'foldt, complex real128')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  c_r4 = [(1,2),(2,4)]
  c_r8 = c_r4
  c_r16 = c_r4

  n = 1
  ntests = 11
  call initialize_tests(tests,ntests)

  tests(n) = assert(head([1_int8,2_int8]) == 1_int8,'head, int8')
  n = n + 1

  tests(n) = assert(head([1_int16,2_int16]) == 1_int16,'head, int16')
  n = n + 1

  tests(n) = assert(head([1_int32,2_int32]) == 1_int32,'head, int32')
  n = n + 1

  tests(n) = assert(head([1_int64,2_int64]) == 1_int64,'head, int64')
  n = n + 1

  tests(n) = assert(head([1._real32,2._real32]) == 1._real32,'head, real32')
  n = n + 1

  tests(n) = assert(head([1._real64,2._real64]) == 1._real64,'head, real64')
  n = n + 1

  tests(n) = assert(head([1._real128,2._real128]) == 1._real128,'head, real128')
  n = n + 1

  tests(n) = assert(head(c_r4) == c_r4(1),'head, complex real32')
  n = n + 1

  tests(n) = assert(head(c_r8) == c_r8(1),'head, complex real64')
  n = n + 1

  tests(n) = assert(head(c_r16) == c_r16(1),'head, complex real128')
  n = n + 1

#ifdef __pgiFortran__
  tests(n) = assert(.true.,'skip head operator, .head.x')
#else
  tests(n) = assert(head([1,2]) == .head.[1,2],'head operator, .head.x')
#endif
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  c_r4 = [(1,2),(2,4)]
  c_r8 = c_r4
  c_r16 = c_r4

  n = 1
  ntests = 13
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(init([1_int8,2_int8]) == [1_int8]),'init, int8')
  n = n + 1

  tests(n) = assert(all(init([1_int16,2_int16]) == [1_int16]),'init, int16')
  n = n + 1

  tests(n) = assert(all(init([1_int32,2_int32]) == [1_int32]),'init, int32')
  n = n + 1

  tests(n) = assert(all(init([1_int64,2_int64]) == [1_int64]),'init, int64')
  n = n + 1

  tests(n) = assert(all(init([1._real32,2._real32]) == [1._real32]),'init, real32')
  n = n + 1

  tests(n) = assert(all(init([1._real64,2._real64]) == [1._real64]),'init, real64')
  n = n + 1

  tests(n) = assert(all(init([1._real128,2._real128]) == [1._real128]),'init, real128')
  n = n + 1

  tests(n) = assert(all(init(c_r4) == [c_r4(1)]),'init, complex real32')
  n = n + 1

  tests(n) = assert(all(init(c_r8) == [c_r8(1)]),'init, complex real64')
  n = n + 1

  tests(n) = assert(all(init(c_r16) == [c_r16(1)]),'init, complex real128')
  n = n + 1

  tests(n) = assert(size(init([1])) == 0,'size(init([1])) == 0')
  n = n + 1

  tests(n) = assert(size(init(init([1]))) == 0,'size(init(init([1]))) == 0')
  n = n + 1

#ifdef __pgiFortran__
  tests(n) = assert(.true.,'skip init operator, .init.x')
#else
  tests(n) = assert(all(init([1,2]) == .init.[1,2]),'init operator, .init.x')
#endif
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 12
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(insert(2_int8,2,[1_int8,3_int8]) == [1,2,3]),&
       'insert, int8')
  n = n + 1

  tests(n) = assert(all(insert(2_int16,2,[1_int16,3_int16]) == [1,2,3]),&
       'insert, int16')
  n = n + 1

  tests(n) = assert(all(insert(2_int32,2,[1_int32,3_int32]) == [1,2,3]),&
       'insert, int32')
  n = n + 1

  tests(n) = assert(all(insert(2_int64,2,[1_int64,3_int64]) == [1,2,3]),&
       'insert, int64')
  n = n + 1

  tests(n) = assert(all(insert(2._real32,2,[1._real32,3._real32]) == [1,2,3]),&
       'insert, real32')
  n = n + 1

  tests(n) = assert(all(insert(2._real64,2,[1._real64,3._real64]) == [1,2,3]),&
       'insert, real64')
  n = n + 1

  tests(n) = assert(all(insert(2._real128,2,[1._real128,3._real128]) == [1,2,3]),&
       'insert, real128')
  n = n + 1

  tests(n) = assert(all(insert(cmplx(2._real32,0._real32, kind=real32),2,&
       [cmplx(1._real32,0._real32, kind=real32),cmplx(3._real32,0._real32, kind=real32)])&
       == arange(cmplx(1._real32,0._real32, kind=real32),cmplx(3._real32,0._real32, kind=real32))),&
       'insert, real32')
  n = n + 1

  tests(n) = assert(all(insert(cmplx(2._real64,0._real64, kind=real64),2,&
       [cmplx(1._real64,0._real64, kind=real64),cmplx(3._real64,0._real64, kind=real64)])&
       == arange(cmplx(1._real64,0._real64, kind=real64),cmplx(3._real64,0._real64, kind=real64))),&
       'insert, real64')
  n = n + 1

  tests(n) = assert(all(insert(cmplx(2._real128,0._real128, kind=real128),2,&
       [cmplx(1._real128,0._real128, kind=real128),cmplx(3._real128,0._real128, kind=real128)])&
       == arange(cmplx(1._real128,0._real128, kind=real128),cmplx(3._real128,0._real128, kind=real128))),&
       'insert, real128')
  n = n + 1

  tests(n) = assert(all(insert(1,1,arange(1,0)) == [1]),&
       'insert into empty array')
  n = n + 1

  tests(n) = assert(all(insert(2,2,[1]) == [1,2]),&
       'insert out of bounds')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 11
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(intersection([1_int8,2_int8],[2_int8,3_int8]) == [2]),&
       'intersection, int8')
  n = n + 1

  tests(n) = assert(all(intersection([1_int16,2_int16],[2_int16,3_int16]) == [2]),&
       'intersection, int16')
  n = n + 1

  tests(n) = assert(all(intersection([1_int32,2_int32],[2_int32,3_int32]) == [2]),&
       'intersection, int32')
  n = n + 1

  tests(n) = assert(all(intersection([1_int64,2_int64],[2_int64,3_int64]) == [2]),&
       'intersection, int64')
  n = n + 1

  tests(n) = assert(all(intersection([1._real32,2._real32],[2._real32,3._real32]) == [2]),&
       'intersection, real32')
  n = n + 1

  tests(n) = assert(all(intersection([1._real64,2._real64],[2._real64,3._real64]) == [2]),&
       'intersection, real64')
  n = n + 1

  tests(n) = assert(all(intersection([1._real128,2._real128],[2._real128,3._real128]) == [2]),&
       'intersection, real128')
  n = n + 1

  tests(n) = assert(all(intersection([cmplx(1._real32,0._real32, kind=real32),cmplx(2._real32,0._real32, kind=real32)],&
       [cmplx(2._real32,0._real32, kind=real32),cmplx(3._real32,0._real32, kind=real32)])&
       == [cmplx(2._real32,0._real32, kind=real32)]),&
       'intersection, complex real32')
  n = n + 1

  tests(n) = assert(all(intersection([cmplx(1._real64,0._real64, kind=real64),cmplx(2._real64,0._real64, kind=real64)],&
       [cmplx(2._real64,0._real64, kind=real64),cmplx(3._real64,0._real64, kind=real64)])&
       == [cmplx(2._real64,0._real64, kind=real64)]),&
       'intersection, complex real64')
  n = n + 1

  tests(n) = assert(all(intersection([cmplx(1._real128,0._real128, kind=real128),cmplx(2._real128,0._real128, kind=real128)],&
       [cmplx(2._real128,0._real128, kind=real128),cmplx(3._real128,0._real128, kind=real128)])&
       == [cmplx(2._real128,0._real128, kind=real128)]),&
       'intersection, complex real128')
  n = n + 1

  tests(n) = assert(all(intersection([1,2],[2,3]) == ([1,2].intersection.[2,3])),&
       'intersection operator, x.intersection.y')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 10
  call initialize_tests(tests,ntests)

  tests(n) = assert(iterfold(isum_i1,0_int8,[1_int8,2_int8,3_int8,4_int8,5_int8]) == 15,&
       'iterfold, int8')
  n = n + 1

  tests(n) = assert(iterfold(isum_i2,0_int16,[1_int16,2_int16,3_int16,4_int16,5_int16]) == 15,&
       'iterfold, int16')
  n = n + 1

  tests(n) = assert(iterfold(isum_i4,0_int32,[1_int32,2_int32,3_int32,4_int32,5_int32]) == 15,&
       'iterfold, int32')
  n = n + 1

  tests(n) = assert(iterfold(isum_i8,0_int64,[1_int64,2_int64,3_int64,4_int64,5_int64]) == 15,&
       'iterfold, int64')
  n = n + 1

  tests(n) = assert(iterfold(isum_r4,0._real32,[1._real32,2._real32,3._real32,4._real32,5._real32]) == 15,&
       'iterfold, real32')
  n = n + 1

  tests(n) = assert(iterfold(isum_r8,0._real64,[1._real64,2._real64,3._real64,4._real64,5._real64]) == 15,&
       'iterfold, real64')
  n = n + 1

  tests(n) = assert(iterfold(isum_r16,0._real128,[1._real128,2._real128,3._real128,4._real128,5._real128]) == 15,&
       'iterfold, real128')
  n = n + 1

  c4 = arange(cmplx(1,0, kind=real32),cmplx(5,0, kind=real32))
  c8 = c4
  c16 = c4

  c8_start = cmplx(0,0, kind=real32)
  c16_start = c8_start

  tests(n) = assert(iterfold(isum_c4,cmplx(0.,0., kind=real32),c4) == cmplx(15,0, kind=real32),&
       'iterfold, complex real32')
  n = n + 1

  tests(n) = assert(iterfold(isum_c8,c8_start,c8) == cmplx(15._real64,0._real64, kind=real64),&
       'iterfold, complex real64')
  n = n + 1

  tests(n) = assert(iterfold(isum_c16,c16_start,c16) == cmplx(15._real128,0._real128, kind=real128),&
       'iterfold, complex real128')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  c_r4 = [(1,2),(2,4)]
  c_r8 = c_r4
  c_r16 = c_r4

  n = 1
  ntests = 11
  call initialize_tests(tests,ntests)

  tests(n) = assert(last([1_int8,2_int8]) == 2_int8,'last, int8')
  n = n + 1

  tests(n) = assert(last([1_int16,2_int16]) == 2_int16,'last, int16')
  n = n + 1

  tests(n) = assert(last([1_int32,2_int32]) == 2_int32,'last, int32')
  n = n + 1

  tests(n) = assert(last([1_int64,2_int64]) == 2_int64,'last, int64')
  n = n + 1

  tests(n) = assert(last([1._real32,2._real32]) == 2._real32,'last, real32')
  n = n + 1

  tests(n) = assert(last([1._real64,2._real64]) == 2._real64,'last, real64')
  n = n + 1

  tests(n) = assert(last([1._real128,2._real128]) == 2._real128,'last, real128')
  n = n + 1

  tests(n) = assert(last(c_r4) == c_r4(2),'last, complex real32')
  n = n + 1

  tests(n) = assert(last(c_r8) == c_r8(2),'last, complex real64')
  n = n + 1

  tests(n) = assert(last(c_r16) == c_r16(2),'last, complex real128')
  n = n + 1

#ifdef __pgiFortran__
  tests(n) = assert(.true.,'skip last operator, .last.x')
#else
  tests(n) = assert(last([1,2]) == .last.[1,2],'last operator, .last.x')
#endif
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 11
  call initialize_tests(tests,ntests)

  tests(n) = assert(limit(2_int8,1_int8,3_int8) == 2_int8,&
       'limit, int8')
  n = n + 1

  tests(n) = assert(limit(2_int16,1_int16,3_int16) == 2_int16,&
       'limit, int16')
  n = n + 1

  tests(n) = assert(limit(2_int32,1_int32,3_int32) == 2_int32,&
       'limit, int32')
  n = n + 1

  tests(n) = assert(limit(2_int64,1_int64,3_int64) == 2_int64,&
       'limit, int64')
  n = n + 1

  tests(n) = assert(limit(2._real32,1._real32,3._real32) == 2._real32,&
       'limit, real32')
  n = n + 1

  tests(n) = assert(limit(2._real64,1._real64,3._real64) == 2._real64,&
       'limit, real32')
  n = n + 1

  tests(n) = assert(limit(2._real128,1._real128,3._real128) == 2._real128,&
       'limit, real128')
  n = n + 1

  tests(n) = assert(limit(cmplx(-0.5,1.5, kind=real32),cmplx(0,0, kind=real32),cmplx(1,1, kind=real32)) == cmplx(0,1, kind=real32),&
       'limit, complex real32')
  n = n + 1

  tests(n) = assert(limit(cmplx(-0.5_real64,1.5_real64, kind=real64),cmplx(0._real64,0._real64, kind=real64),&
       cmplx(1._real64,1._real64, kind=real64)) == cmplx(0._real64,1._real64, kind=real64),&
       'limit, complex real64')
  n = n + 1

  tests(n) = assert(limit(cmplx(-0.5_real128,1.5_real128, kind=real128),cmplx(0._real128,0._real128, kind=real128),&
       cmplx(1._real128,1._real128, kind=real128)) == cmplx(0._real128,1._real128, kind=real128),&
       'limit, complex real128')
  n = n + 1

  tests(n) = assert(all(limit(arange(1,3),2,2) == [2,2,2]),&
       'limit works on arrays')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 10
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(map(xpowx_i1,[1_int8,2_int8,3_int8])&
       == [1_int8,4_int8,27_int8]),'map, int8')
  n = n + 1

  tests(n) = assert(all(map(xpowx_i2,[1_int16,2_int16,3_int16])&
       == [1_int16,4_int16,27_int16]),'map, int16')
  n = n + 1

  tests(n) = assert(all(map(xpowx_i4,[1_int32,2_int32,3_int32])&
       == [1_int32,4_int32,27_int32]),'map, int32')
  n = n + 1

  tests(n) = assert(all(map(xpowx_i8,[1_int64,2_int64,3_int64])&
       == [1_int64,4_int64,27_int64]),'map, int64')
  n = n + 1

  tests(n) = assert(all(map(xpowx_r4,[1._real32,2._real32,3._real32])&
       == [1._real32,4._real32,27._real32]),'map, real32')
  n = n + 1

  tests(n) = assert(all(map(xpowx_r8,[1._real64,2._real64,3._real64])&
       == [1._real64,4._real64,27._real64]),'map, real64')
  n = n + 1

  tests(n) = assert(all(map(xpowx_r16,[1._real128,2._real128,3._real128])&
       == [1._real128,4._real128,27._real128]),'map, real128')
  n = n + 1

  c4 = [cmplx(1.,0., kind=real32),cmplx(2.,0., kind=real32),cmplx(3.,0., kind=real32)]
  c8 = [cmplx(1._real64,0._real64, kind=real64),&
       cmplx(2._real64,0._real64, kind=real64),&
       cmplx(3._real64,0._real64, kind=real64)]
  c16 = [cmplx(1._real128,0._real128, kind=real128),&
       cmplx(2._real128,0._real128, kind=real128),&
       cmplx(3._real128,0._real128, kind=real128)]

  tests(n) = assert(all(map(xpowx_c4,c4) == c4**c4),'map, complex real32')
  n = n + 1

  tests(n) = assert(all(map(xpowx_c8,c8) == c8**c8),'map, complex real64')
  n = n + 1

  tests(n) = assert(all(map(xpowx_c16,c16) == c16**c16),'map, complex real128')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 11
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(reverse(arange(1_int8,3_int8)) == [3,2,1]),&
       'reverse, int8')
  n = n + 1

  tests(n) = assert(all(reverse(arange(1_int16,3_int16)) == [3,2,1]),&
       'reverse, int16')
  n = n + 1

  tests(n) = assert(all(reverse(arange(1_int32,3_int32)) == [3,2,1]),&
       'reverse, int32')
  n = n + 1

  tests(n) = assert(all(reverse(arange(1_int64,3_int64)) == [3,2,1]),&
       'reverse, int64')
  n = n + 1

  tests(n) = assert(all(reverse(arange(1._real32,3._real32)) == [3,2,1]),&
       'reverse, real32')
  n = n + 1

  tests(n) = assert(all(reverse(arange(1._real64,3._real64)) == [3,2,1]),&
       'reverse, real64')
  n = n + 1

  tests(n) = assert(all(reverse(arange(1._real128,3._real128)) == [3,2,1]),&
       'reverse, real128')
  n = n + 1

  tests(n) = assert(all(reverse(arange(cmplx(1._real32,0._real32, kind=real32),&
       cmplx(3._real32,0._real32, kind=real32)))&
       == arange(cmplx(3._real32,0._real32, kind=real32),&
       cmplx(1._real32,0._real32, kind=real32),&
       cmplx(-1._real32,0._real32, kind=real32))),&
       'reverse, complex real32')
  n = n + 1

  tests(n) = assert(all(reverse(arange(cmplx(1._real64,0._real64, kind=real64),&
       cmplx(3._real64,0._real64, kind=real64)))&
       == arange(cmplx(3._real64,0._real64, kind=real64),&
       cmplx(1._real64,0._real64, kind=real64),&
       cmplx(-1._real64,0._real64, kind=real64))),&
       'reverse, complex real64')
  n = n + 1

  tests(n) = assert(all(reverse(arange(cmplx(1._real128,0._real128, kind=real64),&
       cmplx(3._real128,0._real128, kind=real64)))&
       == arange(cmplx(3._real128,0._real128, kind=real64),&
       cmplx(1._real128,0._real128, kind=real64),&
       cmplx(-1._real128,0._real128, kind=real64))),&
       'reverse, complex real128')
  n = n + 1

#ifdef __pgiFortran__
  tests(n) = assert(.true.,'skip reverse operator, .reverse.x')
#else
  tests(n) = assert(all(reverse([1,2,3]) == .reverse.[1,2,3]),&
       'reverse operator, .reverse.x')
#endif
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 12
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(set([1_int8,2_int8,2_int8,3_int8]) == [1,2,3]),&
       'set, int8')
  n = n + 1

  tests(n) = assert(all(set([1_int16,2_int16,2_int16,3_int16]) == [1,2,3]),&
       'set, int16')
  n = n + 1

  tests(n) = assert(all(set([1_int32,2_int32,2_int32,3_int32]) == [1,2,3]),&
       'set, int32')
  n = n + 1

  tests(n) = assert(all(set([1_int64,2_int64,2_int64,3_int64]) == [1,2,3]),&
       'set, int64')
  n = n + 1

  tests(n) = assert(all(set([1._real32,2._real32,2._real32,3._real32]) == [1,2,3]),&
       'set, real32')
  n = n + 1

  tests(n) = assert(all(set([1._real64,2._real64,2._real64,3._real64]) == [1,2,3]),&
       'set, real64')
  n = n + 1

  tests(n) = assert(all(set([1._real128,2._real128,2._real128,3._real128]) == [1,2,3]),&
       'set, real128')
  n = n + 1

  c4 = [cmplx(1,0, kind=real32),cmplx(2,0, kind=real32),cmplx(2,0, kind=real32),cmplx(3,0, kind=real32)]
  c4_res = [cmplx(1,0, kind=real32),cmplx(2,0, kind=real32),cmplx(3,0, kind=real32)]
  tests(n) = assert(all(set(c4) == c4_res),'set, complex real32')
  n = n + 1

  c8 = c4
  c8_res = c4_res
  tests(n) = assert(all(set(c8) == c8_res),'set, complex real64')
  n = n + 1

  c16 = c4
  c16_res = c4_res
  tests(n) = assert(all(set(c16) == c16_res),'set, complex real128')
  n = n + 1

  tests(n) = assert(all(set(arange(1,0)) == arange(1,0)),&
       'set of empty array is an empty array')
  n = n + 1

#ifdef __pgiFortran__
  tests(n) = assert(.true.,'skip set operator, .set.x')
#else
  tests(n) = assert(all(set([1,2,2,3]) == .set.[1,2,2,3]),&
       'set operator, .set.x')
#endif
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 12
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(sort([3_int8,2_int8,1_int8]) == [1,2,3]),&
       'sort, int8')
  n = n + 1

  tests(n) = assert(all(sort([3_int16,2_int16,1_int16]) == [1,2,3]),&
       'sort, int16')
  n = n + 1

  tests(n) = assert(all(sort([3_int32,2_int32,1_int32]) == [1,2,3]),&
       'sort, int32')
  n = n + 1

  tests(n) = assert(all(sort([3_int64,2_int64,1_int64]) == [1,2,3]),&
       'sort, int8')
  n = n + 1

  tests(n) = assert(all(sort([3._real32,2._real32,1._real32]) == [1,2,3]),&
       'sort, real32')
  n = n + 1

  tests(n) = assert(all(sort([3._real64,2._real64,1._real64]) == [1,2,3]),&
       'sort, real64')
  n = n + 1

  tests(n) = assert(all(sort([3._real128,2._real128,1._real128]) == [1,2,3]),&
       'sort, real128')
  n = n + 1

  tests(n) = assert(all(sort(arange(cmplx(3._real32,0._real32, kind=real32),&
       cmplx(1._real32,0._real32, kind=real32),&
       cmplx(-1._real32,0._real32, kind=real32)))&
       == arange(cmplx(1._real32,0._real32, kind=real32),&
       cmplx(3._real32,0._real32, kind=real32))),&
       'sort, complex real32')
  n = n + 1

  tests(n) = assert(all(sort(arange(cmplx(3._real64,0._real64, kind=real64),&
       cmplx(1._real64,0._real64, kind=real64),&
       cmplx(-1._real64,0._real64, kind=real64)))&
       == arange(cmplx(1._real64,0._real64, kind=real64),&
       cmplx(3._real64,0._real64, kind=real64))),&
       'sort, complex real64')
  n = n + 1

  tests(n) = assert(all(sort(arange(cmplx(3._real128,0._real128, kind=real128),&
       cmplx(1._real128,0._real128, kind=real128),&
       cmplx(-1._real128,0._real128, kind=real128)))&
       == arange(cmplx(1._real128,0._real128, kind=real128),&
       cmplx(3._real128,0._real128, kind=real128))),&
       'sort, complex real128')
  n = n + 1

  call random_number(x)
  tests(n) = assert(all(tail(sort(x)) >= init(sort(x))),&
       'all(tail(sort(x)) >= init(sort(x))')
  n = n + 1

#ifdef __pgiFortran__
  tests(n) = assert(.true.,'skip sort operator, .sort.x')
#else
  tests(n) = assert(all(sort(x) == .sort.x),&
       'sort operator, .sort.x')
#endif
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 16
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(split(arange(1_int8,10_int8),1) == arange(1,5)),&
       'split(x,1), int8')
  n = n + 1

  tests(n) = assert(all(split(arange(1_int8,10_int8),2) == arange(6,10)),&
       'split(x,2), int8')
  n = n + 1

  tests(n) = assert(all(split(arange(1_int16,10_int16),1) == arange(1,5)),&
       'split(x,1), int16')
  n = n + 1

  tests(n) = assert(all(split(arange(1_int16,10_int16),2) == arange(6,10)),&
       'split(x,2), int16')
  n = n + 1

  tests(n) = assert(all(split(arange(1_int32,10_int32),1) == arange(1,5)),&
       'split(x,1), int32')
  n = n + 1

  tests(n) = assert(all(split(arange(1_int32,10_int32),2) == arange(6,10)),&
       'split(x,2), int32')
  n = n + 1

  tests(n) = assert(all(split(arange(1_int64,10_int64),1) == arange(1,5)),&
       'split(x,1), int64')
  n = n + 1

  tests(n) = assert(all(split(arange(1_int64,10_int64),2) == arange(6,10)),&
       'split(x,2), int64')
  n = n + 1

  tests(n) = assert(all(split(arange(1._real32,10._real32),1) == arange(1,5)),&
       'split(x,1), real32')
  n = n + 1

  tests(n) = assert(all(split(arange(1._real32,10._real32),2) == arange(6,10)),&
       'split(x,2), real32')
  n = n + 1

  tests(n) = assert(all(split(arange(1._real64,10._real64),1) == arange(1,5)),&
       'split(x,1), real64')
  n = n + 1

  tests(n) = assert(all(split(arange(1._real64,10._real64),2) == arange(6,10)),&
       'split(x,2), real64')
  n = n + 1

  tests(n) = assert(all(split(arange(1._real128,10._real128),1) == arange(1,5)),&
       'split(x,1), real128')
  n = n + 1

  tests(n) = assert(all(split(arange(1._real128,10._real128),2) == arange(6,10)),&
       'split(x,2), real128')
  n = n + 1

  tests(n) = assert(all(split([1],1) == arange(1,0)),'split([1],1) returns an empty array')
  n = n + 1

  tests(n) = assert(all(split([1],2) == [1]),'split([1],2) returns [1]')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 11
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(subscript([1_int8,2_int8,3_int8],[2_int8]) == [2_int8]),&
       'subscript, int8')
  n = n + 1

  tests(n) = assert(all(subscript([1_int16,2_int16,3_int16],[2_int16]) == [2_int16]),&
       'subscript, int16')
  n = n + 1

  tests(n) = assert(all(subscript([1_int32,2_int32,3_int32],[2_int32]) == [2_int32]),&
       'subscript, int32')
  n = n + 1

  tests(n) = assert(all(subscript([1_int64,2_int64,3_int64],[2_int64]) == [2_int64]),&
       'subscript, int64')
  n = n + 1

  tests(n) = assert(all(subscript([1._real32,2._real32,3._real32],[2]) == [2._real32]),&
       'subscript, real32')
  n = n + 1

  tests(n) = assert(all(subscript([1._real64,2._real64,3._real64],[2]) == [2._real64]),&
       'subscript, real64')
  n = n + 1

  tests(n) = assert(all(subscript([1._real128,2._real128,3._real128],[2]) == [2._real128]),&
       'subscript, real128')
  n = n + 1

  tests(n) = assert(size(subscript([1,2,3],[0])) == 0,&
       'subscript out of bounds returns empty array')
  n = n + 1

  tests(n) = assert(all(subscript(arange(cmplx(1._real32,0._real32, kind=real32),&
       cmplx(3._real32,0._real32, kind=real32)),[2])&
       == [cmplx(2._real32,0._real32, kind=real32)]),&
       'subscript, complex real32')
  n = n + 1

  tests(n) = assert(all(subscript(arange(cmplx(1._real64,0._real64, kind=real64),&
       cmplx(3._real64,0._real64, kind=real64)),[2])&
       == [cmplx(2._real64,0._real64, kind=real64)]),&
       'subscript, complex real64')
  n = n + 1

  tests(n) = assert(all(subscript(arange(cmplx(1._real128,0._real128, kind=real128),&
       cmplx(3._real128,0._real128, kind=real128)),[2])&
       == [cmplx(2._real128,0._real128, kind=real128)]),&
       'subscript, complex real128')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  c_r4 = [(1,2),(2,4)]
  c_r8 = c_r4
  c_r16 = c_r4

  n = 1
  ntests = 13
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(tail([1_int8,2_int8]) == [2_int8]),'tail, int8')
  n = n + 1

  tests(n) = assert(all(tail([1_int16,2_int16]) == [2_int16]),'tail, int16')
  n = n + 1

  tests(n) = assert(all(tail([1_int32,2_int32]) == [2_int32]),'tail, int32')
  n = n + 1

  tests(n) = assert(all(tail([1_int64,2_int64]) == [2_int64]),'tail, int64')
  n = n + 1

  tests(n) = assert(all(tail([1._real32,2._real32]) == [2._real32]),'tail, real32')
  n = n + 1

  tests(n) = assert(all(tail([1._real64,2._real64]) == [2._real64]),'tail, real64')
  n = n + 1

  tests(n) = assert(all(tail([1._real128,2._real128]) == [2._real128]),'tail, real128')
  n = n + 1

  tests(n) = assert(all(tail(c_r4) == [c_r4(2)]),'tail, complex real32')
  n = n + 1

  tests(n) = assert(all(tail(c_r8) == [c_r8(2)]),'tail, complex real64')
  n = n + 1

  tests(n) = assert(all(tail(c_r16) == [c_r16(2)]),'tail, complex real128')
  n = n + 1

  tests(n) = assert(size(tail([1._real32])) == 0,'size(tail([1])) == 0')
  n = n + 1

  tests(n) = assert(size(tail(tail([1._real32]))) == 0,'size(tail(tail([1]))) == 0')
  n = n + 1

#ifdef __pgiFortran__
  tests(n) = assert(.true.,'skip tail operator, .tail.x')
#else
  tests(n) = assert(all(tail([1,2]) == .tail.[1,2]),'tail operator, .tail.x')
#endif
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 10
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(unfold(addone_i1,[1_int8],3_int8) == [1,2,3]),&
       'unfold, int8')
  n = n + 1

  tests(n) = assert(all(unfold(addone_i2,[1_int16],3_int16) == [1,2,3]),&
       'unfold, int16')
  n = n + 1

  tests(n) = assert(all(unfold(addone_i4,[1_int32],3_int32) == [1,2,3]),&
       'unfold, int32')
  n = n + 1

  tests(n) = assert(all(unfold(addone_i8,[1_int64],3_int64) == [1,2,3]),&
       'unfold, int64')
  n = n + 1

  tests(n) = assert(all(unfold(addone_r4,[1._real32],3_int32) == [1,2,3]),&
       'unfold, real32')
  n = n + 1

  tests(n) = assert(all(unfold(addone_r8,[1._real64],3_int32) == [1,2,3]),&
       'unfold, real64')
  n = n + 1

  tests(n) = assert(all(unfold(addone_r16,[1._real128],3_int32) == [1,2,3]),&
       'unfold, real128')
  n = n + 1

  tests(n) = assert(all(unfold(addone_c4,[cmplx(1._real32,0._real32, kind=real32)],3)&
       == arange(cmplx(1._real32,0._real32, kind=real32),&
       cmplx(3._real32,0._real32, kind=real32))),&
       'unfold, complex real32')
  n = n + 1

  c8 = [cmplx(1._real64,0._real64, kind=real64)]
  tests(n) = assert(all(unfold(addone_c8,c8,3)&
       == arange(cmplx(1._real64,0._real64, kind=real64),&
       cmplx(3._real64,0._real64, kind=real64))),&
       'unfold, complex real64')
  n = n + 1

  c16 = [cmplx(1._real128,0._real128, kind=real128)]
  tests(n) = assert(all(unfold(addone_c16,c16,3)&
       == arange(cmplx(1._real128,0._real128, kind=real128),&
       cmplx(3._real128,0._real128, kind=real128))),&
       'unfold, complex real128')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  n = 1
  ntests = 14
  call initialize_tests(tests,ntests)

  tests(n) = assert(all(union([1_int8,2_int8],[2_int8,3_int8]) == [1,2,3]),&
       'union, int8')
  n = n + 1

  tests(n) = assert(all(union([1_int16,2_int16],[2_int16,3_int16]) == [1,2,3]),&
       'union, int16')
  n = n + 1

  tests(n) = assert(all(union([1_int32,2_int32],[2_int32,3_int32]) == [1,2,3]),&
       'union, int32')
  n = n + 1

  tests(n) = assert(all(union([1_int64,2_int64],[2_int64,3_int64]) == [1,2,3]),&
       'union, int64')
  n = n + 1

  tests(n) = assert(all(union([1._real32,2._real32],[2._real32,3._real32]) == [1,2,3]),&
       'union, real32')
  n = n + 1

  tests(n) = assert(all(union([1._real64,2._real64],[2._real64,3._real64]) == [1,2,3]),&
       'union, real64')
  n = n + 1

  tests(n) = assert(all(union([1._real128,2._real128],[2._real128,3._real128]) == [1,2,3]),&
       'union, real128')
  n = n + 1

  tests(n) = assert(all(union([cmplx(1._real32,0._real32, kind=real32),cmplx(2._real32,0._real32, kind=real32)],&
       [cmplx(2._real32,0._real32, kind=real32),cmplx(3._real32,0._real32, kind=real32)])&
       == [cmplx(1._real32,0._real32, kind=real32),&
       cmplx(2._real32,0._real32, kind=real32),&
       cmplx(3._real32,0._real32, kind=real32)]),&
       'union, complex real32')
  n = n + 1

  tests(n) = assert(all(union([cmplx(1._real64,0._real64, kind=real64),cmplx(2._real64,0._real64, kind=real64)],&
       [cmplx(2._real64,0._real64, kind=real64),cmplx(3._real64,0._real64, kind=real64)])&
       == [cmplx(1._real64,0._real64, kind=real64),&
       cmplx(2._real64,0._real64, kind=real64),&
       cmplx(3._real64,0._real64, kind=real64)]),&
       'union, complex real64')
  n = n + 1

  tests(n) = assert(all(union([cmplx(1._real128,0._real128, kind=real128),cmplx(2._real128,0._real128, kind=real128)],&
       [cmplx(2._real128,0._real128, kind=real128),cmplx(3._real128,0._real128, kind=real128)])&
       == [cmplx(1._real128,0._real128, kind=real128),&
       cmplx(2._real128,0._real128, kind=real128),&
       cmplx(3._real128,0._real128, kind=real128)]),&
       'union, complex real128')
  n = n + 1

  tests(n) = assert(all(union(arange(1,0),arange(1,0)) == arange(1,0)),&
       'union of empty arrays is an empty array')
  n = n + 1

  tests(n) = assert(all(union([1,2,2,3],arange(1,0)) == set([1,2,2,3])),&
       'union(x,[]) == set(x)')
  n = n + 1

  tests(n) = assert(all(union(arange(1,0),[1,2,2,3]) == set([1,2,2,3])),&
       'union([],x) == set(x)')
  n = n + 1

  tests(n) = assert(all(union([1,2],[3,4]) == ([1,2].union.[3,4])),&
       'union operator, a.union.b')
  n = n + 1

  test_failed = .false.
  call report_tests(tests,test_failed)
  isgood = isgood .and. (.not. test_failed)

  write(*,*) ''
  if (isgood) then
     write(*,*) 'mo_functional ', color('o.k.', c_green)
  else
     write(*,*) 'mo_functional ', color('failed!', c_red)
  endif

end program main


text:
!
! wavy - A spectral ocean wave modeling and development framework
! Copyright (c) 2017, Wavebit Scientific LLC
! All rights reserved.
!
! Licensed under the BSD-3-clause license. See LICENSE for details.

module mod_io

use mod_precision,only:ik => intkind,rk => realkind

implicit none

private
contains
!-------------------------------------------------------------------------------
endmodule mod_io

max_stars_repo_name: dksasaki/wavy

max_issues_repo_name: dksasaki/wavy

max_forks_repo_name: dksasaki/wavy


text:
7Eleven is a Convenience Stores convenience store just like every other 7Eleven from Kathmandu to Kentucky! They offer all kinds of sundries from tampons to Ben & Jerrys ice cream, snacks, sandwiches, hot dogs, Forties beer, & Smoking smokes. They are most widely known for their Frozen Drinks Slurpees (frozen soft drinks). The store also has a fountain drinks soda fountain. There really is a 7Eleven in Kathmandu as well as in Bangkok. When you purchase nachos, you get an unlimited selfserve supply of hot nacho cheese.

Every 0711 (July 11) of each year, 7Eleven offers free 7.11 ounce size cups of Slurpees.

There is a Citibank branded ATMs ATM inside the store which Citibank users can use without paying a surcharge. Since there is no Citibank in Davis, Citibank users can use the ATM here without any fees. It is a Vcom brand ATM with deposit capability. The ATM is also part of the CoOp ATM network and is free for members of all Davis Credit Unions credit unions and virtually all other credit unions across the country as well. There is also a East Davis Phones payphone for public use.

7Eleven celebrated the launch of Slurpee Lite (sugar free versions) with SlurpFREE Day. On May 23, 2012 guests from coast to coast received FREE 7.11 ounce Slurpees.

Twelve 7Eleven stores were temporarily transformed into KwikEMarts in anticipation of The Simpsons Movie. The promotion ended August 1, 2007.

7Eleven stores featured Simpsons inspired products: Squishee drinks, KrustyOs Cereal, Buzz Cola, Pink Movie Donuts, and Radioactive Man Comic Book Edition #711.

As of March 15, 2010, the Davis 7Eleven location has a convenient Redbox DVD rental machine in front of the store.

The store changed hands and replaced most or all of their staff in mid2012.  It was also remodeled. As of August 2012 this location is renting out kegs. The following brands are available as kegs: Coors Kegs, Miller kegs, Budweiser kegs.

Older comments have been archived /archive here.

20110712 01:23:26 nbsp After walking in post Free Slurpee Day all I can say is some people are disgusting, greedy slobs. There was such a mess it was utterly ridiculous... The staff was cleaning for almost two hours trying to recover from the insanity. Users/WesP
  But you probably already knew that from working at Little Prague. Users/hankim

20120728 10:16:29 nbsp I just went into the store, they remodeled and its looking great now plus theres new staff. Theyre also giving out free small slurpees! Users/SoupJill
  Removed the entirety of the old staff, even those who worked their tails off for the company for several years, and threw a new coat of paint on old dingy walls... Store went under new ownership for losing money quarterly, now you continue the trend as I am not setting foot in the locale again. Users/WesP

20120729 19:17:32 nbsp I dont understand what your problem is... I am just conveying my opinion. Users/SoupJill
  And http://daviswiki.org/7Eleven?actiondiff&version267&version156 removing the opinions of seven other people.  It happens, and can be forgiven, but all edits are openly displayed.  ⁓Users/JabberWokky ʝ⍵

20120731 08:36:36 nbsp I am sorry. I figured that since they are so old that they are no longer accurate... Users/SoupJill

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program higherorder
  real :: a,b,c,d,e
  d = f(g,b,c)
  e = f(g,b,c)
contains
    ! need a mechanism to express that all uses of h need to have the same units

  real function f(h,x,y)
    != unit m :: x
    != unit s :: y

    real :: x, y
    real, external :: h
    f = h(x) + h(y)
  end function f

  real function g(y)
    real :: y
    g = y
  end function g
end program higherorder

max_stars_repo_name: raehik/camfort

max_issues_repo_name: apthorpe/camfort

max_forks_repo_name: raehik/camfort


text:
      DOUBLE PRECISION FUNCTION MA02JZ( LTRAN1, LTRAN2, N, Q1, LDQ1, Q2,
     $                                  LDQ2, RES, LDRES )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To compute || Q^H Q - I ||_F for a complex matrix of the form
C
C                       [  op( Q1 )  op( Q2 ) ]
C                  Q =  [                     ],
C                       [ -op( Q2 )  op( Q1 ) ]
C
C     where Q1 and Q2 are N-by-N matrices. This residual can be used to
C     test wether Q is numerically a unitary symplectic matrix.
C
C     FUNCTION VALUE
C
C     MA02JZ  DOUBLE PRECISION
C             The computed residual.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     LTRAN1  LOGICAL
C             Specifies the form of op( Q1 ) as follows:
C             = .FALSE.:  op( Q1 ) = Q1;
C             = .TRUE. :  op( Q1 ) = Q1'.
C
C     LTRAN2  LOGICAL
C             Specifies the form of op( Q2 ) as follows:
C             = .FALSE.:  op( Q2 ) = Q2;
C             = .TRUE. :  op( Q2 ) = Q2'.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrices Q1 and Q2.  N >= 0.
C
C     Q1      (input) COMPLEX*16 array, dimension (LDQ1,N)
C             On entry, the leading N-by-N part of this array must
C             contain the matrix op( Q1 ).
C
C     LDQ1    INTEGER
C             The leading dimension of the array Q1.  LDQ1 >= MAX(1,N).
C
C     Q2      (input) COMPLEX*16 array, dimension (LDQ2,N)
C             On entry, the leading N-by-N part of this array must
C             contain the matrix op( Q2 ).
C
C     LDQ2    INTEGER
C             The leading dimension of the array Q2.  LDQ2 >= MAX(1,N).
C
C     Workspace
C
C     RES     DOUBLE PRECISION array, dimension (LDRES,N)
C
C     LDRES   INTEGER
C             The leading dimension of the array RES.
C             LDRES >= MAX(1,N).
C
C     METHOD
C
C     The routine computes the residual by simple elementary operations.
C
C     CONTRIBUTORS
C
C     V. Sima, Research Institute for Informatics, Bucharest, Romania,
C     Sep. 2012. Based on SLICOT Library routine MA02JD.
C
C     REVISIONS
C
C     -
C
C     KEYWORDS
C
C     Elementary operations.
C
C     ******************************************************************
C
C     .. Parameters ..
      COMPLEX*16        ZERO, ONE
      PARAMETER        ( ZERO = (0.0D0,0.0D0), ONE = (1.0D0,0.0D0) )
      DOUBLE PRECISION  TWO
      PARAMETER         ( TWO = 2.0D0 )
C     .. Scalar Arguments ..
      LOGICAL           LTRAN1, LTRAN2
      INTEGER           LDQ1, LDQ2, LDRES, N
C     .. Array Arguments ..
      COMPLEX*16        Q1(LDQ1,*), Q2(LDQ2,*), RES(LDRES,*)
C     .. Local Scalars ..
      INTEGER           I
      DOUBLE PRECISION  TEMP
C     .. Local Arrays ..
      DOUBLE PRECISION  DUMMY(1)
C     .. External Subroutines ..
      EXTERNAL          ZGEMM
C     .. External Functions ..
      DOUBLE PRECISION  DLAPY2, ZLANGE
      EXTERNAL          DLAPY2, ZLANGE
C     .. Intrinsic Functions ..
      INTRINSIC         SQRT
C
C     .. Executable Statements ..
C
      IF ( LTRAN1 ) THEN
         CALL ZGEMM( 'No Transpose', 'Conj Transpose', N, N, N, ONE, Q1,
     $               LDQ1, Q1, LDQ1, ZERO, RES, LDRES )
      ELSE
         CALL ZGEMM( 'Conj Transpose', 'No Transpose', N, N, N, ONE, Q1,
     $               LDQ1, Q1, LDQ1, ZERO, RES, LDRES )
      END IF
      IF ( LTRAN2 ) THEN
         CALL ZGEMM( 'No Transpose', 'Conj Transpose', N, N, N, ONE, Q2,
     $               LDQ2, Q2, LDQ2, ONE, RES, LDRES )
      ELSE
         CALL ZGEMM( 'Conj Transpose', 'No Transpose', N, N, N, ONE, Q2,
     $               LDQ2, Q2, LDQ2, ONE, RES, LDRES )
      END IF
      DO 10 I = 1, N
         RES(I,I) = RES(I,I) - ONE
   10 CONTINUE
      TEMP = ZLANGE( 'Frobenius', N, N, RES, LDRES, DUMMY )
      IF ( LTRAN1 .AND. LTRAN2 ) THEN
         CALL ZGEMM( 'No Transpose', 'Conj Transpose', N, N, N, ONE, Q2,
     $               LDQ2, Q1, LDQ1, ZERO, RES, LDRES )
         CALL ZGEMM( 'No Transpose', 'Conj Transpose', N, N, N, ONE, Q1,
     $               LDQ1, Q2, LDQ2, -ONE, RES, LDRES )
      ELSE IF ( LTRAN1 ) THEN
         CALL ZGEMM( 'Conj Transpose', 'Conj Transpose', N, N, N, ONE,
     $               Q2, LDQ2, Q1, LDQ1, ZERO, RES, LDRES )
         CALL ZGEMM( 'No Transpose', 'No Transpose', N, N, N, ONE, Q1,
     $               LDQ1, Q2, LDQ2, -ONE, RES, LDRES )
      ELSE IF ( LTRAN2 ) THEN
         CALL ZGEMM( 'No Transpose', 'No Transpose', N, N, N, ONE, Q2,
     $               LDQ2, Q1, LDQ1, ZERO, RES, LDRES )
         CALL ZGEMM( 'Conj Transpose', 'Conj Transpose', N, N, N, ONE,
     $               Q1, LDQ1, Q2, LDQ2, -ONE, RES, LDRES )
      ELSE
         CALL ZGEMM( 'Conj Transpose', 'No Transpose', N, N, N, ONE, Q2,
     $               LDQ2, Q1, LDQ1, ZERO, RES, LDRES )
         CALL ZGEMM( 'Conj Transpose', 'No Transpose', N, N, N, ONE, Q1,
     $               LDQ1, Q2, LDQ2, -ONE, RES, LDRES )
      END IF
      TEMP = DLAPY2( TEMP, ZLANGE( 'Frobenius', N, N, RES, LDRES,
     $                             DUMMY ) )
      MA02JZ = SQRT( TWO )*TEMP
      RETURN
C *** Last line of MA02JZ ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: bnavigator/SLICOT-Reference

max_forks_repo_name: bnavigator/SLICOT-Reference


text: #redirect Jennifer Beeman


text:
      SUBROUTINE FA1PKV (AZ,AMK,AMB,N,E1,CZ,BREF,PI,VEL,IBUF)
C
      INTEGER         IBUF(1),IV(6),TRL(7)
      REAL            E1(5),V(6),E(2),AMK(1),AMB(1),CZ(1)
      COMPLEX         AZ(1),CEIG,EIGEN,EIGZ
      COMMON /SYSTEM/ SYSBUF,NOUT,SPACE(6),NLPP,X(2),LINES
      EQUIVALENCE     (V(1),IV(1)),(EIGEN,E(1))
      DATA    ISCR  / 301/, IPASS /0/
C
      EIGZ = (0.0,0.0)
      IF (N .LT. 2) GO TO 1000
      E(1) = E1(1)
      E(2) = E1(2)
      IF (IPASS .NE. 0) GO TO 5
      CALL OPEN (*1000,ISCR,IBUF,1)
      GO TO 9
    5 CALL OPEN (*1000,ISCR,IBUF,3)
    9 IPASS = IPASS + 1
C
C     BUILD A = IP2 + M-1B P + M-1K
C
      CEIG = EIGEN*EIGEN
      K = 0
      DO 10 I = 1,N
      DO 10 J = 1,N
      K = K + 1
      AZ(K) = -AMB(K)*EIGEN - AMK(K)
      IF (I .EQ. J) AZ(K) = AZ(K) + CEIG
   10 CONTINUE
C
C     CORE FOR EGNVCT
C
      N2 = N*2
      NA = 1  + N2*N
      NB = NA + N2
      NC = NB + N2
      ND = NC + N2
      CALL EGNVCT (AZ,CZ(NA),EIGZ,CZ(NB),CZ(NC),CZ(ND),N)
C
C     BUILD ON SCR1 DATA FOR VECTOR OUTPUT
C
      IV(1) = IPASS
      IV(2) = IPASS
      V (3) = E1(1)
      V (4) = E1(2)
      IF (E1(2) .EQ. 0.0) GO TO 20
      V(5)  = E1(3)
      V(6)  = E1(5)
      GO TO 22
   20 V(5)  = 0.0
      V(6)  = (BREF/(.34657*VEL))*E1(1)
   22 CALL WRITE (ISCR,IV,6,1)
      CALL WRITE (ISCR,CZ(NB),N2,1)
C
C     VECTOR IS IN CZ(NB)
C
      LINES = NLPP
      K = 0
      DO 30 I = 1,N
      IF (LINES .LT. NLPP) GO TO 25
      CALL PAGE1
      WRITE  (NOUT,21) EIGEN
   21 FORMAT (1H0,47X,30HEIGENVECTOR FROM THE PK METHOD, /3X,
     1        13HEIGENVALUE = ,1P,E15.5,1P,E15.5, //3X,11HEIGENVECTOR)
      LINES = LINES + 5
   25 LINES = LINES + 1
      WRITE  (NOUT,26) CZ(NB+K),CZ(NB+K+1)
   26 FORMAT (16X,1P,E15.5,1P,E15.5)
      K = K + 2
   30 CONTINUE
      TRL(1) = ISCR
      TRL(2) = 1
      CALL WRTTRL (TRL)
 1000 CALL CLOSE (ISCR,3)
      RETURN
      END