text:
      SUBROUTINE SETPCL (PCCTAB)

C  SETPCL: loads the pcctab table with the names of the picture
C  change categories.

      CHARACTER  PCCTAB(17)*34

      PCCTAB(01) = 'polyline bundle representation'
      PCCTAB(02) = 'polymarker bundle representation'
      PCCTAB(03) = 'text bundle representation'
      PCCTAB(04) = 'interior bundle representation'
      PCCTAB(05) = 'edge bundle representation'
      PCCTAB(06) = 'pattern representation'
      PCCTAB(07) = 'colour representation'
      PCCTAB(08) = 'view representation'
      PCCTAB(09) = 'workstation transformation'
      PCCTAB(10) = 'highlighting filter'
      PCCTAB(11) = 'invisibility filter'
      PCCTAB(12) = 'HLHSR mode'
      PCCTAB(13) = 'structure content'
      PCCTAB(14) = 'post structure'
      PCCTAB(15) = 'unpost structure'
      PCCTAB(16) = 'delete structure'
      PCCTAB(17) = 'reference modification'

      END


text:
	SUBROUTINE SETVALNB ( LUNIT, TAGPV, NTAGPV, TAGNB, NTAGNB,
     .                        R8VAL, IRET )

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    SETVALNB
C   PRGMMR: J. ATOR          ORG: NCEP       DATE: 2016-07-29
C
C ABSTRACT:  THIS SUBROUTINE SHOULD ONLY BE CALLED WHEN A BUFR FILE IS
C   OPENED FOR OUTPUT, AND A SUBSET DEFINITION MUST ALREADY BE IN SCOPE
C   VIA A PREVIOUS CALL TO BUFR ARCHIVE LIBRARY SUBROUTINE OPENMB OR
C   EQUIVALENT.  THE FUNCTION WILL FIRST SEARCH FOR THE (NTAGPV)th
C   OCCURRENCE OF MNEMONIC TAGPV WITHIN THE OVERALL SUBSET DEFINITION,
C   COUNTING FROM THE BEGINNING OF THE SUBSET.  IF FOUND, IT WILL THEN
C   SEARCH FORWARD (IF NTAGNB IS POSITIVE) OR BACKWARD (IF NTAGNB IS
C   NEGATIVE) FROM THAT POINT WITHIN THE SUBSET FOR THE (NTAGNB)th
C   OCCURRENCE OF MNEMONIC TAGNB AND STORE R8VAL AS THE VALUE
C   CORRESPONDING TO THAT MNEMONIC.
C
C PROGRAM HISTORY LOG:
C 2016-07-29  J. ATOR    -- ORIGINAL AUTHOR; BASED ON GETVALNB
C
C USAGE:    CALL SETVALNB (LUNIT, TAGPV, NTAGPV, TAGNB, NTAGNB,
C                          R8VAL, IRET)
C   INPUT ARGUMENT LIST:
C     LUNIT    - INTEGER: FORTRAN LOGICAL UNIT NUMBER FOR BUFR FILE
C     TAGPV    - CHARACTER*(*): PIVOT MNEMONIC; THE FUNCTION WILL
C                FIRST SEARCH FOR the (NTAGPV)th OCCURRENCE OF THIS
C                MNEMONIC, COUNTING FROM THE BEGINNING OF THE OVERALL
C                SUBSET DEFINITION
C     NTAGPV   - INTEGER: ORDINAL OCCURRENCE OF TAGPV TO SEARCH FOR
C     TAGNB    - CHARACTER*(*): NEARBY MNEMONIC; ASSUMING TAGPV IS
C                SUCCESSFULLY FOUND, THE FUNCTION WILL THEN SEARCH
C                NEARBY FOR THE (NTAGNB)th OCCURRENCE OF TAGNB AND
C                STORE R8VAL AS THE CORRESPONDING VALUE
C     NTAGNB   - INTEGER: ORDINAL OCCURRENCE OF TAGNB TO SEARCH FOR,
C                COUNTING FROM THE LOCATION OF TAGPV WITHIN THE OVERALL
C                SUBSET DEFINITION.  IF NTAGNB IS POSITIVE, THE FUNCTION
C                WILL SEARCH IN A FORWARD DIRECTION FROM THE LOCATION OF
C                TAGPV, OR IF NTAGNB IS NEGATIVE IT WILL INSTEAD SEARCH
C                IN A BACKWARDS DIRECTION.
C     R8VAL    - REAL*8: VALUE TO BE STORED CORRESPONDING TO (NTAGNB)th
C                OCCURRENCE OF TAGNB.
C
C   OUTPUT ARGUMENT LIST:
C     IRET     - INTEGER: RETURN CODE
C                   0 = NORMAL RETURN
C                  -1 = (NTAGNB)th OCCURENCE OF MNEMONIC TAGNB COULD
C                       NOT BE FOUND, OR SOME OTHER ERROR OCCURRED
C
C REMARKS:
C    THIS ROUTINE CALLS:        FSTAG    STATUS
C    THIS ROUTINE IS CALLED BY: None
C                               Normally called only by application
C                               programs
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

	USE MODA_USRINT
	USE MODA_MSGCWD
	USE MODA_TABLES

	INCLUDE 'bufrlib.prm'

	CHARACTER*(*) TAGPV, TAGNB

	REAL*8  R8VAL

C----------------------------------------------------------------------
C----------------------------------------------------------------------

	IRET = -1

C	Get LUN from LUNIT.

	CALL STATUS (LUNIT, LUN, IL, IM )
	IF ( IL .LE. 0 ) RETURN
	IF ( INODE(LUN) .NE. INV(1,LUN) ) RETURN

C	Starting from the beginning of the subset, locate the (NTAGPV)th
C	occurrence of TAGPV.

	CALL FSTAG( LUN, TAGPV, NTAGPV, 1, NPV, IERFT )
	IF ( IERFT .NE. 0 ) RETURN

C	Now, starting from the (NTAGPV)th occurrence of TAGPV, search
C	forward or backward for the (NTAGNB)th occurrence of TAGNB.

	CALL FSTAG( LUN, TAGNB, NTAGNB, NPV, NNB, IERFT )
	IF ( IERFT .NE. 0 ) RETURN

	IRET = 0
	VAL(NNB,LUN) = R8VAL

	RETURN
	END

ext: f


!==============================================================================!
  subroutine Interface_Mod_Create(inter, grid, n_dom)
!------------------------------------------------------------------------------!
!   Create interface between two grids.                                        !
!                                                                              !
!   The subroutine works in the following way:                                 !
!   1. It allocates the entire interface on all processors.                    !
!   2. For both sides of the interface (1 and 2) it stores face coordinates    !
!      (xf_1 to zf_2), cell inside (ic_1, ic_2), bundary cells (ib_1, ib_2).   !
!      and processor ids (ip_1, ip_2).                                         !
!   3. While doint the step above, it takes care of the offsets of each        !
!      processor (off_1 and off_2).                                            !
!   4. Then takes global summs of all the arrays mentioned in step 2.          !
!   5. Sorts the above arrays by their x, y and z coordinates, carrying        !
!      information on cells inside (ic_1, ic), boundary cells (ib_1, ib_2)     !
!      and processor ids (ip_1, ip_2)                                          !
!   6. The sorting from step 5 gives mapping of local cells inside, and local  !
!      boundary cells to the global interface.                                 !
!   7. The mapping information is finally stored in arrays cell_1, cell_2      !
!      (for inside cells), bcel_1, bcel_2 for boundary cells and face_1        !
!      and face_2 for global faces at the interface.                           !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use Work_Mod, only: xf_1 => r_face_01,  &  ! face coordinates ...
                      yf_1 => r_face_02,  &  ! ... on the each side ...
                      zf_1 => r_face_03,  &  ! ... of the interface
                      xf_2 => r_face_04,  &
                      yf_2 => r_face_05,  &
                      zf_2 => r_face_06,  &
                      ic_1 => i_face_01,  &  ! internal cells in domain 1
                      ib_1 => i_face_02,  &  ! boundary cells in domain 1
                      ip_1 => i_face_03,  &  ! processor in domain 1
                      ic_2 => i_face_04,  &  ! internal cells in domain 2
                      ib_2 => i_face_05,  &  ! boundary cells in domain 2
                      ip_2 => i_face_06      ! processor in domain 2
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Interface_Type)    :: inter(MD, MD)
  type(Grid_Type), target :: grid(MD)
  integer                 :: n_dom
!-----------------------------------[Locals]-----------------------------------!
  integer, allocatable :: off_1(:)
  integer, allocatable :: off_2(:)
  integer              :: n, n1, n2, n1_tot, n2_tot, n_tot
  integer              :: pos, c, c1, c2, d1, d2, k, nks, s, p
  logical              :: found
  character(len=80)    :: keys(128)
!==============================================================================!

  ! Allocate memory for offsets
  if(n_proc > 1) then
    allocate(off_1(n_proc)); off_1 = 0
    allocate(off_2(n_proc)); off_2 = 0
  end if

  !--------------------------------------------!
  !                                            !
  !   Browse through all domain combinations   !
  !                                            !
  !--------------------------------------------!
  do d1 = 1, n_dom
    do d2 = 1, n_dom

      ! Innitialize number of faces at the interface
      inter(d1, d2) % n_tot = 0

      ! Store pointers to the grids surrounding the interface
      inter(d1, d2) % pnt_grid1 => grid(d1)
      inter(d1, d2) % pnt_grid2 => grid(d2)

      if(d1 .ne. d2) then

        xf_1(:) = 0;  yf_1(:) = 0; zf_1(:) = 0
        xf_2(:) = 0;  yf_2(:) = 0; zf_2(:) = 0
        ic_1(:) = 0;  ib_1(:) = 0; ip_1(:) = 0
        ic_2(:) = 0;  ib_2(:) = 0; ip_2(:) = 0

        !-----------------------------------------------------------------!
        !   Try to find specified interface condition between d1 and d2   !
        !-----------------------------------------------------------------!
        call Control_Mod_Position_At_Three_Keys('INTERFACE_CONDITION',    &
                                                trim(problem_name(d1)),   &
                                                trim(problem_name(d2)),   &
                                                found,                    &
                                                verbose=.false.)

        !---------------------------------------------------!
        !   Found specification between domains d1 and d2   !
        !---------------------------------------------------!
        if(found) then
          call Control_Mod_Read_Strings_On('BOUNDARY_CONDITIONS',  &
                                           keys, nks, .false.)
          do k = 1, nks
            call To_Upper_Case(keys(k))
          end do

          !---------------------------------------!
          !   Count boundary cells on each side   !
          !---------------------------------------!

          ! Innitialize counter in first domain
          n1 = 0

          ! On physical boundary cells
          do s = 1, grid(d1) % n_faces
            c2 = grid(d1) % faces_c(2,s)
            if(c2 < 0) then
              if(Grid_Mod_Bnd_Cond_Name(grid(d1), c2) .eq. keys(1)) then
                n1 = n1 + 1
              end if
            end if
          end do

          ! On periodic faces of domain 1
          do s = 1, grid(d1) % n_faces
            c1 = grid(d1) % faces_c(1,s)
            if(Grid_Mod_Bnd_Cond_Name(grid(d1), s) .eq. keys(1) .and.  &
               grid(d1) % comm % cell_proc(c1) .eq. this_proc) then
              n1 = n1 + 1
            end if
          end do
          inter(d1, d2) % n1_sub = n1

          ! Innitialize counter in first domain
          n2 = 0

          ! On physical boundary cells
          do s = 1, grid(d2) % n_faces
            c2 = grid(d2) % faces_c(2,s)
            if(c2 < 0) then
              if(Grid_Mod_Bnd_Cond_Name(grid(d2), c2) .eq. keys(2)) then
                n2 = n2 + 1
              end if
            end if
          end do

          ! On periodic faces
          do s = 1, grid(d2) % n_faces
            c1 = grid(d2) % faces_c(1,s)
            if(Grid_Mod_Bnd_Cond_Name(grid(d2), s) .eq. keys(2) .and.  &
               grid(d2) % comm % cell_proc(c1) .eq. this_proc) then
              n2 = n2 + 1
            end if
          end do
          inter(d1, d2) % n2_sub = n2

          !------------------------------------------!
          !     Check if the mapping is conformal    !
          !   (Simply by comparing number of faces   !
          !      on each side of the interface)      !
          !------------------------------------------!
          n1_tot = n1
          n2_tot = n2
          call Comm_Mod_Global_Sum_Int(n1_tot)
          call Comm_Mod_Global_Sum_Int(n2_tot)

          if(n1_tot .ne. n2_tot) then
            if(this_proc < 2) then
              print *, '# Number of cells at the interface between ',  &
                        trim(problem_name(d1)), ' and ',               &
                        trim(problem_name(d2)), ' is not the same!'
              print *, '# Only conformal mappings are supported.  Exiting!'
            end if
            call Comm_Mod_End()
          else
            n_tot = n1_tot
            inter(d1, d2) % n_tot = n_tot
            if(this_proc < 2) then
              print '(5a,i6,a)', ' # Domains ', trim(problem_name(d1)),  &
                                       ' and ', trim(problem_name(d2)),  &
                        ' are connected with ', n_tot, ' interface cells!'
            end if
          end if

          !-----------------------------------------------------!
          !          Store boundary cells on each side          !
          !   (Inside cells in case periodicity is specified)   !
          !-----------------------------------------------------!
          ic_1(1:n_tot) = 0; ic_2(1:n_tot) = 0
          ib_1(1:n_tot) = 0; ib_2(1:n_tot) = 0
          ip_1(1:n_tot) = 0; ip_2(1:n_tot) = 0
          allocate(inter(d1, d2) % phi_1(n_tot, MAX_VARS_INTERFACE))
          allocate(inter(d1, d2) % phi_2(n_tot, MAX_VARS_INTERFACE))
          inter(d1, d2) % phi_1(:,:) = 0.0
          inter(d1, d2) % phi_2(:,:) = 0.0

          ! Work out offsets for each interface for each processor
          if(n_proc > 1) then
            off_1(1:n_proc)  = 0
            off_2(1:n_proc)  = 0
            off_1(this_proc) = n1
            off_2(this_proc) = n2
            call Comm_Mod_Global_Sum_Int_Array(n_proc, off_1)
            call Comm_Mod_Global_Sum_Int_Array(n_proc, off_2)

            do p = n_proc, 2, -1
              off_1(p) = sum(off_1(1:p-1))
              off_2(p) = sum(off_2(1:p-1))
            end do
            off_1(1) = 0
            off_2(1) = 0
          end if

          ! On the side of d1
          n1 = 0
          do s = 1, grid(d1) % n_faces
            c1 = grid(d1) % faces_c(1,s)
            c2 = grid(d1) % faces_c(2,s)
            if(c2 < 0) then
              if(Grid_Mod_Bnd_Cond_Name(grid(d1), c2) .eq. keys(1)) then
                n1 = n1 + 1
                pos = n1
                if(n_proc > 1) pos = pos + off_1(this_proc)
                ic_1(pos) = c1
                ib_1(pos) = c2
                ip_1(pos) = this_proc
                xf_1(pos) = grid(d1) % xf(s)
                yf_1(pos) = grid(d1) % yf(s)
                zf_1(pos) = grid(d1) % zf(s)
              end if
            end if
          end do
          do s = 1, grid(d1) % n_faces
            c1 = grid(d1) % faces_c(1,s)
            c2 = grid(d1) % faces_c(2,s)
            if(Grid_Mod_Bnd_Cond_Name(grid(d1), s) .eq. keys(1) .and.  &
               grid(d1) % comm % cell_proc(c1) .eq. this_proc) then
              n1 = n1 + 1
              pos = n1
              if(n_proc > 1) pos = pos + off_1(this_proc)
              ic_1(pos) = c1
              ib_1(pos) = c2
              ip_1(pos) = this_proc
              xf_1(pos) = grid(d1) % xf(s)
              yf_1(pos) = grid(d1) % yf(s)
              zf_1(pos) = grid(d1) % zf(s)
            end if
          end do

          ! On the side of d2
          n2 = 0
          do s = 1, grid(d2) % n_faces
            c1 = grid(d2) % faces_c(1,s)
            c2 = grid(d2) % faces_c(2,s)
            if(c2 < 0) then
              if(Grid_Mod_Bnd_Cond_Name(grid(d2), c2) .eq. keys(2)) then
                n2 = n2 + 1
                pos = n2
                if(n_proc > 1) pos = pos + off_2(this_proc)
                ic_2(pos) = c1
                ib_2(pos) = c2
                ip_2(pos) = this_proc
                xf_2 (pos) = grid(d2) % xf(s)
                yf_2 (pos) = grid(d2) % yf(s)
                zf_2 (pos) = grid(d2) % zf(s)
              end if
            end if
          end do
          do s = 1, grid(d2) % n_faces
            c1 = grid(d2) % faces_c(1,s)
            c2 = grid(d2) % faces_c(2,s)
            if(Grid_Mod_Bnd_Cond_Name(grid(d2), s) .eq. keys(2) .and.  &
               grid(d2) % comm % cell_proc(c1) .eq. this_proc) then
              n2 = n2 + 1
              pos = n2
              if(n_proc > 1) pos = pos + off_2(this_proc)
              ic_2(pos) = c1
              ib_2(pos) = c2
              ip_2(pos) = this_proc
              xf_2 (pos) = grid(d2) % xf(s)
              yf_2 (pos) = grid(d2) % yf(s)
              zf_2 (pos) = grid(d2) % zf(s)
            end if
          end do

          !----------------------------------------------------------!
          !   Distribute interface coordinates over all processors   !
          !----------------------------------------------------------!
          call Comm_Mod_Global_Sum_Real_Array(n_tot, xf_1(1:n_tot))
          call Comm_Mod_Global_Sum_Real_Array(n_tot, yf_1(1:n_tot))
          call Comm_Mod_Global_Sum_Real_Array(n_tot, zf_1(1:n_tot))
          call Comm_Mod_Global_Sum_Real_Array(n_tot, xf_2(1:n_tot))
          call Comm_Mod_Global_Sum_Real_Array(n_tot, yf_2(1:n_tot))
          call Comm_Mod_Global_Sum_Real_Array(n_tot, zf_2(1:n_tot))
          call Comm_Mod_Global_Sum_Int_Array (n_tot, ic_1(1:n_tot))
          call Comm_Mod_Global_Sum_Int_Array (n_tot, ib_1(1:n_tot))
          call Comm_Mod_Global_Sum_Int_Array (n_tot, ip_1(1:n_tot))
          call Comm_Mod_Global_Sum_Int_Array (n_tot, ic_2(1:n_tot))
          call Comm_Mod_Global_Sum_Int_Array (n_tot, ib_2(1:n_tot))
          call Comm_Mod_Global_Sum_Int_Array (n_tot, ip_2(1:n_tot))

          ! Sort interfaces from domain 1 carrying
          ! information of cells surrounding it along
          call Sort_Mod_3_Real_Carry_3_Int(xf_1(1:n_tot),  &
                                           yf_1(1:n_tot),  &
                                           zf_1(1:n_tot),  &
                                           ic_1(1:n_tot),  &
                                           ib_1(1:n_tot),  &
                                           ip_1(1:n_tot))

          ! Sort interfaces from domain 2 carrying
          ! information of cells surrounding it along
          call Sort_Mod_3_Real_Carry_3_Int(xf_2(1:n_tot),  &
                                           yf_2(1:n_tot),  &
                                           zf_2(1:n_tot),  &
                                           ic_2(1:n_tot),  &
                                           ib_2(1:n_tot),  &
                                           ip_2(1:n_tot))
!         ! Write some debugging information
!         do n = 1, inter(d1, d2) % n_tot
!           write(100*this_proc+d1*10+d2, '(i4, 6f10.5)')   &
!                     n,                                    &
!                     xf_1(n), yf_1(n), zf_1(n),            &
!                     xf_2(n), yf_2(n), zf_2(n)
!         end do

          !------------------------------!
          !   Store buffer information   !
          !------------------------------!
          allocate(inter(d1, d2) % cell_1(inter(d1, d2) % n1_sub))
          allocate(inter(d1, d2) % face_1(inter(d1, d2) % n1_sub))
          allocate(inter(d1, d2) % bcel_1(inter(d1, d2) % n1_sub))
          allocate(inter(d1, d2) % cell_2(inter(d1, d2) % n2_sub))
          allocate(inter(d1, d2) % face_2(inter(d1, d2) % n2_sub))
          allocate(inter(d1, d2) % bcel_2(inter(d1, d2) % n2_sub))
          inter(d1, d2) % cell_1 = 0
          inter(d1, d2) % face_1 = 0
          inter(d1, d2) % bcel_1 = 0
          inter(d1, d2) % cell_2 = 0
          inter(d1, d2) % face_2 = 0
          inter(d1, d2) % bcel_2 = 0

          n1 = 0
          n2 = 0
          do n = 1, inter(d1, d2) % n_tot  ! browse through interface now

            ! Handle domain 1
            if(ip_1(n) .eq. this_proc) then
              n1 = n1 + 1
              do c = 1, grid(d1) % n_cells
                if(c .eq. ic_1(n)) then
                  inter(d1, d2) % cell_1(n1) = c
                  inter(d1, d2) % face_1(n1) = n
                  inter(d1, d2) % bcel_1(n1) = ib_1(n)
                  goto 1
                end if
              end do
1             continue
            end if

            ! Handle domain 2
            if(ip_2(n) .eq. this_proc) then
              n2 = n2 + 1
              do c = 1, grid(d2) % n_cells
                if(c .eq. ic_2(n)) then
                  inter(d1, d2) % cell_2(n2) = c
                  inter(d1, d2) % face_2(n2) = n
                  inter(d1, d2) % bcel_2(n2) = ib_2(n)
                  goto 2
                end if
              end do
2             continue
            end if

          end do  ! next face in the interface

        end if  ! if this interface is found in control mod
      end if

    end do
  end do

  if(n_proc > 1) then
    deallocate(off_1)
    deallocate(off_2)
  end if

  end subroutine


text:
      subroutine phi1_2m(m2,x3,xth,xphi,s3min,p1,p2,p3,wt,*)
c     massive particle p1 decaying into p2 mass m2 and p3 mass-squared s3.
c     with invariant mass of particle three s3 integrated over.
c     s3min is the minimum value of s3.
c     Vectors returned p2 and p3 are in the same frame as p1 is supplied.
c     Expression evaluated is
c     ds3 d^4 p2 d^4 p3 (2 pi)^4 delta(p1-p2-p3)/(2 pi)^6
c     delta(p2^2-m2) delta(p3^2-s3)
      implicit none
      include 'constants.f'
      include 'debug.f'
      include 'zerowidth.f'
      include 'breit.f'
      double precision p1(4),p2(4),p3(4),p3cm(4)
      double precision x3,xth,xphi,costh,sinth,phi,cphi,sphi
      double precision wt,wt0,w3
      double precision s3max,s3min,xx
      double precision m1,m2,m3,s1,s2,s3,lambda
      integer j
      parameter(wt0=one/8d0/pi)

      wt=0d0
      s1=p1(4)**2-p1(1)**2-p1(2)**2-p1(3)**2
      if (s1 .lt. 0d0) return 1
      m1=dsqrt(s1)
      s2=m2**2
      s3max=(m2-m1)**2
      if (s3min .gt. s3max) return 1
      if (n3 .eq. 0) then
         w3=s3max-s3min
         s3=s3max*x3+s3min*(1d0-x3)
         xx=0d0
      elseif (n3 .eq. 1) then
        if ((zerowidth) .and. (s3max .lt. mass3)) return 1
        xx=1d0
        call breitw(x3,s3min,s3max,mass3,width3,s3,w3)
      endif

      m3=dsqrt(s3)
      if (m1-m2-m3.lt. 0d0) return 1

      costh=two*xth-one
      phi=twopi*xphi
      sinth=dsqrt(one-costh**2)
      cphi=dcos(phi)
      sphi=dsin(phi)
      lambda=((s1-s2-s3)**2-4d0*s2*s3)

      if ((m1-m2-m3.lt. 0d0) .or. debug) then
      write(6,*) 'lambda in phi1_2m',lambda
      write(6,*) 's1 in phi1_2m',s1
      write(6,*) 's2 in phi1_2m',s2
      write(6,*) 's3 in phi1_2m',s3
      write(6,*) 'm1 in phi1_2m',m1
      write(6,*) 'm2 in phi1_2m',m2
      write(6,*) 'm3 in phi1_2m',m3
      write(6,*) 'm1-m2-m3 in phi1_2m',m1-m2-m3
      write(6,*) 'xx in phi1_2m',xx
      write(6,*) 'x3 in phi1_2m',x3
      write(6,*) 'n3 in phi1_2m',n3
      write(6,*) 'mass3 in phi1_2m',mass3
      return 1
      endif
      lambda=dsqrt(lambda)

      wt=wt0*w3*lambda/s1

      if(debug) write(6,*) 'wt in phi1_2m',wt

      p3cm(4)=m1/two*(s1+s3-s2)/s1
      p3cm(1)=m1/two*lambda/s1*sinth*sphi
      p3cm(2)=m1/two*lambda/s1*sinth*cphi
      p3cm(3)=m1/two*lambda/s1*costh
      call boost(m1,p1,p3cm,p3)
      do j=1,4
      p2(j)=p1(j)-p3(j)
      enddo

      if (  (p1(4) .lt. 0d0)
     & .or. (p2(4) .lt. 0d0)
     & .or. (p3(4) .lt. 0d0)) then
      return 1
      endif
      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
module kinds_m
	implicit none
	integer, parameter :: MP = kind(0)
	integer, parameter :: SP = kind(0.0e0)
	integer, parameter :: DP = kind(0.0d0)

	integer, parameter :: I2 = 2
	integer, parameter :: I4 = 4
	integer, parameter :: I8 = 8
end module

max_stars_repo_name: juselius/libgetkw

max_issues_repo_name: juselius/libgetkw

max_forks_repo_name: juselius/libgetkw


text:
!>\file module_nst_parameters.f90
!>\defgroup nst_parameters GFS NSST Parameter Module
!! \ingroup GFS_NSST
!!\brief This module includes constants and parameters used in GFS
!! near sea surface temperature scheme.
!! @{
module module_nst_parameters
  use machine, only :  kind_phys &
       ,kind_rad ! for astronomy (date) calculations
  !
  ! air constants and coefficients from the atmospehric model
  use physcons, only: &
       eps =>  con_eps &
       ,cp_a => con_cp &          ! spec heat air @p    (j/kg/k)
       , epsm1 => con_epsm1 &
       , hvap => con_hvap &       ! lat heat h2o cond   (j/kg)
       ,sigma_r => con_sbc  &     ! stefan-boltzmann    (w/m2/k4)
       ,grav => con_g         &   ! acceleration due to gravity (kg/m/s^2)
       ,omega => con_omega    &    ! ang vel of earth    (1/s)
       ,rvrdm1 => con_fvirt &
       ,rd => con_rd &
       ,rocp => con_rocp  &        ! r/cp
       ,pi => con_pi
  !
  ! note: take timestep from here later
  public
  integer :: &
       niter_conv = 5, &
       niter_z_w  = 5, &
       niter_sfs  = 5
  real (kind=kind_phys), parameter :: &
       !
       ! general constants
        sec_in_day=86400.       &
       ,sec_in_hour=3600.       &
       ,solar_time_6am=21600.0  &
       ,const_rot=0.000073      &          !< constant to calculate corioli force
       ,ri_c=0.65               &
       ,ri_g=0.25               &
       ,eps_z_w=0.01            &          !< criteria to finish iterations for z_w
       ,eps_conv=0.01           &          !< criteria to finish iterations for d_conv
       ,eps_sfs=0.01            &          !< criteria to finish iterations for d_sfs
       ,z_w_max=30.0            &          !< max warm layer thickness
       ,z_w_min=0.2             &          !< min warm layer thickness
       ,z_w_ini=0.2             &          !< initial warm layer thickness in dtl_onset
       ,z_c_max=0.01            &          !< maximum of sub-layer thickness (m)
       ,z_c_ini=0.001           &          !< initial value of z_c
       ,ustar_a_min=0.031       &          !< minimum of friction wind speed (m/s): 0.031 ~ 1m/s at 10 m hight
       ,tau_min=0.005           &          !< minimum of wind stress for dtm
       ,exp_const=9.5           &          !< coefficient in exponet profile
       ,delz=0.1                &          !< vertical increment for integral calculation   (m)
       ,von=0.4                 &          !< von karman's "constant"
       ,t0k=273.16              &          !<  celsius to kelvin
       ,gray=0.97               &
       ,sst_max=308.16          &
       ,tw_max=5.0              &
       ,wd_max=2.0              &
       ,omg_m =1.0              &          !< trace factor to apply salinity effect
       ,omg_rot = 1.0           &          !< trace factor to apply rotation effect
       ,omg_sh = 1.0            &          !< trace factor to apply sensible heat due to rainfall effect
       ,visw=1.e-6 &                       !< m2/s kinematic viscosity water
       ,novalue=0 &
       ,smallnumber=1.e-6 &
       ,timestep_oc=sec_in_day/8. &        !< time step in the ocean model (3 hours)
       ,radian=2.*pi/180.       &
       ,rad2deg=180./pi       &
       ,cp_w=4000.   &                     !< specific heat water (j/kg/k )
       ,rho0_w=1022.0 &                    !< density water (kg/m3 ) (or 1024.438)
       ,vis_w=1.e-6  &                     !< kinematic viscosity water (m2/s )
       ,tc_w=0.6    &                      !< thermal conductivity water (w/m/k )
       ,capa_w =3950.0 &                   !< heat capacity of sea water      !
       ,thref =1.0e-3                      !< reference value of specific volume (m**3/kg)

!!$!============================================
!!$
!!$  ,lvapor=2.453e6 &        ! latent heat of vaporization note: make it function of t ????? note the same as hvap
!!$       ,alpha=1 ! thermal expansion coefficient
!!$  ,beta ! saline contraction coefficient
!!$  ,cp=1 !=1 specific heat of sea water
!!$  ,g=1 ! acceleration due to gravity
!!$  ,kw=1 ! thermal conductivity of water
!!$  ,nu=1    !kinematic wiscosity
!!$  ,rho_w=1 !water density
!!$  ,rho_a=1 !air density
!!$  ,l_vapr=2.453e6
!!$  ,novalue=--1.0e+10
!!$
!!$c factors
!!$      beta=1.2     !given as 1.25 in fairall et al.(1996)
!!$      von=0.4      ! von karman's "constant"
!!$c      fdg=1.00     ! fairall's lkb rr to von karman adjustment
!!$      fdg=1.00     !based on results from flux workshop august 1995
!!$      tok=273.16   ! celsius to kelvin
!!$      twopi=3.14159*2.
!!$
!!$c air constants and coefficients
!!$      rgas=287.1                  !j/kg/k     gas const. dry air
!!$      xlv=(2.501-0.00237*ts)*1e+6  !j/kg  latent heat of vaporization at ts
!!$      cpa=1004.67                 !j/kg/k specific heat of dry air (businger 1982)
!!$      cpv=cpa*(1+0.84*q)          !moist air - currently not used (businger 1982)
!!$      rhoa=p*100./(rgas*(t+tok)*(1.+.61*q)) !kg/m3  moist air density ( " )
!!$      visa=1.326e-5*(1+6.542e-3*t+8.301e-6*t*t-4.84e-9*t*t*t)   !m2/s
!!$          !kinematic viscosity of dry air - andreas (1989) crrel rep. 89-11
!!$c
!!$c cool skin constants
!!$      al=2.1e-5*(ts+3.2)**0.79     !water thermal expansion coefft.
!!$      be=0.026                     !salinity expansion coefft.
!!$      cpw=4000.                    !j/kg/k specific heat water
!!$      rhow=1022.                   !kg/m3  density water
!!$      visw=1.e-6                   !m2/s kinematic viscosity water
!!$      tcw=0.6                      !w/m/k   thermal conductivity water
!!$      bigc=16.*grav*cpw*(rhow*visw)**3/(tcw*tcw*rhoa*rhoa)
!!$      wetc=0.622*xlv*qs/(rgas*(ts+tok)**2) !correction for dq;slope of sat. vap.
!!$
!!$!
!!$! functions
!!$
!!$
!!$  real,    parameter :: timestep=86400.    !integration time step, second
!!$
!!$  real, parameter    :: grav =9.81         !gravity, kg/m/s^2
!!$  real, parameter    :: capa =3950.0       !heat capacity of sea water
!!$  real, parameter    :: rhoref = 1024.438  !sea water reference density, kg/m^3
!!$  real   , parameter :: hslab=50.0         !slab ocean depth
!!$  real   , parameter :: bad=-1.0e+10
!!$  real   , parameter :: tmin=2.68e+02
!!$  real   , parameter :: tmax=3.11e+02
!!$
!!$  real, parameter :: grav =9.81           !gravity, kg/m/s^2
!!$  real, parameter :: capa =3950.0         !heat capacity of sea water
!!$  real, parameter :: rhoref = 1024.438    !sea water reference density, kg/m^3
!!$  real, parameter :: tmin=2.68e+02        !normal minimal temp
!!$  real, parameter :: tmax=3.11e+02        !normal max temp
!!$  real, parameter :: smin=1.0             !normal minimal salt
!!$  real, parameter :: smax=50.             !normal maximum salt
!!$  real, parameter :: visct=1.e-5          !viscocity for temperature diffusion
!!$  real, parameter :: viscs=1.e-5          !viscocity for salt diffusion
!!$
!!$
end module module_nst_parameters
!! @}

max_stars_repo_name: cacraigucar/ccpp-physics

max_issues_repo_name: cacraigucar/ccpp-physics

max_forks_repo_name: cacraigucar/ccpp-physics


!-------------------------------------------------------------------------------
!@(#) setting color pairs using ncurses(3c) from Fortran
!(LICENSE:PD)
program color_set_example
    use M_ncurses
    implicit none
    integer(C_SHORT) :: f,b
    integer          :: i,ierr
    stdscr=initscr()
    ierr=cbreak()
    ierr=noecho()
    if (has_colors()) then
        ierr=start_color()
        ierr=pair_content(0_C_SHORT,f,b)
        ierr=printw("pair 0 contains (%d,%d)"//C_NEW_LINE//C_NULL_CHAR,int(f),int(b))
        ierr=getch()
        ierr=printw("Initializing pair 1 to red/black"//C_NEW_LINE//C_NULL_CHAR)
        ierr=init_pair(1_C_SHORT, int(COLOR_RED,C_SHORT), int(COLOR_BLACK,C_SHORT))
        i = color_set(1_C_SHORT, C_NULL_PTR)
        ierr=printw("RED/BLACK (%s)"//C_NEW_LINE//C_NULL_CHAR, SHOW(i)//C_NULL_CHAR)
        ierr=getch()
        ierr=printw("Initializing pair 2 to white/blue"//C_NEW_LINE//C_NULL_CHAR)
        ierr=init_pair(2_C_SHORT, int(COLOR_WHITE,C_SHORT), int(COLOR_BLUE,C_SHORT))
        i = color_set(2_C_SHORT, C_NULL_PTR)
        ierr=printw("WHITE/BLUE (%s)"//C_NEW_LINE//C_NULL_CHAR, SHOW(i)//C_NULL_CHAR)
        ierr=getch()
        ierr=printw("Resetting colors to pair 0"//C_NEW_LINE//C_NULL_CHAR)
        i = color_set(0_C_SHORT, C_NULL_PTR)
        ierr=printw("Default Colors (%s)"//C_NEW_LINE//C_NULL_CHAR, SHOW(i)//C_NULL_CHAR)
        ierr=getch()
        ierr=printw("Resetting colors to pair 1"//C_NEW_LINE//C_NULL_CHAR)
        i = color_set(1_C_SHORT, C_NULL_PTR)
        ierr=printw("RED/BLACK (%s)"//C_NEW_LINE//C_NULL_CHAR, SHOW(i)//C_NULL_CHAR)
        !ierr=getch()
    else
        ierr=printw("This demo requires a color terminal"//C_NEW_LINE//C_NULL_CHAR)
        ierr=getch()
    endif
    ierr=endwin()
    write(*,*)'GOT TO END'
contains
!-------------------------------------------------------------------------------
character(len=3) function SHOW(n)
   integer :: n
   if(n.eq.0)then
      SHOW="OK"
   else
      SHOW="ERR"
   endif
end function SHOW
!-------------------------------------------------------------------------------
end program color_set_example
!-------------------------------------------------------------------------------


text:
subroutine bash
  call system('ls -1t *.f90 | head -n 1 | xargs cat')
end subroutine

max_stars_repo_name: ashwinvis/code-golf

max_issues_repo_name: ashwinvis/code-golf

max_forks_repo_name: ashwinvis/code-golf


text:
program main

  !======= Inclusions ===========
  use iso_fortran_env, only: wp => real64

  use KuraNetwork
  use LoadNetwork
  use iso_varying_string, only: VARYING_STRING, var_str
  use DopriWrap
  use dopri5Interface
  !======= Declarations =========
  implicit none
  ! local variables
  real(wp) :: tstart                   ! initial time
  real(wp) :: tend                     ! final time
  real(wp) :: rtol, atol               ! relative and absolute tolerance
  real(wp) :: dtout                    ! output time interval
  real(wp) :: tout                     ! output time
  real(wp) :: tcur                  ! current time
  integer :: nout                     ! number of outputs
  integer :: outstep                  ! output loop counter
  integer :: neq                  ! output loop counter
  type(dopri5_t)    :: prop
  logical :: status_ok
  integer :: idid
  integer,dimension(:), allocatable  :: icomp
  character(len=:), allocatable :: filename,buffer
  integer   :: testcase
  integer, dimension(:),allocatable  :: n,c
  character, dimension(:),allocatable:: nt
  character, dimension(:),allocatable:: cd
  type(VARYING_STRING) :: database
  real(wp)                          :: start_time, stop_time
  real(wp)                          :: rate,theta
  real(wp), dimension(:), allocatable :: y, time
  real(wp), dimension(:,:), allocatable :: output
  integer                           :: start_clock, stop_clock,i,j
  type(system_table)                :: stable
  real(wp) :: tload, tintegrate, ttotal
  integer                           :: refeid
  real(wp) :: evref

    allocate(character(512) :: filename)
    CALL get_environment_variable("KURABENCH_DB", filename)
    filename = trim(filename)

    if(command_argument_count().gt.0) then
        allocate(character(5) :: buffer)
        call get_command_argument(1,buffer)
        buffer = trim(buffer)
        read(buffer,'(I4.1)') testcase
        deallocate(buffer)
    else
        testcase=1
    end if
        tstart=0.0_wp
    if(command_argument_count().gt.1) then
        allocate(character(5) :: buffer)
        call get_command_argument(2,buffer)
        buffer = trim(buffer)
        read(buffer,*) tend
        deallocate(buffer)
    else
        tend=1.0_wp
    end if
    if(command_argument_count().gt.2) then
        allocate(character(10) :: buffer)
        call get_command_argument(3,buffer)
        buffer = trim(buffer)
        read(buffer,*) atol
        deallocate(buffer)
    else
        atol=1.0d-6
    end if
    if(command_argument_count().gt.3) then
        allocate(character(10) :: buffer)
        call get_command_argument(4,buffer)
        buffer = trim(buffer)
        read(buffer,*) rtol
        deallocate(buffer)
    else
        rtol=1.0d-6
    end if

  tcur   = tstart
  tout   = tstart
  dtout  = 2.0d0
  nout   = ceiling(tend/dtout)
  call cpu_time(start_time)
  call dopri_loadNet(filename,testcase)
  neq = ln_table%nnodes
  allocate(y(neq))
  allocate(icomp(neq))

  icomp = (/(i, i=1,neq, 1)/)
  call load_initial_conditions(y)
  call cpu_time(stop_time)
  tload = stop_time - start_time
  prop = initialize(fcn=dopri_RhsFn,n=neq,solout=solout,iout=0)

  call cpu_time(start_time)
     !now, perform the integration:
  call prop%integrate(tcur,y,tend,rtol,atol,idid=idid)

  call cpu_time(stop_time)
  tintegrate=stop_time - start_time
  !call write_experiment_states()
  call get_errorvref(refeid,evref,y)
  call write_experiment_header(atol, rtol, tstart, tend, &
                              tload, tintegrate, tload+tintegrate,&
                              refeid, evref)
  call write_final_state(y,tcur)
  call close_database_conn()
end program main

max_stars_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks

max_issues_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks

max_forks_repo_name: PIK-ICoNe/NetworkDynamicsBenchmarks


text:
doch jetzt ist der Held gefallen .
neue Modelle werden erprobt .
doch fehlen uns neue Ressourcen .

max_stars_repo_name: hthu/mt

max_issues_repo_name: hthu/mt

max_forks_repo_name: hthu/mt


text: Ergo sum, a journal of bioethics and the philosophy of science that is published by AS Papers.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!!
!> @brief      fundamental constants
!!
!! List of fundamental physical and mathematical constants. Unlike most other modules, all
!! parameters are public.
!!
!! When possible, constants contain more digits than double precision, so that
!! they are rounded correctly. Single letter constants contain underscore so
!! that they do not clash with user variables ("e" and "i" are frequently used as
!! loop variables)
!!
!! @author     Rodrigo Navarro Perez
!!
module constants
use precisions, only: dp
implicit none

public

real(dp), parameter :: pi = 3.1415926535897932384626433832795_dp !<  \f$\pi\f$
real(dp), parameter :: e_ = 2.7182818284590452353602874713527_dp !< base of the natural log, \f$ e \f$
real(dp), parameter :: euler_mascheroni = 0.577216_dp!< Euler Mascheroni constant, typically denoted by \f$ \gamma \f$ ! 0.5772156649015328606065120900824_dp
complex(dp), parameter :: i_ = (0, 1) !< \f$ i = \sqrt{-1} \f$

! Values from National Institute of Standards and Technology (NIST CODATA)
! https://www.nist.gov/pml/fundamental-physical-constants
! Retrieved on April 8th 2021
real(dp), parameter :: hbar_c    = 197.3269804_dp!< \f$\hbar c\f$ in units of MeV fm
real(dp), parameter :: proton_mass  = 938.27208816_dp!< proton mass in units of MeV
real(dp), parameter :: neutron_mass = 939.56542052_dp!< neutron mass in units of MeV
real(dp), parameter :: electron_mass = 0.51099895000_dp!< electron mass in units of MeV
real(dp), parameter :: mu_proton = 2.79284734463_dp!< proton magnetic moment in units of the nuclear magneton \f$ \mu_n \f$
real(dp), parameter :: mu_neutron = -1.91304273_dp!< neutron magnetic moment in units of the nuclear magneton \f$ \mu_n \f$
real(dp), parameter :: alpha = 1/137.035999084_dp!< fine structure constant, dimensionless.

! Values from Particle Data Group (PDG)
! http://pdg.lbl.gov/
! P.A. Zyla et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2020, 083C01 (2020)
real(dp), parameter :: pion_c_mass =  139.57039_dp!< charged pion_mass in units of MeV
real(dp), parameter :: pion_0_mass =  134.9768_dp!< neutral pion_mass in units of MeV
real(dp), parameter :: pion_mass = (2*pion_c_mass + pion_0_mass)/3  !< average pion mass in units of MeV

! Historic, charge independent, recommended value.
real(dp), parameter :: f_pi_n_2 = 0.075_dp !< pion nucleon coupling constant \f$ f^2 \f$. Dimensionless

! ! Original AV18 values
! real(dp), parameter :: hbar_c    = 197.327053_dp!< \f$\hbar c\f$ in units of MeV fm
! real(dp), parameter :: proton_mass  = 938.27231_dp!< proton mass in units of MeV
! real(dp), parameter :: neutron_mass = 939.56563_dp!< neutron mass in units of MeV
! real(dp), parameter :: electron_mass = 0.510999_dp!< electron mass in units of MeV
! real(dp), parameter :: mu_proton = 2.7928474_dp!< proton magnetic moment in units of the nuclear magneton \f$ \mu_n \f$
! real(dp), parameter :: mu_neutron = -1.9130427_dp!< neutron magnetic moment in units of the nuclear magneton \f$ \mu_n \f$
! real(dp), parameter :: alpha = 1/137.035989_dp!< fine structure constant, dimensionless.
! real(dp), parameter :: pion_c_mass =  139.5675_dp!< charged pion_mass in units of MeV
! real(dp), parameter :: pion_0_mass =  134.9739_dp!< neutral pion_mass in units of MeV
! real(dp), parameter :: pion_mass = 138.0363_dp!< average pion mass in units of MeV
! real(dp), parameter :: f_pi_n_2 = 0.075_dp !< pion nucleon coupling constant \f$ f^2 \f$. Dimensionless

end module

max_stars_repo_name: rnavarroperez/nn-scattering-fit

max_issues_repo_name: rnavarroperez/nn-scattering-fit

max_forks_repo_name: rnavarroperez/nn-scattering-fit


text:
C
C	$Id$
C
      SUBROUTINE CSGRCOEF (SIGMA, D,SD)
      DOUBLE PRECISION SIGMA, D, SD
C
C***********************************************************
C
C                                              From SSRFPACK
C                                            Robert J. Renka
C                                  Dept. of Computer Science
C                                       Univ. of North Texas
C                                           renka@cs.unt.edu
C                                                   11/21/96
C
C   This subroutine computes factors involved in the linear
C systems solved by Subroutines CSGRADG and CSSMSGS.
C
C On input:
C
C       SIGMA = Nonnegative tension factor associated with a
C               triangulation arc.
C
C SIGMA is not altered by this routine.
C
C On output:
C
C       D = Diagonal factor.  D = SIG*(SIG*Coshm(SIG) -
C           Sinhm(SIG))/E where E = SIG*Sinh(SIG) - 2*
C           Coshm(SIG).  D > 0, and D = 4 at SIG = 0.
C
C       SD = Off-diagonal factor.  SD = SIG*Sinhm(SIG)/E.
C            SD > 0, and SD = 2 at SIG = 0.
C
C SSRFPACK module required by CSGRCOEF:  CSSNHCSH
C
C Intrinsic function called by CSGRCOEF:  EXP
C
C***********************************************************
C
      DOUBLE PRECISION COSHM, COSHMM, E, EMS, SCM, SIG,
     .                 SINHM, SSINH, SSM
      SIG = SIGMA
      IF (SIG .LT. 1.D-9) THEN
C
C Cubic function:
C
        D = 4.
        SD = 2.
      ELSEIF (SIG .LE. .5) THEN
C
C 0 < SIG .LE. .5.
C
C Use approximations designed to avoid cancellation error
C   in the hyperbolic functions when SIGMA is small.
C
        CALL CSSNHCSH (SIG, SINHM,COSHM,COSHMM)
        E = SIG*SINHM - COSHMM - COSHMM
        D = SIG*(SIG*COSHM-SINHM)/E
        SD = SIG*SINHM/E
      ELSE
C
C SIG > .5.
C
C Scale SINHM, COSHM, and E by 2*EXP(-SIG) in order to
C   avoid overflow when SIGMA is large.
C
        EMS = EXP(-SIG)
        SSINH = 1. - EMS*EMS
        SSM = SSINH - 2.*SIG*EMS
        SCM = (1.-EMS)*(1.-EMS)
        E = SIG*SSINH - SCM - SCM
        D = SIG*(SIG*SCM-SSM)/E
        SD = SIG*SSM/E
      ENDIF
      RETURN
      END

max_stars_repo_name: xylar/cdat

max_issues_repo_name: xylar/cdat

max_forks_repo_name: CDAT/uvcdat


text:
	real x (41), y (41), a (41, 41), b (41, 41)
	real xf1 (82), yf1 (41), rmax (500), rmin (500)
c
	x (1) = -2.
	y (1) = -2.
	do 10 i = 2, 41
		x (i) = x (i-1) + 0.1
		y (i) = x (i)
10              continue
	do 60 i = 1, 41
		do 50 j = 1, 41
			d = x(i) * x(i) + y(j) * y(j) - 4.
			a (i, j) = 0.
			if (d .lt. 0.) a (i, j) = sqrt (-d)
			b (i, j) = - a (i, j)
50                      continue
60              continue
	call setpen (4)
c       call setmtf ('sphere.mtf', 10)
	call page (30., 22., 'SPHERE', 6, 1)
	call setpen (3)
	call region (5., 3., 20., 16., 0, 0, 0)
	call mnmx (-50., -75., 40., x, y, a, 41, 41, 1, 41, 1, 41, s)
	call mnmx (-50., -75., 40., x, y, b, 41, 41, -1, 41, 1, 41, s)
	call prsp (-50., -75., 40., x, y, a, 41, 41, 1, 41, -1, 41,
     *          xf1, yf1, rmax, rmin, 500)
	call prsp (-50., -75., 40., x, y, b, 41, 41, -1, 41, -1, 41,
     *          xf1, yf1, rmax, rmin, 500)
	call endpg ('    ')
	end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


text:
program testburn

  use amrex_constants_module
  use amrex_fort_module, only : rt => amrex_real
  use microphysics_module
  use network
  use eos_module
  use actual_burner_module
  use runtime_init_module, only : runtime_pretty_print

  implicit none

  real(rt) :: dens, temp, dt
  real(rt), dimension(nspec) :: Xin
  type(burn_t) :: state_in, state_out

  integer :: n
  integer :: of

  call microphysics_init()

  dens = 1.e7_rt
  temp = 367821127.2_rt

  Xin(ihe4)  = 0.8_rt
  Xin(ic12)  = 0.1_rt
  Xin(io16)  = 0.1_rt
  Xin(ife56) = 0.0_rt

  dt = 0.001_rt

  open(newunit=of, file="testburn-params.out", status="replace", action="write")
  call runtime_pretty_print(of)
  close(unit=of)

  print *, 'calling the burner...'

  state_in % rho = dens
  state_in % T = temp
  state_in % e = ZERO
  state_in % xn(:) = Xin(:)

  call actual_burner(state_in, state_out, dt, ZERO)

  print *, 'done!'

1000 format(g25.15, 1x, g25.15)

  print *, 'Xin / Xout:'
  do n = 1, nspec
     print 1000, state_in % xn(n), state_out % xn(n)
  enddo

  print *, 'rho_Hnuc: ', dens * (state_out % e - state_in % e) / dt

  print *, "number of RHS evaluations", state_out % n_rhs

  call microphysics_finalize()

end program testburn

max_stars_repo_name: doreenfan/Microphysics

max_issues_repo_name: Youhichka/Microphysics

max_forks_repo_name: Youhichka/Microphysics


text:
        subroutine iddr_rsvd0(m,n,matvect,p1t,p2t,p3t,p4t,
     1                        matvec,p1,p2,p3,p4,krank,u,v,s,ier,
     2                        list,proj,col,work)
c
c       routine iddr_rsvd serves as a memory wrapper
c       for the present routine (please see routine iddr_rsvd
c       for further documentation).
c
        implicit none
        integer m,n,krank,list(n),ier,k
        real*8 p1t,p2t,p3t,p4t,p1,p2,p3,p4,u(m,krank),v(n,krank),
     1         s(krank),proj(krank*(n-krank)),col(m*krank),
     2         work((krank+1)*(m+3*n)+26*krank**2)
        external matvect,matvec
c
c
c       ID a.
c
        call iddr_rid(m,n,matvect,p1t,p2t,p3t,p4t,krank,list,work)
c
c
c       Retrieve proj from work.
c
        do k = 1,krank*(n-krank)
          proj(k) = work(k)
        enddo ! k
c
c
c       Collect together the columns of a indexed by list into col.
c
        call idd_getcols(m,n,matvec,p1,p2,p3,p4,krank,list,col,work)
c
c
c       Convert the ID to an SVD.
c
        call idd_id2svd(m,krank,col,n,list,proj,u,v,s,ier,work)
c
c
        return
        end


text:
module magma_bmm_driver_module
    use, intrinsic :: iso_fortran_env, only : &
            real32, &
            real64, &
            int64

    use, intrinsic :: iso_c_binding, only : &
            c_ptr, &
            c_int, &
            c_float, &
            c_double

    use :: util_api, only : string

    use :: data_api, only : stream

    use :: tensor_api, only : &
            tensor_c_pointer_converter, &
            tensor_builder, &
            tensor, &
            scalar, &
            vector

    use :: pointer_data_bmm_driver_module, only : pointer_data_bmm_driver

    use :: magma_module, only : &
            magma, &
            magma_sgemm_vbatched, &
            magma_dgemm_vbatched, &
            magma_cgemm_vbatched, &
            magma_zgemm_vbatched

    implicit none
    private

    public :: magma_bmm_driver

    type, extends(pointer_data_bmm_driver) :: magma_bmm_driver
        type(magma) :: magma
    contains
        procedure :: bmm_single_real64 => bmm_single_real64
        procedure :: bmm_single_real32 => bmm_single_real32
        procedure :: bmm_single_complex64 => bmm_single_complex64
        procedure :: bmm_single_complex128 => bmm_single_complex128
        procedure :: get_padding_for_sizes => get_padding_for_sizes
        procedure :: initialize => initialize
        procedure :: synchronize => synchronize
        procedure :: cleanup => cleanup
    end type magma_bmm_driver

    interface magma_bmm_driver
        module procedure constructor_empty
        module procedure constructor
    end interface magma_bmm_driver

contains
    function constructor_empty() result(this)
        type(magma_bmm_driver) :: this

        call this%initialize()
    end function constructor_empty

    function constructor(converter, builder) result(this)
        type(tensor_c_pointer_converter), intent(in) :: converter
        type(tensor_builder) :: builder
        type(magma_bmm_driver) :: this

        this = magma_bmm_driver()
        this%converter = converter
        this%builder = builder
    end function constructor

    subroutine initialize(this)
        class(magma_bmm_driver), intent(inout) :: this

        this%magma = magma()
    end subroutine initialize

    subroutine bmm_single_real32(this, batch_count, dst, left, right, m, n, k, &
                transa, transb, alpha, beta, astream)
        class(magma_bmm_driver), intent(in) :: this
        integer, intent(in) :: batch_count
        type(c_ptr), intent(inout) :: dst
        type(c_ptr), intent(in) :: left, right
        type(c_ptr), intent(in) :: m, n, k
        logical, intent(in) :: transa, transb
        real(c_float), intent(in) :: alpha, beta
        type(stream), intent(in), optional :: astream

        integer(c_int) :: my_transa, my_transb
        type(c_ptr) :: stream_queue

        my_transa = 111; my_transb = 111
        if ( transa) my_transa = 112
        if ( transb ) my_transb = 112

        call this%magma%populate_a_queue(stream_queue, astream)

        call magma_sgemm_vbatched( my_transa, my_transb, m, n, k, alpha, left, &
                m, right, k, beta, dst, m, batch_count, stream_queue)

        call this%magma%destroy_queue_if_necessary(stream_queue, astream)
    end subroutine bmm_single_real32

    subroutine bmm_single_real64(this, batch_count, dst, left, right, m, n, k, &
                transa, transb, alpha, beta, astream)
        class(magma_bmm_driver), intent(in) :: this
        integer, intent(in) :: batch_count
        type(c_ptr), intent(inout) :: dst
        type(c_ptr), intent(in) :: left, right
        type(c_ptr), intent(in) :: m, n, k
        logical, intent(in) :: transa, transb
        real(c_double), intent(in) :: alpha, beta
        type(stream), intent(in), optional :: astream

        integer(c_int) :: my_transa, my_transb
        type(c_ptr) :: stream_queue

        my_transa = 111; my_transb = 111
        if ( transa) my_transa = 112
        if ( transb ) my_transb = 112

        call this%magma%populate_a_queue(stream_queue, astream)

        call magma_dgemm_vbatched( my_transa, my_transb, m, n, k, alpha, left, &
                m, right, k, beta, dst, m, batch_count, stream_queue)

        call this%magma%destroy_queue_if_necessary(stream_queue, astream)
    end subroutine bmm_single_real64

    subroutine bmm_single_complex64(this, batch_count, dst, left, right, m, n, k, &
                transa, transb, alpha, beta, astream)
        class(magma_bmm_driver), intent(in) :: this
        integer, intent(in) :: batch_count
        type(c_ptr), intent(inout) :: dst
        type(c_ptr), intent(in) :: left, right
        type(c_ptr), intent(in) :: m, n, k
        logical, intent(in) :: transa, transb
        complex(c_float), intent(in) :: alpha, beta
        type(stream), intent(in), optional :: astream

        integer(c_int) :: my_transa, my_transb
        type(c_ptr) :: stream_queue

        my_transa = 111; my_transb = 111
        if ( transa) my_transa = 112
        if ( transb ) my_transb = 112

        call this%magma%populate_a_queue(stream_queue, astream)

        call magma_cgemm_vbatched( my_transa, my_transb, m, n, k, alpha, left, &
                m, right, k, beta, dst, m, batch_count, stream_queue)

        call this%magma%destroy_queue_if_necessary(stream_queue, astream)
    end subroutine bmm_single_complex64

    subroutine bmm_single_complex128(this, batch_count, dst, left, right, m, n, k, &
                transa, transb, alpha, beta, astream)
        class(magma_bmm_driver), intent(in) :: this
        integer, intent(in) :: batch_count
        type(c_ptr), intent(inout) :: dst
        type(c_ptr), intent(in) :: left, right
        type(c_ptr), intent(in) :: m, n, k
        logical, intent(in) :: transa, transb
        complex(c_double), intent(in) :: alpha, beta
        type(stream), intent(in), optional :: astream

        integer(c_int) :: my_transa, my_transb
        type(c_ptr) :: stream_queue

        my_transa = 111; my_transb = 111
        if ( transa) my_transa = 112
        if ( transb ) my_transb = 112

        call this%magma%populate_a_queue(stream_queue, astream)

        call magma_zgemm_vbatched( my_transa, my_transb, m, n, k, alpha, left, &
                m, right, k, beta, dst, m, batch_count, stream_queue)

        call this%magma%destroy_queue_if_necessary(stream_queue, astream)
    end subroutine bmm_single_complex128

    integer function get_padding_for_sizes(this)
        class(magma_bmm_driver), intent(in) :: this

        get_padding_for_sizes = 1
    end function get_padding_for_sizes

    subroutine synchronize(this, astream)
        class(magma_bmm_driver), intent(in) :: this
        type(stream), intent(in), optional :: astream

        integer(c_int) :: error

        error = this%magma%sync(astream)
    end subroutine synchronize

    subroutine cleanup(this)
        class(magma_bmm_driver), intent(inout) :: this

        call this%magma%cleanup()
        call this%converter%cleanup()
        call this%builder%cleanup()
    end subroutine cleanup
end module magma_bmm_driver_module

max_stars_repo_name: cheshyre/ntcl-algorithms

max_issues_repo_name: cheshyre/ntcl-algorithms

max_forks_repo_name: cheshyre/ntcl-algorithms


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!> @brief Utility module for defining vector types.
!>
!> This module exposes all of the modules which provide VectorType
!> functionality, as enabled at configure time. This includes the Native vector
!> types, and optionally PETSc and/or Trilinos vector types.
!> It also provides abstract factory routines to aid in constructing VectorType
!> objects, as well as interfaces to the BLAS algorithms
!>
!> The objects are initialized with a parameter list. For valid reference lists
!> see @ref VectorTypes::VectorTypes_Declare_ValidParams
!> "VectorTypes_Declare_ValidParams".
!>
!> @par EXAMPLES
!> @code
!> PROGRAM ExampleVector
!>   TYPE(RealVectorType) :: vector
!>   CLASS(VectorType),POINTER :: vec_p
!>   TYPE(ParamType) :: params
!>
!>   CALL params%add("VectorType->n",36_SIK)
!>
!>   CALL vector%init(params)
!>   CALL params%clear()
!>   CALL vector%set(1,10._SRK)
!>   value=vector%get(1)
!>   CALL vector%clear()
!>
!>   ! Create a Trilinos vector using the abstract factory
!>   CALL params%add("VectorType->n",36_SIK)
!>   CALL params%add("VectorType->nlocal",6_SIK)
!>   CALL params%add("VectorType->engine",VM_TRILINOS)
!>   CALL params%add("VectorType->MPI_Comm_ID",MPI_COMM_WORLD)
!>
!>   CALL VectorFactory(vec_p,params)
!>
!>   ! Clean up
!>   CALL vec_p%clear()
!>   DEALLOCATE(vec_p)
!>
!> ENDPROGRAM ExampleVector
!> @endcode
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
MODULE VectorTypes
USE IntrType
USE Allocs
USE ParameterLists
USE VectorTypes_Base
USE VectorTypes_Native
USE VectorTypes_PETSc
USE VectorTypes_Trilinos
USE trilinos_interfaces
USE BLAS1,           ONLY: BLAS1_asum  => BLAS_asum,  &
                           BLAS1_axpy  => BLAS_axpy,  &
                           BLAS1_copy  => BLAS_copy,  &
                           BLAS1_dot   => BLAS_dot,   &
                           BLAS1_iamax => BLAS_iamax, &
                           BLAS1_iamin => BLAS_iamin, &
                           BLAS1_nrm2  => BLAS_nrm2,  &
                           BLAS1_scal  => BLAS_scal,  &
                           BLAS1_swap  => BLAS_swap

#ifdef FUTILITY_HAVE_PETSC
#include <petscversion.h>
#if (((PETSC_VERSION_MAJOR>=3) && (PETSC_VERSION_MINOR>6)) || (PETSC_VERSION_MAJOR>=4))
USE PETSCVEC
#endif
#endif

IMPLICIT NONE

#ifdef FUTILITY_HAVE_PETSC
#include <petscversion.h>
#if (((PETSC_VERSION_MAJOR>=3) && (PETSC_VERSION_MINOR>=6)) || (PETSC_VERSION_MAJOR>=4))
#include <petsc/finclude/petsc.h>
#else
#include <finclude/petsc.h>
#endif
#undef IS
#else
#ifdef HAVE_MPI
#include <mpif.h>
#endif
#endif

PRIVATE
!
! List of public members
PUBLIC :: eVectorType
PUBLIC :: VectorFactory
PUBLIC :: VectorResemble
PUBLIC :: VectorResembleAlloc
PUBLIC :: VectorType
PUBLIC :: NativeVectorType
PUBLIC :: DistributedVectorType
PUBLIC :: RealVectorType
#ifdef FUTILITY_HAVE_PETSC
PUBLIC :: PETScVectorType
#endif
#ifdef FUTILITY_HAVE_Trilinos
PUBLIC :: TrilinosVectorType
#endif
#ifdef HAVE_MPI
PUBLIC :: NativeDistributedVectorType
#endif
! Matrix structure enumerations
PUBLIC :: REAL_NATIVE,DISTRIBUTED_NATIVE
!> Enumerated matrix-vector engines
INTEGER(SIK),PARAMETER,PUBLIC :: VM_PETSC=0,VM_TRILINOS=1,VM_NATIVE=2
PUBLIC :: VectorType_Declare_ValidParams
PUBLIC :: VectorType_Clear_ValidParams
PUBLIC :: BLAS_asum
PUBLIC :: BLAS_axpy
PUBLIC :: BLAS_copy
PUBLIC :: BLAS_dot
PUBLIC :: BLAS_iamax
PUBLIC :: BLAS_iamin
PUBLIC :: BLAS_nrm2
PUBLIC :: BLAS_scal
PUBLIC :: BLAS_swap

!> @brief Allows for copying of an input allocatable VectorType structure/size
INTERFACE VectorResembleAlloc
  !> @copybrief VectorTypes::VectorResembleAllocArray
  !> @copydetails VectorTypes::VectorResembleAllocArray
  MODULE PROCEDURE VectorResembleAllocArray
  !> @copybrief VectorTypes::VectorResembleAllocScal
  !> @copydetails VectorTypes::VectorResembleAllocScal
  MODULE PROCEDURE VectorResembleAllocScal
ENDINTERFACE VectorResembleAlloc

!> @brief Adds to the @ref BLAS1::BLAS_asum "BLAS_asum" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_asum
  !> @copybrief VectorTypes::asum_VectorType
  !> @copydetails VectorTypes::asum_VectorType
  MODULE PROCEDURE asum_VectorType
ENDINTERFACE BLAS_asum

!> @brief Adds to the @ref BLAS1::BLAS_axpy "BLAS_axpy" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_axpy
  !> @copybrief VectorTypes::axpy_scalar_VectorType
  !> @copydetails VectorTypes::axpy_scalar_VectorType
  MODULE PROCEDURE axpy_scalar_VectorType
  !> @copybrief VectorTypes::axpy_vector_VectorType
  !> @copydetails VectorTypes::axpy_vector_VectorType
  MODULE PROCEDURE axpy_vector_VectorType
ENDINTERFACE BLAS_axpy

!> @brief Adds to the @ref BLAS1::BLAS_copy "BLAS_copy" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_copy
  !> @copybrief VectorTypes::copy_VectorType
  !> @copydetails VectorTypes::copy_VectorType
  MODULE PROCEDURE copy_VectorType
ENDINTERFACE BLAS_copy

!> @brief Adds to the @ref BLAS1::BLAS_dot "BLAS_dot" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_dot
  !> @copybrief VectorTypes::dot_VectorType
  !> @copydetails VectorTypes::dot_VectorType
  MODULE PROCEDURE dot_VectorType
ENDINTERFACE BLAS_dot

!> @brief Adds to the @ref BLAS1::BLAS_iamax "BLAS_iamax" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_iamax
  !> @copybrief VectorTypes::iamax_VectorType
  !> @copydetails VectorTypes::iamax_VectorType
  MODULE PROCEDURE iamax_VectorType
ENDINTERFACE BLAS_iamax

!> @brief Adds to the @ref BLAS1::BLAS_iamin "BLAS_iamin" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_iamin
  !> @copybrief VectorTypes::iamin_VectorType
  !> @copydetails VectorTypes::iamin_VectorType
  MODULE PROCEDURE iamin_VectorType
ENDINTERFACE BLAS_iamin

!> @brief Adds to the @ref BLAS1::BLAS_nrm2 "BLAS_nrm2" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_nrm2
  !> @copybrief VectorTypes::nrm2_VectorType
  !> @copydetails VectorTypes::nrm2_VectorType
  MODULE PROCEDURE nrm2_VectorType
ENDINTERFACE BLAS_nrm2

!> @brief Adds to the @ref BLAS1::BLAS_scal "BLAS_scal" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_scal
  !> @copybrief VectorTypes::scal_scalar_VectorType
  !> @copydetails VectorTypes::scal_scalar_VectorType
  MODULE PROCEDURE scal_scalar_VectorType
  !> @copybrief VectorTypes::scal_vector_VectorType
  !> @copydetails VectorTypes::scal_vector_VectorType
  MODULE PROCEDURE scal_vector_VectorType
ENDINTERFACE BLAS_scal

!> @brief Adds to the @ref BLAS1::BLAS_swap "BLAS_swap" interface so that
!> the vector types defined in this module are also supported.
INTERFACE BLAS_swap
  !> @copybrief VectorTypes::swap_VectorType
  !> @copydetails VectorTypes::swap_VectorType
  MODULE PROCEDURE swap_VectorType
ENDINTERFACE BLAS_swap

#ifdef FUTILITY_HAVE_PETSC
!> Scratch variable for petsc error code.
!> It is an integer type.
PetscErrorCode  :: iperr
#endif

!> Name of module
CHARACTER(LEN=*),PARAMETER :: modName='VECTORTYPES'

!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
!> @brief Abstract factory for all enabled VectorTypes
!>
!> @param vector the vector pointer to allocate and initialize. Should be NULL
!> @param params the parameter list to use in determining the type of and
!> initializing the vector.
!>
!> This is an abstract factory routine for the base VectorType. It will use the
!> 'engine' parameter in the passed parameter list to determine which type of
!> Vector to allocate the pointer to. The parameter list is then used to
!> initialize the vector.
SUBROUTINE VectorFactory(vector, params)
  CHARACTER(LEN=*),PARAMETER :: myName="VectorFactory"
  CLASS(VectorType),POINTER,INTENT(INOUT) :: vector
  CLASS(ParamType),INTENT(IN) :: params
  !
  INTEGER(SIK) :: engine,vecType

  IF(ASSOCIATED(vector)) THEN
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Vector pointer is already allocated")
    RETURN
  ENDIF

  ! Default to the native engine
  engine=VM_NATIVE
  ! Default to RealVectorType
  vecType=REAL_NATIVE

  IF(params%has("VectorType->engine")) THEN
    CALL params%get("VectorType->engine",engine)
  ENDIF

  IF(params%has("VectorType->vecType")) THEN
    CALL params%get("VectorType->vecType",vecType)
  ENDIF

  SELECTCASE(engine)
  CASE(VM_NATIVE)
    SELECTCASE(vecType)
    CASE(REAL_NATIVE)
      ALLOCATE(RealVectorType :: vector)
#ifdef HAVE_MPI
    CASE(DISTRIBUTED_NATIVE)
      ALLOCATE(NativeDistributedVectorType :: vector)
#endif
    CASE DEFAULT
      CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
          "Unrecognized vector type requested")
    ENDSELECT
  CASE(VM_PETSC)
#ifdef FUTILITY_HAVE_PETSC
    ALLOCATE(PETScVectorType :: vector)
#endif
  CASE(VM_TRILINOS)
#ifdef FUTILITY_HAVE_Trilinos
    ALLOCATE(TrilinosVectorType :: vector)
#endif
  CASE DEFAULT
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Unrecognized vector engine requested")
  ENDSELECT

  CALL vector%init(params)

ENDSUBROUTINE VectorFactory
!
!-------------------------------------------------------------------------------
!> @brief Create a new vector of compatible size and type to the input vector pointer and return a pointer
!>
!> @param dest Pointer to the destination VectorType to allocate and construct.
!> This should be NULL
!> @param source Pointer to source vector type to use in determining the type and parameters
!> of the dest vector
!> @param params the parameters to use in overriding settings from the source
!> vector
!>
!> For now, the source vector shall be initialized, though in the future it
!> might be nice to support uninitialized vectors. Unlike the corresponding routine for MatrixTypes,
!> this routine will attempt to adopt parameters from the source vector to initialize the dest vector.
!> This is only done for required parameters that are not provided in the passed
!> parameter list. Providing the parameters on the parameter list will override
!> the corresponding parameters from the source matrix. This behavior will be
!> consistent for all VectorResemble routines
SUBROUTINE VectorResemble(dest, source, params)
  CHARACTER(LEN=*),PARAMETER :: myName="VectorResemble"
  CLASS(VectorType),POINTER,INTENT(INOUT) :: dest
  CLASS(VectorType),POINTER,INTENT(IN) :: source
  CLASS(ParamType),INTENT(INOUT) :: params

  IF(.NOT. ASSOCIATED(source)) THEN
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Source vector is not associated")
    RETURN
  ENDIF

  IF(.NOT. source%isInit) THEN
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Source vector is not initialized")
  ENDIF

  IF(ASSOCIATED(dest)) THEN
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Destination vector is already associated")
    RETURN
  ENDIF

  IF(.NOT. params%has("VectorType->n")) THEN
    CALL params%add("VectorType->n",source%n)
  ENDIF
  SELECT TYPE(source)
  CLASS IS(DistributedVectorType)
    IF(.NOT. params%has("VectorType->nlocal")) THEN
      CALL params%add("VectorType->nlocal",source%nlocal)
    ENDIF
    IF(.NOT. params%has("VectorType->MPI_Comm_Id")) THEN
      CALL params%add("VectorType->MPI_Comm_Id",source%comm)
    ENDIF
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    IF(.NOT. params%has("VectorType->chunkSize")) THEN
      CALL params%add("VectorType->chunkSize",source%chunkSize)
    ENDIF
    IF(.NOT. params%has("VectorType->MPI_Comm_Id")) THEN
      CALL params%add("VectorType->MPI_Comm_Id",source%comm)
    ENDIF
#endif
  ENDSELECT

  SELECT TYPE(source)
  TYPE IS(RealVectorType)
    ALLOCATE(RealVectorType :: dest)
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    ALLOCATE(NativeDistributedVectorType :: dest)
#endif
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    ALLOCATE(PETScVectorType :: dest)
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosVectorType)
    ALLOCATE(TrilinosVectorType :: dest)
#endif
  ENDSELECT
  CALL dest%init(params)

  CALL params%clear()
ENDSUBROUTINE VectorResemble
!
!-------------------------------------------------------------------------------
!> @brief Create a new vector of compatible size and type to the input vector
!>
!> @param dest the destination VectorType array to allocate and initialize.
!> @param source the vector type to use in determining the type and parameters
!> of the dest vector
!> @param inparams the parameters to use in overriding settings from the source
!> vector
!>
SUBROUTINE VectorResembleAllocScal(dest, source, inparams)
  CLASS(VectorType),ALLOCATABLE,INTENT(INOUT) :: dest
  CLASS(VectorType),INTENT(IN) :: source
  CLASS(ParamType),INTENT(INOUT),OPTIONAL :: inparams

  CHARACTER(LEN=*),PARAMETER :: myName="VectorResembleAllocScal"
  TYPE(ParamType) :: params

  IF(.NOT. source%isInit) THEN
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Source vector is not initialized")
  ENDIF

  IF(ALLOCATED(dest)) THEN
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Destination vector is already allocated")
    RETURN
  ENDIF

  IF(PRESENT(inparams)) params=inparams
  IF(.NOT. params%has("VectorType->n")) THEN
    CALL params%add("VectorType->n",source%n)
  ENDIF
  SELECT TYPE(source)
  CLASS IS(DistributedVectorType)
    IF(.NOT. params%has("VectorType->nlocal")) THEN
      CALL params%add("VectorType->nlocal",source%nlocal)
    ENDIF
    IF(.NOT. params%has("VectorType->MPI_Comm_Id")) THEN
      CALL params%add("VectorType->MPI_Comm_Id",source%comm)
    ENDIF
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    IF(.NOT. params%has("VectorType->nlocal")) THEN
      CALL params%add("VectorType->nlocal",SIZE(source%b))
    ENDIF
    IF(.NOT. params%has("VectorType->chunkSize")) THEN
      CALL params%add("VectorType->chunkSize",source%chunkSize)
    ENDIF
    IF(.NOT. params%has("VectorType->MPI_Comm_Id")) THEN
      CALL params%add("VectorType->MPI_Comm_Id",source%comm)
    ENDIF
#endif
  ENDSELECT

  SELECT TYPE(source)
  TYPE IS(RealVectorType)
    ALLOCATE(RealVectorType :: dest)
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    ALLOCATE(NativeDistributedVectorType :: dest)
#endif
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    ALLOCATE(PETScVectorType :: dest)
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosVectorType)
    ALLOCATE(TrilinosVectorType :: dest)
#endif
  ENDSELECT
  CALL dest%init(params)

  IF(PRESENT(inparams)) CALL inparams%clear()
  CALL params%clear()

ENDSUBROUTINE VectorResembleAllocScal
!
!-------------------------------------------------------------------------------
!> @brief Create a new vector of compatible size and type to the input vector
!>
!> @param dest the destination VectorType array to allocate and initialize.
!> @param source the vector type to use in determining the type and parameters
!> of the dest vector
!> @param inparams the parameters to use in overriding settings from the source
!> vector
!> @param nvec length to allocate array of vectors to
!>
SUBROUTINE VectorResembleAllocArray(dest, source, nvec, inparams)
  CLASS(VectorType),ALLOCATABLE,INTENT(INOUT) :: dest(:)
  CLASS(VectorType),INTENT(IN) :: source
  INTEGER(SIK),INTENT(IN) :: nvec
  CLASS(ParamType),INTENT(INOUT),OPTIONAL :: inparams

  CHARACTER(LEN=*),PARAMETER :: myName="VectorResembleAllocArray"
  INTEGER(SIK) :: i
  TYPE(ParamType) :: params

  IF(.NOT. source%isInit) THEN
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Source vector is not initialized")
  ENDIF

  IF(ALLOCATED(dest)) THEN
    CALL eVectorType%raiseError(modName//"::"//myName//" - "// &
        "Destination vector is already allocated")
    RETURN
  ENDIF

  IF(PRESENT(inparams)) params=inparams
  IF(.NOT. params%has("VectorType->n")) THEN
    CALL params%add("VectorType->n",source%n)
  ENDIF
  SELECT TYPE(source)
  CLASS IS(DistributedVectorType)
    IF(.NOT. params%has("VectorType->nlocal")) THEN
      CALL params%add("VectorType->nlocal",source%nlocal)
    ENDIF
    IF(.NOT. params%has("VectorType->MPI_Comm_Id")) THEN
      CALL params%add("VectorType->MPI_Comm_Id",source%comm)
    ENDIF
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    IF(.NOT. params%has("VectorType->nlocal")) THEN
      CALL params%add("VectorType->nlocal",SIZE(source%b))
    ENDIF
    IF(.NOT. params%has("VectorType->chunkSize")) THEN
      CALL params%add("VectorType->chunkSize",source%chunkSize)
    ENDIF
    IF(.NOT. params%has("VectorType->MPI_Comm_Id")) THEN
      CALL params%add("VectorType->MPI_Comm_Id",source%comm)
    ENDIF
#endif
  ENDSELECT

  SELECT TYPE(source)
  TYPE IS(RealVectorType)
    ALLOCATE(RealVectorType :: dest(nvec))
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    ALLOCATE(NativeDistributedVectorType :: dest(nvec))
#endif
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    ALLOCATE(PETScVectorType :: dest(nvec))
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosVectorType)
    ALLOCATE(TrilinosVectorType :: dest(nvec))
#endif
  ENDSELECT
  DO i=1,nvec
    CALL dest(i)%init(params)
  ENDDO

  IF(PRESENT(inparams)) CALL inparams%clear()
  CALL params%clear()

ENDSUBROUTINE VectorResembleAllocArray
!
!-------------------------------------------------------------------------------
!> @brief Function provides an interface to vector absolute value summation
!> of a vector (x).
!> @param thisVector derived vector type
!> @param n the size of the vector @c x
!> @param incx the increment to use when looping over elements in @c x
!> @return r the sum of the absolute values of @c x
!>
FUNCTION asum_VectorType(thisVector,n,incx) RESULT(r)
  CHARACTER(LEN=*),PARAMETER :: myName='asum_VectorType'
  CLASS(VectorType),INTENT(INOUT)     :: thisVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  REAL(SRK) :: r

  SELECT TYPE(thisVector)
  TYPE IS(RealVectorType)
    IF(PRESENT(n) .AND. PRESENT(incx)) THEN
      r=BLAS1_asum(n,thisVector%b,incx)
    ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
      r=BLAS1_asum(n,thisVector%b)
    ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx)) THEN
      r=BLAS1_asum(thisVector%b,incx)
    ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
      r=BLAS1_asum(thisVector%b)
    ENDIF
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble(iperr)
    IF(iperr == 0) CALL VecNorm(thisVector%b,NORM_1,r,iperr)
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosVectorType)
    CALL ForPETRA_VecSUM(thisVector%b,r)
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT
ENDFUNCTION asum_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to compute the result of a vector (y)
!> plus a vector (x) times a scalar (a).
!> @param thisVector derived vector type
!> @param newVector resulting derived vector type
!> @param a the constant to multiply with @c x
!> @param n the size of the vectors @c x and @c y
!> @param incx the increment to use when looping over elements in @c x
!> @param incy the increment to use when looping over elements in @c y
!>
SUBROUTINE axpy_scalar_VectorType(thisVector,newVector,a,n,incx,incy)
  CHARACTER(LEN=*),PARAMETER :: myName='axpy_scalar_VectorType'
  CLASS(VectorType),INTENT(INOUT)  :: thisVector
  CLASS(VectorType),INTENT(INOUT)  :: newVector
  REAL(SRK),INTENT(IN),OPTIONAL :: a
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incy
  REAL(SRK) :: alpha

  alpha=1._SRK
  IF(PRESENT(a)) alpha=a

  SELECT TYPE(thisVector); TYPE IS(RealVectorType)
    SELECT TYPE(newVector); TYPE IS(RealVectorType)
      IF(PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_axpy(n,alpha,thisVector%b,incx,newVector%b,incy)
      ELSEIF(.NOT. PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_axpy(SIZE(thisVector%b),alpha,thisVector%b,incx,newVector%b,incy)
      ELSEIF(.NOT. PRESENT(n) .AND. PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_axpy(alpha,thisVector%b,newVector%b,incx)
      ELSEIF(.NOT. PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_axpy(alpha,thisVector%b,newVector%b,incy)
      ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_axpy(n,alpha,thisVector%b,newVector%b,incy)
      ELSEIF(PRESENT(n) .AND. PRESENT(incx) .AND. .NOT. PRESENT(incy)) THEN
        CALL BLAS1_axpy(n,alpha,thisVector%b,newVector%b,incx)
      ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_axpy(n,alpha,thisVector%b,newVector%b)
      ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_axpy(alpha,thisVector%b,newVector%b)
      ENDIF
    ENDSELECT
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    SELECT TYPE(newVector); TYPE IS(NativeDistributedVectorType)
      IF(PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_axpy(SIZE(thisVector%b),alpha,thisVector%b,incx,newVector%b,incy)
      ELSEIF(PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_axpy(alpha,thisVector%b,newVector%b,incx)
      ELSEIF(.NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_axpy(alpha,thisVector%b,newVector%b,incy)
      ELSEIF(.NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_axpy(alpha,thisVector%b,newVector%b)
      ENDIF
    ENDSELECT
#endif
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    SELECT TYPE(newVector)
    TYPE IS(PETScVectorType)
      IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble(iperr)
      IF(iperr == 0) CALL VecAXPY(newVector%b,alpha,thisVector%b,iperr)
    ENDSELECT
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosVectorType)
    SELECT TYPE(newVector)
    TYPE IS(TrilinosVectorType)
      IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble()
      CALL ForPETRA_VecAXPY(newVector%b,thisVector%b,alpha,1.0_SRK)
    ENDSELECT
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT
ENDSUBROUTINE axpy_scalar_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to compute the result of a vector (y)
!> plus a vector (x) times a vector (aVector).
!> @param thisVector derived vector type
!> @param newVector resulting derived vector type
!> @param a the constant to multiply with @c x
!> @param n the size of the vectors @c x and @c y
!> @param incx the increment to use when looping over elements in @c x
!> @param incy the increment to use when looping over elements in @c y
!>
!> TIBWSFB: This is inefficient. It makes a copy of the vectors, which
!> is a pretty big waste. The one upside to this is the algorithm can be
!> implemented using a type-agnostic approach. But we still select type
!> everything for some reason. Worst of both worlds
SUBROUTINE axpy_vector_VectorType(thisVector,newVector,aVector,n,incx,incy)
  CHARACTER(LEN=*),PARAMETER :: myName='axpy_vector_VectorType'
  CLASS(VectorType),INTENT(INOUT)  :: thisVector
  CLASS(VectorType),INTENT(INOUT)  :: newVector
  CLASS(VectorType),INTENT(INOUT)  :: aVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incy
  REAL(SRK),ALLOCATABLE :: tmpthis(:),tmpnew(:),tmpa(:)

  SELECT TYPE(thisVector)
  TYPE IS(RealVectorType)
    SELECT TYPE(newVector)
    TYPE IS(RealVectorType)
      SELECT TYPE(aVector)
      TYPE IS(RealVectorType)
        ALLOCATE(tmpthis(thisVector%n))
        ALLOCATE(tmpnew(newVector%n))
        ALLOCATE(tmpa(aVector%n))
        CALL thisVector%get(tmpthis)
        CALL newVector%get(tmpnew)
        CALL aVector%get(tmpa)
      ENDSELECT
    ENDSELECT
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    SELECT TYPE(newVector)
    TYPE IS(PETScVectorType)
      SELECT TYPE(aVector)
      TYPE IS(PETScVectorType)
        ALLOCATE(tmpthis(thisVector%n))
        ALLOCATE(tmpnew(newVector%n))
        ALLOCATE(tmpa(aVector%n))
        CALL thisVector%get(tmpthis)
        CALL newVector%get(tmpnew)
        CALL aVector%get(tmpa)
      ENDSELECT
    ENDSELECT
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT

  IF(PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
    CALL BLAS1_axpy(n,tmpa,tmpthis,incx,tmpnew,incy)
  ELSEIF(PRESENT(n) .AND. PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
    CALL BLAS1_axpy(n,tmpa,tmpthis,tmpnew,incx)
  ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
    CALL BLAS1_axpy(n,tmpa,tmpthis,tmpnew,incy)
  ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
    CALL BLAS1_axpy(n,tmpa,tmpthis,tmpnew)
  ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
    CALL BLAS1_axpy(tmpa,tmpthis,tmpnew)
  ENDIF

  CALL newVector%set(tmpnew)

  DEALLOCATE(tmpthis)
  DEALLOCATE(tmpnew)
  DEALLOCATE(tmpa)
ENDSUBROUTINE axpy_vector_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to copy a vector (x) to another
!> vector (y).
!> @param x derived vector type
!> @param y resulting derived vector type
!> @param n the size of the vectors @c x and @c y
!> @param incx the increment to use when looping over elements in @c x
!> @param incy the increment to use when looping over elements in @c y
!>
SUBROUTINE copy_VectorType(thisVector,newVector,n,incx,incy)
  CHARACTER(LEN=*),PARAMETER :: myName='copy_VectorType'
  CLASS(VectorType),INTENT(INOUT)  :: thisVector
  CLASS(VectorType),INTENT(INOUT)  :: newVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incy

  SELECT TYPE(thisVector); TYPE IS(RealVectorType)
    SELECT TYPE(newVector); TYPE IS(RealVectorType)
      IF(PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_copy(n,thisVector%b,incx,newVector%b,incy)
      ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_copy(thisVector%b,incx,newVector%b,incy)
      ELSEIF(PRESENT(n) .AND. PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_copy(n,thisVector%b,newVector%b,incx)
      ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_copy(n,thisVector%b,newVector%b,incy)
      ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_copy(thisVector%b,newVector%b,incx)
      ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_copy(thisVector%b,newVector%b,incy)
      ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_copy(n,thisVector%b,newVector%b)
      ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_copy(thisVector%b,newVector%b)
      ENDIF
    ENDSELECT
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    SELECT TYPE(newVector); TYPE IS(NativeDistributedVectorType)
      IF(PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_copy(thisVector%b,incx,newVector%b,incy)
      ELSEIF(PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_copy(thisVector%b,newVector%b,incx)
      ELSEIF(.NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_copy(thisVector%b,newVector%b,incy)
      ELSEIF(.NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_copy(thisVector%b,newVector%b)
      ENDIF
    ENDSELECT
#endif
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    SELECT TYPE(newVector)
    TYPE IS(PETScVectorType)
      IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble(iperr)
      IF(iperr == 0) CALL VecCopy(thisVector%b,newVector%b,iperr)
    ENDSELECT
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosVectorType)
    SELECT TYPE(newVector)
    TYPE IS(TrilinosVectorType)
      IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble()
      CALL ForPETRA_VecCopy(newVector%b,thisVector%b)
    ENDSELECT
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT
ENDSUBROUTINE copy_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to compute dot product of two
!> vectors (x and y).
!> @param thisVector derived vector type.
!> @param thatVector derived vector type.
!> @param n the size of the vectors @c x and @c y
!> @param incx the increment to use when looping over elements in @c x
!> @param incy the increment to use when looping over elements in @c y
!> @return r the dot product of @c x and @c y
!>
FUNCTION dot_VectorType(thisVector,thatVector,n,incx,incy)  RESULT(r)
  CHARACTER(LEN=*),PARAMETER :: myName='dot_VectorType'
  CLASS(VectorType),INTENT(INOUT)     :: thisVector
  CLASS(VectorType),INTENT(INOUT)     :: thatVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incy
  REAL(SRK) :: r
#ifdef HAVE_MPI
  INTEGER(SIK) :: mpierr
#endif
  SELECT TYPE(thisVector); TYPE IS(RealVectorType)
    SELECT TYPE(thatVector); TYPE IS(RealVectorType)
      IF(PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
        r=BLAS1_dot(n,thisVector%b,incx,thatVector%b,incy)
      ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
        r=BLAS1_dot(thisVector%b,incx,thatVector%b,incy)
      ELSEIF(PRESENT(n) .AND. PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        r=BLAS1_dot(n,thisVector%b,thatVector%b,incx)
      ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        r=BLAS1_dot(n,thisVector%b,thatVector%b,incy)
      ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        r=BLAS1_dot(n,thisVector%b,thatVector%b)
      ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        r=BLAS1_dot(thisVector%b,thatVector%b,incx)
      ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        r=BLAS1_dot(thisVector%b,thatVector%b,incy)
      ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        r=BLAS1_dot(thisVector%b,thatVector%b)
      ENDIF
    ENDSELECT
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    SELECT TYPE(thatVector); TYPE IS(NativeDistributedVectorType)
      IF(PRESENT(incx) .AND. PRESENT(incy)) THEN
        r=BLAS1_dot(thisVector%b,incx,thatVector%b,incy)
      ELSEIF(PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        r=BLAS1_dot(thisVector%b,thatVector%b,incx)
      ELSEIF(.NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        r=BLAS1_dot(thisVector%b,thatVector%b,incy)
      ELSEIF(.NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        r=BLAS1_dot(thisVector%b,thatVector%b)
      ENDIF
      CALL MPI_Allreduce(MPI_IN_PLACE,r,1,MPI_DOUBLE_PRECISION,MPI_SUM,thisVector%comm,mpierr)
    ENDSELECT
#endif
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    SELECT TYPE(thatVector)
    TYPE IS(PETScVectorType)
      IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble(iperr)
      IF(iperr == 0) CALL VecTDot(thisVector%b,thatVector%b,r,iperr)
    ENDSELECT
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
       modName//'::'//myName//' - This interface is not available.')
  ENDSELECT

ENDFUNCTION dot_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to compute the inex of the absolute
!> maximum of a vector (x).
!> @param thisVector derived vector type.
!> @param n the size of the vectors @c x
!> @param incx the increment to use when looping over elements in @c x
!> @return imax index of the absolute max of @c y
!>
FUNCTION iamax_VectorType(thisVector,n,incx)  RESULT(imax)
  CHARACTER(LEN=*),PARAMETER :: myName='iamax_VectorType'
  CLASS(VectorType),INTENT(INOUT) :: thisVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  REAL(SRK),ALLOCATABLE :: tmpthis(:)
  INTEGER(SIK) :: imax

  SELECT TYPE(thisVector)
  TYPE IS(RealVectorType)
    ALLOCATE(tmpthis(thisVector%n))
    CALL thisVector%get(tmpthis)
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    ALLOCATE(tmpthis(thisVector%n))
    CALL thisVector%get(tmpthis)
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT

  IF(PRESENT(n) .AND. PRESENT(incx)) THEN
    imax=BLAS1_iamax(n,tmpthis,incx)
  ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx)) THEN
    imax=BLAS1_iamax(tmpthis,incx)
  ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
    imax=BLAS1_iamax(n,tmpthis)
  ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
    imax=BLAS1_iamax(tmpthis)
  ENDIF
  DEALLOCATE(tmpthis)
ENDFUNCTION iamax_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to compute the index of the absolute
!> minimum of a vector (x).
!> @param thisVector derived vector type.
!> @param n the size of the vectors @c x
!> @param incx the increment to use when looping over elements in @c x
!> @return imin index of the absolute min of @c x
!>

FUNCTION iamin_VectorType(thisVector,n,incx)  RESULT(imin)
  CHARACTER(LEN=*),PARAMETER :: myName='iamin_VectorType'
  CLASS(VectorType),INTENT(INOUT) :: thisVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  REAL(SRK),ALLOCATABLE :: tmpthis(:)
  INTEGER(SIK) :: imin

  SELECT TYPE(thisVector)
  TYPE IS(RealVectorType)
    ALLOCATE(tmpthis(thisVector%n))
    CALL thisVector%get(tmpthis)
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    ALLOCATE(tmpthis(thisVector%n))
    CALL thisVector%get(tmpthis)
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT

  IF(PRESENT(n) .AND. PRESENT(incx)) THEN
    imin=BLAS1_iamin(n,tmpthis,incx)
  ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx)) THEN
    imin=BLAS1_iamin(tmpthis,incx)
  ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
    imin=BLAS1_iamin(n,tmpthis)
  ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
    imin=BLAS1_iamin(tmpthis)
  ENDIF
  DEALLOCATE(tmpthis)
ENDFUNCTION iamin_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to compute the 2-norm of a
!> vector (x).
!> @param thisVector derived vector type
!> @param n the size of the vectors @c x
!> @param incx the increment to use when looping over elements in @c x
!> @return norm2 the 2-norm of @c x
!>
FUNCTION nrm2_VectorType(thisVector,n,incx)  RESULT(norm2)
  CHARACTER(LEN=*),PARAMETER :: myName='nrm2_VectorType'
  CLASS(VectorType),INTENT(INOUT) :: thisVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  REAL(SRK) :: norm2
#ifdef HAVE_MPI
  INTEGER(SIK) :: mpierr
#endif

  SELECT TYPE(thisVector); TYPE IS(RealVectorType)
    IF(PRESENT(n) .AND. PRESENT(incx)) THEN
      norm2=BLAS1_nrm2(n,thisVector%b,incx)
    ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx)) THEN
      norm2=BLAS1_nrm2(thisVector%b,incx)
    ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
      norm2=BLAS1_nrm2(n,thisVector%b)
    ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
      norm2=BLAS1_nrm2(thisVector%b)
    ENDIF
#ifdef HAVE_MPI
  TYPE IS(NativeDistributedVectorType)
    IF(PRESENT(incx)) THEN
      norm2=BLAS1_dot(thisVector%b,thisVector%b)
    ELSE
      norm2=BLAS1_dot(thisVector%b,thisVector%b)
    ENDIF
    CALL MPI_Allreduce(MPI_IN_PLACE,norm2,1,MPI_DOUBLE_PRECISION,MPI_SUM,thisVector%comm,mpierr)
    norm2 = SQRT(norm2)
#endif
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble(iperr)
    IF(iperr == 0) CALL VecNorm(thisVector%b,NORM_2,norm2,iperr)
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosVectorType)
    IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble()
    CALL ForPETRA_VecNorm2(thisVector%b,norm2)
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
       modName//'::'//myName//' - This interface is not available.')
  ENDSELECT

ENDFUNCTION nrm2_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to scale a vector (x) by a
!> scalar (a).
!> @param thisVector derived vector type
!> @param a the constant to multiply with @c x
!> @param n the size of the vectors @c x
!> @param incx the increment to use when looping over elements in @c x
!>
SUBROUTINE scal_scalar_VectorType(thisVector,a,n,incx)
  CHARACTER(LEN=*),PARAMETER :: myName='scal_scalar_VectorType'
  CLASS(VectorType),INTENT(INOUT) :: thisVector
  REAL(SRK),INTENT(IN) :: a
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx

  SELECT TYPE(thisVector)
  TYPE IS(RealVectorType)
    IF(PRESENT(n) .AND. PRESENT(incx)) THEN
      CALL BLAS1_scal(n,a,thisVector%b,incx)
    ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx)) THEN
      CALL BLAS1_scal(a,thisVector%b,incx)
    ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
      CALL BLAS1_scal(n,a,thisVector%b)
    ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
      CALL BLAS1_scal(a,thisVector%b)
    ENDIF
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble(iperr)
    IF(iperr == 0) CALL VecScale(thisVector%b,a,iperr)
#endif
#ifdef FUTILITY_HAVE_Trilinos
  TYPE IS(TrilinosVectorType)
    CALL ForPETRA_VecScale(thisVector%b,a)
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT
ENDSUBROUTINE scal_scalar_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine provides an interface to scale a vector (x) by another
!> vector (aVector).
!> @param thisVector derived vector type
!> @param a the constant to multiply with @c x
!> @param n the size of the vectors @c x
!> @param incx the increment to use when looping over elements in @c x
!>
SUBROUTINE scal_vector_VectorType(thisVector,aVector,n,incx)
  CHARACTER(LEN=*),PARAMETER :: myName='scal_vector_VectorType'
  CLASS(VectorType),INTENT(INOUT) :: thisVector
  CLASS(VectorType),INTENT(INOUT) :: aVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  REAL(SRK),ALLOCATABLE :: tmpthis(:),tmpa(:)

  SELECT TYPE(thisVector)
  TYPE IS(RealVectorType)
    SELECT TYPE(aVector)
    TYPE IS(RealVectorType)
      ALLOCATE(tmpthis(thisVector%n))
      ALLOCATE(tmpa(aVector%n))
      CALL thisVector%get(tmpthis)
      CALL aVector%get(tmpa)
    ENDSELECT
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    SELECT TYPE(aVector)
    TYPE IS(PETScVectorType)
      ALLOCATE(tmpthis(thisVector%n))
      ALLOCATE(tmpa(aVector%n))
      CALL thisVector%get(tmpthis)
      CALL aVector%get(tmpa)
    ENDSELECT
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT

  IF(PRESENT(n) .AND. PRESENT(incx)) THEN
    CALL BLAS1_scal(n,tmpa,tmpthis,incx)
  ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx)) THEN
    CALL BLAS1_scal(tmpa,tmpthis,incx)
  ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
    CALL BLAS1_scal(n,tmpa,tmpthis)
  ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx)) THEN
    CALL BLAS1_scal(tmpa,tmpthis)
  ENDIF

  CALL thisVector%set(tmpthis)

  DEALLOCATE(tmpthis)
  DEALLOCATE(tmpa)
ENDSUBROUTINE scal_vector_VectorType
!
!-------------------------------------------------------------------------------
!> @brief Subroutine swaps a vector @c x with a vector @c y
!> @param thisVector derived vector type
!> @param thatVector derived vector type
!> @param n the size of the vectors @c x and @c y
!> @param incx the increment to use when looping over elements in @c x
!> @param incy the increment to use when looping over elements in @c y
!>
SUBROUTINE swap_VectorType(thisVector,thatVector,n,incx,incy)
  CHARACTER(LEN=*),PARAMETER :: myName='swap_VectorType'
  CLASS(VectorType),INTENT(INOUT) :: thisVector
  CLASS(VectorType),INTENT(INOUT) :: thatVector
  INTEGER(SIK),INTENT(IN),OPTIONAL :: n
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incx
  INTEGER(SIK),INTENT(IN),OPTIONAL :: incy

  SELECT TYPE(thisVector)
  TYPE IS(RealVectorType)
    SELECT TYPE(thatVector)
    TYPE IS(RealVectorType)
      IF(PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_swap(n,thisVector%b,incx,thatVector%b,incy)
      ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_swap(thisVector%b,incx,thatVector%b,incy)
      ELSEIF(PRESENT(n) .AND. PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_swap(n,thisVector%b,thatVector%b,incx)
      ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_swap(n,thisVector%b,thatVector%b,incy)
      ELSEIF(.NOT.PRESENT(n) .AND. PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_swap(thisVector%b,thatVector%b,incx)
      ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. PRESENT(incy)) THEN
        CALL BLAS1_swap(thisVector%b,thatVector%b,incy)
      ELSEIF(PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_swap(n,thisVector%b,thatVector%b)
      ELSEIF(.NOT.PRESENT(n) .AND. .NOT.PRESENT(incx) .AND. .NOT.PRESENT(incy)) THEN
        CALL BLAS1_swap(thisVector%b,thatVector%b)
      ENDIF
    ENDSELECT
#ifdef FUTILITY_HAVE_PETSC
  TYPE IS(PETScVectorType)
    SELECT TYPE(thatVector)
    TYPE IS(PETScVectorType)
      IF(.NOT.thisVector%isAssembled) CALL thisVector%assemble(iperr)
      IF(.NOT.thatVector%isAssembled) CALL thatVector%assemble(iperr)
      IF(iperr == 0) CALL VecSwap(thisVector%b,thatVector%b,iperr)
    ENDSELECT
#endif
  CLASS DEFAULT
    CALL eVectorType%raiseFatalError('Incorrect call to '// &
        modName//'::'//myName//' - This interface is not available.')
  ENDSELECT
ENDSUBROUTINE swap_VectorType
!
ENDMODULE VectorTypes

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: picmc/Futility


text: #redirect Highway 50