text:
program matrix_main
  use matrix_mod
  implicit none
! create one new object image
 type(matrix),pointer     :: m
 type(pgmImage),pointer   :: image

 !// Holding the number of command line arguments
  INTEGER                                          :: argc
  character(len=255)                               :: filename
  character(len=255)                               :: typeOfImage

  argc = COMMAND_ARGUMENT_COUNT()

  if (argc == 2) THEN
    call get_command_argument(1,filename)
    call get_command_argument(2,typeOfImage)

    SELECT CASE(typeOfImage)
     CASE('PGM')
      allocate(image)
      image = pgmImage()
      call image%load(filename) ! not need to pass 'this' argument
      call image%dump()
      deallocate(image)
     CASE('MATRIX')
      allocate(m)
      m = matrix()
      call m%load(filename)
      call m%dump()
      deallocate(m)

    END SELECT
  else
    write(*,*) 'A filename is needed!'
  endif

end program matrix_main

max_stars_repo_name: annefou/Fortran


text:
program poshalo
  implicit none
  integer::ncpu,ndim,npart,i,j,icpu,ipos,nstar,ngroup
  integer::ncpu2,npart2,ndim2,id_group
  integer,dimension(:),allocatable::group_id,n_group
  real(kind=8),dimension(:),allocatable::xref_group,yref_group,zref_group
  logical,dimension(:),allocatable::flag_group
  real(KIND=8)::mtot,mcut=1000000,contamine
  real(KIND=8)::period=1,xmin=0,xmax=1,ymin=0,ymax=1,zmin=0,zmax=1,rvir
  real(kind=8),dimension(:),allocatable::m_group,x_group,y_group,z_group,u_group,v_group,w_group,mpure_group
  real(KIND=8),dimension(:,:),allocatable::x,v
  real(KIND=8),dimension(:),allocatable::m,age
  integer,dimension(:),allocatable::id
  character(LEN=5)::nstring,ncharcpu
  character(LEN=80)::directory,file_groupe_in
  character(LEN=128)::nomfich,repository
  logical::ok

  call read_params

  ! Reading groups
  file_groupe_in=TRIM(directory)//'.tag'
  Open(15,file=trim(file_groupe_in),FORM='UNFORMATTED')
  Read(15)npart, ngroup
  Write(*,*)'npart=',npart
  Write(*,*)'ngroup=',ngroup
  allocate(group_id(1:npart))
  allocate(n_group(1:ngroup),m_group(1:ngroup),mpure_group(1:ngroup))
  allocate(flag_group(1:ngroup))
  allocate(xref_group(1:ngroup),yref_group(1:ngroup),zref_group(1:ngroup))
  allocate(x_group(1:ngroup),y_group(1:ngroup),z_group(1:ngroup))
  allocate(u_group(1:ngroup),v_group(1:ngroup),w_group(1:ngroup))
  Read(15)(group_id(j),j=1,npart)
  Close(15)
  x_group=0.0
  y_group=0.0
  z_group=0.0
  u_group=0.0
  v_group=0.0
  w_group=0.0
  m_group=0.0
  mpure_group=0.0
  n_group=0.0
  flag_group=.true.  ! True if group reference particle is UNSET
  xref_group=0.0
  yref_group=0.0
  zref_group=0.0
  group_id=group_id+1  ! WARNING: convert from 0 to ng-1 to 1 to ng

  !-----------------------------------------------
  ! Lecture du fichier particules au format RAMSES
  !-----------------------------------------------
  ipos=INDEX(repository,'output_')
  nstring=repository(ipos+7:ipos+13)
  nomfich=TRIM(repository)//'/part_'//TRIM(nstring)//'.out00001'
  inquire(file=nomfich, exist=ok) ! verify input file
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
     stop
  endif

  nomfich=TRIM(repository)//'/info_'//TRIM(nstring)//'.txt'
  inquire(file=nomfich, exist=ok) ! verify input file
  if ( .not. ok ) then
     print *,TRIM(nomfich)//' not found.'
     stop
  endif
  open(unit=10,file=nomfich,form='formatted',status='old')
  read(10,'(13X,I11)')ncpu
  read(10,'(13X,I11)')ndim
  close(10)

  npart=0
  do icpu=1,ncpu
     call title(icpu,ncharcpu)
     nomfich=TRIM(repository)//'/part_'//TRIM(nstring)//'.out'//TRIM(ncharcpu)
     open(unit=1,file=nomfich,status='old',form='unformatted')
     read(1)ncpu2
     read(1)ndim2
     read(1)npart2
     read(1)
     read(1)nstar
     close(1)
     npart=npart+npart2
  end do
  write(*,*)'Found ',npart,' particles.'
  if(nstar>0)then
     write(*,*)'Discard star particles.'
  endif

  mtot=0.0d0
  npart=0
  do icpu=1,ncpu
     call title(icpu,ncharcpu)
     nomfich=TRIM(repository)//'/part_'//TRIM(nstring)//'.out'//TRIM(ncharcpu)
     open(unit=1,file=nomfich,status='old',form='unformatted')
     !     write(*,*)'Processing file '//TRIM(nomfich)
     read(1)ncpu2
     read(1)ndim2
     read(1)npart2
     read(1)
     read(1)
     read(1)
     read(1)
     read(1)
     allocate(m(1:npart2))
     if(nstar>0)then
        allocate(age(1:npart2))
        allocate(id(1:npart2))
     endif
     allocate(x(1:npart2,1:ndim2))
     allocate(v(1:npart2,1:ndim2))
     ! Read position
     do i=1,ndim
        read(1)m
        x(1:npart2,i)=m
     end do
     ! Read velocity
     do i=1,ndim
        read(1)m
        v(1:npart2,i)=m
     end do
     ! Read mass
     read(1)m
     if(nstar>0)then
        read(1)id
        read(1) ! Skip level
        read(1)age
     endif
     close(1)

     do i=1,npart2
        if(nstar>0) then
           if(age(i)/=0.or.id(i)<1) cycle
        end if
        npart=npart+1
        id_group=group_id(npart)
        if(id_group>0.and.id_group.LE.ngroup)then
           if(flag_group(id_group)) then
              flag_group(id_group)=.false.
              xref_group(id_group)=x(i,1)
              yref_group(id_group)=x(i,2)
              zref_group(id_group)=x(i,3)
           end if
           n_group(id_group)=n_group(id_group)+1
           m_group(id_group)=m_group(id_group)+m(i)
           if(m(i)<=mcut)mpure_group(id_group)=mpure_group(id_group)+m(i)
           x_group(id_group)=x_group(id_group)+m(i)*shift1d(x(i,1),xref_group(id_group),period)
           y_group(id_group)=y_group(id_group)+m(i)*shift1d(x(i,2),yref_group(id_group),period)
           z_group(id_group)=z_group(id_group)+m(i)*shift1d(x(i,3),zref_group(id_group),period)
           u_group(id_group)=u_group(id_group)+m(i)*v(i,1)
           v_group(id_group)=v_group(id_group)+m(i)*v(i,2)
           w_group(id_group)=w_group(id_group)+m(i)*v(i,3)
        endif
     enddo
     deallocate(x,m,v)
     if(nstar>0)deallocate(age,id)
  end do

  open(18,file=TRIM(directory)//'.pos')
  write(*,'(A100)')'   #   npart       mass  cont.frac         xc         yc         zc         uc         vc         wc'
  write(18,'(A100)')'   #   npart       mass  cont.frac         xc         yc         zc         uc         vc         wc'
  do i=1,ngroup
     x_group(i)=shift1d(x_group(i)/m_group(i),0.5d0,period) ! Ensure group center is within bounds
     y_group(i)=shift1d(y_group(i)/m_group(i),0.5d0,period)
     z_group(i)=shift1d(z_group(i)/m_group(i),0.5d0,period)
     u_group(i)=u_group(i)/m_group(i)
     v_group(i)=v_group(i)/m_group(i)
     w_group(i)=w_group(i)/m_group(i)
     if(         x_group(i).gt.xmin.and.x_group(i).lt.xmax &
          & .and.y_group(i).gt.ymin.and.y_group(i).lt.ymax &
          & .and.z_group(i).gt.zmin.and.z_group(i).lt.zmax)then
        rvir=(m_group(i)/200./(4./3.*3.1415926))**(1./3.)
        contamine=(m_group(i)-mpure_group(i))/m_group(i)
        if(i<=10)write(*,'(I5,A,I7,A,1PE10.3,A,1E10.3,A,1E10.3,A,1E10.3,A,1E10.3,A,1E10.3,A,1E10.3,A,1E10.3)')&
             & i,' ',n_group(i),' ',m_group(i),' ',contamine,' ',x_group(i),' ',y_group(i),' ',z_group(i),' ', &
             & u_group(i),' ',v_group(i),' ',w_group(i)
                write(18,'(I5,A,I7,A,1PE10.3,A,1E10.3,A,1E10.3,A,1E10.3,A,1E10.3,A,1E10.3,A,1E10.3,A,1E10.3)') &
             & i,' ',n_group(i),' ',m_group(i),' ',contamine,' ',x_group(i),' ',y_group(i),' ',z_group(i),' ', &
             & u_group(i),' ',v_group(i),' ',w_group(i)
     endif
  end do
  close(18)

contains

  real(kind=8) function shift1d(x, xref, per)
      !> Shifts x by a multiple of per such that abs(x-xref) is minimal
      real(kind=8) :: x, xref, per
      integer :: ishift
      if(per>0d0) then
         ishift = floor((xref-x)/per+0.5)
         shift1d = x+ishift*per
      else
         shift1d = x
      end if
  end function shift1d

  subroutine read_params

      implicit none

      integer       :: i,n
      integer       :: iargc
      character(len=4)   :: opt
      character(len=128) :: arg

      n = iargc()
      if (n < 4) then
         print *, 'usage: poshalo -inp ramses_input -pre hop_prefix'
         print *, 'ex: poshalo -inp output_00001 -pre zregroup'
         stop
      end if

      do i = 1,n,2
         call getarg(i,opt)
         if (i == n) then
            print '("option ",a4," has no argument")', opt
            stop 2
         end if
         call getarg(i+1,arg)
         select case (opt)
         case ('-pre')
            directory = trim(arg)
         case ('-inp')
            repository = trim(arg)
         case ('-xmi')
            read (arg,*) xmin
         case ('-xma')
            read (arg,*) xmax
         case ('-ymi')
            read (arg,*) ymin
         case ('-yma')
            read (arg,*) ymax
         case ('-zmi')
            read (arg,*) zmin
         case ('-cut')
            read (arg,*) mcut
         case ('-zma')
            read (arg,*) zmax
         case ('-per')
            read (arg,*) period
         case default
            print '("unknown option ",a2," ignored")', opt
         end select
      end do

      return

    end subroutine read_params
end program poshalo

subroutine title(n,nstring)
  implicit none
  integer::n
  character*5::nstring

  character*1::nchar1
  character*2::nchar2
  character*3::nchar3
  character*4::nchar4
  character*5::nchar5

  if(n.ge.10000)then
     write(nchar5,'(i5)') n
     nstring = nchar5
  elseif(n.ge.1000)then
     write(nchar4,'(i4)') n
     nstring = '0'//nchar4
  elseif(n.ge.100)then
     write(nchar3,'(i3)') n
     nstring = '00'//nchar3
  elseif(n.ge.10)then
     write(nchar2,'(i2)') n
     nstring = '000'//nchar2
  else
     write(nchar1,'(i1)') n
     nstring = '0000'//nchar1
  endif
end subroutine title

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rknop/amuse


text:
      PROGRAM TRIPLE_NESTED
      IMPLICIT NONE

      INTEGER MONTH, DATE, DAY

      MONTH = 1

      DO WHILE (MONTH <= 12)
        DATE = 1
        DO WHILE (DATE <= 7)
          DAY = 1
          DO WHILE (DAY <= DATE)
            IF (DAY == 1) THEN
              PRINt *, "MONTH: ", MONTH, " DAY: ", DAY, ", SUNDAY"
            ELSE IF (DAY == 2) THEN
              PRINt *, "MONTH: ", MONTH, " DAY: ", DAY, ", MONDAY"
            ELSE IF (DAY == 3) THEN
              PRINt *, "MONTH: ", MONTH, " DAY: ", DAY, ", TUESDAY"
            ELSE IF (DAY == 4) THEN
              PRINt *, "MONTH: ", MONTH, " DAY: ", DAY, ", WEDNESDAY"
            ELSE IF (DAY == 5) THEN
              PRINt *, "MONTH: ", MONTH, " DAY: ", DAY, ", THURSDAY"
            ELSE IF (DAY == 6) THEN
              PRINt *, "MONTH: ", MONTH, " DAY: ", DAY, ", FRIDAY"
            ELSE
              PRINt *, "MONTH: ", MONTH, " DAY: ", DAY, ", SATURDAY"
            END IF
            DAY = DAY + 1
          END DO
          DATE = DATE +1
        END DO
        MONTH = MONTH + 1
      END DO
      END PROGRAM TRIPLE_NESTED

ext: f


text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! 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.
!
! 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.
!
!--------------------------------------------------------------------
!
!     This module defines data structures for Bueno-Orovio cellular
!     activation model for cardiac electrophysiology.
!
!--------------------------------------------------------------------

      MODULE BOMOD
      USE TYPEMOD
      IMPLICIT NONE

      PRIVATE :: ISZERO, STEP, DELTA

      CONTAINS
!--------------------------------------------------------------------
      SUBROUTINE BO_INIT(nX, X)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: nX
      REAL(KIND=RKIND), INTENT(INOUT) :: X(nX)

      INCLUDE "PARAMS_BO.f"

      X(1) = Voffset
      X(2) = 1._RKIND
      X(3) = 1._RKIND
      X(4) = 0._RKIND

      RETURN
      END SUBROUTINE BO_INIT
!--------------------------------------------------------------------
!     Time integration performed using Forward Euler method
      SUBROUTINE BO_INTEGFE(imyo, nX, X, Ts, Ti, Istim, Ksac, RPAR)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: imyo, nX
      REAL(KIND=RKIND), INTENT(IN) :: Ts, Ti, Istim, Ksac
      REAL(KIND=RKIND), INTENT(INOUT) :: X(nX), RPAR(5)

      INCLUDE "PARAMS_BO.f"

      REAL(KIND=RKIND) :: t, dt, f(nX), fext, Isac

      t    = Ts / Tscale
      dt   = Ti / Tscale

      Isac = Ksac * (Vrest - X(1))
      fext = (Istim + Isac) * Tscale / Vscale

      X(1) = (X(1) - Voffset)/Vscale

      CALL BO_GETF(imyo, nX, X, f, fext, RPAR)
      X(:) = X(:) + dt*f(:)

      X(1) = X(1)*Vscale + Voffset

      RETURN
      END SUBROUTINE BO_INTEGFE
!--------------------------------------------------------------------
!     Time integration performed using 4th order Runge-Kutta method
      SUBROUTINE BO_INTEGRK(imyo, nX, X, Ts, Ti, Istim, Ksac, RPAR)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: imyo, nX
      REAL(KIND=RKIND), INTENT(IN) :: Ts, Ti, Istim, Ksac
      REAL(KIND=RKIND), INTENT(INOUT) :: X(nX), RPAR(5)

      INCLUDE "PARAMS_BO.f"

      REAL(KIND=RKIND) :: t, dt, dt6, fext, Isac, Xrk(nX), frk(nX,4)

      t    = Ts / Tscale
      dt   = Ti / Tscale
      dt6  = dt/6._RKIND

      Isac = Ksac * (Vrest - X(1))
      fext = (Istim + Isac) * Tscale / Vscale
      X(1) = (X(1) - Voffset)/Vscale

!     RK4: 1st pass
      Xrk  = X
      CALL BO_GETF(imyo, nX, Xrk, frk(:,1), fext, RPAR)

!     RK4: 2nd pass
      Xrk  = X + 0.5_RKIND*dt*frk(:,1)
      CALL BO_GETF(imyo, nX, Xrk, frk(:,2), fext, RPAR)

!     RK4: 3rd pass
      Xrk  = X + 0.5_RKIND*dt*frk(:,2)
      CALL BO_GETF(imyo, nX, Xrk, frk(:,3), fext, RPAR)

!     RK4: 4th pass
      Xrk  = X + dt*frk(:,3)
      CALL BO_GETF(imyo, nX, Xrk, frk(:,4), fext, RPAR)

      X = X + dt6*(frk(:,1) + 2._RKIND*(frk(:,2) + frk(:,3)) + frk(:,4))

      X(1) = X(1)*Vscale + Voffset

      RETURN
      END SUBROUTINE BO_INTEGRK
!--------------------------------------------------------------------
!     Time integration performed using Crank-Nicholson method
      SUBROUTINE BO_INTEGCN2(imyo, nX, Xn, Ts, Ti, Istim, Ksac, IPAR,
     2   RPAR)
      USE MATFUN
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: imyo, nX
      INTEGER(KIND=IKIND), INTENT(INOUT) :: IPAR(2)
      REAL(KIND=RKIND), INTENT(IN) :: Ts, Ti, Istim, Ksac
      REAL(KIND=RKIND), INTENT(INOUT) :: Xn(nX), RPAR(5)

      INCLUDE "PARAMS_BO.f"

      REAL(KIND=RKIND), PARAMETER :: eps = EPSILON(eps)

      INTEGER(KIND=IKIND) :: i, k, itMax
      LOGICAL :: l1, l2, l3
      REAL(KIND=RKIND) :: t, dt, fext, atol, rtol, Xk(nX), fn(nX),
     2   fK(nX), rK(nX), Im(nX,nX), JAC(nX,nX), rmsA, rmsR, Isac

      itMax = IPAR(1)
      atol  = RPAR(1)
      rtol  = RPAR(2)

      t     = Ts / Tscale
      dt    = Ti / Tscale
      Isac  = Ksac * (Vrest - Xn(1))
      fext  = (Istim + Isac) * Tscale / Vscale

      Xn(1) = (Xn(1) - Voffset)/Vscale
      Im    = MAT_ID(nX)

      CALL BO_GETF(imyo, nX, Xn, fn, fext, RPAR)

      k  = 0
      Xk = Xn
      l1 = .FALSE.
      l2 = .FALSE.
      l3 = .FALSE.
      t  = Ts + dt
      DO
         k = k + 1
         CALL BO_GETF(imyo, nX, Xk, fk, fext, RPAR)
         rK(:) = Xk(:) - Xn(:) - 0.5_RKIND*dt*(fk(:) + fn(:))

         rmsA = 0._RKIND
         rmsR = 0._RKIND
         DO i=1, nX
            rmsA = rmsA + rK(i)**2._RKIND
            rmsR = rmsR + ( rK(i) / (Xk(i)+eps) )**2._RKIND
         END DO
         rmsA = SQRT(rmsA/REAL(nX, KIND=RKIND))
         rmsR = SQRT(rmsR/REAL(nX, KIND=RKIND))

         l1   = k .GT. itMax
         l2   = rmsA .LE. atol
         l3   = rmsR .LE. rtol
         IF (l1 .OR. l2 .OR. l3) EXIT

         CALL BO_GETJ(imyo, nX, Xk, JAC)
         JAC   = Im - 0.5_RKIND*dt*JAC
         JAC   = MAT_INV(JAC, nX)
         rK(:) = MATMUL(JAC, rK)
         Xk(:) = Xk(:) - rK(:)
      END DO
      Xn(:) = Xk(:)
      CALL BO_GETF(imyo, nX, Xn, fn, fext, RPAR)
      Xn(1) = Xn(1)*Vscale + Voffset

      IF (.NOT.l2 .AND. .NOT.l3) IPAR(2) = IPAR(2) + 1

      RETURN
      END SUBROUTINE BO_INTEGCN2
!--------------------------------------------------------------------
      SUBROUTINE BO_GETF(i, n, X, f, fext, RPAR)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: i, n
      REAL(KIND=RKIND), INTENT(IN) :: X(n), fext
      REAL(KIND=RKIND), INTENT(OUT) :: f(n)
      REAL(KIND=RKIND), INTENT(INOUT) :: RPAR(5)

      INCLUDE "PARAMS_BO.f"

      REAL(KIND=RKIND) :: u, v, w, s, H_uv, H_uw, H_umv, H_uo, taum_v,
     2   taum_w, tau_so, tau_s, tau_o, v_inf, w_inf, I_fi, I_so, I_si

!     Create local copies of the 4 state variables
      u = X(1)
      v = X(2)
      w = X(3)
      s = X(4)

!     Define step functions
      H_uv  = STEP(u - theta_v(i))
      H_uw  = STEP(u - theta_w(i))
      H_umv = STEP(u - thetam_v(i))
      H_uo  = STEP(u - theta_o(i))

!     Define additional constants
      taum_v = (1._RKIND-H_umv)*taum_v1(i) + H_umv*taum_v2(i)
      taum_w = taum_w1(i) + 0.5_RKIND*(taum_w2(i)-taum_w1(i))*
     2   (1._RKIND + TANH(km_w(i)*(u-um_w(i))))
      tau_so = tau_so1(i) + 0.5_RKIND*(tau_so2(i)-tau_so1(i))*
     2   (1._RKIND+DTANH(k_so(i)*(u-u_so(i))))
      tau_s  = (1._RKIND-H_uw)*tau_s1(i) + H_uw*tau_s2(i)
      tau_o  = (1._RKIND-H_uo)*tau_o1(i) + H_uo*tau_o2(i)
      v_inf  = (1._RKIND-H_umv)
      w_inf  = (1._RKIND-H_uo)*(1._RKIND - u/tau_winf(i)) +
     2   H_uo*ws_inf(i)

!     Compute RHS of state variable equations
      I_fi = -v*H_uv*(u-theta_v(i))*(u_u(i) - u)/tau_fi(i)
      I_so =  (u-u_o(i))*(1._RKIND-H_uw)/tau_o + H_uw/tau_so
      I_si = -H_uw*w*s/tau_si(i)

      f(1) = -(I_fi + I_so + I_si + fext)

      f(2) = (1._RKIND-H_uv)*(v_inf-v)/taum_v - H_uv*v/taup_v(i)

      f(3) = (1._RKIND-H_uw)*(w_inf-w)/taum_w - H_uw*w/taup_w(i)

      f(4) = (0.5_RKIND*(1._RKIND + TANH(k_s(i)*(u-u_s(i))))-s)/tau_s

      RPAR(3) = I_fi
      RPAR(4) = I_so
      RPAR(5) = I_si

      RETURN
      END SUBROUTINE BO_GETF
!--------------------------------------------------------------------
      SUBROUTINE BO_GETJ(i, n, X, JAC)
      IMPLICIT NONE
      INTEGER(KIND=IKIND), INTENT(IN) :: i, n
      REAL(KIND=RKIND), INTENT(IN) :: X(n)
      REAL(KIND=RKIND), INTENT(OUT) :: JAC(n,n)

      INCLUDE "PARAMS_BO.f"

      REAL(KIND=RKIND) :: u, v, w, s, H_uv, H_uw, H_umv, H_uo, D_uw,
     2   D_uv, taum_v, taum_w, tau_so, tau_s, tau_o, v_inf, w_inf, n1,
     3   n2, n3

!     Create local variables
      u = X(1)
      v = X(2)
      w = X(3)
      s = X(4)

!     Define step functions
      H_uv  = STEP(u - theta_v(i))
      H_uw  = STEP(u - theta_w(i))
      H_umv = STEP(u - thetam_v(i))
      H_uo  = STEP(u - theta_o(i))

!     Define delta functions
      D_uw  = DELTA(u - theta_w(i))
      D_uv  = DELTA(u - theta_v(i))

!     Define additional constants
      taum_v = (1._RKIND-H_umv)*taum_v1(i) + H_umv*taum_v2(i)
      taum_w = taum_w1(i) + 0.5_RKIND*(taum_w2(i)-taum_w1(i))*
     2   (1._RKIND + TANH(km_w(i)*(u-um_w(i))))
      tau_so = tau_so1(i) + 0.5_RKIND*(tau_so2(i)-tau_so1(i))*
     2   (1._RKIND+DTANH(k_so(i)*(u-u_so(i))))
      tau_s  = (1._RKIND-H_uw)*tau_s1(i) + H_uw*tau_s2(i)
      tau_o  = (1._RKIND-H_uo)*tau_o1(i) + H_uo*tau_o2(i)
      v_inf  = (1._RKIND-H_umv)
      w_inf  = (1._RKIND-H_uo)*(1._RKIND - u/tau_winf(i)) +
     2   H_uo*ws_inf(i)

!     Define Jacobian
      JAC(:,:) = 0._RKIND

      n1 = v*H_uv*(u_u(i) + theta_v(i) - 2._RKIND*u)/tau_fi(i)
      n2 = -(1._RKIND-H_uw)/tau_fi(i)
      n3 = (-1._RKIND/tau_so + (theta_w(i)-u_o(i))/tau_o +
     2   w*s/tau_si(i))*D_uw
      JAC(1,1) = n1 + n2 + n3

      JAC(1,2) = H_uv*(u-theta_v(i))*(u_u(i)-u)/tau_fi(i)

      n1 = H_uw/tau_si(i)
      JAC(1,3) = n1*s
      JAC(1,4) = n1*w

      n1 = -1._RKIND/taum_v
      n2 = -1._RKIND/taup_v(i)
      JAC(2,1) = ((v_inf-v)*n1 + v*n2)*D_uv
      JAC(2,2) = (1._RKIND-H_uv)*n1 + H_uv*n2

      n1 = -1._RKIND/taum_w
      n2 = -1._RKIND/taup_w(i)
      JAC(3,1) = ((w_inf-w)*n1 + w*n2)*D_uw
      JAC(3,3) = (1._RKIND-H_uw)*n1 + H_uw*n2

      n1 = COSH(k_s(i)*(u-u_s(i)))
      n2 = 1._RKIND/(n1*n1)
      n3 = 1._RKIND/tau_s
      JAC(4,1) = 0.5_RKIND*k_s(i)*n2*n3
      JAC(4,4) = -n3

      RETURN
      END SUBROUTINE BO_GETJ
!--------------------------------------------------------------------
!     Compute activation force for electromechanics based on active
!     stress model
      SUBROUTINE BO_ACTVSTRS(X, dt, Tact, epsX)
      IMPLICIT NONE
      REAL(KIND=RKIND), INTENT(IN) :: X, dt
      REAL(KIND=RKIND), INTENT(OUT) :: epsX
      REAL(KIND=RKIND), INTENT(INOUT) :: Tact

      INCLUDE "PARAMS_BO.f"

      REAL(KIND=RKIND) :: nr

      epsX = EXP(-EXP(-xi_T*(X - Vcrit)))
      epsX = eps_0 + (eps_i - eps_0)*epsX
      nr   = Tact + epsX*dt*eta_T*(X - Vrest)
      Tact = nr / (1._RKIND + epsX*dt)

      RETURN
      END SUBROUTINE BO_ACTVSTRS
!--------------------------------------------------------------------
!     Compute macroscopic fiber strain based on sacromere force-length
!     relationship and slow inward current variable (s)
      SUBROUTINE BO_ACTVSTRN(c, I4f, dt, gf)
      IMPLICIT NONE
      REAL(KIND=RKIND), INTENT(IN) :: c, I4f, dt
      REAL(KIND=RKIND), INTENT(INOUT) :: gf

      INCLUDE "PARAMS_BO.f"

      REAL(KIND=RKIND) :: SL, Fa, rtmp

!     fiber length
      SL = I4f * SL0

!     Sacromere force-length relationship
      IF (SL.GE.SLmin .AND. SL.LE.SLmax) THEN
         SL = 0.5_RKIND*f0 + fc1*COS(SL) + fs1*SIN(SL) +
     2      fc2*COS(2._RKIND*SL) + fs2*SIN(2._RKIND*SL)  +
     3      fc3*COS(3._RKIND*SL) + fs3*SIN(3._RKIND*SL)
      ELSE
         SL = 0._RKIND
      END IF

!     Active force
      Fa   = alFa * (c-c0)*(c-c0) * SL

      rtmp = 2._RKIND*I4f*(1._RKIND/(1._RKIND+gf)**3._RKIND - 1._RKIND)
      gf = gf + dt*(Fa + rtmp)/(mu_C * c * c)

      RETURN
      END SUBROUTINE BO_ACTVSTRN
!--------------------------------------------------------------------
      FUNCTION STEP(r)
      IMPLICIT NONE
      REAL(KIND=RKIND), INTENT(IN) :: r
      REAL(KIND=RKIND) STEP

      IF (r .LT. 0) THEN
         STEP = 0._RKIND
      ELSE
         STEP = 1._RKIND
      END IF

      RETURN
      END FUNCTION STEP
!--------------------------------------------------------------------
      FUNCTION DELTA(r)
      IMPLICIT NONE
      REAL(KIND=RKIND), INTENT(IN) :: r
      REAL(KIND=RKIND) DELTA

      DELTA = 0._RKIND
      IF (ISZERO(r)) DELTA = 1._RKIND

      RETURN
      END FUNCTION DELTA
!--------------------------------------------------------------------
      FUNCTION ISZERO(ia)
      IMPLICIT NONE
      REAL(KIND=RKIND), INTENT(IN) :: ia
      LOGICAL ISZERO

      REAL(KIND=RKIND), PARAMETER :: epsil = EPSILON(epsil)
      REAL(KIND=RKIND) a, b, nrm

      a   = ABS(ia)
      b   = 0._RKIND
      nrm = MAX(a,epsil)

      ISZERO = .FALSE.
      IF ((a-b)/nrm .LT. 1._RKIND*epsil) ISZERO = .TRUE.

      RETURN
      END FUNCTION ISZERO
!--------------------------------------------------------------------
      END MODULE BOMOD
!####################################################################

max_stars_repo_name: ktbolt/svFSI

max_issues_repo_name: ktbolt/svFSI

max_forks_repo_name: ktbolt/svFSI


text:
      program main

c*********************************************************************72
c
cc MAIN is the main program for CG_PRB.
c
c  Discussion:
c
c    CG_PRB tests the CG library.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    By Sourangshu Ghosh
c
c  Author:
c
c    Sourangshu Ghosh
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CG_PRB'
      write ( *, '(a)' ) '  FORTRAN77 version'
      write ( *, '(a)' ) '  Test the CG library.'

      call test01 ( )
      call test02 ( )
      call test023 ( )
      call test025 ( )
      call test03 ( )
      call test04 ( )
      call test05 ( )
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'CG_PRB'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ' '
      call timestamp ( )

      stop
      end
      subroutine test01 ( )

c*********************************************************************72
c
cc TEST01 tests R8GE_CG for a full storage matrix.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    By Sourangshu Ghosh
c
c  Author:
c
c    Sourangshu Ghosh
c
      implicit none

      integer n
      parameter ( n = 10 )

      double precision a(n,n)
      double precision b(n)
      double precision e_norm
      integer i
      double precision r(n)
      double precision r_norm
      double precision r8vec_diff_norm
      double precision r8vec_norm
      integer seed
      double precision x1(n)
      double precision x2(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'TEST01'
      write ( *, '(a)' )
     &  '  Test R8GE_CG, applying CG to a full storage matrix.'
c
c  Choose a random positive definite symmetric matrix A.
c
      seed = 123456789
      call pds_random ( n, seed, a )
c
c  Choose a random solution.
c
      seed = 123456789
      call r8vec_uniform_01 ( n, seed, x1 )
c
c  Compute the corresponding right hand side.
c
      call r8ge_mv ( n, n, a, x1, b )
c
c  Call the CG routine.
c
      do i = 1, n
        x2(i) = 1.0D+00
      end do

      call r8ge_cg ( n, a, b, x2 )
c
c  Compute the residual.
c
      call r8ge_resid ( n, n, a, x2, b, r )
      r_norm = r8vec_norm ( n, r )
c
c  Compute the error.
c
      e_norm = r8vec_diff_norm ( n, x1, x2 )
c
c  Report.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Number of variables N = ', n
      write ( *, '(a,g14.6)' ) '  Norm of residual ||Ax-b|| = ', r_norm
      write ( *, '(a,g14.6)' ) '  Norm of error ||x1-x2|| = ', e_norm

      return
      end
      subroutine test02 ( )

c*********************************************************************72
c
cc TEST02 tests R83_CG.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    By Sourangshu Ghosh
c
c  Author:
c
c    Sourangshu Ghosh
c
      implicit none

      integer n
      parameter ( n = 10 )

      double precision a(3,n)
      double precision b(n)
      double precision e_norm
      integer i
      double precision r(n)
      double precision r_norm
      double precision r8vec_diff_norm
      double precision r8vec_norm
      integer seed
      double precision x1(n)
      double precision x2(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'TEST02'
      write ( *, '(a)' ) '  Test R83_CG, applying CG to an R83 matrix.'
c
c  Let A be the -1 2 -1 matrix.
c
      call r83_dif2 ( n, n, a )
c
c  Choose a random solution.
c
      seed = 123456789
      call r8vec_uniform_01 ( n, seed, x1 )
c
c  Compute the corresponding right hand side.
c
      call r83_mv ( n, n, a, x1, b )
c
c  Call the CG routine.
c
      do i = 1, n
        x2(i) = 1.0D+00
      end do

      call r83_cg ( n, a, b, x2 )
c
c  Compute the residual.
c
      call r83_resid ( n, n, a, x2, b, r )
      r_norm = r8vec_norm ( n, r )
c
c  Compute the error.
c
      e_norm = r8vec_diff_norm ( n, x1, x2 )
c
c  Report.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Number of variables N = ', n
      write ( *, '(a,g14.6)' ) '  Norm of residual ||Ax-b|| = ', r_norm
      write ( *, '(a,g14.6)' ) '  Norm of error ||x1-x2|| = ', e_norm

      return
      end
      subroutine test023 ( )

c*********************************************************************72
c
cc TEST023 tests R83S_CG.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    09 July 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 10 )

      double precision a(3)
      double precision b(n)
      double precision e_norm
      integer i
      double precision r(n)
      double precision r_norm
      double precision r8vec_diff_norm
      double precision r8vec_norm
      integer seed
      double precision x1(n)
      double precision x2(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'TEST023'
      write ( *, '(a)' )
     &  '  Test R83S_CG, applying CG to an R83S matrix.'
c
c  Let A be the -1 2 -1 matrix.
c
      call r83s_dif2 ( n, n, a )
c
c  Choose a random solution.
c
      seed = 123456789
      call r8vec_uniform_01 ( n, seed, x1 )
c
c  Compute the corresponding right hand side.
c
      call r83s_mv ( n, n, a, x1, b )
c
c  Call the CG routine.
c
      do i = 1, n
        x2(i) = 1.0D+00
      end do

      call r83s_cg ( n, a, b, x2 )
c
c  Compute the residual.
c
      call r83s_resid ( n, n, a, x2, b, r )
      r_norm = r8vec_norm ( n, r )
c
c  Compute the error.
c
      e_norm = r8vec_diff_norm ( n, x1, x2 )
c
c  Report.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Number of variables N = ', n
      write ( *, '(a,g14.6)' ) '  Norm of residual ||Ax-b|| = ', r_norm
      write ( *, '(a,g14.6)' ) '  Norm of error ||x1-x2|| = ', e_norm

      return
      end
      subroutine test025 ( )

c*********************************************************************72
c
cc TEST025 tests R83T_CG.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    18 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 10 )

      double precision a(3,n)
      double precision b(n)
      double precision e_norm
      integer i
      double precision r(n)
      double precision r_norm
      double precision r8vec_diff_norm
      double precision r8vec_norm
      integer seed
      double precision x1(n)
      double precision x2(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'TEST025'
      write ( *, '(a)' )
     &  '  Test R83T_CG, applying CG to an R83T matrix.'
c
c  Let A be the -1 2 -1 matrix.
c
      call r83t_dif2 ( n, n, a )
c
c  Choose a random solution.
c
      seed = 123456789
      call r8vec_uniform_01 ( n, seed, x1 )
c
c  Compute the corresponding right hand side.
c
      call r83t_mv ( n, n, a, x1, b )
c
c  Call the CG routine.
c
      do i = 1, n
        x2(i) = 1.0D+00
      end do

      call r83t_cg ( n, a, b, x2 )
c
c  Compute the residual.
c
      call r83t_resid ( n, n, a, x2, b, r )
      r_norm = r8vec_norm ( n, r )
c
c  Compute the error.
c
      e_norm = r8vec_diff_norm ( n, x1, x2 )
c
c  Report.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Number of variables N = ', n
      write ( *, '(a,g14.6)' ) '  Norm of residual ||Ax-b|| = ', r_norm
      write ( *, '(a,g14.6)' ) '  Norm of error ||x1-x2|| = ', e_norm

      return
      end
      subroutine test03 ( )

c*********************************************************************72
c
cc TEST03 tests R8PBU_CG.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    02 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer mu
      parameter ( mu = 1 )
      integer n
      parameter ( n = 10 )

      double precision a(mu+1,n)
      double precision b(n)
      double precision e_norm
      integer i
      double precision r(n)
      double precision r_norm
      double precision r8vec_diff_norm
      double precision r8vec_norm
      integer seed
      double precision x1(n)
      double precision x2(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'TEST03'
      write ( *, '(a)' )
     &  '  Test R8PBU_CG, applying CG to an R8PBU matrix.'
c
c  Let A be the -1 2 -1 matrix.
c
      call r8pbu_dif2 ( n, n, mu, a )
c
c  Choose a random solution.
c
      seed = 123456789
      call r8vec_uniform_01 ( n, seed, x1 )
c
c  Compute the corresponding right hand side.
c
      call r8pbu_mv ( n, n, mu, a, x1, b )
c
c  Call the CG routine.
c
      do i = 1, n
        x2(i) = 1.0D+00
      end do

      call r8pbu_cg ( n, mu, a, b, x2 )
c
c  Compute the residual.
c
      call r8pbu_resid ( n, n, mu, a, x2, b, r )
      r_norm = r8vec_norm ( n, r )
c
c  Compute the error.
c
      e_norm = r8vec_diff_norm ( n, x1, x2 )
c
c  Report.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Number of variables N = ', n
      write ( *, '(a,g14.6)' ) '  Norm of residual ||Ax-b|| = ', r_norm
      write ( *, '(a,g14.6)' ) '  Norm of error ||x1-x2|| = ', e_norm

      return
      end
      subroutine test04 ( )

c*********************************************************************72
c
cc TEST04 tests R8SD_CG.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    02 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 10 )
      integer ndiag
      parameter ( ndiag = 2 )

      double precision a(n,ndiag)
      double precision b(n)
      double precision e_norm
      integer i
      integer offset(ndiag)
      double precision r(n)
      double precision r_norm
      double precision r8vec_diff_norm
      double precision r8vec_norm
      integer seed
      double precision x1(n)
      double precision x2(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'TEST04'
      write ( *, '(a)' )
     &  '  Test R8SD_CG, applying CG to an R8SD matrix.'

      offset(1) = 0
      offset(2) = 1
c
c  Set A to the [-1 2 -1] matrix.
c
      call r8sd_dif2 ( n, n, ndiag, offset, a )
c
c  Choose a random solution.
c
      seed = 123456789
      call r8vec_uniform_01 ( n, seed, x1 )
c
c  Compute the corresponding right hand side.
c
      call r8sd_mv ( n, n, ndiag, offset, a, x1, b )
c
c  Call the CG routine.
c
      do i = 1, n
        x2(i) = 1.0D+00
      end do

      call r8sd_cg ( n, ndiag, offset, a, b, x2 )
c
c  Compute the residual.
c
      call r8sd_resid ( n, n, ndiag, offset, a, x2, b, r )
      r_norm = r8vec_norm ( n, r )
c
c  Compute the error.
c
      e_norm = r8vec_diff_norm ( n, x1, x2 )
c
c  Report.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Number of variables N = ', n
      write ( *, '(a,g14.6)' ) '  Norm of residual ||Ax-b|| = ', r_norm
      write ( *, '(a,g14.6)' ) '  Norm of error ||x1-x2|| = ', e_norm

      return
      end
      subroutine test05 ( )

c*********************************************************************72
c
cc TEST05 tests R8SP_CG.
c
c  Licensing:
c
c    This code is distributed under the GNU LGPL license.
c
c  Modified:
c
c    02 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      integer n
      parameter ( n = 10 )
      integer nz_num
      parameter ( nz_num = 3 * n - 2 )

      double precision a(nz_num)
      double precision b(n)
      integer col(nz_num)
      double precision e_norm
      integer i
      double precision r(n)
      integer row(nz_num)
      double precision r_norm
      double precision r8vec_diff_norm
      double precision r8vec_norm
      integer seed
      double precision x1(n)
      double precision x2(n)

      write ( *, '(a)' ) ''
      write ( *, '(a)' ) 'TEST05'
      write ( *, '(a)' )
     &  '  Test R8SP_CG, applying CG to an R8SP matrix.'
c
c  Set A to the [-1 2 -1] matrix.
c
      call r8sp_dif2 ( n, n, nz_num, row, col, a )
c
c  Choose a random solution.
c
      seed = 123456789
      call r8vec_uniform_01 ( n, seed, x1 )
c
c  Compute the corresponding right hand side.
c
      call r8sp_mv ( n, n, nz_num, row, col, a, x1, b )
c
c  Call the CG routine.
c
      do i = 1, n
        x2(i) = 1.0D+00
      end do

      call r8sp_cg ( n, nz_num, row, col, a, b, x2 )
c
c  Compute the residual.
c
      call r8sp_resid ( n, n, nz_num, row, col, a, x2, b, r )
      r_norm = r8vec_norm ( n, r )
c
c  Compute the error.
c
      e_norm = r8vec_diff_norm ( n, x1, x2 )
c
c  Report.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a,i4)' ) '  Number of variables N = ', n
      write ( *, '(a,g14.6)' ) '  Norm of residual ||Ax-b|| = ', r_norm
      write ( *, '(a,g14.6)' ) '  Norm of error ||x1-x2|| = ', e_norm

      return
      end


text:
C  *********************************************************
C  *                                                       *
C  *    SUBROUTINE 03.05/setdif                            *
C  *                                                       *
C  *    PHIGS Validation Tests, produced by NIST           *
C  *                                                       *
C  *********************************************************

      SUBROUTINE SETDIF (BASSIZ, BASEL, SUBSIZ, SUBEL, ANSSIZ, ANSEL)

C  Setdif performs set difference:
C     ansel = basel - subel
C
C  Input parameters:
C     basel  : source set
C     bassiz : size of basel
C     subel  : subtracted set
C     sunsiz : size of subel
C  Output parameter:
C     ansel  : result set
C     anssiz : size of ansel

      INTEGER BASSIZ, BASEL(*), SUBSIZ, SUBEL(*), ANSSIZ, ANSEL(*),
     1        INDX1, INDX2

      ANSSIZ = 0
      DO 100, INDX1 = 1, BASSIZ
         DO 200, INDX2 = 1, SUBSIZ
            IF (BASEL(INDX1) .EQ. SUBEL(INDX2)) THEN
               GOTO 100
            ENDIF
200      CONTINUE
         ANSSIZ = ANSSIZ + 1
         ANSEL(ANSSIZ) = BASEL(INDX1)
100   CONTINUE
      END


text:
      program master

      include 'p4f.h'

      integer i,slaves,type,from,start,done,time,retcde,recvlen
	  real a(600000)
      integer tagcnt, tagdat, tagend
      parameter (tagcnt = 10)
      parameter (tagdat = 20)
      parameter (tagend = 30)

      call p4init()
      call p4crpg()
      call p4setavlbuf(1,100)
      print *,'Master has started'

      n = 50
	  print *,'Enter interation count:'
	  read(*,*) n
      slaves = p4ntotids() - 1
      length = 8

      do 10 i = 1,slaves
         print *,'MASTER sending control message to ',i
      call p4sendr(tagcnt,i,n,length,retcde)
10    continue

C Initialize socket connections between the slaves
      call p4sendr(tagcnt,1,a,length,retcde)
      type = tagcnt
      from = -1
      call p4recv(type,from,a,length,recvlen,retcde)

      print *,'Number of slaves = ',slaves
      print *,' '
	  print *,'SEND to one slave only'
	  do 55 i = 6,19
         length = 2**i
		 time = 0
	     do 51 k = 1,n
      start = p4clock()
      call p4sendr(i,1,a,length,retcde)
      done = p4clock()
            time = done - start + time
51       continue
         time = time / n
         print *,'Total time to send message = ',time,
     &    '  bytes = ',length
C        print *,'Buffer pools------------------'
C        call p4avlbufs()
55    continue

      print *,' '
	  print *,'SENDR to one slave only'
	  do 60 i = 6,19
         length = 2**i
		 time = 0
	     do 61 k = 1,n
      start = p4clock()
      call p4sendr(i,1,a,length,retcde)
      done = p4clock()
            time = done - start + time
61       continue
	     time = time / n
         print *,'Total time to send and acknowledge message = ',time,
     &    '  bytes = ',length
60    continue

      print *,' '
	  print *,'SEND and RECEIVE sent message from one slave only'
	  do 65 i = 6,19
         length = 2**i
         time = 0
	     do 66 k = 1,n
      start = p4clock()
      call p4send(i,1,a,length,retcde)
      type = i
      from = -1
      call p4recv(type,from,a,length,recvlen,retcde)
      done = p4clock()
            time = done - start + time
66       continue
         time = time / n
         print *,'Total time for send and receive = ',time,
     &	 '  bytes = ',length
65    continue

      print *,' '
	  print *,'SEND message and RECEIVE time from one slave'
	  do 67 i = 6,19
         length = 2**i
		 time = 0
		 do 68 k = 1,n
      start = p4clock()
      call p4send(i,1,a,length,retcde)
      type = i
      from = -1
      call p4recv(type,from,done,4,recvlen,retcde)
            time = done - start + time
68       continue
		 time = time / n
         print *,'Total time for send and receive slave = ',time,
     &	 '  bytes = ',length
67    continue

      print *,' '
	  print *,'SEND message around a ring of slaves'
	  do 70 i = 6,19
         length = 2**i
		 time = 0
		 do 72 k = 1,n
      start = p4clock()
      call p4send(i,1,a,length,retcde)
      type = i
      from = -1
      call p4recv(type,from,a,length,recvlen,retcde)
      done = p4clock()
            time = done - start + time
72       continue
         time = time / n
         print *,'Total time for broadcast = ',time,'  bytes = ',length
70    continue

      print *,' '
	  print *,'SENDR message to each slave sequentially'
	  do 80 i = 6,19
         length = 2**i
		 time = 0
		 do 81 k = 1,n
      start = p4clock()
	        do 82 j = 1,slaves
      call p4sendr(i,j,a,length,retcde)
82          continue
      done = p4clock()
         time = done - start + time
81       continue
		 time = time / n
         print *,'Total time for broadcast = ',time,'  bytes = ',length
80    continue

      length = 0
      do 90 i = 1,slaves
	  print *,'Master sending end message to ',i
      call p4sendr(tagend,i,a,length,retcde)
90    continue

999     call p4cleanup()
      print *,'Master exiting and has cleaned up'
      end

max_stars_repo_name: RWTH-OS/MP-MPICH

max_issues_repo_name: RWTH-OS/MP-MPICH

max_forks_repo_name: RWTH-OS/MP-MPICH


text:
*----------------------------------------------------------------------*
      subroutine transform(
     &     label_me_in,label_me_out,label_form,init,
     &     op_info,form_info,str_info,strmap_info,orb_info)
*----------------------------------------------------------------------*
*
*     transform each record on list with label "label_me_in"
*     to a corresp. record on list with label "label_me_out"
*     using given formula (must be factorized)
*
*     op_info:  operator definitions and files
*     str_info: string information (to be passed to subroutines)
*     orb_info: orbital space information (to be passed)
*
*     andreas, Oct. 2014
*
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'
      include 'stdunit.h'
      include 'ioparam.h'
      include 'mdef_operator_info.h'
      include 'def_graph.h'
      include 'def_strinf.h'
      include 'def_strmapinf.h'
      include 'def_orbinf.h'
      include 'def_contraction.h'
      include 'def_formula_item.h'
      include 'mdef_formula_info.h'
      include 'def_dependency_info.h'
      include 'ifc_memman.h'
      include 'ifc_input.h'
      include 'par_opnames_gen.h'

      integer, parameter ::
     &     ntest = 000
      character(len=9), parameter ::
     &     i_am = 'transform'

      character(*) ::
     &     label_form, label_me_in, label_me_out
      logical ::
     &     init
      type(formula_info) ::
     &     form_info
      type(operator_info) ::
     &     op_info
      type(strinf) ::
     &     str_info
      type(strmapinf) ::
     &     strmap_info
      type(orbinf) ::
     &     orb_info

      integer ::
     &     ifree, nout, iout, idx,
     &     idxmel, recst_in, recnd_in, recst_out, recnd_out,
     &     idx_in, idx_out

      type(formula), pointer ::
     &     f_eval
      type(formula_item) ::
     &     fl_eval
      type(dependency_info) ::
     &     depend

      type(me_list), pointer ::
     &     me_in, me_out

      real(8), pointer ::
     &     xret(:)

      integer, external ::
     &     idx_formlist, idx_mel_list

      ifree = mem_setmark('transform')

      if (ntest.ge.100) then
        call write_title(lulog,wst_dbg_subr,'entered transform')
      end if

      ! get lists
      idxmel = idx_mel_list(label_me_in,op_info)
      if (idxmel.le.0)
     &     call quit(1,i_am,'list not found: '//trim(label_me_in))
      me_in => op_info%mel_arr(idxmel)%mel
      idxmel = idx_mel_list(label_me_out,op_info)
      if (idxmel.le.0)
     &     call quit(1,i_am,'list not found: '//trim(label_me_out))
      me_out => op_info%mel_arr(idxmel)%mel

      ! get formula
      idx = idx_formlist(label_form,form_info)
      if (idx.le.0)
     &     call quit(1,i_am,
     &     'did not find formula '//trim(label_form))
      f_eval => form_info%form_arr(idx)%form

      ! read formula
      call read_form_list(f_eval%fhand,fl_eval,.true.)

c dbg
c     call print_form_list_short(lulog,fl_eval,op_info)
c dbg
      ! set dependency info for submitted formula list
      call set_formula_dependencies(depend,fl_eval,op_info)

      nout = depend%ntargets
      allocate(xret(nout))

      recst_in = me_in%fhand%active_records(1)
      recnd_in = me_in%fhand%active_records(2)
      recst_out = me_out%fhand%active_records(1)
      recnd_out = me_out%fhand%active_records(2)

      if (recnd_in-recst_in.ne.recnd_out-recst_out)
     &     call quit(1,i_am,'incompatible number of records')

      idx_out = recst_out-1
      do idx_in = recst_in, recnd_in
        idx_out = idx_out+1
        call switch_mel_record(me_in,idx_in)
        call switch_mel_record(me_out,idx_out)

        ! call the scheduler
        call frm_sched(xret,fl_eval,depend,0,0,
     &               init,.false.,op_info,str_info,strmap_info,orb_info)

        if (iprlvl.ge.10) then
          call write_title(lulog,wst_title,
     &         'norms/values of output operators')
          do iout = 1, nout
            write(lulog,'(">>>",1x,i4," --> ",g20.14)') iout, xret(iout)
          end do
        end if

        ! this was take from "evaluate" without adaption ...
        if (ntest.ge.1000) then
          do iout = 1, nout
            idx = depend%idxlist(iout)
            write(lulog,*) 'dump of result for ',
     &           trim(op_info%mel_arr(idx)%mel%label)
            call wrt_mel_file(lulog,5,
     &           op_info%mel_arr(idx)%mel,
     &           1,op_info%mel_arr(idx)%mel%op%n_occ_cls,
     &           str_info,orb_info)
          end do
        end if

      end do

      deallocate(xret)
      call dealloc_formula_list(fl_eval)
      call clean_formula_dependencies(depend)

      ifree = mem_flushmark()

      return
      end

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
      subroutine clawpack46_qinit(maxmx,maxmy, meqn,mbc,mx,my,
     &      xlower,ylower,dx,dy,q,maux,aux)
      implicit none

      integer meqn, mbc, mx, my, maux, maxmx, maxmy
      double precision xlower, ylower, dx, dy
      double precision q(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, meqn)
      double precision aux(1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, maux)

      integer i, j, mq, blockno, fc2d_clawpack46_get_block
      double precision xc,yc, xp, yp, zp, xlow, ylow, w, tol
      double precision dxc,xm,ym

      double precision q0_physical

      integer*8 cont, get_context
      logical fclaw2d_map_is_used

      integer ii,jj

      integer initchoice
      common /initchoice_comm/ initchoice

      cont = get_context()

      blockno = fc2d_clawpack46_get_block()

      do i = 1-mbc,mx+mbc
          xc = xlower + (i-0.5d0)*dx
          do j = 1-mbc,my+mbc
              yc = ylower + (j-0.5d0)*dy

              if (initchoice .eq. 0) then
c                 # Discontinuous solution
                  xlow = xlower + (i-1)*dx
                  ylow = ylower + (j-1)*dy
                  call cellave2(blockno,xlow,ylow,dx,dy,w)
                  q(i,j,1) = w
              elseif (initchoice .eq. 1) then
                  call fclaw2d_map_c2m(cont,blockno,xc,yc,xp,yp,zp)
                  q(i,j,1) = q0_physical(xp,yp,zp)
              elseif (initchoice .eq. 2) then
                  q(i,j,1) = 1.d0
              elseif (initchoice .eq. 3) then
c                 # Not sure what this one is about
                  if (fclaw2d_map_is_used(cont)) then
                      call fclaw2d_map_c2m(cont,
     &                           blockno,xc,yc,xp,yp,zp)
                  else
                      xp = xc
                      yp = yc
                  endif
                  dxc = 1.d0/8.d0
                  xm = 0.25
                  ym = 0.25
                  w = dxc/2.d0
                  q(i,j,mq) = 0
                  if (blockno .eq.2) then
                      do jj=7,8
                          do ii=1,2
                              q(ii,jj,1) = 1.d0
                          enddo
                      enddo
                  endif
                  if (blockno .eq. 1) then
                      do ii=5,7
c                         q(ii,my,1) = 1.d0
                      enddo
                  endif
              endif
          enddo
      enddo

      return
      end

max_stars_repo_name: scivision/forestclaw

max_issues_repo_name: scivision/forestclaw

max_forks_repo_name: scivision/forestclaw


!
!                Parallel Sparse BLAS  version 3.5
!      (C) Copyright 2006, 2010, 2015, 2017
!        Salvatore Filippone    Cranfield University
!        Alfredo Buttari        CNRS-IRIT, Toulouse
!
!    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. The name of the PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific 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 PSBLAS GROUP OR ITS 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.
!
!
!
! package: psb_c_csr_mat_mod
!
! This module contains the definition of the psb_c_csr_sparse_mat type
! which implements an actual storage format (the CSR in this case) for
! a sparse matrix as well as the related methods (those who are
! specific to the type and could not be defined higher in the
! hierarchy). We are at the bottom level of the inheritance chain.
!
! Please refere to psb_c_base_mat_mod for a detailed description
! of the various methods, and to psb_c_csr_impl for implementation details.
!
module psb_c_csr_mat_mod

  use psb_c_base_mat_mod

  !> \namespace  psb_base_mod  \class  psb_c_csr_sparse_mat
  !! \extends psb_c_base_mat_mod::psb_c_base_sparse_mat
  !!
  !! psb_c_csr_sparse_mat type and the related methods.
  !! This is a very common storage type, and is the default for assembled
  !! matrices in our library
  type, extends(psb_c_base_sparse_mat) :: psb_c_csr_sparse_mat

    !> Pointers to beginning of rows in JA and VAL.
    integer(psb_ipk_), allocatable :: irp(:)
    !> Column indices.
    integer(psb_ipk_), allocatable :: ja(:)
    !> Coefficient values.
    complex(psb_spk_), allocatable :: val(:)

  contains
    procedure, pass(a) :: is_by_rows  => c_csr_is_by_rows
    procedure, pass(a) :: get_size    => c_csr_get_size
    procedure, pass(a) :: get_nzeros  => c_csr_get_nzeros
    procedure, nopass  :: get_fmt     => c_csr_get_fmt
    procedure, pass(a) :: sizeof      => c_csr_sizeof
    procedure, pass(a) :: csmm        => psb_c_csr_csmm
    procedure, pass(a) :: csmv        => psb_c_csr_csmv
    procedure, pass(a) :: inner_cssm  => psb_c_csr_cssm
    procedure, pass(a) :: inner_cssv  => psb_c_csr_cssv
    procedure, pass(a) :: scals       => psb_c_csr_scals
    procedure, pass(a) :: scalv       => psb_c_csr_scal
    procedure, pass(a) :: maxval      => psb_c_csr_maxval
    procedure, pass(a) :: spnmi       => psb_c_csr_csnmi
    procedure, pass(a) :: rowsum      => psb_c_csr_rowsum
    procedure, pass(a) :: arwsum      => psb_c_csr_arwsum
    procedure, pass(a) :: colsum      => psb_c_csr_colsum
    procedure, pass(a) :: aclsum      => psb_c_csr_aclsum
    procedure, pass(a) :: reallocate_nz => psb_c_csr_reallocate_nz
    procedure, pass(a) :: allocate_mnnz => psb_c_csr_allocate_mnnz
    procedure, pass(a) :: cp_to_coo   => psb_c_cp_csr_to_coo
    procedure, pass(a) :: cp_from_coo => psb_c_cp_csr_from_coo
    procedure, pass(a) :: cp_to_fmt   => psb_c_cp_csr_to_fmt
    procedure, pass(a) :: cp_from_fmt => psb_c_cp_csr_from_fmt
    procedure, pass(a) :: mv_to_coo   => psb_c_mv_csr_to_coo
    procedure, pass(a) :: mv_from_coo => psb_c_mv_csr_from_coo
    procedure, pass(a) :: mv_to_fmt   => psb_c_mv_csr_to_fmt
    procedure, pass(a) :: mv_from_fmt => psb_c_mv_csr_from_fmt
    procedure, pass(a) :: csput_a     => psb_c_csr_csput_a
    procedure, pass(a) :: get_diag    => psb_c_csr_get_diag
    procedure, pass(a) :: csgetptn    => psb_c_csr_csgetptn
    procedure, pass(a) :: csgetrow   => psb_c_csr_csgetrow
    procedure, pass(a) :: get_nz_row  => c_csr_get_nz_row
    procedure, pass(a) :: reinit      => psb_c_csr_reinit
    procedure, pass(a) :: trim        => psb_c_csr_trim
    procedure, pass(a) :: print       => psb_c_csr_print
    procedure, pass(a) :: free        => c_csr_free
    procedure, pass(a) :: mold        => psb_c_csr_mold

  end type psb_c_csr_sparse_mat

  private :: c_csr_get_nzeros, c_csr_free,  c_csr_get_fmt, &
       & c_csr_get_size, c_csr_sizeof, c_csr_get_nz_row, &
       & c_csr_is_by_rows

  !> \memberof psb_c_csr_sparse_mat
  !| \see psb_base_mat_mod::psb_base_reallocate_nz
  interface
    subroutine  psb_c_csr_reallocate_nz(nz,a)
      import :: psb_ipk_, psb_c_csr_sparse_mat
      integer(psb_ipk_), intent(in) :: nz
      class(psb_c_csr_sparse_mat), intent(inout) :: a
    end subroutine psb_c_csr_reallocate_nz
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !| \see psb_base_mat_mod::psb_base_reinit
  interface
    subroutine psb_c_csr_reinit(a,clear)
      import :: psb_ipk_, psb_c_csr_sparse_mat
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      logical, intent(in), optional :: clear
    end subroutine psb_c_csr_reinit
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !| \see psb_base_mat_mod::psb_base_trim
  interface
    subroutine  psb_c_csr_trim(a)
      import :: psb_ipk_, psb_c_csr_sparse_mat
      class(psb_c_csr_sparse_mat), intent(inout) :: a
    end subroutine psb_c_csr_trim
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !| \see psb_base_mat_mod::psb_base_mold
  interface
    subroutine psb_c_csr_mold(a,b,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_base_sparse_mat, psb_long_int_k_
      class(psb_c_csr_sparse_mat), intent(in)                  :: a
      class(psb_c_base_sparse_mat), intent(inout), allocatable :: b
      integer(psb_ipk_), intent(out)                           :: info
    end subroutine psb_c_csr_mold
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !| \see psb_base_mat_mod::psb_base_allocate_mnnz
  interface
    subroutine  psb_c_csr_allocate_mnnz(m,n,a,nz)
      import :: psb_ipk_, psb_c_csr_sparse_mat
      integer(psb_ipk_), intent(in) :: m,n
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      integer(psb_ipk_), intent(in), optional :: nz
    end subroutine psb_c_csr_allocate_mnnz
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_print
  interface
    subroutine psb_c_csr_print(iout,a,iv,head,ivr,ivc)
      import :: psb_ipk_, psb_c_csr_sparse_mat
      integer(psb_ipk_), intent(in)               :: iout
      class(psb_c_csr_sparse_mat), intent(in) :: a
      integer(psb_ipk_), intent(in), optional     :: iv(:)
      character(len=*), optional        :: head
      integer(psb_ipk_), intent(in), optional     :: ivr(:), ivc(:)
    end subroutine psb_c_csr_print
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_cp_to_coo
  interface
    subroutine psb_c_cp_csr_to_coo(a,b,info)
      import :: psb_ipk_, psb_c_coo_sparse_mat, psb_c_csr_sparse_mat
      class(psb_c_csr_sparse_mat), intent(in) :: a
      class(psb_c_coo_sparse_mat), intent(inout) :: b
      integer(psb_ipk_), intent(out)            :: info
    end subroutine psb_c_cp_csr_to_coo
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_cp_from_coo
  interface
    subroutine psb_c_cp_csr_from_coo(a,b,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_coo_sparse_mat
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      class(psb_c_coo_sparse_mat), intent(in)    :: b
      integer(psb_ipk_), intent(out)               :: info
    end subroutine psb_c_cp_csr_from_coo
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_cp_to_fmt
  interface
    subroutine psb_c_cp_csr_to_fmt(a,b,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_base_sparse_mat
      class(psb_c_csr_sparse_mat), intent(in)   :: a
      class(psb_c_base_sparse_mat), intent(inout) :: b
      integer(psb_ipk_), intent(out)                       :: info
    end subroutine psb_c_cp_csr_to_fmt
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_cp_from_fmt
  interface
    subroutine psb_c_cp_csr_from_fmt(a,b,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_base_sparse_mat
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      class(psb_c_base_sparse_mat), intent(in)   :: b
      integer(psb_ipk_), intent(out)                        :: info
    end subroutine psb_c_cp_csr_from_fmt
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_mv_to_coo
  interface
    subroutine psb_c_mv_csr_to_coo(a,b,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_coo_sparse_mat
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      class(psb_c_coo_sparse_mat), intent(inout)   :: b
      integer(psb_ipk_), intent(out)            :: info
    end subroutine psb_c_mv_csr_to_coo
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_mv_from_coo
  interface
    subroutine psb_c_mv_csr_from_coo(a,b,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_coo_sparse_mat
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      class(psb_c_coo_sparse_mat), intent(inout) :: b
      integer(psb_ipk_), intent(out)                        :: info
    end subroutine psb_c_mv_csr_from_coo
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_mv_to_fmt
  interface
    subroutine psb_c_mv_csr_to_fmt(a,b,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_base_sparse_mat
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      class(psb_c_base_sparse_mat), intent(inout)  :: b
      integer(psb_ipk_), intent(out)                        :: info
    end subroutine psb_c_mv_csr_to_fmt
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_mv_from_fmt
  interface
    subroutine psb_c_mv_csr_from_fmt(a,b,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_c_base_sparse_mat
      class(psb_c_csr_sparse_mat), intent(inout)  :: a
      class(psb_c_base_sparse_mat), intent(inout) :: b
      integer(psb_ipk_), intent(out)                         :: info
    end subroutine psb_c_mv_csr_from_fmt
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_cp_from
  interface
    subroutine psb_c_csr_cp_from(a,b)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      type(psb_c_csr_sparse_mat), intent(in)   :: b
    end subroutine psb_c_csr_cp_from
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_mv_from
  interface
    subroutine psb_c_csr_mv_from(a,b)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(inout)  :: a
      type(psb_c_csr_sparse_mat), intent(inout) :: b
    end subroutine psb_c_csr_mv_from
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_csput_a
  interface
    subroutine psb_c_csr_csput_a(nz,ia,ja,val,a,imin,imax,jmin,jmax,info,gtl)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      complex(psb_spk_), intent(in)      :: val(:)
      integer(psb_ipk_), intent(in)             :: nz,ia(:), ja(:),&
           &  imin,imax,jmin,jmax
      integer(psb_ipk_), intent(out)            :: info
      integer(psb_ipk_), intent(in), optional   :: gtl(:)
    end subroutine psb_c_csr_csput_a
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_base_mat_mod::psb_base_csgetptn
  interface
    subroutine psb_c_csr_csgetptn(imin,imax,a,nz,ia,ja,info,&
         & jmin,jmax,iren,append,nzin,rscale,cscale)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      integer(psb_ipk_), intent(in)                  :: imin,imax
      integer(psb_ipk_), intent(out)                 :: nz
      integer(psb_ipk_), allocatable, intent(inout)  :: ia(:), ja(:)
      integer(psb_ipk_),intent(out)                  :: info
      logical, intent(in), optional        :: append
      integer(psb_ipk_), intent(in), optional        :: iren(:)
      integer(psb_ipk_), intent(in), optional        :: jmin,jmax, nzin
      logical, intent(in), optional        :: rscale,cscale
    end subroutine psb_c_csr_csgetptn
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_csgetrow
  interface
    subroutine psb_c_csr_csgetrow(imin,imax,a,nz,ia,ja,val,info,&
         & jmin,jmax,iren,append,nzin,rscale,cscale)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      integer(psb_ipk_), intent(in)                  :: imin,imax
      integer(psb_ipk_), intent(out)                 :: nz
      integer(psb_ipk_), allocatable, intent(inout)  :: ia(:), ja(:)
      complex(psb_spk_), allocatable,  intent(inout)    :: val(:)
      integer(psb_ipk_),intent(out)                  :: info
      logical, intent(in), optional        :: append
      integer(psb_ipk_), intent(in), optional        :: iren(:)
      integer(psb_ipk_), intent(in), optional        :: jmin,jmax, nzin
      logical, intent(in), optional        :: rscale,cscale
    end subroutine psb_c_csr_csgetrow
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_csgetblk
  interface
    subroutine psb_c_csr_csgetblk(imin,imax,a,b,info,&
       & jmin,jmax,iren,append,rscale,cscale)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_, psb_c_coo_sparse_mat
      class(psb_c_csr_sparse_mat), intent(in) :: a
      class(psb_c_coo_sparse_mat), intent(inout) :: b
      integer(psb_ipk_), intent(in)                  :: imin,imax
      integer(psb_ipk_),intent(out)                  :: info
      logical, intent(in), optional        :: append
      integer(psb_ipk_), intent(in), optional        :: iren(:)
      integer(psb_ipk_), intent(in), optional        :: jmin,jmax
      logical, intent(in), optional        :: rscale,cscale
    end subroutine psb_c_csr_csgetblk
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_cssv
  interface
    subroutine psb_c_csr_cssv(alpha,a,x,beta,y,info,trans)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      complex(psb_spk_), intent(in)          :: alpha, beta, x(:)
      complex(psb_spk_), intent(inout)       :: y(:)
      integer(psb_ipk_), intent(out)                :: info
      character, optional, intent(in)     :: trans
    end subroutine psb_c_csr_cssv
  end interface
  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_cssm
  interface
    subroutine psb_c_csr_cssm(alpha,a,x,beta,y,info,trans)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      complex(psb_spk_), intent(in)          :: alpha, beta, x(:,:)
      complex(psb_spk_), intent(inout)       :: y(:,:)
      integer(psb_ipk_), intent(out)                :: info
      character, optional, intent(in)     :: trans
    end subroutine psb_c_csr_cssm
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_csmv
  interface
    subroutine psb_c_csr_csmv(alpha,a,x,beta,y,info,trans)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      complex(psb_spk_), intent(in)          :: alpha, beta, x(:)
      complex(psb_spk_), intent(inout)       :: y(:)
      integer(psb_ipk_), intent(out)                :: info
      character, optional, intent(in)     :: trans
    end subroutine psb_c_csr_csmv
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_csmm
  interface
    subroutine psb_c_csr_csmm(alpha,a,x,beta,y,info,trans)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      complex(psb_spk_), intent(in)          :: alpha, beta, x(:,:)
      complex(psb_spk_), intent(inout)       :: y(:,:)
      integer(psb_ipk_), intent(out)                :: info
      character, optional, intent(in)     :: trans
    end subroutine psb_c_csr_csmm
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_maxval
  interface
    function psb_c_csr_maxval(a) result(res)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      real(psb_spk_)         :: res
    end function psb_c_csr_maxval
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_csnmi
  interface
    function psb_c_csr_csnmi(a) result(res)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      real(psb_spk_)         :: res
    end function psb_c_csr_csnmi
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_rowsum
  interface
    subroutine psb_c_csr_rowsum(d,a)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      complex(psb_spk_), intent(out)              :: d(:)
    end subroutine psb_c_csr_rowsum
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_arwsum
  interface
    subroutine psb_c_csr_arwsum(d,a)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      real(psb_spk_), intent(out)              :: d(:)
    end subroutine psb_c_csr_arwsum
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_colsum
  interface
    subroutine psb_c_csr_colsum(d,a)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      complex(psb_spk_), intent(out)              :: d(:)
    end subroutine psb_c_csr_colsum
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_aclsum
  interface
    subroutine psb_c_csr_aclsum(d,a)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      real(psb_spk_), intent(out)              :: d(:)
    end subroutine psb_c_csr_aclsum
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_get_diag
  interface
    subroutine psb_c_csr_get_diag(a,d,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(in) :: a
      complex(psb_spk_), intent(out)     :: d(:)
      integer(psb_ipk_), intent(out)            :: info
    end subroutine psb_c_csr_get_diag
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_scal
  interface
    subroutine psb_c_csr_scal(d,a,info,side)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      complex(psb_spk_), intent(in)      :: d(:)
      integer(psb_ipk_), intent(out)            :: info
      character, intent(in), optional :: side
    end subroutine psb_c_csr_scal
  end interface

  !> \memberof psb_c_csr_sparse_mat
  !! \see psb_c_base_mat_mod::psb_c_base_scals
  interface
    subroutine psb_c_csr_scals(d,a,info)
      import :: psb_ipk_, psb_c_csr_sparse_mat, psb_spk_
      class(psb_c_csr_sparse_mat), intent(inout) :: a
      complex(psb_spk_), intent(in)      :: d
      integer(psb_ipk_), intent(out)            :: info
    end subroutine psb_c_csr_scals
  end interface

contains

  ! == ===================================
  !
  !
  !
  ! Getters
  !
  !
  !
  !
  !
  ! == ===================================

  function c_csr_is_by_rows(a) result(res)
    implicit none
    class(psb_c_csr_sparse_mat), intent(in) :: a
    logical  :: res
    res = .true.

  end function c_csr_is_by_rows

  function c_csr_sizeof(a) result(res)
    implicit none
    class(psb_c_csr_sparse_mat), intent(in) :: a
    integer(psb_long_int_k_) :: res
    res = 8
    res = res + (2*psb_sizeof_sp)  * psb_size(a%val)
    res = res + psb_sizeof_int * psb_size(a%irp)
    res = res + psb_sizeof_int * psb_size(a%ja)

  end function c_csr_sizeof

  function c_csr_get_fmt() result(res)
    implicit none
    character(len=5) :: res
    res = 'CSR'
  end function c_csr_get_fmt

  function c_csr_get_nzeros(a) result(res)
    implicit none
    class(psb_c_csr_sparse_mat), intent(in) :: a
    integer(psb_ipk_) :: res
    res = a%irp(a%get_nrows()+1)-1
  end function c_csr_get_nzeros

  function c_csr_get_size(a) result(res)
    implicit none
    class(psb_c_csr_sparse_mat), intent(in) :: a
    integer(psb_ipk_) :: res

    res = -1

    if (allocated(a%ja)) then
      res = size(a%ja)
    end if
    if (allocated(a%val)) then
      if (res >= 0) then
        res = min(res,size(a%val))
      else
        res = size(a%val)
      end if
    end if

  end function c_csr_get_size

  function  c_csr_get_nz_row(idx,a) result(res)

    implicit none

    class(psb_c_csr_sparse_mat), intent(in) :: a
    integer(psb_ipk_), intent(in)                  :: idx
    integer(psb_ipk_) :: res

    res = 0

    if ((1<=idx).and.(idx<=a%get_nrows())) then
      res = a%irp(idx+1)-a%irp(idx)
    end if

  end function c_csr_get_nz_row

  ! == ===================================
  !
  !
  !
  ! Data management
  !
  !
  !
  !
  !
  ! == ===================================

  subroutine  c_csr_free(a)
    implicit none

    class(psb_c_csr_sparse_mat), intent(inout) :: a

    if (allocated(a%irp)) deallocate(a%irp)
    if (allocated(a%ja)) deallocate(a%ja)
    if (allocated(a%val)) deallocate(a%val)
    call a%set_null()
    call a%set_nrows(izero)
    call a%set_ncols(izero)

    return

  end subroutine c_csr_free

end module psb_c_csr_mat_mod


text:
subroutine do_burn() bind (C)

  use network
  use eos_module
  use eos_type_module
  use burner_module
  use actual_burner_module

  use amrex_fort_module, only : rt => amrex_real
  implicit none

  type (burn_t) :: state_in, state_out

  real(rt)         :: time = 0.0_rt, dt = 1.0e-6_rt

  type (eos_t) :: eos_state

  state_in % rho    = 2.0e7_rt
  state_in % T      = 8.0e9_rt

  state_in % xn(:)  = ONE / nspec

  call burn_to_eos(state_in, eos_state)
  call eos(eos_input_rt, eos_state)
  call eos_to_burn(eos_state, state_in)

  print *, "rho_in: ", state_in % rho
  print *, "T_in: ", state_in % T
  print *, "X_in: ", state_in % xn

  state_out = state_in

  call actual_burner(state_in, state_out, dt, time)

  print *, 'done!'

  print *, "rho_out: ", state_out % rho
  print *, "T_out: ", state_out % T
  print *, "X_out: ", state_out % xn

  print *, "Energy change: ", state_out % e - state_in % e

end subroutine do_burn

max_stars_repo_name: ChrisDeGrendele/Microphysics

max_issues_repo_name: Youhichka/Microphysics

max_forks_repo_name: Youhichka/Microphysics


text:
MODULE m_potdis
CONTAINS
  SUBROUTINE potdis(stars,vacuum,atoms,sphhar, input,cell,sym)
    !
    !     *****************************************************
    !     calculates the root-mean-square distance between input
    !     and output potentials (averaged over unit cell volume)
    !                                 based on code by   c.l.fu
    !     *****************************************************
    !
    USE m_intgr, ONLY : intgr3, intgz0
    USE m_constants
    USE m_loddop
    USE m_cfft
    USE m_types
    IMPLICIT NONE
    TYPE(t_input),INTENT(IN)   :: input
    TYPE(t_vacuum),INTENT(IN)  :: vacuum
    TYPE(t_sym),INTENT(IN)     :: sym
    TYPE(t_stars),INTENT(IN)   :: stars
    TYPE(t_cell),INTENT(IN)    :: cell
    TYPE(t_sphhar),INTENT(IN)  :: sphhar
    TYPE(t_atoms),INTENT(IN)   :: atoms
    !     .. Local Scalars ..
    REAL fact,facv,rhs,sumis,sumz
    COMPLEX phase
    INTEGER i,i1,i2,i3,id2,id3,io,ip,iter,ivac,j,k1,k2,k3,lh,n,&
                 nk12,npz,nt,num,na
    LOGICAL tail
    !     ..
    !     .. Local Arrays ..
    COMPLEX rhpw(stars%ng3,2,2),rhv1(vacuum%nmzxyd,stars%ng2-1,2,2,2)
    REAL rhsp(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,2,2),rhv0(vacuum%nmzd,2,2,2)
    REAL dis(2),disz(vacuum%nmzd),rh(atoms%jmtd)
    REAL af3((2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1),input%jspins),&
         bf3((2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1),input%jspins), &
         bf2((2*stars%mx1+1)*(2*stars%mx2+1),input%jspins),&
         af2((2*stars%mx1+1)*(2*stars%mx2+1),input%jspins)
    REAL pdis(0:4,0:atoms%ntype,2)
    !     ..
    !

    dis = 0.
    pdis=0.0

    tail = .TRUE.
    npz = vacuum%nmz + 1
    nt = (2*stars%mx1+1)*(2*stars%mx2+1)*(2*stars%mx3+1)
    nk12 = (2*stars%mx1+1)*(2*stars%mx2+1)
    fact = cell%omtil/REAL(nt)
    facv = 1.0
    !     ---> reload potentials
    OPEN (9,file='nrp',form='unformatted',status='unknown')
    REWIND 9
    DO  io = 1,2
       CALL loddop(stars,vacuum,atoms,sphhar, input,sym,&
            9, iter,rhsp(:,0:,:,:,io),rhpw(:,:,io), rhv0(:,:,:,io),rhv1(:,:,:,:,io))
    ENDDO
    CLOSE (9)
    IF (input%jspins.EQ.1) THEN
       !       ---> total potential difference
       CALL priv_cdndif(1,1,1,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
    ELSE
       !       ---> total input potential
       CALL priv_cdndif(1,1,2,1,1,1,1.0,1.0,input%film, rhsp,rhpw,rhv0,rhv1)
       !       ---> total output potential
       CALL priv_cdndif(1,1,2,2,2,2,1.0,1.0,input%film, rhsp,rhpw,rhv0,rhv1)
       !       ---> input 'spin' potential
       CALL priv_cdndif(2,1,2,1,1,1,1.0,-2.0,input%film, rhsp,rhpw,rhv0,rhv1)
       !       ---> output 'spin' potential
       CALL priv_cdndif(2,1,2,2,2,2,1.0,-2.0,input%film, rhsp,rhpw,rhv0,rhv1)
       !       ---> potential difference
       CALL priv_cdndif(1,1,1,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
       !       ---> spin potential difference
       CALL priv_cdndif(2,2,2,1,1,2,1.0,-1.0,input%film, rhsp,rhpw,rhv0,rhv1)
    END IF
    DO  num = 1,input%jspins
       !     ----> m.t. part
       na = 1
       DO n = 1,atoms%ntype
          rh(:atoms%jri(n))=rhsp(:atoms%jri(n),0,n,num,1)*rhsp(:atoms%jri(n),0,n,num,1) *fpi_const
          CALL intgr3(rh,atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),pdis(0,n,num))
          pdis(4,n,num)=pdis(0,n,num)
          DO lh = 1,sphhar%nlh(sym%ntypsy(na))
             rh(:atoms%jri(n))=rh(:atoms%jri(n))+rhsp(:atoms%jri(n),lh,n,num,1)*&
                  rhsp(:atoms%jri(n),lh,n,num,1)* atoms%rmsh(:atoms%jri(n),n)*atoms%rmsh(:atoms%jri(n),n)
             CALL intgr3(rh,atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),rhs)
             IF (lh<4) pdis(lh,n,num)=rhs
             pdis(4,n,num)=pdis(4,n,num)+rhs
          ENDDO
          pdis(:,n,num) = pdis(:,n,num)*atoms%neq(n)
          dis(num)=dis(num)+SUM(pdis(:,n,num))
          pdis(:,n,num) = SQRT(pdis(:,n,num)/atoms%volmts(n))*1000.
          na = na + atoms%neq(n)
       ENDDO
       !     ----> interstitial part
       !         ---> create density in the real space
       i = 0
       DO  i3 = 0,2*stars%mx3
          k3 = i3
          IF (k3.GT.stars%mx3) k3 = k3 - 2*stars%mx3-1
          DO  i2 = 0,2*stars%mx2
             k2 = i2
             IF (k2.GT.stars%mx2) k2 = k2 - 2*stars%mx2-1
             DO  i1 = 0,2*stars%mx1
                k1 = i1
                IF (k1.GT.stars%mx1) k1 = k1 - 2*stars%mx3-1
                i = i + 1
                id3 = stars%ig(k1,k2,k3)
                phase = stars%rgphs(k1,k2,k3)
                IF (id3.EQ.0) THEN
                   af3(i,1) = 0.
                   bf3(i,1) = 0.
                ELSE
                   af3(i,1) =  REAL(rhpw(id3,num,1))*REAL(phase) &
                           - AIMAG(rhpw(id3,num,1))*AIMAG(phase)
                   bf3(i,1) = AIMAG(rhpw(id3,num,1))*REAL(phase) &
                           +  REAL(rhpw(id3,num,1))*AIMAG(phase)
                END IF
             ENDDO
          ENDDO
       ENDDO

       CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx1*2+1,stars%mx1*2+1,1)
       CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx2*2+1,nk12,1)
       CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx3*2+1,nt,1)
       !         ---> form the dot product
       i = 0
       DO  i3 = 0,stars%mx3*2
          DO  i2 = 0,stars%mx2*2
             DO  i1 = 0,stars%mx3*2
                i = i + 1
                af3(i,1) = af3(i,1)*af3(i,1)
                bf3(i,1) = 0.
             ENDDO
          ENDDO
       ENDDO
       !         ---> back fft
       CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx1*2+1,stars%mx1*2+1,-1)
       CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx2*2+1,nk12,-1)
       CALL cfft(af3(1,1),bf3(1,1),nt,stars%mx3*2+1,nt,-1)
       sumis = 0.
       i = 0
       DO  i3 = 0,stars%mx3*2
          k3 = i3
          IF (k3.GT.stars%mx3) k3 = k3 - stars%mx3*2-1
          DO  i2 = 0,stars%mx2*2
             k2 = i2
             IF (k2.GT.stars%mx2) k2 = k2 - stars%mx2*2-1
             DO  i1 = 0,stars%mx1*2
                k1 = i1
                IF (k1.GT.stars%mx1) k1 = k1 - stars%mx1*2-1
                i = i + 1
                phase = stars%rgphs(k1,k2,k3)
                id3 = stars%ig(k1,k2,k3)
                IF (id3.NE.0) THEN
                   sumis = sumis + REAL(CMPLX(af3(i,1),bf3(i,1))* CONJG(stars%ustep(id3)))*phase*fact
                END IF
             ENDDO
          ENDDO
       ENDDO
       pdis(0,0,num)=SQRT(sumis/cell%volint)*1000.
       dis(num) = dis(num) + sumis
       IF (input%film) THEN
          !     ----> vacuum part
          DO  ivac = 1,vacuum%nvac
             DO  ip = 1,vacuum%nmzxy
                !         ---> create density in the real space
                i = 0
                DO  i2 = 0,stars%mx2*2
                   k2 = i2
                   IF (k2.GT.stars%mx2) k2 = k2 - stars%mx2*2-1
                   DO  i1 = 0,stars%mx1*2
                      k1 = i1
                      IF (k1.GT.stars%mx1) k1 = k1 - stars%mx1*2-1
                      i = i + 1
                      id3 = stars%ig(k1,k2,0)
                      IF (id3.EQ.0) THEN
                         af2(i,1) = 0.
                         bf2(i,1) = 0.
                      ELSE
                         id2 = stars%ig2(id3)
                         phase = stars%rgphs(k1,k2,0)
                         IF (id2.EQ.1) THEN
                            af2(i,1) = rhv0(ip,ivac,num,1)
                            bf2(i,1) = 0.
                         ELSE
                            af2(i,1) = REAL(rhv1(ip,id2-1,ivac,num,1))
                            bf2(i,1) =AIMAG(rhv1(ip,id2-1,ivac,num,1))
                         END IF
                      END IF
                   ENDDO
                ENDDO
                CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx1*2+1,stars%mx1*2+1,1)
                CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx2*2+1,nk12,1)
                !         ---> form dot product
                i = 0
                DO  i2 = 0,stars%mx2*2
                   DO  i1 = 0,stars%mx1*2
                      i = i + 1
                      af2(i,1) = af2(i,1)*af2(i,1)
                      bf2(i,1) = 0.
                   ENDDO
                ENDDO
                !         ---> back fft
                CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx1*2+1,stars%mx1*2+1,-1)
                CALL cfft(af2(1,1),bf2(1,1),nk12,stars%mx2*2+1,nk12,-1)
                disz(npz-ip) = cell%area*af2(1,1)/REAL(nk12)
             ENDDO
             !         ---> beyond warping region
             DO  ip = vacuum%nmzxy + 1,vacuum%nmz
                disz(npz-ip) = cell%area*rhv0(ip,ivac,num,1)
             ENDDO
             CALL intgz0(disz,vacuum%delz,vacuum%nmz,sumz,tail)
             IF (sym%zrfs .OR. sym%invs) facv = 2.0
             dis(num) = dis(num) + facv*sumz
          ENDDO
       END IF
    ENDDO

    dis(:) = SQRT(dis(:)/cell%vol)*1000.

    WRITE (oUnit,FMT=8000) iter,dis(1)
8000 FORMAT (/,'----> distance of  the potential for it=',i3,':',f11.6, ' mhtr/bohr**3')
    WRITE(oUnit,*) "Details of potential differences for each atom type"
    WRITE(oUnit,*) "Atom: total difference: difference of first sphhar"
    DO n=1,atoms%ntype
       WRITE(oUnit,"(i5,' : ',f10.6,' : ',4f10.6)") n,pdis(4,n,1),pdis(0:3,n,1)
    ENDDO
    WRITE(oUnit,*) "Difference of interstitial:",pdis(0,0,1)
    IF (input%jspins.EQ.2) THEN
       WRITE (oUnit,FMT=8010) iter,dis(2)
8010   FORMAT (/,'----> distance of spin potential for it=',i3,':', f11.6,' mhtr/bohr**3')
       WRITE(oUnit,*) "Details of potential differences for each atom type"
       WRITE(oUnit,*) "Atom: total difference: difference of first sphhar"
       DO n=1,atoms%ntype
          WRITE(oUnit,"(i5,' : ',f10.6,' : ',4f10.6)") n,pdis(4,n,2),pdis(0:3,n,2)
       ENDDO
       WRITE(oUnit,*) "Difference of interstitial:",pdis(0,0,2)
    END IF
  CONTAINS
    SUBROUTINE priv_cdndif(m1,m2,m3,n1,n2,n3,s1,s2,film,rhsp,rhpw,rhv0,rhv1)
      IMPLICIT NONE
      INTEGER, INTENT (IN) :: m1,m2,m3,n1,n2,n3
      REAL,    INTENT (IN) :: s1,s2
      COMPLEX, INTENT (INOUT) :: rhpw(:,:,:)
      COMPLEX, INTENT (INOUT) :: rhv1(:,:,:,:,:)
      REAL,    INTENT (INOUT) :: rhsp(:,:,:,:,:)
      REAL,    INTENT (INOUT) :: rhv0(:,:,:,:)
      LOGICAL,INTENT (IN)     :: film

      rhsp(:,:,:,m1,n1)=s1*rhsp(:,:,:,m2,n2) +  s2*rhsp(:,:,:,m3,n3)
      rhpw(:,m1,n1)    =s1*rhpw(:,m2,n2)     +  s2*rhpw(:,m3,n3)
      IF (film) THEN
         rhv0(:,:,m1,n1)  =s1*rhv0(:,:,m2,n2)   + s2*rhv0(:,:,m3,n3)
         rhv1(:,:,:,m1,n1)=s1*rhv1(:,:,:,m2,n2) + s2*rhv1(:,:,:,m3,n3)
      ENDIF
    END SUBROUTINE priv_cdndif

  END SUBROUTINE potdis
END MODULE m_potdis

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR