MODULE mixing_module

  use mixing_broyden_module
  use mixing_pulay_module
  use io_tools_module
  use rsdft_mpi_module
  use parallel_module, only: RSDFT_MPI_COMPLEX16
  use maxloc_err_module

  implicit none

  PRIVATE
  PUBLIC :: imix, beta &
           ,init_mixing,read_mixing,perform_mixing &
           ,finalize_mixing, restart_mixing
  PUBLIC :: calc_sqerr_mixing

  include 'mpif.h'

  integer :: imix = 10
  integer :: mmix = 4
  real(8) :: beta = 1.1d0
  integer :: iomix = 0
  integer :: iochk(2) = 0

  real(8),allocatable :: Xold(:,:,:)
  complex(8),allocatable :: Xin(:,:,:),Xou(:,:,:)
  real(8) :: beta0
  complex(8),parameter :: zero=(0.d0,0.d0)
  integer :: mmix_count=0
  real(8) :: diff(2),dif0(2),dif1(2)

  integer :: ML, ML0, MF0, MSP
  integer :: comm_grid, comm_spin
  real(8) :: dV
  logical :: disp_switch

  integer :: unit=2
  character(16),parameter :: file_name="vrho_mixing.dat1"
  integer,allocatable :: ir(:),id(:)
  integer :: myrank
  logical :: first_time=.true.

CONTAINS

  SUBROUTINE read_mixing
    implicit none
    logical :: disp
    call write_border( 0, " read_mixing(start)" )
    call IOTools_readIntegerKeyword(  "IMIX",  imix )
    call IOTools_readIntegerKeyword(  "MMIX",  mmix )
    call IOTools_readIntegerKeyword( "IOMIX", iomix )
    call IOTools_readIntegerKeyword( "IC", iochk(1) )
    call IOTools_readIntegerKeyword( "OC", iochk(2) )
    call IOTools_readReal8Keyword( "BETA", beta )
    call check_disp_switch( disp, 0 )
    if ( disp ) then
       if ( mmix < 1 ) then
          mmix=1
          write(*,*) "mmix is replaced to 1 : mmix=",mmix
       end if
       write(*,*) "beta       =",beta
       write(*,*) "iomix,IC,OC=",iomix,iochk(1:2)
    end if
    call write_border( 0, " read_mixing(end)" )
  END SUBROUTINE read_mixing

  SUBROUTINE init_mixing(ML0_in,MSP_in,nf1,nf2,comm_grid_in,comm_spin_in &
       ,dV_in,f,g,ir_in,id_in,myrank_in)
    implicit none
    integer,intent(IN) :: ML0_in,MSP_in,nf1,nf2,comm_grid_in,comm_spin_in
    integer,intent(IN) :: ir_in(0:),id_in(0:),myrank_in
    real(8),intent(IN) :: dV_in,f(ML0_in,nf1:nf2), g(ML0_in,nf1:nf2)
    integer :: m,ierr

    call write_border( 0, " init_mixing(start)" )

    myrank = myrank_in

    beta0      = 1.0d0-beta
    mmix_count = 0

    comm_grid = comm_grid_in
    comm_spin = comm_spin_in

    ML0 = ML0_in
    MSP = MSP_in
    MF0 = MSP*2
    call mpi_allreduce( ML0, ML, 1,MPI_INTEGER,MPI_SUM,comm_grid,ierr )

    dV = dV_in

    if ( .not.allocated(ir) ) then
       m=size(ir_in)
       allocate( ir(0:m-1) ) ; ir=0
       allocate( id(0:m-1) ) ; id=0
       ir(:)=ir_in(:)
       id(:)=id_in(:)
    end if

    if ( allocated(Xou)  ) deallocate( Xou )
    if ( allocated(Xin)  ) deallocate( Xin )
    if ( allocated(Xold) ) deallocate( Xold )
    allocate( Xold(ML0,MF0,mmix) ) ; Xold=0.0d0
    allocate(  Xin(ML0,MSP,mmix) ) ;  Xin=zero
    allocate(  Xou(ML0,MSP,mmix) ) ;  Xou=zero

    if ( iomix == 1 .and. iochk(1) == 3 ) then

       call restart_mixing

    else

       m=ML0*(nf2-nf1+1)
       call mpi_allgather( f(1,nf1),m,MPI_REAL8 &
            ,Xold(1,1,1),m,MPI_REAL8,comm_spin,ierr )
       call mpi_allgather( g(1,nf1),m,MPI_REAL8 &
            ,Xold(1,MSP+1,1),m,MPI_REAL8,comm_spin,ierr )

       if ( imix <= 29 ) then
          if ( mod(imix,2) == 0 ) then
             Xin(:,:,mmix) = Xold(:,1:MSP,1)
          else if ( mod(imix,2) == 1 ) then
             Xin(:,:,mmix) = Xold(:,MSP+1:2*MSP,1)
          end if
       else if ( imix >= 30 .and. imix <= 39 ) then
          if ( mod(imix,2) == 0 ) then
             call set_init_pulay( Xold(:,1:MSP,1), Xin(:,:,mmix) )
          else if ( mod(imix,2) == 1 ) then
             call set_init_pulay( Xold(:,MSP+1:2*MSP,1), Xin(:,:,mmix) )
          end if
       end if

       dif0(:) = 0.0d0

    end if

    call write_border( 0, " init_mixing(end)" )

  END SUBROUTINE init_mixing

  SUBROUTINE perform_mixing( m, n1, n2, f_io, g_io, disp_sw_in )
    implicit none
    integer,intent(IN)    :: m,n1,n2
    real(8),intent(INOUT) :: f_io(m,n1:n2),g_io(m,n1:n2)
    logical,optional,intent(IN) :: disp_sw_in
    real(8) :: err0(2),err(2),sum0(2),beta_bak,beta_min,dif_min(2)
    integer :: n,ierr,loop,max_loop,mmix_count_bak,i,i0
    complex(8),allocatable :: Xou_bak(:,:,:),Xin_bak(:,:,:)
    real(8),allocatable :: f(:,:),g(:,:),h(:,:),h_old(:,:),h_bak(:,:)

    call write_border( 1, " perform_mixing(start)" )

    disp_switch = .false.
    if ( present(disp_sw_in) ) disp_switch = disp_sw_in

!    if ( disp_switch ) write(*,'("----- mixing")')

    allocate( f(ML0,MSP) ) ; f=0.0d0
    allocate( g(ML0,MSP) ) ; g=0.0d0

    n=m*(n2-n1+1)
    call mpi_allgather( f_io(1,n1),n,MPI_REAL8,f,n,MPI_REAL8,comm_spin,ierr)
    call mpi_allgather( g_io(1,n1),n,MPI_REAL8,g,n,MPI_REAL8,comm_spin,ierr)

    allocate( h(ML0,MSP)     ) ; h=0.0d0
    allocate( h_old(ML0,MSP) ) ; h_old=0.0d0

    if ( mod(imix,2) == 0 ) then
       i0=0
       h_old(:,:) = Xold(:,i0+1:i0+MSP,1)
       h(:,:) = f(:,:)
    else if ( mod(imix,2) == 1 ) then
       i0=MSP
       h_old(:,:) = Xold(:,i0+1:i0+MSP,1)
       h(:,:) = g(:,:)
    end if

    allocate( Xin_bak(ML0,MSP,mmix) ) ; Xin_bak=zero
    allocate( Xou_bak(ML0,MSP,mmix) ) ; Xou_bak=zero
    allocate( h_bak(ML0,MSP)        ) ; h_bak=0.0d0

    Xin_bak(:,:,:) = Xin(:,:,:)
    Xou_bak(:,:,:) = Xou(:,:,:)
    h_bak(:,:)     = h(:,:)
    beta_bak       = beta
    mmix_count_bak = mmix_count
    beta_min       = beta
    dif_min(1:MSP) = 1.d10

    max_loop = 1

    if ( beta > 1.0d0 ) then
       if ( all(dif0==0.0d0) ) then
          max_loop = 1
          beta = 0.05d0
       else
          max_loop = 11
          beta = 1.0d0
       end if
    end if

    do loop=1,max_loop

       select case(imix)
       case default
          call simple_mixing( ML0, MSP, h_old, h )
       case(10:19)
          call pulay_mixing &
               ( ML0, MSP, comm_grid, mmix_count, mmix, beta, h, Xin, Xou )
!          call pulay_r2_mixing( ML0, MSP, h )
       case(20:29)
          call broyden_mixing &
               ( ML0, MSP, comm_grid, mmix_count, mmix, beta, h, Xin, Xou )
       case(30:39)
          call pulay_g_mixing &
               ( ML0, MSP, comm_grid, mmix_count, mmix, beta, h, Xin, Xou )
       end select

       sum0(:)=0.0d0
       do i=1,MSP
          sum0(i) = sum( (h(:,i)-Xold(:,i0+i,1))**2 )/ML
       end do
       call mpi_allreduce(sum0,dif1,MSP,MPI_REAL8,MPI_SUM,comm_grid,ierr)

       diff(:)=0.0d0
       do i=1,MSP
          if ( dif0(i) /= 0.0d0 ) diff(i) = dif1(i)/dif0(i)
       end do

!       if ( present(disp_sw_in) ) then
!          if ( disp_switch ) then
!             write(*,'(1x,"diff=",i4,2x,2(4g14.6,2x))') &
!                  loop,(diff(i),dif1(i),dif0(i),beta,i=1,MSP)
!          end if
!       end if

       if ( all( dif1(1:MSP) < dif_min(1:MSP) ) ) then
          beta_min = beta
          dif_min(1:MSP) = dif1(1:MSP)
       end if

       if ( loop < max_loop ) then

          if ( any( diff(1:MSP) > 0.5d0 ) ) then

             Xin(:,:,:) = Xin_bak(:,:,:)
             Xou(:,:,:) = Xou_bak(:,:,:)
             h(:,:)     = h_bak(:,:)
             mmix_count = mmix_count_bak
             beta       = beta*0.5d0
             beta0      = 1.0d0 - beta

          else

             exit

          end if

          if ( loop == max_loop-1 ) then
             beta  = beta_min
             beta0 = 1.0d0-beta
          end if

       end if

    end do ! loop

    beta        = beta_bak
    beta0       = 1.0d0 - beta
    dif0(1:MSP) = dif1(1:MSP)

    if ( mod(imix,2) == 0 ) then
       call chk_chrg( ML0,MSP,h )
       f(:,:) = h(:,:)
    else if ( mod(imix,2) == 1 ) then
       g(:,:) = h(:,:)
    end if

    f_io(:,n1:n2) = f(:,n1:n2)
    g_io(:,n1:n2) = g(:,n1:n2)

    Xold(:,1:MSP,1)       = f(:,:)
    Xold(:,MSP+1:2*MSP,1) = g(:,:)

    deallocate( h_bak )
    deallocate( Xou_bak )
    deallocate( Xin_bak )
    deallocate( h_old )
    deallocate( h,g,f )

    call write_border( 1, " perform_mixing(end)" )

  END SUBROUTINE perform_mixing

  SUBROUTINE calc_sqerr_mixing( m, n1, n2, f_io, g_io, sqerr_out )
    implicit none
    integer,intent(IN)    :: m,n1,n2
    real(8),intent(INOUT) :: f_io(m,n1:n2),g_io(m,n1:n2)
    real(8),intent(OUT)   :: sqerr_out(4)
    integer :: n,ierr
    real(8),allocatable :: f(:,:),g(:,:)
    call write_border( 1, " calc_sqerr_mixing(start)" )
    allocate( f(m,MSP) ) ; f=0.0d0
    allocate( g(m,MSP) ) ; g=0.0d0
    n=m*(n2-n1+1)
    call mpi_allgather( f_io(1,n1),n,MPI_REAL8,f,n,MPI_REAL8,comm_spin,ierr)
    call mpi_allgather( g_io(1,n1),n,MPI_REAL8,g,n,MPI_REAL8,comm_spin,ierr)
    call calc_sqerr( m, MSP, f, g, sqerr_out )
    deallocate( g,f )
    call write_border( 1, " calc_sqerr_mixing(end)" )
    return
  END SUBROUTINE calc_sqerr_mixing

  SUBROUTINE calc_sqerr(m,n,f,g,sqerr_out)
    implicit none
    integer,intent(IN)  :: m,n
    real(8),intent(IN)  :: f(m,n),g(m,n)
    real(8),intent(OUT) :: sqerr_out(4)
    real(8) :: err0(2*n),err(2*n)
    integer :: i,ierr,loc

    !flag_maxloc_err = .true.

    do i=1,n
       err0(i  )=sum( ( f(:,i)-Xold(:,i  ,1) )**2 )/ML
       err0(i+n)=sum( ( g(:,i)-Xold(:,i+n,1) )**2 )/ML
       call show_maxloc_err( f(:,i), Xold(:,i  ,1), 3 )
       call show_maxloc_err( g(:,i), Xold(:,i+n,1), 3 )
    end do
    call mpi_allreduce(err0,err,2*n,MPI_REAL8,MPI_SUM,comm_grid,ierr)

    sqerr_out(:)       = 0.0d0
    sqerr_out(1:n)     = err(1:n)
    sqerr_out(n+1:2*n) = err(n+1:2*n)

  END SUBROUTINE calc_sqerr

  SUBROUTINE chk_chrg( m,n,f )
    implicit none
    integer,intent(IN)    :: m,n
    real(8),intent(INOUT) :: f(m,n)
    real(8) :: chrg(2,n)
    integer :: i,ierr

    chrg=0.0d0
    do i=1,n
       chrg(1,i) = sum( f(:,i) )*dV
       chrg(2,i) = sum( f(:,i),mask=(f(:,i)<0.0d0) )*dV
    end do
    call rsdft_allreduce_sum( chrg, comm_grid )

    if ( any( chrg(2,1:n) < 0.0d0 ) ) then
!       if ( disp_switch ) then
!          write(*,'(1x,"NEGATIVE density exist (tot,neg)=",4g18.9)') &
!               ( chrg(1,i),chrg(2,i), i=1,n )
!       end if
    end if

!    where( f < 0.0d0 )
!       f = 0.0d0
!    end where
    f(:,:) = abs( f(:,:) )

  END SUBROUTINE chk_chrg

  SUBROUTINE simple_mixing(m,n,fold,f)
    implicit none
    integer,intent(IN)    :: m,n
    real(8),intent(IN)    :: fold(m,n)
    real(8),intent(INOUT) :: f(m,n)
    f(:,:) = beta0*fold(:,:) + beta*f(:,:)
  END SUBROUTINE simple_mixing

  SUBROUTINE pulay_r_mixing(m,n,f)
    implicit none
    integer,intent(IN)    :: m,n
    real(8),intent(INOUT) :: f(m,n)
    integer :: s,mmix0,ierr,i,i0,j0,mm,j
    integer,allocatable :: ipiv(:)
    real(8),allocatable :: rwork(:,:)
    complex(8) :: zc
    complex(8),allocatable :: A0(:,:),A1(:,:),b1(:),X(:),Y(:)

    Xou(:,:,mmix) = f(:,:)

    mmix_count = mmix_count + 1

    mmix0 = min( mmix_count, mmix )

    if ( mmix == 1 .or. mmix0 < mmix ) then
       do i=2,mmix
          Xin(:,:,i-1)=Xin(:,:,i)
          Xou(:,:,i-1)=Xou(:,:,i)
       end do
       Xin(:,:,mmix) = Xin(:,:,mmix) + beta*( Xou(:,:,mmix)-Xin(:,:,mmix) )
       f(:,:) = Xin(:,:,mmix)
       return
    end if

    allocate( ipiv(mmix0) )
    allocate( X(m),Y(m) )
    allocate( b1(mmix0) )
    allocate( A1(mmix0,mmix0) )
    allocate( A0(mmix0,mmix0) )

    b1(:)   = zero
    A1(:,:) = zero
    A0(:,:) = zero
    mm      = size(A1)

    do s=1,n

       do j0=1 ,mmix0
       do i0=j0,mmix0
          i=mmix-mmix0+i0
          j=mmix-mmix0+j0
          A0(i0,j0)=sum( conjg(Xou(:,s,i)-Xin(:,s,i)) &
                             *(Xou(:,s,j)-Xin(:,s,j)) )
          A0(j0,i0)=conjg( A0(i0,j0) )
       end do
       end do

       call mpi_allreduce(A0,A1,mm,RSDFT_MPI_COMPLEX16,mpi_sum,comm_grid,ierr)

       b1(1:mmix0) = (1.d0,0.d0)
       A0(:,:)     = A1(:,:)

       call zgesv(mmix0,1,A1,mmix0,ipiv,b1,mmix0,ierr)

       zc=1.d0/sum( b1(1:mmix0) )
       b1(1:mmix0)=zc*b1(1:mmix0)

       X(:)=zero
       Y(:)=zero

       do i0=1,mmix0
          i=mmix-mmix0+i0
          X(:)=X(:)+b1(i0)*Xin(:,s,i)
          Y(:)=Y(:)+b1(i0)*Xou(:,s,i)
       end do

       do i=max(1,mmix-mmix0),mmix-1
          Xin(:,s,i)=Xin(:,s,i+1)
          Xou(:,s,i)=Xou(:,s,i+1)
       end do

       Xin(:,s,mmix) = X(:) + beta*( Y(:)-X(:) )

    end do ! s

    f(:,:) = real( Xin(:,:,mmix) )

    deallocate( A0,A1,b1,Y,X,ipiv )
    return

  END SUBROUTINE pulay_r_mixing

  SUBROUTINE pulay_r2_mixing(m,n,f)
    implicit none
    integer,intent(IN)    :: m,n
    real(8),intent(INOUT) :: f(m,n)
    integer :: s,mmix0,ierr,i,i0,j0,mm,j
    integer,allocatable :: ipiv(:)
    real(8),allocatable :: rwork(:,:)
    complex(8) :: zc
    complex(8),allocatable :: A0(:,:),A1(:,:),b1(:),X(:),Y(:)

    call write_border( 1, " pulay_r2_mixing(start)" )

    Xou(:,:,mmix) = f(:,:)

    mmix_count = mmix_count + 1

    mmix0 = min( mmix_count, mmix )

    if ( mmix == 1 .or. mmix0 < mmix ) then
       do i=2,mmix
          Xin(:,:,i-1)=Xin(:,:,i)
          Xou(:,:,i-1)=Xou(:,:,i)
       end do
       Xin(:,:,mmix) = Xin(:,:,mmix) + beta*( Xou(:,:,mmix)-Xin(:,:,mmix) )
       f(:,:) = Xin(:,:,mmix)
       call write_border( 1, " pulay_r2_mixing(return)" )
       return
    end if

    allocate( ipiv(mmix0) )
    allocate( X(m),Y(m) )
    allocate( b1(mmix0) )
    allocate( A1(mmix0,mmix0) )
    allocate( A0(mmix0,mmix0) )

    b1(:)   = zero
    A1(:,:) = zero
    A0(:,:) = zero
    mm      = size(A1)

    do j0=1 ,mmix0
    do i0=j0,mmix0
       i=mmix-mmix0+i0
       j=mmix-mmix0+j0
       A0(i0,j0)=sum( conjg(Xou(:,:,i)-Xin(:,:,i)) &
                          *(Xou(:,:,j)-Xin(:,:,j)) )
       A0(j0,i0)=conjg( A0(i0,j0) )
    end do
    end do

    call mpi_allreduce(A0,A1,mm,RSDFT_MPI_COMPLEX16,mpi_sum,comm_grid,ierr)

    b1(1:mmix0) = (1.d0,0.d0)
    A0(:,:)     = A1(:,:)

    call zgesv(mmix0,1,A1,mmix0,ipiv,b1,mmix0,ierr)

    zc=1.d0/sum( b1(1:mmix0) )
    b1(1:mmix0)=zc*b1(1:mmix0)

    do s=1,n

       X(:)=zero
       Y(:)=zero

       do i0=1,mmix0
          i=mmix-mmix0+i0
          X(:)=X(:)+b1(i0)*Xin(:,s,i)
          Y(:)=Y(:)+b1(i0)*Xou(:,s,i)
       end do

       do i=max(1,mmix-mmix0),mmix-1
          Xin(:,s,i)=Xin(:,s,i+1)
          Xou(:,s,i)=Xou(:,s,i+1)
       end do

       Xin(:,s,mmix) = X(:) + beta*( Y(:)-X(:) )

    end do ! s

    f(:,:) = real( Xin(:,:,mmix) )

    deallocate( A0,A1,b1,Y,X,ipiv )

    call write_border( 1, " pulay_r2_mixing(end)" )

    return
  END SUBROUTINE pulay_r2_mixing

  SUBROUTINE finalize_mixing
    implicit none
    integer :: np,mr,n,m,s,ierr
    real(8) :: data_size
    real(8),allocatable :: d(:)
    complex(8),allocatable :: z(:)

    if ( iomix == 0 .or. iochk(2) /= 3 ) return

    data_size = 8.0d0*size(Xold) + 16.0d0*size(Xin) + 16.0d0*size(Xou)
    if ( myrank == 0 ) then
       write(*,'(a60," finalize_mmixg")') repeat("-",60)
       write(*,*) "data_size(MB)=",data_size/1024.d0**2
    end if
    if ( myrank == 0 ) then
       open( unit, file=file_name, form='unformatted' )
       write(unit) ML,MSP,MF0,mmix,mmix_count
       write(unit) beta,beta0,dif0
    end if
    allocate( d(ML) ) ; d=0.0d0
    do m=1,mmix
    do s=1,MF0
       call mpi_allgatherv( Xold(1,s,m),ML0,MPI_REAL8,d,ir,id, &
            MPI_REAL8,comm_grid,ierr )
       if ( myrank == 0 ) write(unit) d
    end do
    end do
    deallocate( d )
    allocate( z(ML) ) ; z=(0.0d0,0.0d0)
    do m=1,mmix
    do s=1,MSP
       call mpi_allgatherv( Xin(1,s,m),ML0,RSDFT_MPI_COMPLEX16,z,ir,id, &
            RSDFT_MPI_COMPLEX16,comm_grid,ierr )
       if ( myrank == 0 ) write(unit) z
    end do
    end do
    do m=1,mmix
    do s=1,MSP
       call mpi_allgatherv( Xou(1,s,m),ML0,RSDFT_MPI_COMPLEX16,z,ir,id, &
            RSDFT_MPI_COMPLEX16,comm_grid,ierr )
       if ( myrank == 0 ) write(unit) z
    end do
    end do
    deallocate( z )
    if ( myrank == 0 ) close( unit )
  END SUBROUTINE finalize_mixing

  SUBROUTINE restart_mixing
    implicit none
    integer :: m,s,ierr,itmp(5)
    real(8) :: dtmp(4)
    real(8),allocatable :: d(:)
    complex(8),allocatable :: z(:)
    if ( myrank == 0 ) then
       write(*,'(a60," restart_mmixg")') repeat("-",60)
    end if
    if ( myrank == 0 ) then
       open( unit, file=file_name, form='unformatted' )
       read(unit) itmp(1:5)
       read(unit) dtmp(1:4)
    end if
    call mpi_bcast(itmp,5,MPI_INTEGER,0,MPI_COMM_WORLD,ierr)
    call mpi_bcast(dtmp,4,MPI_REAL8,0,MPI_COMM_WORLD,ierr)
    ML=itmp(1)
    MSP=itmp(2)
    MF0=itmp(3)
    mmix=itmp(4)
    mmix_count=itmp(5)
    beta=dtmp(1)
    beta0=dtmp(2)
    dif0(1:2)=dtmp(3:4)
    allocate( d(ML) ) ; d=0.0d0
    do m=1,mmix
    do s=1,MF0
       if ( myrank == 0 ) read(unit) d
       call mpi_scatterv &
            (d,ir,id,MPI_REAL8,Xold(1,s,m),ML0,MPI_REAL8,0,comm_grid,ierr)
    end do
    end do
    deallocate( d )
    allocate( z(ML) ) ; z=(0.0d0,0.0d0)
    do m=1,mmix
    do s=1,MSP
       if ( myrank == 0 ) read(unit) z
       call mpi_scatterv(z,ir,id,RSDFT_MPI_COMPLEX16 &
                        ,Xin(1,s,m),ML0,RSDFT_MPI_COMPLEX16,0,comm_grid,ierr)
    end do
    end do
    do m=1,mmix
    do s=1,MSP
       if ( myrank == 0 ) read(unit) z
       call mpi_scatterv(z,ir,id,RSDFT_MPI_COMPLEX16 &
                        ,Xou(1,s,m),ML0,RSDFT_MPI_COMPLEX16,0,comm_grid,ierr)
    end do
    end do
    deallocate( z )
    if ( myrank == 0 ) close( unit )
  END SUBROUTINE restart_mixing

END MODULE mixing_module


! { dg-do compile }
! testcase from PR 17583
module bidon

 interface
  subroutine drivexc(nspden,rho_updn)
   integer,  intent(in) :: nspden
   integer, intent(in) :: rho_updn(nspden)
  end subroutine drivexc
 end interface

end module bidon

 subroutine nonlinear(nspden)

 use bidon

 integer,intent(in) :: nspden

 end subroutine nonlinear


text: The Buena Vista Hotel was widely considered to be the best hotels hotel in Davis, even though there were hog pens only 60 feet away and cow barns 200 feet away. The essence of Davis would have been particularly strong here, owing to the fact that manure was allowed to collect for months on end in the barnyards. The hotel burned down in 1918.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module type_harmonic_balance
#include<messenger.h>
    use messenger,                      only: write_line
    use mod_kinds,                      only: rk,ik
    use mod_constants,                  only: NO_FACE, NO_ID, ZERO, dQ_DIFF
    use mod_spatial,                    only: update_space

    use type_time_integrator_spectral,  only: time_integrator_spectral_t
    use type_system_assembler,          only: system_assembler_t

    use type_chidg_data,                only: chidg_data_t
    use type_nonlinear_solver,          only: nonlinear_solver_t
    use type_linear_solver,             only: linear_solver_t
    use type_preconditioner,            only: preconditioner_t
    use type_rvector,                   only: rvector_t
    use mod_HB_matrices,                only: calc_pseudo_spectral_operator
    use type_seed,                      only: seed_t
    use type_element_info,              only: element_info_t, element_info
    use type_chidg_vector
    use DNAD_D

    implicit none
    private

    !>
    !!
    !!  @author Mayank Sharma + Matteo Ugolotti
    !!  @date   2/13/2017
    !!
    !---------------------------------------------------------------------------------
    type, extends(time_integrator_spectral_t), public    :: harmonic_balance_t

    contains
        procedure   :: init
        procedure   :: step
    end type harmonic_balance_t
    !*********************************************************************************

    !>  Object for assembling the implicit system
    !!
    !!  @author Matteo Ugolotti + Mayank Sharma
    !!  @date 2/13/2017
    !!
    !!
    !---------------------------------------------------------------------------------
    type , extends(system_assembler_t), public :: assemble_harmonic_balance_t

    contains
        procedure   :: assemble
    end type assemble_harmonic_balance_t
    !*********************************************************************************

contains

    !>  Initialize the harmonic_balance_t time integrator
    !!
    !!  Main activity here, crate the assembler and attach it to the time_integrator
    !!  object so it can be passed to the nonlinear solver.
    !!
    !!  @author Matteo Ugolotti + Mayank Sharma
    !!  @date 2/13/2017
    !!
    !!  @author Mayank Sharma
    !!  @date   2/26/2017
    !!
    !---------------------------------------------------------------------------------
    subroutine init(self)
        class(harmonic_balance_t),  intent(inout)   :: self

        integer(ik)                         :: ierr
        type(assemble_harmonic_balance_t)   :: assemble_harmonic_balance

        call self%set_name('Harmonic Balance')

        if (allocated(self%system)) deallocate(self%system)
        allocate(self%system, source=assemble_harmonic_balance, stat=ierr)
        if (ierr /= 0) call AllocationError

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

    !>  Solving a steady-like equation set. Time levels are included in the LHS
    !!  by means of pseudo-spectral operator
    !!
    !!  Solving R(Q) = 0
    !!
    !!  @author Matteo Ugolotti + Mayank Sharma
    !!  @date 2/13/2017
    !!
    !!
    !---------------------------------------------------------------------------------
    subroutine step(self,data,nonlinear_solver,linear_solver,preconditioner)
        class(harmonic_balance_t),                      intent(inout)   :: self
        type(chidg_data_t),                             intent(inout)   :: data
        class(nonlinear_solver_t),   optional,          intent(inout)   :: nonlinear_solver
        class(linear_solver_t),      optional,          intent(inout)   :: linear_solver
        class(preconditioner_t),     optional,          intent(inout)   :: preconditioner

        ! Simply solve the nonlinear system. No iteration in time.
        call nonlinear_solver%solve(data,self%system,linear_solver,preconditioner)

        ! Store end residual from nonlinear solver.
        call self%residual_norm%push_back(nonlinear_solver%residual_norm%at(nonlinear_solver%residual_norm%size()))

    end subroutine step
    !*********************************************************************************

    !>  Assemble the system for the 'Harmonic Balance' equations with temporal contributions
    !!
    !!  @author Mayank Sharma + Matteo Ugolotti
    !!  @date   2/13/2017
    !!
    !---------------------------------------------------------------------------------
    subroutine assemble(self,data,differentiate,timing)
        class(assemble_harmonic_balance_t),   intent(inout)               :: self
        type(chidg_data_t),                   intent(inout)               :: data
        integer(ik),                          intent(in)                  :: differentiate
        real(rk),                             intent(inout), optional     :: timing

        integer(ik)             :: itime_outer, itime_inner, idom, ielem, ifield, ierr, ntime, &
                                   irow_start, irow_end, icol_start, icol_end
        real(rk),   allocatable :: temp_1(:), temp_2(:), D(:,:), temp_mat(:,:), hb_mat(:,:)
        type(chidg_vector_t)    :: rhs_tmp
        type(seed_t)            :: seed
        type(element_info_t)    :: elem_info

        associate ( rhs => data%sdata%rhs, lhs => data%sdata%lhs, q => data%sdata%q )
        call rhs%clear()
        if (differentiate == dQ_DIFF) call lhs%clear()

        ! Set local variables equal to the values set in time_manager
        ntime = size(data%time_manager%times)
        D = data%time_manager%D

        do itime_outer = 1,ntime

            ! Spatial update needed
            data%time_manager%itime = itime_outer
            data%time_manager%t     = data%time_manager%times(itime_outer)
            call update_space(data,differentiate,timing)

            do idom = 1,data%mesh%ndomains()

                if (allocated(temp_1) .and. allocated(temp_2)) deallocate(temp_1,temp_2)
                allocate(temp_1(data%mesh%domain(idom)%nterms_s),temp_2(data%mesh%domain(idom)%nterms_s), stat=ierr)
                if (ierr /= 0) call AllocationError

                do ielem = 1,data%mesh%domain(idom)%nelem
                    elem_info = data%mesh%get_element_info(idom,ielem)

                    do itime_inner = 1,ntime

                        ! LHS HB contribution for each variable
                        temp_mat = D(itime_outer,itime_inner)*data%mesh%domain(idom)%elems(ielem)%mass
                        if (allocated(hb_mat)) deallocate(hb_mat)
                        allocate(hb_mat(elem_info%nterms_s*elem_info%nfields,elem_info%nterms_s*elem_info%nfields), stat=ierr)
                        if (ierr /= 0) call AllocationError
                        hb_mat = ZERO

                        ! Accumulate variable contributions
                        do ifield = 1,data%eqnset(elem_info%eqn_ID)%prop%nprimary_fields()

                            ! Temporary variables for computing the temporal rhs contributions
                            temp_1 = D(itime_outer,itime_inner)*matmul(data%mesh%domain(idom)%elems(ielem)%mass, &
                                                            q%get_field(elem_info,ifield,itime_inner))

                            call rhs%add_field(temp_1,elem_info,ifield,itime_outer)

                            ! Accumulate hb contribution for the variable
                            irow_start = 1 + (ifield-1)*elem_info%nterms_s
                            irow_end   = irow_start + (elem_info%nterms_s-1)
                            icol_start = 1 + (ifield-1)*elem_info%nterms_s
                            icol_end   = icol_start + (elem_info%nterms_s-1)
                            hb_mat(irow_start:irow_end,icol_start:icol_end) = temp_mat

                        end do  ! ifield

                        ! LHS Contribution
                        call seed%init(elem_info%idomain_g,  &
                                       elem_info%idomain_l,  &
                                       elem_info%ielement_g, &
                                       elem_info%ielement_l, &
                                       elem_info%nfields,    &
                                       elem_info%nterms_s,   &
                                       elem_info%nterms_c,   &
                                       IRANK,                &
                                       itime_inner,          &
                                       elem_info%dof_start,  &
                                       elem_info%xdof_start, &
                                       NO_ID,                &
                                       NO_ID,                &
                                       NO_ID)

                        ! Store HB contribution for all fields to lhs at one time
                        if (itime_inner /= itime_outer) then
                            if (differentiate == dQ_DIFF) call lhs%store_element(hb_mat,elem_info,seed,itime_outer)
                        end if

                    end do  ! itime_inner
                end do  ! ielem
            end do  ! idom
        end do  ! itime_outer

        end associate

        ! Assemble RHS, LHS after HB contributions
        call data%sdata%rhs%assemble()
        if (differentiate == dQ_DIFF) call data%sdata%lhs%assemble()

    end subroutine assemble
    !*********************************************************************************

end module type_harmonic_balance

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
! ##################################################################################################################################
! Begin MIT license text.
! _______________________________________________________________________________________________________

! Copyright 2019 Dr William R Case, Jr (dbcase29@gmail.com)

! Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
! associated documentation files (the "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to
! the following conditions:

! The above copyright notice and this permission notice shall be included in all copies or substantial
! portions of the Software and documentation.

! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
! THE SOFTWARE.
! _______________________________________________________________________________________________________

! End MIT license text.

      SUBROUTINE WRITE_ELEM_ENGR_FORCE ( JSUB, NUM, IHDR )

! Writes blocks of element engineering force output for one element type, one subcase. Elements that can have engineering force
! output are the ones enumerated below fin the IF(TYPE == ???).

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  WRT_ERR, WRT_LOG, ANS, ERR, F04, F06
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, INT_SC_NUM, NDOFR, NUM_CB_DOFS, NVEC, SOL_NAME
      USE TIMDAT, ONLY                :  TSEC
      USE DEBUG_PARAMETERS, ONLY      :  DEBUG
      USE NONLINEAR_PARAMS, ONLY      :  LOAD_ISTEP
      USE LINK9_STUFF, ONLY           :  EID_OUT_ARRAY, OGEL
      USE SUBR_BEGEND_LEVELS, ONLY    :  WRITE_ELEM_ENGR_FORCE_BEGEND
      USE MODEL_STUF, ONLY            :  ELEM_ONAME, LABEL, SCNUM, STITLE, TITLE, TYPE
      USE WRITE_ELEM_ENGR_FORCE_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'WRITE_ELEM_ENGR_FORCE'
      CHARACTER(LEN=*), INTENT(IN)    :: IHDR              ! Indicator of whether to write an output header
      CHARACTER(128*BYTE)             :: FILL              ! Padding for output format
      CHARACTER(LEN=LEN(ELEM_ONAME))  :: ONAME             ! Element name to write out in F06 file

      INTEGER(LONG), INTENT(IN)       :: JSUB              ! Solution vector number
      INTEGER(LONG), INTENT(IN)       :: NUM               ! The number of rows of OGEL to write out
      INTEGER(LONG)                   :: BDY_COMP          ! Component (1-6) for a boundary DOF in CB analyses
      INTEGER(LONG)                   :: BDY_GRID          ! Grid for a boundary DOF in CB analyses
      INTEGER(LONG)                   :: BDY_DOF_NUM       ! DOF number for BDY_GRID/BDY_COMP
      INTEGER(LONG)                   :: I,J,J1            ! DO loop indices or counters
      INTEGER(LONG)                   :: NUM_TERMS         ! Number of terms to write out for shell elems
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = WRITE_ELEM_ENGR_FORCE_BEGEND

      REAL(DOUBLE)                    :: ABS_ANS(8)        ! Max ABS for all grids output for each of the 6 disp components
      REAL(DOUBLE)                    :: MAX_ANS(8)        ! Max for all grids output for each of the 6 disp components
      REAL(DOUBLE)                    :: MIN_ANS(8)        ! Min for all grids output for each of the 6 disp components

! **********************************************************************************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9001) SUBR_NAME,TSEC
 9001    FORMAT(1X,A,' BEGN ',F10.3)
      ENDIF

! **********************************************************************************************************************************
      FILL(1:) = ' '

! Get element output name

      ONAME(1:) = ' '
      CALL GET_ELEM_ONAME ( ONAME )

! Write output headers.

headr:IF (IHDR == 'Y') THEN

!--- Subcase num, TITLE, SUBT, LABEL:

         WRITE(F06,*)
         WRITE(F06,*)
         IF    ((SOL_NAME(1:7) == 'STATICS') .OR. (SOL_NAME(1:8) == 'NLSTATIC')) THEN

            WRITE(F06,101) SCNUM(JSUB)

         ELSE IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 1)) THEN

            WRITE(F06,101) SCNUM(JSUB)

         ELSE IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 2)) THEN

            WRITE(F06,102) JSUB

         ELSE IF (SOL_NAME(1:5) == 'MODES') THEN

            WRITE(F06,102) JSUB

         ELSE IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN   ! Write info on what CB DOF the output is for

            IF ((JSUB <= NDOFR) .OR. (JSUB >= NDOFR+NVEC)) THEN
               IF (JSUB <= NDOFR) THEN
                  BDY_DOF_NUM = JSUB
               ELSE
                  BDY_DOF_NUM = JSUB-(NDOFR+NVEC)
               ENDIF
               CALL GET_GRID_AND_COMP ( 'R ', BDY_DOF_NUM, BDY_GRID, BDY_COMP  )
            ENDIF

            IF       (JSUB <= NDOFR) THEN
               WRITE(F06,103) JSUB, NUM_CB_DOFS, 'acceleration', BDY_GRID, BDY_COMP
            ELSE IF ((JSUB > NDOFR) .AND. (JSUB <= NDOFR+NVEC)) THEN
               WRITE(F06,104) JSUB, NUM_CB_DOFS, JSUB-NDOFR
            ELSE
               WRITE(F06,103) JSUB, NUM_CB_DOFS, 'displacement', BDY_GRID, BDY_COMP
            ENDIF

         ENDIF

         IF (TITLE(INT_SC_NUM)(1:)  /= ' ') THEN
            WRITE(F06,201) TITLE(INT_SC_NUM)
         ENDIF

         IF (STITLE(INT_SC_NUM)(1:) /= ' ') THEN
            WRITE(F06,201) STITLE(INT_SC_NUM)
         ENDIF

         IF (LABEL(INT_SC_NUM)(1:)  /= ' ') THEN
            WRITE(F06,201) LABEL(INT_SC_NUM)
         ENDIF

         WRITE(F06,*)

!--- 1st 2 lines of element specific headers - general info on what type of output:

         IF      (TYPE(1:3) == 'BAR') THEN
            IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
               WRITE(F06,302) FILL(1:33)
            ELSE
               WRITE(F06,301) FILL(1:39)
            ENDIF
            WRITE(F06,401) FILL(1:45), ONAME

         ELSE IF (TYPE(1:4) == 'BUSH') THEN
            IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
               WRITE(F06,302) FILL(1:19)
            ELSE
               WRITE(F06,301) FILL(1:24)
            ENDIF
            WRITE(F06,401) FILL(1:29), ONAME

         ELSE IF (TYPE(1:4) == 'ELAS') THEN
            IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
               WRITE(F06,302) FILL(1:27)
            ELSE
               WRITE(F06,301) FILL(1:33)
            ENDIF
            WRITE(F06,401) FILL(1:37), ONAME

         ELSE IF (TYPE(1:3) == 'ROD') THEN
            IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
               WRITE(F06,302) FILL(1:27)
            ELSE
               WRITE(F06,301) FILL(1:33)
            ENDIF
            WRITE(F06,401) FILL(1:37), ONAME

         ELSE IF (TYPE(1:5) == 'SHEAR') THEN
            IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
               WRITE(F06,302) FILL(1:22)
            ELSE
               WRITE(F06,301) FILL(1:29)
            ENDIF
            WRITE(F06,401) FILL(1:32), ONAME

         ELSE IF((TYPE(1:5) == 'TRIA3') .OR. (TYPE(1:5) == 'QUAD4')) THEN
            IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
               WRITE(F06,302) FILL(1:33)
            ELSE
               WRITE(F06,301) FILL(1:39)
            ENDIF
            WRITE(F06,401) FILL(1:43), ONAME

         ENDIF

!--- Header lines describing columns of output for an element type:

         IF      (TYPE(1:3) == 'BAR'  ) THEN
            WRITE(F06,1101) FILL(1: 0), FILL(1: 0)

         ELSE IF (TYPE(1:4) == 'ELAS') THEN
            WRITE(F06,1201) FILL(1: 0), FILL(1: 0)

         ELSE IF (TYPE(1:3) == 'ROD') THEN
            WRITE(F06,1301) FILL(1: 0), FILL(1: 0)

         ELSE IF (TYPE(1:5) == 'SHEAR') THEN
            WRITE(F06,1401) FILL(1: 0), FILL(1: 0)

         ELSE IF ((TYPE(1:5) == 'TRIA3') .OR. (TYPE(1:5) == 'QUAD4')) THEN
            WRITE(F06,1501) FILL(1: 0), FILL(1: 0), FILL(1: 0)

         ELSE IF (TYPE(1:4) == 'BUSH') THEN
            WRITE(F06,1601) FILL(1: 0), FILL(1: 0)

         ENDIF

         IF (DEBUG(200) > 0) CALL WRITE_ANS ( 'Headers' )

      ENDIF headr

! Write element force output

      IF      (TYPE == 'BAR     ') THEN

         DO I=1,NUM
            WRITE(F06,1102) FILL(1: 0), EID_OUT_ARRAY(I,1),(OGEL(I,J),J=1,8)
         ENDDO
         CALL GET_MAX_MIN_ABS ( 1, 8 )
         WRITE(F06,1103) FILL(1: 0), FILL(1: 0), (MAX_ANS(J),J=1,8), FILL(1: 0), (MIN_ANS(J),J=1,8), FILL(1: 0),                   &
                                                 (ABS_ANS(J),J=1,8), FILL(1: 0)
         IF (DEBUG(200) > 0) CALL WRITE_ANS ( 'BAR' )

      ELSE IF (TYPE(1:4) == 'ELAS') THEN                   ! Engr force for ELAS was put into OGEL(I,1)

!xx      WRITE(F06,1202) FILL(1: 0), (EID_OUT_ARRAY(I,1),OGEL(I,1),I=1,NUM)
         J1 = 1
         DO I=1,NUM,5
            IF (J1+4 <= NUM) THEN
               WRITE(F06,1202) FILL(1: 0), (EID_OUT_ARRAY(J,1), OGEL(J,1), J=J1,J1+4)
               J1 = J1 + 5
            ELSE
               WRITE(F06,1202) FILL(1: 0), (EID_OUT_ARRAY(J,1), OGEL(J,1), J=J1,NUM)
            ENDIF
         ENDDO
         CALL GET_MAX_MIN_ABS ( 1, 1 )
         WRITE(F06,1203) FILL(1: 0), FILL(1: 0), (MAX_ANS(J),J=1,1), FILL(1: 0), (MIN_ANS(J),J=1,1), FILL(1: 0),                   &
                                                 (ABS_ANS(J),J=1,1), FILL(1: 0)
         IF (DEBUG(200) > 0) CALL WRITE_ANS ( 'ELAS' )

      ELSE IF (TYPE == 'ROD     ') THEN

         J1 = 1
         DO I=1,NUM,3
            IF (J1+2 <= NUM) THEN
               WRITE(F06,1302) FILL(1: 0), (EID_OUT_ARRAY(J,1), OGEL(J,7), OGEL(J,8), J=J1,J1+2)
               J1 = J1 + 3
            ELSE
               WRITE(F06,1302) FILL(1: 0), (EID_OUT_ARRAY(J,1), OGEL(J,7), OGEL(J,8), J=J1,NUM)
            ENDIF
         ENDDO
         CALL GET_MAX_MIN_ABS ( 7, 8 )
         WRITE(F06,1303) FILL(1: 0), FILL(1: 0), (MAX_ANS(J),J=7,8), FILL(1: 0), (MIN_ANS(J),J=7,8), FILL(1: 0),                   &
                                                 (ABS_ANS(J),J=7,8), FILL(1: 0)
         IF (DEBUG(200) > 0) CALL WRITE_ANS ( 'ROD' )

      ELSE IF (TYPE == 'SHEAR   ') THEN

         J1 = 1
         DO I=1,NUM,2
            IF      (J1+1 <= NUM) THEN
               WRITE(F06,1402) FILL(1: 0), (EID_OUT_ARRAY(J,1), OGEL(J,1), OGEL(J,2), OGEL(J,3), J=J1,J1+1)
               J1 = J1 + 2
            ELSE
               WRITE(F06,1402) FILL(1: 0), (EID_OUT_ARRAY(J,1), OGEL(J,1), OGEL(J,2), OGEL(J,3), J=J1,NUM)
            ENDIF
         ENDDO
         CALL GET_MAX_MIN_ABS ( 1, 3 )
         WRITE(F06,1403) FILL(1: 0), FILL(1: 0), (MAX_ANS(J),J=1,3), FILL(1: 0), (MIN_ANS(J),J=1,3), FILL(1: 0),                   &
                                                 (ABS_ANS(J),J=1,3), FILL(1: 0)
         IF (DEBUG(200) > 0) CALL WRITE_ANS ( 'SHEAR' )

      ELSE IF ((TYPE == 'TRIA3K  ') .OR. (TYPE == 'QUAD4K  ')) THEN

         DO I=1,NUM
            WRITE(F06,1512) FILL(1: 0), EID_OUT_ARRAY(I,1),(OGEL(I,J),J=1,6)
         ENDDO
         CALL GET_MAX_MIN_ABS ( 1, 8 )
         WRITE(F06,1513) FILL(1: 0), FILL(1: 0), (MAX_ANS(J),J=1,6), FILL(1: 0), (MIN_ANS(J),J=1,6), FILL(1: 0),                   &
                                                 (ABS_ANS(J),J=1,6), FILL(1: 0)
         NUM_TERMS = 6
         IF (DEBUG(200) > 0) CALL WRITE_ANS ( 'SHELL' )

      ELSE IF ((TYPE == 'TRIA3   ') .OR. (TYPE == 'QUAD4   ')) THEN

         DO I=1,NUM
            WRITE(F06,1522) FILL(1: 0), EID_OUT_ARRAY(I,1),(OGEL(I,J),J=1,8)
         ENDDO
         CALL GET_MAX_MIN_ABS ( 1, 8 )
         WRITE(F06,1523) FILL(1: 0), FILL(1: 0), (MAX_ANS(J),J=1,8), FILL(1: 0), (MIN_ANS(J),J=1,8), FILL(1: 0),                   &
                                                 (ABS_ANS(J),J=1,8), FILL(1: 0)
         NUM_TERMS = 8
         IF (DEBUG(200) > 0) CALL WRITE_ANS ( 'SHELL' )

      ELSE IF (TYPE(1:4) == 'BUSH') THEN                   ! Engr force for BUSH was put into OGEL(I,1-6)

         DO I=1,NUM
            WRITE(F06,1602) FILL(1: 0), EID_OUT_ARRAY(I,1),(OGEL(I,J),J=1,6)
         ENDDO
         CALL GET_MAX_MIN_ABS ( 1, 6 )
         WRITE(F06,1603) FILL(1: 0), FILL(1: 0), (MAX_ANS(J),J=1,6), FILL(1: 0), (MIN_ANS(J),J=1,6), FILL(1: 0),                   &
                                                 (ABS_ANS(J),J=1,6), FILL(1: 0)

         IF (DEBUG(200) > 0) CALL WRITE_ANS ( 'BUSH' )

      ENDIF

! **********************************************************************
      IF (WRT_LOG >= SUBR_BEGEND) THEN
         CALL OURTIM
         WRITE(F04,9002) SUBR_NAME,TSEC
 9002    FORMAT(1X,A,' END  ',F10.3)
      ENDIF

      RETURN

! **********************************************************************************************************************************
  101 FORMAT(' OUTPUT FOR SUBCASE ',I8)

  102 FORMAT(' OUTPUT FOR EIGENVECTOR ',I8)

  103 FORMAT(' OUTPUT FOR CRAIG-BAMPTON DOF ',I8,' OF ',I8,' (boundary ',A,' for grid',I8,' component',I2,')')

  104 FORMAT(' OUTPUT FOR CRAIG-BAMPTON DOF ',I8,' OF ',I8,' (modal acceleration for mode ',I8,')')

  201 FORMAT(1X,A)

  301 FORMAT(16X,A,'E L E M E N T   E N G I N E E R I N G   F O R C E S')

  302 FORMAT(16X,A,'C B   E L E M E N T   E N G I N E E R I N G   F O R C E   O T M')

  401 FORMAT(16X,A,'F O R   E L E M E N T   T Y P E   ',A11)

! BAR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1101 FORMAT(16X,A,' Element       Bend-Moment End A           Bend-Moment End B              - Shear -              Axial'        &
          ,'         Torque'  &
          ,/,16X,A,'    ID       Plane 1       Plane 2       Plane 1       Plane 2      Plane 1       Plane 2        Force')

 1102 FORMAT(16X,A,I8,8(1ES14.6))

 1103 FORMAT(1X,A,'         ------------- ------------- ------------- ------------- ------------- ------------- -------------',    &
                        ' -------------',/,                                                                                &
             16X,A,'MAX* :  ',8(ES14.6),/,                                                                                &
             16X,A,'MIN* :  ',8(ES14.6),//,                                                                                &
             16X,A,'ABS* :  ',8(ES14.6),/,                                                                                &
             16X,A,'*for output set')

! ELAS >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1201 FORMAT(16X,A,' Element     Force      Element     Force      Element     Force      Element     Force      Element     Force'&
          ,/,16X,A,'    ID                     ID                     ID                     ID                     ID')

 1202 FORMAT(16X,A,5(I8,1ES14.6))

 1203 FORMAT(16X,A,'         -------------',/,                                                                                &
             16X,A,'MAX* :  ',1(ES14.6),/,                                                                                &
             16X,A,'MIN* :  ',1(ES14.6),//,                                                                                &
             16X,A,'ABS* :  ',1(ES14.6),/,                                                                                &
             16X,A,'*for output set')

! ROD >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1301 FORMAT(16X,A,' Element     Axial        Torque      Element     Axial        Torque      Element     Axial        Torque'    &
          ,/,16X,A,'    ID       Force                       ID       Force                       ID       Force')

 1302 FORMAT(16X,A,3(I8,1ES14.6,1ES14.6))

 1303 FORMAT(16X,A,'         ------------- -------------',/,                                                                       &
             16X,A,'MAX* :  ',2(1ES14.6),/,                                                                                &
             16X,A,'MIN* :  ',2(1ES14.6),//,                                                                                &
             16X,A,'ABS* :  ',2(1ES14.6),/,                                                                                &
             16X,A,'*for output set')

! SHEAR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1401 FORMAT(16X,A,' Element        N o r m a l   F o r c e s           Element        N o r m a l   F o r c e s          '        &
          ,/,16X,A,'    ID        Nxx           Nyy           Nxy          ID        Nxx           Nyy           Nxy')

 1402 FORMAT(1X,A,2(I8,3(1ES14.6),1X))

 1403 FORMAT(16X,A,'         ------------- ------------- -------------',/,                                                         &
             16X,A,'MAX* :  ',3ES14.6,/,                                                                                &
             16X,A,'MIN* :  ',3ES14.6,//,                                                                                &
             16X,A,'ABS* :  ',3ES14.6,/,                                                                                &
             16X,A,'*for output set')

! SHELL >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1501 FORMAT(16X,A,' Element         N o r m a l   F o r c e s                       M o m e n t s'                                &
          ,19X,'T r a n s v e r s e',/16X,A,'    ID', 89X,'S h e a r   F o r c e s'                                                &
          ,/,16X,A,'              Nxx           Nyy           Nxy           Mxx           Myy           Mxy            Qx         '&
          ,'  Qy')

!            WRITE(F06,1501) FILL(1: 0), FILL(1: 0), FILL(1: 0)
 1512 FORMAT(16X,A,I8,6(1ES14.6))

 1513 FORMAT(16X,A,'          ------------- ------------- ------------- ------------- ------------- -------------',/,              &
             16X,A,'MAX* :  ',6(ES14.6),/,                                                                                &
             16X,A,'MIN* :  ',6(ES14.6),//,                                                                                &
             16X,A,'ABS* :  ',6(ES14.6),/,                                                                                &
             16X,A,'*for output set')

 1522 FORMAT(16X,A,I8,8(1ES14.6))

 1523 FORMAT(16X,A,'          ------------- ------------- ------------- ------------- ------------- ------------- -------------',  &
                        ' -------------',/,                                                                                &
             16X,A,'MAX* :  ',8(ES14.6),/,                                                                                &
             16X,A,'MIN* :  ',8(ES14.6),//,                                                                                &
             16X,A,'ABS* :  ',8(ES14.6),/,                                                                                &
             16X,A,'*for output set')

! BUSH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1601 FORMAT(16X,A,' Element      Force         Force         Force        Moment        Moment        Moment'                     &
          ,/,16X,A,'    ID         XE            YE            ZE            XE            YE            ZE')

 1602 FORMAT(16X,A,I8,6(1ES14.6))

 1603 FORMAT(16X,A,'          ------------- ------------- ------------- ------------- ------------- ------------- ',/,             &
             16X,A,'MAX* :  ',6(ES14.6),/,                                                                                &
             16X,A,'MIN* :  ',6(ES14.6),//,                                                                                &
             16X,A,'ABS* :  ',6(ES14.6),/,                                                                                &
             16X,A,'*for output set')

! **********************************************************************************************************************************

      CONTAINS

! ##################################################################################################################################

      SUBROUTINE WRITE_ANS ( WHICH )

      USE PENTIUM_II_KIND, ONLY       :  BYTE

      IMPLICIT NONE

      CHARACTER(LEN=*), INTENT(IN)    :: WHICH             ! Which of the below to write during this call

      INTEGER(LONG)                   :: II,JJ             ! DO loop indices or counters

! **********************************************************************************************************************************
      IF (WHICH == 'Headers') THEN

         IF (DEBUG(200) > 0) THEN
            WRITE(ANS,*)
            WRITE(ANS,*)
            IF    ((SOL_NAME(1:7) == 'STATICS') .OR. (SOL_NAME(1:8) == 'NLSTATIC')) THEN

               WRITE(ANS,101) SCNUM(JSUB)

         ELSE IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 1)) THEN

            WRITE(F06,101) SCNUM(JSUB)

         ELSE IF ((SOL_NAME(1:8) == 'BUCKLING') .AND. (LOAD_ISTEP == 2)) THEN

            WRITE(F06,102) JSUB

            ELSE IF (SOL_NAME(1:5) == 'MODES') THEN

               WRITE(ANS,102) JSUB

            ELSE IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN   ! Write info on what CB DOF the output is for

               IF ((JSUB <= NDOFR) .OR. (JSUB >= NDOFR+NVEC)) THEN
                  IF (JSUB <= NDOFR) THEN
                     BDY_DOF_NUM = JSUB
                  ELSE
                     BDY_DOF_NUM = JSUB-(NDOFR+NVEC)
                  ENDIF
                  CALL GET_GRID_AND_COMP ( 'R ', BDY_DOF_NUM, BDY_GRID, BDY_COMP  )
               ENDIF

               IF       (JSUB <= NDOFR) THEN
                  WRITE(ANS,103) JSUB, NUM_CB_DOFS, 'acceleration', BDY_GRID, BDY_COMP
               ELSE IF ((JSUB > NDOFR) .AND. (JSUB <= NDOFR+NVEC)) THEN
                  WRITE(ANS,104) JSUB, NUM_CB_DOFS, JSUB-NDOFR
               ELSE
                  WRITE(ANS,103) JSUB, NUM_CB_DOFS, 'displacement', BDY_GRID, BDY_COMP
               ENDIF

            ENDIF

            WRITE(ANS,*)

            IF ((TYPE(1:4) /= 'HEXA') .AND. (TYPE(1:5) /= 'PENTA') .AND. (TYPE(1:5) /= 'TETRA')) THEN
            IF (SOL_NAME(1:12) == 'GEN CB MODEL') THEN
               WRITE(ANS,302) FILL(1:16)
            ELSE
               WRITE(ANS,301) FILL(1:16)
            ENDIF
               WRITE(ANS,401) FILL(1:16), ONAME
            ENDIF

            IF      (TYPE(1:4) == 'ELAS') THEN
               WRITE(ANS,1111) FILL(1:16), FILL(1:16)

            ELSE IF (TYPE == 'BAR     ') THEN
               WRITE(ANS,1211) FILL(1:16), FILL(1:16)

            ELSE IF (TYPE == 'ROD     ') THEN
               WRITE(ANS,1311) FILL(1:16), FILL(1:16)

            ELSE IF (TYPE == 'SHEAR   ') THEN
               WRITE(ANS,1411) FILL(1:16), FILL(1:16)

            ELSE IF ((TYPE(1:5)== 'TRIA3') .OR. (TYPE(1:5) == 'QUAD4')) THEN
               WRITE(ANS,1511) FILL(1:16), FILL(1:16), FILL(1:16)

            ELSE IF (TYPE == 'BUSH    ') THEN
               WRITE(ANS,1611) FILL(1:16), FILL(1:16)

            ENDIF

         ENDIF

      ELSE IF (WHICH == 'ELAS' ) THEN
         WRITE(ANS,1112) FILL(1:16), (EID_OUT_ARRAY(II,1),OGEL(II,1),II=1,NUM)
         WRITE(ANS,1113) (MAX_ANS(JJ),JJ=1,1),(MIN_ANS(JJ),JJ=1,1),(ABS_ANS(JJ),JJ=1,1)

      ELSE IF (WHICH == 'BAR'  ) THEN
         DO II=1,NUM
            WRITE(ANS,1212) FILL(1:16), EID_OUT_ARRAY(II,1),(OGEL(II,JJ),JJ=1,8)
         ENDDO
         WRITE(ANS,1213) (MAX_ANS(JJ),JJ=1,8),(MIN_ANS(JJ),JJ=1,8),(ABS_ANS(JJ),JJ=1,8)

      ELSE IF (WHICH == 'ROD'  ) THEN
         WRITE(ANS,1312) (FILL(1:16), EID_OUT_ARRAY(II,1), OGEL(II,7), OGEL(II,8),II=1,NUM)
         WRITE(ANS,1313) (MAX_ANS(JJ),JJ=7,8),(MIN_ANS(JJ),JJ=7,8),(ABS_ANS(JJ),JJ=7,8)

      ELSE IF (WHICH == 'SHEAR') THEN
         DO II=1,NUM
            WRITE(ANS,1412) FILL(1:16), EID_OUT_ARRAY(II,1),(OGEL(II,JJ),JJ=1,3)
         ENDDO
         WRITE(ANS,1413) (MAX_ANS(JJ),JJ=1,8),(MIN_ANS(JJ),JJ=1,8),(ABS_ANS(JJ),JJ=1,3)

      ELSE IF (WHICH == 'SHELL') THEN
         DO II=1,NUM
            WRITE(ANS,1512) FILL(1:16), EID_OUT_ARRAY(II,1),(OGEL(II,JJ),JJ=1,NUM_TERMS)
         ENDDO
         WRITE(ANS,1513) (MAX_ANS(JJ),JJ=1,8),(MIN_ANS(JJ),JJ=1,8),(ABS_ANS(JJ),JJ=1,NUM_TERMS)

      ELSE IF (WHICH == 'BUSH' ) THEN
         DO II=1,NUM
            WRITE(ANS,1612) FILL(1:16), EID_OUT_ARRAY(II,1),(OGEL(II,JJ),JJ=1,6)
         ENDDO
         WRITE(ANS,1613) (MAX_ANS(JJ),JJ=1,6),(MIN_ANS(JJ),JJ=1,6),(ABS_ANS(JJ),JJ=1,6)

      ENDIF

      RETURN

! **********************************************************************************************************************************
  101 FORMAT(' OUTPUT FOR SUBCASE ',I8)

  102 FORMAT(' OUTPUT FOR EIGENVECTOR ',I8)

  103 FORMAT(' OUTPUT FOR CRAIG-BAMPTON DOF ',I8,' OF ',I8,' (boundary ',A,' for grid',I8,' component',I2,')')

  104 FORMAT(' OUTPUT FOR CRAIG-BAMPTON DOF ',I8,' OF ',I8,' (modal acceleration for mode ',I8,')')

  301 FORMAT(39X,A,'E L E M E N T   E N G I N E E R I N G   F O R C E S')

  302 FORMAT(33X,A,'C B   E L E M E N T   E N G I N E E R I N G   F O R C E   O T M')

  401 FORMAT(44X,A,'F O R   E L E M E N T   T Y P E   ',A11)

! ELAS >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1111 FORMAT(2X,A,'Element     Force'                                                                                &
          ,/,2X,A,'   ID')

 1112 FORMAT(A,5(I8,1ES14.6))

 1113 FORMAT(11X,'              -------------',/,                                                                                &
             1X,'MAX (for output set):  ',1(ES14.6),/,                                                                             &
             1X,'MIN (for output set):  ',1(ES14.6),//,                                                                            &
             1X,'ABS (for output set):  ',1(ES14.6))

! BAR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1211 FORMAT(2X,A,'Element       Bend-Moment End A           Bend-Moment End B              - Shear -              Axial'          &
          ,'         Torque'  &
          ,/,2X,A,'   ID       Plane 1       Plane 2       Plane 1       Plane 2      Plane 1       Plane 2        Force')

 1212 FORMAT(A,I8,8(1ES14.6))

 1213 FORMAT(11X,'              ------------- ------------- ------------- ------------- ------------- ------------- -------------',&
                             ' -------------',/,                                                                                &
             1X,'MAX (for output set):  ',8(ES14.6),/,                                                                             &
             1X,'MIN (for output set):  ',8(ES14.6),//,                                                                            &
             1X,'ABS (for output set):  ',8(ES14.6))

! ROD >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1311 FORMAT(2X,A,'Element     Axial        Torque'                                                                                &
          ,/,2X,A,'   ID       Force')

 1312 FORMAT(A,I8,1ES14.6,1ES14.6)

 1313 FORMAT(11X,'              ------------- -------------',/,                                                                    &
             1X,'MAX (for output set):  ',2(ES14.6),/,                                                                             &
             1X,'MIN (for output set):  ',2(ES14.6),//,                                                                            &
             1X,'ABS (for output set):  ',2(ES14.6))

! SHEAR >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1411 FORMAT(1X,A,'                                        Element        N o r m a l   F o r c e s'                               &
          ,/,1X,A,'                                           ID        Nxx           Nyy           Nxy')

 1412 FORMAT(A,40X,I8,3(1ES14.6))

 1413 FORMAT(1X,A,'                                                ------------- ------------- -------------',/,                   &
             1X,A,'                                        MAX* : ',3(ES14.6),/,                                                   &
             1X,A,'                                        MIN* : ',3(ES14.6),//,                                                  &
             1X,A,'                                        ABS* : ',3(ES14.6),/,                                                   &
             1X,A,'                                        *for output set')

! SHELL >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1511 FORMAT(2X,A,'Element        N o r m a l   F o r c e s                       M o m e n t s'                                   &
          ,19X,'T r a n s v e r s e',/,A,95x,'S h e a r   F o r c e s'                                                             &
          ,/,2X,A,'   ID       Nxx           Nyy           Nxy           Mxx           Myy           Mxy            Qx          '  &
          ,'  Qy')

 1512 FORMAT(A,I8,8(1ES14.6))

 1513 FORMAT(11X,'              ------------- ------------- ------------- ------------- ------------- ------------- -------------',&
                             ' -------------',/,                                                                                &
             1X,'MAX (for output set):  ',8(ES14.6),/,                                                                             &
             1X,'MIN (for output set):  ',8(ES14.6),//,                                                                            &
             1X,'ABS (for output set):  ',8(ES14.6))

! BUSH >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
 1611 FORMAT(1X,A,' Element      Force         Force         Force         Force         Force         Force' &
          ,/,1X,A,'    ID         TX            TY            TZ            RX            RY            RZ')

 1612 FORMAT(A,I8,6(1ES14.6))

 1613 FORMAT(1X,'                        ------------- ------------- ------------- ------------- ------------- ------------- ',/,  &
             1X,'MAX (for output set):  ',6(ES14.6),/,                                                                             &
             1X,'MIN (for output set):  ',6(ES14.6),//,                                                                            &
             1X,'ABS (for output set):  ',6(ES14.6))

! **********************************************************************************************************************************

      END SUBROUTINE WRITE_ANS

! ##################################################################################################################################

      SUBROUTINE GET_MAX_MIN_ABS ( BEG_COL, END_COL )

      USE PENTIUM_II_KIND, ONLY       :  LONG
      USE CONSTANTS_1, ONLY           :  ZERO
      USE MACHINE_PARAMS, ONLY        :  MACH_LARGE_NUM

      IMPLICIT NONE

      INTEGER(LONG), INTENT(IN)       :: BEG_COL           ! Col number in OGEL where to beg for averaging to get max, min, abs
      INTEGER(LONG), INTENT(IN)       :: END_COL           ! Col number in OGEL where to end for averaging to get max, min, abs
      INTEGER(LONG)                   :: II,JJ             ! DO loop indices or counters

! **********************************************************************************************************************************
! Get MAX, MIN, ABS values

      DO JJ=BEG_COL,END_COL
         MAX_ANS(JJ) = -MACH_LARGE_NUM
      ENDDO

      DO II=1,NUM
         DO JJ=BEG_COL,END_COL
            IF (OGEL(II,JJ) > MAX_ANS(JJ)) THEN
               MAX_ANS(JJ) = OGEL(II,JJ)
            ENDIF
         ENDDO
      ENDDO

      DO JJ=BEG_COL,END_COL
         MIN_ANS(JJ) = MAX_ANS(JJ)
      ENDDO

      DO II=1,NUM
         DO JJ=BEG_COL,END_COL
            IF (OGEL(II,JJ) < MIN_ANS(JJ)) THEN
               MIN_ANS(JJ) = OGEL(II,JJ)
            ENDIF
         ENDDO
      ENDDO

      DO II=BEG_COL,END_COL
         ABS_ANS(II) = MAX( DABS(MAX_ANS(II)), DABS(MIN_ANS(II)) )
      ENDDO

      END SUBROUTINE GET_MAX_MIN_ABS

      END SUBROUTINE WRITE_ELEM_ENGR_FORCE

max_stars_repo_name: ceanwang/MYSTRAN

max_issues_repo_name: ceanwang/MYSTRAN

max_forks_repo_name: ceanwang/MYSTRAN


text:
program P17626
    implicit none
    integer :: n, i, j, sq
    integer, dimension(0:50000) :: save
    read *, n
    save(0) = 0
    do i = 1, n
        save(i) = i
        j = 1
        sq = 1
        do while (sq <= i)
            save(i) = min(save(i), save(i - sq) + 1)
            j = j + 1
            sq = j * j
        end do
    end do
    write (*, '(I0)') save(n)
end program P17626

max_stars_repo_name: daily-boj/kiwiyou

max_issues_repo_name: daily-boj/kiwiyou

max_forks_repo_name: daily-boj/kiwiyou


text:
+incdir+${TVIP_AXI_HOME}/sample/env
${TVIP_AXI_HOME}/sample/env/tvip_axi_sample_pkg.sv

max_stars_repo_name: azarmadr/tvip-axi

max_issues_repo_name: azarmadr/tvip-axi

max_forks_repo_name: azarmadr/tvip-axi