SUBROUTINE read_inflow_dgfem(inflow_file)

!-------------------------------------------------------------
!
! Reads fixed inflow BC
!
!-------------------------------------------------------------
USE invar
use precision_module, only: dp
IMPLICIT NONE

CHARACTER(30), INTENT(IN) :: inflow_file
INTEGER :: octant,n
INTEGER :: ix,iy,iz,jx,jy,jz,dir,l
INTEGER :: sgn_mu,sgn_eta,sgn_xi
real(kind=dp) :: sgx,sgy,sgz

ALLOCATE(  tfrbc(nz,nx,0:lambda,0:lambda,apo,4),&
           tbabc(nz,nx,0:lambda,0:lambda,apo,4),&
           tlebc(ny,nz,0:lambda,0:lambda,apo,4),&
           tribc(ny,nz,0:lambda,0:lambda,apo,4),&
           tbobc(nx,ny,0:lambda,0:lambda,apo,4),&
           ttobc(nx,ny,0:lambda,0:lambda,apo,4) )

IF (solvertype == "LD" .or. solvertype == "DENSE") THEN
  ALLOCATE( frbc(dofpc,nz,nx,apo,4,1) )
  ALLOCATE( babc(dofpc,nz,nx,apo,4,1) )
  ALLOCATE( lebc(dofpc,ny,nz,apo,4,1) )
  ALLOCATE( ribc(dofpc,ny,nz,apo,4,1) )
  ALLOCATE( bobc(dofpc,nx,ny,apo,4,1) )
  ALLOCATE( tobc(dofpc,nx,ny,apo,4,1) )
ELSE IF (solvertype == "LAGRANGE") THEN
  ALLOCATE( frbc(ordcb,nz,nx,apo,4,1) )
  ALLOCATE( babc(ordcb,nz,nx,apo,4,1) )
  ALLOCATE( lebc(ordcb,ny,nz,apo,4,1) )
  ALLOCATE( ribc(ordcb,ny,nz,apo,4,1) )
  ALLOCATE( bobc(ordcb,nx,ny,apo,4,1) )
  ALLOCATE( tobc(ordcb,nx,ny,apo,4,1) )
END IF

OPEN(UNIT=12, FILE=inflow_file,STATUS = "OLD", ACTION = "READ",FORM='UNFORMATTED')
DO octant=1,8
   DO n=1,apo
      READ(12) sgn_mu,sgn_eta,sgn_xi
      ! >> front/back
      IF(sgn_xi>0 .and. sgn_mu>0) THEN
         dir=1
      ELSE IF (sgn_xi>0 .and. sgn_mu<0) THEN
         dir=2
      ELSE IF (sgn_xi<0 .and. sgn_mu>0) THEN
         dir=3
      ELSE
         dir=4
      END IF
      !
      IF (sgn_eta>0) THEN
         DO iz=1,nz
            DO ix=1,nx
               DO jz=0,lambda
                  DO jx=0,lambda
                     READ(12) tfrbc(iz,ix,jz,jx,n,dir)
                  END DO
                END DO
            END DO
         END DO
      ELSE
         DO iz=1,nz
            DO ix=1,nx
               DO jz=0,lambda
                  DO jx=0,lambda
                     READ(12) tbabc(iz,ix,jz,jx,n,dir)
                  END DO
                END DO
            END DO
         END DO
      END IF
      ! >> left/right
      IF(sgn_eta>0 .and. sgn_xi>0) THEN
         dir=1
      ELSE IF (sgn_eta>0 .and. sgn_xi<0) THEN
         dir=2
      ELSE IF (sgn_eta<0 .and. sgn_xi>0) THEN
         dir=3
      ELSE
         dir=4
      END IF
      !
      IF (sgn_mu>0) THEN
        DO iy=1,ny
          DO iz=1,nz
            DO jy=0,lambda
              DO jz=0,lambda
                READ(12) tlebc(iy,iz,jy,jz,n,dir)
              END DO
            END DO
          END DO
        END DO
      ELSE
        DO iy=1,ny
          DO iz=1,nz
            DO jy=0,lambda
              DO jz=0,lambda
                READ(12) tribc(iy,iz,jy,jz,n,dir)
              END DO
            END DO
          END DO
        END DO
      END IF
      ! >> bottom/top
      IF(sgn_mu>0 .and. sgn_eta>0) THEN
         dir=1
      ELSE IF (sgn_mu>0 .and. sgn_eta<0) THEN
         dir=2
      ELSE IF (sgn_mu<0 .and. sgn_eta>0) THEN
         dir=3
      ELSE
         dir=4
      END IF
      !
      IF (sgn_xi>0) THEN
        DO ix=1,nx
          DO iy=1,ny
            DO jx=0,lambda
              DO jy=0,lambda
                READ(12) tbobc(ix,iy,jx,jy,n,dir)
              END DO
            END DO
          END DO
        END DO
      ELSE
        DO ix=1,nx
          DO iy=1,ny
            DO jx=0,lambda
              DO jy=0,lambda
                READ(12) ttobc(ix,iy,jx,jy,n,dir)
              END DO
            END DO
          END DO
        END DO
      END IF
      !
      !
   END DO
END DO

CLOSE(UNIT=12)

      !
      ! translate standard bc format to lagrange-type dgfem format.
      do dir=1,4
         do n=1,apo
            !
            ! front/back
            !
            ! space
            ! >> front/back
            if      (dir.eq.1) then
               sgz=1.0d0
               sgx=1.0d0
            else if (dir.eq.2) then
               sgz= 1.0d0
               sgx=-1.0d0
            else if (dir.eq.3) then
               sgz=-1.0d0
               sgx= 1.0d0
            else if (dir.eq.4) then
               sgz=-1.0d0
               sgx=-1.0d0
            end if
            do ix=1,nx
               do iz=1,nz
                  ! moments
                  do jx=0,lambda
                    do jy=0,lambda
                      do jz=0,lambda
                        IF (solvertype == "LD" .or. solvertype == "DENSE") THEN
                         IF(jx+jy+jz .le. lambda) THEN
                           l=jz+1-jy*(-3+2*jx+jy-2*lambda)/2+jx*(11+jx**2-3*jx*(2+lambda)+3*lambda*(4+lambda))/6
                           frbc(l,iz,ix,n,dir,1)=-ang(n,2)*(-1.0d0)**jy*dx(ix)*dz(iz)*sgz**jz*sgx**jx*tfrbc(iz,ix,jz,jx,n,dir)
                           babc(l,iz,ix,n,dir,1)=-ang(n,2)*(-1.0d0)**jy*dx(ix)*dz(iz)*sgz**jz*sgx**jx*tbabc(iz,ix,jz,jx,n,dir)
                         END IF
                        ELSE IF (solvertype == "LAGRANGE") THEN
                         l=jz+1+(lambda+1)*jy+(lambda+1)**2*jx
                         frbc(l,iz,ix,n,dir,1)=-ang(n,2)*(-1.0d0)**jy*dx(ix)*dz(iz)*sgz**jz*sgx**jx*tfrbc(iz,ix,jz,jx,n,dir)
                         babc(l,iz,ix,n,dir,1)=-ang(n,2)*(-1.0d0)**jy*dx(ix)*dz(iz)*sgz**jz*sgx**jx*tbabc(iz,ix,jz,jx,n,dir)
                        END IF
                      end do
                    end do
                  end do
                  !
               end do
            end do
            !
            ! left/right
            !
            ! space
            if      (dir.eq.1) then
               sgy=1.0d0
               sgz=1.0d0
            else if (dir.eq.2) then
               sgy= 1.0d0
               sgz=-1.0d0
            else if (dir.eq.3) then
               sgy=-1.0d0
               sgz= 1.0d0
            else if (dir.eq.4) then
               sgy=-1.0d0
               sgz=-1.0d0
            end if
            do iz=1,nz
               do iy=1,ny
                  ! moments
                  do jx=0,lambda
                    do jy=0,lambda
                      do jz=0,lambda
                        IF (solvertype == "LD" .or. solvertype == "DENSE") THEN
                         IF(jx+jy+jz .le. lambda) THEN
                           l=jz+1-jy*(-3+2*jx+jy-2*lambda)/2+jx*(11+jx**2-3*jx*(2+lambda)+3*lambda*(4+lambda))/6
                           lebc(l,iy,iz,n,dir,1)=-ang(n,1)*(-1.0d0)**jx*dy(iy)*dz(iz)*sgy**jy*sgz**jz*tlebc(iy,iz,jy,jz,n,dir)
                           ribc(l,iy,iz,n,dir,1)=-ang(n,1)*(-1.0d0)**jx*dy(iy)*dz(iz)*sgy**jy*sgz**jz*tribc(iy,iz,jy,jz,n,dir)
                         END IF
                        ELSE IF (solvertype == "LAGRANGE") THEN
                         l=jz+1+(lambda+1)*jy+(lambda+1)**2*jx
                         lebc(l,iy,iz,n,dir,1)=-ang(n,1)*(-1.0d0)**jx*dy(iy)*dz(iz)*sgy**jy*sgz**jz*tlebc(iy,iz,jy,jz,n,dir)
                         ribc(l,iy,iz,n,dir,1)=-ang(n,1)*(-1.0d0)**jx*dy(iy)*dz(iz)*sgy**jy*sgz**jz*tribc(iy,iz,jy,jz,n,dir)
                        END IF
                      end do
                    end do
                  end do
                  !
               end do
            end do
            !
            ! bottom/top
            !
            ! space
            if      (dir.eq.1) then
               sgx=1.0d0
               sgy=1.0d0
            else if (dir.eq.2) then
               sgx= 1.0d0
               sgy=-1.0d0
            else if (dir.eq.3) then
               sgx=-1.0d0
               sgy= 1.0d0
            else if (dir.eq.4) then
               sgx=-1.0d0
               sgy=-1.0d0
            end if
            do iy=1,ny
               do ix=1,nx
                  ! moments
                  do jx=0,lambda
                    do jy=0,lambda
                      do jz=0,lambda
                        IF (solvertype == "LD" .or. solvertype == "DENSE") THEN
                         IF(jx+jy+jz .le. lambda) THEN
                           l=jz+1-jy*(-3+2*jx+jy-2*lambda)/2+jx*(11+jx**2-3*jx*(2+lambda)+3*lambda*(4+lambda))/6
                           tobc(l,ix,iy,n,dir,1)=-ang(n,3)*(-1.0d0)**jz*dx(ix)*dy(iy)*sgx**jx*sgy**jy*ttobc(ix,iy,jx,jy,n,dir)
                           bobc(l,ix,iy,n,dir,1)=-ang(n,3)*(-1.0d0)**jz*dx(ix)*dy(iy)*sgx**jx*sgy**jy*tbobc(ix,iy,jx,jy,n,dir)
                         END IF
                        ELSE IF (solvertype == "LAGRANGE") THEN
                         l=jz+1+(lambda+1)*jy+(lambda+1)**2*jx
                         tobc(l,ix,iy,n,dir,1)=-ang(n,3)*(-1.0d0)**jz*dx(ix)*dy(iy)*sgx**jx*sgy**jy*ttobc(ix,iy,jx,jy,n,dir)
                         bobc(l,ix,iy,n,dir,1)=-ang(n,3)*(-1.0d0)**jz*dx(ix)*dy(iy)*sgx**jx*sgy**jy*tbobc(ix,iy,jx,jy,n,dir)
                        END IF
                      end do
                    end do
                  end do
                  !
               end do
            end do
            !
         end do
      end do
      !
      ! deallocate temporary bc arrays
      !
      deallocate( tfrbc,tbabc,tlebc,tribc,tbobc,ttobc )
      !
END SUBROUTINE read_inflow_dgfem


text:
!-----------------------------------------------------------------------
      program test

      include 'ceedf.h'

      integer ceed,err
      integer x,n
      real*8 a(10)
      real*8 b(10)
      real*8 diff
      integer*8 aoffset,boffset
      character arg*32

      call getarg(1,arg)

      call ceedinit(trim(arg)//char(0),ceed,err)

      n=10

      call ceedvectorcreate(ceed,n,x,err)

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

      aoffset=0
      call ceedvectorsetarray(x,ceed_mem_host,ceed_use_pointer,a,aoffset,err)
      call ceedvectorgetarray(x,ceed_mem_host,b,boffset,err)
      b(boffset+3) = -3.14
      call ceedvectorrestorearray(x,b,boffset,err)
      diff=a(3)+3.14
      if (abs(diff)>1.0D-15) then
! LCOV_EXCL_START
        write(*,*) 'Error writing array a(3)=',a(3)
! LCOV_EXCL_STOP
      endif

      call ceedvectordestroy(x,err)
      call ceeddestroy(ceed,err)

      end
!-----------------------------------------------------------------------

max_stars_repo_name: jrwrigh/libCEED

max_issues_repo_name: jrwrigh/libCEED

max_forks_repo_name: jrwrigh/libCEED


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1993  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     #############################################################
c     ##                                                         ##
c     ##  subroutine esolv2  --  atom-by-atom solvation Hessian  ##
c     ##                                                         ##
c     #############################################################
c
c
c     "esolv2" calculates second derivatives of the implicit
c     solvation energy for surface area, generalized Born,
c     generalized Kirkwood and Poisson-Boltzmann solvation models
c
c     note this version does not contain the chain rule terms
c     for derivatives of Born radii with respect to coordinates
c
c
      subroutine esolv2 (i)
      implicit none
      include 'sizes.i'
      include 'potent.i'
      include 'solute.i'
      include 'warp.i'
      integer i
      real*8 probe
c     real*8, allocatable :: aes(:)
c     real*8, allocatable :: des(:,:)
c
c
c     set a value for the solvent molecule probe radius
c
      probe = 1.4d0
c
c     perform dynamic allocation of some local arrays
c
c     allocate (aes(n))
c     allocate (des(3,n))
c
c     compute the surface area-based solvation energy term
c
c     call surface1 (es,aes,des,rsolv,asolv,probe)
c
c     perform deallocation of some local arrays
c
c     deallocate (aes)
c     deallocate (des)
c
c     get the generalized Born term for GB/SA solvation
c
      if (use_born .and. solvtyp(1:2).ne.'GK') then
         if (use_smooth) then
            call egb2b (i)
         else
            call egb2a (i)
         end if
      end if
      return
      end
c
c
c     ###############################################################
c     ##                                                           ##
c     ##  subroutine egb2a  --  atom-by-atom GB solvation Hessian  ##
c     ##                                                           ##
c     ###############################################################
c
c
c     "egb2a" calculates second derivatives of the generalized
c     Born energy term for the GB/SA solvation models
c
c
      subroutine egb2a (i)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'hessn.i'
      include 'shunt.i'
      include 'solute.i'
      integer i,j,k,kk
      real*8 e,de,d2e
      real*8 fi,fik
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 r,r2,r3,r4
      real*8 r5,r6,r7
      real*8 dwater,rb2,rm2
      real*8 expterm,shift
      real*8 d2edx,d2edy,d2edz
      real*8 taper,dtaper,d2taper
      real*8 trans,dtrans,d2trans
      real*8 fgb,fgb2,dfgb
      real*8 dfgb2,d2fgb
      real*8 term(3,3)
      character*6 mode
c
c
c     first see if the atom of interest carries a charge
c
      do k = 1, nion
         if (iion(k) .eq. i) then
            fi = pchg(k)
            goto 10
         end if
      end do
      return
   10 continue
c
c     store the coordinates of the atom of interest
c
      xi = x(i)
      yi = y(i)
      zi = z(i)
c
c     set the solvent dielectric and energy conversion factor
c
      dwater = 78.3d0
      fi = -electric * (1.0d0 - 1.0d0/dwater) * fi
c
c     set cutoff distances and switching function coefficients
c
      mode = 'CHARGE'
      call switch (mode)
c
c     calculate GB polarization energy Hessian elements
c
      do kk = 1, nion
         k = iion(kk)
         if (i .ne. k) then
            xr = xi - x(k)
            yr = yi - y(k)
            zr = zi - z(k)
            r2 = xr*xr + yr*yr + zr*zr
            if (r2 .le. off2) then
               r = sqrt(r2)
               fik = fi * pchg(kk)
c
c     compute chain rule terms for Hessian matrix elements
c
               rb2 = rborn(i) * rborn(k)
               expterm = exp(-0.25d0*r2/rb2)
               fgb2 = r2 + rb2*expterm
               fgb = sqrt(fgb2)
               dfgb = (1.0d0-0.25d0*expterm) * r / fgb
               dfgb2 = dfgb * dfgb
               d2fgb = -dfgb2/fgb + dfgb/r
     &                    + 0.125d0*(r2/rb2)*expterm/fgb
               de = -fik * dfgb / fgb2
               d2e = -fik * (d2fgb-2.0d0*dfgb2/fgb) / fgb2
c
c     use energy switching if near the cutoff distance
c
               if (r2 .gt. cut2) then
                  e = fik / fgb
                  rm2 = (0.5d0 * (off+cut))**2
                  shift = fik / sqrt(rm2 + rb2*exp(-0.25d0*rm2/rb2))
                  e = e - shift
                  r3 = r2 * r
                  r4 = r2 * r2
                  r5 = r2 * r3
                  r6 = r3 * r3
                  r7 = r3 * r4
                  taper = c5*r5 + c4*r4 + c3*r3 + c2*r2 + c1*r + c0
                  dtaper = 5.0d0*c5*r4 + 4.0d0*c4*r3
     &                        + 3.0d0*c3*r2 + 2.0d0*c2*r + c1
                  d2taper = 20.0d0*c5*r3 + 12.0d0*c4*r2
     &                         + 6.0d0*c3*r + 2.0d0*c2
                  trans = fik * (f7*r7 + f6*r6 + f5*r5 + f4*r4
     &                            + f3*r3 + f2*r2 + f1*r + f0)
                  dtrans = fik * (7.0d0*f7*r6 + 6.0d0*f6*r5
     &                            + 5.0d0*f5*r4 + 4.0d0*f4*r3
     &                            + 3.0d0*f3*r2 + 2.0d0*f2*r + f1)
                  d2trans = fik * (42.0d0*f7*r5 + 30.0d0*f6*r4
     &                             + 20.0d0*f5*r3 + 12.0d0*f4*r2
     &                             + 6.0d0*f3*r + 2.0d0*f2)
                  d2e = e*d2taper + 2.0d0*de*dtaper
     &                     + d2e*taper + d2trans
                  de = e*dtaper + de*taper + dtrans
               end if
c
c     form the individual Hessian element components
c
               de = de / r
               d2e = (d2e-de) / r2
               d2edx = d2e * xr
               d2edy = d2e * yr
               d2edz = d2e * zr
               term(1,1) = d2edx*xr + de
               term(1,2) = d2edx*yr
               term(1,3) = d2edx*zr
               term(2,1) = term(1,2)
               term(2,2) = d2edy*yr + de
               term(2,3) = d2edy*zr
               term(3,1) = term(1,3)
               term(3,2) = term(2,3)
               term(3,3) = d2edz*zr + de
c
c     increment diagonal and non-diagonal Hessian elements
c
               do j = 1, 3
                  hessx(j,i) = hessx(j,i) + term(1,j)
                  hessy(j,i) = hessy(j,i) + term(2,j)
                  hessz(j,i) = hessz(j,i) + term(3,j)
                  hessx(j,k) = hessx(j,k) - term(1,j)
                  hessy(j,k) = hessy(j,k) - term(2,j)
                  hessz(j,k) = hessz(j,k) - term(3,j)
               end do
            end if
         end if
      end do
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine egb2b  --  GB solvation Hessian for smoothing  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "egb2b" calculates second derivatives of the generalized
c     Born energy term for the GB/SA solvation models for use with
c     potential smoothing methods
c
c
      subroutine egb2b (i)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'charge.i'
      include 'chgpot.i'
      include 'hessn.i'
      include 'math.i'
      include 'solute.i'
      include 'warp.i'
      integer i,j,k,kk
      real*8 fi,fik,de,d2e
      real*8 xi,yi,zi
      real*8 xr,yr,zr
      real*8 dwater,width
      real*8 r,r2,rb2
      real*8 fgb,fgb2
      real*8 dfgb,dfgb2,d2fgb
      real*8 d2edx,d2edy,d2edz
      real*8 sterm,expterm
      real*8 erf,erfterm
      real*8 term(3,3)
      external erf
c
c
c     first see if the atom of interest carries a charge
c
      do k = 1, nion
         if (iion(k) .eq. i) then
            fi = pchg(k)
            goto 10
         end if
      end do
      return
   10 continue
c
c     store the coordinates of the atom of interest
c
      xi = x(i)
      yi = y(i)
      zi = z(i)
c
c     set the solvent dielectric and energy conversion factor
c
      dwater = 78.3d0
      fi = -electric * (1.0d0 - 1.0d0/dwater) * fi
c
c     set the extent of smoothing to be performed
c
      sterm = 0.5d0 / sqrt(diffc)
c
c     calculate GB polarization energy Hessian elements
c
      do kk = 1, nion
         k = iion(kk)
         if (i .ne. k) then
            xr = xi - x(k)
            yr = yi - y(k)
            zr = zi - z(k)
            r2 = xr*xr + yr*yr + zr*zr
            r = sqrt(r2)
            fik = fi * pchg(kk)
c
c     compute chain rule terms for Hessian matrix elements
c
            rb2 = rborn(i) * rborn(k)
            expterm = exp(-0.25d0*r2/rb2)
            fgb2 = r2 + rb2*expterm
            fgb = sqrt(fgb2)
            dfgb = (1.0d0-0.25d0*expterm) * r / fgb
            dfgb2 = dfgb * dfgb
            d2fgb = -dfgb2/fgb + dfgb/r
     &                 + 0.125d0*(r2/rb2)*expterm/fgb
            de = -fik * dfgb / fgb2
            d2e = -fik * (d2fgb-2.0d0*dfgb2/fgb) / fgb2
c
c     use a smoothable GB analogous to the Coulomb solution
c
            if (deform .gt. 0.0d0) then
               width = deform + 0.15d0*rb2*exp(-0.006d0*rb2/deform)
               width = sterm / sqrt(width)
               erfterm = erf(width*fgb)
               expterm = width * exp(-(width*fgb)**2) / sqrtpi
               de = de * (erfterm-2.0d0*expterm*fgb)
               d2e = d2e*erfterm + 2.0d0*fik*expterm
     &                  * (d2fgb/fgb-2.0d0*dfgb2*(1.0d0/fgb2+width**2))
            end if
c
c     form the individual Hessian element components
c
            de = de / r
            d2e = (d2e-de) / r2
            d2edx = d2e * xr
            d2edy = d2e * yr
            d2edz = d2e * zr
            term(1,1) = d2edx*xr + de
            term(1,2) = d2edx*yr
            term(1,3) = d2edx*zr
            term(2,1) = term(1,2)
            term(2,2) = d2edy*yr + de
            term(2,3) = d2edy*zr
            term(3,1) = term(1,3)
            term(3,2) = term(2,3)
            term(3,3) = d2edz*zr + de
c
c     increment diagonal and non-diagonal Hessian elements
c
            do j = 1, 3
               hessx(j,i) = hessx(j,i) + term(1,j)
               hessy(j,i) = hessy(j,i) + term(2,j)
               hessz(j,i) = hessz(j,i) + term(3,j)
               hessx(j,k) = hessx(j,k) - term(1,j)
               hessy(j,k) = hessy(j,k) - term(2,j)
               hessz(j,k) = hessz(j,k) - term(3,j)
            end do
         end if
      end do
      return
      end

ext: f

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


subroutine dsyr2m (uplo,transa,transb,m,n, alpha,a,lda,b,ldb,beta,c,ldc)
  !
  !  Purpose
  !  =======
  !
  !  dsyr2m performs one of the matrix-matrix operations
  !
  !     op1(C) := alpha*op2( A )*op2( B ) + beta*op1(C),
  !
  !  where  op1( X ) means that only the UPPER or only the LOWER
  !  part of X is accessed
  !
  !  and  op2( X ) is one of
  !
  !     op2( X ) = X   or   op2( X ) = X',
  !
  !  alpha and beta are scalars, and A, B and C are matrices, with
  !  op2( A ) an m by n matrix,
  !  op2( B ) an n by m matrix, and
  !  C an m by m matrix.
  !
  !     uplo: 'U'pper or 'L'ower
  !     transa,transb: op2( X ) = X   or   op2( X ) = X', with X=A,B
  !     m: number of rows of op2(A) and op2(B) and of
  !        rows and columns of C
  !     n: number of columns of op2(A) and op2(B)
  !     alpha
  !     op2(A): mXn matrix
  !     lda: leading dimension of a
  !     op2(B): mXn matrix
  !     ldb: leading dimension of b
  !     beta
  !     c: symmetric mXm matrix,
  !        stored either above (uplo='U') or below (uplo='L') the
  !        diagonal
  !     ldc: leading dimension of c
  !
  implicit none
  character*1 uplo,transa,transb
  integer m,n,lda,ldb,ldc
  double precision alpha,beta
  double precision a(lda,*), b(ldb,*), c(ldc,m)

  integer nblock
  parameter (nblock=128)
  double precision one,zero
  parameter (one=1.0d0, zero=0.0d0)

  integer j,ii,jj,imax,jmax
  logical upper,atrans,btrans

  logical lsame
  external lsame

  upper=lsame(uplo, 'U')
  atrans=lsame(transa, 'T')
  btrans=lsame(transb, 'T')

  if(.not.upper .and. .not.atrans .and. btrans) then
     do ii=1,m,nblock
        imax=min(m,ii+nblock-1)
        do j=ii,imax
           call dgemv('n', imax-j+1, n, alpha, a(j,1), lda, b(j,1), ldb, beta, c(j,j), 1)
        end do
        if(imax.lt.m) call dgemm('n','t',m-imax,imax-ii+1,n, alpha,a(imax+1,1),lda,b(ii,1),ldb, beta,c(imax+1,ii),ldc)
     end do
  else if(.not.upper .and. atrans .and. .not.btrans) then
     do ii=1,m,nblock
        imax=min(m,ii+nblock-1)
        do j=ii,imax
           call dgemv('t', n, imax-j+1, alpha, a(1,j), lda, b(1,j), 1, beta, c(j,j), 1)
        end do
        if(imax.lt.m) &
             call dgemm('t','n',m-imax,imax-ii+1,n,alpha,a(1,imax+1),lda,b(1,ii),ldb,beta,c(imax+1,ii),ldc)
     end do
  else if(upper .and. atrans .and. .not.btrans) then
     do jj=1,m,nblock
        jmax=min(m,jj+nblock-1)
        do j=jj,jmax
           call dgemv('t', n, j-jj+1, alpha, a(1,jj), lda,b(1,j), 1, beta, c(jj,j), 1)
        end do
        if(jj.ne.1) call dgemm('t','n',jj-1,jmax-jj+1,n,alpha,a(1,1),lda,b(1,jj),ldb,beta,c(1,jj),ldc)
     end do
  else
     stop 'wrong parameter in dsyr2m'
  end if
end subroutine dsyr2m


text:
! RUN: %python %S/test_folding.py %s %flang_fc1
! Tests folding of UNPACK (valid cases)
module m
  integer, parameter :: vector(*) = [1, 2, 3, 4]
  integer, parameter :: field(2,3) = reshape([(-j,j=1,6)], shape(field))
  logical, parameter :: mask(*,*) = reshape([.false., .true., .true., .false., .false., .true.], shape(field))
  logical, parameter :: test_unpack_1 = all(unpack(vector, mask, 0) == reshape([0,1,2,0,0,3], shape(mask)))
  logical, parameter :: test_unpack_2 = all(unpack(vector, mask, field) == reshape([-1,1,2,-4,-5,3], shape(mask)))
end module

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm