text:
!=======================================================================!
! Copyright (c) Intel Corporation - All rights reserved.                !
! This file is part of the LIBXSMM library.                             !
!                                                                       !
! For information on the license, see the LICENSE file.                 !
! Further information: https://github.com/hfp/libxsmm/                  !
! SPDX-License-Identifier: BSD-3-Clause                                 !
!=======================================================================!
! Hans Pabst (Intel Corp.)
!=======================================================================!
      PROGRAM dispatch_udt
        USE, INTRINSIC :: ISO_C_BINDING,  ONLY: C_PTR, C_LOC,           &
     &                                          C_ASSOCIATED,           &
     &                                          C_F_POINTER
        USE :: LIBXSMM, ONLY: LIBXSMM_BLASINT_KIND,                     &
     &                        LIBXSMM_MMFUNCTION => LIBXSMM_DMMFUNCTION,&
     &                        libxsmm_mmdispatch => libxsmm_dmmdispatch,&
     &                        libxsmm_mmcall => libxsmm_dmmcall,        &
     &                        libxsmm_xregister, libxsmm_xdispatch
        IMPLICIT NONE
        INTEGER, PARAMETER :: T = KIND(0D0)
        INTEGER :: batchsize = 1000, i
        INTEGER(LIBXSMM_BLASINT_KIND) :: j, ki, nrepeat = 100
        INTEGER(LIBXSMM_BLASINT_KIND) :: m = 13, n = 5, k = 7
        REAL(T), ALLOCATABLE :: a(:,:,:), b(:,:,:), c(:,:)
        TYPE(LIBXSMM_MMFUNCTION), TARGET  :: xmm(2) ! array of kernels
        TYPE(LIBXSMM_MMFUNCTION), POINTER :: udt(:)
        INTEGER(LIBXSMM_BLASINT_KIND), TARGET :: key(3)
        TYPE(C_PTR) :: ptr

        ALLOCATE(a(m,k,batchsize), b(k,n,batchsize), c(m,n))
        ! initialize input
        DO i = 1, batchsize
          DO ki = 1, k
            DO j = 1, m
              a(j,ki,i) = REAL(1, T) / REAL(MOD(i+j+ki, 25), T)
            END DO
            DO j = 1, n
              b(ki,j,i) = REAL(7, T) / REAL(MOD(i+j+ki, 75), T)
            END DO
          END DO
        END DO
        c(:,:) = 0

        ! repeat inner part to exercise libxsmm_xdispatch
        DO j = 1, nrepeat
          key = (/m, n, k/) ! setup key
          ! query associated value using key
          ptr = libxsmm_xdispatch(                                      &
     &      C_LOC(key), SIZE(key) * LIBXSMM_BLASINT_KIND)

          IF (C_ASSOCIATED(ptr)) THEN ! value was already registered
            ! convert C-ptr to Fortran POINTER
            CALL C_F_POINTER(ptr, udt, (/SIZE(xmm)/))
          ELSE ! no value registered yet
            ! generate and dispatch a series of kernels
            CALL libxsmm_mmdispatch(xmm(1), m, n, k,                    &
     &        alpha=REAL(1, T), beta=REAL(1, T))
            CALL libxsmm_mmdispatch(xmm(2), m, n, k + 2,                &
     &        alpha=REAL(1, T), beta=REAL(1, T))
            ! register an entry that contains all kernels from above
            ptr = libxsmm_xregister(                                    &
     &        C_LOC(key), SIZE(key) * LIBXSMM_BLASINT_KIND,             &
     &        SIZE(xmm) * 8, C_LOC(xmm))
            ! point udt to xmm (below code uses udt to refer to kernels
            udt => xmm ! alternatively, use C_F_POINTER
          END IF

          ! here we executed libxsmm_xdispatch one time (for this round)
          ! all kernels have been dispatched at once (udt)
          DO i = 1, batchsize
            CALL libxsmm_mmcall(udt(1), a(:,:,i), b(:,:,i), c)
          END DO
        END DO
        DEALLOCATE(a, b, c)
      END PROGRAM

ext: f

max_stars_repo_name: abhisek-kundu/libxsmm

max_issues_repo_name: abhisek-kundu/libxsmm

max_forks_repo_name: abhisek-kundu/libxsmm


text:
      SUBROUTINE IDPLOT (IDX)
C
      COMMON /OUTPUT/ SKPOUT(32,6),ID(32)
      COMMON /PLTDAT/ SKPPLT(2),XYMIN(2),XYMAX(2),AXYMAX(2),EDGE(12)
     1,      SKPA(3),CNTX,CNTY,SKPB(4),PLTYPE
      INTEGER PLTYPE
C
      INTEGER BLANK
      REAL SAVE(2,4)
      DATA BLANK,LINSIZ / 1H ,3 /
C
C     DOES A PLOT ID EXIST AT ALL
C
      IDX = 1
      DO 101 I = 1,20
      IF (ID(I).NE.BLANK)  GO TO 102
  101 CONTINUE
      IDX = 0
      GO TO 200
C
  102 DO 103 I = 1,2
      SAVE(I,1) = XYMIN(I)
      XYMIN(I) = 0.
      SAVE(I,2) = XYMAX(I)
      XYMAX(I) = AXYMAX(I)+EDGE(I)
      SAVE(I,3) = AXYMAX(I)
      AXYMAX(I) = XYMAX(I)
      SAVE(I,4) = EDGE(I)
      EDGE(I) = 0.
  103 CONTINUE
      NLINES = (AXYMAX(2)-7.*CNTY) / FLOAT(2*LINSIZ) + .1
      IF (IABS(PLTYPE).NE.1)  GO TO 122
C
C     FILL TOP HALF OF PLOT WITH X-AXIS LINES ALL THE WAY ACROSS.
C
      CALL AXIS (0,0,0,0,0,-1)
      DO 111 I = 1,NLINES
      Y = XYMAX(2) - FLOAT((I-1)*LINSIZ)
      CALL AXIS (XYMIN(1),Y,XYMAX(1),Y,1,0)
  111 CONTINUE
C
C     PRINT THE PLOT ID 2 TIMES IN THE MIDDLE OF THE PLOT.
C
      CALL PRINT (0,0,0,0,0,-1)
      X = XYMIN(1) + AMAX1(0.,(AXYMAX(1)-80.*CNTX)/2.)
      YY = Y-CNTY
      DO 116 I = 1,2
      Y = YY - CNTY*FLOAT(I-1)
      CALL PRINT (X,Y,1,ID,20,0)
  116 CONTINUE
C
C     FILL BOTTOM HALF OF PLOT WITH X-AXIS LINES ALL THE WAY ACROSS.
C
      CALL AXIS (0,0,0,0,0,-1)
      DO 121 I = 1,NLINES
      Y = XYMIN(2) + FLOAT((I-1)*LINSIZ)
      CALL AXIS (XYMIN(1),Y,XYMAX(1),Y,1,0)
  121 CONTINUE
      CALL AXIS (0,0,0,0,0,1)
      GO TO 125
C
C     NOT A CRT PLOTTER. TYPE THE ID ONCE AT THE BOTTOM OF THE PAPER.
C
  122 CALL PRINT (0,0,0,0,0,-1)
      X = XYMIN(1) + AMAX1(0.,(AXYMAX(1)-80.*CNTX)/2.)
      Y = 0.
      IF (PLTYPE.LT.0)  Y=CNTY/2.
      CALL PRINT (X,Y,1,ID,20,0)
C
C     END OF ID PLOT. PUT BLANKS IN THE PLOT ID.
C
  125 CALL PRINT (0,0,0,0,0,1)
      DO 126 I = 1,20
      ID(I) = BLANK
  126 CONTINUE
      DO 127 I = 1,2
      XYMIN(I) = SAVE(I,1)
      XYMAX(I) = SAVE(I,2)
      AXYMAX(I) = SAVE(I,3)
      EDGE(I) = SAVE(I,4)
  127 CONTINUE
C
  200 RETURN
      END


module pic

    use iso_c_binding
    implicit none

contains

pure integer function mod1( x, y)

integer(c_int32_t), intent(in) :: x
integer(c_int32_t), intent(in) :: y

mod1 = modulo(x-1, y) + 1

end function mod1

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

pure subroutine push_v( nbpart, p, dt, nx, ny, dx, dy, ex, ey, bz) bind(C, name="push_v")

integer(c_int32_t), intent(in) :: nbpart
real(c_double), intent(inout) :: p(5,nbpart)
real(c_double), intent(in) :: dt
integer(c_int32_t), intent(in) :: nx, ny
real(c_double), intent(in) :: dx, dy
real(c_double), intent(in) :: ex(nx+1,ny+1)
real(c_double), intent(in) :: ey(nx+1,ny+1)
real(c_double), intent(in) :: bz(nx+1,ny+1)

real(c_double) :: tantheta, sintheta
real(c_double) :: hdt, v1, v2, e1, e2, b3
real(c_double) :: a1, a2, a3, a4
real(c_double) :: xp, yp, dxp, dyp
integer(c_int32_t) :: i, j, ipart

hdt = 0.5 * dt

do ipart = 1, nbpart

   v1 = p(3,ipart)
   v2 = p(4,ipart)

   xp = p(1,ipart) / dx
   yp = p(2,ipart) / dy

   i = floor( xp ) + 1
   j = floor( yp ) + 1

   dxp = xp - i + 1
   dyp = yp - j + 1

   a1 = (1-dxp) * (1-dyp)
   a2 = dxp * (1-dyp)
   a3 = dxp * dyp
   a4 = (1-dxp) * dyp

   e1 = a1 * ex(i,j) + a2 * ex(i+1,j) + a3 * ex(i+1,j+1) + a4 * ex(i,j+1)
   e2 = a1 * ey(i,j) + a2 * ey(i+1,j) + a3 * ey(i+1,j+1) + a4 * ey(i,j+1)
   b3 = a1 * bz(i,j) + a2 * bz(i+1,j) + a3 * bz(i+1,j+1) + a4 * bz(i,j+1)

   v1 = v1 + hdt * e1
   v2 = v2 + hdt * e2

   tantheta = hdt * b3
   sintheta = 2.0 * tantheta / ( 1. + tantheta*tantheta)

   v1 = v1 + v2 * tantheta
   v2 = v2 - v1 * sintheta
   v1 = v1 + v2 * tantheta

   p(3,ipart) = v1 + hdt * e1
   p(4,ipart) = v2 + hdt * e2

end do

end subroutine push_v

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

pure subroutine push_x( nbpart, dimx, dimy, p, dt ) bind(C, name="push_x")

integer(c_int32_t), intent(in) :: nbpart
real(c_double), intent(in) :: dimx, dimy, dt
real(c_double), intent(inout) :: p(5,nbpart)
real(c_double) ::  p1, p2
integer :: ipart

do ipart=1,nbpart
   p1 = p(1,ipart) + p(3,ipart) * dt
   p2 = p(2,ipart) + p(4,ipart) * dt
   if ( p1 > dimx ) p1 = p1 - dimx
   if ( p1 < 0.0  ) p1 = p1 + dimx
   if ( p2 > dimy ) p2 = p2 - dimy
   if ( p2 < 0.0  ) p2 = p2 + dimy
   p(1,ipart) = p1
   p(2,ipart) = p2
end do

end subroutine push_x

pure subroutine deposition( nbpart, nx, ny, dx, dy, p, jx, jy ) bind(C, name="deposition")

integer(c_int32_t), intent(in) :: nbpart, nx, ny
real(c_double), intent(in)  :: dx, dy
real(c_double), intent(in)  :: p(5,nbpart)
real(c_double), intent(out) :: jx(nx+1,ny+1)
real(c_double), intent(out) :: jy(nx+1,ny+1)

real(c_double) :: a1, a2, a3, a4, w1, w2, xp, yp, dxp, dyp, factor
integer(c_int32_t) :: ipart, i, j

jx = 0.0
jy = 0.0

factor = real(nx * ny) / nbpart

do ipart=1,nbpart

   xp = p(1,ipart) / dx
   yp = p(2,ipart) / dy

   i = floor(xp) + 1
   j = floor(yp) + 1

   dxp = xp - i + 1
   dyp = yp - j + 1

   a1 = (1-dxp) * (1-dyp)
   a2 = dxp * (1-dyp)
   a3 = dxp * dyp
   a4 = (1-dxp) * dyp

   w1 = p(3,ipart) * factor
   w2 = p(4,ipart) * factor

   jx(i,j)     = jx(i,j)     + a1*w1
   jx(i+1,j)   = jx(i+1,j)   + a2*w1
   jx(i+1,j+1) = jx(i+1,j+1) + a3*w1
   jx(i,j+1)   = jx(i,j+1)   + a4*w1

   jy(i,j)     = jy(i,j)     + a1*w2
   jy(i+1,j)   = jy(i+1,j)   + a2*w2
   jy(i+1,j+1) = jy(i+1,j+1) + a3*w2
   jy(i,j+1)   = jy(i,j+1)   + a4*w2

end do

end subroutine deposition

end module pic


SUBROUTINE invert(matrix)
!
! This subroutine inverts a small square matrix onto itself.
!
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 REAL(iwp),INTENT(IN OUT)::matrix(:,:)
 REAL(iwp)::det,j11,j12,j13,j21,j22,j23,j31,j32,j33,con
 INTEGER::ndim,i,k
 ndim=UBOUND(matrix,1)
 IF(ndim==2)THEN
   det=matrix(1,1)*matrix(2,2)-matrix(1,2)*matrix(2,1)
   j11=matrix(1,1)
   matrix(1,1)=matrix(2,2)
   matrix(2,2)=j11
   matrix(1,2)=-matrix(1,2)
   matrix(2,1)=-matrix(2,1)
   matrix=matrix/det
 ELSE IF(ndim==3)THEN
   det=matrix(1,1)*(matrix(2,2)*matrix(3,3)-matrix(3,2)*matrix(2,3))
   det=det-matrix(1,2)*(matrix(2,1)*matrix(3,3)-matrix(3,1)*matrix(2,3))
   det=det+matrix(1,3)*(matrix(2,1)*matrix(3,2)-matrix(3,1)*matrix(2,2))
   j11=matrix(2,2)*matrix(3,3)-matrix(3,2)*matrix(2,3)
   j21=-matrix(2,1)*matrix(3,3)+matrix(3,1)*matrix(2,3)
   j31=matrix(2,1)*matrix(3,2)-matrix(3,1)*matrix(2,2)
   j12=-matrix(1,2)*matrix(3,3)+matrix(3,2)*matrix(1,3)
   j22=matrix(1,1)*matrix(3,3)-matrix(3,1)*matrix(1,3)
   j32=-matrix(1,1)*matrix(3,2)+matrix(3,1)*matrix(1,2)
   j13=matrix(1,2)*matrix(2,3)-matrix(2,2)*matrix(1,3)
   j23=-matrix(1,1)*matrix(2,3)+matrix(2,1)*matrix(1,3)
   j33=matrix(1,1)*matrix(2,2)-matrix(2,1)*matrix(1,2)
   matrix(1,1)=j11
   matrix(1,2)=j12
   matrix(1,3)=j13
   matrix(2,1)=j21
   matrix(2,2)=j22
   matrix(2,3)=j23
   matrix(3,1)=j31
   matrix(3,2)=j32
   matrix(3,3)=j33
   matrix=matrix/det
 ELSE
   DO k=1,ndim
     con=matrix(k,k)
     matrix(k,k)=1.0_iwp
     matrix(k,:)=matrix(k,:)/con
     DO i=1,ndim
       IF(i/=k)THEN
         con=matrix(i,k)
         matrix(i,k)=0.0_iwp
         matrix(i,:)=matrix(i,:)-matrix(k,:)*con
       END IF
     END DO
   END DO
 END IF
RETURN
END SUBROUTINE invert


text:
Gunrock is the name of the official blue horses horsebased school mascot of UC Davis.  It frequently makes appearances at UC Davis sporting events, acting as the face of UCD Athletics.

Students have described it as an evil trollish looking thing. Others have called it a giddy cow in a horses costume. It stands in as mascot at UCD sports games, as we have yet to discover specifically what an Glossary aggie really is.

 Whered Gunrock Come From?
For a long time there was just the loose concept of the Mustang representing the school, and UC Davis was therefore without an actual costumed mascot for many decades until 2003 when a blue mustang came trotting onto the campus. He took on the name Gunrock after an online internet poll showed that the UCD community overwhelmingly supported this appellation.

But where did the name Gunrock come from?

Gunrock is a very important name in the history of UC Davis, according to UCD Sports Information and Aggie Pack.  In the 1920s, the efforts of choosing a school mascot led many Cal Davis University Farm students to consider their agricultural and veterinary school roots.  Eventually, they chose the glorious mustang... a strongwilled, enduring, intelligent beast exhibiting exactly 1 horsepower of strapping vigor.

But what of the name?

In the 1920s, the University Farm was also used as a breeding facility for supplying highquality stock thoroughbreds for the United States Army Remount Service. Around this time of mascot namecalling, a thoroughbred named Gun Rock was brought to the campus for the Army breeding program. The offspring of Triple Crown winner wiki:wikipedia:Rock_Sand Rock Sand and race mare Gunfire, it turns out that Gun Rock was a relative of the famed racehorse, Man o War, and coincidentally, was also a relative of another famous racehorse named Sea Biscuit, repopularized by the 2003 Summerhit movie of the same name.

In honor of this handsome creature, the students of the 20s chose Gunrock as the name for the official school mascot, which was reaffirmed by the UCD community in 2003.  In two separate votes of the ASUCD during the 1990s, students voted to change the mascot to a Cows cow.

Sometime in 19992000, one of the free rewards that you got for attending an Aggie Pack event was a plush blue horse stuffed animal. You could say that this was probably the first appearance of the mascot in some kind of physical form.

Pictures

 Visiting Team Mascots

When Gunrock trots out onto the field what sort of mascot will he be up against?

  Stanford Tree
  UC Irvine Anteater?

Any way we could find a list of Gunrocks real identities?  I believe the current one just got his PhD, and that made me wonder who it really was.  Users/AllisonEriksen
    I know the current Gunrock I think he would prefer it to be kept a secret.  Users/RobRoy
    A gunrock just recently (March 2006) passed on the gunrocktorch, so to speak.  We ought to figure out who all the past gunrocks are.  There was a Davis Enterprise article about the old gunrock during this time, if someone wants to try and find it in their archives.
    The old Gunrock is friends with Michael Leahy of KDVS.  I wont say more than that.  Users/JeremyNewton

20081227 18:18:23 nbsp United States Army Calvary?
Its CAVALRY.
Calvary is the hill where Jesus was crucified.
Cavalry are mouted soldiers. Users/rswhitaker
   Why are you complaining after you fixed it?  Somebody made an honest mistake and everybody else missed the error until you caught it.  No need to mock everybody involved.  Even your comment has a benign error: Cavalry are mounted soldiers, not mouted soldiers.  Nobody is perfect... the process of helping out each other with editing is how the wiki works. Users/JabberWokky Evan JabberWokky Edwards

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      program testh
c...............................................................|
c
c   program used to study performance:  test host
c
c   code:  erlebacher
c          - calculation of variable derivatives
c          - tri-diagonal solver
c          - sequential for HYPERCUBE target
c
c   Thomas M. Eidson
c   ICASE
c   March, 1991
c...............................................................|
c
      parameter (ioprt=6)
c
      common /grids/ dx,dy,dz
      common /times/ dt,dtx,dty,dtz
      double precision dt,dtx,dty,dtz
      include 'erle.h'

      call setgrid
      call dotest
C     call prntout

      call checksum(dux)
      call checksum(duy)
      call checksum(duz)

      stop 0
      end

c===============================================================|
      subroutine checksum(ar)
      include 'par.h'
      double precision ar(N,N,N)

      double precision sum
      sum = 0.0
      do i = 1, N
         do j = 1,N
            do k = 1,N
              sum = sum + ar(k, j, i)
            enddo
         enddo
      enddo
      print *, sum
      end

c===============================================================|
      subroutine tridvpk(a,b,c,d,e,tot)

      include 'erle.h'
      double precision tot(N,N)
      double precision a(N),b(N),c(N),d(N),e(N)

c
c...............................................................|
c
c     perform forward substitution
c  //@; BEGIN(looplist1=CODE.LoopNestList)
      do 10 j=1,N
      do 10 i=1,N
         duz(i,j,1) = duz(i,j,1)*b(1)
   10 continue

      do  k=2,N-1
      do  j=1,N
      do  i=1,N
         duz(i,j,k)=(duz(i,j,k)-a(k)*duz(i,j,k-1))*b(k)
      enddo
      enddo
      enddo

   20 continue
c
      do 30 j=1,N
      do 30 i=1,N
         tot(i,j) = 0.
   30 continue
c
      do 40 k=1,N-1
      do 40 j=1,N
      do 40 i=1,N
            tot(i,j) = tot(i,j) + d(k)*duz(i,j,k)
   40 continue
c
      do 50 j=1,N
      do 50 i=1,N
         duz(i,j,N) = (duz(i,j,N) - tot(i,j))*b(N)
   50 continue
c
c     perform backward substitution
c
      do 60 j=1,N
      do 60 i=1,N
         duz(i,j,N-1)=duz(i,j,N-1) - e(N-1)*duz(i,j,N)
   60 continue
c
      do  k=N-2,1,-1
      do  j=1,N
      do  i=1,N
         duz(i,j,k) = duz(i,j,k) - c(k)*duz(i,j,k+1) - e(k)*duz(i,j,N)
      enddo
      enddo
      enddo

   70 continue
c
      return
      end
c===============================================================|
      subroutine tridvpj(a,b,c,d,e,tot)

      include 'erle.h'
        double precision tot(N,N)
        double precision a(N),b(N),c(N),d(N),e(N)

c
c...............................................................|
c
c     perform forward substitution
c
c  //@; BEGIN(looplist2=CODE.LoopNestList)
      do 10 k=1,N
      do 10 i=1,N
         duy(i,1,k) = duy(i,1,k)*b(1)
   10 continue
c
      do 20 k=1,N
      do 20 j=2,N-1
      do 20 i=1,N
         duy(i,j,k)=(duy(i,j,k)-a(j)*duy(i,j-1,k))*b(j)
   20 continue
c
      do 30 k=1,N
      do 30 i=1,N
         tot(i,k) = 0.
   30 continue
c
      do 40 k=1,N
      do 40 j=1,N-1
      do 40 i=1,N
         tot(i,k) = tot(i,k) + d(j)*duy(i,j,k)
   40 continue
c
      do 50 k=1,N
      do 50 i=1,N
         duy(i,N,k) = (duy(i,N,k) - tot(i,k))*b(N)
   50 continue
c
c     perform backward substitution
c
      do 60 k=1,N
      do 60 i=1,N
         duy(i,N-1,k)=duy(i,N-1,k) - e(N-1)*duy(i,N,k)
   60 continue
c
      do 70 k=1,N
      do 70 j=N-2,1,-1
      do 70 i=1,N
         duy(i,j,k) = duy(i,j,k) - c(j)*duy(i,j+1,k) - e(j)*duy(i,N,k)
   70 continue
c
      return
      end
c===============================================================|
      subroutine tridvpi(a,b,c,d,e,tot)

      include 'erle.h'
        double precision tot(N,N)
        double precision a(N),b(N),c(N),d(N),e(N)

c
c...............................................................|
c
c     perform forward substitution
c
c  //@; BEGIN(looplist3=CODE.LoopNestList)
      do 10 k=1,N
      do 10 j=1,N
         dux(1,j,k) = dux(1,j,k)*b(1)
   10 continue
c
      do 20 k=1,N
      do 20 j=1,N
      do 20 i=2,N-1
         dux(i,j,k)=(dux(i,j,k)-a(i)*dux(i-1,j,k))*b(i)
   20 continue
c
      do 30 k=1,N
      do 30 j=1,N
         tot(j,k) = 0.
   30 continue
c
      do 40 k=1,N
      do 40 j=1,N
      do 40 i=1,N-1
         tot(j,k) = tot(j,k) + d(i)*dux(i,j,k)
   40 continue
c
      do 50 k=1,N
      do 50 j=1,N
         dux(N,j,k) = (dux(N,j,k) - tot(j,k))*b(N)
   50 continue
c
c     perform backward substitution
c
      do 60 k=1,N
      do 60 j=1,N
         dux(N-1,j,k)=dux(N-1,j,k) - e(N-1)*dux(N,j,k)
   60 continue
c
      do 70 k=1,N
      do 70 j=1,N
      do 70 i=N-2,1,-1
         dux(i,j,k) = dux(i,j,k) - c(i)*dux(i+1,j,k) - e(i)*dux(N,j,k)
   70 continue
c
      return
      end
c===============================================================|
      subroutine tridvpf(a,b,c,d,e)
c...............................................................|
c
c   routine:
c      factorization stage of vectorized tridiagonal solver
c      equations are periodic
c      all equations have the same coefficients
c
c...............................................................|
c
      include 'erle.h'
      double precision a(N),b(N),c(N),d(N),e(N)
c
c...............................................................|
c
c     generate first elements of LU
c
      b(1) = 1.0/b(1)
      c(1) = c(1)*b(1)
      e(1) = a(1)*b(1)
      d(1) = c(N)
c
c     generate n=2 to n=N-2 elements of LU
c
      do 10 k=2,N-2
         b(k) = b(k) - a(k)*c(k-1)
         b(k) = 1.0/b(k)
         c(k) = c(k)*b(k)
         e(k) = -a(k)*e(k-1)*b(k)
         d(k) = -d(k-1)*c(k-1)
   10 continue
c
c     generate N-1 elements
c
      b(N-1) = b(N-1) - a(N-1)*c(N-2)
      b(N-1) = 1.0/b(N-1)
      e(N-1) = (c(N-1) - a(N-1)*e(N-2))*b(N-1)
      d(N-1) = a(N) - d(N-2)*c(N-2)
c
c     generate the n-th element
c
      tot = 0.
      do 20 k=1,N-1
         tot = tot + d(k)*e(k)
   20 continue
      b(N) = b(N) - tot
      b(N) = 1.0/b(N)
c
      return
      end

      subroutine prnt(ar)
      include 'erle.h'
      double precision ar(N,N,N)

c
c
      parameter (size = N)
      do i = 1, size
	 do j = 1, size
	    print *, "ROW ", i , " COL ", j
	    do k = 1, size
	       print *, ar(i,j,k)
            enddo
         enddo
      enddo
      end
c===============================================================|
      subroutine dz3d6p(dz)
c...............................................................|
c
c   routine:
c      computes the z-derivatives of u storing the result in ud
c      uses a 6th order Pade scheme (Periodic domain)
c      periodic in z (evenly spaced grid)
c      third index corresponds to z direction
c
c...............................................................|
c

      parameter (pi=3.1415926535898)
c
      include 'erle.h'
      double precision a(N),b(N),c(N),d(N),e(N)

c
c...............................................................|
c
      nij=N*N
c
c  first derivative
c
         c13    = 1./3.
         c79dz  = 7./(9.*dz)
         c136dz = 1./(36.*dz)
c
c  set up the left hand side and factor the periodic matrix
c
         do 30 i=1,N
            a(i) = c13
            b(i) = 1.
            c(i) = c13
   30    continue
         call tridvpf(a,b,c,d,e)
c
c  set up right hand side and do forward/backward substitution
c
   91    do 41 j=1,N
         do 41 i=1,N
            duz(i,j,1) = c79dz*(u(i,j,2) - u(i,j,N))
     1                 + c136dz*(u(i,j,3) - u(i,j,N-1))
            duz(i,j,2) = c79dz*(u(i,j,3) - u(i,j,1))
     1                 + c136dz*(u(i,j,4) - u(i,j,N))
   41    continue

   92  do 40 j=1,N
         do 40 i=1,N
           duz(i,j,N-1) = c79dz*(u(i,j,N) - u(i,j,N-2))
     1                    + c136dz*(u(i,j,1) - u(i,j,N-3))
            duz(i,j,N)   = c79dz*(u(i,j,1) - u(i,j,N-1))
     1                    + c136dz*(u(i,j,2) - u(i,j,N-2))
   40    continue

         do 50 k=3,N-2
         do 50 j=1,N
         do 50 i=1,N
              duz(i,j,k) = c79dz*(u(i,j,k+1) - u(i,j,k-1))
     1                   + c136dz*(u(i,j,k+2) - u(i,j,k-2))
   50    continue
         call tridvpk(a,b,c,d,e,wrk)
c
      return
      end
c===============================================================|
      subroutine setgrid
c
      common /grids/ dx,dy,dz

      include 'erle.h'
c
c...............................................................|
c
c calculate boundary indices for data section on each node
c

c
c calculate grid
c
      pi=acos(-1.0)
      dx=2.*pi/N
      dy=2.*pi/N
      dz=2.*pi/N
c
      do 5 i=1,N
         x(i)=(i-1)*dx
    5 continue
c
      do 6 j=1,N
         y(j)=(j-1)*dy
    6 continue
c
      do 7 k=1,N
         z(k)=(k-1)*dz
    7 continue
c
      return
      end
c===============================================================|
      subroutine dotest

      common /grids/ dx,dy,dz
      common /times/ dt,dtx,dty,dtz
      double precision dt,dtx,dty,dtz

c...............................................................|
c
c*** generate data
c
      call genvar
c
c...............................................................|
c
c*** time derivative calculation
c
      call dx3d6p(dx)
      call dy3d6p(dy)
      call dz3d6p(dz)

      return
      end
c===============================================================|
      subroutine genvar
c
      common /grids/ dx,dy,dz
      include 'erle.h'

      do 10 kk=1,N
         do 10 jj=1,N
            do 10 ii=1,N
               u(ii,jj,kk)   =sin(x(ii))*cos(y(jj))*cos(3.*z(kk))
   10 continue
c
      return
      end
c===============================================================|
      subroutine dx3d6p(dx)
c...............................................................|
c
c   routine:
c      computes the x-derivatives of u storing the result in ud
c      uses a 6th order Pade scheme (Periodic domain)
c      periodic in x (evenly spaced grid)
c      first index corresponds to x direction
c
c...............................................................|
c
      parameter (pi=3.1415926535898)
c
      include 'erle.h'
        double precision a(N),b(N),c(N),d(N),e(N)

c...............................................................|
c
c  first derivative
c
         c13    = 1./3.
         c79dx  = 7./(9.*dx)
         c136dx = 1./(36.*dx)
c
c  set up the left hand side and factor the periodic matrix
c
         do 30 i=1,N
            a(i) = c13
            b(i) = 1.
            c(i) = c13
   30    continue
         call tridvpf(a,b,c,d,e)
c
c  set up right hand side and do forward/backward substitution
c
         do 40 k=1,N
         do 40 j=1,N
            dux(1,j,k) = c79dx*(u(2,j,k) - u(N,j,k))
     &               + c136dx*(u(3,j,k) - u(N-1,j,k))
            dux(2,j,k) = c79dx*(u(3,j,k) - u(1,j,k))
     &               + c136dx*(u(4,j,k) - u(N,j,k))
   40    continue
c
         do 41 k=1,N
         do 41 j=1,N
            dux(N-1,j,k) = c79dx*(u(N,j,k) - u(N-2,j,k))
     &                     + c136dx*(u(1,j,k)    - u(N-3,j,k))
            dux(N,j,k) = c79dx*(u(1,j,k) - u(N-1,j,k))
     &                   + c136dx*(u(2,j,k) - u(N-2,j,k))
   41    continue
c
         do 50 k=1,N
         do 50 j=1,N
         do 50 i=3,N-2
              dux(i,j,k) = c79dx*(u(i+1,j,k) - u(i-1,j,k))
     &                 + c136dx*(u(i+2,j,k) - u(i-2,j,k))
   50    continue
         call tridvpi(a,b,c,d,e,wrk)
      return
      end
c===============================================================|
      subroutine dy3d6p(dy)
c...............................................................|
c
c   routine:
c      computes the y-derivatives of u storing the result in ud
c      uses a 6th order Pade scheme (Periodic domain)
c      periodic in y (evenly spaced grid)
c      second index corresponds to y direction
c
c...............................................................|
c

      parameter (pi=3.1415926535898)
c
      include 'erle.h'

        double precision a(N),b(N),c(N),d(N),e(N)
c...............................................................|
c
c  first derivative
c
         c13    = 1./3.
         c79dy  = 7./(9.*dy)
         c136dy = 1./(36.*dy)
c
c  set up the left hand side and factor the periodic matrix
c
         do 30 i=1,N
            a(i) = c13
            b(i) = 1.
            c(i) = c13
  30     continue
         call tridvpf(a,b,c,d,e)
c
c  set up right hand side and do forward/backward substitution
c
         do 40 k=1,N
         do 40 i=1,N
            duy(i,1,k) = c79dy*(u(i,2,k) - u(i,N,k))
     1               + c136dy*(u(i,3,k) - u(i,N-1,k))
            duy(i,2,k) = c79dy*(u(i,3,k) - u(i,1,k))
     1               + c136dy*(u(i,4,k) - u(i,N,k))
  40     continue

         do 41 k=1,N
         do 41 i=1,N
            duy(i,N-1,k) = c79dy*(u(i,N,k)
     1                          - u(i,N-2,k))
     2                     + c136dy*(u(i,1,k) - u(i,N-3,k))
            duy(i,N,k) = c79dy*(u(i,1,k) - u(i,N-1,k))
     1                   + c136dy*(u(i,2,k) - u(i,N-2,k))
  41     continue
c
         do 50 k=1,N
         do 50 j=3,N-2
         do 50 i=1,N
              duy(i,j,k) = c79dy*(u(i,j+1,k) - u(i,j-1,k))
     1                 + c136dy*(u(i,j+2,k) - u(i,j-2,k))
  50     continue
         call tridvpj(a,b,c,d,e,wrk)
c
      return
      end

max_stars_repo_name: Artoriasviel/POET

max_issues_repo_name: Artoriasviel/POET

max_forks_repo_name: DeveloperPOET/POET


text:
program int1list_example
  use type_int1list
  implicit none

  type(Int1List), allocatable :: myList

  ! Explicit allocation and deallocation to test finalization
  allocate(myList)
  call Int1List_init(myList)
  print *, 'Creating list [[1, 2], [2, 3, 4], [5, 6, 7, 8]]'
  call myList%append([1, 2])
  call myList%append([2, 3, 4])
  call myList%append([5, 6, 7, 8])
  print *, 'Size of the list:', size(myList)   ! Should yield 3
  print *, 'Second element:', myList%get(2)  ! Should yield 2, 3, 4
  print *, 'Position of item with value [2, 3, 4]:', myList%find([2, 3, 4])  ! Should yield 2
  print *, 'Position of item with value [1, 1]:', myList%find([1, 1])  ! Should yield 0

  print *, 'Inserting [9, 10] at position 2'
  call myList%insert(2, [9, 10])
  print *, 'List size:', size(myList)  ! 4
  print *, 'Elements 1:3:'
  print *, '1:', myList%get(1)  ! 1, 2
  print *, '2:', myList%get(2)  ! 9, 10
  print *, '3:', myList%get(3)  ! 2, 3, 4
  deallocate(myList)

end program int1list_example

ext: fpp

max_stars_repo_name: aradi/fycollections

max_issues_repo_name: aradi/fycollections

max_forks_repo_name: aradi/fycollections


text:
Describe Users/KelseyCondell here.
Hi Kelsey, and Welcome to the Wiki/Business Owner welcome to the Wiki!  Thanks for editing under your real name.  Identity is a big plus for interacting with others on the wiki.  Users/TomGarberson

20110523 13:55:29 nbsp Thanks! I actually feel like a shmuck for editing in my real name, because most people dont seem to.  Users/KelseyCondell

20110707 11:20:48 nbsp Hi Kelsey, I was wondering what the internet is REALLY like at Cambridge House. Theres a lot of haters on the wiki, and its hard to extract the truth from it. You said it worked well for you, but were you streaming movies off it, or what? Thanks! Users/RyanPhillips

20120218 20:14:29 nbsp aw, RyanPhillips, I just got your message! Doh!! The internet worked great for me... but they screwed me over at move out. surprise surprise! Users/KelseyCondell


!***********************************/
!	Name:物体まわりの点を並び変えて格納する
!	Alias:UMakeInternalObject
!	Description:
!	Type:UnstructuredGrid
!	Input:
!	Output:
!	Note:
!	Author:Akitaka Toyota
!	Date:2018.02.08
!	Update:
!	Other:
!***********************************/
subroutine UMakeInternalObject(UG)
    use StructVar_Mod
    use LoopVar_Mod
    use FrequentOperation
    implicit none
    type(UnstructuredGrid), intent(inout) :: UG
    integer, allocatable :: Point2(:),Point(:)
    double precision, allocatable :: Angle(:),Norm(:)
    double precision, allocatable :: RelativeVector(:)
    double precision :: ErrorNorm, MinYCD, LastMinYCD
    integer :: iCount, iProcess, iEndFlag, iCheckMin

    iCount = 0
    do iCell = UG%GI%RealCells+1, UG%GI%AllCells
        if(UG%VC%Type(iCell) == 2) then
            iCount = iCount + 1
        end if
    end do

    allocate(Point(iCount))
    allocate(Point2(iCount*2))

    iCount = 0
    do iCell = UG%GI%RealCells+1, UG%GI%AllCells
        if(UG%VC%Type(iCell) == 2) then
            do iLoop=1,2
                iCount = iCount + 1
                Point2(iCount) = UG%Line%Point(UG%VC%Edge(iCell),iLoop)

            end do
        end if
    end do

    iEndFlag = 0
    iProcess = 0
    do while(iEndFlag == 0)
        iEndFlag = 1
        iProcess = iProcess + 1
        do iPoint=1, ubound(Point2,1)-iProcess

            if(Point2(iPoint) > Point2(iPoint+1)) then
                iEndFlag = 0
                call SwapInteger(Point2(iPoint),Point2(iPoint+1))
            end if
        end do
    end do

    MinYCD = 10**30

    do iPoint = 1, ubound(Point,1)
        Point(iPoint) = Point2(2*iPoint)

        LastMinYCD = MinYCD
        MinYCD = min(MinYCD,UG%CD%Point(Point(iPoint),2))
        if(LastMinYCD /= MinYCD) iCheckMin = Point(iPoint)
    end do
        if(LastMinYCD /= MinYCD) iCheckMin = Point(ubound(Point,1)) !Final Step

    deallocate(Point2)
    allocate(Angle(ubound(Point,1)),Norm(ubound(Point,1)))
    allocate(RelativeVector(3))

    do iPoint = 1, ubound(Point,1)
        if(iPoint /= iCheckMin) then
            RelativeVector = UG%CD%Point(Point(iPoint),1:3) - UG%CD%Point(iCheckMin,1:3)
            Angle(iPoint) = atan2(RelativeVector(2),RelativeVector(1))
            Norm(iPoint) = dot_product(RelativeVector(1:3),RelativeVector(1:3))

        else
            Angle(iPoint) = 0.0d0
            Norm(iPoint) = 0.0d0
        end if
    end do

    ErrorNorm = 10**(-10)

    iEndFlag = 0
    iProcess = 0
    do while(iEndFlag == 0)
        iEndFlag = 1
        iProcess = iProcess + 1
        do iPoint=1, ubound(Point,1)-iProcess
            if(Angle(iPoint)-ErrorNorm < Angle(iPoint) .and.  Angle(iPoint) < Angle(iPoint)+ErrorNorm) then
                if(Norm(iPoint) > Norm(iPoint + 1)) then
                    iEndFlag = 0
                    call swapAll

                end if

            else if(Angle(iPoint) < Angle(iPoint+1)) then
                iEndFlag = 0
                call swapAll
            end if
        end do
    end do

    allocate(UG%IO%PointNum(ubound(Point,1)))

    UG%IO%iTotal = ubound(Point,1)
    UG%IO%PointNum = Point

    return
contains
    subroutine swapAll
    implicit none
        call SwapInteger(Point(iPoint),Point(iPoint+1))
        call SwapDoublePrecision(Angle(iPoint),Angle(iPoint+1))
        call SwapDoublePrecision(Norm(iPoint),Norm(iPoint+1))
        return
    end subroutine swapAll

end subroutine UMakeInternalObject


text:
************************************************************************
* This file implements basic composed geometries as a set of
* simple geometries. It provides the basic functionality like
* in the case of simple geometries but tries to augment it in this
* more general case.
*
* Like in the case of simple geometries the routines provided here
* separate into the following types of routines:
*
* 1.) Test if a point is inside of the geometry
* 2.) Write an approximation of the boundary into a GMV-file.
* 3.) Calculate the volume of the geometry
* 4.) Calculate the length of the boundary interface ("arc length")
* 5.) Calculate the outer normal vector of the geometry in a point on
*     the boundary interface
* 6.) Project a point onto the boundary
* 7.) Calculate the minimum distance of an arbitrary point to the
*     object.
*
* Each routine accepts (like described in SCOMPGEOMETRIES.INC)
* an index-integer-array and a double-array describing the composed
* geometry, as well as a wrapper function that calls the appropriate
* simple-geometry-routine according to the object type.
* It's acceptable to specify the appropriate GxxWRP as wrapper routines
* for a basic functionality of simple-geometry wrappers.
************************************************************************

************************************************************************
* GxxWRP Wrapper-routines
*
* The following wrapper routines are designed to map the TYPE
* specifier in the object structure to a call to the appropriate
* routine defined in GEOMETRIES.F. They share the same syntax
* as the "simple"-routines and are fully interchangeable with them.
************************************************************************

*-----------------------------------------------------------------------
* Test if the point (X,Y) is inside of the geometry xxx.
*-----------------------------------------------------------------------

      INTEGER FUNCTION GINWRP (GEOSIM,CORSYS, X, Y)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      DOUBLE PRECISION GEOSIM(*),X,Y
      DOUBLE PRECISION CORSYS(SZCORS)
      INTEGER TP

      INTEGER  GINCIR,GINELL,GINSQR,GINREC,GINSWV,GINPSL
      EXTERNAL GINCIR,GINELL,GINSQR,GINREC,GINSWV,GINPSL

      TP = GEOSIM(OTYPE)

      GINWRP = 0
      IF (TP.EQ.SGTCIR) GINWRP = GINCIR(GEOSIM,CORSYS,X,Y)
      IF (TP.EQ.SGTELL) GINWRP = GINELL(GEOSIM,CORSYS,X,Y)
      IF (TP.EQ.SGTSQR) GINWRP = GINSQR(GEOSIM,CORSYS,X,Y)
      IF (TP.EQ.SGTREC) GINWRP = GINREC(GEOSIM,CORSYS,X,Y)
      IF (TP.EQ.SGTSWV) GINWRP = GINSWV(GEOSIM,CORSYS,X,Y)
      IF (TP.EQ.SGTPSL) GINWRP = GINPSL(GEOSIM,CORSYS,X,Y)

      END

*-----------------------------------------------------------------------
* Calculate the volume of a geometry object.
*-----------------------------------------------------------------------

      DOUBLE PRECISION FUNCTION GVLWRP (GEOSIM,CORSYS)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      DOUBLE PRECISION GEOSIM(*)
      DOUBLE PRECISION CORSYS(SZCORS)
      INTEGER TP

      DOUBLE PRECISION GVLCIR,GVLELL,GVLSQR,GVLREC,GVLSWV,GVLPSL
      EXTERNAL GVLCIR,GVLELL,GVLSQR,GVLREC,GVLSWV,GVLPSL

      TP = GEOSIM(OTYPE)

      GVLWRP = -1
      IF (TP.EQ.SGTCIR) GVLWRP = GVLCIR(GEOSIM,CORSYS)
      IF (TP.EQ.SGTELL) GVLWRP = GVLELL(GEOSIM,CORSYS)
      IF (TP.EQ.SGTSQR) GVLWRP = GVLSQR(GEOSIM,CORSYS)
      IF (TP.EQ.SGTREC) GVLWRP = GVLREC(GEOSIM,CORSYS)
      IF (TP.EQ.SGTSWV) GVLWRP = GVLSWV(GEOSIM,CORSYS)
      IF (TP.EQ.SGTPSL) GVLWRP = GVLPSL(GEOSIM,CORSYS)

      END

*-----------------------------------------------------------------------
* Calculate the length of the boundary interface
*-----------------------------------------------------------------------

      DOUBLE PRECISION FUNCTION GALWRP (GEOSIM,CORSYS)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      DOUBLE PRECISION GEOSIM(*)
      DOUBLE PRECISION CORSYS(SZCORS)
      INTEGER TP

      DOUBLE PRECISION GALCIR,GALELL,GALSQR,GALREC,GALSWV,GALPSL
      EXTERNAL GALCIR,GALELL,GALSQR,GALREC,GALSWV,GALPSL

      TP = GEOSIM(OTYPE)

      GALWRP = -1
      IF (TP.EQ.SGTCIR) GALWRP = GALCIR(GEOSIM,CORSYS)
      IF (TP.EQ.SGTELL) GALWRP = GALELL(GEOSIM,CORSYS)
      IF (TP.EQ.SGTSQR) GALWRP = GALSQR(GEOSIM,CORSYS)
      IF (TP.EQ.SGTREC) GALWRP = GALREC(GEOSIM,CORSYS)
      IF (TP.EQ.SGTSWV) GALWRP = GALSWV(GEOSIM,CORSYS)
      IF (TP.EQ.SGTPSL) GALWRP = GALPSL(GEOSIM,CORSYS)

      END

*-----------------------------------------------------------------------
* Write an approximation of the boundary interface into a GMV-file
* identified by the file handle MUNIT.
* Use material number MMAT to identify the coloring in the GMV-file.
*-----------------------------------------------------------------------

      SUBROUTINE GGMWRP(GEOSIM,CORSYS, MUNIT, MMAT)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      DOUBLE PRECISION GEOSIM(*)
      DOUBLE PRECISION CORSYS(SZCORS)
      INTEGER MUNIT,MMAT
      INTEGER TP

      TP = GEOSIM(OTYPE)

      IF (TP.EQ.SGTCIR) CALL GGMCIR(GEOSIM,CORSYS, MUNIT, MMAT)
      IF (TP.EQ.SGTELL) CALL GGMELL(GEOSIM,CORSYS, MUNIT, MMAT)
      IF (TP.EQ.SGTSQR) CALL GGMSQR(GEOSIM,CORSYS, MUNIT, MMAT)
      IF (TP.EQ.SGTREC) CALL GGMREC(GEOSIM,CORSYS, MUNIT, MMAT)
      IF (TP.EQ.SGTSWV) CALL GGMSWV(GEOSIM,CORSYS, MUNIT, MMAT)
      IF (TP.EQ.SGTPSL) CALL GGMPSL(GEOSIM,CORSYS, MUNIT, MMAT)

      END

*-----------------------------------------------------------------------
* Calculate the normalized outer normal vector in the boundary-point
* (X,Y) on the interface of the geometry object GEOSIM.
*-----------------------------------------------------------------------

      SUBROUTINE GNMWRP (GEOSIM,CORSYS, X,Y, XN,YN)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      DOUBLE PRECISION GEOSIM(*),X,Y,XN,YN
      DOUBLE PRECISION CORSYS(SZCORS)
      INTEGER TP

      TP = GEOSIM(OTYPE)

      IF (TP.EQ.SGTCIR) CALL GNMCIR(GEOSIM,CORSYS, X,Y, XN,YN)
      IF (TP.EQ.SGTELL) CALL GNMELL(GEOSIM,CORSYS, X,Y, XN,YN)
      IF (TP.EQ.SGTSQR) CALL GNMSQR(GEOSIM,CORSYS, X,Y, XN,YN)
      IF (TP.EQ.SGTREC) CALL GNMREC(GEOSIM,CORSYS, X,Y, XN,YN)
      IF (TP.EQ.SGTSWV) CALL GNMSWV(GEOSIM,CORSYS, X,Y, XN,YN)
      IF (TP.EQ.SGTPSL) CALL GNMPSL(GEOSIM,CORSYS, X,Y, XN,YN)

      END

*-----------------------------------------------------------------------
* Calculate the distance of the point (X,Y) to the boundary interface
* of the object.
*-----------------------------------------------------------------------

      INTEGER FUNCTION GDSWRP (GEOSIM,CORSYS, X,Y, DIST)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      DOUBLE PRECISION GEOSIM(*),X,Y,DIST
      DOUBLE PRECISION CORSYS(SZCORS)
      INTEGER TP

      INTEGER  GDSCIR,GDSELL,GDSSQR,GDSREC,GDSSWV,GDSPSL
      EXTERNAL GDSCIR,GDSELL,GDSSQR,GDSREC,GDSSWV,GDSPSL

      TP = GEOSIM(OTYPE)

      GDSWRP = 0
      IF (TP.EQ.SGTCIR) GDSWRP = GDSCIR(GEOSIM,CORSYS, X,Y, DIST)
      IF (TP.EQ.SGTELL) GDSWRP = GDSELL(GEOSIM,CORSYS, X,Y, DIST)
      IF (TP.EQ.SGTSQR) GDSWRP = GDSSQR(GEOSIM,CORSYS, X,Y, DIST)
      IF (TP.EQ.SGTREC) GDSWRP = GDSREC(GEOSIM,CORSYS, X,Y, DIST)
      IF (TP.EQ.SGTSWV) GDSWRP = GDSSWV(GEOSIM,CORSYS, X,Y, DIST)
      IF (TP.EQ.SGTPSL) GDSWRP = GDSPSL(GEOSIM,CORSYS, X,Y, DIST)

      END

*-----------------------------------------------------------------------
* Project a point onto the boundary interface of the geometry object.
*-----------------------------------------------------------------------

      INTEGER FUNCTION GPRWRP (GEOSIM,CORSYS, X,Y, XP,YP)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      DOUBLE PRECISION GEOSIM(*),X,Y,XP,YP
      DOUBLE PRECISION CORSYS(SZCORS)
      INTEGER TP

      INTEGER  GPRCIR,GPRELL,GPRSQR,GPRREC,GPRSWV,GPRPSL
      EXTERNAL GPRCIR,GPRELL,GPRSQR,GPRREC,GPRSWV,GPRPSL

      TP = GEOSIM(OTYPE)

      GPRWRP = 0
      IF (TP.EQ.SGTCIR) GPRWRP = GPRCIR(GEOSIM,CORSYS, X,Y, XP,YP)
      IF (TP.EQ.SGTELL) GPRWRP = GPRELL(GEOSIM,CORSYS, X,Y, XP,YP)
      IF (TP.EQ.SGTSQR) GPRWRP = GPRSQR(GEOSIM,CORSYS, X,Y, XP,YP)
      IF (TP.EQ.SGTREC) GPRWRP = GPRREC(GEOSIM,CORSYS, X,Y, XP,YP)
      IF (TP.EQ.SGTSWV) GPRWRP = GPRSWV(GEOSIM,CORSYS, X,Y, XP,YP)
      IF (TP.EQ.SGTPSL) GPRWRP = GPRPSL(GEOSIM,CORSYS, X,Y, XP,YP)

      END

************************************************************************
* The following set of routines now define the functionality
* of the composed-geometry objects. The same set as in the case of
* simple geometries is provided, but the calling conventions are
* different.
************************************************************************

************************************************************************
* INTEGER FUNCTION GINCMP (GEOIDX, GEOCOM, GINFCT, X, Y)
*
* Test if the point (X,Y) is inside of the geometry xxx.
*
* In:
*  X,Y    - coordinate of the point
*  GEOIDX - the GEOCOMIDX integer array describing the index positions
*           of the objects in GEOCOM
*  GEOCOM - the GEOCOM double array describing the sub-objects
*  GINFCT - a wrapper function that provides calling the correct
*           GINxxx-routine. GINWRP is acceptable.
*
* Out:
*  Result > 0, if the point is inside of the geometry object.
*              If the group-object is not completely inverted (INVER=1),
*              the number describes the number of sub-objects
*              that contain this point.
*         = 0, if the point is not inside of the geometry object
*         < 0, if the point is inside of the geometry object and
*              the object is inverted
************************************************************************

      INTEGER FUNCTION GINCMP (GEOIDX, GEOCOM, GINFCT, X, Y)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      INCLUDE 'scompgeometries.inc'
      DOUBLE PRECISION GEOCOM(*),X,Y
      INTEGER GEOIDX(*),GINFCT
      EXTERNAL GINFCT

      INTEGER I,CNT

      CNT = 0

C Loop through all sub-objects

      DO I=1,GEOIDX(OGCCNT)

C Add the results of the GINxxx-routines. Negative values
C decrease the number until it's <= 0 - that indicates
C that there are so many objects substracted from the first
C object, that the point is considered to be "outside".

        CNT = CNT + GINFCT(GEOCOM(GEOIDX(OGCIDX+I-1)),
     *      GEOCOM(OGCMAIN+OCSYS-1), X, Y)

      END DO

      GINCMP = CNT

C Switch the result if we are in "inverted" mode

      IF (GEOCOM(OINVER).NE.0) THEN
        GINCMP = -GINCMP
      END IF

      END

************************************************************************
* DOUBLE PRECISION FUNCTION GVLCMP (GEOIDX, GEOCOM, GVLFCT)
*
* Calculate sum of all volume of a geometry object.
*
* In:
*  GEOIDX - the GEOCOMIDX integer array describing the index positions
*           of the objects in GEOCOM
*  GEOCOM - the GEOCOM double array describing the sub-objects
*  GVLFCT - a wrapper function that provides calling the correct
*           GVLxxx-routine. GVLWRP is acceptable.
*
* Out:
*  Result = sum of all volumes of the sub-objects
*         = -1, if the volume cannot be calculated
************************************************************************

      DOUBLE PRECISION FUNCTION GVLCMP (GEOIDX, GEOCOM, GVLFCT)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      INCLUDE 'scompgeometries.inc'
      DOUBLE PRECISION GEOCOM(*)
      INTEGER GEOIDX(*)
      DOUBLE PRECISION GVLFCT
      EXTERNAL GVLFCT

      DOUBLE PRECISION VOL
      INTEGER I

      VOL = 0D0

C Loop through all sub-objects

      DO I=1,GEOIDX(OGCCNT)

C sum up the volumes

        VOL = VOL + GVLFCT(GEOCOM(GEOIDX(OGCIDX+I-1)),
     *                     GEOCOM(OGCMAIN+OCSYS-1) )

      END DO

      GVLCMP = VOL

      END

************************************************************************
* DOUBLE PRECISION FUNCTION GALCMP (GEOIDX, GEOCOM, GALFCT)
*
* Calculate the length of the boundary interface
*
* In:
*  GEOIDX - the GEOCOMIDX integer array describing the index positions
*           of the objects in GEOCOM
*  GEOCOM - the GEOCOM double array describing the sub-objects
*  GALFCT - a wrapper function that provides calling the correct
*           GALxxx-routine. GALWRP is acceptable.
*
* Out:
*  Result = sum of all lenths of the interfaces of the sub-objects
*         = -1, if the length cannot be calculated
************************************************************************

      DOUBLE PRECISION FUNCTION GALCMP (GEOIDX, GEOCOM, GALFCT)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      INCLUDE 'scompgeometries.inc'
      DOUBLE PRECISION GEOCOM(*)
      INTEGER GEOIDX(*)
      DOUBLE PRECISION GALFCT
      EXTERNAL GALFCT

      DOUBLE PRECISION ALE
      INTEGER I

      ALE = 0D0

C Loop through all sub-objects

      DO I=1,GEOIDX(OGCCNT)

C sum up the arc lengths

        ALE = ALE + GALFCT(GEOCOM(GEOIDX(OGCIDX+I-1)),
     *                     GEOCOM(OGCMAIN+OCSYS-1) )

      END DO

      GALCMP = ALE

      END

************************************************************************
* SUBROUTINE GGMCMP (GEOIDX, GEOCOM, GGMFCT, GDSFCT, MUNIT, MMAT)
*
* Write an approximation of the boundary interface into a GMV-file
* identified by the file handle MUNIT.
* Use material number MMAT to identify the coloring in the GMV-file.
*
* The caller of this routine is responsible for the administration
* of the sections in the GMV file. This routine only has to write out
* the content for the POLYGON's section to approximate the boundary
* components. The header/footer of this section is written to the
* file by the caller.
*
* In:
*  GEOIDX - the GEOCOMIDX integer array describing the index positions
*           of the objects in GEOCOM
*  GEOCOM - the GEOCOM double array describing the sub-objects
*  GGMFCT - a wrapper function that provides calling the correct
*           GGMxxx-routine. GGMWRP is acceptable.
*  MUNIT  - File unit where to write the GMV output to.
*  MMAT   - The first material number to use.
************************************************************************

      SUBROUTINE GGMCMP (GEOIDX, GEOCOM, GGMFCT, MUNIT, MMAT)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      INCLUDE 'scompgeometries.inc'
      DOUBLE PRECISION GEOCOM(*)
      INTEGER GEOIDX(*),MUNIT, MMAT
      EXTERNAL GGMFCT

      INTEGER I

C Loop through all sub-objects

      DO I=1,GEOIDX(OGCCNT)

C tell the sub-object to write itself to disc

        CALL GGMFCT(GEOCOM(GEOIDX(OGCIDX+I-1)),
     *              GEOCOM(OGCMAIN+OCSYS-1),
     *              MUNIT, MMAT)

      END DO

      END

************************************************************************
* SUBROUTINE GNMCMP (GEOIDX, GEOCOM, GNMFCT, X,Y, XN,YN)
*
* Calculate the normalized outer normal vector in the boundary-point
* (X,Y) on the interface of the geometry object GEOSIM.
*
* In:
*  GEOIDX - the GEOCOMIDX integer array describing the index positions
*           of the objects in GEOCOM
*  GEOCOM - the GEOCOM double array describing the sub-objects
*  GNMFCT - a wrapper function that provides calling the correct
*           GGMxxx-routine. GNMWRP is acceptable.
*  GDSFCT - a wrapper function that provides calling the correct
*           GDSxxx-routine for distance calculation. GDSWRP is
*           acceptable.
*  X,Y    - coordinates of a point on the interface where to calculate
*           the outer normal vector
*
* Out:
*  XN,YN  - Outer normal vector in the point (X,Y),
*           normalized to length = 1
*           =(0,0) if the calculation is not possible.
*
* This routine only works completely correct for points outside of
* the geometry. Points inside of the geometry may be projected
* wrong if they are lying in more than one sub-component!
************************************************************************

      SUBROUTINE GNMCMP (GEOIDX, GEOCOM, GNMFCT, GDSFCT, X,Y, XN,YN)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      INCLUDE 'scompgeometries.inc'
      DOUBLE PRECISION GEOCOM(*),X,Y,XN,YN
      INTEGER GEOIDX(*),GDSFCT
      EXTERNAL GNMFCT,GDSCMP,GDSFCT
      INTEGER GDSCMP

      DOUBLE PRECISION DIST
      INTEGER I

      XN = 0D0
      YN = 0D0

C Calculate the sub-object with the minimum distance

      I = GDSCMP (GEOIDX, GEOCOM, GDSFCT, X,Y, DIST)

C Calculate the normal vector by that calculation routine

      IF (I.NE.0) THEN
        CALL GNMFCT(GEOCOM(GEOIDX(OGCIDX+I-1)),
     *              GEOCOM(OGCMAIN+OCSYS-1),
     *              X,Y, XN,YN)
      END IF

C Switch the result if we are in "inverted" mode

      IF (GEOCOM(OINVER).NE.0) THEN
        XN = -XN
        YN = -YN
      END IF

      END

************************************************************************
* INTEGER FUNCTION GDSCMP (GEOIDX, GEOCOM, GDSFCT, X,Y, DIST)
*
* Calculate the distance of the point (X,Y) to the boundary interface
* of the object.
*
* In:
*  GEOIDX - the GEOCOMIDX integer array describing the index positions
*           of the objects in GEOCOM
*  GEOCOM - the GEOCOM double array describing the sub-objects
*  GDSFCT - a wrapper function that provides calling the correct
*           GGMxxx-routine. GDSWRP is acceptable.
*  X,Y    - the point whose distance should be calculated
*
* Out:
*  DIST   - min. distance of (X,Y) to the boundary interfaces; the sign
*           identifies the position of the point:
*           > 0: the point is outside of the geometry object
*           <=0: the point is inside of the geometry object
*
*  Result: > 0, if the distance was calculated successfully.
*               The number describes the number of the sub-component
*               that served for calculating the distance.
*          <=0, if the computation was not possible
************************************************************************

      INTEGER FUNCTION GDSCMP (GEOIDX, GEOCOM, GDSFCT, X,Y, DIST)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      INCLUDE 'scompgeometries.inc'
      DOUBLE PRECISION GEOCOM(*),X,Y,DIST
      INTEGER GEOIDX(*)
      INTEGER GDSFCT
      EXTERNAL GDSFCT

      DOUBLE PRECISION DST
      INTEGER I,J

      GDSCMP = 0

C Loop through all sub-objects

      DO I=1,GEOIDX(OGCCNT)

C Calculate the distance to all objects that allow calculating the
C distance.

        J = GDSFCT(GEOCOM(GEOIDX(OGCIDX+I-1)),
     *             GEOCOM(OGCMAIN+OCSYS-1),
     *             X,Y,DST)

        IF (J.GT.0) THEN

C Build the minimum of the distances.
C Note: For points outside of the geometry this is the distance of the
C  point to the geometry. For points inside the geometry, the minimum
C  is build from the *negative* distances: This results in the maximum
C  distance to all objects containing the point. This is of course not
C  the real distance to the boundary interface of the composed geometry,
C  but somehow an approximation to it.

          IF ((I.EQ.1).OR.(DST.LT.DIST)) THEN
            GDSCMP = I
            DIST = DST
          END IF

        END IF

      END DO

C Switch the result if we are in "inverted" mode

      IF (GEOCOM(OINVER).NE.0) THEN
        DIST = -DIST
      END IF

      END

************************************************************************
* INTEGER FUNCTION GPRCMP (GEOIDX, GEOCOM, GPRFCT, X,Y, XP,YP)
*
* Project a point onto the boundary interface of the geometry object.
*
* In:
*  GEOIDX - the GEOCOMIDX integer array describing the index positions
*           of the objects in GEOCOM
*  GEOCOM - the GEOCOM double array describing the sub-objects
*  GPRFCT - a wrapper function that provides calling the correct
*           GPRxxx-routine. GPRWRP is acceptable.
*  GDSFCT - a wrapper function that provides calling the correct
*           GDSxxx-routine for distance calculation. GDSWRP is
*           acceptable.
*  X,Y    - the point that should be projected
*
* Out:
*  XP,YP  - the projected point on the interface of GEOSIM
*
*  Result: > 0, if the projection was calculated successfully.
*               The number represents the number of the sub-component
*               onto which the projection was done.
*          <=0, if the computation was not possible
*
* This routine only works completely correct for points outside of
* the geometry. Points inside of the geometry may be projected
* wrong if they are lying in more than one sub-component!
************************************************************************

      INTEGER FUNCTION GPRCMP (GEOIDX, GEOCOM, GPRFCT, GDSFCT,
     *                         X,Y, XP,YP)
      IMPLICIT NONE
      INCLUDE 'sgeometries.inc'
      INCLUDE 'scompgeometries.inc'
      DOUBLE PRECISION GEOCOM(*),X,Y,XP,YP
      INTEGER GEOIDX(*)
      INTEGER GDSCMP,GPRFCT,GDSFCT
      EXTERNAL GPRFCT,GDSCMP,GDSFCT

      DOUBLE PRECISION DIST
      INTEGER I

C Calculate the sub-object with the minimum distance

      I = GDSCMP (GEOIDX, GEOCOM, GDSFCT, X,Y, DIST)

C Project onto that boundary.

      IF (I.NE.0) THEN
        I = GPRFCT(GEOCOM(GEOIDX(OGCIDX+I-1)),
     *             GEOCOM(OGCMAIN+OCSYS-1),
     *             X,Y, XP,YP)
      END IF

      GPRCMP = I

      END

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


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

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

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

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

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

! End MIT license text.

      SUBROUTINE SURFACE_FIT ( NUM_FITS, NUM_COEFFS, XI, YI, WI, XO, YO, WO, DEB, MESSAGE, OUNT, POLY_PCT_ERR, PCT_ERR_MAX, IERR )

! Fits a 2D surface polynomial to a set of input values. The subroutine is coded to fit a surface of up to 2nd order:

!              WF(X,Y) = B(0) + B(1)*X + B(2)*Y + B(3)*XY + B(4)*X^2 + B(5)*Y^2

      USE PENTIUM_II_KIND, ONLY       :  BYTE, LONG, DOUBLE
      USE IOUNT1, ONLY                :  ERR, F04, F06, WRT_LOG
      USE SCONTR, ONLY                :  BLNK_SUB_NAM, FATAL_ERR
      USE TIMDAT, ONLY                :  TSEC
      USE CONSTANTS_1, ONLY           :  ZERO, ONE
      USE SUBR_BEGEND_LEVELS, ONLY    :  SURFACE_FIT_BEGEND
      USE LSQ_MYSTRAN

      USE SURFACE_FIT_USE_IFs

      IMPLICIT NONE

      CHARACTER(LEN=LEN(BLNK_SUB_NAM)):: SUBR_NAME = 'SURFACE_FIT'
      CHARACTER(LEN=*), INTENT(IN)    :: MESSAGE               ! Message printed if debug output is printed

      INTEGER(LONG), INTENT(IN)       :: DEB                   ! If > 0 then DEB_SURFACE_FIT will be run
      INTEGER(LONG), INTENT(IN)       :: NUM_FITS              ! Number of data points to fit
      INTEGER(LONG), INTENT(IN)       :: NUM_COEFFS            ! Number of coefficients in the fitting polynomial
      INTEGER(LONG), INTENT(IN)       :: OUNT(2)               ! Output units for SURFACE_FIT
      INTEGER(LONG), INTENT(OUT)      :: IERR                  ! Error indicator
      INTEGER(LONG)                   :: I,J                   ! DO loop indices
      INTEGER(LONG)                   :: IFAULT                ! Return code from subr REGCF
      INTEGER(LONG), PARAMETER        :: MAX_COEFFS = 6        ! Maximum number of coefficients coded for ther polynomial fit
      INTEGER(LONG), PARAMETER        :: SUBR_BEGEND = SURFACE_FIT_BEGEND

      LOGICAL                         :: LINDEP(MAX_COEFFS)

      REAL(DOUBLE), INTENT(IN)        :: WI(NUM_FITS)          ! Values of the function to fit at the input data points
      REAL(DOUBLE), INTENT(IN)        :: XI(NUM_FITS)          ! X coords of the input  data points
      REAL(DOUBLE), INTENT(IN)        :: YI(NUM_FITS)          ! Y coords of the input  data points
      REAL(DOUBLE), INTENT(IN)        :: XO(NUM_FITS)          ! X coords of the output data points
      REAL(DOUBLE), INTENT(IN)        :: YO(NUM_FITS)          ! Y coords of the output data points
      REAL(DOUBLE), INTENT(OUT)       :: POLY_PCT_ERR(NUM_FITS)! % difference between fitted & actual data (norm'd to max act data)
      REAL(DOUBLE), INTENT(OUT)       :: PCT_ERR_MAX           ! Max value from array POLY_PCT_ERR (max error at any input data pt)
      REAL(DOUBLE), INTENT(OUT)       :: WO(NUM_FITS)          ! Values of the function to fit at the output data points
      REAL(DOUBLE)                    :: DEN                   ! Intermediate value in a calculation
      REAL(DOUBLE)                    :: B(0:MAX_COEFFS-1)     ! Coefficients in the polynomial fit
      REAL(DOUBLE)                    :: RES(NUM_FITS)         ! Actual minus fitted data
      REAL(DOUBLE)                    :: RES_ABS               ! Max abs of the residual
      REAL(DOUBLE)                    :: WF(NUM_FITS)          ! Values at the input XI, YI coords of the fitted polynomial
      REAL(DOUBLE)                    :: WI_MAX                ! Max abs value of WI

      REAL(DOUBLE)                    :: XROW(0:MAX_COEFFS-1)! Array of XYI values for one data point

                                                             ! Values of XI, YI and their products to go into the poly fit eqn
      REAL(DOUBLE)                    :: XYI(NUM_FITS,0:MAX_COEFFS-1)

                                                             ! Values of XO, YO and their products to go into the poly fit eqn
      REAL(DOUBLE)                    :: XYO(NUM_FITS,0:MAX_COEFFS-1)

      REAL(DOUBLE), PARAMETER         :: WT = ONE            ! Parameter

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

! **********************************************************************************************************************************
      IERR = 0

! Make sure that we have not requested an impossible situation with NUM_COEFFS

      IF ((NUM_COEFFS < 1) .OR. (NUM_COEFFS > MAX_COEFFS)) THEN
         FATAL_ERR = FATAL_ERR + 1
         WRITE(OUNT(1),958) SUBR_NAME, NUM_COEFFS, MAX_COEFFS
         IF (OUNT(2) /= OUNT(1)) THEN
            WRITE(OUNT(2),958) SUBR_NAME, NUM_COEFFS, MAX_COEFFS
         ENDIF
         IERR = 1
         RETURN
      ENDIF

! Initialize

      DO I=1,NUM_FITS
         DO J=1,MAX_COEFFS
            XYI(I,J-1) = ZERO
            XROW(J-1)  = ZERO
            LINDEP(J) = .FALSE.
         ENDDO
      ENDDO

! Calc polynomial fit to the input WI at its coords, XI, YI

      CALL STARTUP (NUM_COEFFS-1, .TRUE.)

      WI_MAX = ZERO
      DO I=1,NUM_FITS
         IF (DABS(WI(I)) > WI_MAX) WI_MAX = WI(I)
         IF (NUM_COEFFS >= 1) THEN
            XYI(I,0) = ONE
         ENDIF
         IF (NUM_COEFFS >= 2) THEN
            XYI(I,1) = XI(I)
         ENDIF
         IF (NUM_COEFFS >= 3) THEN
            XYI(I,2) = YI(I)
         ENDIF
         IF (NUM_COEFFS >= 4) THEN
            XYI(I,3) = XI(I)*XI(I)
         ENDIF
         IF (NUM_COEFFS >= 5) THEN
            XYI(I,4) = XI(I)*YI(I)
         ENDIF
         IF (NUM_COEFFS >= 6) THEN
            XYI(I,5) = YI(I)*YI(I)
         ENDIF
         XROW(1:NUM_COEFFS-1) = XYI(I,1:NUM_COEFFS-1)
         XROW(0) = ONE
         CALL INCLUD ( WT, XROW, WI(I) )
      ENDDO

      CALL TOLSET ()

      CALL SING ( LINDEP, IFAULT )

      CALL REGCF ( B, NUM_COEFFS, IFAULT )

! Calc percent errors in polynomial fit to input data

      RES_ABS = ZERO
      DO I=1,NUM_FITS
         WF(I) = ZERO
         DO J=0,NUM_COEFFS-1
            WF(I) = WF(I) + B(J)*XYI(I,J)
         ENDDO
         RES(I) = WI(I) - WF(I)
         IF (DABS(RES(I)) > RES_ABS) THEN
            RES_ABS = DABS(RES(I))
         ENDIF
      ENDDO

      PCT_ERR_MAX = ZERO
      DO I=1,NUM_FITS
         IF (DABS(WI_MAX) > 0.0) THEN
            DEN = DABS(WI_MAX)
         ELSE
            DEN = ONE
         ENDIF
         POLY_PCT_ERR(I) = 100*RES(I)/DEN
         IF (DABS(PCT_ERR_MAX) < DABS(POLY_PCT_ERR(I))) THEN
            PCT_ERR_MAX = POLY_PCT_ERR(I)
         ENDIF
      ENDDO

! Calc fit to output coords, XO, YO

      DO I=1,NUM_FITS
         IF (NUM_COEFFS >= 1) THEN
            XYO(I,0)   = ONE
         ENDIF
         IF (NUM_COEFFS >= 2) THEN
            XYO(I,1)   = XO(I)
         ENDIF
         IF (NUM_COEFFS >= 3) THEN
            XYO(I,2)   = YO(I)
         ENDIF
         IF (NUM_COEFFS >= 4) THEN
            XYO(I,3)   = XO(I)*XO(I)
         ENDIF
         IF (NUM_COEFFS >= 5) THEN
            XYO(I,4)   = XO(I)*YO(I)
         ENDIF
         IF (NUM_COEFFS >= 6) THEN
            XYO(I,5)   = YO(I)*YO(I)
         ENDIF
         WO(I) = ZERO
         DO J=0,NUM_COEFFS-1
            WO(I) = WO(I) + B(J)*XYO(I,J)
         ENDDO
      ENDDO

      IF (DEB > 0) CALL DEB_SURFACE_FIT

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

      RETURN

! **********************************************************************************************************************************
  958 FORMAT(' *ERROR   958: PROGRAMMING ERROR IN SUBROUTINE ',A                                                                   &
                    ,/,14X,' THE VX VECTOR FOR ',A,' ELEMENT ',I8,' WAS LEFT UNSORTED. IT MUST BE SORTED TO DETERMINE VY, VZ')

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

      CONTAINS

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

      SUBROUTINE DEB_SURFACE_FIT

      IMPLICIT NONE

      INTEGER(LONG)                   :: II                ! DO loop index

! **********************************************************************************************************************************
      WRITE(OUNT(2),2000)
      WRITE(OUNT(2),'(A)') ' DEBUG 175 results from subr SURFACE_FIT'//MESSAGE//':'
      WRITE(OUNT(2),*)     ' --------------------------------------'
      WRITE(OUNT(2),2001) IFAULT
      WRITE(OUNT(2),'(A,1ES10.2,A)') ' PCT_ERR_MAX = ',PCT_ERR_MAX,'%'
      WRITE(OUNT(2),*)

      IF      (NUM_COEFFS == 1) THEN
         WRITE(OUNT(2),2011) NUM_FITS, NUM_COEFFS
         WRITE(OUNT(2),2021) (B(I), I=0,NUM_COEFFS-1)
      ELSE IF (NUM_COEFFS == 2) THEN
         WRITE(OUNT(2),2012) NUM_FITS, NUM_COEFFS
         WRITE(OUNT(2),2022) (B(I), I=0,NUM_COEFFS-1)
      ELSE IF (NUM_COEFFS == 3) THEN
         WRITE(OUNT(2),2013) NUM_FITS, NUM_COEFFS
         WRITE(OUNT(2),2023) (B(I), I=0,NUM_COEFFS-1)
      ELSE IF (NUM_COEFFS == 4) THEN
         WRITE(OUNT(2),2014) NUM_FITS, NUM_COEFFS
         WRITE(OUNT(2),2024) (B(I), I=0,NUM_COEFFS-1)
      ELSE IF (NUM_COEFFS == 5) THEN
         WRITE(OUNT(2),2014) NUM_FITS, NUM_COEFFS
         WRITE(OUNT(2),2025) (B(I), I=0,NUM_COEFFS-1)
      ELSE IF (NUM_COEFFS == 6) THEN
         WRITE(OUNT(2),2014) NUM_FITS, NUM_COEFFS
         WRITE(OUNT(2),2026) (B(I), I=0,NUM_COEFFS-1)
      ENDIF
      WRITE(OUNT(2),*)

      IF (RES_ABS > 0) THEN
         WRITE(OUNT(2),2031)
      ELSE
         WRITE(OUNT(2),2032)
      ENDIF

      DO II=1,NUM_FITS
         IF (DABS(WI_MAX) > 0.0) THEN
            WRITE(OUNT(2),2051) II, XI(II), YI(II), WI(II), WF(II), RES(II), POLY_PCT_ERR(II)
         ELSE
            WRITE(OUNT(2),2052) II, XI(II), YI(II), WI(II), WF(II), RES(II)
         ENDIF
      ENDDO

      WRITE(OUNT(2),*)
      WRITE(OUNT(2),2061)
      DO II=1,NUM_FITS
         WRITE(OUNT(2),2062) II, XO(II), YO(II), WO(II)
      ENDDO

! **********************************************************************************************************************************
 2000 FORMAT(' ******************************************************************************************************************',&
             '******************')

 2001 FORMAT(' Return code from polynomial fit: IFAULT = ',I3,/)

 2011 FORMAT(' Fitted quadratic surface for ',I2,' data points using a ',I2,'st order polynomial:',/,                              &
             ' -------------------------------------------------------------------------')

 2012 FORMAT(' Fitted quadratic surface for ',I2,' data points using a ',I2,'nd order polynomial:',/,                              &
             ' -------------------------------------------------------------------------')

 2013 FORMAT(' Fitted quadratic surface for ',I2,' data points using a ',I2,'rd order polynomial:',/,                              &
             ' -------------------------------------------------------------------------')

 2014 FORMAT(' Fitted quadratic surface for ',I2,' data points using a ',I2,'th order polynomial:',/,                              &
             ' -------------------------------------------------------------------------')

 2021 FORMAT('  Y = B0                                              ',//,                                                          &
             '      where: B0 = ',1ES14.6,/)

 2022 FORMAT('  Y = B0 + B1*X1                                      ',//,                                                          &
             '      where: B0 = ',1ES14.6,/,                                                                                &
             '             B1 = ',1ES14.6,/)

 2023 FORMAT('  Y = B0 + B1*X1 + B2*X2                              ',//,                                                          &
             '      where: B0 = ',1ES14.6,/,                                                                                &
             '             B1 = ',1ES14.6,/,                                                                                &
             '             B2 = ',1ES14.6,/)

 2024 FORMAT('  Y = B0 + B1*X1 + B2*X2 + B3*X1^2                    ',//,                                                          &
             '      where: B0 = ',1ES14.6,/,                                                                                &
             '             B1 = ',1ES14.6,/,                                                                                &
             '             B2 = ',1ES14.6,/,                                                                                &
             '             B3 = ',1ES14.6,/)

 2025 FORMAT('  Y = B0 + B1*X1 + B2*X2 + B3*X1^2 +B4*X1*Y1          ',//,                                                          &
             '      where: B0 = ',1ES14.6,/,                                                                                &
             '             B1 = ',1ES14.6,/,                                                                                &
             '             B2 = ',1ES14.6,/,                                                                                &
             '             B3 = ',1ES14.6,/,                                                                                &
             '             B4 = ',1ES14.6,/)

 2026 FORMAT('  Y = B0 + B1*X1 + B2*X2 +B3*X1^2 + B4*X1*Y1 + B5*Y2^2',//,                                                          &
             '      where: B0 = ',1ES14.6,/,                                                                                &
             '             B1 = ',1ES14.6,/,                                                                                &
             '             B2 = ',1ES14.6,/,                                                                                &
             '             B3 = ',1ES14.6,/,                                                                                &
             '             B4 = ',1ES14.6,/,                                                                                &
             '             B5 = ',1ES14.6,/)

 2031 FORMAT(' Comparison of actual data to fitted curve:',/,' -----------------------------------------',/,                       &
             '     I      XI(I)         YI(I)         WI(I)         WF(I)       WI - WF       Err')

 2032 FORMAT(' Comparison of actual data to fitted curve:',/,' -----------------------------------------',/,                       &
             '     I      XI(I)         YI(I)         WI(I)         WF(I)       WI - WF')

 2051 FORMAT(I6,5(1ES14.6),F8.3,'%')

 2052 FORMAT(I6,5(1ES14.6))

 2061 FORMAT(' Polynomial values at output data points:',/,' -----------------------------------------',/,                       &
             '     I      XO(I)         YO(I)         WO(I)')

 2062 FORMAT(I6,3(1ES14.6))

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

      END SUBROUTINE DEB_SURFACE_FIT

      END SUBROUTINE SURFACE_FIT


text:
program combsimple
 integer :: i,nbod1,nbod2,nsph1,nsph2,nstar1,nstar2,ioer
 real :: pmass,x(3),v(3)
 character*200 :: filename1,filename2,outfile

print*,'** combine to <simple> snapshot files **'
print*,'file 1?'
read*, filename1
print*,'file 2?'
read*,filename2
print*,'output file?'
read*,outfile

 open (1, file=filename1, form='unformatted')
 read (1,iostat=ioer) nbod1,nsph1,nstar1
   if(ioer) then
    rewind(1)
    read (1,iostat=ioer) nbod1
    nsph1=0
    nstar1=0
   endif
 print*, 'file1 has:',nbod1,nsph1,nstar1

 open (2, file=filename2, form='unformatted')
 read (2,iostat=ioer) nbod2,nsph2,nstar2
   if(ioer) then
    rewind(2)
    read (2,iostat=ioer) nbod2
    nsph2=0
    nstar2=0
   endif
 print*, 'file2 has:',nbod2,nsph2,nstar2

 open (3, file=outfile, form='unformatted')
 write(3), nbod1+nbod2,nsph1+nsph2,nstar1+nstar2

 do i=1,nsph1
  read(1) pmass,x,v
  write(3) pmass,x,v
 enddo
 do i=1,nsph2
  read(2) pmass,x,v
  write(3) pmass,x,v
 enddo
 do i=1,nbod1-nsph1-nstar1
  read(1) pmass,x,v
  write(3) pmass,x,v
 enddo
 do i=1,nbod2-nsph2-nstar2
  read(2) pmass,x,v
  write(3) pmass,x,v
 enddo
 do i=1,nstar1
  read(1) pmass,x,v
  write(3) pmass,x,v
 enddo
 do i=1,nstar2
  read(2) pmass,x,v
  write(3) pmass,x,v
 enddo

 close(1)
 close(2)
 close(3)

end program

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse