text:
subroutine inidis(lundia    ,error     ,runid     ,cyclic    ,timnow    , &
                & itdis     ,itstrt    ,itfinish  ,sferic    ,grdang    , &
                & nsrc      ,nsrcd     ,lstsc     ,j         ,nmmaxj    , &
                & icx       ,icy       ,namsrc    ,disint    ,dismmt    , &
                & namcon    ,mnksrc    ,alfas     ,disch     , &
                & disch0    ,disch1    ,rint      ,rint0     ,rint1     , &
                & umdis     ,umdis0    ,umdis1    ,vmdis     ,vmdis0    , &
                & vmdis1    ,bubble    ,kmax      ,kspu      ,kspv      , &
                & upwsrc    ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: inidis.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/inichk/inidis.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Reads the time dependent data from file for the first time
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use globaldata
    use string_module
    use dfparall
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer  , pointer :: nxbub
    integer  , pointer :: itdate
    integer  , pointer :: lundis
    real(fp) , pointer :: tstop
    real(fp) , pointer :: dt
    real(fp) , pointer :: scalef
    real(fp) , pointer :: timscl
!
! Global variables
!
    integer                                                                     , intent(in)  :: icx      !!  Increment in the X-dir., if ICX= NMAX then computation proceeds in the X-dir. If icx=1 then computation proceeds in the Y-dir.
    integer                                                                     , intent(in)  :: icy      !  Increment in the Y-dir. (see ICX)
    integer                                                                     , intent(in)  :: itfinish !  Description and declaration in inttim.igs
    integer                                                                     , intent(in)  :: itstrt   !  Description and declaration in inttim.igs
    integer                                                                                :: j        !!  Begin pointer for arrays which have been transformed into 1D arrays. Due to the shift in the 2nd (M-) index, J = -2*NMAX + 1
    integer                                                                     , intent(in)  :: kmax     !  Description and declaration in esm_alloc_int.f90
    integer                                                                     , intent(in)  :: lstsc    !  Description and declaration in dimens.igs
    integer                                                                     , intent(in)  :: lundia   !  Description and declaration in inout.igs
    integer                                                                     , intent(in)  :: nmmaxj   !  Description and declaration in dimens.igs
    integer                                                                     , intent(in)  :: nsrc     !  Description and declaration in esm_alloc_int.f90
    integer                                                                     , intent(in)  :: nsrcd    !  Description and declaration in dimens.igs
    integer                                                                     , intent(in)  :: upwsrc   !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(5, nsrcd)                                          , intent(out) :: itdis    !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(7, nsrc)                                           , intent(in)  :: mnksrc   !  Description and declaration in esm_alloc_int.f90
    logical                                                                     , intent(in)  :: bubble   !  Description and declaration in procs.igs
    logical                                                                     , intent(in)  :: cyclic   !!  Flag = TRUE if cyclic system assumed
    logical                                                                     , intent(out) :: error    !!  Flag=TRUE if an erroris encountered
    logical                                                                     , intent(in)  :: sferic   !  Description and declaration in tricom.igs
    real(fp)                                                                                :: grdang   !  Description and declaration in tricom.igs
    real(fp)                                                                    , intent(in)  :: timnow
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)                             , intent(in)  :: alfas    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(lstsc, nsrcd)                                                    :: rint     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(lstsc, nsrcd)                                                    :: rint0    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(lstsc, nsrcd)                                                    :: rint1    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: disch    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: disch0   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: disch1   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: umdis    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: umdis0   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: umdis1   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: vmdis    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: vmdis0   !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(nsrcd)                                                           :: vmdis1   !  Description and declaration in esm_alloc_real.f90
    character(*)                                                                , intent(in)  :: runid
    character(1) , dimension(nsrcd)                                                           :: disint   !  Description and declaration in esm_alloc_char.f90
    character(1) , dimension(nsrcd)                                                           :: dismmt   !  Description and declaration in esm_alloc_char.f90
    character(20), dimension(lstsc)                                             , intent(in)  :: namcon   !  Description and declaration in esm_alloc_char.f90
    character(20), dimension(nsrc)                                                            :: namsrc   !  Description and declaration in esm_alloc_char.f90
    integer      , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax)                :: kspu     !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nlb:gdp%d%nub, gdp%d%mlb:gdp%d%mub, 0:kmax)                :: kspv     !  Description and declaration in esm_alloc_int.f90
!
!> Local variables
!
    integer                                    :: ddb
    integer                                    :: icxy        ! MAX value of ICX and ICY
    integer                                    :: iend
    integer                                    :: ifound
    integer                                    :: iocond
    integer                                    :: irecrd      ! Counter of records if input file is a direct access file
    integer, dimension(:), allocatable         :: irecs
    integer                                    :: isrc        ! Index number of discharge location
    integer                                    :: isrc_nodup  ! isrc, skipping duplicate discharges, introduced by bubble screens
    integer                                    :: istart
    integer                                    :: ja_upw      ! Default only around momentum discharge locations upwind is expected = 1 For no upwind JA_UPW = -1
    integer                                    :: l           ! Loop counter over LSTSC
    integer                                    :: lrec        ! Record length of direct access file
    integer                                    :: lrid        ! Length of character string runid
    integer                                    :: m
    integer                                    :: md
    integer                                    :: n
    integer                                    :: nd
    integer                                    :: newlun
    integer                                    :: nm          ! N,M index for discharge location
    integer                                    :: npara       ! Number of parameter records in time dependent direct access file
    integer                                    :: nparrd      ! NR. of parameter records actual read
    integer                                    :: nst_nobub   !< Number of discharges excluding bubble points
    integer                                    :: ntimrd
    real(fp)                                   :: alpha       !< linear interpolation factor
    real(fp)                                   :: rtdis0      !< Previous  read time for discharge time-series
    real(fp)                                   :: rtdis1      !< Following read time for discharge time-series
    logical                                    :: access      ! Flag to read file as direct access or sequential
    logical                                    :: opend       ! Help flag = TRUE when file is still open (Delft3D) and
    character(1)                               :: dumchr      ! Dummy character (#) in first record of direct access file
    character(10)                              :: cntent      ! String with <contents> input
    character(max(300,(4+lstsc)*16))           :: record      ! Record for DIS file (125 or 153)
    character(20)                              :: chlp20
    character(20)                              :: namhlp      ! Name of NAMCON(L) in small characters
    character(20)                              :: namhlp_old  ! Previous value of namhlp
    character(256)                             :: filnam      ! Help var. for file name
    character(36), dimension(:), allocatable   :: parnam      ! Number of parameter records in time dependent direct access files for DIS
    character(36), dimension(3)                :: defpar      ! Default parameters
!
!! executable statements -------------------------------------------------------
!
    nxbub   => gdp%d%nxbub
    timscl  => gdp%gdinidis%timscl
    lundis  => gdp%gdluntmp%lundis
    itdate  => gdp%gdexttim%itdate
    tstop   => gdp%gdexttim%tstop
    dt      => gdp%gdexttim%dt
    scalef  => gdp%gdupddis%scalef
    !
    allocate(parnam(lstsc + 4))
    if (parll) then
       allocate(irecs(nsrcd))
       irecs(:) = -1
    endif
    !
    ddb       = gdp%d%ddbound
    icxy      = max(icx, icy)
    nst_nobub = nsrcd - nxbub
    defpar(1) = 'flux/discharge rate '
    defpar(2) = 'flow magnitude      '
    defpar(3) = 'flow direction      '
    !
    ! define length of RUNID
    !
    call remove_leading_spaces(runid     ,lrid      )
    !
    ! Time dependent discharges
    !
    filnam = 'TMP_' // runid(:lrid) // '.dis'
    !
    ! Test if file is already opened (multi entry Delft3D)
    !
    inquire (file = filnam(:8 + lrid), opened = opend)
    if (.not.opend) then
       lundis = newlun(gdp)
       open (lundis, file = filnam(:8 + lrid), form = 'formatted',              &
            & status = 'old')
       read (lundis, '(a1,i5)', iostat = iocond) dumchr, lrec
       !
       ! erroror EOF (IOCOND <> 0), not allowed
       !
       if (iocond/=0) then
          call prterr(lundia    ,'G007'    ,filnam(:8 + lrid)    )
          error= .true.
          goto 9999
       endif
       close (lundis)
       !
       ! for parallel runs, for discharges in this partion: find locations of discharge in input file
       !
       if (parll) then
          irecrd = 0
          !
          ! file not open as direct access!
          !
          open (lundis, file = filnam(:8 + lrid), form = 'formatted')
   10     continue
          irecrd = irecrd + 1
          read (lundis, '(a)', end=20) record(:lrec - 1)
          iend   = len(record)
          istart = 1
          call srckey(record    ,istart    ,iend      ,ifound    ,gdp       )
          if (ifound == 3) then
             !
             ! location keyword found
             !
             call keyinp(record(istart:iend)  ,chlp20   )
             call small(chlp20    ,len(chlp20)          )
             isrc_nodup = 0
             namhlp_old = ' '
             do isrc = 1, nsrc
                namhlp = namsrc(isrc)
                call small(namhlp, len(namhlp))
                if (namhlp /= namhlp_old) then
                   isrc_nodup = isrc_nodup + 1
                endif
                namhlp_old = namhlp
                if (chlp20 == namhlp) then
                   !
                   ! Discharge record found in the input file
                   ! Store the location of this record table-name in irecs
                   !
                   irecs(isrc_nodup) = irecrd - 2
                   exit
                endif
             enddo
          endif
          goto 10
   20     close (lundis)
       endif
       !
       ! Open file as direct access
       !
       open (lundis, file = filnam(:8 + lrid), form = 'formatted',              &
            & access = 'direct', recl = lrec)
       !
       ! Initialize ITDIS array
       !
       irecrd = 2
       access = .true.
       !
       do isrc = 1, nsrcd
          itdis(1, isrc) = -1
          itdis(2, isrc) = -1
          itdis(3, isrc) = -1
          itdis(4, isrc) = -1
          itdis(5, isrc) = -1
       enddo
    endif
    !
    ! Loop over nsrcd discharge points
    !
    do isrc = 1, nsrcd
       !
       ! No bubble point: mnksrc(6,i) can be used to check whether this point is inside this partition
       ! Bubble point   : mnksrc mismatches with isrc. This is repaired in cnvbub
       !
       if (isrc<=nst_nobub .and. mnksrc(6,isrc)==-1) cycle
       if (parll) then
          irecrd = irecs(isrc)
       endif
       !
       ! Always read umdis/vmdis values from file
       !
       itdis(1, isrc) = -1
       itdis(2, isrc) = -1
       if (itdis(3, isrc)== -1) then
          !
          ! Read initial values from file
          ! IRECRD = 2 as input and start of time varying data as
          ! output
          !
          read (lundis, '(a)', rec = irecrd) record(:lrec - 1)
          npara = 2 + lstsc + 2
          call flhnew(lundis    ,lundia    ,error     ,record(:lrec - 1)    ,access    , &
                    & irecrd    ,namsrc(isrc)         ,cntent    ,disint(isrc)         ,itdate    , &
                    & timscl    ,ntimrd    ,parnam    ,npara     ,nparrd    , &
                    & bubble    ,gdp       )
          if (error) then
             exit
          endif
          !
          ! Test string CNTENT and define momentum
          !
          if (cntent(:10)=='regular   ' .or. cntent(:10)=='walking   ' .or.     &
            & cntent(:10)=='inoutlet  ' .or. cntent(:10)=='power     ' .or.     &
            & cntent(:10)=='culvert   ') then
             if (nparrd/=2 + lstsc) then
                call prterr(lundia    ,'V097'    ,' '       )
                error= .true.
                exit
             endif
             dismmt(isrc) = 'N'
          elseif (cntent(:8)=='momentum') then
             if (nparrd/=4 + lstsc) then
                call prterr(lundia    ,'V097'    ,' '       )
                error= .true.
                exit
             endif
             dismmt(isrc) = 'Y'
          else
             call prterr(lundia    ,'V096'    ,cntent    )
             error= .true.
             exit
          endif
          !
          ! Test name of constituent conform read parameter name
          ! Only first 20 characters are of significance
          !
          if (parnam(2)(:20)/=defpar(1)(:20)) then
             call prterr(lundia    ,'V096'    ,parnam(2)(:20)       )
             error= .true.
             exit
          endif
          !
          do l = 1, lstsc
             namhlp = namcon(l)
             call small(namhlp    ,20        )
             if (parnam(2 + l)(:20)/=namhlp) then
                call prterr(lundia    ,'V096'    ,parnam(2 + l)(:20)   )
                error= .true.
                goto 9999
             endif
          enddo
          !
          if (dismmt(isrc)=='Y') then
             do l = 1, 2
                if (parnam(2 + lstsc + l)(:20)/=defpar(1 + l)(:20)) then
                   call prterr(lundia    ,'V096'    ,parnam(2 + lstsc + l)(:20)      )
                   error= .true.
                   goto 9999
                endif
             enddo
          endif
          !
          ! Define ITDIS (3>5,ISRC) record values
          !
          itdis(3, isrc) = irecrd
          itdis(4, isrc) = irecrd + ntimrd - 1
          itdis(5, isrc) = itdis(3, isrc) - 1
          !
          ! Define start record for next information records
          !
          irecrd = irecrd + ntimrd
       endif
       itdis(5, isrc) = itdis(3, isrc) - 1
       !
       select case (mnksrc(7,isrc))
       case (3,4,5,7)
          !
          ! culvert: initialize parameters on zero and cycle
          !
          disch(isrc) = 0.0_fp
          umdis(isrc) = 0.0_fp
          vmdis(isrc) = 0.0_fp
          do l = 1, lstsc
              rint(l, isrc) = 0.0_fp
          enddo
          cycle
       case default
          !
          ! nothing
          !
       end select
       !
       ! Define discharge location
       ! Test over KCS is not necessary. Discharge location always
       ! in point with KCS = 1 (see CHKDIS)
       !
       ! WARNING: nm is wrong for bubble points (due to extending mnksrc)
       ! This does not harm because momentum is not used for bubble points
       !
       nm = (mnksrc(5, isrc) + ddb) + ((mnksrc(4, isrc) - 1) + ddb)*icxy
       !
       ! Read first time step(s) and discharge data
       !
       call upddis(lundis    ,lundia    ,sferic    ,itdis     , &
                 & isrc      ,nm        ,grdang    ,timnow    ,dt        , &
                 & itfinish  ,timscl    ,nsrcd     ,lstsc     ,j         , &
                 & nmmaxj    ,dismmt    ,alfas     , &
                 & disch0    ,disch1    ,rint0     ,rint1     , &
                 & umdis0    ,umdis1    ,vmdis0    ,vmdis1    ,gdp       )
       if (disint(isrc)=='Y') then
            if (itdis(1,isrc) == itdis(2,isrc)) then
               alpha = 0.0_fp
            else
               rtdis0 = real(itdis(1,isrc),fp)
               rtdis1 = real(itdis(2,isrc),fp)
               alpha  = (timnow-rtdis0) / (rtdis1-rtdis0)
            endif
            disch(isrc) = (1.0_fp-alpha)*disch0(isrc) + alpha*disch1(isrc)
            umdis(isrc) = (1.0_fp-alpha)*umdis0(isrc) + alpha*umdis1(isrc)
            vmdis(isrc) = (1.0_fp-alpha)*vmdis0(isrc) + alpha*vmdis1(isrc)
            do l = 1, lstsc
                rint(l, isrc) = (1.0_fp-alpha)*rint0(l, isrc) + alpha*rint1(l, isrc)
            enddo
       else
            disch(isrc) = disch0(isrc)
            umdis(isrc) = umdis0(isrc)
            vmdis(isrc) = vmdis0(isrc)
            do l = 1, lstsc
                rint(l, isrc) = rint0(l, isrc)
            enddo
       endif
       !
       ! Test and fill KSPU/V(nm,0) array
       ! This does not work for bubble points, because n and m are wrong
       ! (due to extending mnksrc)
       ! This does not harm because a structure is not allowed at bubble points
       !
       if (isrc <= nst_nobub) then
          if (upwsrc == 0) then
               if (dismmt(isrc) == 'Y') then
                   ja_upw = 1
               else
                   ja_upw = -1
               endif
          else
               ja_upw = upwsrc
          endif
          !
          ! Check whether a structure is positioned at intake/outlet
          ! Floating structure is allowed
          !
          ! Intake
          !
          m  = mnksrc(1, isrc)
          n  = mnksrc(2, isrc)
          md = max(1, m - 1)
          nd = max(1, n - 1)
          if ( abs(kspu(n, m, 0))>2 .or. abs(kspu(n , md, 0))>2 .or.  &
               abs(kspv(n, m, 0))>2 .or. abs(kspv(nd, m , 0))>2 ) then
             call prterr(lundia, 'V057', namsrc(isrc))
          endif
          kspu(n , m , 0) = ja_upw
          kspu(n , md, 0) = ja_upw
          kspv(n , m , 0) = ja_upw
          kspv(nd, m , 0) = ja_upw
          !
          ! Outfall
          !
          m  = mnksrc(4, isrc)
          n  = mnksrc(5, isrc)
          md = max(1, m - 1)
          nd = max(1, n - 1)
          if ( abs(kspu(n, m, 0))>2 .or. abs(kspu(n , md, 0))>2 .or.  &
               abs(kspv(n, m, 0))>2 .or. abs(kspv(nd, m , 0))>2 ) then
             call prterr(lundia, 'V057', namsrc(isrc))
          endif
          kspu(n , m , 0) = ja_upw
          kspu(n , md, 0) = ja_upw
          kspv(n , m , 0) = ja_upw
          kspv(nd, m , 0) = ja_upw
       endif
    enddo
    !
 9999 continue
    deallocate(parnam)
    if (parll) deallocate(irecs)
end subroutine inidis

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
module foo_module
  type foo
  contains
    procedure, nopass :: bar
  end type
  interface
    module function bar() result(x)
      real x
    end function
  end interface
end module

module foobar_module
  interface
    subroutine foobar()
      use foo_module, only : foo
    end subroutine
  end interface
end module

max_stars_repo_name: rouson/bug-reports

max_issues_repo_name: rouson/bug-reports

max_forks_repo_name: rouson/bug-reports


! 10. Write a program that prints the next 20 leap years.

    program Exercises
    implicit none
    integer :: currentYear
    integer :: count

    integer dateTimeResult(8)
    character(len=12) unused(3)

    ! Intrinsic subroutine
    call date_and_time(unused(1), unused(2), unused(3), dateTimeResult)

    currentYear = dateTimeResult(1)
    count = 0

10&

    if (mod(currentYear, 4) .ne. 0) then
        ! Not a leap year

    else if (mod(currentYear, 100) .ne. 0) then
        ! It is a leap year
        print *, currentYear
        count = count + 1

    else if (mod(currentYear, 400) .ne. 0) then
        ! Not a leap year

    else
        ! It is a leap year
        print *, currentYear
        count = count + 1
    end if

    if (count .lt. 20) then
        currentYear = currentYear + 1
        goto 10
    end if

    end Program Exercises


text:
GIS Day is a global event originally started by National Geographic and ESRI (The leading vendor in GIS software) and associates in 1998. Its held the third Wednesday of each November, during Geography Awareness Week, and is focused on educating people around the world about GIS technologies and their daily applications.  It is a part of the larger Geography Awareness Week that was started by the National Geographic Society in 1987.

UC Davis has a long history of hosting various events on Geographic Information Systems GIS day that include guest speakers, poster exhibits, map collection tours and demonstrations of technology.

Nov 18, 2009

Events
All Outdoor events will move to Hunt 158 in the event of Rain.

Poster/Map Exhibit
10:004:00  Hunt 158
  Check out maps and posters created by geographers across the campus! Posters and Maps created by various researchers from across campus will be on display showing a wide array of ways to use GIS and related tools in research and careers. Geography Club members will be around to discuss and answer questions.

Demo Computers
10:0012:00  Hunt 158
12:002:00  Tables 2122 by the MU
  Ever wanted to try making map on a computer? Stop by and have a look at some of the options out there including free GIS software.  Geography Club members will be there to show you how it works and answer questions about how you could use GIS.
  Some of the Software that will be available:
   Quantum GIS (Free)
   ArcGIS
   GRASS GIS (Free)
   GPS Drive (Free)

Tabling on the Quad
12:002:00  Tables 2122 by the MU
  Global Positioning System (GPS) units will be available for hands on fun along with people who can teach you all about them.
  Demo Computers for hands on experience, see above for more details
  Geography Quiz

Map Collection Tour
35pm Basement of Shields Library
  An Open House of the Shields Library/Map Collection Map Collection
  Tours of the map collection will be available beginning at 3 PM. Find out more about the 400,000 sheet collection of maps and aerial photographs available for your use.  See what GIS resources the library has to aid you in your research. Meet Kathy Stroud, the Map and GIS librarian.
If you can’t stop by on Wednesday, please drop by at another time. The Map Collection is open 1 PM – 5 PM Monday through Friday and morning access for research purposes may be arranged by emailing mapquestions@lib.ucdavis.edu or calling 5307526196.

Brought to you by
  GLOBAL GLOBAL Geography Club
  Information Center for the Environment http://ice.ucdavis.edu
  UC Davis Shields Library Shields Library/Map Collection Map Collection
  Informatics Center http://ic.ucdavis.edu
  Geography Graduate Group

(If this event had been posted on the Events Board people may have actually known about it... (The California Aggie Daily Calendar would be another good place)...

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
program fme

  use fsystem
  use genoutput
  use paramlist

  use fme2d

  implicit none

  type(t_parlist) :: rparam
  integer :: nDim

  ! Initialise system-wide settings:
  call sys_init()

  ! Initialise the output system.
  call output_init()

  ! Initialise the FEAT 2.0 storage management:
  call storage_init(999, 100)

  ! Read in parameter list
  call parlst_init(rparam)
  call parlst_readfromfile(rparam, './data/master.dat')

  ! fetch parameters
  call parlst_getvalue_int(rparam, '', 'NDIM', nDim, 0)

  ! call exporter
  select case(nDim)
  case (2)
    call fme_2d(rparam)
  end select

  ! Release parameter list
  call parlst_done (rparam)

  ! Print out heap statistics
  call output_lbrk ()
  call storage_info(.true.)

  ! Clean up the storage management, finish
  call storage_done()

end program

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


module szcmb
  use precision1, only : dp
  use constants
  use utils ,only : vec2ang,ang2vec,rotate
  use pix_tools, only : pix2ang_ring,nside2npix
  use alm_tools, only : map2alm,alm2cl
  implicit none

contains

    subroutine cmbcal(nside,dtt,iwrite)
!
!   Calculate CMB sky at observer
!   and return (T0-<T0>)/<T0> as a healpix map
!
        use cosmo_params, only : r_obs,theta_obs,model
        use sznullgeo, only : rayshoot_dr
        implicit none
        integer, intent(in) :: nside
        integer, intent(in) :: iwrite
        real(dp), dimension(0:12*nside*nside-1), intent(out) :: dtt
        real(dp), dimension(0:12*nside*nside-1) :: zp1_inv !1/(z+1)
        real(dp) :: tht_ls,phi_ls,z,zp1_inv_ave,obs_pos(3)
        integer :: npix,ipix
        character(len=28) :: filename

        obs_pos=(/r_obs,theta_obs,0.543*PI/)
! sky appears the same for any azimuthal angle (phi arbitrary)

        npix=nside2npix(nside)
!$omp parallel do private(z,tht_ls,phi_ls) &
!$omp shared(obs_pos,zp1_inv) &
!$omp schedule(guided)
        do ipix=0,npix-1
            call pix2ang_ring(nside,ipix,theta=tht_ls,phi=phi_ls)
            call rayshoot_dr(obs_pos,tht_ls,phi_ls,z,iexit=0,val_exit=-400._dp)
            zp1_inv(ipix)=1._dp/(1._dp+z)
        end do
!$omp end parallel do

        zp1_inv_ave=sum(zp1_inv)/real(npix,kind=dp)
        dtt=zp1_inv/zp1_inv_ave-1._dp

        if (iwrite == 1) then
            write(filename,'(a,a,a,i4.4,a)') &
                    'map_',model,'_Nside',nside,'.txt'
            open(unit=44, file=filename, status='replace')
            do ipix=1,npix
                write(44,*) zp1_inv(ipix), dtt(ipix)
            end do
            close(44)
        end if
        return
    end subroutine cmbcal

    subroutine cmbstats(nside,map)
!
!   Compute dipole and quadrupole
!
        use cosmo_params, only : T0CMB
        implicit none
        real(dp), dimension(:), intent(in) :: map
        integer, intent(in) :: nside
        integer, parameter :: lmax=2,mmax=2
        real(dp), dimension(1:2*nside,1) :: dw8
        complex(dp), dimension(1,0:lmax,0:mmax) :: alm
        real(dp), dimension(0:lmax,1) :: cl
        real(dp) :: dipole,quadrupole,map_max,map_min
        real(dp) :: z,zbounds(2)
        dw8=1._dp
        z=1._dp
        zbounds=(/-z,z/) !(/-1,1/) means no mask
        call map2alm(nside,lmax,mmax,map,alm,zbounds,dw8)
        call alm2cl(lmax,mmax,alm,cl)

        map_max=maxval(map)
        map_min=minval(map)

        write(*,*) 'C_0 :', cl(0,1)
        write(*,*) 'C_1 :', cl(1,1)
        write(*,*) 'C_2 :', cl(2,1)
        write(*,*) 'DT/T dpl :', 1.5_dp*sqrt(cl(1,1)/PI)
        write(*,*) 'map maxval :', map_max
        write(*,*) 'map minval :', map_min
        dipole=1.5_dp*sqrt(cl(1,1)/PI)*T0CMB
        quadrupole=(3._dp/PI)*cl(2,1)*T0CMB**2

        write(*,*) '---------------------------------------------'
        write(*,*) 'Dipole [mK] =', dipole * 1e3_dp
        write(*,*) 'Quad [uK^2] =', quadrupole * 1e12_dp
        write(*,*) '---------------------------------------------'

    end subroutine cmbstats

    function chi_squared(nside)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! This routine computes the chi^2 of the composite sample for the    !
! input Szekeres model:                                              !
!                                                                    !
!     chi^2 = \sum_i [(di - dL) / \sigma di]^2                       !
!                                                                    !
! where di is the observed luminosity distance, dL is the            !
! model luminosity distance and depends on the model parameters, and !
! \sigma di is the error.                                            !
! The luminosity distance as a power series is given by              !
!                                                                    !
!     dL = c/H_0  [z + 1/2 (1-q0)z^2 + ...]                          !
!                                                                    !
! where q0 is the deceleration parameter see Weinberg 2008 (1.4.9),  !
! Visser 2004 [gr-qc/0309109].                                       !
! Then the error (to leading order) is                               !
!                                                                    !
!    \sigma dL \sim \sigma (c z / H_0) = \sigma_v / H_0              !
!                                                                    !
! where cz is the velocity and \sigma_v = c \sigma z.                !
!                                                                    !
! This routine computes dL for the Szekeres model such that the      !
! computed redshift coincides with the observed redshift.            !
! The variables tht \in [0,pi] and phi \in [0,2pi] are the local sky !
! angular coordinates related to the observer's celestial sphere     !
! NOT spacetime theta,phi.                                           !
!                                                                    !
! Composite_LG.txt has 4534 lines and 6 columns corresponding to     !
!                                                                    !
!   v[km/s], dL[Mpc/h], v_pec[km/s], sigma_v[km/s], l[deg], b[deg]   !
!                                                                    !
! where v=cz                                                         !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        use cosmo_params, only : r_obs,theta_obs,FLRW, &
                                 ldipole,bdipole,alpha,amp,r0
        use sznullgeo, only : rayshoot_dr
        implicit none
        real(dp) :: chi_squared
        integer,  intent(in) :: nside
        integer,  parameter  :: Ntot=4534
        real(dp), dimension(0:12*nside*nside-1) :: dtt
        real(dp), dimension(1:Ntot) :: chi2_arr
        real(dp), dimension(1:3)    :: dpl_vec
        real(dp), dimension(1:3)    :: obs_pos
        real(dp) :: vi,di,vpeci,sig_vi,li,bi,sig_di
        real(dp) :: tht_dpl,phi_dpl,tht_dpl0,phi_dpl0
        real(dp) :: tht,phi,dL,z_exit,a1,a2,b1,b2
        real(dp) :: eps,Tdpl
        integer  :: i
        real(dp), parameter :: H_0=100._dp*FLRW%h !km/s/Mpc
        real(dp), parameter :: Tdpl0=5.64_dp,Tdpl_sig=0.03_dp
        real(dp), parameter :: &
                VEL2RED=1._dp/(SPEED_OF_LIGHT*1.0e-3_dp)
! observed dipole angle
        tht_dpl0=HALFPI-bdipole*DEG2RAD
        phi_dpl0=ldipole*DEG2RAD

        call cmbcal(nside,dtt,iwrite=0)
        call dipole_vec(nside,dtt,dpl_vec,Tdpl)
        eps=abs(Tdpl-Tdpl0)
        chi_squared=(eps/Tdpl_sig)**2
        print *, Tdpl
        call vec2ang(dpl_vec,tht_dpl,phi_dpl)

        obs_pos=(/ r_obs, theta_obs, 0.543*PI /)

        a1=phi_dpl0
        a2=HALFPI-tht_dpl0
        b1=HALFPI-tht_dpl
        b2=phi_dpl

        open(unit=77,file='Composite_LG.txt',status='OLD',action='READ')

!$omp parallel do private(vi,di,vpeci,sig_vi,li,bi,sig_di,tht,phi,dL,z_exit) &
!$omp shared(obs_pos,chi2_arr,a1,a2,b1,b2)
        do i=1,Ntot
            read(77,*) vi,di,vpeci,sig_vi,li,bi
            di=di/FLRW%h !evaluate h to get in units Mpc
            sig_di=sig_vi/H_0
            z_exit=vi*VEL2RED
            call galactic2szekeres(li,bi,tht,phi,a1,a2,b1,b2)
            call rayshoot_dr(obs_pos,tht,phi,dL,iexit=1,val_exit=z_exit)
            chi2_arr(i)=((di-dL)/sig_di)**2
        end do
!$omp end parallel do

        chi_squared=chi_squared+sum(chi2_arr)

        close(77)
        return
    end function chi_squared

    subroutine dipole_vec(nside,map,dpl_vec,Tdipole)
        use cosmo_params, only : T0CMB
        implicit none
        integer, intent(in) :: nside
        real(dp), dimension(:), intent(in) :: map
        real(dp), dimension(3), intent(out) :: dpl_vec
        real(dp), dimension(1:2*nside,1) :: dw8
        real(dp), dimension(2) :: zbounds=(/-1._dp,1._dp/)
        integer, parameter :: lmax=2,mmax=2
        complex(dp), dimension(1,0:lmax,0:mmax) :: alm
        real(dp), dimension(0:lmax,1) :: cl
        real(dp) :: a10,a11r,a11i,norm,Tdipole
        dw8=1._dp
        call map2alm(nside,lmax,mmax,map,alm,zbounds,dw8)
        a10=real(alm(1,1,0))
        a11r=real(alm(1,1,1))
        a11i=aimag(alm(1,1,1))
        dpl_vec(1)=-SQRT2*a11r
        dpl_vec(2)=+SQRT2*a11i
        dpl_vec(3)=+a10
        norm=sqrt(a10*a10+2._dp*(a11r**2+a11i**2))
        dpl_vec=dpl_vec/norm
! compute Temperature dipole (magnitude)
        call alm2cl(lmax,mmax,alm,cl)
        Tdipole=1.5_dp*sqrt(cl(1,1)/PI)*T0CMB*1e3_dp
    end subroutine dipole_vec

    subroutine galactic2szekeres(l,b,theta,phi,a1,a2,b1,b2) !formerly invrotmap
!
!   Maps galactic coordinates (l,b) [deg] to szekeres local
!   sky coordinates (theta,phi) (NOT spacetime theta,phi)
!
        implicit none
        real(dp), intent(in) :: l,b
        real(dp), intent(out) :: theta,phi
        real(dp), intent(in) :: a1,a2,b1,b2
        real(dp), dimension(3) :: vector
        theta=HALFPI-b*DEG2RAD
        phi=l*DEG2RAD
        call ang2vec(theta,phi,vector)
! rotate dipole vector onto ex=(1,0,0)
        call rotate(vector,psi=-a1,axis='z')
        call rotate(vector,psi=-a2,axis='y')
! rotate ex onto the found dipole in szekeres
        call rotate(vector,psi=+b1,axis='y')
        call rotate(vector,psi=+b2,axis='z')
        call vec2ang(vector,theta,phi)
    end subroutine galactic2szekeres

end module szcmb


!------------------------------------------------------------------------------
! BSD 2-Clause License
!
! Copyright (c) 2021, Science and Technology Facilities Council
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! * Redistributions of source code must retain the above copyright notice, this
!   list of conditions and the following disclaimer.
!
! * Redistributions in binary form must reproduce the above copyright notice,
!   this list of conditions and the following disclaimer in the documentation
!   and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!------------------------------------------------------------------------------
! Author: S. Siso, STFC Daresbury Laboratory

! Mock device implementation using separate 1D host memory
module virtual_device
    use field_mod
    use kind_params_mod, only: go_wp
    use iso_c_binding, only: c_loc, c_ptr, c_f_pointer
    implicit none

contains

    subroutine init_device_memory(field)
        type(r2d_field), intent(inout) :: field
        real(go_wp), dimension(:), allocatable, target :: device_memory

        allocate(device_memory(size(field%data,1) * size(field%data,2)))
        field%device_ptr = C_LOC(device_memory)
        field%data_on_device = .true.
        field%read_from_device_f => read_from_device_impl
        field%write_to_device_f => write_to_device_impl
    end subroutine

    subroutine read_from_device_impl(from, to, startx, starty, nx, ny, blocking)
        type(c_ptr), intent(in) :: from
        real(go_wp), dimension(:,:), target, intent(inout) :: to
        integer, intent(in) :: startx, starty, nx, ny
        logical, intent(in) :: blocking
        real(go_wp), dimension(:), pointer :: device_memory
        integer :: next_offset, i

        write(*, *) "Read operation", startx, starty, nx, ny, blocking
        call C_F_POINTER(from, device_memory, [size(to,1) * size(to,2)])

        next_offset = ((starty - 1) * size(to,1) + (startx-1)) + 1

        do i = starty, starty + ny - 1
            ! Copy next contiguous chunk
            to(startx:startx+nx-1, i) = device_memory(next_offset:next_offset+nx-1)
            next_offset = next_offset + size(to,1)
        enddo

    end subroutine read_from_device_impl

    subroutine write_to_device_impl(from, to, startx, starty, nx, ny, blocking)
        real(go_wp), dimension(:,:), target, intent(in) :: from
        type(c_ptr), intent(in) :: to
        integer, intent(in) ::  startx, starty, nx, ny
        logical, intent(in) :: blocking
        real(go_wp), dimension(:), pointer :: device_memory
        integer :: i, next_offset

        write(*, *) "Write operation", startx, starty, nx, ny, blocking
        call C_F_POINTER(to, device_memory, [size(from,1) * size(from,2)])

        next_offset = ((starty - 1) * size(from,1) + (startx-1)) + 1

        do i = starty, starty + ny - 1
            ! Copy next contiguous chunk
            device_memory(next_offset:next_offset+nx-1) = from(startx:startx+nx-1, i)
            next_offset = next_offset + size(from,1)
        enddo

    end subroutine write_to_device_impl

    subroutine simulate_device_computation(buffer, total_size)
        type(c_ptr), intent(in) :: buffer
        integer, intent(in) :: total_size
        real(go_wp), dimension(:), pointer :: device_memory

        write(*, *) "Device computation"
        call C_F_POINTER(buffer, device_memory, [total_size])
        device_memory = device_memory * 2

    end subroutine simulate_device_computation

end module virtual_device

!> Tests for the field read and write to device functionality.
program test_device_io
    use kind_params_mod
    use parallel_mod
    use grid_mod
    use field_mod
    use gocean_mod
    use virtual_device
    implicit none
    ! Total size of the model domain
    integer :: jpiglo = 5, jpjglo = 5
    ! (Uniform) grid resolution
    real(go_wp) :: dx = 1.0
    real(go_wp) :: dy = 1.0
    !> The grid on which our fields are defined
    type(grid_type), target :: model_grid
    !> An example field
    type(r2d_field) :: test_field
    ! Local definition of the T-point mask which defines whether T points are
    ! wet (1), dry (0) or outside (-1) the simulation domain.
    integer, allocatable :: tmask(:,:)
    integer :: ierr
    integer :: i

    call gocean_initialise()

    !> Create our grid object
    model_grid = grid_type(GO_ARAKAWA_C, &
                           (/GO_BC_EXTERNAL,GO_BC_EXTERNAL,GO_BC_NONE/), &
                           GO_OFFSET_NE)

    !> Generate a domain decomposition. This automatically uses the number
    !! of MPI ranks available.
    call model_grid%decompose(jpiglo, jpjglo)

    !> Create a T-point mask describing the (local) domain
    allocate(tmask(model_grid%subdomain%global%nx,  &
             model_grid%subdomain%global%ny), Stat=ierr)
    if(ierr /= 0)then
        call gocean_stop('Failed to allocate T-mask')
    end if
    ! To keep things simple we set all points to be wet and within the domain
    tmask = 1

    !> Complete the initialisation of the grid using the T-mask and
    !! grid resolution
    call grid_init(model_grid, dx, dy, tmask)

    !> Create fields on U,V,T,F-points of the grid
    test_field = r2d_field(model_grid, GO_U_POINTS)
    call init_device_memory(test_field)

    test_field%data = 0
    call test_field%write_to_device()  ! All device data is 0
    test_field%data = 1
    call test_field%write_to_device(2, 2, 5, 5)  ! A 5x5 block starting at 2,2 is 1
    call simulate_device_computation(test_field%device_ptr, 8*8) ! Double device values
    call test_field%read_from_device(5, 5, 4, 4) ! Read the bottom-right quadrant

    write(*,*) "Resulting array:"
    do i=1, 8
        write(*, "(10f5.1)") test_field%data(:, i)
    enddo

    if ( &
        ! Data not read back must be 1.0
        test_field%data(1, 1) /= 1.0 .or. &
        ! Data written, doubled and read is 2.0
        test_field%data(5, 5) /= 2.0 .or. &
        ! Data not written but read back is still 0.0
        test_field%data(7, 7) /= 0.0 ) then
        stop "Error - Results are incorrect"
    endif

    write(*,*) "Test passed successfuly! Note that this test should still ", &
               "succeed when changing the DL_ESM_ALIGNMENT environment variable."
    call gocean_finalise()

end program test_device_io


text:
!     This software is Copyright (c) 2012-2015 The Regents of the
!     University of California. All Rights Reserved.
!
!     Permission to copy and modify this software and its documentation
!     for educational, research and non-profit purposes, without fee,
!     and without a written agreement is hereby granted, provided that
!     the above copyright notice, this paragraph and the following three
!     paragraphs appear in all copies.
!
!     Permission to make commercial use of this software may be obtained
!     by contacting:
!
!     Technology Transfer Office
!     9500 Gilman Drive, Mail Code 0910
!     University of California
!     La Jolla, CA 92093-0910
!     (858) 534-5815
!     invent@ucsd.edu
!
!     This software program and documentation are copyrighted by The
!     Regents of the University of California. The software program and
!     documentation are supplied "as is", without any accompanying
!     services from The Regents. The Regents does not warrant that the
!     operation of the program will be uninterrupted or error-free. The
!     end-user understands that the program was developed for research
!     purposes and is advised not to rely exclusively on the program for
!     any reason.
!
!     IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY
!     PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
!     DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS
!     SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
!     CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!     THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY
!     WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
!     OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE
!     SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE
!     UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
!     MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

      SUBROUTINE CGRAD(nFaces, dof, nNo, nnz, mynNo, commu, cS, face,
     2   ls, rowPtr, colPtr, D, G, L, R)

      INCLUDE "svLS_STD.h"

      INTEGER, INTENT(IN) :: nFaces, dof, nNo, nnz, mynNo
      TYPE(svLS_commuType), INTENT(IN) :: commu
      TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks)
      TYPE(svLS_faceType), INTENT(IN) :: face(nFaces)
      TYPE(svLS_subLsType), INTENT(INOUT) :: ls
      INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz)
      REAL*8, INTENT(IN) :: D(dof,nnz), G(dof,nnz), L(nnz)
      REAL*8, INTENT(INOUT) :: R(nNo)

      INTEGER i
      REAL*8 errO, err, alpha, eps, time
      REAL*8 CPUT, NORMS, DOTS
      REAL*8, ALLOCATABLE :: X(:), P(:), SP(:), DGP(:), GP(:,:),
     2   unCondU(:,:)

      ALLOCATE(X(nNo), P(nNo), SP(nNo), DGP(nNo), GP(dof,nNo),
     2   unCondU(dof,nNo))

      time     = CPUT()
      ls%suc   = .FALSE.
      errO     = NORMS(mynNo, commu, R)
      ls%iNorm = errO
      eps      = MAX(ls%absTol,ls%relTol*errO)
      eps      = eps*eps
      X        = 0D0

      err      = errO
      err      = err*err
      P        = R
      DO i=1, ls%mItr
         IF (err .LT. eps) THEN
            ls%suc = .TRUE.
            EXIT
         END IF
         errO = err
         CALL SPARMULSV(dof, nNo, nnz, commu, cS, rowPtr, colPtr,
     2      G, P, GP)

         IF (ANY(face%coupledFlag)) THEN
            unCondU = GP
            CALL ADDBCMUL(BCOP_TYPE_PRE, nFaces, dof, nNo, mynNo, commu,
     2         face, unCondU, GP)
         END IF

         CALL SPARMULVS(dof, nNo, nnz, commu, cS, rowPtr, colPtr,
     2      D, GP, DGP)

         CALL SPARMULSS(     nNo, nnz, commu, cS, rowPtr, colPtr,
     2      L, P, SP)

         SP    = SP - DGP
         alpha = errO/DOTS(mynNo, commu, P, SP)
         X     = X + alpha*P
         R     = R - alpha*SP
         err   = NORMS(mynNo, commu, R)
         err   = err*err
         P     = R + err/errO*P
      END DO
      R        = X
      ls%fNorm = SQRT(err)
      ls%callD = CPUT() - time + ls%callD
      ls%dB    = 5D0*LOG(err/errO)
      ls%itr   = ls%itr + i - 1

      DEALLOCATE(X, P, SP, DGP, GP)

      RETURN
      END SUBROUTINE CGRAD

!====================================================================

      SUBROUTINE CGRADV(dof, nNo, nnz, mynNo, commu, cS, ls,
     2   rowPtr, colPtr, K, R)

      INCLUDE "svLS_STD.h"

      INTEGER, INTENT(IN) :: dof, nNo, nnz, mynNo
      TYPE(svLS_commuType), INTENT(IN) :: commu
      TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks)
      TYPE(svLS_subLsType), INTENT(INOUT) :: ls
      INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz)
      REAL*8, INTENT(IN) :: K(dof*dof,nnz)
      REAL*8, INTENT(INOUT) :: R(dof,nNo)

      INTEGER i
      REAL*8 errO, err, alpha, eps
      REAL*8 CPUT, NORMV, DOTV
      REAL*8, ALLOCATABLE :: P(:,:), KP(:,:), X(:,:)

      ALLOCATE(P(dof,nNo), KP(dof,nNo), X(dof,nNo))

      ls%callD = CPUT()
      ls%suc   = .FALSE.
      err      = NORMV(dof, mynNo, commu, R)
      ls%iNorm = err
      eps      = MAX(ls%absTol,ls%relTol*err)
      eps      = eps*eps
      err      = err*err
      X        = 0D0
      P        = R

      DO i=1, ls%mItr
         IF (err .LT. eps) THEN
            ls%suc = .TRUE.
            EXIT
         END IF
         errO = err
         CALL SPARMULVV(dof, nNo, nnz, commu, cS, rowPtr, colPtr,
     2      K, P, KP)

         alpha = errO/DOTV(dof, mynNo, commu, P, KP)
         X     = X + alpha*P
         R     = R - alpha*KP
         err   = NORMV(dof, mynNo, commu, R)
         err   = err*err
         P = R + err/errO*P
      END DO

      R        = X
      ls%itr   = i - 1
      ls%fNorm = SQRT(err)
      ls%callD = CPUT() - ls%callD
      ls%dB    = 5D0*LOG(err/errO)
      IF (i .GT. ls%mItr) ls%itr = ls%mItr

      RETURN
      END SUBROUTINE CGRADV

!====================================================================

      SUBROUTINE CGRADS(nNo, nnz, mynNo, commu, cS, ls,
     2   rowPtr, colPtr, K, R)

      INCLUDE "svLS_STD.h"

      INTEGER, INTENT(IN) :: nNo, nnz, mynNo
      TYPE(svLS_commuType), INTENT(IN) :: commu
      TYPE(svLS_cSType), INTENT(IN) :: cS(commu%nTasks)
      TYPE(svLS_subLsType), INTENT(INOUT) :: ls
      INTEGER, INTENT(IN) :: rowPtr(2,nNo), colPtr(nnz)
      REAL*8, INTENT(IN) :: K(nnz)
      REAL*8, INTENT(INOUT) :: R(nNo)

      INTEGER i
      REAL*8 errO, err, alpha, eps
      REAL*8 CPUT, NORMS, DOTS
      REAL*8, ALLOCATABLE :: P(:), KP(:), X(:)

      ALLOCATE(P(nNo), KP(nNo), X(nNo))

      ls%callD = CPUT()
      ls%suc   = .FALSE.
      err      = NORMS(mynNo, commu, R)
      ls%iNorm = err
      eps      = MAX(ls%absTol,ls%relTol*err)
      eps      = eps*eps
      err      = err*err
      X        = 0D0
      P        = R

      DO i=1, ls%mItr
         IF (err .LT. eps) THEN
            ls%suc = .TRUE.
            EXIT
         END IF
         errO = err
         CALL SPARMULSS(nNo, nnz, commu, cS, rowPtr, colPtr, K, P, KP)
         alpha = errO/DOTS(mynNo, commu, P, KP)
         X     = X + alpha*P
         R     = R - alpha*KP
         err   = NORMS(mynNo, commu, R)
         err   = err*err
         P = R + err/errO*P
      END DO

      R        = X
      ls%itr   = i - 1
      ls%fNorm = SQRT(err)
      ls%callD = CPUT() - ls%callD
      ls%dB    = 5D0*LOG(err/errO)
      IF (i .GT. ls%mItr) ls%itr = ls%mItr

      RETURN
      END SUBROUTINE CGRADS


text:
! @@name:	atomic.2f
! @@type:	F-fixed
! @@compilable:	yes
! @@linkable:	no
! @@expect:	success
       function atomic_read(p)
       integer :: atomic_read
       integer, intent(in) :: p
! Guarantee that the entire value of p is read atomically. No part of
! p can change during the read operation.

!$omp atomic read
       atomic_read = p
       return
       end function atomic_read

       subroutine atomic_write(p, value)
       integer, intent(out) :: p
       integer, intent(in) :: value
! Guarantee that value is stored atomically into p. No part of p can change
! until after the entire write operation is completed.
!$omp atomic write
       p = value
       end subroutine atomic_write

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
module thermopack_constants
  ! Constants.
  implicit none
  save
  public

  real, parameter :: Tmax=2000.0 !< K
  real, parameter :: Tmin=50.0 !< K
  real, parameter :: kB_const=1.380649e-23 !< J/K (Boltzmanns const)
  real, parameter :: h_const=6.626069311e-34  !< Js (Plancks constant)
  real, parameter :: N_Avogadro = 6.02214076e23
  real, parameter :: Rgas_default = N_Avogadro*kB_const
  real :: Rgas=Rgas_default !< J/mol/K
  real :: kRgas=1000.0*Rgas_default !< J/kmol/K
  real :: tpTmax=999.0 !< K
  real :: tpTmin=80.0 !< K
  real :: tpPmax=1.0e8 !< Pa
  real :: tpPmin=1.0e1 !< Pa
  !> Apply special treatement to phases with low consentrations of electrolytes
  real, parameter :: elEps = 1.0e-25
  !> Trace components, and zero components
  real, parameter :: traceEps = 1.0e-20
  !< Minimum mol number for apparent composition
  real, parameter :: min_mol_num = 1.0e-150
  !> String length
  integer, parameter :: clen=2048
  !> Control output of debug information
  logical :: verbose = .false.
  !> Phase identifiers
  integer, parameter :: TWOPH=0,LIQPH=1,VAPPH=2,MINGIBBSPH=3,&
       SINGLEPH=4,SOLIDPH=5,FAKEPH=6,VAPSOLPH=7
  !> Liquid phase type identifiers
  integer, parameter :: NONWATER=-1, WATER=-2
  !> Library used to solve EoS
  integer, parameter :: THERMOPACK=1, TREND=2
  !> Method to discriminate between liquid and vapor in the event of an
  !> undefined single phase
  integer, parameter :: PSEUDO_CRIT_ZFAC=1, PSEUDO_CRIT_MOLAR_VOLUME=2, &
       VOLUME_COVOLUME_RATIO=3
  !> Ignore components where z <= zLimit
  real, parameter :: zLimit = 0.0
  !> Continue on error?
  logical :: continueOnError = .false.
  !> Test type
  integer, parameter :: GRID = 1, ENVELOPE_PL = 2, &
       BINARY_PL = 3, BINARY_VLLE_PL = 4, META_LIMIT_PL = 5, SOLIDENVELOPE_PL = 6
  !> Property type
  integer, parameter :: PROP_OVERALL=0, PROP_RESIDUAL=1, PROP_IDEAL=2
  !> String lengths
  integer, parameter :: uid_len = 20
  integer, parameter :: ref_len = 40
  integer, parameter :: comp_name_len = 40
  integer, parameter :: formula_len = 20
  integer, parameter :: bibref_len = 80
  integer, parameter :: eosid_len = 20
  integer, parameter :: mix_len = 20
  integer, parameter :: eos_name_len = 30
  integer, parameter :: short_label_len = 20
  integer, parameter :: label_len = 100

contains

  subroutine set_Rgas(RgasIn)
    implicit none
    real, intent(in) :: RgasIn
    Rgas = RgasIn
    kRgas = 1000.0*Rgas
  end subroutine set_Rgas

  !----------------------------------------------------------------------
  subroutine get_templimits(Tmin, Tmax)
    !> Get EoSlib-specific max/min supported temperature
    !>
    !> \author EA, 2014-05
    implicit none
    ! Output:
    real,     intent(out) :: Tmin !< Minimum supported temperature (K)
    real,     intent(out) :: Tmax !< Maximum supported temperature (K)

    Tmin = tpTmin
    Tmax = tpTmax
  end subroutine get_templimits

  subroutine get_presslimits(Pmin, Pmax)
    !> Get EoSlib-specific max/min supported pressure
    implicit none
    ! Output:
    real,     intent(out) :: Pmin !< Minimum supported temperature (Pa)
    real,     intent(out) :: Pmax !< Maximum supported temperature (Pa)

    Pmin = tpPmin
    Pmax = tpPmax
  end subroutine get_presslimits

  !----------------------------------------------------------------------
  subroutine phaseIntToName(phase,phaseName)
    integer, intent(in) :: phase
    character(len=*), intent(out) :: phaseName
    !
    if (len(phaseName) >= 10) then
      select case(phase)
      case (TWOPH)
        phaseName = 'TWOPH'
      case (LIQPH)
        phaseName = 'LIQPH'
      case (VAPPH)
        phaseName = 'VAPPH'
      case (MINGIBBSPH)
        phaseName = 'MINGIBBSPH'
      case (SINGLEPH)
        phaseName = 'SINGLEPH'
      case (SOLIDPH)
        phaseName = 'SOLIDPH'
      case (FAKEPH)
        phaseName = 'FAKEPH'
      end select
    endif
  end subroutine phaseIntToName

end module thermopack_constants

max_stars_repo_name: SINTEF/Thermopack

max_issues_repo_name: SINTEF/Thermopack

max_forks_repo_name: SINTEF/Thermopack


module precision
  !  This module will assign fdp (floating point double precision)
  !  isp (integer single precision), idp (integer double precision)
  !  csp (complex single precision), cdp (complex double precision)
  !  fqp (floating point quadratic precision)
  implicit none
  integer, parameter :: fqp=selected_real_kind(20,200)
  integer, parameter :: fdp=selected_real_kind(10,200)
  integer, parameter :: fsp=selected_real_kind(5,50)
  integer, parameter :: isp=selected_int_kind(6)
  integer, parameter :: idp=selected_int_kind(9)
  integer, parameter :: cdp=fdp
  integer, parameter :: csp=fsp
end module precision


PROGRAM LA_CGESDD_ET_EXAMPLE
!
!  -- LAPACK95 interface driver routine (version 3.0) --
!     UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK
!     September, 2000
!
!  .. USE STATEMENTS
   USE LA_PRECISION, ONLY: WP => SP
   USE F95_LAPACK, ONLY: LA_GESDD
!  .. IMPLICIT STATEMENT ..
   IMPLICIT NONE
!  .. PARAMETERS ..
   CHARACTER(LEN=*), PARAMETER :: FMT = '(4(1X,1H(,F7.3,1H,,F7.3,1H):))'
   INTEGER, PARAMETER :: NIN=5, NOUT=6
!  .. LOCAL SCALARS ..
   INTEGER :: I, INFO, M, N, MN
!  .. LOCAL ARRAYS ..
   CHARACTER(LEN=1) :: JOB
   REAL(WP), ALLOCATABLE :: AA(:,:), S(:), WW(:)
   COMPLEX(WP), ALLOCATABLE :: A(:,:), U(:,:), VT(:,:), DUMMY(:,:)
!  .. EXECUTABLE STATEMENTS ..
   WRITE (NOUT,*) 'CGESDD ET_Example Program Results.'
   READ ( NIN, * )   ! SKIP HEADING IN DATA FILE
   READ ( NIN, * ) M, N
   PRINT *, 'M = ', M, ' N = ', N
   MN = MIN(M,N)
   ALLOCATE ( A(M,N), AA(M,N), S(MN), U(M,M), VT(N,N), WW(1:MN-1) )
!
      READ (NIN, *) AA
      A=AA
      WRITE(NOUT,*) 'The matrix A:'
      DO I = 1, M; WRITE (NOUT,*) 'I = ', I; WRITE (NOUT,FMT) A(I,:); ENDDO
!
   WRITE ( NOUT, * )'---------------------------------------------------------'
   WRITE ( NOUT, * )
   WRITE ( NOUT, * )'Details of LA_CGESDD LAPACK Subroutine Results.'
   WRITE ( NOUT, * )
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U, VT, WW, ''N'', INFO )'
   A=AA
   CALL LA_GESDD( A, S, U, VT, WW, 'N', INFO )
   WRITE(NOUT,*) 'INFO = ', INFO, ' Singular values:'
   WRITE(NOUT,FMT) S
   WRITE(NOUT,*) 'Orthogonal/unitary matrix U):'
   DO I = 1, M; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) U(:,I); END DO
   WRITE(NOUT,*) 'Orthogonal/unitary matrix V**H):'
   DO I = 1, N; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) VT(:,I); END DO
   IF( INFO > 0 )THEN
      WRITE(NOUT,*) 'Unconverged superdiagonal elements:'
      WRITE(NOUT,FMT) WW
   ENDIF
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U, VT, WW, ''U'', INFO )'
   A=AA
   CALL LA_GESDD( A, S, U, VT, WW, 'U', INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U, VT, WW, ''V'', INFO )'
   A=AA
   CALL LA_GESDD( A, S, U, VT, WW, 'V', INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U, VT, WW, ''B'', INFO )'
   A=AA
   CALL LA_GESDD( A, S, U, VT, WW, 'B', INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, VT=VT, WW=WW, JOB=''U'', INFO=INFO )'
   A=AA
   CALL LA_GESDD( A, S, VT=VT, WW=WW, JOB='U', INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U=U, WW=WW, JOB=''V'', INFO=INFO )'
   A=AA
   CALL LA_GESDD( A, S, U=U, WW=WW, JOB='V', INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, WW=WW, JOB=''V'', INFO=INFO )'
   A=AA
   CALL LA_GESDD( A, S, WW=WW, JOB='V', INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, WW=WW, JOB=''U'', INFO=INFO )'
   A=AA
   CALL LA_GESDD( A, S, WW=WW, JOB='U', INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, WW=WW, JOB=''N'', INFO=INFO )'
   A=AA
   CALL LA_GESDD( A, S, WW=WW, JOB='N', INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U, VT )'
   A=AA
   CALL LA_GESDD( A, S, U, VT )
   WRITE(NOUT,*) ' Singular values:'
   WRITE(NOUT,FMT) S
   WRITE(NOUT,*) 'Orthogonal/unitary matrix U):'
   DO I = 1, M; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) U(:,I); END DO
   WRITE(NOUT,*) 'Orthogonal/unitary matrix V**H):'
   DO I = 1, N; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) VT(:,I); END DO
!
     WRITE(NOUT,*)
     WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U(:,1:MN), VT(1:MN,:) )'
     A=AA; U = HUGE(1.0_WP); VT = HUGE(1.0_WP)
     CALL LA_GESDD( A, S, U(:,1:MN), VT(1:MN,:) )
     WRITE(NOUT,*) ' Singular values:'
     WRITE(NOUT,FMT) S
     WRITE(NOUT,*) 'Orthogonal/unitary matrix U):'
     DO I = 1, M; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) U(:,I); END DO
       WRITE(NOUT,*) 'Orthogonal/unitary matrix V**H):'
       DO I = 1, N; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) VT(:,I); END DO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S )'
   A=AA
   CALL LA_GESDD( A, S )
   WRITE(NOUT,*) ' Singular values:'
   WRITE(NOUT,FMT) S
! STARTING ERROR TESTS
! ERROR 1
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( DUMMY, S, INFO=INFO )'
   A=AA
   CALL LA_GESDD( DUMMY, S, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 2
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S(1:MAX(M,N)+1), INFO=INFO )'
   A=AA
!   CALL LA_GESDD( A, S(1:MAX(M,N)+1), INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 3
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U(1:M-1,:), VT, WW, ''N'',INFO )'
   A=AA
!   CALL LA_GESDD( A, S, U(1:M-1,:), VT, WW, 'N', INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 3
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U(:,1:MN-1), VT, WW, ''N'',INFO )'
   A=AA
   CALL LA_GESDD( A, S, U(:,1:MN-1), VT, WW, 'N', INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 4
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U, VT(1:MN-1,:), WW, ''N'',INFO )'
   A=AA
   CALL LA_GESDD( A, S, U, VT(1:MN-1,:), WW, 'N', INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 4
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, U, VT(:,1:N-1), WW, ''N'',INFO )'
   A=AA
   CALL LA_GESDD( A, S, U, VT(:,1:N-1), WW, 'N', INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 6
   JOB = 'T'
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_GESDD( A, S, JOB=JOB, INFO=INFO )'
   A=AA
   CALL LA_GESDD( A, S,JOB=JOB, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
END PROGRAM LA_CGESDD_ET_EXAMPLE


text:
program gemv_test
    use :: util_mod
    implicit none
    integer, parameter :: rows = 3, cols = 2
    integer :: i
    real, dimension(rows, cols) :: A = reshape( [ (real(i), i=1, rows*cols) ], [ rows, cols ])
    real, dimension(cols) :: v = [ 1.0, 2.0 ]
    real, dimension(rows) :: u
    character(len=12) :: fmt_str = '(*(F6.1))'

    call print_matrix(A, 'A', fmt_str)
    call print_vector(v, 'v', fmt_str)
    call sgemv('n', rows, cols, 1.0, A, rows, v, 1, 0.0, u, 1)
    call print_vector(u, 'u', fmt_str)

end program gemv_test

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: gjbex/CMake-intro

max_forks_repo_name: gjbex/CMake-intro


text:
      program unstr

c     test of entry and exit node extraction by Ronan's restructurer
c     plus unstructured decomposition

 200  continue
      j = 3
      print *, j
      if(j.lt.3) go to 200

      print *, j

      end


text:
*DECK SVOUT
      SUBROUTINE SVOUT (N, SX, IFMT, IDIGIT)
C***BEGIN PROLOGUE  SVOUT
C***SUBSIDIARY
C***PURPOSE  Subsidiary to SPLP
C***LIBRARY   SLATEC
C***TYPE      SINGLE PRECISION (SVOUT-S, DVOUT-D)
C***AUTHOR  (UNKNOWN)
C***DESCRIPTION
C
C     SINGLE PRECISION VECTOR OUTPUT ROUTINE.
C
C  INPUT..
C
C  N,SX(*) PRINT THE SINGLE PRECISION ARRAY SX(I),I=1,...,N, ON
C          OUTPUT UNIT LOUT. THE HEADING IN THE FORTRAN FORMAT
C          STATEMENT IFMT(*), DESCRIBED BELOW, IS PRINTED AS A FIRST
C          STEP. THE COMPONENTS SX(I) ARE INDEXED, ON OUTPUT,
C          IN A PLEASANT FORMAT.
C  IFMT(*) A FORTRAN FORMAT STATEMENT. THIS IS PRINTED ON OUTPUT
C          UNIT LOUT WITH THE VARIABLE FORMAT FORTRAN STATEMENT
C                WRITE(LOUT,IFMT)
C  IDIGIT  PRINT AT LEAST ABS(IDIGIT) DECIMAL DIGITS PER NUMBER.
C          THE SUBPROGRAM WILL CHOOSE THAT INTEGER 4,6,10 OR 14
C          WHICH WILL PRINT AT LEAST ABS(IDIGIT) NUMBER OF
C          PLACES.  IF IDIGIT.LT.0, 72 PRINTING COLUMNS ARE UTILIZED
C          TO WRITE EACH LINE OF OUTPUT OF THE ARRAY SX(*). (THIS
C          CAN BE USED ON MOST TIME-SHARING TERMINALS). IF
C          IDIGIT.GE.0, 133 PRINTING COLUMNS ARE UTILIZED. (THIS CAN
C          BE USED ON MOST LINE PRINTERS).
C
C  EXAMPLE..
C
C  PRINT AN ARRAY CALLED (COSTS OF PURCHASES) OF LENGTH 100 SHOWING
C  6 DECIMAL DIGITS PER NUMBER. THE USER IS RUNNING ON A TIME-SHARING
C  SYSTEM WITH A 72 COLUMN OUTPUT DEVICE.
C
C     DIMENSION COSTS(100)
C     N = 100
C     IDIGIT = -6
C     CALL SVOUT(N,COSTS,'(''1COSTS OF PURCHASES'')',IDIGIT)
C
C***SEE ALSO  SPLP
C***ROUTINES CALLED  I1MACH
C***REVISION HISTORY  (YYMMDD)
C   811215  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   891107  Added comma after 1P edit descriptor in FORMAT
C           statements.  (WRB)
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900328  Added TYPE section.  (WRB)
C***END PROLOGUE  SVOUT
      DIMENSION SX(*)
      CHARACTER IFMT*(*)
C
C     GET THE UNIT NUMBER WHERE OUTPUT WILL BE WRITTEN.
C***FIRST EXECUTABLE STATEMENT  SVOUT
      J=2
      LOUT=I1MACH(J)
      WRITE(LOUT,IFMT)
      IF(N.LE.0) RETURN
      NDIGIT = IDIGIT
      IF(IDIGIT.EQ.0) NDIGIT = 4
      IF(IDIGIT.GE.0) GO TO 80
C
      NDIGIT = -IDIGIT
      IF(NDIGIT.GT.4) GO TO 20
C
      DO 10 K1=1,N,5
      K2 = MIN(N,K1+4)
      WRITE(LOUT,1000) K1,K2,(SX(I),I=K1,K2)
   10 CONTINUE
      RETURN
C
   20 CONTINUE
      IF(NDIGIT.GT.6) GO TO 40
C
      DO 30 K1=1,N,4
      K2 = MIN(N,K1+3)
      WRITE(LOUT,1001) K1,K2,(SX(I),I=K1,K2)
   30 CONTINUE
      RETURN
C
   40 CONTINUE
      IF(NDIGIT.GT.10) GO TO 60
C
      DO 50 K1=1,N,3
      K2=MIN(N,K1+2)
      WRITE(LOUT,1002) K1,K2,(SX(I),I=K1,K2)
   50 CONTINUE
      RETURN
C
   60 CONTINUE
      DO 70 K1=1,N,2
      K2 = MIN(N,K1+1)
      WRITE(LOUT,1003) K1,K2,(SX(I),I=K1,K2)
   70 CONTINUE
      RETURN
C
   80 CONTINUE
      IF(NDIGIT.GT.4) GO TO 100
C
      DO 90 K1=1,N,10
      K2 = MIN(N,K1+9)
      WRITE(LOUT,1000) K1,K2,(SX(I),I=K1,K2)
   90 CONTINUE
      RETURN
C
  100 CONTINUE
      IF(NDIGIT.GT.6) GO TO 120
C
      DO 110 K1=1,N,8
      K2 = MIN(N,K1+7)
      WRITE(LOUT,1001) K1,K2,(SX(I),I=K1,K2)
  110 CONTINUE
      RETURN
C
  120 CONTINUE
      IF(NDIGIT.GT.10) GO TO 140
C
      DO 130 K1=1,N,6
      K2 = MIN(N,K1+5)
      WRITE(LOUT,1002) K1,K2,(SX(I),I=K1,K2)
  130 CONTINUE
      RETURN
C
  140 CONTINUE
      DO 150 K1=1,N,5
      K2 = MIN(N,K1+4)
      WRITE(LOUT,1003) K1,K2,(SX(I),I=K1,K2)
  150 CONTINUE
      RETURN
 1000 FORMAT(1X,I4,' - ',I4,1P,10E12.3)
 1001 FORMAT(1X,I4,' - ',I4,1X,1P,8E14.5)
 1002 FORMAT(1X,I4,' - ',I4,1X,1P,6E18.9)
 1003 FORMAT(1X,I4,' - ',I4,1X,1P,5E24.13)
      END


text:
One persons garbage is another persons treasure.

Dumpster diving is another Davis form of recycling. Many businesses and school organizations throw out wonderful things. Also on moving day in Davis dumpster diving is practically a civic duty. Dumpster diving early and often can help furnish your apartment! Dont pay for furniture that might get ruined quickly. Items as large as couches and as small as wall art are often found in dumpsters at Moving Day end of August/beginning of September. Also see Freshmen Move Out Day.

The The California Aggie Aggie did a front page story called http://www.californiaaggie.com/media/storage/paper981/news/2005/05/16/FrontPage/Diving.For.Dinner1320576.shtml Diving for Dinner on May 15, 2005.

Kind and considerate folks know to leave their usable items besides the dumpster, rather than throw them in, so that others who want it may find them in good shape. Reuse is by far the most efficient form of recycling.

Here is how you do it. Find a dumpster, climb in, and search for something useful. Some amazing things can be found in dumpsters in this town. Many of the grocery stores regularly throw away perfectly good food (some of which is then eaten at Whole Earth Festival meetings).

A notable dumpster diver named Gill used to be seen rummaging the towns dumpsters.

Dumpster Diving is also a classic hacker technique for finding information that leads to an exploit.

dumpsterbating: whereupon an individual finds useful things in their own trash.

 Trash Temples
 Your neighbors dumpster
 Apartment complexes are best because at any time somebody is probably cleaning a closet or moving; middle class neighborhoods seem to generate more stuff than the rich folk.
 Gourmet diving
 Eating food daily is essential, whether youve got the bucks or not. Diving around places that actually produce/resell food versus restaurants will probably yield the best quality, e.g. bakeries, grocery stores. Fast food places regularly toss anything that doesnt sell that day, candy stores sell a product that doesnt really go bad, and movie theaters regularly toss giant bags of popcorn. And if youre diving the burbs, you can score when freezers get cleaned or canned products get tossed on moving day. Diving for food can be risky, so you could instead try hunting for other sources of free food in Davis.
 Bookstores
 When paperbacks and magazines dont sell, the covers get ripped off and sent back to the publisher for a return credit. If youre looking for material for your toilet/library, check behind these places. And dont forget that music stores carry music mags that can be enjoyable potty reading.
 Discount stores
 Cheapie places like The Dollar Tree toss lots of empty boxes from which you can construct a threeseason shelter. Most students opt to hit up mom/dad for rent money, but in a pinch... These places, including Longs, regularly toss slightly damaged/returned goods and plants that theyve very nearly killed.
 Toy/novelty/thrift stores
 Success in this arena is a long shot in Davis, but dont forego the dumpsters behind Davis few high end baby, toy, and novelty stores... if you can drive a real car thats been beat to sht, you can certainly play with a beat up toy car! On this note, dont forget your neighborhood thrift stores or consignment shops  after a while, some things just arent worth the real estate cost for retailers.
 Florist shops/funeral home
 Perhaps little more than a cut above stealing them from the Davis Cemetery or somebodys yard, check out the dumpsters behind one of Davis florist shops for slightly damaged and wilted flowers. If you arent too afraid, look behind the funeral home after services to see whats gotten tossed once the funeral party has gone.
 Copy stores
 Copy shops occasionally toss lots of slightly defective office supplies such as boxes of envelopes that dont quite seal.
 Photo developers
 These dumpsters just make good entertainment. Imperfect and often risqué photos are tossed in damp, chemicalladen bags. Find out what your neighbors have been up to when the curtains are shut!
 Organizations of Evil
 I guess the University might fall within this category, and it certainly does provide ripe pickins for books, interesting papers, old computers/software/lab equipment, and the like. But perhaps more importantly, you can find some great information in dumpsters of organizations such as Davis College Republicans, Womens Center, International House, Sword and Sandals meetings, fraternities, etc.
 Dumpster Ethics If a dumpster is locked, theres probably a reason.

 Good Dumpsters

http://fallingfruit.org/?z13&y38.54635&x121.71932&mtrue&troadmap&lfalse&f843 FallingFruit.org: The global freegan / foraging map now includes all the dumpsters listed below. The map can be edited, much like a wiki, to pinpoint other choice locations that folks might want to share.

  Anderson Plaza dumpsters often have expired food and occasional videos. Dinner and a date!
  Campus dumpsters provide good pickins.
  http://pom.ucdavis.edu/ Department of Pomologys compost bin between Wickson Hall and Wellman Hall is usually filled to the brim with fresh fruits and nuts. Most of it has been experimented on so beware.
    Hardcore, man.  Users/ChrisTakemura ct
  Cost Plus World Market Cost Plus requires its employees to throw out any damaged item, even if its just a scratch, so you may be lucky to find a chair or something out there.
 The Lexington Apartments on Olive Drive have their dumpsters dived really on an hourly basis. This is due to the low income population (actually not the students this time) that live down the block.
  Liquor store dumpsters are worth a look, although they may already be popular with homeless people. Many treasures there to be found for the cashstrapped student or underager (I do not endorse underage drinking). But FYI liquor store dumpsters are also known as urinals.
  The Davis Food CoOp had a good dumpster, now it is under lock and key due to Yolo County health regulations.
  The dumpsters along Drake Drive have had Jennys Furniture free furniture.
  The dumpsters at the North Davis Safeway on Covell generally have better food than the South Davis Safeway on Cowell.

Departed Dumpsters
  Borders dumpster often has stripped books.
    If youre dumpster diving for pornography magazines, youre definitely cheap and pathetic.
      Though it is an interesting way of not financially supporting a questionable industry Users/ArielaHaro
          Are you cheap and pathetic if you did it when you were well under the age of 18 and legally unable to purchase it? If memory serves correctly they used to cycle the girlie mags around the 18th >:D
               Nah, at that point youre just pathetic :)
               Myself, I actually found a foottall stack of porn when I was 12. Angels were singing, I swear.

 Opportunity Knocks
In August and early September of 2004 Sycamore Lane Apartments was gutting all of their apartments in an early phase of the renovation still visible today (2/14/05).  Because it was a furnished apartment complex, this meant they were also scrapping the furniture from all of the apartments.  Their giant commercial dumpsters overflowed with couches, tables, lamps, filing cabinets, and chairs (many of which were still in great condition).  The overflow was not only scattered around the dumpsters, but nearly filled the parking lots as well.  Picture frames, glass cutting boards, mattresses, boxsprings, bed framesthe longer one looked, the more one found.  Once word was out that everything was open game, the ensuing freeforall was a sight to see.  People of all types came from far and wide and took was they could.  In cars, in trucks, on bicycles, and on Walking Distance foot, home furnishings spread throughout the Davis city in all directions.
  But it happened at the expense of people like me that were living there and got thrown out before the lease was up and had nowhere to put stuff.  I would have grabbed some tables and hauled them home but that was home and I had nowhere to put it...  It may have been a nice place to go diving, but it wont be again for a while since people are going to be moving in and not out.  Oh, and dont live there, the management may be nice but the owners are Shadiness Factor SHADY.  Users/MarieHuynh

A friend of mine once found a live, large tarantula still in its plastic tank. Some jerk apparently decided to just throw the furry guy away. Users/DanMasiel (This was a sad find, but we saved it, and named it Tranny Users/MikeyNolan)

Can someone explain to me why a Davis College Republicans DCR dumpster is good? That joke went way over my head and i dont get it. Did someone find something good on DCR from a dumpster?  Users/GeorgejAndrews
   I think it has something to do with gathering secret information on DCR.  Not that very much is very secret about you all, its just that liberals tend to assume that Republicans are all making secret evil plans. Users/BrentLaabs
     Cough, cough... goes the other way, too.  ;)  Ever hear a neocon talk show ranting about the Islamicfacist threat or the Feminazis? (No idea if those are spelled right, not a topic I read about often.  Heh.) Users/JabberWokky JabberWokky the Not Evil, Really

I went to Borders to check out the booktossing miracle. Indeed there are hundreds of coverless books. Unfortunately, however, it seems they didnt just remove the covers. Theyve removed the first 1030 pages. Bummer. Couldnt we just send a few copies of every book to the Davis libraries or perhaps to my house? Users/MichaelGiardina
     The covers are torn off because Borders is destroying them as unsold.  The author and publisher dont get paid for them.  They do not count as sold books  its just more expensive than they are worth to ship them back.  You may feel different, but IMO, if the author isnt getting paid and the publisher is being told they have been destroyed, the ethical thing is to respect the wishes of the author and the actions of the publisher and retailer.  All authors, publishers and book stores do this with paperbacks and magazines.  Its a standard aspect of the agreement between them. Users/JabberWokky
     Do they always tear off some of the first pages, though? I thought it was just the cover, but that doesnt seem to be the case. Also, it appears the magazines are still in perfect order. No rippage. Users/MichaelGiardina
      The amount of damage to the stripped out book/magazine varies. I worked for a major bookstore chain for 5 years and used to do all the receiving/stripping.  While I was doing it, I used to try to be careful that the books and magazines were still in usable shape. I felt bad that it was all just going to end up in a landfill and hoped that some industrious diver could find them.  After a couple of years, a new manager saw what I was doing and had me start tearing the magazines in half and also tearing out the first 2030 pages of the paperbacks. You would be shocked to see the volume of product that even a mid size bookstore destroys just because its too expensive to mail back the whole book. The employees werent allowed to take any of the product home....so it all went in the dumpster.  Its not quite as bad as book burnings...but I still hate to see a book destroyed.  Users/PeterRorvig
     My boyfriend works at Borders and the Mgmt. discovered that books were being taken.They decided to remedy this by ripping the covers off the books and THEN ripping them in half. However, management has recently changed, so the book mutilation may no longer be a policy.Users/StephPolizzi
     What time is best to go diving for Borders?  I definitely would like to check out what kind of magazines are avaliable.Users/AliciaHall
     Im curious about what they do with old Bibles and other Holy Books. I mean, in theory, if Bibles were sitting on a shelf, and were really old, and not sold, would they dare destroy them? Do they ship them back in this case? Do they donate them?Users/MatthewTom
     I currently work at Borders, and can answer your questions.  The only books that are stripped are scifi, romance, and mystery books.  Others are shipped back (including Bibles).  The cover of these books are usually all that is torn off, but sometimes we grab a little too much.  As for the publisher/authors wishes, it really is best to let these books go to the book farm up north.  We have enough trouble with textbook publishers justifying costs, and if these books were recouped, it would only be a matter of time before other publishers caught on an started bringing up book prices.

On August 28th I stopped by a dumpster at one of those apartment complexes off of Alvarado Avenue, but I was beat to the punch by a guy driving a nice Corvette.  As I stepped out of my car he quickly raised his head from inside the dumpster telling me, man, Ive already got this one cleaned out.  Youd better move on.  I had the feeling he was funding his Corvette (or its insurance) by selling his newfound dumpster treasures.

Although dumpster diving for furniture, clothes and other nonperishable items is all good, I find the idea of dumpster diving for food dangerous. If the food has been there for a while, scraping off the fuzzy mold doesnt make the food good. The fuzzy part is only the spores and dispersal structures, the vast majority of the mold is the hyphae that has probably worked its way all the way through the foodstuff. Its also notable that decomposing molds produce very potent toxins, for example, the mold that grows on the shells of peanuts and the outsides of other legumes and nuts, is the source of the second most potent carcinogen known to humans, the number one being dioxin. In addition, flies, such as Calliphora and Muscidae flies, home in on fresh organic matter, and have very promiscuous feeding habits. Calliphorids and Muscids tend to be born and raised in feces, and have been shown to be competent vectors of a wide variety of pathogens, such as dysentery, typhus and other gastrointestinal illnesses. And if that doesnt discourage you, whatever microbes and germs that are already in the dumpster can easily transfer themselves to any foodstuffs inside, and especially if youre diving in the suburbs, realize that people not only toss food, but soiled diapers, used tampons and plenty of household chemicals not fit for consumption in the same trash can. This isnt to say that diving for dinner is going to kill you. Just that it may make you very sick.

my apartment complex has concrete dumpster areas, in this area i found a large drawer to a desk that was equally huge but broken completely in half. Since the drawer was solid wood, not pressboard or anything i dragged it upstairs to my apartment. It also had a wooden drawer divider. I later cleaned it off, sanded and painted it white. It now stands on its side, and makes a lovely night stand. Users/KirstenHaney

Thank you so much Kirsten for pointing that out. I think dumpster diving for (most) clothes, furniture, and other nonfoodstuffs is definitely safe  though some things like soiled clothes are questionable; my mom always chided me on that one when we shopped at thrift stores. People make it sound like dumpster diving for food is all good, when it probably isnt. Businesses probably throw away most food for a reason, like maybe certain foods have expired or has grown mold on it. Dumpster diving for some foodstuffs that are still in its package and are manufactured to have a very long shelflife, like candy, will probably be your safest bet. Other food items are more questionable and I think people should be more discerning. Hey, I know were a wasteful society that makes way too much food for its people and then throws it away  but that does not mean we should eat food that could cause serious illness. If you want to be resourceful with food, there are safer, better ways of doing it. Check out the Free Food page for some better tips. Just because our society is stupid does not mean you have to be stupid. Users/ArielaHaro

I find it disturbing that people think of Manson when they think of Dumpster Divers.  Im insulted. Users/TobinJones

As both a business and dumpster owner and an unrepentant scrounge, I strongly suggest that after a dumpster diving adventure you clean up the area.  Most business owners dont mind a little dumpster recycling but they do mind a big mess around the dumpster.  Put the trash back into the dumpster, put the dumpster back where it was and close the lids.  If you leave junk piled up around the dumpster, it attracts more junk from passersby and the owner has to clean it up.  Nothing encourages an owner to put a lock on a dumpster like having to clean up a mess a couple of times.  Users/grumpyoldgeek

wiki:Ithaca participates in a program called Dump and Run http://www.dumpandrun.org/ where students can leave their good stuff out (at certain locations) and volunteers pick it up.  Then they have a huge rummage sale in the late summer when the new students come back, and all the proceeds go to local charities.  Although I love dumpster diving (and especially loved it in Davis where furniture can sit out for weeks in the summer without fear of rain damage), I am certain that there is lots of great stuff that is not found by divers and goes to the dump...especially after Freshmen Move Out Day.  Dont know how a program like this gets started, but I think Davis could definately benefit!  Users/JanelleAlvstadMattson Janelle

Its great that the area is so lax/supportive of dumpster diving, which is, legally speaking, generally illegal. I am a total diver, and just read this article about dumpster diving where most all of the comments condemned dumpster diving, and it was terrible. I love the amount of recycling going on around Davis though.
Users/DChang

What are the best times for dumpster diving before and during moving day in Davis?
Users/elizabeth2014
 Anytime is best time! more stuff crops up after the move out day but Im already seeing trucks full of obvious scavenging so people are getting out in force ahead of time Daubert


text:
      SUBROUTINE setgmn(meanv,covm,ldcovm,p,parm,ierr)
C      SUBROUTINE setgmn(meanv,covm,p,parm)
C     JJV changed this routine to take leading dimension of COVM
C     JJV argument and pass it to SPOFA, making it easier to use
C     JJV if the COVM which is used is contained in a larger matrix
C     JJV and to make the routine more consistent with LINPACK.
C     JJV Changes are in comments, declarations, and the call to SPOFA.
C**********************************************************************
C
C     SUBROUTINE SETGMN( MEANV, COVM, LDCOVM, P, PARM)
C            SET Generate Multivariate Normal random deviate
C
C
C                              Function
C
C
C      Places P, MEANV, and the Cholesky factoriztion of COVM
C      in PARM for GENMN.
C
C
C                              Arguments
C
C
C     MEANV --> Mean vector of multivariate normal distribution.
C                                        DOUBLE PRECISION MEANV(P)
C
C     COVM   <--> (Input) Covariance   matrix    of  the  multivariate
C                 normal distribution.  This routine uses only the
C                 (1:P,1:P) slice of COVM, but needs to know LDCOVM.
C
C                 (Output) Destroyed on output
C                                        DOUBLE PRECISION COVM(LDCOVM,P)
C
C     LDCOVM --> Leading actual dimension of COVM.
C                                        INTEGER LDCOVM
C
C     P     --> Dimension of the normal, or length of MEANV.
C                                        INTEGER P
C
C     PARM <-- Array of parameters needed to generate multivariate
C                normal deviates (P, MEANV and Cholesky decomposition
C                of COVM).
C                1 : 1                - P
C                2 : P + 1            - MEANV
C                P+2 : P*(P+3)/2 + 1  - Cholesky decomposition of COVM
C                                             DOUBLE PRECISION PARM(P*(P+3)/2 + 1)
C
C**********************************************************************
C     .. Scalar Arguments ..
C      INTEGER p
c      include '../stack.h'
      INTEGER p, ldcovm
C     ..
C     .. Array Arguments ..
C      DOUBLE PRECISION covm(p,p),meanv(p),parm(p* (p+3)/2+1)
      DOUBLE PRECISION covm(ldcovm,p),meanv(p),parm(p* (p+3)/2+1)
C     ..
C     .. Local Scalars ..
      INTEGER i,icount,info,j
C     ..
C     .. External Subroutines ..
      EXTERNAL spofa
C     ..
C     .. Executable Statements ..
C
C
C     TEST THE INPUT
C
   10 parm(1) = p
C
C     PUT P AND MEANV INTO PARM
C
      DO 20,i = 2,p + 1
          parm(i) = meanv(i-1)
   20 CONTINUE
C
C      Cholesky decomposition to find A s.t. trans(A)*(A) = COVM
C
C      CALL spofa(covm,p,p,info)
      CALL spofa(covm,ldcovm,p,info)
      ierr=0
      IF (.NOT. (info.NE.0)) GO TO 30
      call basout(io,wte,"Rand: COV not positive definite")
      ierr=1
      return
   30 icount = p + 1
C
C     PUT UPPER HALF OF A, WHICH IS NOW THE CHOLESKY FACTOR, INTO PARM
C          COVM(1,1) = PARM(P+2)
C          COVM(1,2) = PARM(P+3)
C                    :
C          COVM(1,P) = PARM(2P+1)
C          COVM(2,2) = PARM(2P+2)  ...
C
      DO 50,i = 1,p
          DO 40,j = i,p
              icount = icount + 1
              parm(icount) = covm(i,j)
   40     CONTINUE
   50 CONTINUE
      RETURN
C
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
C     Naive implementation of backpropagation  to aid understanding of the principles.

      program derp
      implicit none

      integer nx, nh, ny, ntrain
      parameter ( nx=2, nh=2, ny=1, ntrain=4 )
      real        x(nx), h(nh), y(ny)

C     a(i,j): weight from node x(j) to node h(i).
C     b(i,j): weight from node h(j) to node y(i).
C     same convention as bprod.m in NNSYMID package (minus bias unit):
C     http://www.iau.dtu.dk/research/control/nnsysid.html

      real a(nh,nx), b(ny,nh)
      real agrad(size(a,1),size(a,2)), bgrad(size(b,1),size(b,2))
      real ydelta(ny), hdelta(nh)
      real xtrain(nx,ntrain), ytrain(ny,ntrain)

C     Setup training examples

      data xtrain / 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0 /
      data ytrain / 0.0, 1.0, 1.0, 0.0 /

      integer i, j, iter, niter
      real    ecost
      real    tdelta, t0, t1
      real    alpha
      parameter ( niter = 100000 )

C     Initialize weights to something random.

      print*,'=================================================='

      print*, 'Initializing random weights.'

      call random_number(a)
      call random_number(b)

      print*,'=================================================='

      call cost(ecost)

      call printmat('cost', reshape([ecost], [1,1]))

      print*,'=================================================='

      print*, 'Numerical gradient with central differences.'

      agrad = 0.0
      bgrad = 0.0

      call ngrad()

      call printmat('dE/da', agrad)
      call printmat('dE/db', bgrad)

      print*,'=================================================='

      print*, 'Gradient with backpropagation.'

      agrad = 0.0
      bgrad = 0.0

      call bpgrad()

      call printmat('dE/da', agrad)
      call printmat('dE/db', bgrad)

      print*,'=================================================='

      print*, 'Train a xor function to a feed forward neural network'

      call random_number(a)
      call random_number(b)

      alpha = 2.0
      iter = 1

      do
        call cost(ecost)
C       call bpgrad()
        call ngrad()
        a = a - alpha*agrad
        b = b - alpha*bgrad

        if (mod(iter, 1000) == 0) print*, 'iter:', iter,
     $                                    'cost:', ecost
        if (ecost < 1e-3) exit
        iter = iter + 1
      end do
      print*, 'iter:', iter, 'cost:', ecost

      print*,'=================================================='

      do i = 1, ntrain
        x = xtrain(:,i)
        call forward
        call printmat('x', reshape(x,[1,nx]))
        call printmat('y', reshape([y],[1,1]))
      end do

      print*,'=================================================='

      call cpu_time(t0)

      do iter = 1, niter
        call ngrad()
      end do

      call cpu_time(t1)
      tdelta = t1 - t0
      print*, 'Central differences ', niter, ' iterations:', tdelta

      call cpu_time(t0)

      do iter = 1, niter
        call bpgrad()
      end do

      call cpu_time(t1)
      tdelta = t1 - t0
      print*, 'Backpropagation ', niter, ' iterations:', tdelta

      print*,'=================================================='

      contains

      subroutine printmat ( name, smat )
        character(len=*)  name
        real              smat(:,:)
        integer           i
        print*, trim(name), ':'
        do i = 1, size(smat,1)
          print*, smat(i,:)
        end do
      end

      real elemental function sigmoid ( x )
        real, intent(in) :: x
        sigmoid = 1.0 / ( 1.0 + exp(-x) )
      end

      subroutine bnorm ( x )
        real        x(:), xmean, xvar, eps
        integer     m
        parameter ( eps = 100.0 * epsilon(1.0) )
        m = size(x)
        xmean = sum(x) / m
        xvar = sum( ( x - xmean )**2 ) / m
        x = ( x - xmean ) / sqrt(xvar + eps)
        x = 0.8*sqrt(xvar)*x + 3.0*xmean
      end

      subroutine forward ()
        h = matmul(a,x)
        call bnorm(h)
        h = sigmoid(h)
        y = matmul(b,h)
        call bnorm(y)
        y = sigmoid(y)
      end

      subroutine cost ( ecost )
        real              ecost
        integer           k
        ecost = 0.0
        do k = 1, ntrain
          x = xtrain(:,k)
          call forward
          ecost = ecost + (0.5 / ntrain) * sum((y-ytrain(:,k))**2)
        end do
      end

C     compute a numerical approximation of the gradient

      subroutine ngrad ()
        real              delta, aij, bij, ep, en
        integer           i, j
        parameter       ( delta = sqrt(epsilon(1.0)) )
C       Approximate dE/da with central differences.
        do j = 1, size(a,2)
          do i = 1, size(a,1)
            aij = a(i,j)
            a(i,j) = aij + delta
            call cost(ep)
            a(i,j) = aij - delta
            call cost(en)
            agrad(i,j) = (ep-en)/(2*delta)
            a(i,j) = aij
          end do
        end do
C       Approximate dE/db with central differences.
        do j = 1, size(b,2)
          do i = 1, size(b,1)
            bij = b(i,j)
            b(i,j) = bij + delta
            call cost(ep)
            b(i,j) = bij - delta
            call cost(en)
            bgrad(i,j) = (ep-en)/(2*delta)
            b(i,j) = bij
          end do
        end do
      end

C     compute gradient with backpropagation NO UPDATE

      subroutine bpgrad ()
        integer           i, j, k
        agrad = 0.0
        bgrad = 0.0
        do k = 1, ntrain
          x = xtrain(:,k)
          call forward
          do i = 1, ny
            ydelta(i) = ( y(i) - ytrain(i,k) ) * y(i) * ( 1.0 - y(i) )
     $                * (1.0 / ntrain)
          end do
C     compute dE/db
          do j = 1, size(b,2)
            do i = 1, size(b,1)
              bgrad(i,j) = bgrad(i,j) + ydelta(i) * h(j)
            end do
          end do
C     compute hdelta
          hdelta = 0.0
          do j = 1, nh
            hdelta(j) = h(j)*(1.0-h(j))*dot_product(ydelta,b(:,j))
C           hdelta(j) = 0.0
C           do i = 1, ny
C             hdelta(j) = hdelta(j)
C    $                  + ydelta(i) * h(j) * ( 1.0 - h(j) ) * b(i,j)
C           end do
          end do
C     compute dE/da
C       This is a rank-1 update of agrad. Could be optimized with BLAS.
          do j = 1, size(a,2)
            do i = 1, size(a,1)
              agrad(i,j) = agrad(i,j) + hdelta(i) * x(j) ! <<< this is right
            end do
          end do
C     end compute dE/da
        end do
      end

      subroutine stocastic
      end

      subroutine batch
      end

      end program derp


text:
subroutine par(xl,xu,pi)
!
!  Subroutine par sets the parameters for the linear
!  PDE problem
!
!  Double precision coding is used
implicit double precision(a-h,o-z)
!
!  Problem parameters
xl=0.0d0
xu=1.0d0
pi=4.0d0*datan(1.0d0)
return
!
!  End of par
end

max_stars_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN

max_issues_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN

max_forks_repo_name: cunyizju/Runge-Kutta-ODE-FORTRAN