text:
subroutine rp1(maxmx, num_eqn, num_waves, num_aux, num_ghost, mx, ql, qr,      &
               auxl, auxr, fwave, s, amdq, apdq)
!======================================================================
!
! Solves normal Riemann problems for the 2D SHALLOW WATER equations
!     with topography:
!     #        h_t + (hu)_x + (hv)_y = 0                           #
!     #        (hu)_t + (hu^2 + 0.5gh^2)_x + (huv)_y = -ghb_x      #
!     #        (hv)_t + (huv)_x + (hv^2 + 0.5gh^2)_y = -ghb_y      #
!! On input, ql contains the state vector at the left edge of each cell
!     qr contains the state vector at the right edge of each cell
!
! This data is along a slice in the x-direction if ixy=1
!     or the y-direction if ixy=2.
!!  Note that the i'th Riemann problem has left state qr(i-1,:)
!     and right state ql(i,:)
!  From the basic clawpack routines, this routine is called with
!     ql = qr
!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                           !
!      # This Riemann solver is for the shallow water equations.            !
!                                                                           !
!       It allows the user to easily select a Riemann solver in             !
!       riemannsolvers_geo.f. this routine initializes all the variables    !
!       for the shallow water equations, accounting for wet dry boundary    !
!       dry cells, wave speeds etc.                                         !
!                                                                           !
!           David George, Vancouver WA, Feb. 2009                           !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    ! Input
    integer, intent(in) :: maxmx, num_eqn, num_waves, num_ghost, num_aux, mx
    real(kind=8), intent(in out) :: ql(num_eqn,1-num_ghost:maxmx+num_ghost)
    real(kind=8), intent(in out) :: qr(num_eqn,1-num_ghost:maxmx+num_ghost)
    real(kind=8), intent(in out) :: auxl(num_aux,1-num_ghost:maxmx+num_ghost)
    real(kind=8), intent(in out) :: auxr(num_aux,1-num_ghost:maxmx+num_ghost)

    ! Output
    real(kind=8), intent(in out) :: s(num_waves, 1-num_ghost:maxmx+num_ghost)
    real(kind=8), intent(in out) :: fwave(num_eqn, num_waves, 1-num_ghost:maxmx+num_ghost)
    real(kind=8), intent(in out) :: amdq(num_eqn, 1-num_ghost:maxmx+num_ghost)
    real(kind=8), intent(in out) :: apdq(num_eqn, 1-num_ghost:maxmx+num_ghost)

    ! Local storage
    integer :: i, n, m

    real(kind=8) :: h_l, h_r, hu_l, hu_r, u_l, u_r, b_l, b_r, phi_l, phi_r

    real(kind=8) :: s1m, s2m, rare1, rare2, h_star

    real(kind=8) :: s_l, s_r, u_hat, c_hat, s_roe_l, s_roe_r, s_einfeldt_l
    real(kind=8) :: s_einfeldt_r

    real(kind=8) :: wall(num_waves), sw(num_waves), fw(num_eqn, num_waves)

    ! Parameters
    integer, parameter :: MAX_ITERATIONS = 1

    ! Common block
    integer :: mcapa
    real(kind=8) :: g, dry_tolerance, earth_radius, deg2rad
    common /cparam/ g, dry_tolerance, earth_radius, deg2rad, mcapa

    fwave = 0.d0
    s = 0.d0
    amdq = 0.d0
    apdq = 0.d0

    ! Primary loop over cell edges
    do i = 2 - num_ghost, mx + num_ghost

        ! Extract state
        h_l = qr(1, i-1)
        h_r = ql(1, i)
        hu_l = qr(2, i-1)
        hu_r = ql(2, i)
        b_l = auxr(1, i-1)
        b_r = auxl(1, i)

        ! Check for negative values
        if (h_l < 0.d0 .or. h_r < 0.d0) then
            print *, "Negative input detected: "
            print *, "    (h_l, h_r, i) = (", h_l, ", ", h_r, ", ", i, ")"
            stop
        end if

        ! Skip problem if completely dry
        if (h_l < dry_tolerance .and. h_r < dry_tolerance) then
            cycle
        end if

        ! Extract derived quantities
        if (h_l > dry_tolerance) then
            u_l = hu_l / h_l
            phi_l = 0.5d0 * g * h_l**2 + hu_l**2 / h_l
        else
            h_l = 0.d0
            u_l = 0.d0
            phi_l = 0.d0
        end if
        if (h_r > dry_tolerance) then
            u_r = hu_r / h_r
            phi_r = 0.5d0 * g * h_r**2 + hu_r**2 / h_r
        else
            h_r = 0.d0
            u_r = 0.d0
            phi_r = 0.d0
        end if

        ! Setup different dry-state problems
        wall = 1.d0
        if (h_r <= dry_tolerance) then
            ! Test to see if the current wave would overcome the wall dry-state
            call riemanntype(h_l, h_l, u_l, -u_l, h_star, s1m, s2m, rare1,    &
                                rare2, 1, dry_tolerance, g)

            if (max(h_l, h_star) + b_l < b_r) then
                ! Right state should become ghost values that mirror left state
                ! so that a wall reflection can happen
                wall(2:) = 0.d0
                h_r = h_l
                hu_r = -hu_l
                b_r = b_l
                phi_r = phi_l
                u_r = -u_l
            else if (h_l + b_l < b_r) then
                ! Special case of the wall bc?
                b_r = h_l + b_l
            end if

        else if (h_l <= dry_tolerance) then
            ! Test to see if the current wave would overcome the wall dry-state
            call riemanntype(h_r, h_r, u_r, -u_r, h_star, s1m, s2m, rare1,    &
                                rare2, 1, dry_tolerance, g)

            if (max(h_r, h_star) + b_r < b_l) then
                ! Right state should become ghost values that mirror left state
                ! so that a wall reflection can happen
                wall(2:) = 0.d0
                h_L = h_r
                hu_L = -hu_r
                b_L = b_r
                phi_L = phi_r
                u_L = -u_r
            else if (h_r + b_r < b_l) then
                ! Special case of the wall bc?
                b_l = h_r + b_r
            end if

        end if

        ! Calculate relevant wave speeds estimates
        ! Basic wave bounds
        s_l = u_l - sqrt(g * h_l)
        s_r = u_r - sqrt(g * h_r)

        ! Roe speeds
        u_hat = (sqrt(g * h_l) * u_l + sqrt(g * h_r) * u_r) / (sqrt(g * h_r) + sqrt(g * h_l))
        c_hat = sqrt(g * 0.5d0 * (h_r + h_l))
        s_roe_l = u_hat - c_hat
        s_roe_r = u_hat + c_hat

        ! Einfeldt speeds
        s_einfeldt_l = min(s_l, s_roe_l)
        s_einfeldt_r = max(s_r, s_roe_r)

        ! Solve Riemann problem
        call riemann_aug_JCP(MAX_ITERATIONS, num_eqn, num_waves, h_l, h_r,  &
                             hu_l, hu_r, 0.d0, 0.d0, b_l, b_r, u_l, u_r,    &
                             0.d0, 0.d0, phi_l, phi_r, s_einfeldt_l,       &
                             s_einfeldt_r, dry_tolerance, g, sw, fw)

        ! Eliminate ghost fluxes for wall
        do m = 1, num_waves
            s(m, i) = sw(m) * wall(m)
            do n = 1, num_eqn
                fwave(n, m, i) = fw(n, m) * wall(m)
            end do
        end do

    end do

    ! Compuate fluctuations
    do i = 2 - num_ghost, mx + num_ghost
        do m = 1, num_waves
            if (s(m, i) < 0.d0) then
                amdq(:, i) = amdq(:, i) + fwave(:, m, i)
            else if (s(m, i) > 0.d0) then
                apdq(:, i) = apdq(:, i) + fwave(:, m, i)
            else
                amdq(:, i) = amdq(:, i) + 0.5d0 * fwave(:, m, i)
                apdq(:, i) = apdq(:, i) + 0.5d0 * fwave(:, m, i)
            end if
        end do
    end do

end subroutine rp1

max_stars_repo_name: mandli/coastal

max_issues_repo_name: mandli/coastal

max_forks_repo_name: mandli/coastal


program test

  ! Import dependencies
  use lapack, only : dgetrf
  use linear_algebra
  use direct_linear_solve, only : dlufactor, dluppfactor

  implicit none

  integer, parameter :: sizes(6) = [200, 400, 800, 1000, 2000, 4000]
  integer :: i, npts

  ! Run for all problem sizes
  do i = 1, size(sizes)

     npts = sizes(i)

     call gesolve(npts)
     call geppsolve(npts)
     !call exact(npts)

  end do

!  call time_report(.false.)
!  call time_report(.true.)

contains

  subroutine geppsolve(npts)

    ! Problem setup
    integer, intent(in) :: npts

    ! Matrices and vectors
    real(8), allocatable, dimension(:,:) :: A, L, U, P
    real(8), allocatable, dimension(:)   :: x, xtmp, b, y

    ! Physics parameters
    real(8), parameter :: lambdas(2) = [0.0d0, 10.0d0]
    real(8) :: lambda
    integer :: flag, i, k

    ! Filename
    character(len=50) :: filename
    character(len=50) :: strlambda, strnpts

    ! Run for each lambda
    do k = 1, size(lambdas)

       lambda = lambdas(k)

       ! Create filename and open
       write(strlambda,*) int(lambda)
       write(strnpts,*) int(npts)
       filename = 'gepp-solution-lam-'// trim(adjustl(strlambda)) &
            & //"-npts-" // trim(adjustl(strnpts)) //'.dat'
       open(11,file=filename)
       write(11,*) "x ", "u"

       ! Allocate matrices
       allocate(A(npts, npts))
       allocate(L, U, P, source = A)
       A = 0.0d0
       L = 0.0d0
       U = 0.0d0
       P = 0.0d0

       ! Allocate vectors
       allocate(x(npts))
       allocate(xtmp, b, y, source = x)
       x    = 0.0d0
       xtmp = 0.0d0
       b    = 0.0d0

       ! Solve x
       call assemble_system(0.0d0, 1.0d0, npts, A, b, x, lambda)
       call dluppfactor(A, L, U, P, flag)
       b = matmul(P,b)
       call fwdsub(L, b, y, flag)
       call backsub(U, y, x, flag)

       ! Write output
       do i = 1, npts
          write(11, *) dble(i)/dble(npts), x(i)
       end do

       close(11)

       deallocate(A,L,U,P,x,xtmp,b,y)

    end do

  end subroutine geppsolve

  pure real(8) function uhat(x, lambda)

    real(8), intent(in) :: x
    real(8), intent(in) :: lambda

    real(8) :: alpha

    alpha = sqrt(lambda)
    uhat =  2.d0*exp(-alpha*x)*(exp(2.d0*alpha)-exp(2.d0*alpha*x))/(exp(2.d0*alpha)-1.d0)

  end function uhat

  subroutine gesolve(npts)

    ! Problem setup
    integer, intent(in) :: npts

    ! Matrices and vectors
    real(8), allocatable, dimension(:,:) :: A, L, U, P
    real(8), allocatable, dimension(:)   :: x, xtmp, b, y

    ! Physics parameters
    real(8), parameter :: lambdas(2) = [0.0d0, 10.0d0]
    real(8) :: lambda
    integer :: flag, i, k

    ! Filename
    character(len=50) :: filename
    character(len=50) :: strlambda, strnpts

    ! Run for each lambda
    do k = 1, size(lambdas)

       lambda = lambdas(k)

       ! Create filename and open
       write(strlambda,*) int(lambda)
       write(strnpts,*) int(npts)
       filename = 'ge-solution-lam-'// trim(adjustl(strlambda)) &
            & //"-npts-" // trim(adjustl(strnpts)) //'.dat'
       open(11,file=filename)
       write(11,*) "x ", "u"

       ! Allocate matrices
       allocate(A(npts, npts))
       allocate(L, U, P, source = A)
       A = 0.0d0
       L = 0.0d0
       U = 0.0d0
       P = 0.0d0

       ! Allocate vectors
       allocate(x(npts))
       allocate(xtmp, b, y, source = x)
       x    = 0.0d0
       xtmp = 0.0d0
       b    = 0.0d0

       ! Solve x
       call assemble_system(0.0d0, 1.0d0, npts, A, b, x, lambda)
       call dlufactor(A, L, U, flag)
       call fwdsub(L, b, y, flag)
       call backsub(U, y, x, flag)

       ! Write output
       do i = 1, npts
          write(11, *) dble(i)/dble(npts), x(i)
       end do

       close(11)

       deallocate(A,L,U,P,x,xtmp,b,y)

    end do

  end subroutine gesolve

  subroutine time_report(pivot)

    ! Problem setup
    integer :: npts
    logical :: pivot

    ! Matrices and vectors
    real(8), allocatable, dimension(:,:) :: A, L, U, P
    real(8), allocatable, dimension(:)   :: x, xtmp, b, y

    ! Timing variables
    real(8) :: total_start, total_finish, total_time
    real(8) :: assembly_start, assembly_finish, assembly_time
    real(8) :: factor_start, factor_finish, factor_time
    real(8) :: elim_start, elim_finish, elim_time

    integer, parameter :: sizes(6) = [200, 400, 800, 1000, 2000, 4000]
    real(8), parameter :: lambdas(2) = [0.0d0, 10.0d0]
    real(8) :: lambda
    integer :: flag, i, k

    ! Run for each lambda
    do k = 1, size(lambdas)

       lambda = lambdas(k)

       ! Write the time for each linear solve
       if (pivot .eqv. .false.) then
          if (k .eq. 1) then
             open(11, file='time_study_ge_lam0.dat')
          else
             open(11, file='time_study_ge_lam10.dat')
          end if
       else
          if (k .eq. 1) then
             open(11, file='time_study_gepp_lam0.dat')
          else
             open(11, file='time_study_gepp_lam10.dat')
          end if
       end if
       write(11, *) "size ", "factor_time ", "elim_time ", "total_time"

       ! Run for all problem sizes
       do i = 1, size(sizes)

          npts = sizes(i)

          ! Allocate matrices
          allocate(A(npts, npts))
          allocate(L, U, P, source = A)
          A = 0.0d0
          L = 0.0d0
          U = 0.0d0
          P = 0.0d0

          ! Allocate vectors
          allocate(x(npts))
          allocate(xtmp, b, y, source = x)
          x    = 0.0d0
          xtmp = 0.0d0
          b    = 0.0d0

          call cpu_time(total_start)

          ! Assemble linear system
          call cpu_time(assembly_start)
          call assemble_system(0.0d0, 1.0d0, npts, A, b, x, lambda)
          call cpu_time(assembly_finish)
          assembly_time = assembly_finish-assembly_start
          print *, "assembled system in ", assembly_time

          call cpu_time(factor_start)
          if (pivot .eqv. .false.) then
             ! Ax = b  is now LUx = b
             call dlufactor(A, L, U, flag)
          else
             call dluppfactor(A, L, U, P, flag)
             b = matmul(P,b)
          end if
          call cpu_time(factor_finish)
          factor_time = factor_finish-factor_start
          print *, "factorized matrix in ", factor_time

          call cpu_time(elim_start)

          ! Solve Ly = b
          call fwdsub(L, b, y, flag)

          ! Solve Ux = y
          call backsub(U, y, x, flag)

          call cpu_time(elim_finish)
          elim_time = elim_finish-elim_start
          print *, "elimination completed in ", elim_time

          call cpu_time(total_finish)
          total_time = total_finish - total_start
          print '("Time = ",f6.3," seconds.")', total_time

          ! Write output data
          write(11, *) npts, factor_time, elim_time, total_time

          deallocate(A,L,U,P,x,xtmp,b,y)

       end do

       close(11)

    end do

  end subroutine time_report

  ! Model problem to solve
  subroutine assemble_system(a, b, npts, V, rhs, u, lambda, homogenous)

    use linear_algebra, only : inv

    real(8), intent(in)  :: a, b ! bound of domain
    integer              :: npts ! number of points
    real(8), intent(out) :: V(npts,npts)
    real(8), intent(out) :: rhs(npts)
    real(8), intent(out) :: u(npts)
    real(8)              :: h
    real(8), intent(in)  :: lambda
    real(8), parameter   :: PI = 3.141592653589793d0
    integer              :: m, n, i, j
    logical, intent(in), optional :: homogenous

    h = (b-a)/dble(npts+1)
    V = 0.0d0

    m = npts
    n = npts
    do i = 1, m
       do j = 1, n
          if (i .eq. j-1) then
             V(i,j) = -1.0d0
          else if (i .eq. j+1) then
             V(i,j) = -1.0d0
          else if (i .eq. j) then
             V(i,i) = 2.0d0 + lambda*h*h
          else
             ! skip zeros
          end if
       end do
    end do

    if(present(homogenous)) then
       ! Assemble the RHS
       rhs = 0.0d0
       rhs(1)    = 2.0d0
       rhs(npts) = 0.0d0
    else
       ! Assemble the RHS
       rhs = 0.0d0
       rhs(1)    = 2.0d0 + h*h*(3.0d0*(dble(i)*h)**2.0d0 - 0.5d0)
       rhs(npts) = 0.0d0 + h*h*(3.0d0*(dble(npts)*h)**2.0d0 - 0.5d0)
    end if

    ! Initial solution profile (straight line)
    do i = 1, npts
       u(i) = 2.0d0 + ((0.0d0+2.0d0)/(b-a))*(dble(i)*h)
    end do

  end subroutine assemble_system

end program test


text:
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 2015 Adrian Martinez Vargas                            %
!                                                                      %
! This software may be modified and distributed under the terms        %
! of the MIT license.  See the LICENSE.txt file for details.           %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
! Junior University.  All rights reserved.                             %
!                                                                      %
! The programs in GSLIB are distributed in the hope that they will be  %
! useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
! responsibility to anyone for the consequences of using them or for   %
! whether they serve any particular purpose or work at all, unless he  %
! says so in writing.  Everyone is granted permission to copy, modify  %
! and redistribute the programs in GSLIB, but only under the condition %
! that this notice and the above copyright notice remain intact.       %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------
!
! The functions and subroutines below were modified from the
! version 2.0 of the gslib code written in fortran 77
!
! The objective is to add functionality to GSLIB and to link
! this code with python using f2py. It uses f90 code convention.
! The arrays are dynamic and externally declared (from python)
!
! the code was converted from Fortran 77 to Fortran 90 using F2F.pl
!
! for more information please refer to:
! - gslib77 source code: http://www.statios.com/software/gslib77_ls.tar.gz
! - GSLIB: Geostatistical Software Library and User's Guide. Second edition
!   by Clayton V. Deutsch, Andre G. Journel, 1997.
! - F2PY Users Guide: http://docs.scipy.org/doc/numpy-dev/f2py/
! - F2F https://bitbucket.org/lemonlab/f2f
!-----------------------------------------------------------------------

subroutine cova3(x1,y1,z1,x2,y2,z2,ivarg,nst,c0,it,cc,aa, &
    irot,maxrot,rotmat,cmax,cova)
    !-----------------------------------------------------------------------

    !                    Covariance Between Two Points
    !                    *****************************

    ! This subroutine calculated the covariance associated with a variogram
    ! model specified by a nugget effect and nested varigoram structures.
    ! The anisotropy definition can be different for each nested structure.

    ! INPUT VARIABLES:

    !   x1,y1,z1         coordinates of first point
    !   x2,y2,z2         coordinates of second point
    !   nst              number of nested structures (same number for all variograms)
    !   ivarg            variogram number (set to 1 unless doing cokriging
    !                       or indicator kriging and in this case use same number of structures in all variograms)
    !   c0(ivarg)        isotropic nugget constant
    !   it(i)            type of each nested structure:
    !                      1. spherical model of range a;
    !                      2. exponential model of parameter a;
    !                           i.e. practical range is 3a
    !                      3. gaussian model of parameter a;
    !                           i.e. practical range is a*sqrt(3)
    !                      4. power model of power a (a must be gt. 0  and
    !                           lt. 2).  if linear model, a=1,c=slope.
    !                      5. hole effect model
    !   cc(i)            multiplicative factor of each nested structure.
    !                      (sill-c0) for spherical, exponential,and gaussian
    !                      slope for linear model.
    !   aa(i)            parameter "a" of each nested structure.
    !   irot             index of the rotation matrix for the first nested
    !                    structure (the second nested structure will use
    !                    irot+1, the third irot+2, and so on)
    !   maxrot           size of rotation matrix arrays
    !   rotmat           rotation matrices
    !
    !  Note that aa, it and cc are 1D arrays with size (MXVARG*MAXNST).
    !  MAXNST was removed from this code and recalculated as MAXROT=MAXNST+1
    !  MXVARG is equal to ivarg
    !
    ! OUTPUT VARIABLES:

    !   cmax             maximum covariance
    !   cova             covariance between (x1,y1,z1) and (x2,y2,z2)

    ! EXTERNAL REFERENCES: sqdist    computes anisotropic squared distance
    !                      rotmat    computes rotation matrix for distance
    !-----------------------------------------------------------------------

    implicit none
    !external references
    real*8 :: sqdist

    ! input
    real, intent(in) :: x1,y1,z1,x2,y2,z2
    integer, intent(in) :: ivarg, irot,maxrot
    integer, intent(in) :: nst
    real, intent(in), dimension(ivarg) :: c0
    integer, intent(in),dimension(nst*ivarg) :: it
    real, intent(in), dimension(nst*ivarg) :: cc, aa
    real*8, intent(in), dimension(maxrot,3,3) :: rotmat

    ! output
    real, intent(out) :: cmax, cova

    ! internal variables
    real ::    hsqd, h, hr
    integer :: ir, is, ist, istart

    !parameters
    real :: DEG2RAD,EPSLON,PI,PMX

    DEG2RAD=3.141592654/180.0
    EPSLON=1.e-20
    PI=3.14159265
    PMX=999.
    EPSLON=1.e-10

    ! Calculate the maximum covariance value (used for zero distances and
    ! for power model covariance):

    istart = 1 + (ivarg-1)*nst
    cmax   = c0(ivarg)
    do is=1,nst
        ist = istart + is - 1
        if(it(ist) == 4) then
            cmax = cmax + PMX
        else
            cmax = cmax + cc(ist)
        endif
    end do

    ! Check for "zero" distance, return with cmax if so:

    hsqd = real(sqdist(x1,y1,z1,x2,y2,z2,irot,maxrot,rotmat))
    if(dble(hsqd) < EPSLON) then
        cova = cmax
        return
    endif

    ! Loop over all the structures:

    cova = 0.0
    do is=1,nst
        ist = istart + is - 1

    ! Compute the appropriate distance:

        if(ist /= 1) then
            ir = min((irot+is-1),maxrot)
            hsqd=real(sqdist(x1,y1,z1,x2,y2,z2,ir,maxrot,rotmat))
        end if
        h = real(sqrt(hsqd))

    ! Spherical Variogram Model?

        if(it(ist) == 1) then
            hr = h/aa(ist)
            if(hr < 1.) cova=cova+cc(ist)*(1.-hr*(1.5-.5*hr*hr))

        ! Exponential Variogram Model?

        else if(it(ist) == 2) then
            cova = cova + cc(ist)*exp(-3.0*h/aa(ist))

        ! Gaussian Variogram Model?

        else if(it(ist) == 3) then
            cova = cova + cc(ist)*exp(-(3.0*h/aa(ist)) &
            *(3.0*h/aa(ist)))

        ! Power Variogram Model?

        else if(it(ist) == 4) then
            cova = cova + cmax - cc(ist)*(h**aa(ist))

        ! Hole Effect Model?

        else if(it(ist) == 5) then
        !                 d = 10.0 * aa(ist)
        !                 cova = cova + cc(ist)*exp(-3.0*h/d)*cos(h/aa(ist)*PI)
            cova = cova + cc(ist)*cos(h/aa(ist)*PI)
        endif
    end do

    ! Finished:

    return

