text:
c
c
c    common routines for generating and processing
c     a level restricted quad tree in 2D
c
c
c
c
      subroutine tree_refine_boxes(irefinebox,nboxes,
     1  ifirstbox,nbloc,centers,bs,nbctr,nlctr,
     2  ilevel,iparent,nchild,ichild)
      implicit none
      integer nboxes,nbloc,nbctr,nlctr
      real *8 centers(3,nboxes),bs
      integer ilevel(nboxes),iparent(nboxes)
      integer ichild(8,nboxes),nchild(nboxes)
      integer irefinebox(nbloc)
      integer ifirstbox
      integer, allocatable :: isum(:)
      integer ii

      integer i,ibox,nel0,j,l,jbox,nel1,nbl,jj

      allocate(isum(nbloc))
      if(nbloc.gt.0) call cumsum(nbloc,irefinebox,isum)

C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,ibox,nbl,j,jbox,ii,jj,l)
      do i = 1,nbloc
        ibox = ifirstbox + i-1
        if(irefinebox(i).eq.1) then
          nbl = nbctr + (isum(i)-1)*8

          nchild(ibox) = 8
          do j=1,8
            ii = 2
            jj = 2
            if(j.eq.1.or.j.eq.2.or.j.eq.5.or.j.eq.6) ii = 1
            if(j.lt.5) jj = 1
            jbox = nbl+j
            centers(1,jbox) = centers(1,ibox)+(-1)**j*bs/2
            centers(2,jbox) = centers(2,ibox)+(-1)**ii*bs/2
            centers(3,jbox) = centers(3,ibox)+(-1)**jj*bs/2
            iparent(jbox) = ibox
            nchild(jbox) = 0
            do l=1,8
              ichild(l,jbox) = -1
            enddo
            ichild(j,ibox) = jbox
            ilevel(jbox) = nlctr
          enddo
        endif
      enddo
C$OMP END PARALLEL DO

      if(nbloc.gt.0) nbctr = nbctr + isum(nbloc)*8

      return
      end
c
c
c
c
c
c
       subroutine tree_copy(nb,centers,ilevel,iparent,nchild,ichild,
     1              centers2,ilevel2,iparent2,nchild2,ichild2)

       implicit none
       integer nd,nb,npb
       real *8 centers(3,nb),centers2(3,nb)
       integer ilevel(nb),ilevel2(nb)
       integer iparent(nb),iparent2(nb)
       integer nchild(nb),nchild2(nb)
       integer ichild(8,nb),ichild2(8,nb)

       integer i,j,nel

C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j)
       do i=1,nb
         centers2(1,i) = centers(1,i)
         centers2(2,i) = centers(2,i)
         centers2(3,i) = centers(3,i)
         ilevel2(i) = ilevel(i)
         iparent2(i) = iparent(i)
         nchild2(i) = nchild(i)
         ichild2(1,i) = ichild(1,i)
         ichild2(2,i) = ichild(2,i)
         ichild2(3,i) = ichild(3,i)
         ichild2(4,i) = ichild(4,i)
         ichild2(5,i) = ichild(5,i)
         ichild2(6,i) = ichild(6,i)
         ichild2(7,i) = ichild(7,i)
         ichild2(8,i) = ichild(8,i)
       enddo
C$OMP END PARALLEL DO

       return
       end
c
c
c
c
c

      subroutine computecoll(nlevels,nboxes,laddr,boxsize,
     1                       centers,iparent,nchild,ichild,
     2                       iper,nnbors,nbors)

c     This subroutine computes the colleagues for an adaptive
c     pruned tree. box j is a colleague of box i, if they share a
c     vertex or an edge and the two boxes are at the same
c     level in the tree
c
c     INPUT arguments
c     nlevels     in: integer
c                 Number of levels
c
c     nboxes      in: integer
c                 Total number of boxes
c
c     laddr       in: integer(2,0:nlevels)
c                 indexing array providing access to boxes at
c                 each level.
c                 the first box on level i is laddr(1,i)
c                 the last box on level i is laddr(2,i)
c
c     boxsize     in: double precision(0:nlevels)
c                 Array of boxsizes
c
c     centers     in: double precision(3,nboxes)
c                 array of centers of boxes
c
c     iparent     in: integer(nboxes)
c                 iparent(i) is the box number of the parent of
c                 box i
c
c     nchild      in: integer(nboxes)
c                 nchild(i) is the number of children of box i
c
c     ichild      in: integer(8,nboxes)
c                 ichild(j,i) is the box id of the jth child of
c                 box i
c
c     iper        in: integer
c                 flag for periodic implementations.
c                 Currently not used. Feature under construction.
c
c----------------------------------------------------------------
c     OUTPUT
c     nnbors      out: integer(nboxes)
c                 nnbors(i) is the number of colleague boxes of
c                 box i
c
c     nbors       out: integer(27,nboxes)
c                 nbors(j,i) is the box id of the jth colleague
c                 box of box i
c---------------------------------------------------------------
      implicit none
      integer nlevels,nboxes,iper
      integer laddr(2,0:nlevels)
      double precision boxsize(0:nlevels)
      double precision centers(3,nboxes)
      integer iparent(nboxes), nchild(nboxes), ichild(8,nboxes)
      integer nnbors(nboxes)
      integer nbors(27,nboxes)

c     Temp variables
      integer ilev,ibox,jbox,kbox,dad
      integer i,j,ifirstbox,ilastbox

c     Setting parameters for level = 0
      nnbors(1) = 1
      nbors(1,1) = 1
      do ilev = 1,nlevels
c        Find the first and the last box at level ilev
         ifirstbox = laddr(1,ilev)
         ilastbox = laddr(2,ilev)
c        Loop over all boxes to evaluate neighbors, list1 and updating
c        hunglists of targets

C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP$PRIVATE(ibox,dad,i,jbox,j,kbox)
         do ibox = ifirstbox,ilastbox
c           Find the parent of the current box
            dad = iparent(ibox)
c           Loop over the neighbors of the parent box
c           to find out list 1 and list 2
            do i=1,nnbors(dad)
                jbox = nbors(i,dad)
                do j=1,8
c               ichild(j,jbox) is one of the children of the
c               neighbors of the parent of the current
c               box
                   kbox = ichild(j,jbox)
                   if(kbox.gt.0) then
c               Check if kbox is a nearest neighbor or in list 2
                      if((abs(centers(1,kbox)-centers(1,ibox)).le.
     1                   1.05*boxsize(ilev)).and.
     2                   (abs(centers(2,kbox)-centers(2,ibox)).le.
     3                   1.05*boxsize(ilev)).and.
     2                   (abs(centers(3,kbox)-centers(3,ibox)).le.
     3                   1.05*boxsize(ilev))) then

                         nnbors(ibox) = nnbors(ibox)+1
                         nbors(nnbors(ibox),ibox) = kbox
                      endif
                   endif
                enddo
            enddo
c           End of computing colleagues of box i
         enddo
C$OMP END PARALLEL DO
      enddo

      return
      end
c
c
c
c
c--------------------------------------------------------------------
      subroutine updateflags(curlev,nboxes,nlevels,laddr,nchild,ichild,
     1                    nnbors,nbors,centers,boxsize,iflag)

c      This subroutine is to check the boxes flagged as flag++
c      and determine which of the boxes need refinement. The flag
c      of the box which need refinement is updated to iflag(box)=1
c      and that of the boxes which do not need refinement is
c      updated to iflag(box) = 0
c
c      INPUT arguments
c      curlev         in: integer
c                     the level for which boxes need to be processed
c
c      nboxes         in: integer
c                     total number of boxes
c
c      nlevels        in: integer
c                     total number of levels
c
c      laddr          in: integer(2,0:nlevels)
c                     boxes from laddr(1,ilev) to laddr(2,ilev)
c                     are at level ilev
c
c      nchild         in: integer(nboxes)
c                     nchild(ibox) is the number of children
c                     of box ibox
c
c      ichild         in: integer(8,nboxes)
c                     ichild(j,ibox) is the box id of the jth
c                     child of box ibox
c
c      nnbors         in: integer(nboxes)
c                     nnbors(ibox) is the number of colleagues
c                     of box ibox
c
c      nbors          in: integer(27,nboxes)
c                     nbors(j,ibox) is the jth colleague of box
c                     ibox
c
c      centers        in: double precision(2,nboxes)
c                     x and y coordinates of the box centers
c
c      boxsize        in: double precision(0:nlevels)
c                     boxsize(i) is the size of the box at level i
c
c      iflag          in/out: integer(nboxes)
c                     iflag(ibox)=3 if it is flag++. iflag(ibox) =1
c                     or 0 at the end of routine depending on
c                     whether box needs to be subdivided or not
c
      implicit none
c     Calling sequence variables
      integer curlev, nboxes, nlevels
      integer laddr(2,0:nlevels),nchild(nboxes),ichild(8,nboxes)
      integer nnbors(nboxes), nbors(27,nboxes)
      integer iflag(nboxes)
      double precision centers(3,nboxes),boxsize(0:nlevels)

c     Temporary variables
      integer i,j,k,l,ibox,jbox,kbox,lbox, ict
      double precision distest,xdis,ydis,zdis

      distest = 1.05d0*(boxsize(curlev) + boxsize(curlev+1))/2.0d0
c     Loop over all boxes at the current level

C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ibox,i,jbox,j,kbox,xdis,ydis)
C$OMP$PRIVATE(zdis,ict)
      do ibox = laddr(1,curlev),laddr(2,curlev)
         if(iflag(ibox).eq.3) then
            iflag(ibox) = 0
c           Loop over colleagues of the current box
            do i=1,nnbors(ibox)
c              Loop over colleagues of flag++ box
               jbox = nbors(i,ibox)

c              Loop over the children of the colleague box
c              Note we do not need to exclude self from
c              the list of colleagues as a self box which
c              is flag++ does not have any children
c              and will not enter the next loop
               do j=1,8
                  kbox = ichild(j,jbox)
                  if(kbox.gt.0) then
                     if(nchild(kbox).gt.0) then
                        xdis = centers(1,kbox) - centers(1,ibox)
                        ydis = centers(2,kbox) - centers(2,ibox)
                        zdis = centers(3,kbox) - centers(3,ibox)
                        ict = 0
                        if(abs(xdis).le.distest) ict = ict + 1
                        if(abs(ydis).le.distest) ict = ict + 1
                        if(abs(zdis).le.distest) ict = ict + 1
                        if(ict.eq.3) then
                           iflag(ibox) = 1
                           goto 1111
                        endif
                     endif
                  endif
c                 End of looping over the children of the child
c                 of the colleague box
               enddo
c              End of looping over the children of the colleague box
            enddo
c           End of looping over colleagues
 1111       continue
         endif
c        End of testing if the current box needs to checked for
      enddo
c     End of looping over boxes at the current level
C$OMP END PARALLEL DO

      return
      end
c
c
c
c
c
c

      subroutine tree_refine_boxes_flag(iflag,nboxes,
     1  ifirstbox,nbloc,centers,bs,nbctr,nlctr,
     2  ilevel,iparent,nchild,ichild)
      implicit none
      integer nboxes,nbloc,nbctr,nlctr
      real *8 centers(3,nboxes),bs
      integer ilevel(nboxes),iparent(nboxes)
      integer ichild(8,nboxes),nchild(nboxes)
      integer iflag(nboxes)
      integer ifirstbox
      integer, allocatable :: isum(:),itmp(:)

      integer i,ibox,nel0,j,l,jbox,nel1,nbl
      integer ii,jj

      allocate(isum(nbloc),itmp(nbloc))

C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ibox,i)
      do i=1,nbloc
        ibox = ifirstbox+i-1
        itmp(i) = 0
        if(iflag(ibox).gt.0) itmp(i) = 1
      enddo
C$OMP END PARALLEL DO
      if(nbloc.gt.0) call cumsum(nbloc,itmp,isum)

C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,ibox,nbl,j,jbox,ii,jj,l)
      do i = 1,nbloc
        ibox = ifirstbox + i-1
        if(iflag(ibox).gt.0) then
          nbl = nbctr + (isum(i)-1)*8

          nchild(ibox) = 8
          do j=1,8
            ii = 2
            jj = 2
            if(j.eq.1.or.j.eq.2.or.j.eq.5.or.j.eq.6) ii = 1
            if(j.lt.5) jj = 1
            jbox = nbl+j
            centers(1,jbox) = centers(1,ibox)+(-1)**j*bs/2
            centers(2,jbox) = centers(2,ibox)+(-1)**ii*bs/2
            centers(3,jbox) = centers(3,ibox)+(-1)**jj*bs/2
            iparent(jbox) = ibox
            nchild(jbox) = 0
            do l=1,8
              ichild(l,jbox) = -1
            enddo
            ichild(j,ibox) = jbox
            ilevel(jbox) = nlctr+1
            if(iflag(ibox).eq.1) iflag(jbox) = 3
            if(iflag(ibox).eq.2) iflag(jbox) = 0
          enddo
        endif
      enddo
C$OMP END PARALLEL DO

      if(nbloc.gt.0) nbctr = nbctr + isum(nbloc)*8

      return
      end
c
c
c
c
c
      subroutine computemnlists(nlevels,nboxes,laddr,boxsize,
     1                   centers,iparent,nchild,
     2                   ichild,isep,nnbors,mnbors,nbors,iper,mnlist1,
     3                   mnlist2,mnlist3,mnlist4)
c     Compute max nuber of boxes in list1,list2,list3,list4
      implicit none
      integer nlevels,nboxes
      integer iper
      integer laddr(2,0:nlevels)
      double precision boxsize(0:nlevels)
      double precision centers(3,nboxes)
      integer iparent(nboxes),nchild(nboxes),ichild(8,nboxes)
      integer mnbors,isep
      integer nnbors(nboxes),nbors(mnbors,nboxes)
      integer mnlist1,mnlist2,mnlist3,mnlist4
      integer nlist1(nboxes),nlist2(nboxes),nlist3(nboxes)
      integer nlist4(nboxes)

c     Temp variables
      integer ilev,ibox,jbox,kbox,i,j,k,l
      integer firstbox,lastbox,dad
      double precision xdis,ydis,zdis,distest

C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
      do i=1,nboxes
         nlist1(i) = 0
         nlist2(i) = 0
         nlist3(i) = 0
         nlist4(i) = 0
      enddo
C$OMP END PARALLEL DO

      nlist1(1) = 1
      nlist2(1) = 0
      nlist3(1) = 0
      nlist4(1) = 0

      do ilev = 1,nlevels
         firstbox = laddr(1,ilev)
         lastbox = laddr(2,ilev)
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP$PRIVATE(ibox,dad,i,jbox,j,kbox,distest,xdis,ydis,zdis)
         do ibox = firstbox,lastbox
            dad = iparent(ibox)
            do i=1,nnbors(dad)
               jbox = nbors(i,dad)
               do j=1,8
                  kbox = ichild(j,jbox)
                  if(kbox.gt.0) then
                     if((abs(centers(1,kbox)-centers(1,ibox)).ge.
     1                  1.05d0*isep*boxsize(ilev)).or.
     2                  (abs(centers(2,kbox)-centers(2,ibox)).ge.
     3                  1.05d0*isep*boxsize(ilev)).or.
     4                  (abs(centers(3,kbox)-centers(3,ibox)).ge.
     5                  1.05d0*isep*boxsize(ilev))) then

                        nlist2(ibox) = nlist2(ibox) + 1
                     endif
                  endif
               enddo
            enddo
c           Compute list1 and list3 of ibox if it is childless
            if(nchild(ibox).eq.0) then
               do i=1,nnbors(ibox)
                  jbox = nbors(i,ibox)

c
cc                     check for list1 at the same level
c
                  if(nchild(jbox).eq.0) then
                     nlist1(ibox) = nlist1(ibox) + 1
                  endif
c
cc                     check for list1 and list3 at one ilev+1
                  if(nchild(jbox).gt.0) then
                     distest = 1.05d0*(boxsize(ilev)+boxsize(ilev+1))/
     1                         2.0d0*isep
                     do j=1,8
                        kbox = ichild(j,jbox)
                        if(kbox.gt.0) then
                           xdis = dabs(centers(1,kbox)-centers(1,ibox))
                           ydis = dabs(centers(2,kbox)-centers(2,ibox))
                           zdis = dabs(centers(3,kbox)-centers(3,ibox))

                           if(xdis.lt.distest.and.ydis.lt.distest.and.
     1                        zdis.lt.distest) then
                              nlist1(ibox) = nlist1(ibox)+1
                           else
                              nlist3(ibox) = nlist3(ibox)+1
                           endif
                        endif
                     enddo
                  endif
               enddo
c
cc               compute list1 and list4 for boxes at level ilev-1
               do i=1,nnbors(dad)
                   jbox = nbors(i,dad)
                   if(nchild(jbox).eq.0) then
                      distest = 1.05d0*(boxsize(ilev)+boxsize(ilev-1))/
     1                         2.0d0*isep
                      xdis = dabs(centers(1,jbox)-centers(1,ibox))
                      ydis = dabs(centers(2,jbox)-centers(2,ibox))
                      zdis = dabs(centers(3,jbox)-centers(3,ibox))
                      if(xdis.lt.distest.and.ydis.lt.distest.and.
     1                  zdis.lt.distest) then
                         nlist1(ibox) = nlist1(ibox)+1
                      endif
                   endif
               enddo
            endif
c
cc           compute list 4 at level ilev-1
c
            do i=1,nnbors(dad)
               jbox = nbors(i,dad)
               if(nchild(jbox).eq.0) then
                   distest = 1.05d0*(boxsize(ilev)+boxsize(ilev-1))/
     1                 2.0d0*isep
                    xdis = dabs(centers(1,jbox)-centers(1,ibox))
                    ydis = dabs(centers(2,jbox)-centers(2,ibox))
                    zdis = dabs(centers(3,jbox)-centers(3,ibox))
                    if(xdis.gt.distest.or.ydis.gt.distest.or.
     1                 zdis.gt.distest) then
                       nlist4(ibox) = nlist4(ibox)+1
                    endif
               endif
            enddo
         enddo
C$OMP END PARALLEL DO
      enddo

      mnlist1 = 0
      mnlist2 = 0
      mnlist3 = 0
      mnlist4 = 0
C$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i)
C$OMP$REDUCTION(max:mnlist1,mnlist2,mnlist3,mnlist4)
      do i=1,nboxes
         if(nlist1(i).gt.mnlist1) mnlist1 = nlist1(i)
         if(nlist2(i).gt.mnlist2) mnlist2 = nlist2(i)
         if(nlist3(i).gt.mnlist3) mnlist3 = nlist3(i)
         if(nlist4(i).gt.mnlist4) mnlist4 = nlist4(i)
      enddo
C$OMP END PARALLEL DO

      return
      end
c---------------------------------------------------------------
      subroutine computelists(nlevels,nboxes,laddr,boxsize,
     1                   centers,iparent,nchild,
     2                   ichild,isep,nnbors,mnbors,nbors,iper,nlist1,
     3                   mnlist1,list1,nlist2,mnlist2,list2,
     4                   nlist3,mnlist3,list3,nlist4,mnlist4,list4)
c     Compute max nuber of boxes in list1,list2,list3,list4
      implicit none
      integer nlevels,nboxes
      integer iper
      integer laddr(2,0:nlevels)
      double precision boxsize(0:nlevels)
      double precision centers(3,nboxes)
      integer iparent(nboxes),nchild(nboxes),ichild(8,nboxes)
      integer mnbors
      integer nnbors(nboxes),nbors(mnbors,nboxes)
      integer mnlist1,mnlist2,mnlist3,mnlist4,isep
      integer nlist1(nboxes),nlist2(nboxes),nlist3(nboxes)
      integer nlist4(nboxes)
      integer list1(mnlist1,nboxes),list2(mnlist2,nboxes)
      integer list3(mnlist3,nboxes),list4(mnlist4,nboxes)

c     Temp variables
      integer ilev,ibox,jbox,kbox,i,j,k,l
      integer firstbox,lastbox,dad
      double precision xdis,ydis,zdis,distest

C$OMP PARALLEL DO DEFAULT(SHARED)
      do i=1,nboxes
         nlist1(i) = 0
         nlist2(i) = 0
         nlist3(i) = 0
         nlist4(i) = 0
      enddo
C$OMP END PARALLEL DO
      if(nchild(1).eq.0) then
         nlist1(1) = 1
         list1(1,1) = 1
      else
         nlist1(1) = 0
      endif
      nlist2(1) = 0
      nlist3(1) = 0
      nlist4(1) = 0

      do ilev = 1,nlevels
         firstbox = laddr(1,ilev)
         lastbox = laddr(2,ilev)
C$OMP PARALLEL DO DEFAULT(SHARED)
C$OMP$PRIVATE(ibox,dad,i,jbox,j,kbox,xdis,ydis,zdis,distest)
         do ibox = firstbox,lastbox
            dad = iparent(ibox)
            do i=1,nnbors(dad)
               jbox = nbors(i,dad)
               do j=1,8
                  kbox = ichild(j,jbox)
                  if(kbox.gt.0) then
                     if((abs(centers(1,kbox)-centers(1,ibox)).ge.
     1                  1.05d0*isep*boxsize(ilev)).or.
     2                  (abs(centers(2,kbox)-centers(2,ibox)).ge.
     3                  1.05d0*isep*boxsize(ilev)).or.
     4                  (abs(centers(3,kbox)-centers(3,ibox)).ge.
     5                  1.05d0*isep*boxsize(ilev))) then

                        nlist2(ibox) = nlist2(ibox) + 1
                        list2(nlist2(ibox),ibox) = kbox
                     endif
                  endif
               enddo
            enddo
c           Compute list1 and list3 of ibox if it is childless
            if(nchild(ibox).eq.0) then
               do i=1,nnbors(ibox)
                  jbox = nbors(i,ibox)

c
cc                boxes in list 1 at the same level
c

                  if(nchild(jbox).eq.0) then
                     nlist1(ibox) = nlist1(ibox) + 1
                     list1(nlist1(ibox),ibox) = jbox
                  else
c
cc                     boxes in list1 and list3 at level ilev+1
c
                     distest = 1.05d0*(boxsize(ilev)+boxsize(ilev+1))/
     1                         2.0d0*isep
                     do j=1,8
                        kbox = ichild(j,jbox)
                        if(kbox.gt.0) then
                           xdis = dabs(centers(1,kbox)-centers(1,ibox))
                           ydis = dabs(centers(2,kbox)-centers(2,ibox))
                           zdis = dabs(centers(3,kbox)-centers(3,ibox))

                           if(xdis.lt.distest.and.ydis.lt.distest.and.
     1                        zdis.lt.distest) then
                              nlist1(ibox) = nlist1(ibox)+1
                              list1(nlist1(ibox),ibox) = kbox
                           else
                              nlist3(ibox) = nlist3(ibox)+1
                              list3(nlist3(ibox),ibox) = kbox
                           endif
                        endif
                     enddo
                  endif
               enddo
c
cc               compute list1 at level ilev-1
               do i=1,nnbors(dad)
                   jbox = nbors(i,dad)
                   if(nchild(jbox).eq.0) then
                      distest = 1.05d0*(boxsize(ilev)+boxsize(ilev-1))/
     1                         2.0d0*isep
                      xdis = dabs(centers(1,jbox)-centers(1,ibox))
                      ydis = dabs(centers(2,jbox)-centers(2,ibox))
                      zdis = dabs(centers(3,jbox)-centers(3,ibox))
                      if(xdis.lt.distest.and.ydis.lt.distest.and.
     1                  zdis.lt.distest) then
                         nlist1(ibox) = nlist1(ibox)+1
                         list1(nlist1(ibox),ibox) = jbox
                      endif
                   endif
                enddo
            endif
c
cc           compute list 4 at level ilev-1
c
            do i=1,nnbors(dad)
               jbox = nbors(i,dad)
               if(nchild(jbox).eq.0) then
                   distest = 1.05d0*(boxsize(ilev)+boxsize(ilev-1))/
     1                 2.0d0*isep
                    xdis = dabs(centers(1,jbox)-centers(1,ibox))
                    ydis = dabs(centers(2,jbox)-centers(2,ibox))
                    zdis = dabs(centers(3,jbox)-centers(3,ibox))
                    if(xdis.gt.distest.or.ydis.gt.distest.or.
     1                 zdis.gt.distest) then
                       nlist4(ibox) = nlist4(ibox)+1
                       list4(nlist4(ibox),ibox)=jbox
                    endif
               endif
            enddo
         enddo
C$OMP END PARALLEL DO
      enddo

      return
      end
c----------------------------------------------------------------

      subroutine getpwlistall(ibox,bs,nboxes,nnbors,nbors,
     1           nchild,ichild,centers,isep,nuall,uall,ndall,dall,nnall,
     2           nall,nsall,sall,neall,eall,nwall,wall,nu1234,u1234,
     3           nd5678,d5678,nn1256,n1256,ns3478,s3478,ne1357,e1357,
     4           nw2468,w2468,nn12,n12,nn56,n56,ns34,s34,ns78,s78,ne13,
     5           e13,ne57,e57,nw24,w24,nw68,w68,ne1,e1,ne3,e3,ne5,e5,
     6           ne7,e7,nw2,w2,nw4,w4,nw6,w6,nw8,w8)
c-------------------------------------------------------------------
      implicit none
      integer ibox
      double precision boxsize,bs
      integer nboxes,nnbors,nbors(nnbors)
      integer nchild, ichild(8,nboxes)
      double precision centers(3,nboxes)
      integer isep
      integer nuall,ndall,nnall,nsall,neall,nwall,nu1234
      integer nd5678,nn1256,ns3478,ne1357,nw2468
      integer nn12,nn56,ns34,ns78,ne13,ne57,nw24,nw68
      integer ne1,ne3,ne5,ne7,nw2,nw4,nw6,nw8
      integer uall(1),dall(1),nall(1),sall(1),eall(1),wall(1)
      integer u1234(1),d5678(1),n1256(1),s3478(1),e1357(1),w2468(1)
      integer n12(1),n56(1),s34(1),s78(1),e13(1),e57(1),w24(1),w68(1)
      integer e1(1),e3(1),e5(1),e7(1)
      integer w2(1),w4(1),w6(1),w8(1)

      integer jbox,kbox
      integer c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c16
      integer i,j

      nuall = 0
      ndall = 0
      nnall = 0
      nsall = 0
      neall = 0
      nwall = 0
      nu1234 = 0
      nd5678 = 0
      nn1256 = 0
      ns3478 = 0
      ne1357 = 0
      nw2468 = 0
      nn12 = 0
      nn56 = 0
      ns34 = 0
      ns78 = 0
      ne13 = 0
      ne57 = 0
      nw24 = 0
      nw68 = 0
      ne1 = 0
      ne3 = 0
      ne5 = 0
      ne7 = 0
      nw2 = 0
      nw4 = 0
      nw6 = 0
      nw8 = 0
      do i=1,nnbors
         jbox = nbors(i)
         do j=1,8
            kbox = ichild(j,jbox)
            if(kbox.gt.0) then
               c1 = 0
               c2 = 0
               c3 = 0
               c4 = 0
               c5 = 0
               c6 = 0
               c7 = 0
               c8 = 0
               c9 = 0
               c10 = 0
               c11 = 0
               c12 = 0
               if((centers(3,kbox)-centers(3,ibox)).ge.
     1              1.01d0*isep*bs+bs/2.0d0) c1 = 1
               if((centers(3,kbox)-centers(3,ibox)).le.
     1              -1.01d0*isep*bs-bs/2.0d0) c2 = 1
               if((centers(2,kbox)-centers(2,ibox)).ge.
     1              1.01d0*isep*bs+bs/2.0d0) c3 = 1
               if((centers(2,kbox)-centers(2,ibox)).le.
     1              -1.01d0*isep*bs-bs/2.0d0) c4 = 1
               if((centers(1,kbox)-centers(1,ibox)).ge.
     1              1.01d0*isep*bs+bs/2.0d0) c5 = 1
               if((centers(1,kbox)-centers(1,ibox)).le.
     1              -1.01d0*isep*bs-bs/2.0d0) c6 = 1
               if((centers(3,kbox)-centers(3,ibox)).ge.
     1              1.01d0*isep*bs-bs/2.0d0) c7 = 1
               if((centers(3,kbox)-centers(3,ibox)).le.
     1              -1.01d0*isep*bs+bs/2.0d0) c8 = 1
               if((centers(2,kbox)-centers(2,ibox)).ge.
     1              1.01d0*isep*bs-bs/2.0d0) c9 = 1
               if((centers(2,kbox)-centers(2,ibox)).le.
     1              -1.01d0*isep*bs+bs/2.0d0) c10 = 1
               if((centers(1,kbox)-centers(1,ibox)).ge.
     1              1.01d0*isep*bs-bs/2.0d0) c11 = 1
               if((centers(1,kbox)-centers(1,ibox)).le.
     1              -1.01d0*isep*bs+bs/2.0d0) c12 = 1
               if(c1.eq.1) then
                  nuall = nuall + 1
                  uall(nuall) = kbox
               endif

               if(c2.eq.1) then
                  ndall = ndall + 1
                  dall(ndall) = kbox
               endif

               if(c3.eq.1.and.c1.eq.0.and.c2.eq.0) then
                  nnall = nnall + 1
                  nall(nnall) = kbox
               endif

               if(c4.eq.1.and.c1.eq.0.and.c2.eq.0) then
                  nsall = nsall + 1
                  sall(nsall) = kbox
               endif

               if(c5.eq.1.and.c1.eq.0.and.c2.eq.0.and.c3.eq.0.and.
     1             c4.eq.0) then
                  neall = neall + 1
                  eall(neall) = kbox
               endif

               if(c6.eq.1.and.c1.eq.0.and.c2.eq.0.and.c3.eq.0.and.
     1            c4.eq.0) then
                  nwall = nwall + 1
                  wall(nwall) = kbox
               endif

               c16 = c1 + c2 + c3 + c4 + c5 +c6
               if(c16.eq.0.and.c7.eq.1) then
                  nu1234 = nu1234 + 1
                  u1234(nu1234) = kbox
               endif

               if(c16.eq.0.and.c8.eq.1) then
                  nd5678 = nd5678 + 1
                  d5678(nd5678) = kbox
               endif

               if(c16.eq.0.and.c9.eq.1.and.c7.eq.0.and.c8.eq.0) then
                  nn1256 = nn1256 + 1
                  n1256(nn1256) = kbox
               endif

               if(c16.eq.0.and.c10.eq.1.and.c7.eq.0.and.c8.eq.0) then
                  ns3478 = ns3478 + 1
                  s3478(ns3478) = kbox
               endif

               if(c16.eq.0.and.c11.eq.1.and.(c7+c8+c9+c10).eq.0) then
                  ne1357 = ne1357 + 1
                  e1357(ne1357) = kbox
               endif

               if(c16.eq.0.and.c12.eq.1.and.(c7+c8+c9+c10).eq.0) then
                  nw2468 = nw2468 + 1
                  w2468(nw2468) = kbox
               endif

               if(c16.eq.0.and.c8.eq.1.and.c9.eq.1) then
                  nn12 = nn12 + 1
                  n12(nn12) = kbox
               endif

               if(c16.eq.0.and.c7.eq.1.and.c9.eq.1) then
                  nn56 = nn56 + 1
                  n56(nn56) = kbox
               endif

               if(c16.eq.0.and.c8.eq.1.and.c10.eq.1) then
                  ns34 = ns34 + 1
                  s34(ns34) = kbox
               endif

               if(c16.eq.0.and.c7.eq.1.and.c10.eq.1) then
                  ns78 = ns78 + 1
                  s78(ns78) = kbox
               endif

               if(c16.eq.0.and.c8.eq.1.and.c11.eq.1
     1           .and.c9.eq.0.and.c10.eq.0) then
                  ne13 = ne13 + 1
                  e13(ne13) = kbox
               endif

               if(c16.eq.0.and.c7.eq.1.and.c11.eq.1
     1           .and.c9.eq.0.and.c10.eq.0) then
                  ne57 = ne57 + 1
                  e57(ne57) = kbox
               endif

               if(c16.eq.0.and.c8.eq.1.and.c12.eq.1
     1           .and.c9.eq.0.and.c10.eq.0) then
                  nw24 = nw24 + 1
                  w24(nw24) = kbox
               endif

               if(c16.eq.0.and.c7.eq.1.and.c12.eq.1
     1           .and.c9.eq.0.and.c10.eq.0) then
                  nw68 = nw68 + 1
                  w68(nw68) = kbox
               endif

               if(c16.eq.0.and.c7.eq.0.and.c10.eq.1.and.c11.eq.1) then
                  ne1 = ne1 + 1
                  e1(ne1) = kbox
               endif
               if(c16.eq.0.and.c7.eq.0.and.c9.eq.1.and.c11.eq.1) then
                  ne3 = ne3 + 1
                  e3(ne3) = kbox
               endif
               if(c16.eq.0.and.c8.eq.0.and.c10.eq.1.and.c11.eq.1) then
                  ne5 = ne5 + 1
                  e5(ne5) = kbox
               endif
               if(c16.eq.0.and.c8.eq.0.and.c9.eq.1.and.c11.eq.1) then
                  ne7 = ne7 + 1
                  e7(ne7) = kbox
               endif
               if(c16.eq.0.and.c7.eq.0.and.c10.eq.1.and.c12.eq.1) then
                  nw2 = nw2 + 1
                  w2(nw2) = kbox
               endif
               if(c16.eq.0.and.c7.eq.0.and.c9.eq.1.and.c12.eq.1) then
                  nw4 = nw4 + 1
                  w4(nw4) = kbox
               endif
               if(c16.eq.0.and.c8.eq.0.and.c10.eq.1.and.c12.eq.1) then
                  nw6 = nw6 + 1
                  w6(nw6) = kbox
               endif
               if(c16.eq.0.and.c8.eq.0.and.c9.eq.1.and.c12.eq.1) then
                  nw8 = nw8 + 1
                  w8(nw8) = kbox
               endif
            endif
         enddo
      enddo

      return
      end
c--------------------------------------------------------------------

      subroutine getlist3pwlistall(ibox,bs,nboxes,nlist3,list3,
     1           isep,centers,nuall,uall,ndall,dall,nnall,
     2           nall,nsall,sall,neall,eall,nwall,wall)
c-------------------------------------------------------------------
      implicit none
      integer ibox
      integer isep
      integer nboxes,nlist3,list3(nlist3)
      double precision centers(3,nboxes)
      double precision sepdist,bs
      integer nuall,ndall,nnall,nsall,neall,nwall
      integer uall(1),dall(1),nall(1),sall(1),eall(1),wall(1)

      integer jbox
      integer c1,c2,c3,c4,c5,c6
      integer j

      nuall = 0
      ndall = 0
      nnall = 0
      nsall = 0
      neall = 0
      nwall = 0
      sepdist = 1.01d0*isep*bs+bs/2.0d0
      do j=1,nlist3
         jbox = list3(j)
C         print *,"jbox: ",jbox
         if(jbox.gt.0) then
            c1 = 0
            c2 = 0
            c3 = 0
            c4 = 0
            c5 = 0
            c6 = 0
            if((centers(3,jbox)-centers(3,ibox)).ge.sepdist) c1 = 1
            if((centers(3,jbox)-centers(3,ibox)).le.-sepdist) c2 = 1
            if((centers(2,jbox)-centers(2,ibox)).ge.sepdist) c3 = 1
            if((centers(2,jbox)-centers(2,ibox)).le.-sepdist) c4 = 1
            if((centers(1,jbox)-centers(1,ibox)).ge.sepdist) c5 = 1
            if((centers(1,jbox)-centers(1,ibox)).le.-sepdist) c6 = 1

C            print *,c1,c2,c3,c4,c5,c6

            if(c1.eq.1) then
               nuall = nuall + 1
               uall(nuall) = jbox
            endif

            if(c2.eq.1) then
               ndall = ndall + 1
               dall(ndall) = jbox
            endif

            if(c3.eq.1.and.c1.eq.0.and.c2.eq.0) then
               nnall = nnall + 1
               nall(nnall) = jbox
            endif

            if(c4.eq.1.and.c1.eq.0.and.c2.eq.0) then
               nsall = nsall + 1
               sall(nsall) = jbox
            endif

            if(c5.eq.1.and.c1.eq.0.and.c2.eq.0.and.c3.eq.0.and.
     1         c4.eq.0) then
               neall = neall + 1
               eall(neall) = jbox
            endif

            if(c6.eq.1.and.c1.eq.0.and.c2.eq.0.and.c3.eq.0.and.
     1         c4.eq.0) then
               nwall = nwall + 1
               wall(nwall) = jbox
            endif
         endif
      enddo

      return
      end
c
c
c
c
c
c
c--------------------------------------------------------------------

      subroutine getlist4pwdir(dir,censrc,centrg,boxsize)
c-------------------------------------------------------------------
      implicit none
ccc   input/output variables
      integer dir
      double precision censrc(3)
      double precision centrg(3)
      double precision boxsize
ccc   scopded function variables
      double precision sepdist
      double precision ctmp(3)
      ctmp(1) = censrc(1)-0*boxsize/2.0d0
      ctmp(2) = censrc(2)-0*boxsize/2.0d0
      ctmp(3) = censrc(3)-0*boxsize/2.0d0

      sepdist=1.51d0*boxsize

      if(ctmp(3)-centrg(3).ge.sepdist) then
        dir=1
      else if(ctmp(3)-centrg(3).le.-sepdist) then
        dir=2
      else if(ctmp(2)-centrg(2).ge.sepdist) then
        dir=3
      else if(ctmp(2)-centrg(2).le.-sepdist) then
        dir=4
      else if(ctmp(1)-centrg(1).ge.sepdist) then
        dir=5
      else if(ctmp(1)-centrg(1).le.-sepdist) then
        dir=6
      else
        dir=0
      end if
C      dir=0

      return
      end

      subroutine getlist4pwdirtest(dir,censrc,centrg,boxsize)
c-------------------------------------------------------------------
      implicit none
ccc   input/output variables
      integer dir
      double precision censrc(3)
      double precision centrg(3)
      double precision boxsize
ccc   scopded function variables
      double precision sepdist
      double precision ctmp(3)
      ctmp(1) = censrc(1)-0*boxsize/2.0d0
      ctmp(2) = censrc(2)-0*boxsize/2.0d0
      ctmp(3) = censrc(3)-0*boxsize/2.0d0

      sepdist=1.51d0*boxsize

      if((ctmp(3)-centrg(3)).ge.sepdist) then
        dir=1
      else if((ctmp(3)-centrg(3)).le.-sepdist) then
        dir=2
      else if((ctmp(2)-centrg(2)).ge.sepdist) then
        dir=3
      else if((ctmp(2)-centrg(2)).le.-sepdist) then
        dir=4
      else if((ctmp(1)-centrg(1)).ge.sepdist) then
        dir=5
      else if((ctmp(1)-centrg(1)).le.-sepdist) then
        dir=6
      else
        dir=0
        print *,"dir:",dir
      end if

      return
      end
c
c
c
c
c
c--------------------------------------------------------------
c
      subroutine subdividebox(pos,npts,center,boxsize,
     1           isorted,iboxfl,subcenters)
      implicit none
      double precision pos(3,npts)
      double precision center(3)
      double precision subcenters(3,8)
      double precision boxsize
      integer npts
      integer isorted(*)
      integer iboxfl(2,8)

c     Temporary variables
      integer isortedtmp(npts)
      integer i,j,i12,i34,istart,jstart,kstart,ii,iii,nss,nee
      integer jj,irefinebox,ntt
      integer i56, i78, i1234, i5678
      integer ibox,ifirstbox,ilastbox,nbfirst
      integer is,it,ie
      integer nc(8)

      i1234 = 0
      i5678 = 0
      do i = 1,npts
         if(pos(3,i)-center(3).lt.0) then
            i1234 = i1234 + 1
            isorted(i1234) = i
         else
            i5678 = i5678 + 1
            isortedtmp(i5678) = i
         endif
      enddo
      do i=1,i5678
         isorted(i1234+i) = isortedtmp(i)
      enddo

c     Sort i1234 into i12 and i34
      i12 = 0
      i34 = 0
      do i = 1,i1234
         if(pos(2,isorted(i))-center(2).lt.0) then
            i12 = i12 + 1
            isorted(i12) = isorted(i)
         else
            i34 = i34 + 1
            isortedtmp(i34) = isorted(i)
         endif
      enddo
c     Note at the end of the loop, i12 is where the particles
c     in part 12 of the box end
c
c     Reorder sources to include 34 in the array
      do i=1,i34
         isorted(i12+i) = isortedtmp(i)
      enddo

c     sort i5678 into i56 and i78
      i56 = i1234
      i78 = 0
      do i=i1234+1,npts
         if(pos(2,isorted(i))-center(2).lt.0) then
            i56 = i56 + 1
            isorted(i56) = isorted(i)
         else
            i78 = i78 + 1
            isortedtmp(i78) = isorted(i)
         endif
      enddo
c     Reorder sources to include 78 in the array
      do i=1,i78
         isorted(i56+i) = isortedtmp(i)
      enddo
c     End of reordering i5678

      nc = 0
c     Sort into boxes 1 and 2
      do i = 1,i12
         if(pos(1,isorted(i))-center(1).lt.0) then
            nc(1) = nc(1) + 1
            isorted(nc(1)) = isorted(i)
         else
            nc(2) = nc(2) + 1
            isortedtmp(nc(2)) = isorted(i)
         endif
      enddo
c     Reorder sources so that sources in 2 are at the
c     end of this part of the array
      do i=1,nc(2)
         isorted(nc(1)+i) = isortedtmp(i)
      enddo

c     Sort into boxes 3 and 4
      do i = i12+1, i1234
         if(pos(1,isorted(i))-center(1).lt.0) then
            nc(3) = nc(3) + 1
            isorted(i12+nc(3)) = isorted(i)
         else
            nc(4) = nc(4)+1
            isortedtmp(nc(4)) = isorted(i)
         endif
      enddo
c     Reorder sources so that sources in 4 are at the
c     end of this part of the array
      do i=1,nc(4)
         isorted(i12+nc(3)+i) = isortedtmp(i)
      enddo

c     Sort into boxes 5 and 6
      do i = i1234+1,i56
         if(pos(1,isorted(i))-center(1).lt.0) then
            nc(5) = nc(5) + 1
            isorted(i1234+nc(5)) = isorted(i)
         else
            nc(6) = nc(6) + 1
            isortedtmp(nc(6)) = isorted(i)
         endif
      enddo
c     Reorder sources so that sources in 6 are at the
c     end of this part of the array
      do i=1,nc(6)
         isorted(i1234+nc(5)+i) = isortedtmp(i)
      enddo
c     End of sorting sources into boxes 5 and 6

c     Sort into boxes 7 and 8
      do i=i56+1,npts
         if(pos(1,isorted(i))-center(1).lt.0) then
            nc(7) = nc(7) + 1
            isorted(i56+nc(7)) = isorted(i)
         else
            nc(8) = nc(8) + 1
            isortedtmp(nc(8)) = isorted(i)
         endif
      enddo
c     Reorder sources so that sources in 8 are at the
c     end of the array
      do i=1,nc(8)
         isorted(i56+nc(7)+i) = isortedtmp(i)
      enddo

      istart=1
      iboxfl=0
      subcenters=0.0d0
      do i=1,8
         ii = 2
         jj = 2
         if(i.eq.1.or.i.eq.2.or.i.eq.5.or.i.eq.6) ii = 1
         if(i.lt.5) jj = 1
         if(nc(i).gt.0) then
            subcenters(1,i) = center(1)+(-1)**i*boxsize/2.0d0
            subcenters(2,i) = center(2)+(-1)**ii*boxsize/2.0d0
            subcenters(3,i) = center(3)+(-1)**jj*boxsize/2.0d0

            iboxfl(1,i) = istart
            iboxfl(2,i) = istart+nc(i)-1

            istart = istart+nc(i)
          endif
      enddo

      return
      end
c--------------------------------------------------------------------
c
c
c


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions
C of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C
C     * Neither the name of NTESS nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

      SUBROUTINE PCKLAB (CURVE,ACURVE,NCURVE)
C
C     THIS SUBROUTINE PACKS A 8 CHARACTER WORD AND A 'I8' MAXIMUM
C     INTEGER INTO A 16 CHARACTER WORD BY REMOVING
C     INCLUDED BLANKS -- USED TO CREATE VARIABLE 'CURVE' FOR PLTONE
C
      CHARACTER*8 ACURVE
      CHARACTER*16 CURVE
      CHARACTER*1 BLANK
      DATA BLANK/' '/
C
      IF (NCURVE .NE. 0) THEN
         WRITE (CURVE, 10) ACURVE, NCURVE
   10    FORMAT (A8,I8)
       ELSE
         CURVE(1:8)= ACURVE(1:8)
         CURVE(9:16)= '        '
       END IF
      L = 0
      DO 20 J=1,16
         IF (CURVE(J:J).NE.BLANK) THEN
            L = L + 1
            CURVE(L:L)= CURVE(J:J)
          END IF
   20    CONTINUE
      L1 = L + 1
      DO 30 J=L1,16
         CURVE(J:J)= BLANK
   30    CONTINUE
      RETURN
C
      END

ext: f

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief This module provides functions to solve sparse system of
!> \linear equitions in the case of contact analysis using standard
!> \Lagrange multiplier algorithm

module m_solve_LINEQ_contact

  use m_fstr
  use m_solve_LINEQ_mkl
  use m_solve_LINEQ_direct_serial_lag
  use m_solve_LINEQ_MUMPS_contact
  use m_solve_LINEQ_iter_contact
  use m_fstr_mat_resid_contact

  implicit none

  private
  public :: solve_LINEQ_contact_init
  public :: solve_LINEQ_contact

contains

  !> \brief This subroutine
  subroutine solve_LINEQ_contact_init(hecMESH,hecMAT,fstrMAT,is_sym)
    type (hecmwST_local_mesh)                :: hecMESH        !< hecmw mesh
    type (hecmwST_matrix)                    :: hecMAT         !< type hecmwST_matrix
    type (fstrST_matrix_contact_lagrange)    :: fstrMAT        !< type fstrST_matrix_contact_lagrange)
    logical :: is_sym

    if( hecMAT%Iarray(99)==1 )then
      call solve_LINEQ_iter_contact_init(hecMESH,hecMAT,fstrMAT,is_sym)
    else if( hecMAT%Iarray(99)==3 )then
      call solve_LINEQ_mkl_init(hecMAT,fstrMAT,is_sym)
    elseif( hecMAT%Iarray(99)==4 )then
      call solve_LINEQ_serial_lag_hecmw_init(hecMAT,fstrMAT,is_sym)
    elseif( hecMAT%Iarray(99)==5 ) then
      call solve_LINEQ_mumps_contact_init(hecMESH,hecMAT,fstrMAT,is_sym)
    endif
  end subroutine solve_LINEQ_contact_init

  !> \brief This subroutine
  subroutine solve_LINEQ_contact(hecMESH,hecMAT,fstrMAT,istat,rf,conMAT)

    type (hecmwST_local_mesh)                :: hecMESH        !< hecmw mesh
    type (hecmwST_matrix)                    :: hecMAT         !< type hecmwST_matrix
    type (fstrST_matrix_contact_lagrange)    :: fstrMAT        !< type fstrST_matrix_contact_lagrange)
    integer(kind=kint), intent(out)          :: istat
    real(kind=kreal), optional               :: rf
    type (hecmwST_matrix),optional           :: conMAT

    real(kind=kreal)                         :: factor
    real(kind=kreal) :: resid
    real(kind=kreal) :: t1, t2

    factor = 1.0d0
    if( present(rf) )factor = rf

    t1 = hecmw_wtime()

    istat = 0
    if( hecMAT%Iarray(99)==1 )then
      if(paraContactFlag.and.present(conMAT)) then
        call solve_LINEQ_iter_contact(hecMESH,hecMAT,fstrMAT,conMAT)
      else
        call solve_LINEQ_iter_contact(hecMESH,hecMAT,fstrMAT)
      endif
    elseif( hecMAT%Iarray(99)==3 )then
      call solve_LINEQ_mkl(hecMESH,hecMAT,fstrMAT,istat)
    elseif( hecMAT%Iarray(99)==4 )then
      call solve_LINEQ_serial_lag_hecmw(hecMESH,hecMAT,fstrMAT)
    elseif( hecMAT%Iarray(99)==5 ) then
      ! ----  For Parallel Contact with Multi-Partition Domains
      if(paraContactFlag.and.present(conMAT)) then
        call solve_LINEQ_mumps_contact(hecMESH,hecMAT,fstrMAT,istat,conMAT)
      else
        call solve_LINEQ_mumps_contact(hecMESH,hecMAT,fstrMAT,istat)
      endif
    endif

    t2 = hecmw_wtime()
    if (hecmw_mat_get_timelog(hecMAT) .ge. 1) then
      if (myrank==0) write(*,*) ' solve time :', t2 - t1
    endif

    if(paraContactFlag.and.present(conMAT)) then
    else
      resid=fstr_get_resid_max_contact(hecMESH,hecMAT,fstrMAT)
      if (myrank==0) then
        write(*,*) ' maximum residual = ', resid
        if( hecmw_mat_get_solver_type(hecMAT) /= 1 .and. resid >= 1.0d-8) then
          write(*,*) ' ###Maximum residual exceeded 1.0d-8---Direct Solver### '
          !          stop
        endif
      endif
    endif

    hecMAT%X=factor*hecMAT%X

  end subroutine solve_LINEQ_contact

end module m_solve_LINEQ_contact


text:
module mod_monolis_sparse_util
  use mod_monolis_util
  use mod_monolis_stdlib
  use mod_monolis_graph
  implicit none

contains

  subroutine monolis_get_nonzero_pattern(monolis, nnode, nbase_func, ndof, nelem, elem)
    use iso_c_binding
    implicit none
    type(monolis_structure) :: monolis
    integer(kint) :: nnode, nbase_func, ndof, nelem, elem(:,:)
    integer(kint), pointer :: ebase_func(:), connectivity(:)
    integer(c_int), pointer :: index(:), item(:)

    call monolis_convert_mesh_to_connectivity &
      & (nelem, nbase_func, elem, ebase_func, connectivity)

    call monolis_convert_connectivity_to_nodal &
      & (nnode, nelem, ebase_func, connectivity, index, item)

     call monolis_get_nonzero_pattern_by_nodal &
       & (monolis, nnode, ndof, index, item)
  end subroutine monolis_get_nonzero_pattern

  subroutine monolis_get_nonzero_pattern_with_arbitrary_dof &
    (monolis, nnode, nbase_func, n_dof_list, nelem, elem)
    use iso_c_binding
    implicit none
    type(monolis_structure) :: monolis
    integer(kint) :: nnode, nbase_func, n_dof_list(:), nelem, elem(:,:)
    integer(kint), pointer :: ebase_func(:), connectivity(:)
    integer(c_int), pointer :: index(:), item(:)

    call monolis_convert_mesh_to_connectivity &
      & (nelem, nbase_func, elem, ebase_func, connectivity)

    call monolis_convert_connectivity_to_nodal &
      & (nnode, nelem, ebase_func, connectivity, index, item)

     call monolis_get_nonzero_pattern_by_nodal_graph_with_arbitrary_dof &
       & (monolis, nnode, n_dof_list, index, item)
  end subroutine monolis_get_nonzero_pattern_with_arbitrary_dof

  subroutine monolis_get_nonzero_pattern_by_connectivity &
      & (monolis, nnode, ndof, nelem, ebase_func, connectivity)
    use iso_c_binding
    implicit none
    type(monolis_structure) :: monolis
    integer(kint) :: nnode, ndof, nelem
    integer(kint), pointer :: ebase_func(:), connectivity(:)
    integer(c_int), pointer :: index(:), item(:)

    call monolis_convert_connectivity_to_nodal &
      & (nnode, nelem, ebase_func, connectivity, index, item)

     call monolis_get_nonzero_pattern_by_nodal &
      & (monolis, nnode, ndof, index, item)
  end subroutine monolis_get_nonzero_pattern_by_connectivity

  subroutine monolis_get_nonzero_pattern_by_nodal(monolis, nnode, ndof, index, item)
    use iso_c_binding
    implicit none
    type(monolis_structure) :: monolis
    integer(kint) :: nnode, ndof
    integer(kint) :: i, j, nz, jS, jE
    integer(c_int), pointer :: index(:), item(:)

    monolis%MAT%N = nnode
    monolis%MAT%NP = nnode
    monolis%MAT%NDOF = ndof
    allocate(monolis%MAT%X(ndof*nnode), source = 0.0d0)
    allocate(monolis%MAT%B(ndof*nnode), source = 0.0d0)
    allocate(monolis%MAT%index(0:nnode), source = 0)
    do i = 1, nnode
      monolis%MAT%index(i) = index(i+1) + i
    enddo

    nz = monolis%MAT%index(nnode)
    monolis%MAT%NZ = nz
    allocate(monolis%MAT%A(ndof*ndof*nz), source = 0.0d0)
    allocate(monolis%MAT%item(nz), source = 0)
    do i = 1, nnode
      jS = monolis%MAT%index(i-1) + 1
      jE = monolis%MAT%index(i)
      monolis%MAT%item(jS) = i
      do j = jS+1, jE
        monolis%MAT%item(j) = item(j-i)
      enddo
      call monolis_qsort_int(monolis%MAT%item(jS:jE), 1, jE - jS + 1)
    enddo

    allocate(monolis%MAT%indexR(0:nnode), source = 0)
    allocate(monolis%MAT%itemR(nz), source = 0)
    allocate(monolis%MAT%permR(nz), source = 0)

    call monolis_get_CRR_format(monolis%MAT%N, nz, monolis%MAT%index, monolis%MAT%item, &
      & monolis%MAT%indexR, monolis%MAT%itemR, monolis%MAT%permR)

    nullify(index)
    nullify(item)
  end subroutine monolis_get_nonzero_pattern_by_nodal

  subroutine monolis_get_nonzero_pattern_by_nodal_graph_with_arbitrary_dof &
    (monolis, nnode, n_dof_list, index, item)
    use iso_c_binding
    implicit none
    type(monolis_structure) :: monolis
    integer(kint) :: nnode, n_dof_list(:)
    integer(kint) :: i, j, k, nz, jS, jE, kS, kE
    integer(kint) :: total_dof, in, jn, kn, nrow, ncol, l
    integer(kint), allocatable :: n_dof_index(:)
    integer(c_int), pointer :: index(:), item(:)

    total_dof = 0
    do i = 1, nnode
      total_dof = total_dof + n_dof_list(i)
    enddo

    allocate(n_dof_index(nnode), source = 0)
    call monolis_get_n_dof_index(nnode, n_dof_list, n_dof_index)

    monolis%MAT%N = total_dof
    monolis%MAT%NP = total_dof
    monolis%MAT%NDOF = 1
    allocate(monolis%MAT%X(total_dof), source = 0.0d0)
    allocate(monolis%MAT%B(total_dof), source = 0.0d0)
    allocate(monolis%MAT%index(0:total_dof), source = 0)

    !> count nz
    nz = 0
    do i = 1, nnode
      jS = index(i) + 1
      jE = index(i+1)
      nz = nz + n_dof_list(i)*n_dof_list(i)
      do j = jS, jE
        jn = item(j)
        nz = nz +  n_dof_list(i)*n_dof_list(jn)
      enddo
    enddo

    monolis%MAT%NZ = nz
    allocate(monolis%MAT%A(nz), source = 0.0d0)
    allocate(monolis%MAT%item(nz), source = 0)

    !> construct index and item
    in = 0
    ncol = 0
    do i = 1, nnode
      jS = index(i) + 1
      jE = index(i+1)
      do k = 1, n_dof_list(i)
        ncol = ncol + 1
        nrow = n_dof_list(i)
        do j = 1, n_dof_list(i)
          in = in + 1
          monolis%MAT%item(in) = n_dof_index(i) + j
        enddo
        do j = jS, jE
          jn = item(j)
          nrow = nrow + n_dof_list(jn)
          do l = 1, n_dof_list(jn)
            in = in + 1
            monolis%MAT%item(in) = n_dof_index(jn) + l
          enddo
        enddo
        monolis%MAT%index(ncol) = monolis%MAT%index(ncol-1) + nrow

        kS = monolis%MAT%index(ncol-1) + 1
        kE = monolis%MAT%index(ncol)
        call monolis_qsort_int(monolis%MAT%item(kS:kE), 1, kE-kS+1)
      enddo
    enddo

    allocate(monolis%MAT%indexR(0:total_dof), source = 0)
    allocate(monolis%MAT%itemR(nz), source = 0)
    allocate(monolis%MAT%permR(nz), source = 0)

    call monolis_get_CRR_format(monolis%MAT%N, nz, monolis%MAT%index, monolis%MAT%item, &
      & monolis%MAT%indexR, monolis%MAT%itemR, monolis%MAT%permR)

    nullify(index)
    nullify(item)
  end subroutine monolis_get_nonzero_pattern_by_nodal_graph_with_arbitrary_dof

  subroutine monolis_get_n_dof_index(n_node, n_dof_list, n_dof_index)
    implicit none
    integer(kint), intent(in) :: n_node, n_dof_list(:)
    integer(kint) :: i, n_dof_index(:)

    do i = 1, n_node - 1
      n_dof_index(i+1) = n_dof_index(i) + n_dof_list(i)
    enddo
  end subroutine monolis_get_n_dof_index

  !> setter
  subroutine monolis_sparse_matrix_set_value(index, item, A, ndof, ci, cj, csub_i, csub_j, val)
    implicit none
    integer(kint), intent(in) :: ndof
    integer(kint), intent(in) :: index(0:), item(:), ci, cj, csub_i, csub_j
    real(kdouble), intent(inout) :: A(:)
    real(kdouble), intent(in) :: val
    integer(kint) :: j, jn, im, jS, jE, NDOF2

    NDOF2 = ndof*ndof
    if(ndof < csub_i) call monolis_stop_by_submatrix_access(ndof, csub_i)
    if(ndof < csub_j) call monolis_stop_by_submatrix_access(ndof, csub_j)

    jS = index(ci-1) + 1
    jE = index(ci)
    do j = jS, jE
      jn = item(j)
      if(jn == cj)then
        im = NDOF2*(j-1) + ndof*(csub_i-1) + csub_j
        A(im) = val
        return
      endif
      call monolis_stop_by_matrix_assemble(ci, cj)
    enddo
  end subroutine monolis_sparse_matrix_set_value

  !> getter
  subroutine monolis_sparse_matrix_get_value(index, item, A, ndof, ci, cj, csub_i, csub_j, val)
    implicit none
    integer(kint), intent(in) :: ndof
    integer(kint), intent(in) :: index(0:), item(:), ci, cj, csub_i, csub_j
    real(kdouble), intent(inout) :: A(:)
    real(kdouble), intent(out) :: val
    integer(kint) :: j, jn, im, jS, jE, NDOF2

    NDOF2 = ndof*ndof
    if(ndof < csub_i) call monolis_stop_by_submatrix_access(ndof, csub_i)
    if(ndof < csub_j) call monolis_stop_by_submatrix_access(ndof, csub_j)

    jS = index(ci-1) + 1
    jE = index(ci)
    do j = jS, jE
      jn = item(j)
      if(jn == cj)then
        im = NDOF2*(j-1) + ndof*(csub_i-1) + csub_j
        val = A(im)
        return
      endif
      call monolis_stop_by_matrix_assemble(ci, cj)
    enddo
  end subroutine monolis_sparse_matrix_get_value

  !> adder
  subroutine monolis_add_scalar_to_sparse_matrix(monolis, i, j, sub_i, sub_j, val)
    implicit none
    type(monolis_structure) :: monolis
    integer(kint), intent(in) :: i, j, sub_i, sub_j
    real(kdouble), intent(in) :: val

    call monolis_sparse_matrix_add_value(monolis%MAT%index, monolis%MAT%item, monolis%MAT%A, &
      & monolis%MAT%ndof, i, j, sub_i, sub_j, val)
  end subroutine monolis_add_scalar_to_sparse_matrix

  subroutine monolis_add_matrix_to_sparse_matrix(monolis, nbase_func, connectivity, stiff)
    implicit none
    type(monolis_structure) :: monolis
    integer(kint), intent(in) :: nbase_func, connectivity(nbase_func)
    real(kdouble), intent(in) :: stiff(:,:)

    call monolis_sparse_matrix_add_matrix(monolis%MAT%index, monolis%MAT%item, monolis%MAT%A, &
      & nbase_func, nbase_func, monolis%MAT%ndof, connectivity, connectivity, stiff)
  end subroutine monolis_add_matrix_to_sparse_matrix

  subroutine monolis_add_matrix_to_sparse_matrix_offdiag(monolis, n1, n2, c1, c2, stiff)
    implicit none
    type(monolis_structure) :: monolis
    integer(kint), intent(in) :: n1, n2, c1(n1), c2(n2)
    real(kdouble), intent(in) :: stiff(:,:)

    call monolis_sparse_matrix_add_matrix(monolis%MAT%index, monolis%MAT%item, monolis%MAT%A, &
      & n1, n2, monolis%MAT%ndof, c1, c2, stiff)
  end subroutine monolis_add_matrix_to_sparse_matrix_offdiag

  subroutine monolis_sparse_matrix_add_matrix(index, item, A, n1, n2, ndof, e1t, e2t, stiff)
    implicit none
    integer(kint), intent(in) :: n1, n2, ndof
    integer(kint), intent(in) :: index(0:), item(:), e1t(n1), e2t(n2)
    real(kdouble), intent(inout) :: A(:)
    real(kdouble), intent(in) :: stiff(n1*ndof,n2*ndof)
    integer(kint) :: e1(n1), e2(n2)
    integer(kint) :: i, j, k, in, jn, im, jS, jE, i2, j2, i1, j1, NDOF2
    integer(kint) :: eperm1(n1), eperm2(n2)
    real(kdouble) :: temp(n1*ndof,n2*ndof)

    NDOF2 = ndof*ndof

    e1 = e1t
    do i = 1, n1
      eperm1(i) = i
    enddo
    call monolis_qsort_int_with_perm(e1, 1, n1, eperm1)

    e2 = e2t
    do i = 1, n2
      eperm2(i) = i
    enddo
    call monolis_qsort_int_with_perm(e2, 1, n2, eperm2)

    temp = 0.0d0
    do i = 1, n2
      i1 = eperm2(i)
      do j = 1, n1
        j1 = eperm1(j)
        do i2 = 1, ndof
          do j2 = 1, ndof
            temp(ndof*(i-1)+i2, ndof*(j-1)+j2) = stiff(ndof*(j1-1)+j2, ndof*(i1-1)+i2)
          enddo
        enddo
      enddo
    enddo

    do i = 1, n1
      in = e1(i)
      jS = index(in-1) + 1
      jE = index(in)
      aa:do j = 1, n2
        do k = jS, jE
          jn = item(k)
          if(jn == e2(j))then
            do i1 = 1, ndof
            do i2 = 1, ndof
              im = NDOF2*(k-1) + ndof*(i1-1) + i2
              A(im) = A(im) + temp(ndof*(j-1)+i2, ndof*(i-1)+i1)
            enddo
            enddo
            jS = k + 1
            cycle aa
          endif
        enddo
      call monolis_stop_by_matrix_assemble(e1(i), e2(j))
      enddo aa
    enddo
  end subroutine monolis_sparse_matrix_add_matrix

  subroutine monolis_sparse_matrix_add_value(index, item, A, ndof, ci, cj, csub_i, csub_j, val)
    implicit none
    integer(kint), intent(in) :: ndof
    integer(kint), intent(in) :: index(0:), item(:), ci, cj, csub_i, csub_j
    real(kdouble), intent(inout) :: A(:)
    real(kdouble), intent(in) :: val
    integer(kint) :: j, jn, im, jS, jE, NDOF2
    character :: cerr*128

    NDOF2 = ndof*ndof
    if(ndof < csub_i) call monolis_stop_by_submatrix_access(ndof, csub_i)
    if(ndof < csub_j) call monolis_stop_by_submatrix_access(ndof, csub_j)

    jS = index(ci-1) + 1
    jE = index(ci)
    do j = jS, jE
      jn = item(j)
      if(jn == cj)then
        im = NDOF2*(j-1) + ndof*(csub_i-1) + csub_j
        A(im) = A(im) + val
        return
      endif
    enddo

    call monolis_stop_by_matrix_assemble(ci, cj)
  end subroutine monolis_sparse_matrix_add_value

  subroutine monolis_set_Dirichlet_bc(monolis, B, node_id, ndof_bc, val)
    implicit none
    type(monolis_structure) :: monolis
    integer(kint), intent(in) :: node_id, ndof_bc
    real(kdouble), intent(in) :: val
    real(kdouble) :: B(:)

    call monolis_sparse_matrix_add_bc(monolis%MAT%index, monolis%MAT%item, monolis%MAT%A, B, &
      & monolis%MAT%indexR, monolis%MAT%itemR, monolis%MAT%permR, &
      & monolis%MAT%ndof, node_id, ndof_bc, val)
  end subroutine monolis_set_Dirichlet_bc

  subroutine monolis_sparse_matrix_add_bc(index, item, A, B, indexR, itemR, permA, &
    & ndof, nnode, idof, val)
    implicit none
    integer(kint), intent(in) :: nnode, ndof, idof
    integer(kint), intent(in) :: index(0:), item(:), indexR(0:), itemR(:), permA(:)
    real(kdouble), intent(inout) :: A(:), B(:)
    real(kdouble), intent(in) :: val
    integer(kint) :: j, k, jn, kn, jS, jE, NDOF2
    logical :: is_add

    if(ndof < idof) call monolis_stop_by_submatrix_access(ndof, idof)

    is_add = .false.
    NDOF2 = ndof*ndof

    jS = indexR(nnode-1) + 1
    jE = indexR(nnode)
    do j = jS, jE
      jn = itemR(j)
      kn = permA(j)
      do k = 1, ndof
        B(ndof*(jn-1)+k) = B(ndof*(jn-1)+k) - val*A(NDOF2*(kn-1) + ndof*(k-1) + idof)
        A(NDOF2*(kn-1) + ndof*(k-1) + idof) = 0.0d0
      enddo
    enddo

    jS = index(nnode-1) + 1
    jE = index(nnode)
    do j = jS, jE
      do k = 1, ndof
        A(NDOF2*(j-1) + ndof*(idof-1) + k) = 0.0d0
      enddo

      jn = item(j)
      if(jn == nnode)then
        A(NDOF2*(j-1) + (ndof+1)*(idof-1) + 1) = 1.0d0
        is_add = .true.
      endif
    enddo

    if(.not. is_add) stop "error: not find a diagonal element in monolis_sparse_matrix_add_bc"

    B(ndof*nnode-ndof+idof) = val
  end subroutine monolis_sparse_matrix_add_bc

  subroutine monolis_get_CRR_format(N, NZ, index, item, indexR, itemR, permR)
    implicit none
    integer(kint), intent(in) :: N, NZ, index(0:), item(:)
    integer(kint), pointer :: indexR(:), itemR(:), temp(:), permR(:)
    integer(kint) :: i, j, in, jS, jE, m, p

    allocate(temp(N), source = 0)
    do i = 1, N
      jS = index(i-1) + 1
      jE = index(i)
      do j = jS, jE
        in = item(j)
        temp(in) = temp(in) + 1
      enddo
    enddo

    do i = 1, N
      indexR(i) = indexR(i-1) + temp(i)
    enddo

    temp = 0
    do i = 1, N
      jS = index(i-1) + 1
      jE = index(i)
      do j = jS, jE
        in = item(j)
        m = indexR(in-1)
        temp(in) = temp(in) + 1
        p = temp(in)
        itemR(m + p) = i
        permR(m + p) = j
      enddo
    enddo

    do i = 1, N
      jS = indexR(i-1) + 1
      jE = indexR(i)
      call monolis_qsort_int(itemR(jS:jE), 1, jE-jS+1)
    enddo
  end subroutine monolis_get_CRR_format

  subroutine monolis_stop_by_matrix_assemble(ci, cj)
    integer(kint), intent(in) :: ci, cj
    write(*,"(a,i0,a,i0,a)") "error: The non-zero element at (", ci, ", ", cj, &
      & ") is not allocated. The value is not accessible."
    stop
  end subroutine monolis_stop_by_matrix_assemble

  subroutine monolis_stop_by_submatrix_access(ndof, sub_dof)
    integer(kint), intent(in) :: ndof, sub_dof
    write(*,"(a)")   "error: set value greater than the DoF of submatrix."
    write(*,"(a,i8)")"       the DoF of submatrix: ", ndof
    write(*,"(a,i8)")"       value:                ", sub_dof
    stop
  end subroutine monolis_stop_by_submatrix_access

end module mod_monolis_sparse_util

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


text:
Winters Jr Warriors Football & Cheerleading  is a nonprofit organizations nonprofit organization 501(c)3.

Hi,

Winters Jr Warriors is a football and cheerleading organization in Winters, CA, Located in Northern CA. We need to raise money to purchase new equipment so we can replace our old equipment, but our organization has little to no money to purchase it on our own. We like your product and though it would be a good idea to raise funds through partnership through your company. I believe this could benefit both sides of the partnership, as you would receive advertising at a very low cost, while we receive a product a discounted price allowing us to fundraise even more money.

Sponsorship

If you are interested in sponsoring the Winters Jr Warriors and/or one of the football and cheerleading squads, please contact Jason Trojanowski at jmtrojanowski@me.com.
Sponsor the Winters Jr Warriors, and your business logo will be viewable on our team site www.myjrwarriors.com, and your sponsorship will be recognized at home games.
Sponsor an individual team, which your business is the only sponsor on their site, as well as being listed on our team site. Your sponsorship will be recognized at home games, and we may allow a reasonable sized business banner to be displayed on game day.
Sponsor through product donation. We accept money/food/equipment/etc., donations to help support our aging program. By sponsoring our program through product donation and/or discounted products, we will display you business logo on our team site, recognize your sponsorship at home games,  and allow a reasonable sized business banner to be placed on game day.

Winters Jr Warriors gladly accepts donations to support our football and cheerleading organization at:

Winters Jr Warriors Football & Cheerleading
PO Box 1273
Winters, CA 95694
Attn: Donations

Thank you,

Jason Trojanowski
Winters Jr Warriors, Board Member

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine scale ( arrin  , factor , nitem  , nvals  )

!       Deltares Software Centre

!>\file
!>                   Scales an array

!     Created            : March 1988 By M.E. Sileon / L. Postma
!                          May   2011    Leo Postma : Fortran 90 look and feel

!     Logical units      : none

!     Subroutines called : none

      use timers       !   performance timers

      implicit none

!     Parameters

!     kind           function         name                        Descriptipon

      integer  ( 4), intent(in   ) :: nvals                     !< number of values
      integer  ( 4), intent(in   ) :: nitem                     !< number of items
      real     ( 4), intent(inout) :: arrin (nvals,nitem)       !< number of items
      real     ( 4), intent(in   ) :: factor(nvals)             !< scale factors

!     local decalations

      integer        i1, i2       ! loop counters
      real     ( 4)  fact         ! factor
      integer(4) :: ithndl = 0
      if (timon) call timstrt( "scale", ithndl )

      do i1 = 1, nvals
         fact = factor(i1)
         do i2 = 1, nitem
            arrin (i1,i2) = arrin (i1,i2) * fact
         enddo
      enddo

      if (timon) call timstop( ithndl )
      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
Colors Photography is a photography service.  We do family portraits, high school seniors, children, pets, commercial images

   Davis has a wide variety of resources for Photography

20070705 23:55:05 nbsp does this place.. well have a place? I mean it has hours but is there a studio?  Users/DavidPoole


text:
*> \brief \b SGENND
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       LOGICAL FUNCTION SGENND (M, N, A, LDA)
*
*       .. Scalar Arguments ..
*       INTEGER M, N, LDA
*       ..
*       .. Array Arguments ..
*       REAL A( LDA, * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*>    SGENND tests that its argument has a non-negative diagonal.
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows in A.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns in A.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is REAL array, dimension (LDA, N)
*>          The matrix.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          Leading dimension of A.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup single_lin
*
*  =====================================================================
      LOGICAL FUNCTION SGENND (M, N, A, LDA)
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER M, N, LDA
*     ..
*     .. Array Arguments ..
      REAL A( LDA, * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO
      PARAMETER          ( ZERO = 0.0E0 )
*     ..
*     .. Local Scalars ..
      INTEGER I, K
*     ..
*     .. Intrinsics ..
      INTRINSIC MIN
*     ..
*     .. Executable Statements ..
      K = MIN( M, N )
      DO I = 1, K
         IF( A( I, I ).LT.ZERO ) THEN
            SGENND = .FALSE.
            RETURN
         END IF
      END DO
      SGENND = .TRUE.
      RETURN
      END

max_stars_repo_name: sydha/Lapack

max_issues_repo_name: sydha/Lapack

max_forks_repo_name: sydha/Lapack


!> Subroutine for Video 2 of VSCode for Fortran
!!
!! @author Lukas Lamm

subroutine getName(name)
    implicit none

    character(len=6), intent(out) :: name

    write(*,'(A)') "Please type your first name:"
    read(*,'(A)') name

end subroutine getName


text: Davis is home to both Bicycling bicyclists and people who enjoy Your Five Favorite Things In Davis nice things. An example of this intersection is http://www.brookssaddles.com/brooksengland.html Brooks Bicycle Saddles.


text:
    module MiscUtils
    implicit none

    INTERFACE PresentDefault
    module procedure PresentDefault_S, PresentDefault_L, PresentDefault_I, PresentDefault_R, PresentDefault_D
    END INTERFACE PresentDefault

    INTERFACE IfThenElse
    module procedure IfThenElse_S, IfThenElse_L, IfThenElse_I, IfThenElse_R, IfThenElse_D
    END INTERFACE IfThenElse

    INTERFACE IsFloat
    module procedure IsFloat0, IsFloat1, IsFloat2
    END INTERFACE IsFloat

    contains

    function DefaultFalse(S) result(Sout)
    logical, intent(in), optional :: S
    logical :: Sout

    if (present(S)) then
        SOut = S
    else
        SOut = .false.
    end if

    end function DefaultFalse

    function DefaultTrue(S) result(Sout)
    logical, intent(in), optional :: S
    logical :: Sout

    if (present(S)) then
        SOut = S
    else
        SOut = .true.
    end if

    end function DefaultTrue

    function PresentDefault_S(default, S) result(Sout)
    character(LEN=*), intent(in), target :: default
    character(LEN=*), intent(in), target, optional :: S
    character(LEN=:), pointer :: Sout

    if (present(S)) then
        SOut => S
    else
        SOut => default
    end if
    end function PresentDefault_S

    function PresentDefault_L(default, S) result(Sout)
    logical, intent(in) :: default
    logical, intent(in), optional :: S
    logical :: Sout

    if (present(S)) then
        SOut = S
    else
        SOut = default
    end if
    end function PresentDefault_L

    function PresentDefault_I(default, S) result(Sout)
    integer, intent(in) :: default
    integer, intent(in), optional :: S
    integer :: Sout

    if (present(S)) then
        SOut = S
    else
        SOut = default
    end if
    end function PresentDefault_I

    function PresentDefault_R(default, S) result(Sout)
    real, intent(in) :: default
    real, intent(in), optional :: S
    real:: Sout

    if (present(S)) then
        SOut = S
    else
        SOut = default
    end if
    end function PresentDefault_R

    function PresentDefault_D(default, S) result(Sout)
    double precision, intent(in) :: default
    double precision, intent(in), optional :: S
    double precision:: Sout

    if (present(S)) then
        SOut = S
    else
        SOut = default
    end if
    end function PresentDefault_D

    function IfThenElse_S(flag, either, or) result(IfThenElse)
    logical, intent(in) :: flag
    character(LEN=:), pointer :: IfThenElse
    character(LEN=*), target :: either, or

    if (flag) then
        IfThenElse => either
    else
        IfThenElse => or
    end if

    end function

    function IfThenElse_L(flag, either, or) result(IfThenElse)
    logical, intent(in) :: flag
    logical :: IfThenElse, either, or

    if (flag) then
        IfThenElse = either
    else
        IfThenElse = or
    end if

    end function

    function IfThenElse_I(flag, either, or) result(IfThenElse)
    logical, intent(in) :: flag
    integer :: IfThenElse, either, or

    if (flag) then
        IfThenElse = either
    else
        IfThenElse = or
    end if

    end function

    function IfThenElse_R(flag, either, or) result(IfThenElse)
    logical, intent(in) :: flag
    real :: IfThenElse, either, or

    if (flag) then
        IfThenElse = either
    else
        IfThenElse = or
    end if

    end function

    function IfThenElse_D(flag, either, or) result(IfThenElse)
    logical, intent(in) :: flag
    double precision :: IfThenElse, either, or

    if (flag) then
        IfThenElse = either
    else
        IfThenElse = or
    end if

    end function

    logical function isFloat0(R)
    class(*), intent(in) :: R

    select type(R)
    type is (real)
        isFloat0 = .true.
    type is (double precision)
        isFloat0 = .true.
        class default
        isFloat0 = .false.
    end select
    end function isFloat0

    logical function isFloat1(R)
    class(*), intent(in) :: R(:)

    isFloat1 = isFloat0(R(LBOUND(R,1)))

    end function isFloat1

    logical function isFloat2(R)
    class(*), intent(in) :: R(:,:)

    isFloat2 = isFloat0(R(LBOUND(R,1),LBOUND(R,2)))

    end function isFloat2

    end module MiscUtils

max_stars_repo_name: sjoudaki/CosmoJBD

max_issues_repo_name: sjoudaki/CosmoJBD

max_forks_repo_name: sjoudaki/CosmoJBD


text:
!  *******************************************************************
!  *                                                                 *
!  *                     Function dlafs3                             *
!  *                                                                 *
!  *******************************************************************
!  Double Precision Version 1.4
!  Written by Gordon A. Fenton, TUNS, Mar. 15, 1995
!  Latest Update: Jan 11, 2005
!
!  PURPOSE  returns the covariance between two points in a partially
!           separable 3-D fractional Gaussian noise field.
!
!  Returns the covariance between points in a 3-D fractional Gaussian noise
!  process which is isotropic in the X-Y plane. The covariance function is
!  partially separable, and is given by
!
!         B(X,Y,Z) = var*r(X,Y)*r(Z)
!
!  where r(Z) is the correlation function in the z-direction,
!
!               1               2G       2G           2G
!     r(Z) = -------- [ |Z + pb|  -  2|Z|  +  |Z - pb|  ]		(1)
!                 2G
!             2*pb
!
!  while, in the X-Y plane, the isotropic correlation function r(X,Y) is
!  given by
!
!                     1               2H       2H           2H
!  r(X,Y) = r(s) = -------- [ |s + pb|  -  2|s|  +  |s - pb|  ]		(2)
!                       2H
!                   2*pb
!
!  where s = sqrt(X*X + Y*Y), and `pb' is the length over which the fractional
!  Brownian motion is averaged in order to make this, the derivative process,
!  exist. Normally `pb' is selected to be quite small (of the order of the
!  size of the discretization interval). The parameters `G' and `H' are the
!  self-similar parameters, or Hurst exponents, (0.5 < G,H < 1) with
!  G = H = 1 giving a total correlated field and G = H = 0.5 giving a white
!  noise field (for arbitrarily small pb).
!
!  The parameters var, pb, H, and G are brought in through the common block
!  /dparam/. X, Y, and Z are the function arguments and are the elements
!  of the lag vector separating the two points in the field for which the
!  covariance is desired.
!
!  If var < 0, then this function computes the variance of a local average
!  of the random field, averaged over a domain having side dimensions {X,Y,Z}.
!  In this case, the variance function (which is also separable)
!
!     V(X,Y,Z) = V(X,Y)*V(Z)
!
!  is used, and the function returns |var|*V(X,Y,Z). The 1-D
!  variance function corresponding to Eq. (1) is given by
!
!                        r        r      r            r
!                |Z + pb| - 2( |Z| + |pb| ) + |Z - pb|
!       V(Z) = -----------------------------------------		(3)
!                      2                  2H
!                     Z * (2H+1)*(2H+2)*pb
!
!  where r = (2*G+2). Notice that as Z goes to zero the above becomes 0/0, but
!  its limiting value is 1.0.
!
!  The variance function which yields a close approximation to Eq. (2)c  is given by
!
!                      r          r        r              r
!              |D + pb|  -  2( |D|  +  |pb|  )  + |D - pb|
!   V(X,Y) = --------------------------------------------------		(4)
!                          2                  2G
!                         D * (2G+1)*(2G+2)*pb
!
!  where r = (2H+2) and D = |X| + |Y| (the 1-norm of {X,Y}). The actual
!  covariance function corresponding to (4) is equal to (2) along the
!  X,Y coordinate axes, and remains close to (2) for G > 0.8. For G < 0.7,
!  the covariance in the diagonal direction is too high, but since the
!  covariance in such fields dies off very rapidly, the error is
!  negligible.
!
!  Parameters to this process are brought in through the common block
!  /dparam/ and are described as follows;
!
!     var	the point variance of the fGn process.
!
!      pb	the averaging length discussed above.
!
!       H	the Hurst exponent governing the process in the X,Y plane.
!		In this implementation, 0.5 < H < 1; values of H near 1 yield
!		a correlation structure which remains very high (and thus a
!		covariance matrix which may be nearly singular). Values of H
!		near 0.5 yield a band- limited white noise process.
!
!       G	the Hurst exponent governing the process in the Z direction.
!		In this implementation, 0.5 < G < 1; values of G near 1 yield
!		a correlation structure which remains very high (and thus a
!		covariance matrix which may be nearly singular). Values of G
!		near 0.5 yield a band- limited white noise process.
!
!      da	dummy placeholder provided so that the common block /dparam/
!		has the same form for a variety of variance functions.
!
!  Arguments to this function are just the X, Y, and Z dimensions of the
!  lag distance (or the physical averaging region).
!
!  REVISION HISTORY:
!  1.1	added lvarfn flag to end of /dparam/ common block. Set it to
!	true if the variance function is to be returned. Otherwise this
!	function returns the covariance. (Apr 29/00)
!  1.11	minor changes to documentation above (Jul 14/00)
!  1.2	eliminated lvarfn -- now return covariances only if var < 0 (Apr 9/01)
!  1.3	reversed default - now return covariances if var > 0 (Apr 11/01)
!  1.31	revised above documentation to reflect revision 1.3 (May 9/01)
!  1.4	made X-Y the isotropic plane (Z is usually vertical) (Jan 11/05)
!----------------------------------------------------------------------
      real*8 function dlafs3( X, Y, Z )
      implicit real*8 (a-h,o-z)
      common/dparam/ var, pb, H, G, da
      data zero/0.d0/, one/1.d0/, two/2.d0/
      abs(q)  = dabs(q)
      sqrt(q) = dsqrt(q)

      if( var .lt. zero ) then			! return variance function
!						Z-direction variance function
         if( Z .eq. zero ) then
            v1 = one
         else
            e0 = two*G
            e1 = e0 + one
            e2 = e0 + two
            t1 = abs(Z+pb)**e2 - two*(abs(Z)**e2+pb**e2) + abs(Z-pb)**e2
            v1 = t1/(Z*Z*e1*e2*(pb**e0))
         endif
!						X-Y plane variance function
         if( (X .eq. zero) .and. (Y .eq. zero) ) then
            v2 = one
         else
            D  = abs(Y) + abs(X)
            e0 = two*H
            e1 = e0 + one
            e2 = e0 + two
            t1 = (D+pb)**e2 - two*(D**e2 + pb**e2) + abs(D-pb)**e2
            v2 = t1/(D*D*e1*e2*pb**e0)
         endif
         dlafs3 = -var*v1*v2
      else					! return covariance
!						Z-direction covariance
         if( Z .eq. zero ) then
            v1 = one
         else
            e0 = two*G
            t1 = abs(Z+pb)**e0 - two*(abs(Z)**e0) + abs(Z-pb)**e0
            v1 = t1/(two*(pb**e0))
         endif
!						Y-direction covariance
         if( (Y .eq. zero) .and. (X .eq. zero) ) then
            v2 = one
         else
            e0 = two*H
            D  = sqrt(Y*Y + X*X)
            t1 = (D+pb)**e0 - two*(D**e0) + abs(D-pb)**e0
            v2 = t1/(two*(pb**e0))
         endif
         dlafs3 = var*v1*v2
      endif

      return
      end

max_stars_repo_name: leemargetts/PoreFEM

max_issues_repo_name: leemargetts/PoreFEM

max_forks_repo_name: leemargetts/PoreFEM


text:
! { dg-do run }
! { dg-options "-ff2c" }
! Verifies that internal functions are not broken by f2c calling conventions
program test
  real, target :: f
  real, pointer :: q
  real :: g
  f = 1.0
  q=>f
  g = foo(q)
  if (g .ne. 1.0) call abort
contains
function foo (p)
  real, pointer :: foo
  real, pointer :: p
  foo => p
end function
end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
module line_mod
    use :: types_mod, only : sp
    implicit none

    private
    real(kind=sp) :: epsilon = 1.0e-4_sp
    public :: on_line, line, set_epsilon, get_epsilon

contains

    subroutine set_epsilon(eps)
        implicit none
        real(kind=sp), intent(in) :: eps
        epsilon = eps
    end subroutine set_epsilon

    function get_epsilon() result(eps)
        implicit none
        real(kind=sp) :: eps
        eps = epsilon
    end function get_epsilon

    elemental function line(x) result(y)
        implicit none
        real(kind=sp), intent(in) :: x
        real(kind=sp) :: y
        y = 3.0 + 0.5*x
    end function line

    elemental function on_line(x, y) result(is_on_line)
        implicit none
        real(kind=sp), intent(in) :: x, y
        logical :: is_on_line
        real(kind=sp) :: y_exp
        y_exp = line(x)
        is_on_line = equal(y, y_exp)

    contains

        pure function equal(a, b) result(is_equal)
            implicit none
            real(kind=sp), intent(in) :: a, b
            logical :: is_equal
            is_equal = abs(a - b) < epsilon
        end function equal

    end function on_line

end module line_mod

max_stars_repo_name: Gjacquenot/training-material

max_issues_repo_name: Gjacquenot/training-material

max_forks_repo_name: Gjacquenot/training-material


text: integer a (10)

max_stars_repo_name: mullikine/RosettaCodeData

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
!>
! @license Apache-2.0
!
! Copyright (c) 2019 The Stdlib Authors.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!    http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!<

!> Computes the dot product of two vectors.
!
! ## Notes
!
! * Modified verstion of reference BLAS level1 routine (version 3.7.0). Updated to "free form" Fortran 95.
!
! ## Authors
!
! * Univ. of Tennessee
! * Univ. of California Berkeley
! * Univ. of Colorado Denver
! * NAG Ltd.
!
! ## History
!
! * Jack Dongarra, linpack, 3/11/78.
!
!   - modified 12/3/93, array(1) declarations changed to array(*)
!
! ## License
!
! From <http://netlib.org/blas/faq.html>:
!
! > The reference BLAS is a freely-available software package. It is available from netlib via anonymous ftp and the World Wide Web. Thus, it can be included in commercial software packages (and has been). We only ask that proper credit be given to the authors.
! >
! > Like all software, it is copyrighted. It is not trademarked, but we do ask the following:
! >
! > * If you modify the source for these routines we ask that you change the name of the routine and comment the changes made to the original.
! >
! > * We will gladly answer any questions regarding the software. If a modification is done, however, it is the responsibility of the person who modified the routine to provide support.
!
! @param {integer} N - number of values over which to compute the dot product
! @param {Array<double>} sx - first array
! @param {integer} strideX - `sx` stride length
! @param {Array<double>} sy - second array
! @param {integer} strideY - `sy` stride length
! @returns {double} the dot product of `sx` and `sy`
!<
double precision function ddot( N, sx, strideX, sy, strideY )
  implicit none
  ! ..
  ! Scalar arguments:
  integer :: strideX, strideY, N
  ! ..
  ! Array arguments:
  double precision, intent(in) :: sx(*), sy(*)
  ! ..
  ! Local scalars:
  double precision :: stemp
  integer :: mp1, ix, iy, i, m
  ! ..
  ! Intrinsic functions:
  intrinsic mod
  ! ..
  stemp = 0.0e0
  ddot = 0.0e0
  ! ..
  if ( N <= 0 ) then
    return
  end if
  ! ..
  ! If both strides are equal to `1`, use unrolled loops...
  if ( strideX == 1 .AND. strideY == 1 ) then
    m = mod( N, 5 )
   ! ..
    ! If we have a remainder, do a clean-up loop...
    if ( m /= 0 ) then
      do i = 1, m
        stemp = stemp + ( sx( i ) + sy( i ) )
      end do
    end if
    if ( N < M ) then
      ddot = stemp
      return
    end if
    mp1 = m + 1
    do i = mp1, N, 5
      stemp = stemp + &
        ( sx( i ) * sy( i ) ) + &
        ( sx( i+1 ) * sy( i+1 ) ) + &
        ( sx( i+2 ) * sy( i+2 ) ) + &
        ( sx( i+3 ) * sy( i+3 ) ) + &
        ( sx( i+4 ) * sy( i+4 ) )
    end do
  else
    if ( strideX < 0 ) then
      ix = ( 1-N ) * strideX
    else
      ix = 0
    endif
    if ( strideY < 0 ) then
      iy = ( 1-N ) * strideY
    else
      iy = 0
    endif
    do i = 1, N
      stemp = stemp + ( sx( ix ) * sy( iy ) )
      ix = ix + strideX
      iy = iy + strideY
    end do
  endif
  ddot = stemp
  return
end function ddot

max_stars_repo_name: JamesGelok/stdlib

max_issues_repo_name: JamesGelok/stdlib

max_forks_repo_name: JamesGelok/stdlib