! CHECK DIMENSION
function check_dimension_kb_contour_gf(G,params) result(bool)
  type(kb_contour_gf)     :: G
  type(kb_contour_params) :: params
  logical                 :: bool
  integer                 :: N,L
  bool=.false.
  N=params%Ntime              !<== check size at maximum time
  L=params%Ntau
  if( (size(G%less)/=N**2) .OR. (size(G%ret)/=N**2) )&
       stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions less/ret"
  if( size(G%lmix)/=N*(L+1) )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions lmix"
  if( size(G%mats)/=(L+1) )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions mats"
  bool=.true.
end function check_dimension_kb_contour_gf
!
function check_dimension_kb_contour_sigma(S,params) result(bool)
  type(kb_contour_sigma)  :: S
  type(kb_contour_params) :: params
  logical                 :: bool
  integer                 :: N,L
  bool=.false.
  N=params%Ntime              !<== check size at maximum time
  L=params%Ntau
  if( (size(S%reg%less)/=N**2) .OR. (size(S%reg%ret)/=N**2) )&
       stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions less/ret"
  if( size(S%reg%lmix)/=N*(L+1) )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions lmix"
  if( size(S%reg%mats)/=(L+1) )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions mats"
  if( size(S%hf)/=N )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions HF"
  bool=.true.
end function check_dimension_kb_contour_sigma
!
function check_dimension_kb_contour_gf_(G,N,L) result(bool)
  type(kb_contour_gf)     :: G
  logical                 :: bool
  integer                 :: N,L
  bool=.false.
  if( (size(G%less)/=N**2) .OR. (size(G%ret)/=N**2) )&
       stop "ERROR contour_gf/check_dimension_kb_contour_gf: wrong dimensions less/ret"
  if( size(G%lmix)/=N*(L+1) )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions lmix"
  if( size(G%mats)/=(L+1) )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions mats"
  bool=.true.
end function check_dimension_kb_contour_gf_
!
function check_dimension_kb_contour_sigma_(S,N,L) result(bool)
  type(kb_contour_sigma)  :: S
  logical                 :: bool
  integer                 :: N,L
  bool=.false.
  if( (size(S%reg%less)/=N**2) .OR. (size(S%reg%ret)/=N**2) )&
       stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions less/ret"
  if( size(S%reg%lmix)/=N*(L+1) )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions lmix"
  if( size(S%reg%mats)/=(L+1) )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions mats"
  if( size(S%hf)/=N )stop "contour_gf/check_dimension_kb_contour_gf: wrong dimensions HF"
  bool=.true.
end function check_dimension_kb_contour_sigma_
!
function check_dimension_kb_contour_dgf(dG,params) result(bool)
  type(kb_contour_dgf)    :: dG
  type(kb_contour_params) :: params
  logical                 :: bool
  integer                 :: N,L
  bool=.false.
  N=params%Ntime              !<== check size at maximum time
  L=params%Ntau
  if( (size(dG%less)/=N) .OR. (size(dG%ret)/=N) )&
       stop "ERROR contour_gf/check_dimension_kb_contour_dgf: wrong dimensions less/ret"
  if( size(dG%lmix)/=(L+1) )&
       stop "ERROR contour_gf/check_dimension_kb_contour_dgf: wrong dimensions lmix"
  bool=.true.
end function check_dimension_kb_contour_dgf
!
function check_dimension_kb_contour_dgf_(dG,N,L) result(bool)
  type(kb_contour_dgf)    :: dG
  logical                 :: bool
  integer                 :: N,L
  bool=.false.
  if( (size(dG%less)/=N) .OR. (size(dG%ret)/=N) )&
       stop "ERROR contour_gf/check_dimension_kb_contour_dgf: wrong dimensions less/ret"
  if( size(dG%lmix)/=(L+1) )&
       stop "ERROR contour_gf/check_dimension_kb_contour_dgf: wrong dimensions lmix"
  bool=.true.
end function check_dimension_kb_contour_dgf_

! SAVE
subroutine save_kb_contour_gf(G,file)
  type(kb_contour_gf) :: G
  character(len=*)    :: file
  integer             :: unit,N,L,Lf
  if(.not.G%status)stop "contour_gf/save_kb_contour_gf: G not allocated"
  unit = free_unit()
  open(unit,file=reg(file)//"_dimension.data.neqipt")
  N=size(G%less,1)
  L=size(G%lmix,2)-1
  Lf=size(G%iw)
  write(unit,"(A1,3A4,3X)")"#","N","L","Lf"
  write(unit,"(1X,3I4,3X)")N,L,Lf
  close(unit)
  call save_array(reg(file)//"_less.data.neqipt",G%less(:,:))
  call save_array(reg(file)//"_ret.data.neqipt", G%ret(:,:))
  call save_array(reg(file)//"_lmix.data.neqipt",G%lmix(:,0:))
  call save_array(reg(file)//"_mats.data.neqipt",G%mats(0:))
  call save_array(reg(file)//"_iw.data.neqipt",G%iw(:))
end subroutine save_kb_contour_gf
!
subroutine save_kb_contour_sigma(S,file)
  type(kb_contour_sigma) :: S
  character(len=*)       :: file
  integer                :: unit,i
  if(.not.S%status)stop "contour_gf/save_kb_contour_gf: S not allocated"
  unit = free_unit()
  open(unit,file=reg(file)//"_hfb.data.neqipt")
  do i=1,size(S%hf(:))
     write(unit,*)S%hf(i)
  enddo
  close(unit)
  call save_kb_contour_gf(S%reg,file)
end subroutine save_kb_contour_sigma
!
subroutine save_kb_contour_dgf(dG,file)
  type(kb_contour_dgf) :: dG
  character(len=*)     :: file
  integer              :: unit,N,L
  if(.not.dG%status)stop "contour_gf/save_kb_contour_dgf: dG not allocated"
  unit = free_unit()
  N=size(dG%less)
  L=size(dG%lmix)-1
  open(unit,file=reg(file)//"_dimension.data.neqipt")
  write(unit,"(A1,2A4,2X)")"#","N","L"
  write(unit,"(1X,2I4,2X)")N,L
  close(unit)
  call save_array(reg(file)//"_less.data.neqipt",dG%less(:))
  call save_array(reg(file)//"_ret.data.neqipt", dG%ret(:))
  call save_array(reg(file)//"_lmix.data.neqipt",dG%lmix(0:))
end subroutine save_kb_contour_dgf

! READ
subroutine read_kb_contour_gf(G,file)
  type(kb_contour_gf)  :: G
  character(len=*)     :: file
  logical              :: check
  check = inquire_kb_contour_gf(file)
  if(.not.G%status.OR..not.check)stop "contour_gf/read_kb_contour_gf: G not allocated"
  call read_array(trim(file)//"_less.data.neqipt",G%less(:,:))
  call read_array(trim(file)//"_ret.data.neqipt",G%ret(:,:))
  call read_array(trim(file)//"_lmix.data.neqipt",G%lmix(:,0:))
  call read_array(trim(file)//"_mats.data.neqipt",G%mats(0:))
  call read_array(trim(file)//"_iw.data.neqipt",G%iw(:))
end subroutine read_kb_contour_gf
!
subroutine read_kb_contour_sigma(S,file)
  type(kb_contour_sigma) :: S
  character(len=*)       :: file
  logical                :: check
  integer                :: unit,i,Len
  check = inquire_kb_contour_sigma(file)
  if(.not.S%status.OR..not.check)stop "contour_gf/read_kb_contour_sigma: S not allocated"
  Len  = file_length(reg(file)//"_hfb.data.neqipt")
  if( Len > size(S%hf) ) stop "contour_gf/read_kb_contour_sigma: length(file) > size(S%hf)"
  unit = free_unit()
  open(unit,file=reg(file)//"_hfb.data.neqipt")
  do i=1,Len
     read(unit,*)S%hf(i)
  enddo
  close(unit)
  call read_kb_contour_gf(S%reg,file)
end subroutine read_kb_contour_sigma
!
subroutine read_kb_contour_dgf(dG,file)
  type(kb_contour_dgf) :: dG
  character(len=*)     :: file
  logical              :: check
  check = inquire_kb_contour_dgf(file)
  if(.not.dG%status.OR..not.check)stop "contour_gf/read_kb_contour_dgf: dG not allocated"
  call read_array(trim(file)//"_less.data.neqipt",dG%less(:))
  call read_array(trim(file)//"_ret.data.neqipt",dG%ret(:))
  call read_array(trim(file)//"_lmix.data.neqipt",dG%lmix(0:))
end subroutine read_kb_contour_dgf
!
!
!aux: inquire:
function inquire_kb_contour_gf(file) result(check)
  integer          :: i
  logical          :: check,bool(5)
  character(len=*) :: file
  character(len=16),dimension(5)  :: ctype=([ character(len=5) :: 'less','ret','lmix','mats','iw'])
  check=.true.
  do i=1,5
     inquire(file=reg(file)//"_"//reg(ctype(i))//".data.neqipt",exist=bool(i))
     if(.not.bool(i))inquire(file=reg(file)//"_"//reg(ctype(i))//".data.neqipt.gz",exist=bool(i))
     check=check.AND.bool(i)
  enddo
end function inquire_kb_contour_gf
!
function inquire_kb_contour_sigma(file) result(check)
  integer          :: i
  logical          :: check,bool(6)
  character(len=*) :: file
  character(len=16),dimension(5)  :: ctype=([ character(len=5) :: 'less','ret','lmix','mats','iw'])
  check=.true.
  do i=1,5
     inquire(file=reg(file)//"_"//reg(ctype(i))//".data.neqipt",exist=bool(i))
     if(.not.bool(i))inquire(file=reg(file)//"_"//reg(ctype(i))//".data.neqipt.gz",exist=bool(i))
     check=check.AND.bool(i)
  enddo
  inquire(file=reg(file)//"_hfb.data.neqipt",exist=bool(6))
  check=check.AND.bool(6)
end function inquire_kb_contour_sigma
!
function inquire_kb_contour_dgf(file) result(check)
  integer          :: i
  logical          :: check,bool(3)
  character(len=*) :: file
  character(len=16),dimension(3)  :: ctype=([character(len=5) :: 'less','ret','lmix'])
  check=.true.
  do i=1,3
     inquire(file=reg(file)//"_"//reg(ctype(i))//".data.neqipt",exist=bool(i))
     if(.not.bool(i))inquire(file=reg(file)//"_"//reg(ctype(i))//".data.neqipt.gz",exist=bool(i))
     check=check.AND.bool(i)
  enddo
end function inquire_kb_contour_dgf

! PLOT
subroutine plot_kb_contour_gf(file,G,params)
  character(len=*)        :: file
  type(kb_contour_gf)     :: G
  type(kb_contour_params) :: params
  integer                 :: Nt,i
  if(.not.G%status)stop "contour_gf/plot_kb_contour_gf: G is not allocated"
  Nt=params%Ntime
  call splot3d(reg(file)//"_less_t_t.data.neqipt",params%t(:Nt),params%t(:Nt),G%less(:Nt,:Nt))
  call splot3d(reg(file)//"_ret_t_t.data.neqipt",params%t(:Nt),params%t(:Nt),G%ret(:Nt,:Nt))
  call splot3d(reg(file)//"_lmix_t_tau.data.neqipt",params%t(:Nt),params%tau(0:),G%lmix(:Nt,0:))
  call splot(reg(file)//"_mats_tau.data.neqipt",params%tau(0:),G%mats(0:))
  call splot(reg(file)//"_mats_iw.data.neqipt",params%wm(:),G%iw(:))
end subroutine plot_kb_contour_gf
!
subroutine plot_kb_contour_sigma(file,S,params)
  character(len=*)        :: file
  type(kb_contour_sigma)  :: S
  type(kb_contour_params) :: params
  integer                 :: Nt,i
  if(.not.S%status)stop "contour_gf/plot_kb_contour_sigma: S is not allocated"
  Nt=params%Ntime
  call splot3d(reg(file)//"_less_t_t.data.neqipt",params%t(:Nt),params%t(:Nt),S%reg%less(:Nt,:Nt))
  call splot3d(reg(file)//"_ret_t_t.data.neqipt",params%t(:Nt),params%t(:Nt),S%reg%ret(:Nt,:Nt))
  call splot3d(reg(file)//"_lmix_t_tau.data.neqipt",params%t(:Nt),params%tau(0:),S%reg%lmix(:Nt,0:))
  call splot(reg(file)//"_mats_tau.data.neqipt",params%tau(0:),S%reg%mats(0:))
  call splot(reg(file)//"_mats_iw.data.neqipt",params%wm(:),S%reg%iw(:)+S%hf(1))
  call splot(reg(file)//"_hfb_t.data.neqipt",params%t(:Nt),S%hf(:Nt))
end subroutine plot_kb_contour_sigma
!
subroutine plot_kb_contour_dgf(file,dG,params)
  character(len=*)        :: file
  type(kb_contour_dgf)    :: dG
  type(kb_contour_params) :: params
  integer                 :: Nt
  if(.not.dG%status)stop "contour_gf/plot_kb_contour_gf: G is not allocated"
  Nt=params%Ntime
  call splot(reg(file)//"_less_t.data.neqipt",params%t(:Nt),dG%less(:Nt))
  call splot(reg(file)//"_ret_t.data.neqipt",params%t(:Nt),dG%ret(:Nt))
  call splot(reg(file)//"_lmix_tau.data.neqipt",params%tau(0:),dG%lmix(0:))
end subroutine plot_kb_contour_dgf

! EXTRAPOLATION
!extrapolate a function from a given time to the next one:
subroutine extrapolate_kb_contour_gf(g,params)
  type(kb_contour_gf)     :: g
  type(kb_contour_params) :: params
  integer                 :: i,j,k,N,L
  if(.not.g%status)     stop "extrapolate_kb_contour_gf: g is not allocated"
  if(.not.params%status)stop "extrapolate_kb_contour_gf: params is not allocated"
  N = params%Nt
  L = params%Ntau
  select case(N)
  case(1)
     return
  case(2)
     !GUESS G AT THE NEXT STEP, GIVEN THE INITIAL CONDITIONS
     do j=1,N
        g%ret(N,j) =g%ret(1,1)
        g%less(N,j)=g%less(1,1)
     end do
     do i=1,N-1
        g%less(i,N)=g%less(1,1)
     end do
     do j=0,L
        g%lmix(N,j)=g%lmix(1,j)
     end do
  case default
     !EXTEND G FROM THE [N-1,N-1] TO THE [N,N] SQUARE TO START DMFT
     !USING QUADRATIC EXTRAPOLATION
     do k=1,N-1
        g%less(N,k)=2.d0*g%less(N-1,k)-g%less(N-2,k)
        g%less(k,N)=2.d0*g%less(k,N-1)-g%less(k,N-2)
     end do
     g%less(N,N)=2.d0*g%less(N-1,N-1)-g%less(N-2,N-2)
     !
     do k=0,L
        g%lmix(N,k)=2.d0*g%lmix(N-1,k)-g%lmix(N-2,k)
     end do
     !
     g%ret(N,N)=-xi
     do k=1,N-2
        g%ret(N,k)=2.d0*g%ret(N-1,k)-g%ret(N-2,k)
     end do
     g%ret(N,N-1)=0.5d0*(g%ret(N,N)+g%ret(N,N-2))
  end select
end subroutine extrapolate_kb_contour_gf

subroutine extrapolate_kb_contour_sigma(s,params)
  type(kb_contour_sigma)  :: s
  type(kb_contour_params) :: params
  integer                 :: i,j,k,N,L
  if(.not.s%status)     stop "extrapolate_kb_contour_gf: g is not allocated"
  if(.not.params%status)stop "extrapolate_kb_contour_gf: params is not allocated"
  N = params%Nt
  select case(N)
  case(1)
     return
  case(2)
     !GUESS Shf AT THE NEXT STEP, GIVEN THE INITIAL CONDITIONS
     s%hf(N)=s%hf(1)
  case default
     !EXTEND Shf FROM THE [N-1,N-1] TO THE [N,N] SQUARE TO START DMFT
     !USING QUADRATIC EXTRAPOLATION
     S%hf(N)=2.d0*S%hf(N-1)-S%hf(N-2)
  end select
  call extrapolate_kb_contour_gf(s%reg,params)
end subroutine extrapolate_kb_contour_sigma


! This file was generated by gpufort

!
! Hints:
! Device variables in scope:
!       character                                   :: uplo, jobz

!       character                                   :: uplo, jobz

!       integer                                     :: n, nz, lda, lwork, istat, info

!       integer                                     :: n, nz, lda, lwork, istat, info

!       integer                                     :: n, nz, lda, lwork, istat, info

!       integer                                     :: n, nz, lda, lwork, istat, info

!       integer                                     :: n, nz, lda, lwork, istat, info

!       integer                                     :: n, nz, lda, lwork, istat, info

!       integer                                     :: lwork_h, liwork_h, ldz_h

!       integer                                     :: lwork_h, liwork_h, ldz_h

!       integer                                     :: lwork_h, liwork_h, ldz_h

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: i, j, k, nb1, nb2, ib, mi, ldt, ldz, il, iu

!       integer                                     :: inde, indtau, indwrk, llwork, llwork_h, indwk2, indwk3, llwrk2

!       integer                                     :: inde, indtau, indwrk, llwork, llwork_h, indwk2, indwk3, llwrk2

!       integer                                     :: inde, indtau, indwrk, llwork, llwork_h, indwk2, indwk3, llwrk2

!       integer                                     :: inde, indtau, indwrk, llwork, llwork_h, indwk2, indwk3, llwrk2

!       integer                                     :: inde, indtau, indwrk, llwork, llwork_h, indwk2, indwk3, llwrk2

!       integer                                     :: inde, indtau, indwrk, llwork, llwork_h, indwk2, indwk3, llwrk2

!       integer                                     :: inde, indtau, indwrk, llwork, llwork_h, indwk2, indwk3, llwrk2

!       integer                                     :: inde, indtau, indwrk, llwork, llwork_h, indwk2, indwk3, llwrk2

!       real(8), parameter                          :: one = 1.0_8

!       type(dim3) :: blocks, threads

!       type(dim3) :: blocks, threads

function test_launch_krnl_e26a05_0_auto()
  ! errorCode > 0 implies that the test has failed
  use iso_c_binding
  use hip
  use dsyevd_gpu_kernels
  implicit none
  integer :: errorCode = 1
  ! TODO fix parameters
  ! - Add missing arguments
  ! - Determine size of arrays (typically indicated by 'type(c_ptr)' type)
  ! - Add target where we need a pointer
  integer(c_int),intent(IN) :: sharedMem
  type(c_ptr),value,intent(IN) :: stream
  integer,value :: n
  TODO declaration not found :: z
  TODO declaration not found :: a
  ! TODO Create initial data on host
  ! TODO Copy data to device ! (dest,src,size,direction)
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  ! ... might be more (or less) than two memcopies
  ! TODO run the test
  CALL launch_krnl_e26a05_0_auto(0,c_null_ptr,sharedMem,stream,n,z,a) ! Modify sharedMem if other than default 0
  CALL hipCheck(hipDeviceSynchronize())
  CALL launch_krnl_e26a05_0_cpu(0,c_null_ptr,sharedMem,stream,n,z,a)

  ! TODO Copy results back to host
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  ! ... might be more (or less) than two memcopies
  ! TODO Compare results
  ! TODO Update error code if the results do not match
  return errorCode
end function

!
! Hints:
! Device variables in scope:
!       integer                               :: n, k, ldv, ldt, ldw

!       integer                               :: n, k, ldv, ldt, ldw

!       integer                               :: n, k, ldv, ldt, ldw

!       integer                               :: n, k, ldv, ldt, ldw

!       integer                               :: n, k, ldv, ldt, ldw

!       real(8), dimension(ldv, k), device    :: v

!       real(8), dimension(k), device         :: tau

!       real(8), dimension(ldt, k), device    :: t

!       real(8), dimension(ldw, k), device    :: w

!       integer                               :: i, j, istat

!       integer                               :: i, j, istat

!       integer                               :: i, j, istat

!       type(dim3)                            :: threads

function test_launch_krnl_b1f342_1_auto()
  ! errorCode > 0 implies that the test has failed
  use iso_c_binding
  use hip
  use dsyevd_gpu_kernels
  implicit none
  integer :: errorCode = 1
  ! TODO fix parameters
  ! - Add missing arguments
  ! - Determine size of arrays (typically indicated by 'type(c_ptr)' type)
  ! - Add target where we need a pointer
  integer(c_int),intent(IN) :: sharedMem
  type(c_ptr),value,intent(IN) :: stream
  type(c_ptr),value :: w
  integer,value :: n
  integer,value :: k
  type(c_ptr),value :: v
  ! TODO Create initial data on host
  ! TODO Copy data to device ! (dest,src,size,direction)
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  ! ... might be more (or less) than two memcopies
  ! TODO run the test
  CALL launch_krnl_b1f342_1_auto(0,c_null_ptr,sharedMem,stream,w,n,k,v) ! Modify sharedMem if other than default 0
  CALL hipCheck(hipDeviceSynchronize())
  CALL launch_krnl_b1f342_1_cpu(0,c_null_ptr,sharedMem,stream,w,n,k,v)

  ! TODO Copy results back to host
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  ! ... might be more (or less) than two memcopies
  ! TODO Compare results
  ! TODO Update error code if the results do not match
  return errorCode
end function

!
! Hints:
! Device variables in scope:
!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: m, n, k, ldv, ldt, ldc, ldw, ldwork, istat

!       integer                               :: i, j

!       integer                               :: i, j

!       real(8), dimension(ldv, k), device    :: v

!       real(8), dimension(ldt, k), device    :: t

!       real(8), dimension(ldw, k), device    :: w

!       real(8), dimension(ldc, n), device    :: c

!       real(8), dimension(ldwork, k), device :: work

function test_launch_krnl_b95769_2_auto()
  ! errorCode > 0 implies that the test has failed
  use iso_c_binding
  use hip
  use dsyevd_gpu_kernels
  implicit none
  integer :: errorCode = 1
  ! TODO fix parameters
  ! - Add missing arguments
  ! - Determine size of arrays (typically indicated by 'type(c_ptr)' type)
  ! - Add target where we need a pointer
  integer(c_int),intent(IN) :: sharedMem
  type(c_ptr),value,intent(IN) :: stream
  integer,value :: m
  type(c_ptr),value :: w
  integer,value :: k
  type(c_ptr),value :: v
  ! TODO Create initial data on host
  ! TODO Copy data to device ! (dest,src,size,direction)
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  ! ... might be more (or less) than two memcopies
  ! TODO run the test
  CALL launch_krnl_b95769_2_auto(0,c_null_ptr,sharedMem,stream,m,w,k,v) ! Modify sharedMem if other than default 0
  CALL hipCheck(hipDeviceSynchronize())
  CALL launch_krnl_b95769_2_cpu(0,c_null_ptr,sharedMem,stream,m,w,k,v)

  ! TODO Copy results back to host
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  ! ... might be more (or less) than two memcopies
  ! TODO Compare results
  ! TODO Update error code if the results do not match
  return errorCode
end function

! Fortran implementation:
!     attributes(global) subroutine finish_T_block_kernel(N, T, ldt, tau)
!        implicit none
!        integer, value                     :: N, ldt
!        real(8), dimension(ldt, K), device :: T
!        real(8), dimension(K), device      :: tau
!        ! T_s contains only lower triangular elements of T in linear array, by row
!        real(8), dimension(2080), shared   :: T_s
!        ! (i,j) --> ((i-1)*i/2 + j)
! #define IJ2TRI(i,j) (ISHFT((i-1)*i,-1) + j)
!
!        integer     :: tid, tx, ty, i, j, k, diag
!        complex(8)  :: cv
!
!        tx = threadIdx%x
!        ty = threadIdx%y
!        tid = (threadIdx%y - 1)*blockDim%x + tx ! Linear thread id
!
!        ! Load T into shared memory
!        if (tx <= N) then
!           do j = ty, N, blockDim%y
!              cv = tau(j)
!              if (tx > j) then
!                 T_s(IJ2TRI(tx, j)) = -cv*T(tx, j)
!              else if (tx == j) then
!                 T_s(IJ2TRI(tx, j)) = cv
!              endif
!           end do
!        end if
!
!        call syncthreads()
!
!        ! Perform column by column update by first thread column
!        do i = N - 1, 1, -1
!           if (ty == 1) then
!              if (tx > i .and. tx <= N) then
!                 cv = 0.0d0
!                 do j = i + 1, tx
!                    cv = cv + T_s(IJ2TRI(j, i))*T_s(IJ2TRI(tx, j))
!                 end do
!              endif
!
!           endif
!
!           call syncthreads()
!           if (ty == 1 .and. tx > i .and. tx <= N) then
!              T_s(IJ2TRI(tx, i)) = cv
!           endif
!           call syncthreads()
!
!        end do
!
!        call syncthreads()
!
!        ! Write T_s to global
!        if (tx <= N) then
!           do j = ty, N, blockDim%y
!              if (tx >= j) then
!                 T(tx, j) = T_s(IJ2TRI(tx, j))
!              endif
!           end do
!        end if
!
!     end subroutine finish_T_block_kernel
!
function test_launch_finish_t_block_kernel()
  ! errorCode > 0 implies that the test has failed
  use iso_c_binding
  use hip
  use dsyevd_gpu_kernels
  implicit none
  integer :: errorCode = 1
  ! TODO fix parameters
  ! - Add missing arguments
  ! - Determine size of arrays (typically indicated by 'type(c_ptr)' type)
  ! - Add target where we need a pointer
  type(dim3,,intent(IN, :: grid
  type(dim3,,intent(IN, :: block
  integer(c_int,,intent(IN, :: sharedMem
  type(c_ptr,,value,intent(IN, :: stream
  integer,value :: n
  integer,value :: ldt
  type(c_ptr),value :: _t
  type(c_ptr),value :: _tau
  ! TODO Create initial data on host
  ! TODO Copy data to device ! (dest,src,size,direction)
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  CALL hipCheck(hipMemcpy(???,c_loc(???),C_SIZEOF(???),hipMemcpyHostToDevice)) !
  ! ... might be more (or less) than two memcopies
  ! TODO run the test
  CALL launch_finish_t_block_kernel(0,c_null_ptr,grid,block,sharedMem,stream,n,ldt,t,tau) ! Modify sharedMem if other than default 0
  CALL hipCheck(hipDeviceSynchronize())
  CALL launch_finish_t_block_kernel_cpu(0,c_null_ptr,grid,block,sharedMem,stream,n,ldt,t,tau)

  ! TODO Copy results back to host
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  CALL hipCheck(hipMemcpy(c_loc(???),???,C_SIZEOF(???),hipMemcpyDeviceToHost)
  ! ... might be more (or less) than two memcopies
  ! TODO Compare results
  ! TODO Update error code if the results do not match
  return errorCode
end function

program test_dsyevd_gpu_kernels
  implicit none
  integer :: globalErrorCode = 0, errorCode, fails = 0, tests = 0
  ! declare test functions and return type
  integer :: test_launch_krnl_e26a05_0_auto
  integer :: test_launch_krnl_b1f342_1_auto
  integer :: test_launch_krnl_b95769_2_auto
  integer :: test_launch_finish_t_block_kernel
  write(*,*) "SUITE test_dsyevd_gpu_kernels run ..."
  errorCode = test_launch_krnl_e26a05_0_auto()
  IF (errorCode > 0) THEN
    fails = fails + 1
    write(*,*) "TEST test_launch_krnl_e26a05_0_auto ... FAILURE"
  ELSE
    write(*,*) "TEST test_launch_krnl_e26a05_0_auto ... SUCCESS"
  END IF
  tests = tests + 1
  globalErrorCode = globalErrorCode + errorCode
  errorCode = test_launch_krnl_b1f342_1_auto()
  IF (errorCode > 0) THEN
    fails = fails + 1
    write(*,*) "TEST test_launch_krnl_b1f342_1_auto ... FAILURE"
  ELSE
    write(*,*) "TEST test_launch_krnl_b1f342_1_auto ... SUCCESS"
  END IF
  tests = tests + 1
  globalErrorCode = globalErrorCode + errorCode
  errorCode = test_launch_krnl_b95769_2_auto()
  IF (errorCode > 0) THEN
    fails = fails + 1
    write(*,*) "TEST test_launch_krnl_b95769_2_auto ... FAILURE"
  ELSE
    write(*,*) "TEST test_launch_krnl_b95769_2_auto ... SUCCESS"
  END IF
  tests = tests + 1
  globalErrorCode = globalErrorCode + errorCode
  errorCode = test_launch_finish_t_block_kernel()
  IF (errorCode > 0) THEN
    fails = fails + 1
    write(*,*) "TEST test_launch_finish_t_block_kernel ... FAILURE"
  ELSE
    write(*,*) "TEST test_launch_finish_t_block_kernel ... SUCCESS"
  END IF
  tests = tests + 1
  globalErrorCode = globalErrorCode + errorCode

  IF (globalErrorCode > 0) THEN
    write(*,*) "SUITE test_dsyevd_gpu_kernels ... FAILURE passed:",(tests-fails)," failed:",fails," total:",tests
  ELSE
    write(*,*) "SUITE test_dsyevd_gpu_kernels ... SUCCESS passed:",(tests-fails)," failed:",fails," total:",tests
  END IF
end program test_dsyevd_gpu_kernels


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2020 Guido Dhondt
!
!     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 2);
!
!
!     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, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine jouleheating(ipkon,lakon,kon,co,elcon,nelcon,
     &     mi,ne,sti,ielmat,nelemload,sideload,xload,nload,nload_,
     &     iamload,nam,idefload,ncmat_,ntmat_,
     &     alcon,nalcon,ithermal,vold,t1)
!
!     determines the effect of Joule heating
!
      implicit none
!
      character*8 lakon(*)
      character*20 label,sideload(*)
!
      integer ipkon(*),nelem,kon(*),mi(*),nope,indexe,j,k,null,
     &  mint3d,jj,iflag,ne,nelemload(2,*),iamload(2,*),nload,nload_,
     &  ielmat(mi(3),*),konl(20),idefload(*),iamplitude,isector,nam,
     &  one,nelcon(2,*),nalcon(2,*),ithermal(*),i1,ncmat_,ntmat_,imat
!
      real*8 co(3,*),xl(3,20),xi,et,ze,xsj,shp(4,20),weight,xload(2,*),
     &  sti(6,mi(1),*),alpha(6),heat,elcon(0:ncmat_,ntmat_,*),volume,
     &  elconloc(21),t1l,alcon(0:6,ntmat_,*),vold(0:mi(2),*),t1(*)
!
      include "gauss.f"
!
      data iflag /2/
!
      null=0
      one=1
!
      do nelem=1,ne
         if(ipkon(nelem).lt.0) cycle
!
!        only elements belonging to the A,V-domain experience
!        Joule heating
!
         if(int(elcon(2,1,ielmat(1,nelem))).ne.2) cycle
!
         imat=ielmat(1,nelem)
         indexe=ipkon(nelem)
!
         if(lakon(nelem)(1:5).eq.'C3D8I') then
            nope=11
         elseif(lakon(nelem)(4:4).eq.'2') then
            nope=20
         elseif(lakon(nelem)(4:4).eq.'8') then
            nope=8
         elseif(lakon(nelem)(4:5).eq.'10') then
            nope=10
         elseif(lakon(nelem)(4:4).eq.'4') then
            nope=4
         elseif(lakon(nelem)(4:5).eq.'15') then
            nope=15
         elseif(lakon(nelem)(4:5).eq.'6') then
            nope=6
         endif
!
         do j=1,nope
            konl(j)=kon(indexe+j)
            do k=1,3
               xl(k,j)=co(k,konl(j))
            enddo
         enddo
!
         heat=0.d0
         volume=0.d0
!
         if(lakon(nelem)(4:5).eq.'8R') then
            mint3d=1
         elseif((lakon(nelem)(4:4).eq.'8').or.
     &           (lakon(nelem)(4:6).eq.'20R')) then
            mint3d=8
         elseif(lakon(nelem)(4:4).eq.'2') then
            mint3d=27
         elseif(lakon(nelem)(4:5).eq.'10') then
            mint3d=4
         elseif(lakon(nelem)(4:4).eq.'4') then
            mint3d=1
         elseif(lakon(nelem)(4:5).eq.'15') then
            mint3d=9
         elseif(lakon(nelem)(4:5).eq.'6') then
            mint3d=2
         endif
!
!        loop over the integration points
!
         do jj=1,mint3d
            if(lakon(nelem)(4:5).eq.'8R') then
               xi=gauss3d1(1,jj)
               et=gauss3d1(2,jj)
               ze=gauss3d1(3,jj)
               weight=weight3d1(jj)
            elseif((lakon(nelem)(4:4).eq.'8').or.
     &              (lakon(nelem)(4:6).eq.'20R'))
     &              then
               xi=gauss3d2(1,jj)
               et=gauss3d2(2,jj)
               ze=gauss3d2(3,jj)
               weight=weight3d2(jj)
            elseif(lakon(nelem)(4:4).eq.'2') then
               xi=gauss3d3(1,jj)
               et=gauss3d3(2,jj)
               ze=gauss3d3(3,jj)
               weight=weight3d3(jj)
            elseif(lakon(nelem)(4:5).eq.'10') then
               xi=gauss3d5(1,jj)
               et=gauss3d5(2,jj)
               ze=gauss3d5(3,jj)
               weight=weight3d5(jj)
            elseif(lakon(nelem)(4:4).eq.'4') then
               xi=gauss3d4(1,jj)
               et=gauss3d4(2,jj)
               ze=gauss3d4(3,jj)
               weight=weight3d4(jj)
            elseif(lakon(nelem)(4:5).eq.'15') then
               xi=gauss3d8(1,jj)
               et=gauss3d8(2,jj)
               ze=gauss3d8(3,jj)
               weight=weight3d8(jj)
            else
               xi=gauss3d7(1,jj)
               et=gauss3d7(2,jj)
               ze=gauss3d7(3,jj)
               weight=weight3d7(jj)
            endif
!
            if(nope.eq.20) then
               call shape20h(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.8) then
               call shape8h(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.10) then
               call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.4) then
               call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
            elseif(nope.eq.15) then
               call shape15w(xi,et,ze,xl,xsj,shp,iflag)
            else
               call shape6w(xi,et,ze,xl,xsj,shp,iflag)
            endif
!
!     calculating the temperature
!
            t1l=0.d0
            if(ithermal(1).eq.1) then
               if(lakon(nelem)(4:5).eq.'8 ') then
                  do i1=1,nope
                     t1l=t1l+t1(konl(i1))/8.d0
                  enddo
               elseif(lakon(nelem)(4:6).eq.'20 ')then
                  call linscal(t1,konl,nope,jj,t1l,one)
               elseif(lakon(nelem)(4:6).eq.'10T') then
                  call linscal10(t1,konl,t1l,null,shp)
               else
                  do i1=1,nope
                     t1l=t1l+shp(4,i1)*t1(konl(i1))
                  enddo
               endif
            elseif(ithermal(1).ge.2) then
               if(lakon(nelem)(4:5).eq.'8 ') then
                  do i1=1,nope
                     t1l=t1l+vold(0,konl(i1))/8.d0
                  enddo
               elseif(lakon(nelem)(4:6).eq.'20 ')then
                  call linscal(vold,konl,nope,jj,t1l,mi(2))
               elseif(lakon(nelem)(4:6).eq.'10T') then
                  call linscal10(vold,konl,t1l,mi(2),shp)
               else
                  do i1=1,nope
                     t1l=t1l+shp(4,i1)*vold(0,konl(i1))
                  enddo
               endif
            endif
!
!        material data (electric conductivity and
!        magnetic permeability)
!
            call materialdata_em(elcon,nelcon,alcon,nalcon,
     &           imat,ntmat_,t1l,elconloc,ncmat_,alpha)
!
            heat=heat+weight*xsj*alpha(1)*
     &           (sti(1,jj,nelem)*sti(1,jj,nelem)+
     &           sti(2,jj,nelem)*sti(2,jj,nelem)+
     &           sti(3,jj,nelem)*sti(3,jj,nelem))
            volume=volume+weight*xsj
         enddo
!
         heat=heat/volume
!
!        adding the Joule heating to the distributed loading
!
         label='BF                  '
         iamplitude=0
         isector=0
         call loadadd(nelem,label,heat,nelemload,sideload,
     &      xload,nload,nload_,iamload,iamplitude,nam,isector,
     &     idefload)
      enddo
!
      return
      end

ext: f

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
! { dg-do compile }
!
! PR fortran/32905 - accepts types with ultimate POINTER components
!
MODULE types
  type :: tp3
    real :: x
    integer, pointer :: i
  end type

  type :: tp2
    type(tp3) :: t
  end type

  type :: tp1
    integer :: i
    type(tp2) :: t
  end type
END MODULE

MODULE nml
USE types
   type(tp1) :: t1
   type(tp3) :: t3

   namelist /a/ t1    ! { dg-error "cannot have POINTER components" }
   namelist /b/ t3    ! { dg-error "cannot have POINTER components" }
END MODULE

! { dg-final { cleanup-modules "types nml" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
!  ************************************************************************
!  *                                                                      *
!  *                 S U B R O U T I N E   E C H O                        *
!  *                                                                      *
!  *                                                         gm   1/21/05 *
!  *                                                                      *
!  ************************************************************************
!  Used in reading input data files.  Skips any line in a file which has
!  an "=" sign in the first column!

!       variable definitions:
!       variable  i/o   description
!       --------  ---   -----------

!       external references:  file associated with ===> inp

!       called from: main, input

!       programmer: ken rojas

      subroutine ECHO(inp)

!      implicit double precision (a-h,o-z)

!  Local Variables
      character*80 line
      integer inp

      if (inp .ne. 5) then

10       read (unit=inp,fmt=100,end=9999) line
         if (line(1:1) .eq. '=') go to 10
         backspace (unit=inp)

      endif

9999  return

100   format(a80)

      end

ext: for

max_stars_repo_name: USDA-ARS-WMSRU/phenologymms-science

max_issues_repo_name: USDA-ARS-WMSRU/phenologymms-science

max_forks_repo_name: USDA-ARS-WMSRU/phenologymms-science


text:
REAL*8 function kern_sqexp(r)
implicit none
REAL*8, intent(in) :: r
kern_sqexp = exp(-0.5d0*r**2)
end function
REAL*8 function dkern_sqexp(r)
implicit none
REAL*8, intent(in) :: r
dkern_sqexp = -r*exp(-0.5d0*r**2)
end function
REAL*8 function d2kern_sqexp(r)
implicit none
REAL*8, intent(in) :: r
d2kern_sqexp = (r**2 - 1)*exp(-0.5d0*r**2)
end function
REAL*8 function kern_matern32(r)
implicit none
REAL*8, intent(in) :: r
kern_matern32 = (sqrt(3.0d0)*r + 1)*exp(-1.7320508075688773d0*r)
end function
REAL*8 function dkern_matern32(r)
implicit none
REAL*8, intent(in) :: r
dkern_matern32 = -3*r*exp(-1.7320508075688773d0*r)
end function
REAL*8 function d2kern_matern32(r)
implicit none
REAL*8, intent(in) :: r
d2kern_matern32 = 3*(sqrt(3.0d0)*r - 1)*exp(-1.7320508075688773d0*r)
end function
REAL*8 function kern_matern52(r)
implicit none
REAL*8, intent(in) :: r
kern_matern52 = ((5.0d0/3.0d0)*r + sqrt(5.0d0)*r + 1)*exp( &
      -2.2360679774997897d0*r)
end function
REAL*8 function dkern_matern52(r)
implicit none
REAL*8, intent(in) :: r
dkern_matern52 = -1.0d0/3.0d0*(5*sqrt(5.0d0)*r + 15*r - 5)*exp( &
      -2.2360679774997897d0*r)
end function
REAL*8 function d2kern_matern52(r)
implicit none
REAL*8, intent(in) :: r
d2kern_matern52 = (5.0d0/3.0d0)*(5*r + 3*sqrt(5.0d0)*r - 2*sqrt(5.0d0) - &
      3)*exp(-2.2360679774997897d0*r)
end function
REAL*8 function kern_wend4(r)
implicit none
REAL*8, intent(in) :: r
kern_wend4 = (1 - r)**4*(4*r + 1)
end function
REAL*8 function dkern_wend4(r)
implicit none
REAL*8, intent(in) :: r
dkern_wend4 = 20*r*(r - 1)**3
end function
REAL*8 function d2kern_wend4(r)
implicit none
REAL*8, intent(in) :: r
d2kern_wend4 = (r - 1)**2*(80*r - 20)
end function

max_stars_repo_name: krystophny/profit

max_issues_repo_name: krystophny/profit

max_forks_repo_name: krystophny/profit


text:
The Hindustani Vocal Ensemble is a performance class (MUS 148) offered through the UC Davis Music Department. The class is instructed by renowned musician and Berkeley resident http://www.ritasahai.com/ Rita Sahai. The course focuses on the basics of Hindustani music through theory and practice.

The twounit class meets Mondays from 4 to 6 pm in the Music Building. Some people like the class so much that they keep taking it for no credit or after theyve graduated. The ensemble is open to all backgrounds and musical abilities  there are no prerequisite classes and students dont need previous experience. You wont find the class listed on the Open Courses website or the Course Catalog, but if you just show up on a Monday afternoon, you will be more than welcome to join.

For more information, contact Professor Sandra Graham at sjgraham@ucdavis.edu.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Founded in 1995, Studio66 is a graphic design & marketing firm, specializing in website design, branding, print and green marketing.

In 2005, they were ranked 29th out of the Sacramento Business Journals 100 Fastest Growing Companies.
In 2006, they were ranked 17th out of the same list.
In 2006, they received two silver ADDY Design Awards for http://www.mydtown.com dtown (a local publication) and a postcard campaign for Brooks Painting.

In 2008, they received a gold ADDY Design Award for the Scotts Valley Target Branding and campaign http://www.scottsvalleytarget.com. They also received a silver ADDY Design Award for a self promotion Backend Programming campaign.

In January 2009, they launched the new Davis Chamber of Commerce website.

In March 2009, the studio was purchased by employees Rashael Parker and Users/PeterAnselmo Peter Anselmo.

They serve the greater Sacramento area, and their client list includes the following:
  http://www.blakeshvac.com/ Blakes Heating and Air Conditioning
  Brooks Painting
  Butchs Board Shop
  Castellucci & Associates
  Cannery Park
  Central Park Gardens
  City of Davis
  Coldwell Banker  Doug Arnold Real Estate
  Davis Chamber of Commerce
  Davis Media Access
  Davis Police Officers Association
  Easton Painting
  Foundation for Teaching Economics
  HSP Painting Company
  Ink Monkey Graphics
  Jon Berkley Management
  Marrone Organic Innovations
  Music Matt
  NexGen
  NorthStar General Contractor
  OTECH
  Pacent Learning Solutions
  Plumbing M.D.
  Purves and Associates
  Regency Centers
  Silverado Country Club & Resort
  Snow Landscape
  Spirit of America (nonprofit)
  Target
  The History Project at UC Davis
  UC Davis
  Vacaville Magazine
  wiki:woodland:Woodland Chamber of Commerce
  West Sacramento Chamber of Commerce
  Yme National Breast Cancer Organization, Northern CA

Davis contains many companies providing business to business services.


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
!
!
!    NEMO - Numerical Engine (for) Multiphysics Operators
! Copyright (c) 2007, Stefano Toninel
!                     Gian Marco Bianchi  University of Bologna
!              David P. Schmidt    University of Massachusetts - Amherst
!              Salvatore Filippone University of Rome Tor Vergata
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
!     1. Redistributions of source code must retain the above copyright notice,
!        this list of conditions and the following disclaimer.
!     2. 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.
!     3. Neither the name of the NEMO project nor the names of its contributors
!        may be used to endorse or promote products derived from this software
!        without specific prior written permission.
!
! 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 OWNER 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.
!
!---------------------------------------------------------------------------------
!
! $Id$
!
! Description:
!    Adds to PDE the contribution of the time derivative of FLD * PHI.
!    Remark: FLD is optional.
!
SUBMODULE (op_ddt) scalar_pde_ddt_implementation
    IMPLICIT NONE

    CONTAINS

    MODULE PROCEDURE scalar_pde_ddt
    USE class_psblas, ONLY : psb_dpk_, mypnum_, sw_pde, psb_cd_get_local_rows, abort_psblas, psb_get_loc_to_glob
    USE class_dimensions, ONLY : dimensions, volume_, time_, OPERATOR(/=)
    USE class_mesh, ONLY : mesh, check_mesh_consistency
    USE tools_operators, ONLY : lhs_, pde_sign, size_blk
    IMPLICIT NONE
    !
    CHARACTER(len=*), PARAMETER :: op_name = 'SCALAR_PDE_DDT'
    INTEGER :: i, ic, ic_glob, ifirst, info, ncells, nel, nmax
    INTEGER, ALLOCATABLE :: ia(:), ja(:)
    INTEGER, ALLOCATABLE :: iloc_to_glob(:)
    REAL(psb_dpk_), ALLOCATABLE :: A(:), b(:)
    REAL(psb_dpk_), ALLOCATABLE :: fld_x_old(:)
    REAL(psb_dpk_), ALLOCATABLE :: phi_x_old(:)
    REAL(psb_dpk_) :: dtinv, fact, fsign, side_
    TYPE(dimensions) :: dim, dim_temp
    TYPE(mesh), POINTER :: msh => NULL()
    TYPE(mesh), POINTER :: msh_phi => NULL(), msh_fld => NULL()

    CALL sw_pde%tic()

    IF(mypnum_() == 0) THEN
        WRITE(*,*) '* ', TRIM(pde%name_()), ': applying the Time Derivative ',&
            & 'operator to the ', TRIM(phi%name_()), ' field'
    END IF

    ! Possible reinit of PDE
    CALL pde%reinit_pde()

    ! Is PHI cell-centered?
    IF(phi%on_faces_()) THEN
        WRITE(*,100) TRIM(op_name)
        CALL abort_psblas
    END IF

    ! Is FLD cell-centered?
    IF(PRESENT(fld)) THEN
        IF(fld%on_faces_()) THEN
            WRITE(*,100) TRIM(op_name)
            CALL abort_psblas
        END IF
    END IF

    ! Checks mesh consistency PDE vs. PHI
    CALL pde%get_mesh(msh)
    CALL phi%get_mesh(msh_phi)
    BLOCK
        USE class_mesh, ONLY : check_mesh_consistency
        CALL check_mesh_consistency(msh,msh_phi,op_name)
    END BLOCK

    ! Checks mesh consistency PHI vs. FLD
    IF(PRESENT(fld)) THEN
        CALL fld%get_mesh(msh_fld)
        BLOCK
            USE class_mesh, ONLY : check_mesh_consistency
            CALL check_mesh_consistency(msh_phi,msh_fld,op_name)
        END BLOCK
    END IF

    NULLIFY(msh_fld)
    NULLIFY(msh_phi)

    ! Equation dimensional check
    dim = phi%dim_() * volume_ / time_
    IF(PRESENT(fld)) dim = fld%dim_() * dim
    IF(dim /= pde%dim_()) THEN
        dim_temp=phi%dim_()
        !CALL phi%dim_()%debug_dim()
        CALL dim_temp%debug_dim()
        dim_temp=pde%dim_()
        !CALL pde%dim_()%debug_dim()
        CALL dim_temp%debug_dim()
        CALL dim%debug_dim()
        WRITE(*,200) TRIM(op_name)
        CALL abort_psblas
    END IF

    ! Computes sign factor
    IF(PRESENT(side)) THEN
        side_ = side
    ELSE
        side_ = lhs_ ! Default = LHS
    END IF
    fsign = pde_sign(sign,side_)

    ! Gets PHI "x" internal values
    CALL phi%get_x(phi_x_old)

    ! Gets FLD "x" internal values
    IF(PRESENT(fld)) THEN
        CALL fld%get_x(fld_x_old)
    ELSE
        ncells = SIZE(phi_x_old)
        ALLOCATE(fld_x_old(ncells),stat=info)
        IF(info /= 0) THEN
            WRITE(*,300) TRIM(op_name)
            CALL abort_psblas
        END IF
        fld_x_old(:) = 1.d0
    END IF

    ! Number of strictly local cells
    ncells = psb_cd_get_local_rows(msh%desc_c)

    ! Gets local to global list for cell indices
    CALL psb_get_loc_to_glob(msh%desc_c,iloc_to_glob)

    ! Computes maximum size of blocks to be inserted
    nmax = size_blk(1,ncells)

    ! Checks timestep size
    IF(dt <= 0.d0) THEN
        WRITE(*,400)
        CALL abort_psblas
    END IF

    dtinv = 1.d0 / dt

    ALLOCATE(A(nmax),b(nmax),ia(nmax),ja(nmax),stat=info)
    IF(info /= 0 ) THEN
        WRITE(*,300) TRIM(op_name)
        CALL abort_psblas
    END IF

    ifirst = 1; ic = 0
    insert: DO
        IF(ifirst > ncells) EXIT insert
        nel = size_blk(ifirst,ncells)

        BLOCK: DO i = 1, nel
            ! Local indices
            ic = ic + 1

            fact = fsign * msh%vol(ic) * dtinv * fld_x_old(ic)

            A(i) = fact
            b(i) = fact * phi_x_old(ic)

            ! Global indices in COO format
            ic_glob = iloc_to_glob(ic)
            ia(i)= ic_glob
            ja(i)= ic_glob
        END DO BLOCK

        CALL pde%spins_pde(nel,ia,ja,A)
        CALL pde%geins_pde(nel,ia,b)

        ifirst = ifirst +  nel

    END DO insert

    DEALLOCATE(A,b,ia,ja)
    DEALLOCATE(iloc_to_glob)

    DEALLOCATE(phi_x_old)
    DEALLOCATE(fld_x_old)
    NULLIFY(msh)

    CALL sw_pde%toc()

100 FORMAT(' ERROR! Operands in ',a,' are not cell centered')
200 FORMAT(' ERROR! Dimensional check failure in ',a)
300 FORMAT(' ERROR! Memory allocation failure in ',a)
400 FORMAT(' ERROR! Missing or illegal time advancing parameters')

    END PROCEDURE scalar_pde_ddt

END SUBMODULE scalar_pde_ddt_implementation

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


text:
      subroutine prompt (str, buf, fd)
      integer str(100), buf(100)
      integer fd
      integer isatty
      if (.not.(isatty(fd) .eq. 1))goto 23000
      call putlin (str, fd)
      call rfflus(fd)
23000 continue
      call getlin (buf, fd)
      return
      end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


text:
C Copyright(C) 1999-2020 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C See packages/seacas/LICENSE for details

C=======================================================================
      SUBROUTINE LNPLOT (NPTIMS, NPTIMW, XLN, YLN, ZLN, IXNODE, *)
C=======================================================================

C   --*** LNPLOT *** (PATHLN) Plot the pathlines
C   --   Written by Amy Gilkey - revised 05/27/88
C   --
C   --LNPLOT plots the pathlines.  For 3D meshes, the pathlines are rotated.
C   --
C   --Parameters:
C   --   NPTIMS - IN - the number of points on a history pathline
C   --   NPTIMW - IN - the number of points on a non-history pathline
C   --   XLN, YLN, ZLN - IN - the pathline data
C   --   IXNODE - SCRATCH - size = NPTIMS (3D only)
C   --   * - return statement if the cancel function is active
C   --
C   --Common Variables:
C   --   Uses NLNCRV, ILVID of /LNVARS/
C   --   Uses ROTMAT, ROTCEN of /ROTOPT/

      PARAMETER (NUMSYM = 6, NUMLIN = 6)

      include 'lnvars.blk'
      include 'd3nums.blk'
      include 'rotopt.blk'

      REAL XLN(NPTIMS,NLNCRV), YLN(NPTIMS,NLNCRV), ZLN(NPTIMS,NLNCRV)
      INTEGER IXNODE(NPTIMS)

      LOGICAL GRABRT
      LOGICAL NUMCRV
      CHARACTER TYP
      CHARACTER*8 LABSID

      lintyp = 1
      isytyp = 0
      numcrv = .false.
      labsid = 'NONE'

      IF (IS3DIM) THEN
         DO 100 I = 1, NPTIMS
            IXNODE(I) = I
  100    CONTINUE
      END IF

      DO 110 NP = 1, NLNCRV

         IF (GRABRT()) RETURN 1
         CALL GRCOLR (NP)
         CALL GRSYMB (LINTYP, ISYTYP, NP)

         CALL DBVTYP_BL (ILVID(1,NP), TYP, IDUM)
         IF (TYP .EQ. 'H') THEN
            NPTS = NPTIMS
         ELSE
            NPTS = NPTIMW
         END IF

C      --Bl_Rotate the 3D pathline

         IF (IS3DIM) THEN
            CALL BL_ROTATE (NPTS, IXNODE, ROTMAT, ROTCEN,
     &         XLN, YLN, ZLN, XLN, YLN, ZLN)
         END IF

C      --Plot pathline

         IF (GRABRT()) RETURN 1
         CALL PLTCUR (XLN(1,NP), YLN(1,NP), NPTS)

         IF (NUMCRV) THEN
            IF (GRABRT()) RETURN 1
            CALL GRNCRV (LABSID, NP, NPTS,
     &         XLN(1,NP), YLN(1,NP), (LINTYP .EQ. 0))
         END IF

C      --Finish plot

         CALL PLTFLU
  110 CONTINUE

      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
module SD_ale_volume_operator
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ZERO,ONE,TWO,HALF

    use type_operator,          only: operator_t
    use type_chidg_worker,      only: chidg_worker_t
    use type_properties,        only: properties_t
    use DNAD_D
    use ieee_arithmetic
    implicit none
    private

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   9/14/2016
    !!
    !!
    !!
    !-------------------------------------------------------------------------
    type, extends(operator_t), public :: SD_ale_volume_operator_t

    contains

        procedure   :: init
        procedure   :: compute

    end type SD_ale_volume_operator_t
    !*************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/29/2016
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(SD_ale_volume_operator_t),   intent(inout)      :: self

        ! Set operator name
        call self%set_name('Scalar Diffusion ALE Volume Operator')

        ! Set operator type
        call self%set_operator_type('Volume Diffusive Operator')

        ! Set operator equations
        call self%add_primary_field('u')

    end subroutine init
    !********************************************************************************

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   9/14/2016
    !!
    !!
    !!
    !------------------------------------------------------------------------------------
    subroutine compute(self,worker,prop)
        class(SD_ale_volume_operator_t),    intent(inout)   :: self
        type(chidg_worker_t),           intent(inout)   :: worker
        class(properties_t),            intent(inout)   :: prop

        type(AD_D), allocatable, dimension(:)   ::  &
            flux_1, flux_2, flux_3, mu

        type(AD_D), allocatable, dimension(:,:)   ::  &
            gradu, &
            flux_ref

        !
        ! Interpolate solution to quadrature nodes
        !
        gradu = worker%get_primary_field_grad_ale_element('u','gradient + lift')

        !
        ! Compute scalar coefficient
        !
        mu = worker%get_model_field_element('Scalar Diffusion Coefficient', 'value')

        !
        ! Compute volume flux at quadrature nodes
        !
        flux_1 = -mu*gradu(:,1)
        flux_2 = -mu*gradu(:,2)
        flux_3 = -mu*gradu(:,3)

        !
        ! Multiply by grid jacobian
        !
        flux_ref = worker%post_process_volume_diffusive_flux_ale(flux_1, flux_2, flux_3)

        !
        ! Integrate volume flux
        !
        call worker%integrate_volume('u',flux_ref(:,1),flux_ref(:,2),flux_ref(:,3))

    end subroutine compute
    !****************************************************************************************************

end module SD_ale_volume_operator

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG


text:
20091111 09:32:13 nbsp Actually, I am an official sponsor of the DBC and their race team.  If you were a member, you would see my logo on their race kit on both sides of the shorts.  As for my sponsoring the UCD Cycling team.  I was supposed to be this years Bike supplier, but at the last minute, they went with an online retailer instead.  I still offer the UCD Cycling team a very nice discount in the shop (Better than any other shop in town, including the UCD Bike Barn, their main sponsor).  I also sponsor a grass roots program that has brought up many young riders into the ranks of collegiate and USAC/UCI events.  Many are racing as As, Bs, and Cs this year. My riders have had podium finishes in many local and national events. Please refrain from editing the APEX Cycles page in regards to my business practices. Users/Aaron.Curtin

20091202 18:08:19 nbsp Just so you know... I am now an official sponsor of the UCD cycling team this year.  Last week they requested my logos for their new jerseys for the upcoming season. Users/Aaron.Curtin


text:
C->>> ------------------------------------------> ems_ca_u_lg_dvx_wt <<<
c     Calls the routines to update Devex weights and possibly
c     dual activities.
c
      subroutine ems_ca_u_lg_dvx_wt(
     &     u_du_act,
     &     vr_in_c,
     &     du_act,
     &     pi_v,
     &     pi_ix,
     &     ed_wt)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RSMICS.INC'
      include 'RSMICOM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      logical u_du_act
      integer vr_in_c(-vr_in_c_n_sn:n_c)
      integer pi_ix(0:n_r)
      double precision du_act(0:mx_n_c+n_r)
      double precision pi_v(0:n_r)
      double precision ed_wt(0:mx_n_c+n_r)

      if (vr_t_en_bs .le. 0 .or. vr_t_lv_bs .eq. 0 .or.
     &     vr_t_en_bs .eq. vr_t_lv_bs) then
         if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9900)
     &        vr_t_en_bs, vr_t_lv_bs
         call ems_msg_wr_li(warn_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?         call ems_dump
CM      ENDIF
         go to 7000
      end if
      if (u_du_act) then
         call ems_u_lg_dvx_wt_du_act(
     &        vr_t_en_bs,
     &        vr_in_c,
     &        du_act,
     &        pi_v,
     &        pi_ix,
     &        ed_wt)
      else
         call ems_u_lg_dvx_wt(
     &        vr_t_en_bs,
     &        vr_in_c,
     &        pi_v,
     &        pi_ix,
     &        ed_wt)
      endif
      if (vr_t_lv_bs .gt. mx_n_c) then
c
c     Set the weight for the variable which has just left the basis.
c
         ed_wt(vr_t_lv_bs) = max(one, ed_wt(vr_t_en_bs)/abs(pv))
      endif
 7000 continue
      return
 9900 format('Calling u_lg_dvx_wt with vr_t_en_bs, vr_t_lv_bs = ',
     &     i9, i9)
      end

C->>> ---------------------------------------------> ems_u_lg_dvx_wt <<<
c     Updates Devex weights for logicals.
c
      subroutine ems_u_lg_dvx_wt(
     &     vr_t_en_bs,
     &     vr_in_c,
     &     pi_v,
     &     pi_ix,
     &     ed_wt)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
CM      IF (emsol_tt .EQ. 1) THEN
C?      include 'EMSTT.INC'
CM      ENDIF
      include 'ICTVR.INC'
      integer vr_t_en_bs, vr_in_c(-vr_in_c_n_sn:n_c)
      integer pi_ix(0:n_r)
      double precision pi_v(0:n_r)
      double precision ed_wt(0:mx_n_c+n_r)
      integer c_n, vr_n, ix_n, r_n
      integer c_loop_ln
      double precision pi_v_mu

      pi_v_mu = ed_wt(vr_t_en_bs)
      c_loop_ln = vr_in_c(os_lg_in_c_l_pc_p)
      if (pi_ix(0) .gt. n_r .or.
     &     tbu_r_loop_mode .eq. tbu_r_loop_no .or.
     &     (tbu_r_loop_mode .eq. tbu_r_loop_poss .and.
     &     2*c_loop_ln .le. pi_ix(0))) then
c
c     Update the weights for the nonbasic logicals. See
c     subroutine ems_u_sed_wt for comments on obtaining the
c     reduced cost of the variable which has just left the basis.
c
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1) call ems_tt_rec(u_lg_dvx_wt_dse_tt, n_bs)
CM      ENDIF
         do 10 c_n = 1, vr_in_c(os_lg_in_c_l_pc_p)
            vr_n = vr_in_c(c_n)
            if (pi_v(vr_n-mx_n_c) .ne. zero) ed_wt(vr_n) =
     &           max(ed_wt(vr_n), abs(pi_v(vr_n-mx_n_c))*pi_v_mu)
 10      continue
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1) call ems_tt_rec(-u_lg_dvx_wt_dse_tt, n_bs)
CM      ENDIF
      else
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1) call ems_tt_rec(u_lg_dvx_wt_sps_tt, n_bs)
CM      ENDIF
         do 20, ix_n = 1, pi_ix(0)
            r_n = pi_ix(ix_n)
            vr_n = mx_n_c + r_n
            ed_wt(vr_n) = max(ed_wt(vr_n), abs(pi_v(r_n))*pi_v_mu)
 20      continue
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1) call ems_tt_rec(-u_lg_dvx_wt_sps_tt, n_bs)
CM      ENDIF
      endif
      return
      end

C->>> --------------------------------------> ems_u_lg_dvx_wt_du_act <<<
c     Updates Devex weights and dual activities for logicals.
c
      subroutine ems_u_lg_dvx_wt_du_act(
     &     vr_t_en_bs,
     &     vr_in_c,
     &     du_act,
     &     pi_v,
     &     pi_ix,
     &     ed_wt)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
CM      IF (emsol_tt .EQ. 1) THEN
C?      include 'EMSTT.INC'
CM      ENDIF
      include 'ICTVR.INC'
      integer vr_t_en_bs, vr_in_c(-vr_in_c_n_sn:n_c)
      integer pi_ix(0:n_r)
      double precision du_act(0:mx_n_c+n_r)
      double precision pi_v(0:n_r)
      double precision ed_wt(0:mx_n_c+n_r)
      integer c_n, vr_n, ix_n, r_n
      integer c_loop_ln
      double precision pi_v_mu

      pi_v_mu = ed_wt(vr_t_en_bs)/abs(du_act(vr_t_en_bs))
      c_loop_ln = vr_in_c(os_lg_in_c_l_pc_p)
      if (pi_ix(0) .gt. n_r .or.
     &     tbu_r_loop_mode .eq. tbu_r_loop_no .or.
     &     (tbu_r_loop_mode .eq. tbu_r_loop_poss .and.
     &     2*c_loop_ln .le. pi_ix(0))) then
c
c     Update the weights and dual activities for the nonbasic variables.
c     See subroutine ems_u_sed_wt for comments on obtaining the
c     reduced cost of the variable which has just left the basis.
c
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1) call ems_tt_rec(u_lg_dvx_wt_dse_tt, n_bs)
CM      ENDIF
         do 10 c_n = 1, vr_in_c(os_lg_in_c_l_pc_p)
            vr_n = vr_in_c(c_n)
            if (pi_v(vr_n-mx_n_c) .ne. zero) then
               ed_wt(vr_n) =
     &              max(ed_wt(vr_n), abs(pi_v(vr_n-mx_n_c))*pi_v_mu)
               du_act(vr_n) = du_act(vr_n) - pi_v(vr_n-mx_n_c)
            end if
 10      continue
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1) call ems_tt_rec(-u_lg_dvx_wt_dse_tt, n_bs)
CM      ENDIF
      else
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1) call ems_tt_rec(u_lg_dvx_wt_sps_tt, n_bs)
CM      ENDIF
         do 20, ix_n = 1, pi_ix(0)
            r_n = pi_ix(ix_n)
            vr_n = mx_n_c + r_n
            ed_wt(vr_n) = max(ed_wt(vr_n), abs(pi_v(r_n))*pi_v_mu)
            du_act(vr_n) = du_act(vr_n) - pi_v(r_n)
 20      continue
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1) call ems_tt_rec(-u_lg_dvx_wt_sps_tt, n_bs)
CM      ENDIF
      endif
      return
      end

C->>> -------------------------------------> ems_perm_ca_u_lg_dvx_wt <<<
c     Calls the routines to update Devex weights and possibly
c     dual activities.
c
      subroutine ems_perm_ca_u_lg_dvx_wt(
     &     u_du_act,
     &     vr_in_c,
     &     du_act,
     &     pi_v,
     &     pi_ix,
     &     ed_wt,
     &     og_t_nw_perm,
     &     nw_t_og_perm)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RSMICS.INC'
      include 'RSMICOM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      logical u_du_act
      integer vr_in_c(-vr_in_c_n_sn:n_c)
      integer og_t_nw_perm(0:n_r)
      integer nw_t_og_perm(0:n_r)
      integer pi_ix(0:n_r)
      double precision du_act(0:mx_n_c+n_r)
      double precision pi_v(0:n_r)
      double precision ed_wt(0:mx_n_c+n_r)

      if (vr_t_en_bs .le. 0 .or. vr_t_lv_bs .eq. 0 .or.
     &     vr_t_en_bs .eq. vr_t_lv_bs) then
         if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9900)
     &        vr_t_en_bs, vr_t_lv_bs
         call ems_msg_wr_li(warn_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?         call ems_dump
CM      ENDIF
         go to 7000
      end if
      if (u_du_act) then
         call ems_perm_u_lg_dvx_wt_du_act(
     &        vr_t_en_bs,
     &        vr_in_c,
     &        du_act,
     &        pi_v,
     &        pi_ix,
     &        ed_wt,
     &        og_t_nw_perm,
     &        nw_t_og_perm)
      else
         call ems_perm_u_lg_dvx_wt(
     &        vr_t_en_bs,
     &        vr_in_c,
     &        pi_v,
     &        pi_ix,
     &        ed_wt,
     &        og_t_nw_perm,
     &        nw_t_og_perm)
      endif
      if (vr_t_lv_bs .gt. mx_n_c) then
c
c     Set the weight for the variable which has just left the basis.
c
         ed_wt(vr_t_lv_bs) = max(one, ed_wt(vr_t_en_bs)/abs(pv))
      endif
 7000 continue
      return
 9900 format('Calling u_lg_dvx_wt with vr_t_en_bs, vr_t_lv_bs = ',
     &     i9, i9)
      end

C->>> ----------------------------------------> ems_perm_u_lg_dvx_wt <<<
c     Updates Devex weights for logicals.
c
      subroutine ems_perm_u_lg_dvx_wt(
     &     vr_t_en_bs,
     &     vr_in_c,
     &     pi_v,
     &     pi_ix,
     &     ed_wt,
     &     og_t_nw_perm,
     &     nw_t_og_perm)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
      include 'ICTVR.INC'
      integer vr_t_en_bs, vr_in_c(-vr_in_c_n_sn:n_c)
      integer og_t_nw_perm(0:n_r)
      integer nw_t_og_perm(0:n_r)
      integer pi_ix(0:n_r)
      double precision pi_v(0:n_r)
      double precision ed_wt(0:mx_n_c+n_r)
      integer c_n, vr_n, nw_r_n, ix_n, r_n, og_r_n
      integer c_loop_ln
      double precision pi_v_mu

      pi_v_mu = ed_wt(vr_t_en_bs)
      c_loop_ln = vr_in_c(os_lg_in_c_l_pc_p)
      if (pi_ix(0) .gt. n_r .or.
     &     tbu_r_loop_mode .eq. tbu_r_loop_no .or.
     &     (tbu_r_loop_mode .eq. tbu_r_loop_poss .and.
     &     2*c_loop_ln .le. pi_ix(0))) then
c
c     Update the weights for the nonbasic logicals. See
c     subroutine ems_u_sed_wt for comments on obtaining the
c     reduced cost of the variable which has just left the basis.
c
         do 10, c_n = 1, vr_in_c(os_lg_in_c_l_pc_p)
            vr_n = vr_in_c(c_n)
            nw_r_n = og_t_nw_perm(vr_n-mx_n_c)
            if (pi_v(nw_r_n) .ne. zero) ed_wt(vr_n) =
     &           max(ed_wt(vr_n), abs(pi_v(nw_r_n))*pi_v_mu)
 10      continue
      else
         do 20, ix_n = 1, pi_ix(0)
            r_n = pi_ix(ix_n)
            og_r_n = nw_t_og_perm(r_n)
            vr_n = mx_n_c + og_r_n
            ed_wt(vr_n) = max(ed_wt(vr_n), abs(pi_v(r_n))*pi_v_mu)
 20      continue
      endif
      return
      end

C->>> ---------------------------------> ems_perm_u_lg_dvx_wt_du_act <<<
c     Updates Devex weights and dual activities for logicals.
c
      subroutine ems_perm_u_lg_dvx_wt_du_act(
     &     vr_t_en_bs,
     &     vr_in_c,
     &     du_act,
     &     pi_v,
     &     pi_ix,
     &     ed_wt,
     &     og_t_nw_perm,
     &     nw_t_og_perm)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
      include 'ICTVR.INC'
      integer vr_t_en_bs, vr_in_c(-vr_in_c_n_sn:n_c)
      integer og_t_nw_perm(0:n_r)
      integer nw_t_og_perm(0:n_r)
      integer pi_ix(0:n_r)
      double precision du_act(0:mx_n_c+n_r)
      double precision pi_v(0:n_r)
      double precision ed_wt(0:mx_n_c+n_r)
      integer c_n, vr_n, nw_r_n, ix_n, r_n, og_r_n
      integer c_loop_ln
      double precision pi_v_mu

      pi_v_mu = ed_wt(vr_t_en_bs)/abs(du_act(vr_t_en_bs))
      c_loop_ln = vr_in_c(os_lg_in_c_l_pc_p)
      if (pi_ix(0) .gt. n_r .or.
     &     tbu_r_loop_mode .eq. tbu_r_loop_no .or.
     &     (tbu_r_loop_mode .eq. tbu_r_loop_poss .and.
     &     2*c_loop_ln .le. pi_ix(0))) then
c
c     Update the weights and dual activities for the nonbasic logicals.
c     See subroutine ems_u_sed_wt for comments on obtaining the
c     reduced cost of the variable which has just left the basis.
c
         do 10 c_n = 1, vr_in_c(os_lg_in_c_l_pc_p)
            vr_n = vr_in_c(c_n)
            nw_r_n = og_t_nw_perm(vr_n-mx_n_c)
            if (pi_v(nw_r_n) .ne. zero) then
               ed_wt(vr_n) = max(ed_wt(vr_n), abs(pi_v(nw_r_n))*pi_v_mu)
               du_act(vr_n) = du_act(vr_n) - pi_v(nw_r_n)
            end if
 10      continue
      else
         do 20, ix_n = 1, pi_ix(0)
            r_n = pi_ix(ix_n)
            og_r_n = nw_t_og_perm(r_n)
            vr_n = mx_n_c + og_r_n
            ed_wt(vr_n) = max(ed_wt(vr_n), abs(pi_v(r_n))*pi_v_mu)
            du_act(vr_n) = du_act(vr_n) - pi_v(r_n)
 20      continue
      endif
      return
      end

C->>> ---------------------------------------> ems_ca_u_struc_dvx_wt <<<
c     Calls the routines to update Devex weights and possibly
c     dual activities.
c
      subroutine ems_ca_u_struc_dvx_wt(
     &     u_du_act,
     &     vr_in_c,
     &     du_act,
     &     tbu_r_v,
     &     tbu_r_ix,
     &     ed_wt)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RSMICS.INC'
      include 'RSMICOM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      logical u_du_act
      integer vr_in_c(-vr_in_c_n_sn:n_c)
      double precision du_act(0:mx_n_c+n_r)
      double precision tbu_r_v(0:n_c)
      integer tbu_r_ix(0:n_c)
      double precision ed_wt(0:mx_n_c+n_r)

      if (vr_t_en_bs .le. 0 .or. vr_t_lv_bs .eq. 0 .or.
     &     vr_t_en_bs .eq. vr_t_lv_bs) then
         if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9900)
     &        vr_t_en_bs, vr_t_lv_bs
         call ems_msg_wr_li(warn_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?         call ems_dump
CM      ENDIF
         go to 7000
      end if
      if (u_du_act) then
         call ems_u_struc_dvx_wt_du_act(
     &        vr_in_c,
     &        du_act,
     &        tbu_r_v,
     &        tbu_r_ix,
     &        ed_wt)
      else
         call ems_u_struc_dvx_wt(
     &        vr_in_c,
     &        tbu_r_v,
     &        tbu_r_ix,
     &        ed_wt)
      endif
      if (vr_t_lv_bs .le. n_c) then
c
c     Set the weight for the variable which has just left the basis.
c
         ed_wt(vr_t_lv_bs) = max(one, ed_wt(vr_t_en_bs)/abs(pv))
      endif
      n_dvx_it = n_dvx_it + 1
      if (iand(ck_msk, ze_a_ck_bt) .ne. 0)
     &     call ems_ck_ze_rl_a(n_c, tbu_r_v)
 7000 continue
      return
 9900 format('Calling u_struc_dvx_wt with vr_t_en_bs, vr_t_lv_bs = ',
     &     i9, i9)
      end

C->>> ------------------------------------------> ems_u_struc_dvx_wt <<<
c     Updates Devex weights for structurals.
c
      subroutine ems_u_struc_dvx_wt(
     &     vr_in_c,
     &     tbu_r_v,
     &     tbu_r_ix,
     &     ed_wt)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
CM      IF (emsol_tt .EQ. 1) THEN
C?      include 'EMSTT.INC'
CM      ENDIF
      include 'RSMICOM.INC'
      include 'ICTVR.INC'
      integer vr_in_c(-vr_in_c_n_sn:n_c)
      double precision tbu_r_v(0:n_c)
      integer tbu_r_ix(0:n_c)
      double precision ed_wt(0:mx_n_c+n_r)
      integer c_loop_ln, c_n, vr_n, ix_n
      double precision pi_v_mu
      integer tbu_r_n_nz, tbu_r_n_c

      pi_v_mu = ed_wt(vr_t_en_bs)
      c_loop_ln =
     &     vr_in_c(os_struc_in_c_l_pc_p) -
     &     vr_in_c(os_struc_in_c_f_p_m1)
      tbu_r_n_c = c_loop_ln
c
c     Update the weights for the nonbasic variables. See
c     subroutine ems_u_sed_wt for comments on obtaining the
c     reduced cost of the variable which has just left the basis.
c
      if (tbu_r_ix(0) .gt. n_c .or.
     &     tbu_r_loop_mode .eq. tbu_r_loop_no .or.
     &     (tbu_r_loop_mode .eq. tbu_r_loop_poss .and.
     &     2*c_loop_ln .le. tbu_r_ix(0))) then
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1)
C?     &        call ems_tt_rec(u_struc_dvx_wt_dse_tt, n_bs)
CM      ENDIF
         tbu_r_n_nz = 0
         do 10, c_n = vr_in_c(os_struc_in_c_f_p_m1) + 1,
     &        vr_in_c(os_struc_in_c_l_pc_p)
            vr_n = vr_in_c(c_n)
            if (tbu_r_v(vr_n) .ne. zero) then
               tbu_r_n_nz = tbu_r_n_nz + 1
               ed_wt(vr_n) =
     &              max(ed_wt(vr_n), pi_v_mu*abs(tbu_r_v(vr_n)))
               tbu_r_v(vr_n) = zero
            endif
 10      continue
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1)
C?     &        call ems_tt_rec(-u_struc_dvx_wt_dse_tt, n_bs)
CM      ENDIF
      else
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1)
C?     &        call ems_tt_rec(u_struc_dvx_wt_sps_tt, n_bs)
CM      ENDIF
         tbu_r_n_nz = tbu_r_ix(0)
         do 20, ix_n = 1, tbu_r_ix(0)
            vr_n = tbu_r_ix(ix_n)
            ed_wt(vr_n) = max(ed_wt(vr_n), pi_v_mu*abs(tbu_r_v(vr_n)))
            tbu_r_v(vr_n) = zero
 20      continue
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1)
C?     &        call ems_tt_rec(-u_struc_dvx_wt_sps_tt, n_bs)
CM      ENDIF
      endif
c
c     Get the density of the tableau row (if required).
c
      if (tbu_r_dse .lt. zero) then
         if (tbu_r_n_c .le. 0) then
c
c     This situation---no nonbasic structurals---is pretty unlikely.
c     Setting tbu_r_dse to an illegal value (rather than 0 or 1) seems
c     sensible.
c
            tbu_r_dse = two
         else
            tbu_r_dse = float(tbu_r_n_nz)/float(tbu_r_n_c)
         endif
      endif
      return
      end

C->>> -----------------------------------> ems_u_struc_dvx_wt_du_act <<<
c     Updates Devex weights and dual activities for structurals.
c
      subroutine ems_u_struc_dvx_wt_du_act(
     &     vr_in_c,
     &     du_act,
     &     tbu_r_v,
     &     tbu_r_ix,
     &     ed_wt)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
CM      IF (emsol_da .EQ. 1) THEN
C?      include 'EMSDA.INC'
CM      ENDIF
CM      IF (emsol_tt .EQ. 1) THEN
C?      include 'EMSTT.INC'
CM      ENDIF
      include 'RSMICOM.INC'
      include 'ICTVR.INC'
      integer vr_in_c(-vr_in_c_n_sn:n_c)
      double precision du_act(0:mx_n_c+n_r)
      double precision tbu_r_v(0:n_c)
      integer tbu_r_ix(0:n_c)
      double precision ed_wt(0:mx_n_c+n_r)
      integer c_n, vr_n, c_loop_ln, ix_n
      double precision tbu_r_v_mu
      integer tbu_r_n_nz, tbu_r_n_c

      tbu_r_v_mu = ed_wt(vr_t_en_bs)/abs(du_act(vr_t_en_bs))
      c_loop_ln =
     &     vr_in_c(os_struc_in_c_l_pc_p) -
     &     vr_in_c(os_struc_in_c_f_p_m1)
      tbu_r_n_c = c_loop_ln
c
c     Update the weights and dual activities for the nonbasic variables.
c     See subroutine ems_u_sed_wt for comments on obtaining the
c     reduced cost of the variable which has just left the basis.
c
      if (tbu_r_ix(0) .gt. n_c .or.
     &     tbu_r_loop_mode .eq. tbu_r_loop_no .or.
     &     (tbu_r_loop_mode .eq. tbu_r_loop_poss .and.
     &     2*c_loop_ln .le. tbu_r_ix(0))) then
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1)
C?     &        call ems_tt_rec(u_struc_dvx_wt_dse_tt, n_bs)
CM      ENDIF
         tbu_r_n_nz = 0
         do 10, c_n = vr_in_c(os_struc_in_c_f_p_m1) + 1,
     &        vr_in_c(os_struc_in_c_l_pc_p)
            vr_n = vr_in_c(c_n)
            if (tbu_r_v(vr_n) .ne. zero) then
               tbu_r_n_nz = tbu_r_n_nz + 1
               ed_wt(vr_n) = max(
     &              ed_wt(vr_n), abs(tbu_r_v(vr_n))*tbu_r_v_mu)
               du_act(vr_n) = du_act(vr_n) + tbu_r_v(vr_n)
               tbu_r_v(vr_n) = zero
            end if
 10      continue
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1)
C?     &        call ems_tt_rec(-u_struc_dvx_wt_dse_tt, n_bs)
CM      ENDIF
      else
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1)
C?     &        call ems_tt_rec(u_struc_dvx_wt_sps_tt, n_bs)
CM      ENDIF
         tbu_r_n_nz = tbu_r_ix(0)
         do 20, ix_n = 1, tbu_r_ix(0)
            vr_n = tbu_r_ix(ix_n)
            if (tbu_r_v(vr_n) .ne. zero) then
               ed_wt(vr_n) = max(
     &              ed_wt(vr_n), abs(tbu_r_v(vr_n))*tbu_r_v_mu)
               du_act(vr_n) = du_act(vr_n) + tbu_r_v(vr_n)
               tbu_r_v(vr_n) = zero
            endif
 20      continue
CM      IF (emsol_tt .EQ. 1) THEN
C?         if (ems_tt_pc_lvl1)
C?     &        call ems_tt_rec(-u_struc_dvx_wt_sps_tt, n_bs)
CM      ENDIF
      endif
c
c     Get the density of the tableau row (if required).
c
      if (tbu_r_dse .lt. zero) then
         if (tbu_r_n_c .le. 0) then
c
c     This situation---no nonbasic structurals---is pretty unlikely.
c     Setting tbu_r_dse to an illegal value (rather than 0 or 1) seems
c     sensible.
c
            tbu_r_dse = two
         else
            tbu_r_dse = float(tbu_r_n_nz)/float(tbu_r_n_c)
         endif
      endif
      return
      end

C->>> ----------------------------------------------> ems_iz_dvx_fwk <<<
c     Initialises the weights and the index sets for Devex.
c
      subroutine ems_iz_dvx_fwk(
     &     vr_in_r, vr_in_c,
     &     ed_wt,
     &     dvx_ix)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RSMICS.INC'
      include 'ICTVR.INC'
      integer vr_in_r(0:n_r)
      integer vr_in_c(-vr_in_c_n_sn:n_c)
      integer dvx_ix(0:mx_n_c+n_r)
      double precision ed_wt(0:mx_n_c+n_r)
      integer c_n, vr_n, r_n
c
c     Set R is the reference framework and consists of the set of
c     variables for pricing for this simplex iteration.
c     Set H is the set indices which are used to compute the Devex
c     norm. It consists of the basic variables in R so, initially, is
c     empty.
c     The implementation requires the vector
c
c     dvx_ix     If ix = dvx_ix(vr_n) then
c                   ix = 1 => the variable vr_n is in H (in R & basic)
c                   ix = 0 => the variable vr_n is in R but not basic.
c                   ix =-1 => the variable vr_n is not in R.
c
c
c     Zero the logical edge weights in case values of pi are just
c     used to update them without checking whether the variable is basic
c     or nonbasic---thus causing an unassigned variable error with EPC.
c
      call ems_cp_rl_a(n_r, zero, ed_wt(mx_n_c+1), 0)
c
c     Original initialisation didn't assign weights to variables not being priced
c     but this led to FPR when such weights for logicals were updated when
c     passing through pi---even if the weights were never used. Probably safe to
c     delete loops 10 and 20, but they are commented out for reference.
c
c      do 10, c_n = 1, vr_in_c(os_lg_in_c_l_pc_p)
c         vr_n = vr_in_c(c_n)
c         ed_wt(vr_n) = one
c         dvx_ix(vr_n) = 0
c 10   continue
c      do 20, c_n = vr_in_c(os_struc_in_c_f_p_m1) + 1,
c     &     vr_in_c(os_struc_in_c_l_pc_p)
c         vr_n = vr_in_c(c_n)
c         ed_wt(vr_n) = one
c         dvx_ix(vr_n) = 0
c     20   continue
      do 10, c_n = 1, n_c
         vr_n = vr_in_c(c_n)
         ed_wt(vr_n) = one
         dvx_ix(vr_n) = 0
 10   continue
      do 30, r_n = 1, n_r
         vr_n = vr_in_r(r_n)
         ed_wt(vr_n) = zero
         dvx_ix(vr_n) = -1
 30   continue
      n_dvx_it = 1
      n_dvx_fwk = n_dvx_fwk + 1
c
c     Indicate that the edge weight information is correct.
c
      ml_da_st_msk = ior(ml_da_st_msk, ml_da_st_ed_wt)
      return
      end

C->>> ----------------------------------------------> ems_iz_blk_dvx <<<
c     Sets up a block for the vector required by Devex.
c
      subroutine ems_iz_blk_dvx(is)
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'EMSPM.INC'
      include 'EMSMEM.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      integer is(0:is_n_en_m1)
      integer blk_n, r_cf, c_cf, a_el_cf, cs, n_wo, p_ml_hdl
      integer mem_mgr_rt_cod
      integer rt_cod

      if (iand(ml_blk_st_msk, ml_blk_st_dvx) .ne. 0) goto 8000
      call ems_g_blk_dvx_n_wo(r_cf, c_cf, a_el_cf, cs, n_wo)
      call ems_mem_mgr_ope_blk(mem_mgr_rt_cod, is,
     &     n_wo, ope_blk_anywhere,
     &     cu_ml_n, dvx_blk_id, blk_n)
      if (mem_mgr_rt_cod .ne. mem_mgr_rt_cod_ok) then
         if (mem_mgr_rt_cod .eq. mem_mgr_rt_cod_serious_no_po) then
            ems_msg_cod = ems_msg_lvl_serious
            go to 7000
         else if (mem_mgr_rt_cod .ge. mem_mgr_rt_lvl_serious) then
             ems_msg_cod = ems_msg_lvl_serious
            go to 7000
         endif
      endif
      call ems_mem_mgr_rp_ope_blk(
     &     mem_mgr_rt_cod, is, -1,
     &     r_cf, mx_n_r,
     &     c_cf, mx_n_c,
     &     0, 0,
     &     cs, n_wo, dvx_blk_id)
      if (mem_mgr_rt_cod .ge. mem_mgr_rt_lvl_serious) then
         ems_msg_cod = ems_msg_lvl_serious
         go to 7000
      endif
      p_ml_hdl = p_ml_bs_blk + (cu_ml_n-1)*ml_bs_blk_n_wo +
     &     ml_bs_blk_os_hdl
      call ems_mem_mgr_nw_hdl(mem_mgr_rt_cod, is,
     &     blk_n, 1+mx_n_c+mx_n_r,
     &     i_wo_z, is(p_ml_hdl+ix_hdl_dvx_ix))
      if (ems_msg_cod .ge. ems_msg_lvl_serious) go to 7000
      call ems_g_ml_p(rt_cod, is)
      if (rt_cod .ne. 0) goto 8020
      ml_blk_st_msk = ml_blk_st_msk + ml_blk_st_dvx
 7000 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 8020 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9802)
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 9800 format('Model already has space for DEVEX')
 9802 format('Error in ems_g_ml_p')
      end

C->>> ------------------------------------------> ems_g_blk_dvx_n_wo <<<
      subroutine ems_g_blk_dvx_n_wo(
     &     r_cf, c_cf, a_el_cf, cs, n_wo)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'ICTVR.INC'
      integer r_cf, c_cf, a_el_cf, cs, n_wo
      r_cf = i_wo_z
      c_cf = i_wo_z
      a_el_cf =   0
      cs =   i_wo_z
      n_wo = r_cf*mx_n_r + c_cf*mx_n_c + cs
      return
      end

C->>> ----------------------------------------------> ems_rm_blk_dvx <<<
c     Removes the block for the vector required by Devex.
c
      subroutine ems_rm_blk_dvx(is)
      implicit none
      include 'EMSV.INC'
      include 'EMSMMGR.INC'
      include 'EMSPM.INC'
      include 'EMSMEM.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      integer is(0:is_n_en_m1)
      integer blk_n, p_ml_hdl
      integer mem_mgr_rt_cod

      if (iand(ml_blk_st_msk, ml_blk_st_dvx) .eq. 0) goto 8000
      p_ml_hdl = p_ml_bs_blk + (cu_ml_n-1)*ml_bs_blk_n_wo +
     &     ml_bs_blk_os_hdl
      blk_n = is(p_ml_hdl+ix_hdl_dvx_ix)
      call ems_mem_mgr_rm_blk(mem_mgr_rt_cod, is,
     &     blk_n)
      if (mem_mgr_rt_cod .ge. mem_mgr_rt_lvl_serious) then
         ems_msg_cod = ems_msg_lvl_serious
         goto 7100
      endif
      ml_blk_st_msk = ml_blk_st_msk - ml_blk_st_dvx
 7000 continue
 7100 continue
      return
 8000 continue
      if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9800)
      call ems_msg_wr_li(bug_msg_n)
      goto 7000
 9800 format('Model does not already has space for DEVEX')
      end

C->>> ------------------------------------------------> ems_u_dvx_ix <<<
c     Updates the index sets for Devex.
c
      subroutine ems_u_dvx_ix(dvx_ix)
      implicit none
      include 'EMSV.INC'
      include 'EMSPM.INC'
      include 'RSMICOM.INC'
      include 'ICTVR.INC'
      include 'EMSMSG.INC'
      integer dvx_ix(0:mx_n_c+n_r)

      if (vr_t_en_bs .le. 0 .or. vr_t_lv_bs .eq. 0 .or.
     &     vr_t_en_bs .eq. vr_t_lv_bs) then
         if (ems_msg_no_prt_fm .ge. 1) write(ems_li, 9900)
     &        vr_t_en_bs, vr_t_lv_bs
         call ems_msg_wr_li(warn_msg_n)
CM      IF (emsol_deb .EQ. 1) THEN
C?         call ems_dump
CM      ENDIF
         go to 7000
      endif
c
c     The variable to leave the basis is in R (so will leave H).
c
      if (dvx_ix(vr_t_lv_bs) .gt. 0) dvx_ix(vr_t_lv_bs) = 0
c
c     The variable to enter the basis is in R (and will enter H)
c
      if (dvx_ix(vr_t_en_bs) .eq. 0) dvx_ix(vr_t_en_bs) = 1
 7000 continue
      return
 9900 format('Calling u_dvx_ix with vr_t_en_bs, vr_t_lv_bs = ',
     &     i9, i9)
      end

max_stars_repo_name: ERGO-Code/EMSOL

max_issues_repo_name: ERGO-Code/EMSOL

max_forks_repo_name: ERGO-Code/EMSOL


text:
!-------------------------------------------------------------------------------
! Copyright (c) 2019 FrontISTR Commons
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> This module provides DOF based sparse matrix data structure (CSR and COO)
module m_sparse_matrix
  use hecmw_util
  use m_hecmw_comm_f
  implicit none

  private

  public :: SPARSE_MATRIX_TYPE_CSR
  public :: SPARSE_MATRIX_TYPE_COO

  public :: SPARSE_MATRIX_SYMTYPE_ASYM
  public :: SPARSE_MATRIX_SYMTYPE_SPD
  public :: SPARSE_MATRIX_SYMTYPE_SYM

  public :: sparse_matrix
  public :: sparse_matrix_set_type
  public :: sparse_matrix_init
  public :: sparse_matrix_finalize
  public :: sparse_matrix_dump
  public :: sparse_matrix_gather_rhs
  public :: sparse_matrix_scatter_rhs
  public :: sparse_matrix_is_sym

  integer(kind=kint), parameter :: SPARSE_MATRIX_TYPE_CSR=1
  integer(kind=kint), parameter :: SPARSE_MATRIX_TYPE_COO=2

  integer(kind=kint), parameter :: SPARSE_MATRIX_SYMTYPE_ASYM=0
  integer(kind=kint), parameter :: SPARSE_MATRIX_SYMTYPE_SPD=1
  integer(kind=kint), parameter :: SPARSE_MATRIX_SYMTYPE_SYM=2

  type sparse_matrix
    integer(kind=kint) :: type
    integer(kind=kint) :: symtype
    integer(kind=kint) :: N, N_loc
    integer(kind=kint) :: NZ
    integer(kind=kint), pointer :: IRN(:) => null()
    integer(kind=kint), pointer :: JCN(:) => null()
    real(kind=kreal), pointer :: A(:) => null()
    real(kind=kreal), pointer :: rhs(:)
    integer(kind=kint) :: offset
    integer(kind=kint), pointer :: N_COUNTS(:) => null()
    integer(kind=kint), pointer :: DISPLS(:) => null()
    integer(kind=kint), pointer :: conv_ext(:) => null()
    integer(kind=kint) :: iterlog
    integer(kind=kint) :: timelog
    logical :: is_initialized = .false.
  end type sparse_matrix

contains

  !!! public subroutines

  subroutine sparse_matrix_set_type(spMAT, type, symtype)
    type (sparse_matrix), intent(inout) :: spMAT
    integer(kind=kint), intent(in) :: type
    integer(kind=kint), intent(in) :: symtype
    if (type == SPARSE_MATRIX_TYPE_CSR .or. &
        type == SPARSE_MATRIX_TYPE_COO) then
      spMAT%type = type
    else
      write(*,*) 'ERROR: unknown sparse matrix type'
      stop
    endif
    if (symtype == SPARSE_MATRIX_SYMTYPE_ASYM .or. &
        symtype == SPARSE_MATRIX_SYMTYPE_SPD .or. &
        symtype == SPARSE_MATRIX_SYMTYPE_SYM) then
      spMAT%symtype = symtype
    else
      write(*,*) 'ERROR: unknown symmetry type for sparse matrix'
      stop
    endif
  end subroutine sparse_matrix_set_type

  subroutine sparse_matrix_init(spMAT, N_loc, NZ)
    type (sparse_matrix), intent(inout) :: spMAT
    integer(kind=kint), intent(in) :: N_loc
    integer(kind=kint), intent(in) :: NZ
    if (spMAT%is_initialized) then
      !write(*,*) "WARNING: sparse_matrix_init_prof: already initialized; freeing"
      call sparse_matrix_finalize(spMAT)
    endif
    call sparse_matrix_set_dimension(spMAT, N_loc)
    call sparse_matrix_set_offsets(spMAT)
    call sparse_matrix_allocate_arrays(spMAT, NZ)
    spMAT%is_initialized = .true.
  end subroutine sparse_matrix_init

  subroutine sparse_matrix_finalize(spMAT)
    type (sparse_matrix), intent(inout) :: spMAT
    call sparse_matrix_free_arrays(spMAT)
    call sparse_matrix_clear(spMAT)
    spMAT%is_initialized = .false.
  end subroutine sparse_matrix_finalize

  subroutine sparse_matrix_dump(spMAT)
    type(sparse_matrix), intent(in) :: spMAT
    integer(kind=kint),save :: num_call
    character(len=128) :: fname
    integer(kind=kint) :: i,myrank
    myrank=hecmw_comm_get_rank()
    write(fname,'(A,I0,A,I0)') 'sparse-matrix-',num_call,'.coo.',myrank
    num_call = num_call + 1
    write(*,*) 'Dumping sparse matrix to ', fname
    open(91,file=fname)
    if (spMAT%type == SPARSE_MATRIX_TYPE_CSR) then
      if (spMAT%symtype == SPARSE_MATRIX_SYMTYPE_ASYM) &
        write(91,*) '%SPARSE MATRIX CSR ASYM'
      if (spMAT%symtype == SPARSE_MATRIX_SYMTYPE_SPD) &
        write(91,*) '%SPARSE MATRIX CSR SPD'
      if (spMAT%symtype == SPARSE_MATRIX_SYMTYPE_SYM) &
        write(91,*) '%SPARSE MATRIX CSR SYM'
      write(91,*) spMAT%N_loc, spMAT%N_loc, spMAT%NZ
      do i=1,spMAT%N_loc+1
        write(91,*) spMAT%IRN(i)
      enddo
      do i=1,spMAT%NZ
        write(91,*) spMAT%JCN(i), spMAT%A(i)
      enddo
    elseif (spMAT%type == SPARSE_MATRIX_TYPE_COO) then
      if (spMAT%symtype == SPARSE_MATRIX_SYMTYPE_ASYM) &
        write(91,*) '%SPARSE MATRIX COO ASYM'
      if (spMAT%symtype == SPARSE_MATRIX_SYMTYPE_SPD) &
        write(91,*) '%SPARSE MATRIX COO SPD'
      if (spMAT%symtype == SPARSE_MATRIX_SYMTYPE_SYM) &
        write(91,*) '%SPARSE MATRIX COO SYM'
      write(91,*) spMAT%N_loc, spMAT%N_loc, spMAT%NZ
      do i=1,spMAT%NZ
        write(91,*) spMAT%IRN(i), spMAT%JCN(i), spMAT%A(i)
      enddo
    endif
    close(91)
  end subroutine sparse_matrix_dump

  subroutine sparse_matrix_gather_rhs(spMAT, rhs_all)
    type (sparse_matrix), intent(in) :: spMAT
    real(kind=kreal), intent(out) :: rhs_all(:)
    integer(kind=kint) :: ierr,i
    if (hecmw_comm_get_size() == 1) then
      do i=1,spMAT%N_loc
        rhs_all(i) = spMAT%rhs(i)
      enddo
    else
      call HECMW_GATHERV_REAL(spMAT%rhs, spMAT%N_loc, &
        rhs_all, spMAT%N_COUNTS, spMAT%DISPLS, &
        0, hecmw_comm_get_comm())
    endif
  end subroutine sparse_matrix_gather_rhs

  subroutine sparse_matrix_scatter_rhs(spMAT, rhs_all)
    type (sparse_matrix), intent(inout) :: spMAT
    real(kind=kreal), intent(in) :: rhs_all(:)
    integer(kind=kint) :: ierr,i
    if (hecmw_comm_get_size() == 1) then
      do i=1,spMAT%N_loc
        spMAT%rhs(i) = rhs_all(i)
      enddo
    else
      call HECMW_SCATTERV_REAL( &
        rhs_all, spMAT%N_COUNTS, spMAT%DISPLS, &
        spMAT%rhs, spMAT%N_loc, &
        0, hecmw_comm_get_comm())
    endif
  end subroutine sparse_matrix_scatter_rhs

  function sparse_matrix_is_sym(spMAT)
    type(sparse_matrix), intent(inout) :: spMAT
    logical :: sparse_matrix_is_sym
    sparse_matrix_is_sym = (spMAT%symtype > 0)
  end function sparse_matrix_is_sym

  !!! private subroutines

  subroutine sparse_matrix_set_dimension(spMAT, N_loc)
    type (sparse_matrix), intent(inout) :: spMAT
    integer(kind=kint), intent(in) :: N_loc
    integer(kind=kint) :: ierr
    spMAT%N_loc = N_loc
    if (hecmw_comm_get_size() == 1) then
      spMAT%N = spMAT%N_loc
    else
      call HECMW_ALLREDUCE_INT_1(spMAT%N_loc, spMAT%N, HECMW_SUM, &
        hecmw_comm_get_comm())
    endif
  end subroutine sparse_matrix_set_dimension

  subroutine sparse_matrix_set_offsets(spMAT)
    type (sparse_matrix), intent(inout) :: spMAT
    integer(kind=kint) :: i,ierr,nprocs,myrank
    nprocs=hecmw_comm_get_size()
    myrank=hecmw_comm_get_rank()
    ! OFFSET
    !if (myrank == 0) then
    if (associated(spMAT%N_COUNTS)) deallocate(spMAT%N_COUNTS)
    if (associated(spMAT%DISPLS)) deallocate(spMAT%DISPLS)
    allocate(spMAT%N_COUNTS(nprocs), spMAT%DISPLS(nprocs), stat=ierr)
    if (ierr /= 0) then
      write(*,*) " Allocation error, spMAT%N_COUNTS, spMAT%DISPLS"
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    !endif
    if (nprocs > 1) then
      call HECMW_ALLGATHER_INT_1(spMAT%N_loc, &
        spMAT%N_COUNTS, hecmw_comm_get_comm())
    endif
    spMAT%DISPLS(1)=0
    do i=1,nprocs-1
      spMAT%DISPLS(i+1)=spMAT%DISPLS(i)+spMAT%N_COUNTS(i)
    enddo
    spMAT%offset = spMAT%DISPLS(myrank+1)
  end subroutine sparse_matrix_set_offsets

  subroutine sparse_matrix_allocate_arrays(spMAT, NZ)
    type(sparse_matrix), intent(inout) :: spMAT
    integer(kind=kint), intent(in) :: NZ
    integer(kind=kint) :: N_loc
    integer(kind=kint) :: ierr
    if (associated(spMAT%IRN)) deallocate(spMAT%IRN)
    if (associated(spMAT%JCN)) deallocate(spMAT%JCN)
    if (associated(spMAT%A)) deallocate(spMAT%A)
    ierr = -1
    N_loc=spMAT%N_loc
    if (spMAT%type == SPARSE_MATRIX_TYPE_CSR) then
      allocate(spMAT%IRN(N_loc+1), spMAT%JCN(NZ), spMAT%A(NZ), stat=ierr)
    elseif (spMAT%type == SPARSE_MATRIX_TYPE_COO) then
      allocate(spMAT%IRN(NZ), spMAT%JCN(NZ), spMAT%A(NZ), stat=ierr)
    endif
    if (ierr /= 0) then
      write(*,*) " Allocation error, spMAT%IRN, spMAT%JCN, spMAT%A"
      call hecmw_abort(hecmw_comm_get_comm())
    endif
    spMAT%NZ = NZ
  end subroutine sparse_matrix_allocate_arrays

  subroutine sparse_matrix_free_arrays(spMAT)
    type (sparse_matrix), intent(inout) :: spMAT
    if (associated(spMAT%IRN)) deallocate(spMAT%IRN)
    if (associated(spMAT%JCN)) deallocate(spMAT%JCN)
    if (associated(spMAT%A)) deallocate(spMAT%A)
    if (associated(spMAT%N_COUNTS)) deallocate(spMAT%N_COUNTS)
    if (associated(spMAT%DISPLS)) deallocate(spMAT%DISPLS)
    if (associated(spMAT%conv_ext)) deallocate(spMAT%conv_ext)
  end subroutine sparse_matrix_free_arrays

  subroutine sparse_matrix_clear(spMAT)
    type(sparse_matrix), intent(inout) :: spMAT
    spMAT%type = 0
    spMAT%symtype = 0
    spMAT%N = 0
    spMAT%N_loc = 0
    spMAT%NZ = 0
    spMAT%IRN => null()
    spMAT%JCN => null()
    spMAT%A => null()
    spMAT%offset = 0
    spMAT%N_COUNTS => null()
    spMAT%DISPLS => null()
    spMAT%conv_ext => null()
    spMAT%iterlog = 0
    spMAT%timelog = 0
  end subroutine sparse_matrix_clear

end module m_sparse_matrix

max_stars_repo_name: masae-hayashi/FrontISTR

max_issues_repo_name: masae-hayashi/FrontISTR

max_forks_repo_name: masae-hayashi/FrontISTR


text:
!
! PRECISION BUG IN MUST_REGIONS
!
! first test case is tagged MAY
! bug second loop case is tagged EXACT,
! although it should be the same.
!
      program bugprec
      real a(10)
      integer i, j, n, p

      do i=1, 10
         a(i) = i
      enddo

      read *, n, p, i, j

!ABC if (1.le.p.and.(n.lt.0.or.n.gt.7)) stop
      if (1.le.p) then
         do i=n, n+3
            a(i) = 10
         enddo
      endif

!ABC if (1.le.p.and.(n.lt.0.or.n.gt.7)) stop
      do j=1, p
         do i=n, n+3
            a(i) = 10
         enddo
      enddo

      end


text:
      SUBROUTINE R_ASSERT
     &           (D_LO, D_VALUE, D_HI, VAL_STR, SUB_STR)
C ---------------------------------------------------------------------------
C
C ---------------------------------------------------------------------------
      IMPLICIT NONE
      INCLUDE 'knd_params.inc'
C
      REAL(KIND=RK) D_LO
      REAL(KIND=RK) D_VALUE
      REAL(KIND=RK) D_HI
      CHARACTER*(*) VAL_STR
      CHARACTER*(*) SUB_STR
C
      IF (D_LO .GT. D_VALUE .OR. D_VALUE .GT. D_HI) THEN
         WRITE (*, *) '*** ASSERTION FAILURE ***'
         WRITE (*, *) VAL_STR, ' is out of range ',
     &      D_LO, ' .. ', D_HI
         WRITE (*, *) VAL_STR, ' = ', D_VALUE
         WRITE (*, *) '*** ASSERTION FAILURE ***'
         CALL ERROR (SUB_STR, 'Value out of range')
      END IF
C
      RETURN
      END

max_stars_repo_name: meggermo/waterwaves

max_issues_repo_name: meggermo/waterwaves

max_forks_repo_name: meggermo/waterwaves


text:
Describe Users/kateboesser here.

I can look into it, keep harassing me on the wiki ♥  Users/StevenDaubert

Regarding your request for information on Unique Bikes Seen In Davis, Peter Wagner is who Id recommend locally. If youre in Alaska, any mechanicallyliterate person with access to scrap parts could build such a thing. SteveDavison