end subroutine cova3

subroutine dcova3(x1,y1,z1,x2,y2,z2,ivarg,nst,c0,it,cc,aa, &
    irot,maxrot,rotmat,cmax,cova)
    !-----------------------------------------------------------------------

    !                    Covariance Between Two Points
    !                    *****************************

    ! This subroutine calculated the covariance associated with a variogram
    ! model specified by a nugget effect and nested varigoram structures.
    ! The anisotropy definition can be different for each nested structure.

    ! INPUT VARIABLES:

    !   x1,y1,z1         coordinates of first point
    !   x2,y2,z2         coordinates of second point
    !   nst              number of nested structures (same number for all variograms)
    !   ivarg            variogram number (set to 1 unless doing cokriging
    !                       or indicator kriging and in this case use same number of structures in all variograms)
    !   c0(ivarg)        isotropic nugget constant
    !   it(i)            type of each nested structure:
    !                      1. spherical model of range a;
    !                      2. exponential model of parameter a;
    !                           i.e. practical range is 3a
    !                      3. gaussian model of parameter a;
    !                           i.e. practical range is a*sqrt(3)
    !                      4. power model of power a (a must be gt. 0  and
    !                           lt. 2).  if linear model, a=1,c=slope.
    !                      5. hole effect model
    !   cc(i)            multiplicative factor of each nested structure.
    !                      (sill-c0) for spherical, exponential,and gaussian
    !                      slope for linear model.
    !   aa(i)            parameter "a" of each nested structure.
    !   irot             index of the rotation matrix for the first nested
    !                    structure (the second nested structure will use
    !                    irot+1, the third irot+2, and so on)
    !   maxrot           size of rotation matrix arrays
    !   rotmat           rotation matrices
    !
    !  Note that aa, it and cc are 1D arrays with size (MXVARG*MAXNST).
    !  MAXNST was removed from this code and recalculated as MAXROT=MAXNST+1
    !  MXVARG is equal to ivarg
    !
    ! OUTPUT VARIABLES:

    !   cmax             maximum covariance
    !   cova             covariance between (x1,y1,z1) and (x2,y2,z2)

    ! EXTERNAL REFERENCES: sqdist    computes anisotropic squared distance
    !                      rotmat    computes rotation matrix for distance
    !-----------------------------------------------------------------------

    implicit none
    !external references
    real*8:: dsqdist

    ! input
    real*8, intent(in) :: x1,y1,z1,x2,y2,z2
    integer, intent(in) :: ivarg, irot,maxrot
    integer, intent(in) :: nst
    real*8, intent(in), dimension(ivarg) :: c0
    integer, intent(in),dimension(nst*ivarg) :: it
    real*8, intent(in), dimension(nst*ivarg) :: cc, aa
    real*8, intent(in), dimension(maxrot,3,3) :: rotmat

    ! output
    real*8, intent(out) :: cmax, cova

    ! internal variables
    real*8 ::    hsqd, h, hr
    integer :: ir, is, ist, istart

    !parameters
    real*8 :: DEG2RAD,EPSLON,PI,PMX

    DEG2RAD=3.141592654/180.0
    EPSLON=1.e-20
    PI=3.14159265
    PMX=999.
    EPSLON=1.e-10

    ! Calculate the maximum covariance value (used for zero distances and
    ! for power model covariance):

    istart = 1 + (ivarg-1)*nst
    cmax   = c0(ivarg)
    do is=1,nst
        ist = istart + is - 1
        if(it(ist) == 4) then
            cmax = cmax + PMX
        else
            cmax = cmax + cc(ist)
        endif
    end do

    ! Check for "zero" distance, return with cmax if so:

    hsqd = dsqdist(x1,y1,z1,x2,y2,z2,irot,maxrot,rotmat)
    if(dble(hsqd) < EPSLON) then
        cova = cmax
        return
    endif

    ! Loop over all the structures:

    cova = 0.0
    do is=1,nst
        ist = istart + is - 1

    ! Compute the appropriate distance:

        if(ist /= 1) then
            ir = min((irot+is-1),maxrot)
            hsqd=dsqdist(x1,y1,z1,x2,y2,z2,ir,maxrot,rotmat)
        end if
        h = dble(dsqrt(hsqd))

    ! Spherical Variogram Model?

        if(it(ist) == 1) then
            hr = h/aa(ist)
            if(hr < 1.) cova=cova+cc(ist)*(1.-hr*(1.5-.5*hr*hr))

        ! Exponential Variogram Model?

        else if(it(ist) == 2) then
            cova = cova + cc(ist)*exp(-3.0*h/aa(ist))

        ! Gaussian Variogram Model?

        else if(it(ist) == 3) then
            cova = cova + cc(ist)*exp(-(3.0*h/aa(ist)) &
            *(3.0*h/aa(ist)))

        ! Power Variogram Model?

        else if(it(ist) == 4) then
            cova = cova + cmax - cc(ist)*(h**aa(ist))

        ! Hole Effect Model?

        else if(it(ist) == 5) then
        !                 d = 10.0 * aa(ist)
        !                 cova = cova + cc(ist)*exp(-3.0*h/d)*cos(h/aa(ist)*PI)
            cova = cova + cc(ist)*cos(h/aa(ist)*PI)
        endif
    end do

    ! Finished:

    return

end subroutine dcova3

max_stars_repo_name: cmrajan/pygslib

max_issues_repo_name: kaufmanno/pygslib

max_forks_repo_name: kaufmanno/pygslib


text:
      logical function common_restr(rst_res,rst1,dag1,rst2,dag2,
     &                        hpvxgas,ngas,nspin)

      implicit none

      include 'opdim.h'
      include 'stdunit.h'

      integer, parameter ::
     &     ntest = 00
      character(len=12), parameter ::
     &     i_am = 'common_restr'

      integer, intent(in) ::
     &     ngas, nspin, hpvxgas(ngas,nspin),
     &     rst1(2,ngas,2,2,nspin),
     &     rst2(2,ngas,2,2,nspin)
      logical, intent(in) ::
     &     dag1, dag2
      integer, intent(out) ::
     &     rst_res(2,ngas,2,2,nspin)

      logical ::
     &     ok
      integer ::
     &     ispin, ica, ica1, ica2

      logical, external ::
     &     is_proper_restr_for_hpvx, zero_ivec

      if (ntest.ge.100) then
        call write_title(lulog,wst_dbg_func,i_am)
        write(lulog,*) 'on entry:'
        call wrt_rstr(lulog,rst1,ngas)
        call wrt_rstr(lulog,rst2,ngas)
      end if

      if (nspin.gt.1) call quit(1,i_am,'nspin.gt.1 -> adapt me fully')

      rst_res = 0
      ok = .true.
      do ispin = 1, nspin
        ! mask = 1: paths included in occupation graph
        do ica = 1, 2
          ica1 = ica
          if (dag1) ica1 = 3-ica
          ica2 = ica
          if (dag2) ica2 = 3-ica
          call common_restr_for_hpvx(rst_res(:,:,ica ,1,ispin),
     &                                  rst1(:,:,ica1,1,ispin),
     &                                  rst2(:,:,ica2,1,ispin),
     &                                  IHOLE,hpvxgas(:,ispin),ngas)
          call common_restr_for_hpvx(rst_res(:,:,ica ,1,ispin),
     &                                  rst1(:,:,ica1,1,ispin),
     &                                  rst2(:,:,ica2,1,ispin),
     &                                  IPART,hpvxgas(:,ispin),ngas)
          call common_restr_for_hpvx(rst_res(:,:,ica ,1,ispin),
     &                                  rst1(:,:,ica1,1,ispin),
     &                                  rst2(:,:,ica2,1,ispin),
     &                                  IVALE,hpvxgas(:,ispin),ngas)
          call common_restr_for_hpvx(rst_res(:,:,ica ,1,ispin),
     &                                  rst1(:,:,ica1,1,ispin),
     &                                  rst2(:,:,ica2,1,ispin),
     &                                  IEXTR,hpvxgas(:,ispin),ngas)

          if (ntest.ge.100) then
            write(lulog,*) 'present ',ispin,ica
            call wrt_rstr(lulog,rst_res,ngas)
          end if

          ok = ok.and.is_proper_restr_for_hpvx(rst_res(:,:,ica,1,ispin),
     &                                      IHOLE,hpvxgas(:,ispin),ngas)
          ok = ok.and.is_proper_restr_for_hpvx(rst_res(:,:,ica,1,ispin),
     &                                      IPART,hpvxgas(:,ispin),ngas)
          ok = ok.and.is_proper_restr_for_hpvx(rst_res(:,:,ica,1,ispin),
     &                                      IVALE,hpvxgas(:,ispin),ngas)
          ok = ok.and.is_proper_restr_for_hpvx(rst_res(:,:,ica,1,ispin),
     &                                      IEXTR,hpvxgas(:,ispin),ngas)

        end do

        ! mask = 2: paths excluded in occupation graph
        if (.not.zero_ivec(rst1(:,:,:,2,ispin),2*ngas*2).or.
     &      .not.zero_ivec(rst2(:,:,:,2,ispin),2*ngas*2)) then
          write(lulog,'(1x,a,16i4)') 'rst1 (mask): ',rst1(:,:,:,2,ispin)
          write(lulog,'(1x,a,16i4)') 'rst2 (mask): ',rst2(:,:,:,2,ispin)
          call quit(1,i_am,'adapt me for mask restrictions !')
        end if

      end do

      common_restr = ok

      if (ntest.ge.100) then
        write(lulog,*) 'on exit: ok = ',ok
        call wrt_rstr(lulog,rst_res,ngas)
      end if

      end

ext: f

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
      SUBROUTINE HIVEL(IH)
*
*
*       High-velocity particle search.
*       ------------------------------
*
      INCLUDE 'common6.h'
*
*
*       Copy membership and add 1 for analogy with HARP/GRAPE.
      NHI = LISTV(1) + 1
*
*       Specify square velocity limit in terms of current state.
      VMAX2 = 16.0*ECLOSE
*
*       Check for removal of distant high-velocity particle.
      IF (IH.LT.0) THEN
    1     LH = 0
          DO 2 L = 2,NHI
              I = LISTV(L)
              VI2 = X0DOT(1,I)**2 + X0DOT(2,I)**2 + X0DOT(3,I)**2
              RI2 = (X(1,I) - RDENS(1))**2 + (X(2,I) - RDENS(2))**2 +
     &                                       (X(3,I) - RDENS(3))**2
*       Save index of fast particle outside 3*<R> or VI < VMAX/4.
              IF (RI2.GT.9.0*RSCALE**2.OR.VI2.LT.VMAX2/16.0) THEN
                  LH = L
                  LI = I
                  RL2 = RI2
                  VL2 = VI2
              END IF
    2     CONTINUE
*       Reduce membership and remove any distant member from the list.
          IF (LH.GT.0) THEN
              if(rank.eq.0)then
              WRITE (29,3)  TIME+TOFF, LI, NAME(LI), SQRT(RL2),
     &                      SQRT(VL2)
    3         FORMAT (' HIVEL REMOVE    T I NAM R V ',F12.4,2I6,2F6.1)
              CALL FLUSH(29)
              end if
              NHI = NHI - 1
              LISTV(1) = LISTV(1) - 1
              DO 4 L = LH,NHI
                  LISTV(L) = LISTV(L+1)
    4         CONTINUE
              GO TO 1
          END IF
      END IF
*
*       Set index to known neutron star or terminated KS/chain.
      IF (IH.GT.0) THEN
          I1 = IH
          I2 = IH
      ELSE
*       Include two first single particles or up to five chain members.
          I1 = IFIRST
          I2 = IFIRST + 1
          IF (IPHASE.EQ.8) THEN
              I2 = IFIRST + 4
*       Search all particles after escaper removal (defined by IPHASE = -1).
          ELSE IF (IH.EQ.0.AND.IPHASE.EQ.-1) THEN
              I2 = NTOT
              NHI = 1
              LISTV(1) = 0
          ELSE IF (IPHASE.EQ.1) THEN
*       See whether the first few locations may have been exchanged.
              DO 5 L = 2,NHI
                  IF (LISTV(L).LT.IFIRST + 3) THEN
                      I2 = NTOT
                  END IF
    5         CONTINUE
          END IF
      END IF
*
*       Add any new high-velocity particles (save F**2 > N & STEP < DTMIN).
      NHV = 0
      DO 10 I = I1,I2
          FI2 = F(1,I)**2 + F(2,I)**2 + F(3,I)**2
          VI2 = X0DOT(1,I)**2 + X0DOT(2,I)**2 + X0DOT(3,I)**2
          RI2 = (X(1,I) - RDENS(1))**2 + (X(2,I) - RDENS(2))**2 +
     &                                   (X(3,I) - RDENS(3))**2
*       Form a list of recently ejected candidates.
          IF (FI2.GT.FLOAT(N).OR.STEP(I).LT.DTMIN.OR.
     &        STEPR(I).LT.20.0*DTMIN) THEN
              IF (VI2.GT.VMAX2) THEN
                  NHV = NHV + 1
                  JLIST(NHV) = I
              END IF
              GO TO 10
          END IF
*      Include all VMAX2 members inside 2*RSCALE or with negative velocity.
          RD = X(1,I)*XDOT(1,I) + X(2,I)*XDOT(2,I) + X(3,I)*XDOT(3,I)
          IF (VI2.GT.VMAX2.AND.(RI2.LT.4.0*RSCALE**2.OR.RD.LT.0.)) THEN
              DO 8 L = 2,NHI
                  IF (I.EQ.LISTV(L)) GO TO 10
    8         CONTINUE
*       Check maximum membership and possible ghost particle.
              IF (NHI.GE.MLV.OR.STEP(I).GT.1.0) GO TO 10
              NHI = NHI + 1
              LISTV(1) = LISTV(1) + 1
              NFAST = NFAST + 1
              LISTV(NHI) = I
              if(rank.eq.0)then
              WRITE (29,9)  TIME+TOFF, NHI-1, I, NAME(I), KSTAR(I),
     &                      SQRT(VI2), SQRT(RI2), STEP(I)
    9         FORMAT (' HIVEL ADD    T NHI I NM K* VI R DT ',
     &                               F12.4,I4,2I8,I4,2F6.2,1P,E10.2)
              CALL FLUSH(29)
              end if
          END IF
   10 CONTINUE
*
*       Consider single fast particle or hyperbolic two-body motion.
      IF ((NHV.EQ.1.OR.NHV.EQ.2).AND.NHI.LE.MLV-NHV) THEN
*       Compare any candidates with existing members.
          DO 20 K = 1,NHV
              DO 15 L = 2,NHI
                  IF (JLIST(K).EQ.LISTV(L)) GO TO 30
   15         CONTINUE
   20     CONTINUE
          I1 = JLIST(1)
*       Include single particles without further tests.
          IF (NHV.EQ.1) THEN
              NHI = NHI + 1
              LISTV(1) = LISTV(1) + 1
              NFAST = NFAST + 1
              LISTV(NHI) = I1
              RI2 = (X(1,I1) - RDENS(1))**2 + (X(2,I1) - RDENS(2))**2 +
     &                                        (X(3,I1) - RDENS(3))**2
              VI2 = X0DOT(1,I1)**2 + X0DOT(2,I1)**2 + X0DOT(3,I1)**2
              if(rank.eq.0)
     &        WRITE (29,22)  TTOT, NHI-1, NAME(I1), IPHASE, SQRT(VI2),
     &                       SQRT(RI2), STEP(I1)
   22         FORMAT (' HIVEL ADD    T NHI NM IPH VI R DT ',
     &                               F12.4,I4,I8,I4,2F6.2,1P,E10.2)
              GO TO 30
          END IF
*       Evaluate two-body energy.
          I2 = JLIST(2)
          RIJ2 = 0.0
          VIJ2 = 0.0
          RDOT = 0.0
          DO 25 K = 1,3
              RIJ2 = RIJ2 + (X(K,I1) - X(K,I2))**2
              VIJ2 = VIJ2 + (XDOT(K,I1) - XDOT(K,I2))**2
              RDOT = RDOT + (X(K,I1) - X(K,I2))*(XDOT(K,I1)-XDOT(K,I2))
   25     CONTINUE
          RIJ = SQRT(RIJ2)
          SEMI = 2.0/RIJ - VIJ2/(BODY(I1) + BODY(I2))
*       Accept outwards hyperbolic motion arising from recent interaction.
          IF (SEMI.LT.0.0.AND.RIJ.LT.10.0*RMIN.AND.RDOT.GT.0.0) THEN
              NHI = NHI + 1
              LISTV(NHI) = I1
              NHI = NHI + 1
              LISTV(NHI) = I2
              LISTV(1) = LISTV(1) + 2
              NFAST = NFAST + 2
              if(rank.eq.0)
     &        WRITE (29,28) TTOT, NHI-1, NAME(I1), NAME(I2), IPHASE, RIJ
   28         FORMAT (' HIVEL ADD    T NHI NM IPH RIJ ',
     &                               F12.4,I4,2I8,I4,1P,E10.2)
          END IF
      END IF
*
   30 CONTINUE
*       Increase counter if loop is over all particles.
*     IF (I2.EQ.NTOT) THEN
*         NHIVEL = NHIVEL + 1
*     END IF
*
      RETURN
*
      END

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
!Program to test SCALE intrinsic function.

program test_scale
  call test_real4 (3.0, 2)
  call test_real4 (33.0, -2)
  call test_real4 (-3., 2)
  call test_real4 (0., 3)
  call test_real8 (0._8, 3)
  call test_real8 (3.0_8, 4)
  call test_real8 (33.0_8, -4)
  call test_real8 (-33._8, 4)
end
subroutine test_real4 (orig, i)
  real x,y,orig
  integer i
  x = orig
  y = x * (2.0 ** i)
  x = scale (x, i)
  if (abs (x - y) .gt. abs(x * 1e-6)) call abort
end

subroutine test_real8 (orig, i)
  real*8 x,y,orig
  integer i
  x = orig
  y = x * (2.0 ** i)
  x = scale (x, i)
  if (abs (x - y) .gt. abs(x * 1e-6)) call abort
end

max_issues_repo_name: vidkidz/crossbridge

max_forks_repo_name: vidkidz/crossbridge


text:
! sample compile command to generate *.o object file:
!     gfortran -fopenmp -c Fig_8.3_regPromote.f90

program main
   use omp_lib
   implicit none

   interface
      function doit(A, N, id)
         integer :: N, id
         real*8 :: A(N)
         real*8 :: doit
      end function
   end interface

   integer, parameter :: MAX = 10000
   integer, parameter :: NMAX = 1000
   real, parameter :: TOL = 0.0001

   integer :: iter, N
   real*8 :: A(NMAX)
   real*8 :: conv
   integer :: id, nthrd

   iter = 0
   N = 1000
   A = 0.0
   conv = 0.0

   !$omp parallel shared(A,N,iter) firstprivate(conv) private(id,nthrd)
      id = omp_get_thread_num()
      nthrd = omp_get_num_threads()

      do while (iter < MAX)
         conv = doit(A, N, id)
         if (conv < TOL) exit
         if (id == 0) iter = iter + 1
      end do

   !$omp end parallel

end program main

max_stars_repo_name: zafar-hussain/OmpCommonCore

max_issues_repo_name: zafar-hussain/OmpCommonCore

max_forks_repo_name: zafar-hussain/OmpCommonCore


text:
      SUBROUTINE GMBDRD
C
      IMPLICIT REAL*8(A-H,O-Z)
C
      COMMON/CGIMBD/ELIN(3,3),ELAX(3),ELCG(3),ELMS,ZTZT(3,3)
C
      COMMON/DMBICS/ELEVI,ELEVID,GMUP(2),GMDN(2)
C
      COMMON/DMINTF/GMK1(2),GMK2(2),GMDMP(2),GMSTP(2)
C
      COMMON/DMPRPL/ GMBAZ,GMBAZD,GMBEL,GMBELD
C
      COMMON/ELKDMP/ OMKDMP(3,10),IOMKDM(10)
C
      COMMON/IGIMBD/IGMBL,NELEV,NE1
C
      DIMENSION HEDGM(5)
      DIMENSION HEDIC(5)
      DIMENSION HEDJD(5)
      REAL*4 BUFF(450)
C
      DATA I8/',A8,'/
      DATA HEDGM/'    ELEV',' GIMBLE ','PHYSICAL',' PROPERT','IES     '/
      DATA HEDIC/'GIMBLE I','NITIAL C','ONDITION','S       ','        '/
      DATA HEDJD/'JET DAMP','ING OPTI','ON INVOK','ED      ','        '/
C
      CALL SETUP(8HDELIN   ,8,ELIN,3,3)
      CALL SETUP(8HDELAX   ,8,ELAX,3)
      CALL SETUP(8HDELCG   ,8,ELCG,3)
      CALL SETUP(8HDELMS   ,8,ELMS)
C
      CALL SETUP(8HDLEVI   ,8,ELEVI)
      CALL SETUP(8HDLEVID  ,8,ELEVID)
      CALL SETUP(8HDMUP    ,8,GMUP,2)
      CALL SETUP(8HDMDN    ,8,GMDN,2)
C
      CALL SETUP(8HDMK1    ,8,GMK1,2)
      CALL SETUP(8HDMK2    ,8,GMK2,2)
      CALL SETUP(8HDMDMP   ,8,GMDMP,2)
      CALL SETUP(8HDMSTP   ,8,GMSTP,2)
C
      CALL SETUP(8HIGMBLD  ,4,IGMBL)
C
      CALL SETUP(8HOMKDMP  ,8,OMKDMP,3,10)
      CALL SETUP(8HIOMKDM  ,4,IOMKDM,10)
C
      CALL DAMPCS
C
C
      RETURN
C
C    ***************************************************************
      ENTRY ECHDGM
C    ***************************************************************
C
      IF(IGMBL.EQ.0) GO TO 10
C
      CALL HVAL(HEDGM)
C
      CALL FVAL('ELIN    ',4,ELIN,3,3,2)
      CALL FVAL('ELAX    ',4,ELAX,3,0,1)
      CALL FVAL('ELCG    ',4,ELCG,3,0,1)
      CALL FVAL('ELMS    ',4,ELMS,0,0,0)
C
      CALL HVAL(HEDIC)
C
      CALL FVAL('ELEVI   ',5,ELEVI,0,0,0)
      CALL FVAL('ELEVID  ',6,ELEVID,0,0,0)
C
   10 CONTINUE
C
      RETURN
C
C    ***************************************************************
      ENTRY DMBPLT(BUFF,INDX)
C    ***************************************************************
C
      INDEX=INDX-1
      INDX=INDX+4
C
      IF(IGMBL.EQ.0) RETURN
C
      BUFF(INDEX+1)=GMBAZ
      BUFF(INDEX+2)=GMBAZD
      BUFF(INDEX+3)=GMBEL
      BUFF(INDEX+4)=GMBELD
C
      RETURN
C
C    ***************************************************************
      ENTRY DMPRNT
C    ***************************************************************
C
      IF(IGMBL.EQ.0) RETURN
C
      CALL SET('GMBL AZ ',0,0,GMBAZ,I8)
      CALL SET('GMBL AZD',0,0,GMBAZD,I8)
      CALL SET('GMBL EL ',0,0,GMBEL,I8)
      CALL SET('GMBL ELD',0,0,GMBELD,I8)
C
C
      RETURN
C
C
C
      END

ext: for


text:
!
! Copyright (C) 2001-2007 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
#define ZERO ( 0.D0, 0.D0 )
#define ONE  ( 1.D0, 0.D0 )
!
!----------------------------------------------------------------------------
SUBROUTINE ccgdiagg( npwx, npw, nbnd, npol, psi, e, btype, precondition, &
                     ethr, maxter, reorder, notconv, avg_iter )
  !----------------------------------------------------------------------------
  !
  ! ... "poor man" iterative diagonalization of a complex hermitian matrix
  ! ... through preconditioned conjugate gradient algorithm
  ! ... Band-by-band algorithm with minimal use of memory
  ! ... Calls h_1psi and s_1psi to calculate H|psi> and S|psi>
  ! ... Works for generalized eigenvalue problem (US pseudopotentials) as well
  !
  USE constants,        ONLY : pi
  USE kinds,            ONLY : DP
  USE mp_global,        ONLY : intra_bgrp_comm
  USE mp,               ONLY : mp_sum
  !
  IMPLICIT NONE
  !
  ! ... I/O variables
  !
  INTEGER,     INTENT(IN)    :: npwx, npw, nbnd, npol, maxter
  INTEGER,     INTENT(IN)    :: btype(nbnd)
  REAL(DP),    INTENT(IN)    :: precondition(npwx*npol), ethr
  COMPLEX(DP), INTENT(INOUT) :: psi(npwx*npol,nbnd)
  REAL(DP),    INTENT(INOUT) :: e(nbnd)
  INTEGER,     INTENT(OUT)   :: notconv
  REAL(DP),    INTENT(OUT)   :: avg_iter
  !
  ! ... local variables
  !
  INTEGER                  :: i, j, m, iter, moved
  COMPLEX(DP), ALLOCATABLE :: hpsi(:), spsi(:), lagrange(:), &
                              g(:), cg(:), scg(:), ppsi(:), g0(:)
  REAL(DP)                 :: psi_norm, a0, b0, gg0, gamma, gg, gg1, &
                              cg0, e0, es(2)
  REAL(DP)                 :: theta, cost, sint, cos2t, sin2t
  LOGICAL                  :: reorder
  INTEGER                  :: kdim, kdmx, kdim2
  REAL(DP)                 :: empty_ethr, ethr_m
  !
  ! ... external functions
  !
  REAL (DP), EXTERNAL :: ddot
  !
  !
  CALL start_clock( 'ccgdiagg' )
  !
  empty_ethr = MAX( ( ethr * 5.D0 ), 1.D-5 )
  !
  IF ( npol == 1 ) THEN
     !
     kdim = npw
     kdmx = npwx
     !
  ELSE
     !
     kdim = npwx * npol
     kdmx = npwx * npol
     !
  END IF
  !
  kdim2 = 2 * kdim
  !
  ALLOCATE( spsi( kdmx ) )
  ALLOCATE( scg(  kdmx ) )
  ALLOCATE( hpsi( kdmx ) )
  ALLOCATE( g(    kdmx ) )
  ALLOCATE( cg(   kdmx ) )
  ALLOCATE( g0(   kdmx ) )
  ALLOCATE( ppsi( kdmx ) )
  !
  ALLOCATE( lagrange( nbnd ) )
  !
  avg_iter = 0.D0
  notconv  = 0
  moved    = 0
  !
  ! ... every eigenfunction is calculated separately
  !
  DO m = 1, nbnd
     !
     IF ( btype(m) == 1 ) THEN
        !
        ethr_m = ethr
        !
     ELSE
        !
        ethr_m = empty_ethr
        !
     END IF
     !
     spsi     = ZERO
     scg      = ZERO
     hpsi     = ZERO
     g        = ZERO
     cg       = ZERO
     g0       = ZERO
     ppsi     = ZERO
     lagrange = ZERO
     !
     ! ... calculate S|psi>
     !
     CALL s_1psi( npwx, npw, psi(1,m), spsi )
     !
     ! ... orthogonalize starting eigenfunction to those already calculated
     !
     CALL ZGEMV( 'C', kdim, m, ONE, psi, kdmx, spsi, 1, ZERO, lagrange, 1 )
     !
     CALL mp_sum( lagrange( 1:m ), intra_bgrp_comm )
     !
     psi_norm = DBLE( lagrange(m) )
     !
     DO j = 1, m - 1
        !
        psi(:,m)  = psi(:,m) - lagrange(j) * psi(:,j)
        !
        psi_norm = psi_norm - &
                   ( DBLE( lagrange(j) )**2 + AIMAG( lagrange(j) )**2 )
        !
     END DO
     !
     psi_norm = SQRT( psi_norm )
     !
     psi(:,m) = psi(:,m) / psi_norm
     !
     ! ... calculate starting gradient (|hpsi> = H|psi>) ...
     !
     CALL h_1psi( npwx, npw, psi(1,m), hpsi, spsi )
     !
     ! ... and starting eigenvalue (e = <y|PHP|y> = <psi|H|psi>)
     !
     ! ... NB:  ddot(2*npw,a,1,b,1) = REAL( zdotc(npw,a,1,b,1) )
     !
     e(m) = ddot( kdim2, psi(1,m), 1, hpsi, 1 )
     !
     CALL mp_sum( e(m), intra_bgrp_comm )
     !
     ! ... start iteration for this band
     !
     iterate: DO iter = 1, maxter
        !
        ! ... calculate  P (PHP)|y>
        ! ... ( P = preconditioning matrix, assumed diagonal )
        !
        g(:)    = hpsi(:) / precondition(:)
        ppsi(:) = spsi(:) / precondition(:)
        !
        ! ... ppsi is now S P(P^2)|y> = S P^2|psi>)
        !
        es(1) = ddot( kdim2, spsi(1), 1, g(1), 1 )
        es(2) = ddot( kdim2, spsi(1), 1, ppsi(1), 1 )
        !
        CALL mp_sum( es , intra_bgrp_comm )
        !
        es(1) = es(1) / es(2)
        !
        g(:) = g(:) - es(1) * ppsi(:)
        !
        ! ... e1 = <y| S P^2 PHP|y> / <y| S S P^2|y> ensures that
        ! ... <g| S P^2|y> = 0
        ! ... orthogonalize to lowest eigenfunctions (already calculated)
        !
        ! ... scg is used as workspace
        !
        CALL s_1psi( npwx, npw, g(1), scg(1) )
        !
        CALL ZGEMV( 'C', kdim, ( m - 1 ), ONE, psi, &
                    kdmx, scg, 1, ZERO, lagrange, 1  )
        !
        CALL mp_sum( lagrange( 1:m-1 ), intra_bgrp_comm )
        !
        DO j = 1, ( m - 1 )
           !
           g(:)   = g(:)   - lagrange(j) * psi(:,j)
           scg(:) = scg(:) - lagrange(j) * psi(:,j)
           !
        END DO
        !
        IF ( iter /= 1 ) THEN
           !
           ! ... gg1 is <g(n+1)|S|g(n)> (used in Polak-Ribiere formula)
           !
           gg1 = ddot( kdim2, g(1), 1, g0(1), 1 )
           !
           CALL mp_sum( gg1, intra_bgrp_comm )
           !
        END IF
        !
        ! ... gg is <g(n+1)|S|g(n+1)>
        !
        g0(:) = scg(:)
        !
        g0(:) = g0(:) * precondition(:)
        !
        gg = ddot( kdim2, g(1), 1, g0(1), 1 )
        !
        CALL mp_sum( gg, intra_bgrp_comm )
        !
        IF ( iter == 1 ) THEN
           !
           ! ... starting iteration, the conjugate gradient |cg> = |g>
           !
           gg0 = gg
           !
           cg(:) = g(:)
           !
        ELSE
           !
           ! ... |cg(n+1)> = |g(n+1)> + gamma(n) * |cg(n)>
           !
           ! ... Polak-Ribiere formula :
           !
           gamma = ( gg - gg1 ) / gg0
           gg0   = gg
           !
           cg(:) = cg(:) * gamma
           cg(:) = g + cg(:)
           !
           ! ... The following is needed because <y(n+1)| S P^2 |cg(n+1)>
           ! ... is not 0. In fact :
           ! ... <y(n+1)| S P^2 |cg(n)> = sin(theta)*<cg(n)|S|cg(n)>
           !
           psi_norm = gamma * cg0 * sint
           !
           cg(:) = cg(:) - psi_norm * psi(:,m)
           !
        END IF
        !
        ! ... |cg> contains now the conjugate gradient
        !
        ! ... |scg> is S|cg>
        !
        CALL h_1psi( npwx, npw, cg(1), ppsi(1), scg(1) )
        !
        cg0 = ddot( kdim2, cg(1), 1, scg(1), 1 )
        !
        CALL mp_sum(  cg0 , intra_bgrp_comm )
        !
        cg0 = SQRT( cg0 )
        !
        ! ... |ppsi> contains now HP|cg>
        ! ... minimize <y(t)|PHP|y(t)> , where :
        ! ...                         |y(t)> = cos(t)|y> + sin(t)/cg0 |cg>
        ! ... Note that  <y|P^2S|y> = 1, <y|P^2S|cg> = 0 ,
        ! ...           <cg|P^2S|cg> = cg0^2
        ! ... so that the result is correctly normalized :
        ! ...                           <y(t)|P^2S|y(t)> = 1
        !
        a0 = 2.D0 * ddot( kdim2, psi(1,m), 1, ppsi(1), 1 ) / cg0
        !
        CALL mp_sum(  a0 , intra_bgrp_comm )
        !
        b0 = ddot( kdim2, cg(1), 1, ppsi(1), 1 ) / cg0**2
        !
        CALL mp_sum(  b0 , intra_bgrp_comm )
        !
        e0 = e(m)
        !
        theta = 0.5D0 * ATAN( a0 / ( e0 - b0 ) )
        !
        cost = COS( theta )
        sint = SIN( theta )
        !
        cos2t = cost*cost - sint*sint
        sin2t = 2.D0*cost*sint
        !
        es(1) = 0.5D0 * (   ( e0 - b0 ) * cos2t + a0 * sin2t + e0 + b0 )
        es(2) = 0.5D0 * ( - ( e0 - b0 ) * cos2t - a0 * sin2t + e0 + b0 )
        !
        ! ... there are two possible solutions, choose the minimum
        !
        IF ( es(2) < es(1) ) THEN
           !
           theta = theta + 0.5D0 * pi
           !
           cost = COS( theta )
           sint = SIN( theta )
           !
        END IF
        !
        ! ... new estimate of the eigenvalue
        !
        e(m) = MIN( es(1), es(2) )
        !
        ! ... upgrade |psi>
        !
        psi(:,m) = cost * psi(:,m) + sint / cg0 * cg(:)
        !
        ! ... here one could test convergence on the energy
        !
        IF ( ABS( e(m) - e0 ) < ethr_m ) EXIT iterate
        !
        ! ... upgrade H|psi> and S|psi>
        !
        spsi(:) = cost * spsi(:) + sint / cg0 * scg(:)
        !
        hpsi(:) = cost * hpsi(:) + sint / cg0 * ppsi(:)
        !
     END DO iterate
     !
     IF ( iter >= maxter ) notconv = notconv + 1
     !
     avg_iter = avg_iter + iter + 1
     !
     ! ... reorder eigenvalues if they are not in the right order
     ! ... ( this CAN and WILL happen in not-so-special cases )
     !
     IF ( m > 1 .AND. reorder ) THEN
        !
        IF ( e(m) - e(m-1) < - 2.D0 * ethr_m ) THEN
           !
           ! ... if the last calculated eigenvalue is not the largest...
           !
           DO i = m - 2, 1, - 1
              !
              IF ( e(m) - e(i) > 2.D0 * ethr_m ) EXIT
              !
           END DO
           !
           i = i + 1
           !
           moved = moved + 1
           !
           ! ... last calculated eigenvalue should be in the
           ! ... i-th position: reorder
           !
           e0 = e(m)
           !
           ppsi(:) = psi(:,m)
           !
           DO j = m, i + 1, - 1
              !
              e(j) = e(j-1)
              !
              psi(:,j) = psi(:,j-1)
              !
           END DO
           !
           e(i) = e0
           !
           psi(:,i) = ppsi(:)
           !
           ! ... this procedure should be good if only a few inversions occur,
           ! ... extremely inefficient if eigenvectors are often in bad order
           ! ... ( but this should not happen )
           !
        END IF
        !
     END IF
     !
  END DO
  !
  avg_iter = avg_iter / DBLE( nbnd )
  !
  DEALLOCATE( lagrange )
  DEALLOCATE( ppsi )
  DEALLOCATE( g0 )
  DEALLOCATE( cg )
  DEALLOCATE( g )
  DEALLOCATE( hpsi )
  DEALLOCATE( scg )
  DEALLOCATE( spsi )
  !
  CALL stop_clock( 'ccgdiagg' )
  !
  RETURN
  !
END SUBROUTINE ccgdiagg

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
The MichelNovelozo Gallery is an Art Galleries art gallery in Davis. It displays Photography photographic works in the front room and basement. It replaced the JGlenn Gallery.

APEX Cycles used to be located in their basement but in 2008 moved to their own building on G Street.

20050930 15:16:57 Do they have a website, or email address? Users/GeoffJohnson
   The site listed is from the Davis Art Association brochure, however it is down as of 20051106.

20061026 01:49:19 nbsp I went there once on gallery night and saw some amazing pinholes!  I also saw my high school zoology/botany teacher out with his daughter.  The gallery had cookies that were a whole buncha delish!! Users/MatthiasGropp

20110926 10:39:09 nbsp The Gallery closed in 2008 and Neil Michel moved his Photography business to A private studio off of G Street. Users/Aaron.Curtin

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine cudaclaw5_set_block(blockno)
      implicit none

      integer blockno, blockno_com
      common /comblock/ blockno_com

      blockno_com = blockno
      end

      integer function fc2d_cudaclaw5_get_block()
      implicit none

      integer blockno_com
      common /comblock/ blockno_com

      fc2d_cudaclaw5_get_block = blockno_com
      return
      end

      subroutine cudaclaw5_unset_block()
      implicit none

      integer blockno, blockno_com
      common /comblock/ blockno_com

      blockno_com = -1
      end

max_stars_repo_name: ECLAIRWaveS/ForestClaw

max_issues_repo_name: ECLAIRWaveS/ForestClaw

max_forks_repo_name: ECLAIRWaveS/ForestClaw


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

contains

subroutine taylor(mg,system,info,stencil,srg,tspsi_in,tspsi_out,sshtpsi,   &
                  ppg,V_local,rt)
  use inputoutput, only: n_hamil
  use structures
  use hamiltonian, only: hpsi
  use sendrecv_grid, only: s_sendrecv_grid
  implicit none
  type(s_rgrid),intent(in) :: mg
  type(s_dft_system),intent(in) :: system
  type(s_parallel_info),intent(in) :: info
  type(s_stencil),intent(in) :: stencil
  type(s_sendrecv_grid),intent(inout) :: srg
  type(s_orbital),intent(inout) :: tspsi_in
  type(s_orbital),intent(inout) :: tspsi_out
  type(s_orbital),intent(inout) :: sshtpsi
  type(s_pp_grid),intent(in) :: ppg
  type(s_scalar) ,intent(in) :: V_local(system%nspin)
  type(s_rt),     intent(in) :: rt
  integer :: nn,ix,iy,iz
  integer :: ik,io,is,nspin
  nspin = system%nspin

  do nn=1,n_hamil
    if(mod(nn,2)==1)then
      call hpsi(tspsi_in,sshtpsi,info,mg,V_local,system,stencil,srg,ppg)
      if (nn==1) then
!$OMP parallel do collapse(5) private(ik,io,is,iz,iy,ix)
        do ik=info%ik_s,info%ik_e
        do io=info%io_s,info%io_e
          do is=1,nspin
            do iz=mg%is(3),mg%ie(3)
            do iy=mg%is(2),mg%ie(2)
            do ix=mg%is(1),mg%ie(1)
              tspsi_out%zwf(ix,iy,iz,is,io,ik,1)=tspsi_in%zwf(ix,iy,iz,is,io,ik,1)+ &
                                                   rt%zc(nn)*sshtpsi%zwf(ix,iy,iz,is,io,ik,1)
            end do
            end do
            end do
          end do
        end do
        end do
      else
!$OMP parallel do collapse(5) private(ik,io,is,iz,iy,ix)
        do ik=info%ik_s,info%ik_e
        do io=info%io_s,info%io_e
          do is=1,nspin
            do iz=mg%is(3),mg%ie(3)
            do iy=mg%is(2),mg%ie(2)
            do ix=mg%is(1),mg%ie(1)
              tspsi_out%zwf(ix,iy,iz,is,io,ik,1)=tspsi_out%zwf(ix,iy,iz,is,io,ik,1)+ &
                                                   rt%zc(nn)*sshtpsi%zwf(ix,iy,iz,is,io,ik,1)
            end do
            end do
            end do
          end do
        end do
        end do
      end if
    else
      call hpsi(sshtpsi,tspsi_in,info,mg,V_local,system,stencil,srg,ppg)
!$OMP parallel do collapse(5) private(ik,io,is,iz,iy,ix)
      do ik=info%ik_s,info%ik_e
      do io=info%io_s,info%io_e
        do is=1,nspin
          do iz=mg%is(3),mg%ie(3)
          do iy=mg%is(2),mg%ie(2)
          do ix=mg%is(1),mg%ie(1)
            tspsi_out%zwf(ix,iy,iz,is,io,ik,1)=tspsi_out%zwf(ix,iy,iz,is,io,ik,1)+ &
                                               rt%zc(nn)*tspsi_in%zwf(ix,iy,iz,is,io,ik,1)
          end do
          end do
          end do
        end do
      end do
      end do
    end if
  end do

end subroutine taylor

end module taylor_sub


text:
program blocks
! program to check parsing of block structures
  real x
  integer n
  x = bar(1.0)
  if( x .gt. 0 ) then
     call baz(x)
     n = int(x)
     select case (n)
        case(1)
           call baz(2.0)
        case(2:10)
           call baz(3.0)
        case(-1,11:100,300)
           call baz(4.0)
        case default
           call baz(5.0)
     end select
  else
     k = 0
     loop1: do i=1,100
        loop2:    do j=1,10
           write(*,*) i,j,i*j
           if( i .eq. 50 ) exit loop1
           k = k + j
           if( k .eq. 50) cycle loop2
        end do loop2
        write(*,*) k*i**2
     enddo loop1
     case2: select case (n)
        case(1) case2
           call baz(2.0)
        case(2:10) case2
           call baz(3.0)
        case(-1,11:100,300) case2
           call baz(4.0)
        case default case2
           do i=1,10
              if( i == 7 ) goto 200
              call baz(5.0)
200        end do
     end select case2
  end if
end program blocks
function bar(c)
  real a, b, c ,d
  read *, a
! This block is from section 8.1.2.3 of the F90 standard, except for
! removing space between some keywords
  if ( a .gt. 0 ) then
     b = c/a
     if (b .gt. 0) then
        d = 1.0
     endif
  elseif (c .gt. 0) then
     b = a/c
     d = -1.0
  else
     b = abs (max (a, c))
     d = 0
  endif
  bar = d*b
end function bar
subroutine baz(c)
  real a, b, c ,d
  read *, a
! Same as above but with spaces restored and construct names added
  first_if: if ( a .gt. 0 ) then
     b = c/a
     second_if: if (b .gt. 0) then
        d = 1.0
     end if second_if
  else if (c .gt. 0) then first_if
     b = a/c
     d = -1.0
  else first_if
     b = abs (max (a, c))
     d = 0
  end if first_if
  print *, a, b, c, d
end subroutine baz

max_stars_repo_name: kurtsansom/ftnchek


text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
SUBMODULE (class_exodus) class_exodus_procedures

    USE iso_fortran_env, ONLY : INT64, REAL64
    USE class_psblas
    USE class_cell
    USE class_face
    USE class_vertex
    USE class_mesh
    USE class_connectivity
    USE class_iterating
    USE class_scalar_field
    USE class_vector_field
    USE class_vector
    USE class_output

    IMPLICIT NONE
    INCLUDE 'exodusII.inc'

    !! author: Hari Radhakrishnan, GSE
    !! date 12/12/2019
    !!
    !! This module implements the routines for writing an exodus file with mesh and results
    CONTAINS

        MODULE PROCEDURE write_exo_morfeus

            CHARACTER(len=100) :: filename
            INTEGER, DIMENSION(:), ALLOCATABLE :: cell_ids
            INTEGER, DIMENSION(:), ALLOCATABLE :: v2cconn
            INTEGER, DIMENSION(:), ALLOCATABLE :: icverts
            INTEGER, DIMENSION(:), ALLOCATABLE :: iproc
            REAL(psb_dpk_), DIMENSION(:,:), ALLOCATABLE :: points

            INTEGER(kind = int64) :: ierr, titlelen
            INTEGER(kind = int64) :: exodus_file_id, num_dims, num_nodes, num_elems,  &
              num_elem_blks, num_node_sets, num_side_sets
            INTEGER(kind = int64) :: cpu_ws, io_ws, num_props, prop_value
            REAL :: vers
            REAL(kind = real64), ALLOCATABLE :: xx(:), yy(:), zz(:)
            INTEGER :: info, err_act
            INTEGER :: icontxt, mypnum
            INTEGER :: i, ig, ic, nc2g, nv2c, nconn, ncells, ngc, ngroups
            INTEGER(kind = int64), ALLOCATABLE :: conn(:)
            CHARACTER(len = MXSTLN) :: coord_names(3), ctype
            INTEGER, POINTER :: ic2g(:) => NULL()
            INTEGER, POINTER :: iv2c(:) => NULL()
            INTEGER, ALLOCATABLE  :: igroup(:)
            INTEGER, ALLOCATABLE :: i_loc(:)
            TYPE(cell), ALLOCATABLE :: cells(:)
            TYPE(face), ALLOCATABLE :: faces(:)
            TYPE(vertex), ALLOCATABLE :: verts(:)
            TYPE(connectivity) :: v2f, v2c, f2c
            CHARACTER(len=32) :: path
            REAL(psb_dpk_), ALLOCATABLE :: scalar_local(:), scalar_global(:), f_s(:)
            TYPE(vector), ALLOCATABLE :: vector_local(:), vector_global(:)
            REAL(psb_dpk_), ALLOCATABLE :: f_x(:), f_y(:), f_z(:)
            ! Sets error handling for PSBLAS-2 routines
            CALL psb_erractionsave(err_act)
            mypnum  = mypnum_()
            icontxt = icontxt_()

            ! Sets output path
            ! IF(PRESENT(iter)) CALL out%set_output_path(iter)
            ! path = out%path_()
            filename = 'out_nemo.e' !TRIM(path)
            print *, filename
            ! Global number of cells
            ncells = psb_cd_get_global_cols(msh%desc_c)

            CALL psb_geall(i_loc,msh%desc_c,info)
            CALL psb_check_error(info,'write_mesh','psb_geall',icontxt)

            ALLOCATE(igroup(ncells),iproc(ncells),stat=info)
            IF(info /= 0) THEN
                WRITE(*,100)
                CALL abort_psblas
            END IF

            ! Gathers MESH components on P0
            CALL l2g_vertex(msh%verts,verts,msh%desc_v)
            CALL l2g_face(msh%faces,faces,msh%desc_f,msh%desc_c)
            CALL l2g_cell(msh%cells,cells,msh%desc_c)
            CALL l2g_conn(msh%v2f,v2f,msh%desc_v,msh%desc_f)
            CALL l2g_conn(msh%v2c,v2c,msh%desc_v,msh%desc_c)
            CALL l2g_conn(msh%f2c,f2c,msh%desc_f,msh%desc_c)

            ! Gathers processor IDs
            i_loc(:) = mypnum
            CALL psb_gather(iproc,i_loc,msh%desc_c,info,root=0)
            CALL psb_check_error(info,'write_mesh','psb_gather',icontxt)

            ! Gathers group IDs
            i_loc(:) = 0
            ngroups = msh%c2g%nel_()

            DO ig = 1, ngroups
                CALL msh%c2g%get_ith_conn(ic2g,ig)
                ngc = SIZE(ic2g) ! number of group cells
                DO i = 1, ngc
                    ic = ic2g(i)
                    i_loc(ic) = ig
                END DO
            END DO

            CALL psb_gather(igroup,i_loc,msh%desc_c,info,root=0)
            CALL psb_check_error(info,'write_mesh','psb_gather',icontxt)

            ! Open exodusII file for writing
            cpu_ws = 0
            io_ws = 8
            call exopts (EXVRBS, IERR) !Verbose error reporting of EXODUS routines
            exodus_file_id = excre("output.exo", EXCLOB, cpu_ws, io_ws, ierr)

            num_dims = msh%ncd
            num_elem_blks = ngroups
            num_node_sets = 0
            num_side_sets = msh%nbc
            num_nodes = size(verts)
            num_elems = size(cells)

            ! Write the database parameters to the exodus file
            CALL expini(exodus_file_id, "MORFEUS Output", num_dims, num_nodes, &
              num_elems, num_elem_blks, num_node_sets, num_side_sets, ierr)

            ! Write the mesh coordinates
            ALLOCATE(xx(num_nodes), yy(num_nodes), zz(num_nodes))

            xx(:) = verts(:)%x_()
            yy(:) = verts(:)%y_()
            zz(:) = verts(:)%z_()

            CALL expcor(exodus_file_id, xx, yy, zz, ierr)

            ! Write the mesh coordinate names
            coord_names(1) = "X"
            coord_names(2) = "Y"
            coord_names(3) = "Z"
            CALL expcon(exodus_file_id, coord_names, ierr)

            ! Write the element block parameters
            DO ig = 1, num_elem_blks !Loop over groups
              CALL msh%c2g%get_ith_conn(ic2g,ig)
              nc2g = SIZE(ic2g) !Number of cells in group
              CALL msh%v2c%get_ith_conn(iv2c,1) ! First cell
              nv2c = SIZE(iv2c) ! Number of nodes per cell
              ctype = "HEX8"
              CALL expelb(exodus_file_id, INT(ig,INT64), ctype, INT(nc2g,INT64), INT(nv2c,INT64), 1_INT64, ierr)

              nconn = nv2c*nc2g !Connectivity vector size

              print *, ig, nconn
              IF (ALLOCATED(conn)) DEALLOCATE(conn)
              ALLOCATE(conn(nconn))
              i = 0
              DO ic = 1, nc2g
                CALL msh%v2c%get_ith_conn(iv2c,ic)
                !print *, ic, i, iv2c
                conn(i+1:i+nv2c) = iv2c
                i = i + nv2c
              END DO
              call expelc(exodus_file_id, INT(ig,INT64), conn, ierr)
            END DO

            ! Loop over scalar fields and gather them on the root processor
            DO i = 1, size(sfield)
              ! Is the field cell-centered?
              IF(sfield(i)%on_faces_()) THEN
                  WRITE(*,100)
                  CYCLE
              END IF
              CALL sfield(i)%get_x(scalar_local)

              ALLOCATE(scalar_global(ncells),stat=info)
              IF(info /= 0) THEN
                  WRITE(*,200)
                  CALL abort_psblas
              END IF

              ! Gathers cell-centered values
              CALL psb_gather(scalar_global,scalar_local,msh%desc_c,info,root=0)

              IF(mypnum == 0) THEN
                DO ig = 1, num_elem_blks
                  CALL msh%c2g%get_ith_conn(ic2g,ig)
                  nc2g = SIZE(ic2g) !Number of cells in group
                  IF(ALLOCATED(f_s)) DEALLOCATE(f_s)
                  ALLOCATE(f_s(nc2g))
                  f_s(1:nc2g) = scalar_global(ic2g(1:nc2g))
                  CALL expev(exodus_file_id, 0.0, i, INT(ig, INT64), INT(nc2g, INT64), f_s, ierr)
                END DO
              END IF
            END DO

            ! Loop over vector fields and gather them on the root processor
            DO i = 1, size(vfield)
              ! Is the field cell-centered?
              IF(vfield(i)%on_faces_()) THEN
                  WRITE(*,100)
                  CYCLE
              END IF
              CALL vfield(i)%get_x(vector_local)

              ALLOCATE(vector_global(ncells),stat=info)
              IF(info /= 0) THEN
                  WRITE(*,200)
                  CALL abort_psblas
              END IF

              ! Gathers cell-centered values
              CALL l2g_vector(vector_global,vector_local,msh%desc_c)

              IF(mypnum == 0) THEN
                DO ig = 1, num_elem_blks
                  CALL msh%c2g%get_ith_conn(ic2g,ig)
                  nc2g = SIZE(ic2g) !Number of cells in group
                  IF(ALLOCATED(f_x)) DEALLOCATE(f_x, f_y, f_z)
                  ALLOCATE(f_x(nc2g), f_y(nc2g), f_z(nc2g))
                  f_x(1:nc2g) = vector_global(ic2g(1:nc2g))%x_()
                  f_y(1:nc2g) = vector_global(ic2g(1:nc2g))%y_()
                  f_z(1:nc2g) = vector_global(ic2g(1:nc2g))%z_()
                  CALL expev(exodus_file_id, 0.0, i, INT(ig, INT64), INT(nc2g, INT64), f_s, ierr)
                END DO
              END IF
            END DO

            ! DO i = 1, size(sfield) + 3*size(vfield)
            !   IF (i <= size(sfield)) THEN
            !     !write sfield to exodus file
            !   ELSE

            !   END IF
            ! END DO
            ! Close the exodus file
            CALL exclos(exodus_file_id, ierr)
            STOP

100         FORMAT(' WARNING! Face-centered field in WRITE_EXODUS. Field will not be output.')
200         FORMAT(' ERROR! Memory allocation failure in WRITE_EXODUS')

        END PROCEDURE write_exo_morfeus

END SUBMODULE class_exodus_procedures

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


text:
program concatenate
implicit none
integer, parameter :: n1 = 2, n2 = 3
integer :: imat(n1,n2),ivec(n1),i
integer, allocatable :: jmat(:,:)
imat(1,:) = [10,11,12]
imat(2,:) = [20,21,22]
ivec = 10*imat(:,n2) ! [120,220]
jmat = reshape([imat,ivec,4,9],[n1,n2+2])
do i=1,size(jmat,1)
   print "(*(i5))",jmat(i,:)
end do
end program concatenate
! output:
!   10   11   12  120    4
!   20   21   22  220    9

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


! RUN: %B/test/Semantics/test_errors.sh %s %flang %t
! Check that a basic arithmetic if compiles.

if ( A ) 100, 200, 300
100 CONTINUE
200 CONTINUE
300 CONTINUE
END


text:
Cathy Kennedy ran for Yolo County Supervisor in the 4^th^ District in the 2008 City Council Election Spring 2008 Election.  Her campaign website, now defunct, was {{{http://www.Kennedy4supervisor.com}}}.

 Campaign Literature

20080602 01:50:22 nbsp Kennedy is the Republican candidate for County Supervisor. She registered over to decline to state less than three months ago. While I agree she has a compelling story, she is not the right candidate for supervisor, especially given the credentials of John Ferrera Ferrera and Jim Provenza Provenza. Users/MaxMikalonis

20081115 15:32:11 nbsp Kennedy is a combo; certainly fiscally conservative; but certainly a social liberal. This lady grew up in social poverty, freely admits that she was a welfare child; put herself through school; takes in children; seems to be rather savvy on business and political issues. With the layout of the BOS, she would be a very compelling addition. Yolo has 4 men and 1 woman; 2 lawyers, 2 professional politicians and 1 farmer. Kennedy does not need Ferrera and Provenza credentials  they are weak. Provenza did a bang job on the school board, and Ferrera has nothing but Califonia Budget experience  like thats a bragging point. Kennedy has fire and doesn’t seem to pander to anyone unlike what weve all seen from Ferrera and Provenza. True the race is done  she entered too late. Kennedy should run for something else but, she might be more than little Yolo can comprehend    Users/YoloPeep


text:
      module parkind

      implicit none
      save

!------------------------------------------------------------------
! rrtmg kinds
! Define integer and real kinds for various types.
!
! Initial version: MJIacono, AER, jun2006
! Revised: MJIacono, AER, aug2008
!------------------------------------------------------------------

!
!     integer kinds
!     -------------
!
      integer, parameter :: kind_ib = selected_int_kind(13)  ! 8 byte integer
      integer, parameter :: kind_im = selected_int_kind(6)   ! 4 byte integer
      integer, parameter :: kind_in = kind(1)                ! native integer

!
!     real kinds
!     ----------
!
      integer, parameter :: kind_rb = selected_real_kind(12) ! 8 byte real
      integer, parameter :: kind_rm = selected_real_kind(6)  ! 4 byte real
      integer, parameter :: kind_rn = kind(1.0)              ! native real

      end module parkind

max_stars_repo_name: nfeldl/climlab

max_issues_repo_name: nfeldl/climlab

max_forks_repo_name: nfeldl/climlab


text:
^^wiki:westsac:Users/BlueDot also^^

20090712 06:06:06 nbsp oH Hai n00b!, joo want just ice, oda wiki:wikipedia:William Leete |_337ness?  TY4 50Ul ;|! Users/thedevil


text:
*DECK R9GMIC
      FUNCTION R9GMIC (A, X, ALX)
C***BEGIN PROLOGUE  R9GMIC
C***SUBSIDIARY
C***PURPOSE  Compute the complementary incomplete Gamma function for A
C            near a negative integer and for small X.
C***LIBRARY   SLATEC (FNLIB)
C***CATEGORY  C7E
C***TYPE      SINGLE PRECISION (R9GMIC-S, D9GMIC-D)
C***KEYWORDS  COMPLEMENTARY INCOMPLETE GAMMA FUNCTION, FNLIB, SMALL X,
C             SPECIAL FUNCTIONS
C***AUTHOR  Fullerton, W., (LANL)
C***DESCRIPTION
C
C Compute the complementary incomplete gamma function for A near
C a negative integer and for small X.
C
C***REFERENCES  (NONE)
C***ROUTINES CALLED  ALNGAM, R1MACH, XERMSG
C***REVISION HISTORY  (YYMMDD)
C   770701  DATE WRITTEN
C   890531  Changed all specific intrinsics to generic.  (WRB)
C   890531  REVISION DATE from Version 3.2
C   891214  Prologue converted to Version 4.0 format.  (BAB)
C   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
C   900720  Routine changed from user-callable to subsidiary.  (WRB)
C***END PROLOGUE  R9GMIC
      SAVE EULER, EPS, BOT
      DATA EULER / .5772156649 015329 E0 /
      DATA EPS, BOT / 2*0.0 /
C***FIRST EXECUTABLE STATEMENT  R9GMIC
      IF (EPS.EQ.0.0) EPS = 0.5*R1MACH(3)
      IF (BOT.EQ.0.0) BOT = LOG(R1MACH(1))
C
      IF (A .GT. 0.0) CALL XERMSG ('SLATEC', 'R9GMIC',
     +   'A MUST BE NEAR A NEGATIVE INTEGER', 2, 2)
      IF (X .LE. 0.0) CALL XERMSG ('SLATEC', 'R9GMIC',
     +   'X MUST BE GT ZERO', 3, 2)
C
      MA = A - 0.5
      FM = -MA
      M = -MA
C
      TE = 1.0
      T = 1.0
      S = T
      DO 20 K=1,200
        FKP1 = K + 1
        TE = -X*TE/(FM+FKP1)
        T = TE/FKP1
        S = S + T
        IF (ABS(T).LT.EPS*S) GO TO 30
 20   CONTINUE
      CALL XERMSG ('SLATEC', 'R9GMIC',
     +   'NO CONVERGENCE IN 200 TERMS OF CONTINUED FRACTION', 4, 2)
C
 30   R9GMIC = -ALX - EULER + X*S/(FM+1.0)
      IF (M.EQ.0) RETURN
C
      IF (M.EQ.1) R9GMIC = -R9GMIC - 1.0 + 1.0/X
      IF (M.EQ.1) RETURN
C
      TE = FM
      T = 1.0
      S = T
      MM1 = M - 1
      DO 40 K=1,MM1
        FK = K
        TE = -X*TE/FK
        T = TE/(FM-FK)
        S = S + T
        IF (ABS(T).LT.EPS*ABS(S)) GO TO 50
 40   CONTINUE
C
 50   DO 60 K=1,M
        R9GMIC = R9GMIC + 1.0/K
 60   CONTINUE
C
      SGNG = 1.0
      IF (MOD(M,2).EQ.1) SGNG = -1.0
      ALNG = LOG(R9GMIC) - ALNGAM(FM+1.0)
C
      R9GMIC = 0.0
      IF (ALNG.GT.BOT) R9GMIC = SGNG*EXP(ALNG)
      IF (S.NE.0.0) R9GMIC = R9GMIC + SIGN (EXP(-FM*ALX+LOG(ABS(S)/FM))
     1  , S)
C
      IF (R9GMIC .EQ. 0.0 .AND. S .EQ. 0.0) CALL XERMSG ('SLATEC',
     +   'R9GMIC', 'RESULT UNDERFLOWS', 1, 1)
      RETURN
C
      END


text:
great_ne(v1,b1).
great_ne(t1,bb1).
great_ne(p1,dd1).
great_ne(jj1,h1).
great_ne(aa1,ff1).
great_ne(q1,x1).
great_ne(u1,i1).
great_ne(q1,z1).
great_ne(kk1,bb1).
great_ne(cc1,u1).
great_ne(q1,y1).
great_ne(ee1,kk1).
great_ne(kk1,i1).
great_ne(x1,v1).
great_ne(dd1,f1).
great_ne(z1,ff1).
great_ne(u1,v1).
great_ne(jj1,a1).
great_ne(z1,jj1).
great_ne(k1,v1).
great_ne(hh1,bb1).
great_ne(w1,f1).
great_ne(b1,i1).
great_ne(t1,v1).
great_ne(p1,ff1).
great_ne(u1,x1).
great_ne(bb1,i1).
great_ne(z1,a1).
great_ne(q1,jj1).
great_ne(aa1,c1).
great_ne(m1,w1).
great_ne(y1,b1).
great_ne(jj1,bb1).
great_ne(jj1,y1).

max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp