text:
!** Copyright (c) 2003, 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.

! check whether inlining with optional arguments works
! if the argument is not type scalar integer
! this test: inline with derived type arguments
module m
 type dt
  integer :: m1
  real :: m2
  integer, dimension(3) :: m3
 end type
 contains

 subroutine s1( a, b, c )
  type(dt) :: a
  integer :: b
  type(dt),optional :: c

  a%m1 = b+1
  a%m2 = b*2
  a%m3(1) = b+2
  a%m3(2) = b+3
  a%m3(3) = b+4
  if( present(c) ) then
   c%m1 = b+5
   c%m2 = b*3
   c%m3(1) = b+6
   c%m3(2) = b+7
   c%m3(3) = b+8
  endif
 end subroutine
end module

program p
 use m
 type(dt) a,c
 integer b
 integer result(15)
 integer expect(15)
 data expect /101,200,102,103,104,105, 300,106,107,108,201,400,202,203,204/

 b = 100
 call s1( a, b, c )
 result( 1) = a%m1
 result( 2) = a%m2
 result( 3) = a%m3(1)
 result( 4) = a%m3(2)
 result( 5) = a%m3(3)
 result( 6) = c%m1
 result( 7) = c%m2
 result( 8) = c%m3(1)
 result( 9) = c%m3(2)
 result(10) = c%m3(3)
 call s1( a, 200 )
 result(11) = a%m1
 result(12) = a%m2
 result(13) = a%m3(1)
 result(14) = a%m3(2)
 result(15) = a%m3(3)
 !print *,result
 call check(result,expect,15)
end program

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


!
! libui-mebu.f90
! Fortran から　libui　を使う。
!
! Author Yoji.Hosokawa 2016/4
!
! gfortran -O2 -g -Wall -o libui-window libui-window.f90 -L.-lui `pkg-config gtk+-3.0 --libs`
!

module handlers
  use iso_c_binding
  implicit none

  type(c_ptr) :: mainwin
  type(c_ptr) :: entry1,entry2

  type(c_ptr) :: fileMenu
  type(c_ptr) :: newItem
  type(c_ptr) :: openItem
  type(c_ptr) :: saveItem
  type(c_ptr) :: shouldQuitItem
  type(c_ptr) :: quitItem
  type(c_ptr) :: editMenu
  type(c_ptr) :: undoItem
  type(c_ptr) :: checkItem
  type(c_ptr) :: accelItem
  type(c_ptr) :: prefsItem
  type(c_ptr) :: testMenu
  type(c_ptr) :: enabledItem
  type(c_ptr) :: enableThisItem
  type(c_ptr) :: forceCheckedItem
  type(c_ptr) :: forceUncheckedItem
  type(c_ptr) :: whatWindowItem
  type(c_ptr) :: moreTestsMenu
  type(c_ptr) :: quitEnabledItem
  type(c_ptr) :: prefsEnabledItem
  type(c_ptr) :: aboutEnabledItem
  type(c_ptr) :: checkEnabledItem
  type(c_ptr) :: multiMenu
  type(c_ptr) :: helpMenu
  type(c_ptr) :: helpItem
  type(c_ptr) :: aboutItem

contains

  ! "close mainwindow" is a this app

  subroutine onClosing()
    use iso_c_binding
    use libuif
      call uiQuit()
  end subroutine onClosing

  function onShouldQuit(d) result(ret)
    use iso_c_binding
    use libuif
      type(c_ptr), value :: d
      integer :: ret
        if (uiMenuItemChecked(shouldQuitItem) == 1) then
          call uiControlDestroy(d)
          ret = 1
          return
        end if
      ret = 0
  end function onShouldQuit

  subroutine onOpenFileClicked()
    use iso_c_binding
    use libuif

    type(c_ptr) :: filename_c
    integer :: i
    character(1000, C_CHAR),pointer :: filename_f

    filename_c = uiOpenFile(mainwin)
    if (C_ASSOCIATED(filename_c) .neqv. .true.) then
      call uiEntrySetText(entry1, "(cancelled)"//c_null_char)
    else
      call c_f_pointer(filename_c,filename_f)
      i = index(filename_f,char(0))
      call uiEntrySetText(entry1, filename_f(1:i))
      call uiFreeText(filename_f)
    end if

  end subroutine onOpenFileClicked

  subroutine onSaveFileClicked()
    use iso_c_binding
    use libuif

    type(c_ptr) :: filename_c
    integer :: i
    character(1000, C_CHAR),pointer :: filename_f

    filename_c = uiSaveFile(mainwin)
    if (C_ASSOCIATED(filename_c)) then
      call c_f_pointer(filename_c,filename_f)
      i = index(filename_f,char(0))
      call uiEntrySetText(entry2, filename_f(1:i))
      call uiFreeText(filename_f)
    else
      call uiEntrySetText(entry2, "(cancelled)"//c_null_char)
    end if

  end subroutine onSaveFileClicked

  subroutine onMsgBoxClicked()
    use iso_c_binding
    use libuif
      call uiMsgBox(mainwin,"このアプリについて"//c_null_char, &
&          "これはlibuiとFortranを使って作成した" // char(13) // 'ＧＵＩのサンプルです。'//c_null_char)
  end subroutine onMsgBoxClicked

  subroutine enableItemTest(item, w, d)
    use iso_c_binding
    use libuif
      type(c_ptr),value :: item, w, d

      if (uiMenuItemChecked(item) == 1) then
        call uiMenuItemEnable(d)
      else
        call uiMenuItemDisable(d)
      end if
  end subroutine enableItemTest

  subroutine forceOn()
    use iso_c_binding
    use libuif
      call uiMenuItemSetChecked(enabledItem, 1)
  end subroutine forceOn

  subroutine forceOff()
    use iso_c_binding
    use libuif
      call uiMenuItemSetChecked(enabledItem, 0)
  end subroutine forceOff

  subroutine whatWindow()
    use iso_c_binding
    use libuif
        print *, "menu item clicked on window %p\n"
  end subroutine whatWindow

  subroutine initMenus()
    use iso_c_binding
    use libuif

    type(c_ptr) :: dummy

      fileMenu = uiNewMenu("File"//c_null_char)
      newItem = uiMenuAppendItem(fileMenu, "New"//c_null_char)
      openItem = uiMenuAppendItem(fileMenu, "Open"//c_null_char)
      call uiMenuItemOnClicked(openItem, c_funloc(onOpenFileClicked), c_null_ptr)
      saveItem = uiMenuAppendItem(fileMenu, "保存"//c_null_char)
      call uiMenuItemOnClicked(saveItem, c_funloc(onSaveFileClicked), c_null_ptr)
      call uiMenuAppendSeparator(fileMenu)
      shouldQuitItem = uiMenuAppendCheckItem(fileMenu, "Should Quit"//c_null_char)
      quitItem = uiMenuAppendQuitItem(fileMenu)

      editMenu = uiNewMenu("Edit"//c_null_char)
      undoItem = uiMenuAppendItem(editMenu, "Undo"//c_null_char)
      call uiMenuItemDisable(undoItem)
      call uiMenuAppendSeparator(editMenu)
      checkItem = uiMenuAppendCheckItem(editMenu, "Check Me\tTest"//c_null_char)
      accelItem = uiMenuAppendItem(editMenu, "A&ccele&&rator T_es__t"//c_null_char)
      prefsItem = uiMenuAppendPreferencesItem(editMenu)

      testMenu = uiNewMenu("Test"//c_null_char)
      enabledItem = uiMenuAppendCheckItem(testMenu, "Enable Below Item"//c_null_char)
      call uiMenuItemSetChecked(enabledItem, 1)
      enableThisItem = uiMenuAppendItem(testMenu, "This Will Be Enabled"//c_null_char)
      call uiMenuItemOnClicked(enabledItem, c_funloc(enableItemTest), enableThisItem)
      forceCheckedItem = uiMenuAppendItem(testMenu, "Force Above Checked"//c_null_char)
      call uiMenuItemOnClicked(forceCheckedItem, c_funloc(forceOn), c_null_ptr)
      forceUncheckedItem = uiMenuAppendItem(testMenu, "Force Above Unchecked"//c_null_char)
      call uiMenuItemOnClicked(forceUncheckedItem, c_funloc(forceOff), c_null_ptr)
      call uiMenuAppendSeparator(testMenu)
      whatWindowItem = uiMenuAppendItem(testMenu, "What Window?"//c_null_char)
      call uiMenuItemOnClicked(whatWindowItem, c_funloc(whatWindow), c_null_ptr)

      moreTestsMenu = uiNewMenu("More Tests"//c_null_char)
      quitEnabledItem = uiMenuAppendCheckItem(moreTestsMenu, "Quit Item Enabled"//c_null_char)
      call uiMenuItemSetChecked(quitEnabledItem, 1)
      prefsEnabledItem = uiMenuAppendCheckItem(moreTestsMenu, "Preferences Item Enabled"//c_null_char)
      call uiMenuItemSetChecked(prefsEnabledItem, 1)
      aboutEnabledItem = uiMenuAppendCheckItem(moreTestsMenu, "About Item Enabled"//c_null_char)
      call uiMenuItemSetChecked(aboutEnabledItem, 1)
      call uiMenuAppendSeparator(moreTestsMenu)
      checkEnabledItem = uiMenuAppendCheckItem(moreTestsMenu, "Check Me Item Enabled"//c_null_char)
      call uiMenuItemSetChecked(checkEnabledItem, 1)

      multiMenu = uiNewMenu("Multi"//c_null_char)
      call uiMenuAppendSeparator(multiMenu)
      call uiMenuAppendSeparator(multiMenu)
      dummy = uiMenuAppendItem(multiMenu, "Item && Item && Item"//c_null_char)
      call uiMenuAppendSeparator(multiMenu)
      call uiMenuAppendSeparator(multiMenu)
      dummy = uiMenuAppendItem(multiMenu, "Item __ Item __ Item"//c_null_char)
      call uiMenuAppendSeparator(multiMenu)
      call uiMenuAppendSeparator(multiMenu)

      helpMenu = uiNewMenu("Help"//c_null_char)
      helpItem = uiMenuAppendItem(helpMenu, "Help"//c_null_char)
      aboutItem = uiMenuAppendAboutItem(helpMenu)

      call uiMenuItemOnClicked(quitEnabledItem, c_funloc(enableItemTest), quitItem)
      call uiMenuItemOnClicked(prefsEnabledItem, c_funloc(enableItemTest), prefsItem)
      call uiMenuItemOnClicked(aboutEnabledItem, c_funloc(enableItemTest), aboutItem)
      call uiMenuItemOnClicked(checkEnabledItem, c_funloc(enableItemTest), checkItem)

      call uiMenuItemOnClicked(aboutItem, c_funloc(onMsgBoxClicked), c_null_ptr)

  end subroutine initMenus

end module handlers

program main
    use iso_c_binding
    use libuif
    use handlers

    implicit none
!    type(c_ptr) :: mainwin
    type(c_ptr) :: vbox,ebtn,grid,button

! initialize
  call uiInit(0)

! start procedure

  call initMenus()

  mainwin = uiNewWindow("libui menu"//c_null_char, 400, 100, 1)
  call uiWindowSetMargined(mainwin, 1)
  call uiWindowOnClosing(mainwin, c_funloc(onClosing), c_null_ptr)

  call uiOnShouldQuit(c_funloc(onShouldQuit), mainwin)

  vbox = uiNewVerticalBox()
  call uiWindowSetChild(mainwin, vbox)

  grid = uiNewGrid()
  call uiGridSetPadded(grid, 1)
  call uiBoxAppend(vbox, grid, 0)

  button = uiNewButton("Open File"//c_null_char)
  entry1 = uiNewEntry()
  call uiEntrySetReadOnly(entry1, 1)
  call uiButtonOnClicked(button, c_funloc(onOpenFileClicked), c_null_ptr)
  call uiGridAppend(grid, button,0, 0, 1, 1,0, uiAlignFill, 0, uiAlignFill)
  call uiGridAppend(grid, entry1,1, 0, 1, 1,1, uiAlignFill, 0, uiAlignFill)

  button = uiNewButton("Save File"//c_null_char)
  entry2 = uiNewEntry()
  call uiEntrySetReadOnly(entry2, 1)
  call uiButtonOnClicked(button, c_funloc(onSaveFileClicked), c_null_ptr)
  call uiGridAppend(grid, button,0, 1, 1, 1,0, uiAlignFill, 0, uiAlignFill)
  call uiGridAppend(grid, entry2,1, 1, 1, 1,1, uiAlignFill, 0, uiAlignFill)

  call uiBoxAppend(vbox,uiNewHorizontalSeparator(),0)

  ebtn = uiNewButton("終了"//c_null_char)
  call uiButtonOnClicked(ebtn, c_funloc(onClosing), c_null_ptr)

  call uiBoxAppend(vbox, ebtn, 0)

  call uiControlShow(mainwin)
  call uiMain()

! nomary do not need if call when close window receive sig fault by button clicked
!  call uiUninit()

end program main


text:
program ex2p3
    use MC_randoms
    implicit none
    !Number of bins, shots, tests and expected value
    integer(kind=ik), parameter:: M=100, N=1000000, T=100, E=N/M
    integer(kind=ik), dimension(M):: y=0
    integer(kind=ik):: i, j
    !Variables for results
    real(kind=rk):: chi=0, seed0, rn
    real(kind=rk), dimension(T):: chi_arr

    !We'll use Fortran's RNG to change seed between runs
    call srand(1237654)
    seed0=irand()

    !-----------------------------------------------------------------------------------
    ! LCG
    !-----------------------------------------------------------------------------------
    call seed_lcg(int(seed0,ik))
    do j=1, T
        !Generate a number and place it in an appropriate bin
        !Since the generator gives a number between 0 and 1, we can multiply it by the number of bins
        !and round it to integer to find its correct bin as change to integer rounds down.
        do i=1, N
            rn=rand_lcg()
            y(int(rn*M+1, ik)) = y(int(rn*M+1, ik))+1
        end do

        !Chi-squared for this run
        do i=1, M
            chi=chi+((y(i)-E)**2)/E
        end do

        !Save the result and get a new sequence
        chi_arr(j)=chi
        chi=0
        y=0
        seed0=irand()
        call seed_lcg(int(seed0, ik))
    end do

    !Save the data for analysis in Python
    open(1,action="write",file="data_lcg.txt",status="replace")
    do i=1,100
        write(1,*) chi_arr(i)
    end do
    close(1)

    print*, "LCG döne"

    !-----------------------------------------------------------------------------------
    !Park-Miller
    !-----------------------------------------------------------------------------------
    call seed_pm(int(seed0,ik))
    do j=1, T
        do i=1, N
            rn=rand_pm()
            y(int(rn*M+1, ik)) = y(int(rn*M+1, ik))+1
        end do

        do i=1, M
            chi=chi+((y(i)-E)**2)/E
        end do

        chi_arr(j)=chi
        chi=0
        y=0
        seed0=irand()
        call seed_pm(int(seed0, ik))
    end do

    open(1,action="write",file="data_pm.txt",status="replace")
    do i=1,100
        write(1,*) chi_arr(i)
    end do
    close(1)

    print*, "PM döne"

    !-----------------------------------------------------------------------------------
    !Mersenne Twister
    !-----------------------------------------------------------------------------------
    call init_genrand64(int(seed0,ik))
    do j=1, T
        do i=1, N
            rn=genrand64_real2()
            y(int(rn*M+1, ik)) = y(int(rn*M+1, ik))+1
        end do

        do i=1, M
            chi=chi+((y(i)-E)**2)/E
        end do

        chi_arr(j)=chi
        chi=0
        y=0
        seed0=irand()
        call init_genrand64(int(seed0,ik))
    end do

    open(1,action="write",file="data_mt.txt",status="replace")
    do i=1,100
        write(1,*) chi_arr(i)
    end do
    close(1)

    print*, "MT döne"

    !-----------------------------------------------------------------------------------
    !QCG
    !-----------------------------------------------------------------------------------
    call seed_qcg(int(seed0,ik))
    do j=1, T
        do i=1, N
            rn=rand_qcg()
            !There were strange things happening with some seeds, so I added a condition just to be safe
            !QCG only gives the same 11 values but somehow it managed to overflow
            if(rn<0) then
                rn=rand_qcg()
            end if
            y(int(rn*M+1, ik)) = y(int(rn*M+1, ik))+1
        end do

        do i=1, M
            chi=chi+((y(i)-E)**2)/E
        end do

        chi_arr(j)=chi
        chi=0
        y=0
        seed0=irand()
        call seed_qcg(int(seed0, ik))
    end do

    open(1,action="write",file="data_qcg.txt",status="replace")
    do i=1,100
        write(1,*) chi_arr(i)
    end do
    close(1)

    print*, "QCG döne"

    !-----------------------------------------------------------------------------------

end program ex2p3

max_stars_repo_name: toicca/basics-of-mc

max_issues_repo_name: toicca/basics-of-mc

max_forks_repo_name: toicca/basics-of-mc


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!

! automatic arrays in namelist

SUBROUTINE S1(I, yy, output)
integer :: yy(i), zz(I), aa(I), ii
NAMELIST /NLIST/ yy, zz, ii
character*40 output(8)
ii    = 5
zz(1) = 3
zz(2) = 4
!write(6     ,NML=NLIST)
write(output,NML=NLIST)
END SUBROUTINE S1

program test
integer aa(3)
data aa/1,2,-99/
character*40 output(8)
integer yy(2), zz(2), ii
NAMELIST /NLIST/ yy, zz, ii
integer result(5)
integer expect(5)
data expect/5,1,2,3,4/
call S1(2,aa,output)
!write(6, '(a)') output
read(output, nlist)
!write(6, '(6i4)') yy, zz, ii

result(1) = ii
result(2:3) = yy
result(4:5) = zz
call check(result, expect, 5)
end

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


!> This module writes the parameters in the
!! normalisations namelist to the netcdf output file.
!! These parameters play no part in the simulation,
!! but will typically have been used in calculating
!! the dimensionless inputs, and will typically be used
!! to dimensionalise the outputs.
module diagnostics_normalisations
contains
  subroutine write_normalisations(gnostics)
    use normalisations, only: norms
    use diagnostics_config, only: diagnostics_type
    use diagnostics_create_and_write, only: create_and_write_variable_noread
    implicit none
    type(diagnostics_type), intent(in) :: gnostics
    call create_and_write_variable_noread(gnostics, gnostics%rtype, "mref", "", &
         "Reference mass in atomic mass units ", "a.m.u.", norms%get_value("mref"))
    call create_and_write_variable_noread(gnostics, gnostics%rtype, "zref", "", &
         "Reference charge", "proton charge", norms%get_value("zref"))
    call create_and_write_variable_noread(gnostics, gnostics%rtype, "nref", "", &
         "The reference density ", "m^-3", norms%get_value("nref"))
    call create_and_write_variable_noread(gnostics, gnostics%rtype, "tref", "", &
         "The reference temperature ", "eV", norms%get_value("tref"))
    call create_and_write_variable_noread(gnostics, gnostics%rtype, "aref", "", &
         "The reference length ", "m", norms%get_value("aref"))
    call create_and_write_variable_noread(gnostics, gnostics%rtype, "vref", "", &
         "The reference (thermal) velocity ", "m/s", norms%get_value("vref"))
    call create_and_write_variable_noread(gnostics, gnostics%rtype, "bref", "", &
         "The reference magnetic field ", "Tesla", norms%get_value("bref"))
    call create_and_write_variable_noread(gnostics, gnostics%rtype, "rhoref", "", &
         "The reference larmour radius ", "m", norms%get_value("rhoref"))
  end subroutine write_normalisations

end module diagnostics_normalisations


text:
*     ------------------------------------------------------------
*     I N I T M
*     ------------------------------------------------------------
*
*     Determine the Rydberg constant, RMASS constant, and
*     initialize the radial grid
*
      SUBROUTINE INITM(IPROB)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (NWD=60,NOD=220)
*     MPI stuff ***********************************************
*
        INCLUDE 'mpif.h'
        parameter (MAXPROC=9999)
        common /MPI/ myid, nprocs, ierr, istat(MPI_STATUS_SIZE)
      PARAMETER (NELEMENT=50)

****************************************************************
*
      COMMON /INOUT/IREAD,IWRITE,ISCW,iuc,iuw,ioul,iouj,
     ;             iouhn,iouhz,iouhs,iouhm,iLS,idisk
*
      COMMON /PARAM/ H,H1,H3,CH,EH,RHO,Z,TOL,NO,ND,NWF,MASS,NCFG,IB,IC,
     :ID,D0,D1,D2,D3,D4,D5,D6,D8,D10,D12,D16,D30,FINE,NSCF,NCLOSD,RMASS
*
      COMMON /RADIAL/R(NOD),RR(NOD),R2(NOD),YK(NOD),YR(NOD),X(NOD)
*
      COMMON /RYDBERG/ RYDBRG,ZMU,ENERGY
*
      LOGICAL CORE,YM,YF
      COMMON /COR/COREF,COREG,CORER,CORES,ZM1,ZM2,FZE,FZA,DR2,CORE,YM,YF
*
      CHARACTER YES
      DOUBLE PRECISION ME
      REAL MN(NELEMENT),FZ(NELEMENT-9)
*
      DATA ME/548.579903D-6/,RY/109737.31534/ MN/1,4,7,9,11,12,14,16,
     : 19,20,23,24,27,28,31,32,35,40,39,40,45,48,51,52,55,56,59,58,
     : 63,64,69,74,75,80,79,84,85,88,89,90,93,98,99,102,103,108,
     : 107,114,115,120/
*
      DATA FZ/163.5,198.4,239.6,283.9,332.7,386.1,444.2,507.8,
     : 576.4,650.7,731.2,817.9,911.2,1011.4,1119.3,1234.9,1359.3,
     : 1491.9,1634.4,1786.6,1948.7,2122.7,2308.2,2507.4,2718.5,
     : 2944.3,3187.5,3442.2,3720.3,4015.4,4326.6,4664.8,5022.4,
     : 5405.9,5811.9,6245.0,6706.8,7201.3,7728.7,8296.5,8889.3/
*

*
* ****   DETERMINE THE RYDBERG CONSTANT AND F(Z) PARAMETER
*
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
	if (myid.eq.0) then

      WRITE(iscw,'(/1X,A)') 'Default Rydberg constant (y/n)'
      READ(iread,'(A1)') YES
      IF (YES.EQ.'Y'.OR.YES.EQ.'y'.AND.(INT(Z).LE.NELEMENT)) THEN
         ZMU=MN(INT(Z))
      ELSE
         WRITE(iscw,'(1X,A)') 'Enter the mass of the atom'
         READ(iread,*) ZMU
      END IF
      RYDBRG = RY/(1.+ME/ZMU)
      IF (IPROB.NE.1) THEN
         IF (INT(Z).LT.10) THEN
            FZE=1.566432
         ELSEIF (INT(Z).LE.NELEMENT) THEN
            FZE=FZ(INT(Z-9))
         ELSE
            WRITE(iscw,'(1X,A)') 'Enter the f(Z) value in MHz'
            READ(iread,*) FZE
         ENDIF
       ENDIF
      ENDIF
      call MPI_BCAST(RYDBRG,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,
     $               ierr)
      call MPI_BCAST(ZMU,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,
     $               ierr)
      call MPI_BCAST(FZE,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ierr)

* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

       RMASS = ME/ZMU
*
* ****  INITIALIZE THE RADIAL GRID
*
      RHO=-4.0
      DO 10 J=1,NO
         R(J)=DEXP(RHO)/Z
         RR(J)=R(J)*R(J)
         R2(J)=DSQRT(R(J))
         RHO=RHO+H
10    CONTINUE
      RETURN
      END

ext: f


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1990  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ###############################################################
c     ##                                                           ##
c     ##  subroutine prtpdb  --  output of Protein Data Bank file  ##
c     ##                                                           ##
c     ###############################################################
c
c
c     "prtpdb" writes out a set of Protein Data Bank coordinates
c     to an external disk file
c
c
      subroutine prtpdb (ipdb)
      implicit none
      include 'sizes.i'
      include 'files.i'
      include 'pdb.i'
      include 'sequen.i'
      include 'titles.i'
      integer i,k,ipdb
      integer start,stop
      integer resmax,resnumb
      integer, allocatable :: resid(:)
      real*8 crdmin,crdmax
      logical opened
      logical rename
      logical reformat
      character*1 chnname
      character*1, allocatable :: chain(:)
      character*2 atmc,resc
      character*3 resname
      character*6 crdc
      character*38 fstr
      character*120 pdbfile
c
c
c     set flags for residue naming and large value formatting
c
      rename = .false.
      reformat = .false.
c
c     open the output unit if not already done
c
      inquire (unit=ipdb,opened=opened)
      if (.not. opened) then
         pdbfile = filename(1:leng)//'.pdb'
         call version (pdbfile,'new')
         open (unit=ipdb,file=pdbfile,status='new')
      end if
c
c     write out the header lines and the title
c
      if (ltitle .eq. 0) then
         fstr = '(''HEADER'',/,''COMPND'',/,''SOURCE'')'
         write (ipdb,fstr(1:32))
      else
         fstr = '(''HEADER'',4x,a,/,''COMPND'',/,''SOURCE'')'
         write (ipdb,fstr(1:37))  title(1:ltitle)
      end if
c
c     perform dynamic allocation of some local arrays
c
      allocate (resid(maxres))
      allocate (chain(maxres))
c
c     find the chain name and chain position for each residue
c
      do i = 1, nchain
         start = ichain(1,i)
         stop = ichain(2,i)
         do k = start, stop
            resid(k) = k - start + 1
            chain(k) = chnnam(i)
         end do
      end do
c
c     change some TINKER residue names to match PDB standards
c
      if (rename) then
         do i = 1, npdb
            if (resnam(i) .eq. 'CYX')  resnam(i) = 'CYS'
            if (resnam(i) .eq. 'CYD')  resnam(i) = 'CYS'
            if (resnam(i) .eq. 'TYD')  resnam(i) = 'TYR'
            if (resnam(i) .eq. 'HID')  resnam(i) = 'HIS'
            if (resnam(i) .eq. 'HIE')  resnam(i) = 'HIS'
            if (resnam(i) .eq. 'HIP')  resnam(i) = 'HIS'
            if (resnam(i) .eq. 'ASH')  resnam(i) = 'ASP'
            if (resnam(i) .eq. 'GLH')  resnam(i) = 'GLU'
            if (resnam(i) .eq. 'LYD')  resnam(i) = 'LYS'
         end do
      end if
c
c     set formatting to match the PDB fixed format standard
c
      atmc = 'i5'
      resc = 'i4'
      crdc = '3f8.3 '
c
c     check for large values requiring extended formatting
c
      resmax = 0
      crdmin = 0.0d0
      crdmax = 0.0d0
      do i = 1, npdb
         if (pdbtyp(i) .eq. 'ATOM  ') then
            resmax = max(resmax,resid(resnum(i)))
         else
            resmax = max(resmax,resnum(i))
         end if
         crdmin = min(crdmin,xpdb(i),ypdb(i),zpdb(i))
         crdmax = max(crdmax,xpdb(i),ypdb(i),zpdb(i))
      end do
      if (reformat) then
         if (npdb .ge. 100000)  atmc = 'i6'
         if (resmax .ge. 10000)  resc = 'i5'
         if (resmax .ge. 100000)  resc = 'i6'
         if (crdmin .le. -100.0d0)  crdc = '3f9.3 '
         if (crdmax .ge. 1000.0d0)  crdc = '3f9.3 '
         if (crdmin .le. -1000.0d0)  crdc = '3f10.3'
         if (crdmax .ge. 10000.0d0)  crdc = '3f10.3'
      end if
c
c     write info and coordinates for each PDB atom
c
      fstr = '(a6,'//atmc//',1x,a4,1x,a3,1x,a1,'//resc//
     &          ',4x,'//crdc//')'
      do i = 1, npdb
         resname = resnam(i)
         if (resname(2:3) .eq. '  ')  resname = '  '//resname(1:1)
         if (resname(3:3) .eq. ' ')  resname = ' '//resname(1:2)
         if (pdbtyp(i) .eq. 'ATOM  ') then
            resnumb = resid(resnum(i))
            chnname = chain(resnum(i))
         else
            resnumb = resnum(i)
            chnname = ' '
         end if
         write (ipdb,fstr)  pdbtyp(i),i,atmnam(i),resname,chnname,
     &                      resnumb,xpdb(i),ypdb(i),zpdb(i)
      end do
c
c     perform deallocation of some local arrays
c
      deallocate (resid)
      deallocate (chain)
c
c     write any connectivity records for PDB atoms
c
      fstr = '(''CONECT'',9'//atmc//')'
      do i = 1, npdb
         if (npdb12(i) .ne. 0) then
            write (ipdb,fstr(1:14))  i,(ipdb12(k,i),k=1,npdb12(i))
         end if
      end do
      fstr = '(''END'')'
      write (ipdb,fstr(1:7))
c
c     close the output unit if opened by this routine
c
c     if (.not. opened)  close (unit=ipdb)
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
      SUBROUTINE LINE(DX,DY)
C
C          ------------------------------------------------
C          ROUTINE NO. (  76)   VERSION (A8.4)    24:NOV:86
C          ------------------------------------------------
C
C          THIS DRAWS AN INCREMENTAL LINE OF LENGTH (DX,DY).
C
C
      COMMON /T0PPOS/ XPLOT0,YPLOT0
      COMMON /T3LIMS/ IMAXI,RMAXI,RMINI
C
C
      IF (ABS(DX).LT.RMINI.AND.ABS(DY).LT.RMINI) RETURN
C
      X= XPLOT0+DX
      Y= YPLOT0+DY
      CALL JOIN(X,Y)
C
      RETURN
      END

max_stars_repo_name: ZedThree/GHOST

max_issues_repo_name: ZedThree/GHOST

max_forks_repo_name: ZedThree/GHOST


text:
C
C $Id: mdrgfa.f,v 1.10 2008-09-18 00:42:18 kennison Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE MDRGFA (XCRA,YCRA,NCRA,IAAI,IAGI,NGPS)
C
C This routine is called by ARSCAM, which has been called by MDRGSX,
C which has been called by MDRGSF.  If the polygon defined by the
C arguments is one that should be filled, the polygon is transformed
C from its position on the unit square to a 1-degree square on a
C lat/lon grid and then to its position on the map defined by the
C current state of EZMAP and then filled.
C
        REAL XCRA(NCRA),YCRA(NCRA)
        INTEGER IAAI(NGPS),IAGI(NGPS),NGPS
C
C Declare required EZMAP common blocks.
C
        COMMON /MAPRGD/  ICOL(5),ICSF(5),IDPF,LCRA,NILN,NILT,OLAT,OLON
        INTEGER          ICOL,ICSF,IDPF,LCRA,NILN,NILT
        REAL             OLAT,OLON
        SAVE   /MAPRGD/
C
C Declare local variables.
C
        REAL             UVAL,VVAL,XVPL,XVPR,YVPB,YVPT,XWDL,XWDR,YWDB,
     +                   YWDT
        INTEGER          LNLG
C
C Because of the way MDRGFA is called, we know that the polygon has but
C one area identifier, relative to edge group identifier 1; if that id
C is greater than 0 and the polygon has at least three points, then we
C transform and fill it.
C
        IF (NCRA.GT.2.AND.IAAI(1).GT.0) THEN
          IF (ICSF(IAAI(1)).GE.0) THEN
            CALL GETSET (XVPL,XVPR,YVPB,YVPT,XWDL,XWDR,YWDB,YWDT,LNLG)
            CALL MDPRS
            NCRM=0
            DO 101 I=1,NCRA
              CALL MAPTRN (OLAT+YCRA(I),OLON+XCRA(I),UVAL,VVAL)
              IF (UVAL.NE.1.E12) THEN
                NCRM=NCRM+1
                XCRA(NCRM)=UVAL
                YCRA(NCRM)=VVAL
              END IF
  101       CONTINUE
            IF (NCRM.GT.2) THEN
              CALL GSFACI (ICSF(IAAI(1)))
              CALL GFA    (NCRM,XCRA,YCRA)
            END IF
            CALL SET (XVPL,XVPR,YVPB,YVPT,XWDL,XWDR,YWDB,YWDT,LNLG)
          END IF
        END IF
C
C Done.
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
	SUBROUTINE GG_DLTN  ( lincol, lintyp, linwid, ifreq, dlltln,
     +			      aloc, ifrmat, iret )
C************************************************************************
C* GG_DLTN								*
C*									*
C* This subroutine draws latitude/longitude lines on the graphics	*
C* device.  								*
C*									*
C* GG_DLTN ( LINCOL, LINTYP, LINWID, IFREQ, DLLTLN, ALOC, IFRMAT, IRET )*
C*									*
C* Input parameters:							*
C*	LINCOL		INTEGER		Line color			*
C*	LINTYP		INTEGER		Line type			*
C*	LINWID		INTEGER		Line width			*
C*	IFREQ(2)	INTEGER		Frequency of labels		*
C*	DLLTLN(2)	REAL		Frequency of lat/lon lines	*
C*	ALOC(2)		REAL		Lat/lon label locations		*
C*	IFRMAT		INTEGER		Lat/lon format			*
C*					  1 = +/- value			*
C*					  2 = N,S,E,W addded, no '-'    *
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*				  	  0 = normal return		*
C*				 	-13 = lines not drawn		*
C**									*
C* Log:									*
C* S. Maxwell/GSC	 1/97	Taken from GG_LTLN			*
C* T. Lee/GSC		 9/97	Fixed label freq in graphic mode	*
C* A. Hardy/GSC		12/00   Added lat/lon label locations, format	*
C* T. Piper/GSC		 5/01	Added GSSMTH				*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	INTEGER		ifreq (*)
	REAL		dlltln(*), aloc(*)
C*
	REAL		xlb (LLAXIS), ylb (LLAXIS)
	LOGICAL		polar
C-----------------------------------------------------------------------
C*	Save current line and smooth information.
C
	CALL GQLINE  ( ilto, ilhwo, iwdo, iwdhwo, ier )
	CALL GQSMTH  ( ismtyp, dens, ier )
	CALL GSLINE  ( lintyp, 0, linwid, 0, ier )
	CALL GSSMTH  ( 0, 0., ier )
C
C*      Get mode to see if this is a graph or map.
C
	CALL GQMODE  ( mode, ier )
C
C*	Set grid line color and draw lines.
C
	CALL GSCOLR  ( lincol, ier )
	IF ( mode .eq. 1 )  THEN
	    CALL GDRGRD  ( dlltln(1), dlltln(2), ifreq, aloc(1),
     +                     aloc(2), ifrmat, ier )
	    IF  ( ier .ne. 0 )  THEN
		CALL ER_WMSG  ( 'GEMPLT', ier, ' ', iret )
		iret = -13
		CALL ER_WMSG  ( 'GG', iret, ' ', ier )
	    END IF
	ELSE IF  ( mode .eq. 2 )  THEN
	    CALL GQGRAF  ( ixtyp, iytyp, ysiz, xlm, ylm, xrm, ytm,
     +                     ier )
C
C*	    Check for special case of polar plots.
C
	    IF  ( ixtyp .eq. 5 )  THEN
		polar = .true.
		IF  ( ( xlm .eq. xrm ) .and.
     +                 ( NINT ( ylm ) .eq. 225 ) .and.
     +                 ( NINT ( ytm ) .eq. 45 ) )  THEN
		    xrm = xlm / SQRT ( 2. )
		    xlm = 0.
		    ylm = 0.
		    ytm = 360.
		END IF
	    ELSE
		polar = .false.
	    END IF
C
C*          Get number of points for labels.
C
	    nx = ABS ( (xrm - xlm) / dlltln (1) ) + 1
	    ny = ABS ( (ytm - ylm) / dlltln (2) ) + 1
	    xlb (1) = xlm
	    ylb (1) = ylm
C
C*	    Get multiplier and label values.
C
	    IF  ( xlm .gt. xrm )  THEN
		ixmul = -1
	    ELSE
		ixmul = 1
	    END IF
	    IF  ( ylm .gt. ytm )  THEN
		iymul = -1
	    ELSE
		iymul = 1
	    END IF
C*
	    DO  i = 2, nx
		xlb ( i ) = xlb ( i-1 ) + ( dlltln (1) * ixmul )
	    END DO
	    DO  i = 2, ny
		ylb ( i ) = ylb ( i-1 ) + ( dlltln (2) * iymul )
	    END DO
C
C*	    Draw the x axis.
C
	    CALL GDAXIS  ( 1, ylm, .true., ifreq (1), ifreq (1),
     +                     ifreq (1), 0, nx, xlb, ier )
	    IF  ( ier .ne. 0 ) CALL ER_WMSG ( 'GEMPLT', ier, ' ', ier1 )
C
C*	    Draw the y axis.
C
	    IF  ( polar )  xlm = xrm
	    CALL GDAXIS  ( 2, xlm, .true., ifreq (2), ifreq (2),
     +                     ifreq (2), 0, ny, ylb, ier )
	    IF  ( ier .ne. 0 ) CALL ER_WMSG ( 'GEMPLT', ier, ' ', ier1 )
	ELSE
	    CALL ER_WMSG  ( 'GEMPLT', -2, 'invalid mode', ier )
	END IF
C
C*	Restore current line and smooth type.
C
	CALL GSLINE  ( ilto, 0, iwdo, 0, ier )
	CALL GSSMTH  ( ismtyp, dens, ier )
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


! { dg-do run }
! PR 30865 - passing a subroutine optional argument to size(dim=...)
! used to segfault.
program main
  implicit none
  integer :: a(2,3)
  integer :: ires

  call checkv (ires, a)
  if (ires /= 6) STOP 1
  call checkv (ires, a, 1)
  if (ires /= 2) STOP 2
contains
  subroutine checkv(ires,a1,opt1)
    integer, intent(out) :: ires
    integer :: a1(:,:)
    integer, optional :: opt1

    ires = size (a1, dim=opt1)
  end subroutine checkv
end program main


text:
c-----------------------------------------------------------------------
      subroutine setupcomm(comm,newcomm,newcommg,path_in,session_in)
      include 'mpif.h'
      include 'SIZE'
      include 'PARALLEL'
      include 'TSTEP'
      include 'INPUT'

      integer comm, newcomm, newcommg
      character session_in*(*), path_in*(*)
      logical flag

      common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal

      integer nid_global_root(0:nsessmax-1)
      character*132 session_mult(0:nsessmax-1), path_mult(0:nsessmax-1)

      logical ifhigh
      logical mpi_is_initialized

      integer*8 ntags

      call mpi_initialized(mpi_is_initialized, ierr)
      if (.not.mpi_is_initialized) call mpi_init(ierr)

      call mpi_comm_dup(comm,newcommg,ierr)
      newcomm = newcommg
      nekcomm = newcommg

      call mpi_comm_size(nekcomm,np_global,ierr)
      call mpi_comm_rank(nekcomm,nid_global,ierr)

      ! check upper tag size limit
      call mpi_comm_get_attr(MPI_COMM_WORLD,MPI_TAG_UB,ntags,flag,ierr)
      if (ntags .lt. np_global) then
         if(nid_global.eq.0) write(6,*) 'ABORT: MPI_TAG_UB too small!'
         call exitt
      endif

      ! set defaults
      nid         = nid_global
      ifneknek    = .false.
      ifneknekc   = .false. ! session are uncoupled
      nsessions   = 1

      ierr = 0
      nlin = 0
      if (nid .eq. 0) then
         l = ltrunc(session_in,len(session_in))
         if (l .gt. 0) then
            call blank(session_mult(0),132)
            call chcopy(session_mult(0), session_in, l)
            l = ltrunc(path_in,len(path_in))
            call blank(path_mult(0)   ,132)
            call chcopy(path_mult(0), path_in, l)
         else
           !write(6,*) 'Reading session file ...'
           open (unit=8,file='SESSION.NAME',status='old',err=24)
 21        read (8,*,END=22)
           nlin = nlin + 1
           goto 21
 22        rewind(8)
           if (nlin.gt.2) read(8,*,err=24) nsessions
           if (nsessions.gt.1) read(8,*,err=24) ifneknekc
           do n=0,nsessions-1
              call blank(session_mult(n),132)
              call blank(path_mult(n)   ,132)
              read(8,11,err=24) session_mult(n)
              read(8,11,err=24) path_mult(n)
              if (nsessions.gt.1) read(8,*,err=24)  npsess(n)
           enddo
 11        format(a132)
           close(8)
         endif
         if (nsessions.gt.1)
     $     write(6,*) 'Number of sessions:',nsessions
         goto 23
 24      ierr = 1
      endif
 23   continue
      call err_chk(ierr,' Error while reading SESSION.NAME!$')

      call bcast(nsessions,ISIZE)
      if (nsessions .gt. nsessmax)
     &   call exitti('nsessmax in SIZE too low!$',nsessmax)
      if (nsessions .gt. 1) ifneknek = .true.

      call bcast(ifneknekc,LSIZE)
      do n = 0,nsessions-1
         call bcast(npsess(n),ISIZE)
         call bcast(session_mult(n),132*CSIZE)
         call bcast(path_mult(n),132*CSIZE)
      enddo

      ! single session run
      if (.not.ifneknek) then
         ifneknekc = .false.
         session   = session_mult(0)
         path      = path_mult(0)
         amgfile  = session
         return
      endif

c     Check if specified number of ranks in each session is consistent
c     with the total number of ranks
      npall=0
      do n=0,nsessions-1
         npall=npall+npsess(n)
      enddo
      if (npall.ne.np_global)
     &   call exitti('Number of ranks does not match!$',npall)

c     Assign key for splitting into multiple groups
      nid_global_root_next=0
      do n=0,nsessions-1
         nid_global_root(n)=nid_global_root_next
         nid_global_root_next=nid_global_root(n)+npsess(n)
         if (nid_global.ge.nid_global_root(n).and.
     &       nid_global.lt.nid_global_root_next) idsess = n
      enddo
      call mpi_comm_split(comm,idsess,nid,newcomm,ierr)

      session = session_mult(idsess)
      path    = path_mult   (idsess)

      if (ifneknekc) then
         if (nsessions.gt.2) call exitti(
     &     'More than 2 coupled sessions are currently not supported!$',
     $     nsessions)
      endif

      amgfile  = session

      return
      end
c---------------------------------------------------------------------
      subroutine iniproc()
      include 'SIZE'
      include 'PARALLEL'
      include 'INPUT'
      include 'mpif.h'

      common /nekmpi/ nid_,np_,nekcomm,nekgroup,nekreal

      logical flag
      integer*8 isize_mpi, lb

      nid  = mynode()
      nid_ = nid
      np   = numnodes()
      np_  = np

      nio = -1             ! Default io flag
      if (nid.eq.0) nio=0  ! Only node 0 writes

      if (nid.eq.nio) then
         if (ifneknek) then
           call set_stdout(' ',idsess)
         else
           call set_stdout(' ',-1)
         endif
      endif

      if (wdsize .eq. 4)
     $   call exitti('Single precision mode not supported!',wdsize)

      call MPI_Type_Get_Extent(MPI_DOUBLE_PRECISION,lb,isize_mpi,ierr)
      if (isize_mpi .ne. wdsize) then
         call exitti('MPI real size does not match$',isize_mpi)
      endif

      call MPI_Type_Get_Extent(MPI_INTEGER,lb,isize_mpi,ierr)
       if (isize_mpi .ne. isize) then
         call exitti('MPI integer size does not match$',isize_mpi)
      endif

      call MPI_Type_Get_Extent(MPI_INTEGER8,lb,isize_mpi,ierr)
       if (isize_mpi .ne. isize8) then
         call exitti('MPI integer8 size does not match$',isize_mpi)
      endif

      PID = 0
      NULLPID=0
      NODE0=0
      NODE= NID+1

C     Test timer accuracy
      edif = 0.0
      do i = 1,10
         e1 = dnekclock()
         e2 = dnekclock()
         edif = edif + e2-e1
      enddo
      edif = edif/10.

      call fgslib_crystal_setup(cr_h,nekcomm,np)  ! set cr handle to new instance
      return
      end
c-----------------------------------------------------------------------
      subroutine gop( x, w, op, n)

c     Global vector commutative operation

      include 'CTIMER'

      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal

      real x(n), w(n)
      character*3 op

      if (ifsync) call nekgsync()

#ifdef TIMER
      if (icalld.eq.0) then
        tgop =0.0d0
        ngop =0
        icalld=1
      endif
      ngop = ngop + 1
      etime1=dnekclock()
#endif
c
      if (op.eq.'+  ') then
      call mpi_allreduce(x,w,n,MPI_DOUBLE_PRECISION,mpi_sum,nekcomm,ie)
      elseif (op.EQ.'M  ') then
      call mpi_allreduce(x,w,n,MPI_DOUBLE_PRECISION,mpi_max,nekcomm,ie)
      elseif (op.EQ.'m  ') then
      call mpi_allreduce(x,w,n,MPI_DOUBLE_PRECISION,mpi_min,nekcomm,ie)
      elseif (op.EQ.'*  ') then
      call mpi_allreduce(x,w,n,MPI_DOUBLE_PRECISION,mpi_prod,nekcomm,ie)
      else
      write(6,*) nid,' OP ',op,' not supported.  ABORT in GOP.'
      call exitt
      endif

      call copy(x,w,n)

#ifdef TIMER
      tgop =tgop +(dnekclock()-etime1)
#endif

      return
      end
c-----------------------------------------------------------------------
      subroutine igop( x, w, op, n)
c
c     Global vector commutative operation
c
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal

      integer x(n), w(n)
      character*3 op

      if     (op.eq.'+  ') then
        call mpi_allreduce (x,w,n,mpi_integer,mpi_sum ,nekcomm,ierr)
      elseif (op.EQ.'M  ') then
        call mpi_allreduce (x,w,n,mpi_integer,mpi_max ,nekcomm,ierr)
      elseif (op.EQ.'m  ') then
        call mpi_allreduce (x,w,n,mpi_integer,mpi_min ,nekcomm,ierr)
      elseif (op.EQ.'*  ') then
        call mpi_allreduce (x,w,n,mpi_integer,mpi_prod,nekcomm,ierr)
      else
        write(6,*) nid,' OP ',op,' not supported.  ABORT in igop.'
        call exitt
      endif

      call icopy(x,w,n)

      return
      end
c-----------------------------------------------------------------------
      subroutine i8gop( x, w, op, n)
c
c     Global vector commutative operation
c
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal

      integer*8 x(n), w(n)
      character*3 op

      if     (op.eq.'+  ') then
        call mpi_allreduce (x,w,n,mpi_integer8,mpi_sum ,nekcomm,ierr)
      elseif (op.EQ.'M  ') then
        call mpi_allreduce (x,w,n,mpi_integer8,mpi_max ,nekcomm,ierr)
      elseif (op.EQ.'m  ') then
        call mpi_allreduce (x,w,n,mpi_integer8,mpi_min ,nekcomm,ierr)
      elseif (op.EQ.'*  ') then
        call mpi_allreduce (x,w,n,mpi_integer8,mpi_prod,nekcomm,ierr)
      else
        write(6,*) nid,' OP ',op,' not supported.  ABORT in igop.'
        call exitt
      endif

      call i8copy(x,w,n)

      return
      end
c-----------------------------------------------------------------------
      subroutine csend(mtype,buf,len,jnid,jpid)
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      real*4 buf(1)

      call mpi_send (buf,len,mpi_byte,jnid,mtype,nekcomm,ierr)

      return
      end
c-----------------------------------------------------------------------
      subroutine crecv(mtype,buf,lenm)
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      integer status(mpi_status_size)
C
      real*4 buf(1)
      len = lenm
      jnid = mpi_any_source

      call mpi_recv (buf,len,mpi_byte
     $              ,jnid,mtype,nekcomm,status,ierr)
c
      if (len.gt.lenm) then
          write(6,*) nid,'long message in mpi_crecv:',len,lenm
          call exitt
      endif
c
      return
      end
c-----------------------------------------------------------------------
      subroutine crecv2(mtype,buf,lenm,jnid)
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      integer status(mpi_status_size)
C
      real*4 buf(1)
      len = lenm

      call mpi_recv (buf,len,mpi_byte
     $              ,jnid,mtype,nekcomm,status,ierr)
c
      if (len.gt.lenm) then
          write(6,*) nid,'long message in mpi_crecv:',len,lenm
          call exitt
      endif
c
      return
      end
c-----------------------------------------------------------------------
      subroutine crecv3(mtype,buf,len,lenm)
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      integer status(mpi_status_size)
C
      real*4 buf(1)
      len = lenm
      jnid = mpi_any_source

      call mpi_recv (buf,len,mpi_byte
     $            ,jnid,mtype,nekcomm,status,ierr)
      call mpi_get_count (status,mpi_byte,len,ierr)
c
      if (len.gt.lenm) then
          write(6,*) nid,'long message in mpi_crecv:',len,lenm
          call exitt
      endif
c
      return
      end
c-----------------------------------------------------------------------
      integer function numnodes()
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal

      call mpi_comm_size (nekcomm, numnodes , ierr)

      return
      end
c-----------------------------------------------------------------------
      integer function mynode()
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      integer myid

      call mpi_comm_rank (nekcomm, myid, ierr)
      mynode = myid

      return
      end
c-----------------------------------------------------------------------
      real*8 function dnekclock()
      include 'mpif.h'
c
      dnekclock=mpi_wtime()
c
      return
      end
c-----------------------------------------------------------------------
      real*8 function dnekclock_sync()
      include 'mpif.h'
c
      call nekgsync()
      dnekclock_sync=mpi_wtime()
c
      return
      end
c-----------------------------------------------------------------------
      subroutine lbcast(ifif)
C
C     Broadcast logical variable to all processors.
C
      include 'SIZE'
      include 'PARALLEL'
      include 'mpif.h'

      logical ifif

      if (np.eq.1) return

      item=0
      if (ifif) item=1
      call bcast(item,isize)
      ifif=.false.
      if (item.eq.1) ifif=.true.

      return
      end
c-----------------------------------------------------------------------
      subroutine bcast(buf,len)
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      real*4 buf(1)

      call mpi_bcast (buf,len,mpi_byte,0,nekcomm,ierr)

      return
      end
c-----------------------------------------------------------------------
      subroutine create_comm(inewcomm)
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal

c      call mpi_comm_group (mpi_comm_world,itmp,ierr)
c      call mpi_comm_create (mpi_comm_world,itmp,icomm,ierr)
c      call mpi_group_free (itmp,ierr)

      call mpi_comm_dup(nekcomm,inewcomm,ierr)

      return
      end
c-----------------------------------------------------------------------
      function isend(msgtag,x,len,jnid,jpid)
c
c     Note: len in bytes
c
      integer x(1)
C
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
C
      call mpi_isend (x,len,mpi_byte,jnid,msgtag
     $       ,nekcomm,imsg,ierr)
      isend = imsg
c     write(6,*) nid,' isend:',imsg,msgtag,len,jnid,(x(k),k=1,len/4)
c
      return
      end
c-----------------------------------------------------------------------
      function irecv(msgtag,x,len)
c
c     Note: len in bytes
c
      integer x(1)
C
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
C
      call mpi_irecv (x,len,mpi_byte,mpi_any_source,msgtag
     $       ,nekcomm,imsg,ierr)
      irecv = imsg
c     write(6,*) nid,' irecv:',imsg,msgtag,len
c
c
      return
      end
c-----------------------------------------------------------------------
      subroutine msgwait(imsg)
c
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      integer status(mpi_status_size)
c
c     write(6,*) nid,' msgwait:',imsg
c
      call mpi_wait (imsg,status,ierr)
c
      return
      end
c-----------------------------------------------------------------------
      subroutine nekgsync()

      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal

      call mpi_barrier(nekcomm,ierr)

      return
      end
c-----------------------------------------------------------------------
      subroutine exittr(stringi,rdata,idata)
      character*1 stringi(132)
      character*1 stringo(132)
      character*25 s25
      include 'SIZE'
      include 'TOTAL'
      include 'CTIMER'

      call blank(stringo,132)
      call chcopy(stringo,stringi,132)
      len = indx1(stringo,'$',1)
      write(s25,25) rdata,idata
   25 format(1x,1p1e14.6,i10)
      call chcopy(stringo(len),s25,25)

      if (nid.eq.0) write(6,1) (stringo(k),k=1,len+24)
    1 format('EXIT: ',132a1)

      call exitt

      return
      end
c-----------------------------------------------------------------------
      subroutine exitti(stringi,idata)
      character*1 stringi(132)
      character*1 stringo(132)
      character*11 s11
      include 'SIZE'
      include 'TOTAL'
      include 'CTIMER'

      call blank(stringo,132)
      call chcopy(stringo,stringi,132)
      len = indx1(stringo,'$',1)
      write(s11,11) idata
   11 format(1x,i10)
      call chcopy(stringo(len),s11,11)

      if (nid.eq.0) write(6,1) (stringo(k),k=1,len+10)
    1 format('EXIT: ',132a1)

      call exitt

      return
      end
c-----------------------------------------------------------------------
      subroutine err_chk(ierr,string)
      character*1 string(132)
      character*1 ostring(132)
      character*10 s10
      include 'SIZE'
c     include 'TOTAL'
c     include 'CTIMER'

      ierr = iglsum(ierr,1)
      if(ierr.eq.0) return

      len = indx1(string,'$',1)
      call blank(ostring,132)
      write(s10,11) ierr
   11 format(1x,' ierr=',i3)

      call chcopy(ostring,string,len-1)
      call chcopy(ostring(len),s10,10)

      if (nid.eq.0) write(6,1) (ostring(k),k=1,len+10)
    1 format('ERROR: ',132a1)

      call exitt

      return
      end
c
c-----------------------------------------------------------------------
      subroutine exitt0

      include 'SIZE'
      include 'TOTAL'

      if (nid.eq.0) then
         write(6,*) ' '
         write(6,'(A)') 'run successful: dying ...'
         write(6,*) ' '
      endif

c      if (nid.eq.0) call close_files
      call print_runtime_info
      call nek_die(0)

      return
      end
c-----------------------------------------------------------------------
      subroutine exitt

      include 'SIZE'
      include 'TOTAL'

      if (nid.eq.0) then
         write(6,*) ' '
         write(6,'(A)') 'an error occured: dying ...'
         write(6,*) ' '
      endif

c      call print_stack()
c      if (nid.eq.0) call close_files
c      call print_runtime_info
      call nek_die(1)

      return
      end
c-----------------------------------------------------------------------
      subroutine print_runtime_info
      include 'SIZE'
      include 'TOTAL'
      include 'CTIMER'
      include 'mpif.h'

#ifdef PAPI
      gflops = glsum(dnekgflops(),1)
#endif

      tstop  = dnekclock_sync()
      ttotal = tstop-etimes
      tsol   = max(ttime - tprep,0.0)
      nxyz   = lx1*ly1*lz1

      dtmp4 = glsum(getmaxrss(),1)/1e9

      if (nid.eq.0) then
         dtmp1 = 0
         dtmp2 = 0
         if(istep.gt.0) then
           dgp   = nvtot
           dgp   = max(dgp,1.)*max(istep,1)
           dtmp0 = np*(ttime-tprep)
           dtmp1 = 0
           if (dtmp0.gt.0) dtmp1 = dgp/dtmp0
           dtmp2 = (ttime-tprep)/max(istep,1)
         endif
         write(6,*) ' '
         write(6,'(5(A,1p1e13.5,A,/))')
     &       'total elapsed time             : ',ttotal, ' sec'
     &      ,'total solver time w/o IO       : ',tsol,   ' sec'
     &      ,'time/timestep                  : ',dtmp2 , ' sec'
     &      ,'avg throughput per timestep    : ',dtmp1 , ' gridpts/CPUs'
     &      ,'total max memory usage         : ',dtmp4 , ' GB'
#ifdef PAPI
         write(6,'(1(A,1p1e13.5,/))')
     &      ,'total Gflops/s                 : ',gflops
#endif
      endif
      call flush_io

      return
      end
c-----------------------------------------------------------------------
      subroutine nek_die(ierr)
      include 'SIZE'
      include 'mpif.h'

      call mpi_finalize (ierr_)
      call cexit(ierr)

      return
      end
c-----------------------------------------------------------------------
      subroutine fgslib_userExitHandler(istatus)

      call exitt

      return
      end
c-----------------------------------------------------------------------
      subroutine printHeader

      include 'SIZE'
      include 'PARALLEL'

      include 'HEADER'
      write(6,*) 'Number of MPI ranks :', np
c      WRITE(6,*) 'REAL     wdsize     :',WDSIZE
c      WRITE(6,*) 'INTEGER  wdsize     :',ISIZE
c      WRITE(6,*) 'INTEGER8 wdsize     :',ISIZE8
      WRITE(6,*) ' '

      return
      end
c-----------------------------------------------------------------------
      function igl_running_sum(in)
c
      include 'mpif.h'
      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      integer status(mpi_status_size)
      integer x,w,r

      x = in  ! running sum
      w = in  ! working buff
      r = 0   ! recv buff

      call mpi_scan(x,r,1,mpi_integer,mpi_sum,nekcomm,ierr)
      igl_running_sum = r

      return
      end
c-----------------------------------------------------------------------
      subroutine platform_timer(ivb) ! mxm, ping-pong, and all_reduce timer

      include 'SIZE'
      include 'TOTAL'

      call mxm_test_all(nid,ivb)  ! measure mxm times
c     call exitti('done mxm_test_all$',ivb)

      call comm_test(ivb)         ! measure message-passing and all-reduce times

      return
      end
c-----------------------------------------------------------------------
      subroutine comm_test(ivb) ! measure message-passing and all-reduce times
                                ! ivb = 0 --> minimal verbosity
                                ! ivb = 1 --> fully verbose
                                ! ivb = 2 --> smaller sample set(shorter)

      include 'SIZE'
      include 'PARALLEL'

      call gop_test(ivb)   ! added, Jan. 8, 2008

      log_np=log2(np)
      np2 = 2**log_np
      if (np2.eq.np) call gp2_test(ivb)   ! added, Jan. 8, 2008

      io = 6
      n512 = min(512,np-1)

      do nodeb=1,n512
         call pingpongo(alphas,betas,0,nodeb,.0005,io,ivb)
         if (nid.eq.0) write(6,2) nodeb,np,alphas,betas
    2    format(2i10,1p2e15.7,' alpha betao')
      enddo

      do kk=0,2
      do nodeb=1,n512
         call pingpong (alphas,betas,0,nodeb,.0005,io,ivb,kk)
         if (nid.eq.0) write(6,1) nodeb,np,alphas,betas,kk
    1    format(2i10,1p2e15.7,' alpha beta',i1)
      enddo
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine pingpong(alphas,betas,nodea,nodeb,dt,io,ivb,kk)

      include 'SIZE'
      common /nekmpi/ mid,np,nekcomm,nekgroup,nekreal

      parameter  (lt=lx1*ly1*lz1*lelt)
      parameter (mwd = 3*lt/2)
      common /scrns/ x(mwd),y(mwd),x1(mwd),y1(mwd)

      include 'mpif.h'
      integer status(mpi_status_size)

      character*10 fname

      if (nid.eq.nodea) then
         write(fname,3) np,nodeb
    3    format('t',i4.4,'.',i4.4)
         if (io.ne.6) open (unit=io,file=fname)
      endif

      call nekgsync
      call get_msg_vol(msg_vol,dt,nodea,nodeb) ! Est. msg vol for dt s

      nwds = 0
      if (nid.eq.nodea.and.ivb.gt.0) write(io,*)

      betas = 0  ! Reported inverse bandwidth
      count = 0

      do itest = 1,500

         nloop = msg_vol/(nwds+2)
         nloop = min(nloop,1000)
         nloop = max(nloop,1)

         len   = 8*nwds

         if (kk.eq.0)
     $      call ping_loop (t1,t0,len,nloop,nodea,nodeb,nid,x,y,x1,y1)
         if (kk.eq.1)
     $      call ping_loop1(t1,t0,len,nloop,nodea,nodeb,nid,x,y)
         if (kk.eq.2)
     $      call ping_loop2(t1,t0,len,nloop,nodea,nodeb,nid,x,y)

         if (nid.eq.nodea) then
            tmsg = (t1-t0)/(2*nloop)   ! 2*nloop--> Double Buffer
            tmsg = tmsg / 2.           ! one-way cost = 1/2 round-trip
            tpwd = tmsg                ! time-per-word
            if (nwds.gt.0) tpwd = tmsg/nwds
            if (ivb.gt.0) write(io,1) nodeb,np,nloop,nwds,tmsg,tpwd,kk
    1       format(3i6,i12,1p2e16.8,' pgn',i1)

            if (nwds.eq.1) then
               alphas = tmsg
            elseif (nwds.gt.10000) then   ! "average" beta
               betas = (betas*count + tpwd)/(count+1)
               count = count + 1
            endif
         endif

         if (ivb.eq.2) then
            nwds = (nwds+1)*1.25
         else
            nwds = (nwds+1)*1.016
         endif
         if (nwds.gt.mwd) then
c        if (nwds.gt.1024) then
            if (nid.eq.nodea.and.io.ne.6) close(unit=io)
            call nekgsync
            return
         endif

      enddo

      if (nid.eq.nodea.and.io.ne.6) close(unit=io)
      call nekgsync

      return
      end
c-----------------------------------------------------------------------
      subroutine pingpongo(alphas,betas,nodea,nodeb,dt,io,ivb)

      include 'SIZE'
      common /nekmpi/ mid,np,nekcomm,nekgroup,nekreal

      parameter  (lt=lx1*ly1*lz1*lelt)
      parameter (mwd = 3*lt)
      common /scrns/ x(mwd),y(mwd)

      include 'mpif.h'
      integer status(mpi_status_size)

      character*10 fname

      if (nid.eq.nodea) then
         write(fname,3) np,nodeb
    3    format('t',i4.4,'.',i4.4)
         if (io.ne.6) open (unit=io,file=fname)
      endif

      call nekgsync
      call get_msg_vol(msg_vol,dt,nodea,nodeb) ! Est. msg vol for dt s

      nwds = 0
      if (nid.eq.nodea.and.ivb.gt.0) write(io,*)

      betas = 0  ! Reported inverse bandwidth
      count = 0

      do itest = 1,500
         call nekgsync
         nloop = msg_vol/(nwds+2)
         nloop = min(nloop,1000)
         nloop = max(nloop,1)

         len   = 8*nwds
         jnid = mpi_any_source

         if (nid.eq.nodea) then

            msg  = irecv(itest,y,1)
            call csend(itest,x,1,nodeb,0)   ! Initiate send, to synch.
            call msgwait(msg)

            t0 = mpi_wtime ()
            do i=1,nloop
               call mpi_irecv(y,len,mpi_byte,mpi_any_source,i
     $                        ,nekcomm,msg,ierr)
               call mpi_send (x,len,mpi_byte,nodeb,i,nekcomm,ierr)
               call mpi_wait (msg,status,ierr)
            enddo
            t1 = mpi_wtime ()
            tmsg = (t1-t0)/nloop
            tmsg = tmsg / 2.       ! Round-trip message time = twice one-way
            tpwd = tmsg
            if (nwds.gt.0) tpwd = tmsg/nwds
            if (ivb.gt.0) write(io,1) nodeb,np,nloop,nwds,tmsg,tpwd
    1       format(3i6,i12,1p2e16.8,' pgo')

            if (nwds.eq.1) then
               alphas = tmsg
            elseif (nwds.gt.10000) then
               betas = (betas*count + tpwd)/(count+1)
               count = count + 1
            endif

         elseif (nid.eq.nodeb) then

            call crecv(itest,y,1)           ! Initiate send, to synch.
            call csend(itest,x,1,nodea,0)

            t0 = dnekclock()
            do i=1,nloop
               call mpi_recv (y,len,mpi_byte
     $               ,jnid,i,nekcomm,status,ierr)
               call mpi_send (x,len,mpi_byte,nodea,i,nekcomm,ierr)
            enddo
            t1 = dnekclock()
            tmsg = (t1-t0)/nloop

         endif

         nwds = (nwds+1)*1.016
         if (nwds.gt.mwd) then
            if (nid.eq.nodea.and.io.ne.6) close(unit=io)
            call nekgsync
            return
         endif

      enddo

      if (nid.eq.nodea.and.io.ne.6) close(unit=io)
      call nekgsync

      return
      end
c-----------------------------------------------------------------------
      subroutine get_msg_vol(msg_vol,dt,nodea,nodeb)
      include 'SIZE'
      common /nekmpi/ mid,np,nekcomm,nekgroup,nekreal
      parameter (lt=lx1*ly1*lz1*lelt)
      common /scrns/ x(3*lt),y(3*lt)
!
!     Est. msg vol for dt s
!
      msg_vol = 1000

      nwds  = min(1000,lt)
      nloop = 50

      tmsg = 0.
      call gop(tmsg,t1,'+  ',1)

      len = 8*nwds
      if (nid.eq.nodea) then

         msg  = irecv(1,y,1)
         call csend(1,x,1,nodeb,0)   ! Initiate send, to synch.
         call msgwait(msg)

         t0 = dnekclock()
         do i=1,nloop
            msg  = irecv(i,y,len)
            call csend(i,x,len,nodeb,0)
            call msgwait(msg)
         enddo
         t1   = dnekclock()
         tmsg = (t1-t0)/nloop
         tpwd = tmsg/nwds

      elseif (nid.eq.nodeb) then

         call crecv(1,y,1)           ! Initiate send, to synch.
         call csend(1,x,1,nodea,0)

         t0 = dnekclock()
         do i=1,nloop
            call crecv(i,y,len)
            call csend(i,x,len,nodea,0)
         enddo
         t1   = dnekclock()
         tmsg = (t1-t0)/nloop
         tmsg = 0.

      endif

      call gop(tmsg,t1,'+  ',1)
      msg_vol = nwds*(dt/tmsg)
c     if (nid.eq.nodea) write(6,*) nid,msg_vol,nwds,dt,tmsg,' msgvol'

      return
      end
c-----------------------------------------------------------------------
      subroutine gop_test(ivb)
      include 'SIZE'
      common /nekmpi/ mid,np,nekcomm,nekgroup,nekreal
      include 'mpif.h'
      integer status(mpi_status_size)

      parameter  (lt=lx1*ly1*lz1*lelt)
      parameter (mwd = 3*lt)
      common /scrns/ x(mwd),y(mwd)
      common /scruz/ times(2,500)
      common /scrcg/ nwd(500)

      nwds  = 1
      mtest = 0
      do itest = 1,500
         nwds = (nwds+1)*1.016
         if (nwds.gt.mwd) goto 100
         mtest = mtest+1
         nwd(mtest) = nwds
      enddo
  100 continue

      nwds = 1
      do itest = mtest,1,-1

         tiny = 1.e-27
         call cfill(x,tiny,mwd)
         nwds = nwd(itest)
         call nekgsync

         t0 = mpi_wtime ()
         call gop(x,y,'+  ',nwds)
         call gop(x,y,'+  ',nwds)
         call gop(x,y,'+  ',nwds)
         call gop(x,y,'+  ',nwds)
         call gop(x,y,'+  ',nwds)
         call gop(x,y,'+  ',nwds)
         t1 = mpi_wtime ()

         tmsg = (t1-t0)/6 ! six calls
         tpwd = tmsg
         if (nwds.gt.0) tpwd = tmsg/nwds
         times(1,itest) = tmsg
         times(2,itest) = tpwd

      enddo
  101 continue

      if (nid.eq.0) then
         nwds = 1
         do itest=1,500
            if (ivb.gt.0.or.itest.eq.1)
     $         write(6,1) np,nwds,(times(k,itest),k=1,2)
    1       format(i12,i12,1p2e16.8,' gop')
            nwds = (nwds+1)*1.016
            if (nwds.gt.mwd) goto 102
         enddo
  102    continue
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine gp2_test(ivb)

      include 'SIZE'
      include 'mpif.h'

      common /nekmpi/ mid,np,nekcomm,nekgroup,nekreal
      integer status(mpi_status_size)

      parameter  (lt=lx1*ly1*lz1*lelt)
      parameter (mwd = 3*lt)
      common /scrns/ x(mwd),y(mwd)
      common /scruz/ times(2,500)

      call rzero(x,mwd)

      nwds = 1
      do itest = 1,500
         call gp2(x,y,'+  ',1,nid,np)

         t0 = mpi_wtime ()
         call gp2(x,y,'+  ',nwds,nid,np)
         call gp2(x,y,'+  ',nwds,nid,np)
         call gp2(x,y,'+  ',nwds,nid,np)
         call gp2(x,y,'+  ',nwds,nid,np)
         t1 = mpi_wtime ()

         tmsg = (t1-t0)/4 ! four calls
         tpwd = tmsg
         if (nwds.gt.0) tpwd = tmsg/nwds
         times(1,itest) = tmsg
         times(2,itest) = tpwd

         nwds = (nwds+1)*1.016
         if (nwds.gt.mwd) goto 101
      enddo
  101 continue

      if (nid.eq.0) then
         nwds = 1
         do itest=1,500
            if (ivb.gt.0.or.itest.eq.1)
     $         write(6,1) np,nwds,(times(k,itest),k=1,2)
    1       format(i12,i12,1p2e16.8,' gp2')
            nwds = (nwds+1)*1.016
            if (nwds.gt.mwd) goto 102
         enddo
  102    continue
      endif

      return
      end
c-----------------------------------------------------------------------
      integer function xor(m,n)
c
c  If NOT running on a parallel processor, it is sufficient to
c  have this routine return a value of XOR=1.
c
c  Pick one of the following:
c
c  UNIX 4.2, f77:
       XOR = OR(M,N)-AND(M,N)
c
c  Intel FTN286:
c     XOR = M.NEQV.N
c
c  Ryan-McFarland Fortran
C      XOR = IEOR(M,N)
c
c     XOR = 0
c     IF(M.EQ.1 .OR.  N.EQ.1) XOR=1
c     IF(M.EQ.0 .AND. N.EQ.0) XOR=0
c     IF(M.EQ.1 .AND. N.EQ.1) XOR=0
c     IF(M.GT.1 .OR.N.GT.1 .OR.M.LT.0.OR.N.LT.0) THEN
c        PRINT*,'ERROR IN XOR'
c        STOP
c     ENDIF
C
      return
      end
c-----------------------------------------------------------------------
      subroutine gp2( x, w, op, n, nid, np)
c
c     Global vector commutative operation using spanning tree.
c
c     Std. fan-in/fan-out

      real x(n), w(n)
      character*3 op

      integer bit, bytes, cnt, diff, spsize, i,
     *   parent, troot, xor, root, lnp, log2
      logical ifgot

      integer type
      save    type
      data    type  /998/

      type  = type+100
      if (type.gt.9992) type=type-998
      typer = type-1
      bytes = 8*n

      root    = 0
      troot   = max0((nid/np)*np, root)
      diff    = xor(nid,troot)
      nullpid = 0

c     Accumulate contributions from children, if any
      level2=1
    5 continue
         level=level2
         level2=level+level
         if (mod(nid,level2).ne.0) goto 20
            call crecv(type,w,bytes)
            if (op.eq.'+  ') then
               do i=1,n
                  x(i) = x(i) + w(i)
               enddo
            elseif (op.eq.'*  ') then
               do i=1,n
                  x(i) = x(i) * w(i)
               enddo
            elseif (op.eq.'M  ') then
               do i=1,n
                  x(i) = max(x(i),w(i))
               enddo
            elseif (op.eq.'m  ') then
               do i=1,n
                  x(i) = min(x(i),w(i))
               enddo
            endif
         if (level2.lt.np) goto 5

c     Pass result back to parent
   20 parent = nid-level
      if (nid .ne. 0) call csend(type,x,bytes,parent,nullpid)

c     Await final answer from node 0 via log_2 fan out
      level=np/2
      ifgot=.false.
      if (nid.eq.root) ifgot=.true.

      lnp = log2(np)
      do i=1,lnp
        if (ifgot) then
           jnid=nid+level
           call csend(typer,x,bytes,jnid,nullpid)
        elseif (mod(nid,level).eq.0) then
           call crecv(typer,x,bytes)
           ifgot=.true.
        endif
        level=level/2
      enddo

      return
      end
c-----------------------------------------------------------------------
      subroutine ping_loop1(t1,t0,len,nloop,nodea,nodeb,nid,x,y)

      common /nekmpi/ mid,np,nekcomm,nekgroup,nekreal

      real x(1),y(1)

      include 'mpif.h'
      integer status(mpi_status_size)

      i=0
      if (nid.eq.nodea) then
         call nekgsync
         call mpi_irecv(y,len,mpi_byte,nodeb,i,nekcomm,msg,ierr)    ! 1b
         call mpi_send (x,len,mpi_byte,nodeb,i,nekcomm,ierr)        ! 1a
c        call mpi_rsend(x,len,mpi_byte,nodeb,i,nekcomm,ierr)        ! 1a
         call msgwait(msg)                                          ! 1b

         t0 = mpi_wtime ()
         do i=1,nloop
            call mpi_irecv(y,len,mpi_byte,nodeb,i,nekcomm,msg,ierr) ! 2b
            call mpi_send (x,len,mpi_byte,nodeb,i,nekcomm,ierr)     ! 2a
c           call mpi_rsend(x,len,mpi_byte,nodeb,i,nekcomm,ierr)     ! 2a
            call mpi_wait (msg,status,ierr)                         ! 2b
         enddo
         t1 = mpi_wtime ()

      elseif (nid.eq.nodeb) then

         call mpi_irecv(y,len,mpi_byte,nodea,i,nekcomm,msg,ierr)    ! 1a
         call nekgsync
         call mpi_wait (msg,status,ierr)                            ! 1a

         j=i
         do i=1,nloop
            call mpi_irecv(y,len,mpi_byte,nodea,i,nekcomm,msg,ierr) ! 2a
c           call mpi_rsend(x,len,mpi_byte,nodea,j,nekcomm,ierr)     ! 1b
            call mpi_send (x,len,mpi_byte,nodea,j,nekcomm,ierr)     ! 1b
            call mpi_wait (msg,status,ierr)                         ! 2a
            j=i
         enddo
c        call mpi_rsend(x,len,mpi_byte,nodea,j,nekcomm,ierr)        ! nb
         call mpi_send (x,len,mpi_byte,nodea,j,nekcomm,ierr)        ! nb

      else
         call nekgsync
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine ping_loop2(t1,t0,len,nloop,nodea,nodeb,nid,x,y)

      common /nekmpi/ mid,np,nekcomm,nekgroup,nekreal

      real x(1),y(1)

      include 'mpif.h'
      integer status(mpi_status_size)

      i=0
      if (nid.eq.nodea) then
         call nekgsync
         call mpi_irecv(y,len,mpi_byte,nodeb,i,nekcomm,msg,ierr)    ! 1b
         call mpi_send (x,len,mpi_byte,nodeb,i,nekcomm,ierr)        ! 1a
         call msgwait(msg)                                          ! 1b

         t0 = mpi_wtime ()
         do i=1,nloop
            call mpi_send (x,len,mpi_byte,nodeb,i,nekcomm,ierr)     ! 2a
            call mpi_irecv(y,len,mpi_byte,nodeb,i,nekcomm,msg,ierr) ! 2b
            call mpi_wait (msg,status,ierr)                         ! 2b
         enddo
         t1 = mpi_wtime ()

      elseif (nid.eq.nodeb) then

         call mpi_irecv(y,len,mpi_byte,nodea,i,nekcomm,msg,ierr)    ! 1a
         call nekgsync
         call mpi_wait (msg,status,ierr)                            ! 1a

         j=i
         do i=1,nloop
            call mpi_send (x,len,mpi_byte,nodea,j,nekcomm,ierr)     ! 1b
            call mpi_irecv(y,len,mpi_byte,nodea,i,nekcomm,msg,ierr) ! 2a
            call mpi_wait (msg,status,ierr)                         ! 2a
            j=i
         enddo
         call mpi_send (x,len,mpi_byte,nodea,j,nekcomm,ierr)        ! nb

      else
         call nekgsync
      endif

      return
      end
c-----------------------------------------------------------------------
      subroutine ping_loop(t1,t0,len,nloop,nodea,nodeb,nid,x1,y1,x2,y2)
c     Double Buffer : does 2*nloop timings

      common /nekmpi/ mid,np,nekcomm,nekgroup,nekreal

      real x1(1),y1(1),x2(1),y2(1)

      include 'mpif.h'
      integer status(mpi_status_size)

      itag=1
      if (nid.eq.nodea) then
         call mpi_irecv(y1,len,mpi_byte,nodeb,itag,nekcomm,msg1,ierr)   ! 1b
         call nekgsync

         t0 = mpi_wtime ()
         do i=1,nloop
            call mpi_send (x1,len,mpi_byte,nodeb,itag,nekcomm,ierr)     ! 1a
            call mpi_irecv(y2,len,mpi_byte,nodeb,itag,nekcomm,msg2,ierr)! 2b
            call mpi_wait (msg1,status,ierr)                            ! 1b
            call mpi_send (x2,len,mpi_byte,nodeb,itag,nekcomm,ierr)     ! 2a
            call mpi_irecv(y1,len,mpi_byte,nodeb,itag,nekcomm,msg1,ierr)! 3b
            call mpi_wait (msg2,status,ierr)                            ! 2b
         enddo
         t1 = mpi_wtime ()
         call mpi_send (x1,len,mpi_byte,nodeb,itag,nekcomm,ierr)        ! nb
         call mpi_wait (msg1,status,ierr)                              ! nb

      elseif (nid.eq.nodeb) then

         call mpi_irecv(y1,len,mpi_byte,nodea,itag,nekcomm,msg1,ierr)   ! nb
         call nekgsync

         do i=1,nloop
            call mpi_wait (msg1,status,ierr)                            ! 1a
            call mpi_send (x1,len,mpi_byte,nodea,itag,nekcomm,ierr)     ! 1b
            call mpi_irecv(y2,len,mpi_byte,nodea,itag,nekcomm,msg2,ierr)! 2a
            call mpi_wait (msg2,status,ierr)                            ! 2a
            call mpi_send (x2,len,mpi_byte,nodea,itag,nekcomm,ierr)     ! 2b
            call mpi_irecv(y1,len,mpi_byte,nodea,itag,nekcomm,msg1,ierr)! 3a
         enddo
         call mpi_wait (msg1,status,ierr)                            ! 2a
         call mpi_send (x1,len,mpi_byte,nodea,itag,nekcomm,ierr)        ! nb

      else
         call nekgsync
      endif

      return
      end
c-----------------------------------------------------------------------
      integer*8 function i8gl_running_sum(in)
c
      include 'mpif.h'

      integer*8 in

      common /nekmpi/ nid,np,nekcomm,nekgroup,nekreal
      integer status(mpi_status_size)
      integer*8 x,r

      x = in  ! running sum
      r = 0   ! recv buff

      call mpi_scan(x,r,1,mpi_integer8,mpi_sum,nekcomm,ierr)
      i8gl_running_sum = r

      return
      end
c-----------------------------------------------------------------------
      subroutine close_files
      logical ifopen

      do ii=1,99
         if(ii.ne.5.and.ii.ne.6) then
           inquire(unit=ii,opened=ifopen)
           if(ifopen) close(ii)
         endif
      enddo

      return
      end
c----------------------------------------------------------------------
      subroutine neknekgsync()

      include 'SIZE'
      include 'PARALLEL'

      call mpi_barrier(iglobalcomm,ierr)

      return
      end
c------------------------------------------------------------------------
      subroutine setnekcomm(comm_in)

      include 'SIZE'

      integer comm_in
      common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal

      nekcomm = comm_in
      call mpi_comm_size(nekcomm,mp,ierr)
      np = mp

      return
      end

max_stars_repo_name: RonRahaman/nekRS

max_issues_repo_name: neams-th-coe/nekRS

max_forks_repo_name: neams-th-coe/nekRS


text:
!** Copyright (c) 1989, 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.

!* Tests for runtime library MATMUL routines

program p

  parameter(NbrTests=2112)
  parameter(o_extent=2)
  parameter(n_extent=6)
  parameter(m_extent=4)
  parameter(k_extent=8)

  integer*4, dimension(n_extent,m_extent) :: arr1
  integer*4, dimension(m_extent,k_extent) :: arr2
  integer*4, dimension(n_extent,k_extent) :: arr3
  integer*4, dimension(n_extent,m_extent,o_extent) :: arr4
  integer*4, dimension(n_extent,o_extent,m_extent) :: arr5
  integer*4, dimension(o_extent,n_extent,m_extent) :: arr6

  integer*4, dimension(o_extent,m_extent,k_extent) :: arr7
  integer*4, dimension(m_extent,o_extent,k_extent) :: arr8
  integer*4, dimension(m_extent,k_extent,o_extent) :: arr9

  integer*4, dimension(n_extent,k_extent,o_extent) :: arr10
  integer*4, dimension(n_extent,o_extent,k_extent) :: arr11
  integer*4, dimension(o_extent,n_extent,k_extent) :: arr12

  integer*4, dimension(2:n_extent+1,m_extent) :: arr13
  integer*4, dimension(2:m_extent+1,k_extent) :: arr14
  integer*4, dimension(2:n_extent+1,k_extent) :: arr15
  integer*4, dimension(n_extent,2:m_extent+1) :: arr16
  integer*4, dimension(m_extent,2:k_extent+1) :: arr17
  integer*4, dimension(n_extent,2:k_extent+1) :: arr18
  integer*4, dimension(n_extent,k_extent) :: arr19

  integer*4 :: expect(NbrTests)
  integer*4 :: results(NbrTests)

  integer:: i,j

  data arr1 /0,1,2,3,4,5,6,7,8,9,10,11,			&
             12,13,14,15,16,17,18,19,20,22,22,23/
  data arr2 /0,1,2,3,4,5,6,7,8,9,10,11,			&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31/
  data arr3 /0,1,2,3,4,5,6,7,8,9,10,11, 			&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47/
  data arr4 /0,1,2,3,4,5,6,7,8,9,10,11, 			&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47/
  data arr5 /0,1,2,3,4,5,6,7,8,9,10,11, 			&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47/
  data arr6 /0,1,2,3,4,5,6,7,8,9,10,11, 			&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47/
  data arr7 /0,1,2,3,4,5,6,7,8,9,10,11, 			&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47,		&
             48,49,50,51,52,53,54,55,56,57,58,59,		&
             60,61,62,63/
  data arr8 /0,1,2,3,4,5,6,7,8,9,10,11, 			&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47,		&
             48,49,50,51,52,53,54,55,56,57,58,59,		&
             60,61,62,63/
  data arr9 /0,1,2,3,4,5,6,7,8,9,10,11, 			&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47,		&
             48,49,50,51,52,53,54,55,56,57,58,59,		&
             60,61,62,63/
  data arr10 /0,1,2,3,4,5,6,7,8,9,10,11,	 		&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47,		&
             48,49,50,51,52,53,54,55,56,57,58,59,		&
             60,61,62,63,64,65,66,67,68,69,70,71,		&
             72,73,74,75,76,77,78,79,80,81,82,83,		&
             84,85,86,87,88,89,90,91,92,93,94,95/
  data arr11 /0,1,2,3,4,5,6,7,8,9,10,11,	 		&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47,		&
             48,49,50,51,52,53,54,55,56,57,58,59,		&
             60,61,62,63,64,65,66,67,68,69,70,71,		&
             72,73,74,75,76,77,78,79,80,81,82,83,		&
             84,85,86,87,88,89,90,91,92,93,94,95/
  data arr12 /0,1,2,3,4,5,6,7,8,9,10,11,	 		&
             12,13,14,15,16,17,18,19,20,21,22,23,		&
             24,25,26,27,28,29,30,31,32,33,34,35,		&
             36,37,38,39,40,41,42,43,44,45,46,47,		&
             48,49,50,51,52,53,54,55,56,57,58,59,		&
             60,61,62,63,64,65,66,67,68,69,70,71,		&
             72,73,74,75,76,77,78,79,80,81,82,83,		&
             84,85,86,87,88,89,90,91,92,93,94,95/
  data arr13 /0,1,2,3,4,5,6,7,8,9,10,11,                        &
             12,13,14,15,16,17,18,19,20,22,22,23/
  data arr14 /0,1,2,3,4,5,6,7,8,9,10,11,                        &
             12,13,14,15,16,17,18,19,20,21,22,23,               &
             24,25,26,27,28,29,30,31/
  data arr15 /0,1,2,3,4,5,6,7,8,9,10,11,                        &
             12,13,14,15,16,17,18,19,20,21,22,23,               &
             24,25,26,27,28,29,30,31,32,33,34,35,               &
             36,37,38,39,40,41,42,43,44,45,46,47/
  data arr16 /0,1,2,3,4,5,6,7,8,9,10,11,                        &
             12,13,14,15,16,17,18,19,20,22,22,23/
  data arr17 /0,1,2,3,4,5,6,7,8,9,10,11,                        &
             12,13,14,15,16,17,18,19,20,21,22,23,               &
             24,25,26,27,28,29,30,31/
  data arr18 /0,1,2,3,4,5,6,7,8,9,10,11,                        &
             12,13,14,15,16,17,18,19,20,21,22,23,               &
             24,25,26,27,28,29,30,31,32,33,34,35,               &
             36,37,38,39,40,41,42,43,44,45,46,47/

  data expect /  &
  ! test 1-48
      84, 90, 96, 105, 108, 114, 228, 250, 272, &
      301, 316, 338, 372, 410, 448, 497, 524, 562, &
      516, 570, 624, 693, 732, 786, 660, 730, 800, &
      889, 940, 1010, 804, 890, 976, 1085, 1148, 1234, &
      948, 1050, 1152, 1281, 1356, 1458, 1092, 1210, 1328, &
      1477, 1564, 1682, &
  ! test 49-96
      0, 90, 96, 105, 108, 114, 0, 250, 272, &
      301, 316, 338, 0, 410, 448, 497, 524, 562, &
      0, 570, 624, 693, 732, 786, 0, 730, 800, &
      889, 940, 1010, 0, 890, 976, 1085, 1148, 1234, &
      0, 1050, 1152, 1281, 1356, 1458, 0, 1210, 1328, &
      1477, 1564, 1682, &
  ! test 97-144
      84, 90, 96, 105, 108, 0, 228, 250, 272, &
      301, 316, 0, 372, 410, 448, 497, 524, 0, &
      516, 570, 624, 693, 732, 0, 660, 730, 800, &
      889, 940, 0, 804, 890, 976, 1085, 1148, 0, &
      948, 1050, 1152, 1281, 1356, 0, 1092, 1210, 1328, &
      1477, 1564, 0, &
  ! test 145-192
      84, 90, 96, 105, 108, 114, 228, 246, 264, &
      289, 300, 318, 372, 402, 432, 473, 492, 522, &
      516, 558, 600, 657, 684, 726, 660, 714, 768, &
      841, 876, 930, 804, 870, 936, 1025, 1068, 1134, &
      948, 1026, 1104, 1209, 1260, 1338, 1092, 1182, 1272, &
      1393, 1452, 1542, &
  ! test 193-240
      84, 90, 96, 105, 108, 114, 228, 246, 264, &
      289, 300, 318, 372, 402, 432, 473, 492, 522, &
      516, 558, 600, 657, 684, 726, 660, 714, 768, &
      841, 876, 930, 804, 870, 936, 1025, 1068, 1134, &
      948, 1026, 1104, 1209, 1260, 1338, 1092, 1182, 1272, &
      1393, 1452, 1542, &
  ! test 241-288
      30, 33, 36, 39, 42, 45, 102, 117, 132, &
      147, 162, 177, 174, 201, 228, 255, 282, 309, &
      246, 285, 324, 363, 402, 441, 318, 369, 420, &
      471, 522, 573, 390, 453, 516, 579, 642, 705, &
      462, 537, 612, 687, 762, 837, 534, 621, 708, &
      795, 882, 969, &
  !test 289-336
      30, 33, 36, 0, 0, 0, 102, 117, 132, &
      0, 0, 0, 174, 201, 228, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, &
  ! test 337-384
      0, 0, 0, 0, 0, 0, 0, 246, 264, &
      289, 0, 0, 0, 402, 432, 473, 0, 0, &
      0, 558, 600, 657, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, &
  ! test 385-432
      30, 33, 36, 39, 42, 45, 0, 0, 0, &
      0, 0, 0, 174, 201, 228, 255, 282, 309, &
      0, 0, 0, 0, 0, 0, 318, 369, 420, &
      471, 522, 573, 0, 0, 0, 0, 0, 0, &
      462, 537, 612, 687, 762, 837, 0, 0, 0, &
      0, 0, 0, &
  ! test 433-480
      48, 0, 54, 0, 60, 0, 192, 0, 222, &
      0, 252, 0, 336, 0, 390, 0, 444, 0, &
      480, 0, 558, 0, 636, 0, 624, 0, 726, &
      0, 828, 0, 768, 0, 894, 0, 1020, 0, &
      912, 0, 1062, 0, 1212, 0, 1056, 0, 1230, &
      0, 1404, 0, &
  ! test 481-528
      30, 0, 36, 0, 42, 0, 0, 0, 0, &
      0, 0, 0, 174, 0, 228, 0, 282, 0, &
      0, 0, 0, 0, 0, 0, 318, 0, 420, &
      0, 522, 0, 0, 0, 0, 0, 0, 0, &
      462, 0, 612, 0, 762, 0, 0, 0, 0, &
      0, 0, 0, &
  ! test 529-576
      48, 0, 54, 0, 60, 0, 0, 0, 0, &
      0, 0, 0, 336, 0, 390, 0, 444, 0, &
      0, 0, 0, 0, 0, 0, 624, 0, 726, &
      0, 828, 0, 0, 0, 0, 0, 0, 0, &
      912, 0, 1062, 0, 1212, 0, 0, 0, 0, &
      0, 0, 0, &
  ! test 577-624
      0, 0, 0, 0, 0, 0, 0, 138, 0, &
      174, 0, 210, 0, 0, 0, 0, 0, 0, &
      0, 306, 0, 390, 0, 474, 0, 0, 0, &
      0, 0, 0, 0, 474, 0, 606, 0, 738, &
      0, 0, 0, 0, 0, 0, 0, 642, 0, &
      822, 0, 1002, &
  !test 625-672
      0, 621, 0, 795, 0, 969, 0, 0, 0, &
      0, 0, 0, 0, 453, 0, 579, 0, 705, &
      0, 0, 0, 0, 0, 0, 0, 285, 0, &
      363, 0, 441, 0, 0, 0, 0, 0, 0, &
      0, 117, 0, 147, 0, 177, 0, 0, 0, &
      0, 0, 0, &
  ! test 673-720
      912, 0, 1062, 0, 1212, 0, 0, 0, 0, &
      0, 0, 0, 624, 0, 726, 0, 828, 0, &
      0, 0, 0, 0, 0, 0, 336, 0, 390, &
      0, 444, 0, 0, 0, 0, 0, 0, 0, &
      48, 0, 54, 0, 60, 0, 0, 0, 0, &
      0, 0, 0, &
  ! test 721-768
      0, 0, 0, 0, 0, 0, 0, 138, 0, &
      174, 0, 210, 0, 0, 0, 0, 0, 0, &
      0, 306, 0, 390, 0, 474, 0, 0, 0, &
      0, 0, 0, 0, 474, 0, 606, 0, 738, &
      0, 0, 0, 0, 0, 0, 0, 642, 0, &
      822, 0, 1002, &
  ! test 767-864
      0, 0, 0, 0, 0, 0, 0, 492, 528, &
      564, 0, 0, 0, 804, 864, 924, 0, 0, &
      0, 1116, 1200, 1284, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, &
  ! test 865-960
      102, 105, 108, 111, 114, 117, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 1542, 1593, 1644, &
      1695, 1746, 1797, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 2982, 3081, 3180, 3279, 3378, 3477, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      4422, 4569, 4716, 4863, 5010, 5157, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, &
  !test 961-1056
      0, 3576, 0, 0, 0, 3774, 0, 0, 0, &
      3972, 0, 0, 0, 4008, 0, 0, 0, 4230, &
      0, 0, 0, 4452, 0, 0, 0, 4440, 0, &
      0, 0, 4686, 0, 0, 0, 4932, 0, 0, &
      0, 4872, 0, 0, 0, 5142, 0, 0, 0, &
      5412, 0, 0, 0, 5304, 0, 0, 0, 5598, &
      0, 0, 0, 5892, 0, 0, 0, 5736, 0, &
      0, 0, 6054, 0, 0, 0, 6372, 0, 0, &
      0, 6168, 0, 0, 0, 6510, 0, 0, 0, &
      6852, 0, 0, 0, 6600, 0, 0, 0, 6966, &
      0, 0, 0, 7332, 0, 0, &
  ! test 1057-1152
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 294, 0, 324, 0, 354, 0, &
      0, 0, 0, 0, 0, 0, 1158, 0, 1284, &
      0, 1410, 0, 0, 0, 0, 0, 0, 0, &
      2022, 0, 2244, 0, 2466, 0, 0, 0, 0, &
      0, 0, 0, 2886, 0, 3204, 0, 3522, 0, &
      0, 0, 0, 0, 0, 0, &
  !test 1153-1248
      0, 0, 0, 0, 0, 0, 96, 0, 102, &
      0, 108, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 672, 0, 726, 0, 780, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      1248, 0, 1350, 0, 1452, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 1824, 0, 1974, &
      0, 2124, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, &
  ! test 1249-1344
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 789, 0, 0, 0, &
      867, 0, 0, 0, 945, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 1701, 0, 0, 0, 1875, 0, 0, &
      0, 2049, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 2613, &
      0, 0, 0, 2883, 0, 0, 0, 3153, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 3525, 0, 0, 0, &
      3891, 0, 0, 0, 4257, 0, &
  !test 1355-1440
      0, 2769, 0, 3501, 0, 4233, 0, 0, 0, &
      0, 0, 0, 0, 2409, 0, 3045, 0, 3681, &
      0, 0, 0, 0, 0, 0, 0, 2049, 0, &
      2589, 0, 3129, 0, 0, 0, 0, 0, 0, &
      0, 1689, 0, 2133, 0, 2577, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, &
  !test 1441-1536
      0, 0, 0, 0, 0, 0, 3720, 0, 4332, &
      0, 4944, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 2568, 0, 2988, 0, 3408, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      1416, 0, 1644, 0, 1872, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 264, 0, 300, &
      0, 336, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, &
  !test 1537-1632
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 474, 0, 0, &
      0, 594, 0, 0, 0, 714, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 1194, 0, 0, 0, 1506, 0, &
      0, 0, 1818, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      1914, 0, 0, 0, 2418, 0, 0, 0, 2922, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 2634, 0, 0, &
      0, 3330, 0, 0, 0, 4026, &
  ! test 11633-1680
      84, 90, 96, 105, 108, 114, 228, 250, 272, &
      301, 316, 338, 372, 410, 448, 497, 524, 562, &
      516, 570, 624, 693, 732, 786, 660, 730, 800, &
      889, 940, 1010, 804, 890, 976, 1085, 1148, 1234, &
      948, 1050, 1152, 1281, 1356, 1458, 1092, 1210, 1328, &
      1477, 1564, 1682, &

  ! test 1681-1728
      84, 90, 96, 105, 108, 114, 228, 250, 272, &
      301, 316, 338, 372, 410, 448, 497, 524, 562, &
      516, 570, 624, 693, 732, 786, 660, 730, 800, &
      889, 940, 1010, 804, 890, 976, 1085, 1148, 1234, &
      948, 1050, 1152, 1281, 1356, 1458, 1092, 1210, 1328, &
      1477, 1564, 1682, &

  ! test 1729-1776
      84, 90, 96, 105, 108, 0, 228, 250, 272, &
      301, 316, 0, 372, 410, 448, 497, 524, 0, &
      516, 570, 624, 693, 732, 0, 660, 730, 800, &
      889, 940, 0, 804, 890, 976, 1085, 1148, 0, &
      948, 1050, 1152, 1281, 1356, 0, 1092, 1210, 1328, &
      1477, 1564, 0, &
  ! test 1777-1824
      84, 90, 96, 105, 108, 114, 228, 246, 264, &
      289, 300, 318, 372, 402, 432, 473, 492, 522, &
      516, 558, 600, 657, 684, 726, 660, 714, 768, &
      841, 876, 930, 804, 870, 936, 1025, 1068, 1134, &
      948, 1026, 1104, 1209, 1260, 1338, 1092, 1182, 1272, &
      1393, 1452, 1542, &
  ! test 1825-1872
      84, 90, 96, 105, 108, 114, 228, 246, 264, &
      289, 300, 318, 372, 402, 432, 473, 492, 522, &
      516, 558, 600, 657, 684, 726, 660, 714, 768, &
      841, 876, 930, 804, 870, 936, 1025, 1068, 1134, &
      948, 1026, 1104, 1209, 1260, 1338, 1092, 1182, 1272, &
      1393, 1452, 1542, &
  ! test 1873-1920
      30, 33, 36, 39, 42, 45, 102, 117, 132, &
      147, 162, 177, 174, 201, 228, 255, 282, 309, &
      246, 285, 324, 363, 402, 441, 318, 369, 420, &
      471, 522, 573, 390, 453, 516, 579, 642, 705, &
      462, 537, 612, 687, 762, 837, 534, 621, 708, &
      795, 882, 969, &
  !test 1921-1968
      30, 33, 36, 0, 0, 0, 102, 117, 132, &
      0, 0, 0, 174, 201, 228, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, &
  ! test 1969-2016
      0, 0, 0, 0, 0, 0, 0, 246, 264, &
      289, 0, 0, 0, 402, 432, 473, 0, 0, &
      0, 558, 600, 657, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, 0, 0, 0, 0, 0, 0, &
      0, 0, 0, &
  ! test 2017-2064
      30, 33, 36, 39, 42, 45, 0, 0, 0, &
      0, 0, 0, 174, 201, 228, 255, 282, 309, &
      0, 0, 0, 0, 0, 0, 318, 369, 420, &
      471, 522, 573, 0, 0, 0, 0, 0, 0, &
      462, 537, 612, 687, 762, 837, 0, 0, 0, &
      0, 0, 0, &
  ! test 2065-2112
      48, 0, 54, 0, 60, 0, 192, 0, 222, &
      0, 252, 0, 336, 0, 390, 0, 444, 0, &
      480, 0, 558, 0, 636, 0, 624, 0, 726, &
      0, 828, 0, 768, 0, 894, 0, 1020, 0, &
      912, 0, 1062, 0, 1212, 0, 1056, 0, 1230, &
      0, 1404, 0 /

  !test 1-48
  arr3=0
  arr3 = matmul(arr1,arr2)
  call assign_result(1,48,arr3,results)
  ! print *,"**** arr3 = matmul(arr1,arr2) ***"
  ! print *,arr3

  ! test 49-96
  arr3=0
  arr3(2:n_extent,:) = matmul(arr1(2:n_extent,:),arr2)
  call assign_result(49,96,arr3,results)
  ! print *,"**** arr3(2:n_extent,:) = matmul(arr1(2:n_extent,:),arr2) ***"
  ! print *,arr3

  ! test 97-144
  arr3=0
  arr3(1:n_extent-1,:) = matmul(arr1(1:n_extent-1,:),arr2)
  call assign_result(97,144,arr3,results)
  ! print *,"**** arr3(1:n_extent-1,:) = matmul(arr1(1:n_extent-1,:),arr2) ***"
  ! print *,arr3

  ! test 145-192
  arr3=0
  arr3 = matmul(arr1(:,2:m_extent),arr2(2:m_extent,:))
  call assign_result(145,192,arr3,results)
  ! print *,"**** arr3 = matmul(arr1(:,2:m_extent),arr2(2:m_extent,:)) ***"
  ! print *,arr3

  ! test 193-240
  arr3=0
  arr3 = matmul(arr1(:,2:m_extent),arr2(2:m_extent,:))
  call assign_result(193,240,arr3,results)
  ! print *,"**** arr3 = matmul(arr1(:,2:m_extent),arr2(2:m_extent,:)) ***"
  ! print *,arr3

  ! test 241-288
  arr3=0
  arr3 = matmul(arr1(:,1:m_extent-1),arr2(1:m_extent-1,:))
  call assign_result(241,288,arr3,results)
  ! print *,"**** arr3 = matmul(arr1(:,1:m_extent-1),arr2(1:m_extent-1,:)) ***"
  ! print *,arr3

  !test 289-336
  arr3=0
  arr3(1:3,1:3) = matmul(arr1(1:3,1:3),arr2(1:3,1:3))
  call assign_result(289,336,arr3,results)
  ! print *,"**** arr3(1:3,1:3) = matmul(arr1(1:3,1:3),arr2(1:3,1:3)) ***"
  ! print *,arr3

  ! test 337-384
  arr3=0
  arr3(2:4,2:4) = matmul(arr1(2:4,2:4),arr2(2:4,2:4))
  call assign_result(337,384,arr3,results)
  ! print *,"**** arr3(2:4,2:4) = matmul(arr1(2:4,2:4),arr2(2:4,2:4)) ***"
  ! print *,arr3

  ! test 385-432
  arr3=0
  arr3(:,1:k_extent:2) = matmul(arr1(:,1:m_extent-1),arr2(1:m_extent-1,1:k_extent:2))
  call assign_result(385,432,arr3,results)
  ! print *,"**** arr3(:,1:k_extent:2) = matmul(arr1(:,1:m_extent-1),arr2(1:m_extent-1,1:k_extent:2)) ***"
  ! print *,arr3

  ! test 433-480
  arr3=0
  arr3(1:n_extent:2,:) = matmul(arr1(1:n_extent:2,2:m_extent),arr2(1:m_extent-1,:))
  call assign_result(433,480,arr3,results)
  ! print *,"**** arr3(1:n_extent:2,:) = matmul(arr1(1:n_extent:2,2:m_extent),arr2(1:m_extent-1,:)) ***"
  ! print *,arr3

  ! test 481-528
  arr3=0
  arr3(1:n_extent:2,1:k_extent:2) = matmul(arr1(1:n_extent:2,1:m_extent-1),      &
                                           arr2(1:m_extent-1,1:k_extent:2))
  call assign_result(481,528,arr3,results)
  ! print *,"**** arr3(1:n_extent:2,1:k_extent:2) = matmul(arr1(1:n_extent:2,1:m_extent-1),arr2(1:m_extent-1,1:k_extent:2)) ***"
  ! print *,arr3

  ! test 529-576
  arr3=0
  arr3(1:n_extent-1:2,1:k_extent-1:2) = matmul(arr1(1:n_extent-1:2,2:m_extent),	&
                                               arr2(1:m_extent-1,1:k_extent:2))
  call assign_result(529,576,arr3,results)
  ! print *,"**** arr3(1:n_extent-1:2,1:k_extent-1:2) = matmul(arr1(1:n_extent-1:2,2:m_extent),arr2(1:m_extent-1,1:k_extent:2))***"
  ! print *,arr3

  ! test 577-624
  arr3=0
  arr3(2:n_extent:2,2:k_extent:2) = matmul(arr1(2:n_extent:2,1:m_extent-1),	&
                                               arr2(2:m_extent,2:k_extent:2))
  call assign_result(577,624,arr3,results)
  ! print *,"**** arr3(2:n_extent:2,2:k_extent:2) = matmul(arr1(2:n_extent:2,1:m_extent-1),arr2(2:m_extent,2:k_extent:2)) ***"
  ! print *,arr3

  !test 625-672
  arr3=0
  arr3(n_extent:1:-2,1:k_extent:2) = matmul(arr1(n_extent:1:-2,1:m_extent-1),      &
                                           arr2(1:m_extent-1,k_extent:1:-2))
  call assign_result(625,672,arr3,results)
  ! print *,"**** arr3(n_extent:1:-2,1:k_extent:2) = matmul(arr1(n_extent:1:-2,1:m_extent-1), arr2(1:m_extent-1,k_extent:1:-2)) ***"
  ! print *,arr3

  ! test 673-720
  arr3=0
  arr3(1:n_extent-1:2,k_extent-1:1:-2) = matmul(arr1(1:n_extent-1:2,m_extent:2:-1),	&
                                               arr2(m_extent-1:1:-1,1:k_extent:2))
  call assign_result(673,720,arr3,results)
  ! print *,"**** arr3(1:n_extent-1:2,k_extent-1:1:-2) = matmul(arr1(1:n_extent-1:2,m_extent:2:-1),arr2(m_extent-1:1:-1,1:k_extent:2)) ***"
  ! print *,arr3

  ! test 721-768
  arr3=0
  arr3(n_extent:2:-2,k_extent:2:-2) = matmul(arr1(n_extent:2:-2,m_extent-1:1:-1),	&
                                               arr2(m_extent:2:-1,k_extent:2:-2))
  call assign_result(721,768,arr3,results)
  ! print *,"**** arr3(n_extent:2:-2,k_extent:2:-2) = matmul(arr1(n_extent:2:-2,m_extent-1:1:-1),arr2(m_extent:2:-1,k_extent:2:-2)) ***"
  ! print *,arr3

  ! test 769-864
  arr10=0
  arr10(2:4,2:4:1) = matmul(arr4(2:4,2:4,1),arr7(1,2:4,2:4))
  call assign_result(769,864,arr10,results)
  ! print *,"**** arr10(2:4,2:4:1) = matmul(arr4(2:4,2:4,1),arr7(1,2:4,2:4)) ***"
  ! print *,arr10

  ! test 865-960
  arr11=0
  arr11(:,1,1:k_extent:2) = matmul(arr4(:,1:m_extent-1,2),arr8(1:m_extent-1,1,1:k_extent:2))
  call assign_result(865,960,arr11,results)
  ! print *,"**** arr11(:,1,1:k_extent:2) = matmul(arr4(:,1:m_extent-1,2),arr8(1:m_extent-1,1,1:k_extent:2)) ***"
  ! print *,arr11

  !test 961-1056
  arr12=0
  arr12(2,1:n_extent:2,:) = matmul(arr4(1:n_extent:2,2:m_extent,2),arr9(1:m_extent-1,:,2))
  call assign_result(961,1056,arr12,results)
  ! print *,"**** arr12(2,1:n_extent:2,:) = matmul(arr4(1:n_extent:2,2:m_extent,2),arr9(1:m_extent-1,:,2)) ***"
  ! print *,arr12

  ! test 1057-1152
  arr10=0
  arr10(1:n_extent:2,1:k_extent:2,2) = matmul(arr5(1:n_extent:2,2,1:m_extent-1),      &
                                           arr8(1:m_extent-1,2,1:k_extent:2))
  call assign_result(1057,1152,arr10,results)
  ! print *,"**** arr10(1:n_extent:2,1:k_extent:2,2) = matmul(arr5(1:n_extent:2,2,1:m_extent-1), arr8(1:m_extent-1,2,1:k_extent:2)) ***"
  ! print *,arr10

  !test 1153-1248
  arr11=0
  arr11(1:n_extent-1:2,2,1:k_extent-1:2) = matmul(arr5(1:n_extent-1:2,1,2:m_extent),	&
                                               arr9(1:m_extent-1,1:k_extent:2,1))
  call assign_result(1153,1248,arr11,results)
  ! print *,"**** arr11(1:n_extent-1:2,2,1:k_extent-1:2) = matmul(arr5(1:n_extent-1:2,1,2:m_extent),arr9(1:m_extent-1,1:k_extent:2,1)) ***"
  ! print *,arr11

  ! test 1249-1344
  arr12=0
  arr12(1,2:n_extent:2,2:k_extent:2) = matmul(arr5(2:n_extent:2,2,1:m_extent-1),	&
                                               arr7(2,2:m_extent,2:k_extent:2))
  call assign_result(1249,1344,arr12,results)
  ! print *,"**** arr12(1,2:n_extent:2,2:k_extent:2) = matmul(arr5(2:n_extent:2,2,1:m_extent-1), arr7(2,2:m_extent,2:k_extent:2)) ***"
  ! print *,arr12

  !test 1345-1440
  arr10=0
  arr10(n_extent:1:-2,1:k_extent:2,1) = matmul(arr6(2,n_extent:1:-2,1:m_extent-1),      &
                                           arr9(1:m_extent-1,k_extent:1:-2,2))
  call assign_result(1345,1440,arr10,results)
  ! print *,"**** arr10(n_extent:1:-2,1:k_extent:2,1) = matmul(arr6(2,n_extent:1:-2,1:m_extent-1),arr9(1:m_extent-1,k_extent:1:-2,2)) ***"
  ! print *,arr10

  !test 1441-1536
  arr11=0
  arr11(1:n_extent-1:2,2,k_extent-1:1:-2) = matmul(arr6(1,1:n_extent-1:2,m_extent:2:-1),	&
                                               arr7(2,m_extent-1:1:-1,1:k_extent:2))
  call assign_result(1441,1536,arr11,results)
  ! print *,"**** arr11(1:n_extent-1:2,2,k_extent-1:1:-2) = matmul(arr6(1,1:n_extent-1:2,m_extent:2:-1),arr7(2,m_extent-1:1:-1,1:k_extent:2)) ***"
  ! print *,arr11

  !test 1537-1632
  arr12=0
  arr12(2,n_extent:2:-2,k_extent:2:-2) = matmul(arr6(2,n_extent:2:-2,m_extent-1:1:-1),	&
                                               arr8(m_extent:2:-1,1,k_extent:2:-2))
  call assign_result(1537,1632,arr12,results)
  ! print *,"**** arr12(2,n_extent:2:-2,k_extent:2:-2) = matmul(arr6(2,n_extent:2:-2,m_extent-1:1:-1),arr8(m_extent:2:-1,1,k_extent:2:-2)) ***"
  ! print *,arr12

  arr19 = 0

  ! test 1633-1680
  arr15=0
  arr15 =  arr19 + matmul(arr13,arr14)
  call assign_result(1633,1680,arr15,results)
  ! print *,arr15

  ! test 1681-1728
  arr18=0
  arr18 = arr19 +  matmul(arr16,arr17)
  call assign_result(1681,1728,arr18,results)
  ! print *,arr15

  ! test 1729-1776
  arr15=0
  arr15(2:n_extent,:) = matmul(arr13(2:n_extent,:),arr14)
  call assign_result(1729,1776,arr15,results)
  ! print *,arr15

  ! test 1777-1824
  arr15=0
  arr15 = arr19 + matmul(arr13(:,2:m_extent),arr14(3:m_extent+1,:))
  call assign_result(1777,1824,arr15,results)
  ! print *,arr15

  ! test 1825-1872
  arr15=0
  arr15 = arr19 + matmul(arr13(:,2:m_extent),arr14(3:m_extent+1,:))
  call assign_result(1825,1872,arr15,results)
  ! print *,arr15

  ! test 1873-1920
  arr15=0
  arr15 = arr19 + matmul(arr13(:,1:m_extent-1),arr14(2:m_extent,:))
  call assign_result(1873,1920,arr15,results)
  !print *,arr15

  !test 1921-1968
  arr15=0
  arr15(2:4,1:3) = arr19(2:4,1:3) + matmul(arr13(2:4,1:3),arr14(2:4,1:3))
  call assign_result(1921,1968,arr15,results)
  ! print *,arr15

  ! test 1969-2016
  arr15=0
  arr15(3:5,2:4) = arr19(3:5,2:4) + matmul(arr13(3:5,2:4),arr14(3:5,2:4))
  call assign_result(1969,2016,arr15,results)
  ! print *,arr15

  ! test 2017-2064
  arr15=0
  arr15(:,1:k_extent:2) = arr19(:,1:k_extent:2) +  &
                         matmul(arr13(:,1:m_extent-1),arr14(2:m_extent,1:k_extent:2))
  call assign_result(2017,2064,arr15,results)
  ! print *,arr15

  ! test 2065-2112
  arr15=0
  arr15(2:n_extent+1:2,:) = arr19(2:n_extent+1:2,:) +   &
                           matmul(arr13(2:n_extent+1:2,2:m_extent),arr14(2:m_extent,:))
  call assign_result(2065,2112,arr15,results)
  ! print *,arr15

  call check(results, expect, NbrTests)
end program

subroutine assign_result(s_idx, e_idx , arr, rslt)
  integer*4, dimension(1:e_idx-s_idx+1) :: arr
  integer*4, dimension(e_idx) :: rslt
  integer:: s_idx, e_idx

  rslt(s_idx:e_idx) = arr

end subroutine