text:
! TITLE:
!   Fast sort (and select)
!
!
! PURPOSE:
!  This file contains the following modules and subroutines for quickly
!  sorting (either partially or fully) arrays in Fortran.
!
!   SWAP
!    SWAP_I64
!    SWAP_R64
!
!   FAST_SELECT
!    ARGSELECT_R64
!    ARGSELECT_I64
!
!   FAST_SORT
!    ARGSORT_R64
!    ARGPARTITION_R64
!    INSERTION_ARGSORT_R64
!
! Respectively, the purposes of the three modules are:
!
!   SWAP - Swap the value of two variables through a temporary.
!
!   FAST_SELECT - Do a fast "select" operation on an array of numbers,
!     resulting in the value at index K being as if the entire array
!     were sorted. All values before K will be less or equal, but not
!     sorted. All values after K will be greater or equal, but not
!     sorted. This operation is O(N).
!
!   FAST_SORT - Do a fast "sort" operation that uses methods
!     appropriate for the size of the array that will maximize speed.
!
!
! AUTHOR:
!   Thomas C.H. Lux (thomas.ch.lux@gmail.com)
!
!
! NOTES:
!   The algorithm for SELECT was inspired by the Floyd-Rivest method
!   of performing selection. The algorithm for SORT was inspired by
!   Data Structures and Algorithms Analysis by Dr. Clifford A. Shaffer
!   (Virginia Polytechnic Institute and State University).

! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
MODULE SWAP
  USE ISO_FORTRAN_ENV, ONLY: REAL64, INT64
  IMPLICIT NONE

CONTAINS
  SUBROUTINE SWAP_I64(V1, V2)
    INTEGER(KIND=INT64), INTENT(INOUT) :: V1, V2
    ! Local temp
    INTEGER(KIND=INT64) :: TEMP
    TEMP = V1
    V1 = V2
    V2 = TEMP
  END SUBROUTINE SWAP_I64

  SUBROUTINE SWAP_R64(V1, V2)
    REAL(KIND=REAL64), INTENT(INOUT) :: V1, V2
    ! Local temp
    REAL(KIND=REAL64) :: TEMP
    TEMP = V1
    V1 = V2
    V2 = TEMP
  END SUBROUTINE SWAP_R64
END MODULE SWAP
! ____________________________________________________________________

! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
MODULE FAST_SELECT
  USE ISO_FORTRAN_ENV, ONLY: REAL64, INT64
  USE SWAP, ONLY: SWAP_I64, SWAP_R64
  IMPLICIT NONE

CONTAINS
  ! ------------------------------------------------------------------
  !                       FastSelect method
  !
  ! Given VALUES list of numbers, rearrange the elements of VALUES
  ! such that the element at index K has rank K (holds its same
  ! location as if all of VALUES were sorted). Symmetrically rearrange
  ! array INDICES to keep track of prior indices.
  !
  ! This algorithm uses the same conceptual approach as Floyd-Rivest,
  ! but instead of standard-deviation based selection of bounds for
  ! recursion, a rank-based method is used to pick the subset of
  ! values that is searched. This simplifies the code and improves
  ! interpretability, while achieving the same tunable performance.
  !
  ! Arguments:
  !
  !   VALUES   --  A 1D array of real numbers.
  !   INDICES  --  A 1D array of original indices for elements of VALUES.
  !   K        --  A positive integer for the rank index about which
  !                VALUES should be rearranged.
  ! Optional:
  !
  !   DIVISOR  --  A positive integer >= 2 that represents the
  !                division factor used for large VALUES arrays.
  !   MAX_SIZE --  An integer >= DIVISOR that represents the largest
  !                sized VALUES for which the worst-case pivot value
  !                selection is tolerable. A worst-case pivot causes
  !                O( SIZE(VALUES)^2 ) runtime. This value should be
  !                determined heuristically based on compute hardware.
  !
  ! Output:
  !
  !   The elements of the array VALUES are rearranged such that the
  !   element at position VALUES(K) is in the same location it would
  !   be if all of VALUES were in sorted order. Also known as,
  !   VALUES(K) has rank K.
  !
  RECURSIVE SUBROUTINE ARGSELECT_R64(VALUES, INDICES, K, DIVISOR, MAX_SIZE)
    ! Arguments
    REAL(KIND=REAL64),   INTENT(INOUT), DIMENSION(:) :: VALUES
    INTEGER(KIND=INT64), INTENT(INOUT), DIMENSION(:) :: INDICES
    INTEGER(KIND=INT64), INTENT(IN)                  :: K
    INTEGER(KIND=INT64), INTENT(IN), OPTIONAL        :: DIVISOR, MAX_SIZE
    ! Locals
    INTEGER(KIND=INT64) :: LEFT, RIGHT, L, R, MS, D
    REAL(KIND=REAL64) :: P
    ! Initialize the divisor (for making subsets).
    IF (PRESENT(DIVISOR)) THEN ; D = DIVISOR
    ELSE IF (SIZE(VALUES) .GE. 2**23) THEN ; D = 2**5
    ELSE IF (SIZE(VALUES) .GE. 2**20) THEN ; D = 2**3
    ELSE                                   ; D = 2**2
    END IF
    ! Initialize the max size (before subsets are created).
    IF (PRESENT(MAX_SIZE)) THEN ; MS = MAX_SIZE
    ELSE                        ; MS = 2**10
    END IF
    ! Initialize LEFT and RIGHT to be the entire array.
    LEFT = 1
    RIGHT = SIZE(VALUES)
    ! Loop until done finding the K-th element.
    DO WHILE (LEFT .LT. RIGHT)
       ! Use SELECT recursively to improve the quality of the
       ! selected pivot value for large arrays.
       IF (RIGHT - LEFT .GT. MS) THEN
          ! Compute how many elements should be left and right of K
          ! to maintain the same percentile in a subset.
          L = K - K / D
          R = L + (SIZE(VALUES) / D)
          ! Perform fast select on an array a fraction of the size about K.
          CALL ARGSELECT_R64(VALUES(L:R), INDICES(L:R), K - L + 1, DIVISOR, MAX_SIZE)
       END IF
       ! Pick a partition element at position K.
       P = VALUES(K)
       L = LEFT
       R = RIGHT
       ! Move the partition element to the front of the list.
       CALL SWAP_R64(VALUES(LEFT), VALUES(K))
       CALL SWAP_I64(INDICES(LEFT), INDICES(K))
       ! Pre-swap the left and right elements (temporarily putting a
       ! larger element on the left) before starting the partition loop.
       IF (VALUES(RIGHT) .GT. P) THEN
          CALL SWAP_R64(VALUES(LEFT), VALUES(RIGHT))
          CALL SWAP_I64(INDICES(LEFT), INDICES(RIGHT))
       END IF
       ! Now partition the elements about the pivot value "P".
       DO WHILE (L .LT. R)
          CALL SWAP_R64(VALUES(L), VALUES(R))
          CALL SWAP_I64(INDICES(L), INDICES(R))
          L = L + 1
          R = R - 1
          DO WHILE (VALUES(L) .LT. P) ; L = L + 1 ; END DO
          DO WHILE (VALUES(R) .GT. P) ; R = R - 1 ; END DO
       END DO
       ! Place the pivot element back into its appropriate place.
       IF (VALUES(LEFT) .EQ. P) THEN
          CALL SWAP_R64(VALUES(LEFT), VALUES(R))
          CALL SWAP_I64(INDICES(LEFT), INDICES(R))
       ELSE
          R = R + 1
          CALL SWAP_R64(VALUES(R), VALUES(RIGHT))
          CALL SWAP_I64(INDICES(R), INDICES(RIGHT))
       END IF
       ! adjust left and right towards the boundaries of the subset
       ! containing the (k - left + 1)th smallest element
       IF (R .LE. K) LEFT = R + 1
       IF (K .LE. R) RIGHT = R - 1
    END DO
  END SUBROUTINE ARGSELECT_R64

  ! ------------------------------------------------------------------

  RECURSIVE SUBROUTINE ARGSELECT_I64(VALUES, INDICES, K, DIVISOR, MAX_SIZE)
    ! Arguments
    INTEGER(KIND=INT64), INTENT(INOUT), DIMENSION(:) :: VALUES
    INTEGER(KIND=INT64), INTENT(INOUT), DIMENSION(:) :: INDICES
    INTEGER(KIND=INT64), INTENT(IN)                  :: K
    INTEGER(KIND=INT64), INTENT(IN), OPTIONAL        :: DIVISOR, MAX_SIZE
    ! Locals
    INTEGER(KIND=INT64) :: LEFT, RIGHT, L, R, MS, D, P
    ! Initialize the divisor (for making subsets).
    IF (PRESENT(DIVISOR)) THEN ; D = DIVISOR
    ELSE IF (SIZE(VALUES) .GE. 2**23) THEN ; D = 2**5
    ELSE IF (SIZE(VALUES) .GE. 2**20) THEN ; D = 2**3
    ELSE                                   ; D = 2**2
    END IF
    ! Initialize the max size (before subsets are created).
    IF (PRESENT(MAX_SIZE)) THEN ; MS = MAX_SIZE
    ELSE                        ; MS = 2**10
    END IF
    ! Initialize LEFT and RIGHT to be the entire array.
    LEFT = 1
    RIGHT = SIZE(VALUES)
    ! Loop until done finding the K-th element.
    DO WHILE (LEFT .LT. RIGHT)
       ! Use SELECT recursively to improve the quality of the
       ! selected pivot value for large arrays.
       IF (RIGHT - LEFT .GT. MS) THEN
          ! Compute how many elements should be left and right of K
          ! to maintain the same percentile in a subset.
          L = K - K / D
          R = L + (SIZE(VALUES) / D)
          ! Perform fast select on an array a fraction of the size about K.
          CALL ARGSELECT_I64(VALUES(L:R), INDICES(L:R), K - L + 1, DIVISOR, MAX_SIZE)
       END IF
       ! Pick a partition element at position K.
       P = VALUES(K)
       L = LEFT
       R = RIGHT
       ! Move the partition element to the front of the list.
       CALL SWAP_I64(VALUES(LEFT), VALUES(K))
       CALL SWAP_I64(INDICES(LEFT), INDICES(K))
       ! Pre-swap the left and right elements (temporarily putting a
       ! larger element on the left) before starting the partition loop.
       IF (VALUES(RIGHT) .GT. P) THEN
          CALL SWAP_I64(VALUES(LEFT), VALUES(RIGHT))
          CALL SWAP_I64(INDICES(LEFT), INDICES(RIGHT))
       END IF
       ! Now partition the elements about the pivot value "T".
       DO WHILE (L .LT. R)
          CALL SWAP_I64(VALUES(L), VALUES(R))
          CALL SWAP_I64(INDICES(L), INDICES(R))
          L = L + 1
          R = R - 1
          DO WHILE (VALUES(L) .LT. P) ; L = L + 1 ; END DO
          DO WHILE (VALUES(R) .GT. P) ; R = R - 1 ; END DO
       END DO
       ! Place the pivot element back into its appropriate place.
       IF (VALUES(LEFT) .EQ. P) THEN
          CALL SWAP_I64(VALUES(LEFT), VALUES(R))
          CALL SWAP_I64(INDICES(LEFT), INDICES(R))
       ELSE
          R = R + 1
          CALL SWAP_I64(VALUES(R), VALUES(RIGHT))
          CALL SWAP_I64(INDICES(R), INDICES(RIGHT))
       END IF
       ! adjust left and right towards the boundaries of the subset
       ! containing the (k - left + 1)th smallest element
       IF (R .LE. K) LEFT = R + 1
       IF (K .LE. R) RIGHT = R - 1
    END DO
  END SUBROUTINE ARGSELECT_I64
END MODULE FAST_SELECT
! ____________________________________________________________________

! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
MODULE FAST_SORT
  USE ISO_FORTRAN_ENV, ONLY: REAL64, INT64
  USE SWAP, ONLY: SWAP_I64, SWAP_R64
  IMPLICIT NONE

CONTAINS
  ! ------------------------------------------------------------------
  !                         FastSort
  !
  ! This routine uses a combination of QuickSort (with modestly
  ! intelligent pivot selection) and Insertion Sort (for small arrays)
  ! to achieve very fast average case sort times for both random and
  ! partially sorted data. The pivot is selected for QuickSort as the
  ! median of the first, middle, and last values in the array.
  !
  ! Arguments:
  !
  !   VALUES   --  A 1D array of real numbers.
  !   INDICES  --  A 1D array of original indices for elements of VALUES.
  !
  ! Optional:
  !
  !   MIN_SIZE --  An positive integer that represents the largest
  !                sized VALUES for which a partition about a pivot
  !                is used to reduce the size of a an unsorted array.
  !                Any size less than this will result in the use of
  !                INSERTION_ARGSORT instead of ARGPARTITION.
  !
  ! Output:
  !
  !   The elements of the array VALUES are sorted and all elements of
  !   INDICES are sorted symmetrically (given INDICES = 1, ...,
  !   SIZE(VALUES) beforehand, final INDICES will show original index
  !   of each element of VALUES before the sort operation).
  !
  RECURSIVE SUBROUTINE ARGSORT_R64(VALUES, INDICES, MIN_SIZE)
    REAL(KIND=REAL64),   INTENT(INOUT), DIMENSION(:)            :: VALUES
    INTEGER(KIND=INT64), INTENT(INOUT), DIMENSION(SIZE(VALUES)) :: INDICES
    INTEGER(KIND=INT64), INTENT(IN), OPTIONAL                   :: MIN_SIZE
    ! Local variables
    INTEGER(KIND=INT64) :: I, MS
    IF (PRESENT(MIN_SIZE)) THEN ; MS = MIN_SIZE
    ELSE                        ; MS = 2**6
    END IF
    ! Base case, return.
    IF (SIZE(VALUES) .LT. MS) THEN
       CALL INSERTION_ARGSORT_R64(VALUES, INDICES)
       ! Call this function recursively after pivoting about the median.
    ELSE
       ! ---------------------------------------------------------------
       ! If you are having slow runtime with the selection of pivot values
       ! provided by ARGPARTITION, then consider using ARGSELECT instead.
       I = ARGPARTITION_R64(VALUES, INDICES)
       ! ---------------------------------------------------------------
       ! I = SIZE(VALUES) / 2
       ! CALL ARGSELECT_R64(VALUES, INDICES, I)
       ! ! Requires 'USE FAST_SELECT' at top of subroutine or module.
       ! ---------------------------------------------------------------
       CALL ARGSORT_R64(VALUES(:I-1), INDICES(:I-1), MS)
       CALL ARGSORT_R64(VALUES(I+1:), INDICES(I+1:), MS)
    END IF
  END SUBROUTINE ARGSORT_R64

  ! This function efficiently partitions values based on the median
  ! of the first, middle, and last elements of the VALUES array. This
  ! function returns the index of the pivot.
  FUNCTION ARGPARTITION_R64(VALUES, INDICES) RESULT(LEFT)
    REAL(KIND=REAL64),   INTENT(INOUT), DIMENSION(:)            :: VALUES
    INTEGER(KIND=INT64), INTENT(INOUT), DIMENSION(SIZE(VALUES)) :: INDICES
    INTEGER(KIND=INT64) :: LEFT, MID, RIGHT
    REAL(KIND=REAL64)   :: PIVOT
    ! Use the median of the first, middle, and last element as the
    ! pivot. Place the pivot at the end of the array.
    MID = (1 + SIZE(VALUES)) / 2
    ! Swap the first and last elements (if the last is smaller).
    IF (VALUES(SIZE(VALUES)) < VALUES(1)) THEN
       CALL SWAP_R64(VALUES(1),  VALUES(SIZE(VALUES)))
       CALL SWAP_I64(INDICES(1), INDICES(SIZE(VALUES)))
    END IF
    ! Swap the middle and first elements (if the middle is smaller).
    IF (VALUES(MID) < VALUES(SIZE(VALUES))) THEN
       CALL SWAP_R64(VALUES(MID),  VALUES(SIZE(VALUES)))
       CALL SWAP_I64(INDICES(MID), INDICES(SIZE(VALUES)))
       ! Swap the last and first elements (if the last is smaller).
       IF (VALUES(SIZE(VALUES)) < VALUES(1)) THEN
          CALL SWAP_R64(VALUES(1),  VALUES(SIZE(VALUES)))
          CALL SWAP_I64(INDICES(1), INDICES(SIZE(VALUES)))
       END IF
    END IF
    ! Set the pivot, LEFT index and RIGHT index (skip the smallest,
    ! which is in location 1, and the pivot at the end).
    PIVOT = VALUES(SIZE(VALUES))
    LEFT  = 2
    RIGHT = SIZE(VALUES) - 1
    ! Partition all elements to the left and right side of the pivot
    ! (left if they are smaller, right if they are bigger).
    DO WHILE (LEFT < RIGHT)
       ! Loop left until we find a value that is greater or equal to pivot.
       DO WHILE (VALUES(LEFT) < PIVOT)
          LEFT = LEFT + 1
       END DO
       ! Loop right until we find a value that is less or equal to pivot (or LEFT).
       DO WHILE (RIGHT .NE. LEFT)
          IF (VALUES(RIGHT) .LT. PIVOT) EXIT
          RIGHT = RIGHT - 1
       END DO
       ! Now we know that [VALUES(RIGHT) < PIVOT < VALUES(LEFT)], so swap them.
       CALL SWAP_R64(VALUES(LEFT),  VALUES(RIGHT))
       CALL SWAP_I64(INDICES(LEFT), INDICES(RIGHT))
    END DO
    ! The last swap was done even though LEFT == RIGHT, we need to undo.
    CALL SWAP_R64(VALUES(LEFT),  VALUES(RIGHT))
    CALL SWAP_I64(INDICES(LEFT), INDICES(RIGHT))
    ! Finally, we put the pivot back into its proper location.
    CALL SWAP_R64(VALUES(LEFT),  VALUES(SIZE(VALUES)))
    CALL SWAP_I64(INDICES(LEFT), INDICES(SIZE(VALUES)))
  END FUNCTION ARGPARTITION_R64

  ! Insertion sort (best for small lists).
  SUBROUTINE INSERTION_ARGSORT_R64(VALUES, INDICES)
    REAL(KIND=REAL64),   INTENT(INOUT), DIMENSION(:)            :: VALUES
    INTEGER(KIND=INT64), INTENT(INOUT), DIMENSION(SIZE(VALUES)) :: INDICES
    ! Local variables.
    REAL(KIND=REAL64)   :: TEMP_VAL
    INTEGER(KIND=INT64) :: I, BEFORE, AFTER, TEMP_IND
    ! Return for the base case.
    IF (SIZE(VALUES) .LE. 1) RETURN
    ! Put the smallest value at the front of the list.
    I = MINLOC(VALUES,1)
    CALL SWAP_R64(VALUES(1),  VALUES(I))
    CALL SWAP_I64(INDICES(1), INDICES(I))
    ! Insertion sort the rest of the array.
    DO I = 3, SIZE(VALUES)
       TEMP_VAL = VALUES(I)
       TEMP_IND = INDICES(I)
       ! Search backwards in the list,
       BEFORE = I - 1
       AFTER  = I
       DO WHILE (TEMP_VAL .LT. VALUES(BEFORE))
          VALUES(AFTER)  = VALUES(BEFORE)
          INDICES(AFTER) = INDICES(BEFORE)
          BEFORE = BEFORE - 1
          AFTER  = AFTER - 1
       END DO
       ! Put the value into its place (where it is greater than the
       ! element before it, but less than all values after it).
       VALUES(AFTER)  = TEMP_VAL
       INDICES(AFTER) = TEMP_IND
    END DO
  END SUBROUTINE INSERTION_ARGSORT_R64
  ! ------------------------------------------------------------------

END MODULE FAST_SORT

max_stars_repo_name: thchang/SortVectors

max_issues_repo_name: thchang/SortVectors

max_forks_repo_name: thchang/SortVectors


text:
C***********************************************************************
C    Module:  mousexy.f
C
C    Copyright (C) 2012 Harold Youngren, Mark Drela
C
C    This program is free software; you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation; either version 2 of the License, or
C    (at your option) any later version.
C
C    This program is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C    You should have received a copy of the GNU General Public License
C    along with this program; if not, write to the Free Software
C    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
C
C    Report problems to:    guppy@maine.com
C                        or drela@mit.edu
C***********************************************************************

      program mousexy
c---------------------------------------------------------------
c     Interactive mouse and continuous input test program.
c
c     Lets user specify points in a spline, then allows the user to
c     drag spline points around to change the curve
c---------------------------------------------------------------
c
      parameter (ipx=100)
c
      character*1 cin
c
      common /Rpoints/
     &       ch,
     &       x1(ipx), y1(ipx), s1(ipx), xp1(ipx), yp1(ipx),
     &       x2(ipx), y2(ipx), s2(ipx), xp2(ipx), yp2(ipx)
      common /Ipoints/
     &       iplt, n1, n2
      logical lok
c
      ch = 0.13
      n1 = 0
      n2 = 0
      iplt = 0
c
 1000 format(a)
C
C---Initialize the plot package before we get into plotting...
      CALL PLINITIALIZE

C---put up plot window and refresh
      call pltall('(terminate with q or double point)')

c
ccc      call NEWFACTOR(6.0)
c
 5    write(*,1050)
 1050 format(/' Enter points for spline (q or ret to end)')

C--- Get initial set of points from user
      call NEWCOLORNAME('green')
      call NEWPEN(1)
C
      xlast = -999.
      ylast = -999.
      ii = 0
      do j = 1, 100
        call GETCURSORXY(xx,yy,ikey)
        cin = char(ikey)
        if(cin.EQ.'q' .OR. cin.EQ.'Q') go to 100
C--- check for doubled point to end input
          if( (xx-xlast).EQ.0.0 .AND.
     &        (yy-ylast).EQ.0.0 ) go to 100
        ii = ii + 1
        x1(ii) = xx
        y1(ii) = yy
         if(ii.EQ.1) then
          call plot(xx,yy,3)
          CALL PLSYMB(999.,999.,0.4*CH,2,0.0,0)
         else
          call plot(xlast,ylast,3)
          call plot(xx,yy,2)
          CALL PLSYMB(999.,999.,0.4*CH,2,0.0,0)
          call PLFLUSH
        endif
        xlast = xx
        ylast = yy
      end do
C
 100  n1 = ii
C
      call plend
C--- initialize modified point arrays
       do j = 1, n1
         x2(j) = x1(j)
         y2(j) = y1(j)
       end do
       n2 = n1

      call pltall('Select point to modify ("i" to reset)')
c
      write(*,*)
      write(*,*) 'Move spline points...'
C
 200  call GETCURSORXY(xx,yy,ikey)
cc      call NEWCOLORNAME('green')
cc      call plot(xx,yy,3)
cc      CALL PLSYMB(999.,999.,0.4*CH,2,0.0,0)
      cin = char(ikey)
      write(*,*) 'ikey cin ',ikey,cin

      if(cin.EQ.'q' .OR. cin.EQ.'Q') go to 400
C--- reset modified point arrays
      if(cin.EQ.'i' .OR. cin.EQ.'I') then
       write(*,*) 're-initializing points...'
       do j = 1, n1
         x2(j) = x1(j)
         y2(j) = y1(j)
       end do
       n2 = n1
       call pltall('Select point to modify ("i" to reset)')
       go to 200
      endif

C--- find closest point to cursor
      jp = 0
      dsqmin = 999.
      do j = 1, n1
        dsq = (xx-x2(j))**2 + (yy-y2(j))**2
        if(dsq.LT.dsqmin) then
          dsqmin = dsq
          jp = j
        endif
      end do
c
      write(*,*) 'Cursor at x,y ',xx,yy
      write(*,*) 'Point ',jp,' selected at x,y ',x1(jp),y1(jp)
C--- now read cursor and move point
 300  call GETCURSORXYC(xx,yy,ibtn)
      write(*,*) xx,yy,ibtn
      if(ibtn.eq.0) go to 200
      x2(jp) = xx
      y2(jp) = yy
      n2 = n1
      call pltall('Move point...("i" to reset)')
      go to 300
c
 400  continue
ccc      pause
      stop
      end

      subroutine pltall(msg)
      character*(*) msg
      parameter (ipx=100)
c
      common /Rpoints/
     &       ch,
     &       x1(ipx), y1(ipx), s1(ipx), xp1(ipx), yp1(ipx),
     &       x2(ipx), y2(ipx), s2(ipx), xp2(ipx), yp2(ipx)
      common /Ipoints/
     &       iplt, n1, n2
C
      call PLOPEN(0.8,0,1)
      CALL PLCHAR(0.5,0.5,1.2*CH,'TEST FOR MOUSE READ',0.0,-1)
      CALL PLCHAR(0.3,0.3,1.2*CH,msg,0.0,-1)
ccc      call PLOT(5.5, 4.25, -3)

      if(n1.GT.1) then
C--- spline input curve and plot
      call scalc(x1,y1,s1,n1)
      call splind(x1,xp1,s1,n1,-999.,-999.)
      call splind(y1,yp1,s1,n1,-999.,-999.)
C
      call NEWCOLORNAME('blue')
      call plot(x1(1),y1(1),3)
      CALL PLSYMB(999.,999.,0.4*CH,1,0.0,0)
      ninter = 10
      dfn = 1.0/float(ninter)
      do j = 1, n1-1
        ds = s1(j+1) - s1(j)
        do n = 1, ninter
          ss = s1(j) + ds*float(n)*dfn
          xx = seval(ss,x1,xp1,s1,n1)
          yy = seval(ss,y1,yp1,s1,n1)
          call plot(xx,yy,2)
        end do
        CALL PLSYMB(999.,999.,0.4*CH,1,0.0,0)
      end do
      endif

      if(n2.GT.1) then
C--- spline modified curve and plot
      call scalc(x2,y2,s2,n2)
      call splind(x2,xp2,s2,n2,-999.,-999.)
      call splind(y2,yp2,s2,n2,-999.,-999.)
C
      call NEWCOLORNAME('red')
      call plot(x2(1),y2(1),3)
      CALL PLSYMB(999.,999.,0.4*CH,1,0.0,0)
      ninter = 10
      dfn = 1.0/float(ninter)
      do j = 1, n2-1
        ds = s2(j+1) - s2(j)
        do n = 1, ninter
          ss = s2(j) + ds*float(n)*dfn
          xx = seval(ss,x2,xp2,s2,n2)
          yy = seval(ss,y2,yp2,s2,n2)
          call plot(xx,yy,2)
        end do
        CALL PLSYMB(999.,999.,0.4*CH,1,0.0,0)
      end do
      endif

      call PLFLUSH
      return
      end

C***********************************************************************
C    Module:  spline.f
C
C    Copyright (C) 2000 Mark Drela
C
C    This program is free software; you can redistribute it and/or modify
C    it under the terms of the GNU General Public License as published by
C    the Free Software Foundation; either version 2 of the License, or
C    (at your option) any later version.
C
C    This program is distributed in the hope that it will be useful,
C    but WITHOUT ANY WARRANTY; without even the implied warranty of
C    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C    GNU General Public License for more details.
C
C    You should have received a copy of the GNU General Public License
C    along with this program; if not, write to the Free Software
C    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
C***********************************************************************

      SUBROUTINE SPLINE(X,XS,S,N)
      DIMENSION X(N),XS(N),S(N)
      PARAMETER (NMAX=1000)
      DIMENSION A(NMAX),B(NMAX),C(NMAX)
C-------------------------------------------------------
C     Calculates spline coefficients for X(S).          |
C     Zero 2nd derivative end conditions are used.      |
C     To evaluate the spline at some value of S,        |
C     use SEVAL and/or DEVAL.                           |
C                                                       |
C     S        independent variable array (input)       |
C     X        dependent variable array   (input)       |
C     XS       dX/dS array                (calculated)  |
C     N        number of points           (input)       |
C                                                       |
C-------------------------------------------------------
      IF(N.GT.NMAX) STOP 'SPLINE: array overflow, increase NMAX'
C
      DO 1 I=2, N-1
        DSM = S(I) - S(I-1)
        DSP = S(I+1) - S(I)
        B(I) = DSP
        A(I) = 2.0*(DSM+DSP)
        C(I) = DSM
        XS(I) = 3.0*((X(I+1)-X(I))*DSM/DSP + (X(I)-X(I-1))*DSP/DSM)
    1 CONTINUE
C
C---- set zero second derivative end conditions
      A(1) = 2.0
      C(1) = 1.0
      XS(1) = 3.0*(X(2)-X(1)) / (S(2)-S(1))
      B(N) = 1.0
      A(N) = 2.0
      XS(N) = 3.0*(X(N)-X(N-1)) / (S(N)-S(N-1))
C
C---- solve for derivative array XS
      CALL TRISOL(A,B,C,XS,N)
C
      RETURN
      END ! SPLINE

      SUBROUTINE SPLIND(X,XS,S,N,XS1,XS2)
      DIMENSION X(N),XS(N),S(N)
      PARAMETER (NMAX=1000)
      DIMENSION  A(NMAX),B(NMAX),C(NMAX)
C-------------------------------------------------------
C     Calculates spline coefficients for X(S).          |
C     Specified 1st derivative and/or usual zero 2nd    |
C     derivative end conditions are used.               |
C     To evaluate the spline at some value of S,        |
C     use SEVAL and/or DEVAL.                           |
C                                                       |
C     S        independent variable array (input)       |
C     X        dependent variable array   (input)       |
C     XS       dX/dS array                (calculated)  |
C     N        number of points           (input)       |
C     XS1,XS2  endpoint derivatives       (input)       |
C              If = 999.0, then usual zero second       |
C              derivative end condition(s) are used     |
C              If = -999.0, then zero third             |
C              derivative end condition(s) are used     |
C                                                       |
C-------------------------------------------------------
      IF(N.GT.NMAX) STOP 'SPLIND: array overflow, increase NMAX'
C
      DO 1 I=2, N-1
        DSM = S(I) - S(I-1)
        DSP = S(I+1) - S(I)
        B(I) = DSP
        A(I) = 2.0*(DSM+DSP)
        C(I) = DSM
        XS(I) = 3.0*((X(I+1)-X(I))*DSM/DSP + (X(I)-X(I-1))*DSP/DSM)
    1 CONTINUE
C
      IF(XS1.EQ.999.0) THEN
C----- set zero second derivative end condition
       A(1) = 2.0
       C(1) = 1.0
       XS(1) = 3.0*(X(2)-X(1)) / (S(2)-S(1))
      ELSE IF(XS1.EQ.-999.0) THEN
C----- set zero third derivative end condition
       A(1) = 1.0
       C(1) = 1.0
       XS(1) = 2.0*(X(2)-X(1)) / (S(2)-S(1))
      ELSE
C----- set specified first derivative end condition
       A(1) = 1.0
       C(1) = 0.
       XS(1) = XS1
      ENDIF
C
      IF(XS2.EQ.999.0) THEN
       B(N) = 1.0
       A(N) = 2.0
       XS(N) = 3.0*(X(N)-X(N-1)) / (S(N)-S(N-1))
      ELSE IF(XS2.EQ.-999.0) THEN
       B(N) = 1.0
       A(N) = 1.0
       XS(N) = 2.0*(X(N)-X(N-1)) / (S(N)-S(N-1))
      ELSE
       A(N) = 1.0
       B(N) = 0.
       XS(N) = XS2
      ENDIF
C
      IF(N.EQ.2 .AND. XS1.EQ.-999.0 .AND. XS2.EQ.-999.0) THEN
       B(N) = 1.0
       A(N) = 2.0
       XS(N) = 3.0*(X(N)-X(N-1)) / (S(N)-S(N-1))
      ENDIF
C
C---- solve for derivative array XS
      CALL TRISOL(A,B,C,XS,N)
C
      RETURN
      END ! SPLIND

      SUBROUTINE SPLINA(X,XS,S,N)
      IMPLICIT REAL (A-H,O-Z)
      DIMENSION X(N),XS(N),S(N)
      LOGICAL LEND
C-------------------------------------------------------
C     Calculates spline coefficients for X(S).          |
C     A simple averaging of adjacent segment slopes     |
C     is used to achieve non-oscillatory curve          |
C     End conditions are set by end segment slope       |
C     To evaluate the spline at some value of S,        |
C     use SEVAL and/or DEVAL.                           |
C                                                       |
C     S        independent variable array (input)       |
C     X        dependent variable array   (input)       |
C     XS       dX/dS array                (calculated)  |
C     N        number of points           (input)       |
C                                                       |
C-------------------------------------------------------
C
      LEND = .TRUE.
      DO 1 I=1, N-1
        DS = S(I+1)-S(I)
        IF (DS.EQ.0.) THEN
          XS(I) = XS1
          LEND = .TRUE.
         ELSE
          DX = X(I+1)-X(I)
          XS2 = DX / DS
          IF (LEND) THEN
            XS(I) = XS2
            LEND = .FALSE.
           ELSE
            XS(I) = 0.5*(XS1 + XS2)
          ENDIF
        ENDIF
        XS1 = XS2
    1 CONTINUE
      XS(N) = XS1
C
      RETURN
      END ! SPLINA

      SUBROUTINE TRISOL(A,B,C,D,KK)
      DIMENSION A(KK),B(KK),C(KK),D(KK)
C-----------------------------------------
C     Solves KK long, tri-diagonal system |
C                                         |
C             A C          D              |
C             B A C        D              |
C               B A .      .              |
C                 . . C    .              |
C                   B A    D              |
C                                         |
C     The righthand side D is replaced by |
C     the solution.  A, C are destroyed.  |
C-----------------------------------------
C
      DO 1 K=2, KK
        KM = K-1
        C(KM) = C(KM) / A(KM)
        D(KM) = D(KM) / A(KM)
        A(K) = A(K) - B(K)*C(KM)
        D(K) = D(K) - B(K)*D(KM)
    1 CONTINUE
C
      D(KK) = D(KK)/A(KK)
C
      DO 2 K=KK-1, 1, -1
        D(K) = D(K) - C(K)*D(K+1)
    2 CONTINUE
C
      RETURN
      END ! TRISOL

      FUNCTION SEVAL(SS,X,XS,S,N)
      DIMENSION X(N), XS(N), S(N)
C--------------------------------------------------
C     Calculates X(SS)                             |
C     XS array must have been calculated by SPLINE |
C--------------------------------------------------
      ILOW = 1
      I = N
C
   10 IF(I-ILOW .LE. 1) GO TO 11
C
      IMID = (I+ILOW)/2
      IF(SS .LT. S(IMID)) THEN
       I = IMID
      ELSE
       ILOW = IMID
      ENDIF
      GO TO 10
C
   11 DS = S(I) - S(I-1)
      T = (SS - S(I-1)) / DS
      CX1 = DS*XS(I-1) - X(I) + X(I-1)
      CX2 = DS*XS(I)   - X(I) + X(I-1)
      SEVAL = T*X(I) + (1.0-T)*X(I-1) + (T-T*T)*((1.0-T)*CX1 - T*CX2)
      RETURN
      END ! SEVAL

      FUNCTION DEVAL(SS,X,XS,S,N)
      DIMENSION X(N), XS(N), S(N)
C--------------------------------------------------
C     Calculates dX/dS(SS)                         |
C     XS array must have been calculated by SPLINE |
C--------------------------------------------------
      ILOW = 1
      I = N
C
   10 IF(I-ILOW .LE. 1) GO TO 11
C
      IMID = (I+ILOW)/2
      IF(SS .LT. S(IMID)) THEN
       I = IMID
      ELSE
       ILOW = IMID
      ENDIF
      GO TO 10
C
   11 DS = S(I) - S(I-1)
      T = (SS - S(I-1)) / DS
      CX1 = DS*XS(I-1) - X(I) + X(I-1)
      CX2 = DS*XS(I)   - X(I) + X(I-1)
      DEVAL = X(I) - X(I-1) + (1.-4.0*T+3.0*T*T)*CX1 + T*(3.0*T-2.)*CX2
      DEVAL = DEVAL/DS
      RETURN
      END ! DEVAL

      FUNCTION D2VAL(SS,X,XS,S,N)
      DIMENSION X(N), XS(N), S(N)
C--------------------------------------------------
C     Calculates d2X/dS2(SS)                       |
C     XS array must have been calculated by SPLINE |
C--------------------------------------------------
      ILOW = 1
      I = N
C
   10 IF(I-ILOW .LE. 1) GO TO 11
C
      IMID = (I+ILOW)/2
      IF(SS .LT. S(IMID)) THEN
       I = IMID
      ELSE
       ILOW = IMID
      ENDIF
      GO TO 10
C
   11 DS = S(I) - S(I-1)
      T = (SS - S(I-1)) / DS
      CX1 = DS*XS(I-1) - X(I) + X(I-1)
      CX2 = DS*XS(I)   - X(I) + X(I-1)
      D2VAL = (6.*T-4.)*CX1 + (6.*T-2.0)*CX2
      D2VAL = D2VAL/DS**2
      RETURN
      END ! D2VAL

      FUNCTION CURV(SS,X,XS,Y,YS,S,N)
      DIMENSION X(N), XS(N), Y(N), YS(N), S(N)
C-----------------------------------------------
C     Calculates curvature of splined 2-D curve |
C     at S = SS                                 |
C                                               |
C     S        arc length array of curve        |
C     X, Y     coordinate arrays of curve       |
C     XS,YS    derivative arrays                |
C              (calculated earlier by SPLINE)   |
C-----------------------------------------------
C
      ILOW = 1
      I = N
C
   10 IF(I-ILOW .LE. 1) GO TO 11
C
      IMID = (I+ILOW)/2
      IF(SS .LT. S(IMID)) THEN
       I = IMID
      ELSE
       ILOW = IMID
      ENDIF
      GO TO 10
C
   11 DS = S(I) - S(I-1)
      T = (SS - S(I-1)) / DS
C
      CX1 = DS*XS(I-1) - X(I) + X(I-1)
      CX2 = DS*XS(I)   - X(I) + X(I-1)
      XD = X(I) - X(I-1) + (1.0-4.0*T+3.0*T*T)*CX1 + T*(3.0*T-2.0)*CX2
      XDD = (6.0*T-4.0)*CX1 + (6.0*T-2.0)*CX2
C
      CY1 = DS*YS(I-1) - Y(I) + Y(I-1)
      CY2 = DS*YS(I)   - Y(I) + Y(I-1)
      YD = Y(I) - Y(I-1) + (1.0-4.0*T+3.0*T*T)*CY1 + T*(3.0*T-2.0)*CY2
      YDD = (6.0*T-4.0)*CY1 + (6.0*T-2.0)*CY2
C
      SD = SQRT(XD*XD + YD*YD)
      SD = MAX(SD,0.001*DS)
C
      CURV = (XD*YDD - YD*XDD) / SD**3
C
      RETURN
      END ! CURV

      FUNCTION CURVS(SS,X,XS,Y,YS,S,N)
      DIMENSION X(N), XS(N), Y(N), YS(N), S(N)
C-----------------------------------------------
C     Calculates curvature derivative of        |
C     splined 2-D curve at S = SS               |
C                                               |
C     S        arc length array of curve        |
C     X, Y     coordinate arrays of curve       |
C     XS,YS    derivative arrays                |
C              (calculated earlier by SPLINE)   |
C-----------------------------------------------
C
      ILOW = 1
      I = N
C
   10 IF(I-ILOW .LE. 1) GO TO 11
C
      IMID = (I+ILOW)/2
      IF(SS .LT. S(IMID)) THEN
       I = IMID
      ELSE
       ILOW = IMID
      ENDIF
      GO TO 10
C
   11 DS = S(I) - S(I-1)
      T = (SS - S(I-1)) / DS
C
      CX1 = DS*XS(I-1) - X(I) + X(I-1)
      CX2 = DS*XS(I)   - X(I) + X(I-1)
      XD = X(I) - X(I-1) + (1.0-4.0*T+3.0*T*T)*CX1 + T*(3.0*T-2.0)*CX2
      XDD = (6.0*T-4.0)*CX1 + (6.0*T-2.0)*CX2
      XDDD = 6.0*CX1 + 6.0*CX2
C
      CY1 = DS*YS(I-1) - Y(I) + Y(I-1)
      CY2 = DS*YS(I)   - Y(I) + Y(I-1)
      YD = Y(I) - Y(I-1) + (1.0-4.0*T+3.0*T*T)*CY1 + T*(3.0*T-2.0)*CY2
      YDD = (6.0*T-4.0)*CY1 + (6.0*T-2.0)*CY2
      YDDD = 6.0*CY1 + 6.0*CY2
C
      SD = SQRT(XD*XD + YD*YD)
      SD = MAX(SD,0.001*DS)
C
      BOT = SD**3
      DBOTDT = 3.0*SD*(XD*XDD + YD*YDD)
C
      TOP = XD*YDD - YD*XDD
      DTOPDT = XD*YDDD - YD*XDDD
C
      CURVS = (DTOPDT*BOT - DBOTDT*TOP) / BOT**2
C
      RETURN
      END ! CURVS

      SUBROUTINE SINVRT(SI,XI,X,XS,S,N)
      DIMENSION X(N), XS(N), S(N)
C-------------------------------------------------------
C     Calculates the "inverse" spline function S(X).    |
C     Since S(X) can be multi-valued or not defined,    |
C     this is not a "black-box" routine.  The calling   |
C     program must pass via SI a sufficiently good      |
C     initial guess for S(XI).                          |
C                                                       |
C     XI      specified X value       (input)           |
C     SI      calculated S(XI) value  (input,output)    |
C     X,XS,S  usual spline arrays     (input)           |
C                                                       |
C-------------------------------------------------------
C
      SISAV = SI
C
      DO 10 ITER=1, 10
        RES  = SEVAL(SI,X,XS,S,N) - XI
        RESP = DEVAL(SI,X,XS,S,N)
        DS = -RES/RESP
        SI = SI + DS
        IF(ABS(DS/(S(N)-S(1))) .LT. 1.0E-5) RETURN
   10 CONTINUE
      WRITE(*,*)
     &  'SINVRT: spline inversion failed. Input value returned.'
      SI = SISAV
C
      RETURN
      END ! SINVRT

      SUBROUTINE SCALC(X,Y,S,N)
      DIMENSION X(N), Y(N), S(N)
C----------------------------------------
C     Calculates the arc length array S  |
C     for a 2-D array of points (X,Y).   |
C----------------------------------------
C
      S(1) = 0.
      DO 10 I=2, N
        S(I) = S(I-1) + SQRT((X(I)-X(I-1))**2 + (Y(I)-Y(I-1))**2)
   10 CONTINUE
C
      RETURN
      END ! SCALC

      SUBROUTINE SPLNXY(X,XS,Y,YS,S,N)
      DIMENSION X(N), XS(N), Y(N), YS(N), S(N)
C-----------------------------------------
C     Splines 2-D shape X(S), Y(S), along |
C     with true arc length parameter S.   |
C-----------------------------------------
      PARAMETER (KMAX=32)
      DIMENSION XT(0:KMAX), YT(0:KMAX)
C
      KK = KMAX
      NPASS = 10
C
C---- set first estimate of arc length parameter
      CALL SCALC(X,Y,S,N)
C
C---- spline X(S) and Y(S)
      CALL SEGSPL(X,XS,S,N)
      CALL SEGSPL(Y,YS,S,N)
C
C---- re-integrate true arc length
      DO 100 IPASS=1, NPASS
C
        SERR = 0.
C
        DS = S(2) - S(1)
        DO I = 2, N
          DX = X(I) - X(I-1)
          DY = Y(I) - Y(I-1)
C
          CX1 = DS*XS(I-1) - DX
          CX2 = DS*XS(I  ) - DX
          CY1 = DS*YS(I-1) - DY
          CY2 = DS*YS(I  ) - DY
C
          XT(0) = 0.
          YT(0) = 0.
          DO K=1, KK-1
            T = FLOAT(K) / FLOAT(KK)
            XT(K) = T*DX + (T-T*T)*((1.0-T)*CX1 - T*CX2)
            YT(K) = T*DY + (T-T*T)*((1.0-T)*CY1 - T*CY2)
          ENDDO
          XT(KK) = DX
          YT(KK) = DY
C
          SINT1 = 0.
          DO K=1, KK
            SINT1 = SINT1
     &            + SQRT((XT(K)-XT(K-1))**2 + (YT(K)-YT(K-1))**2)
          ENDDO
C
          SINT2 = 0.
          DO K=2, KK, 2
            SINT2 = SINT2
     &            + SQRT((XT(K)-XT(K-2))**2 + (YT(K)-YT(K-2))**2)
          ENDDO
C
          SINT = (4.0*SINT1 - SINT2) / 3.0
C
          IF(ABS(SINT-DS) .GT. ABS(SERR))  SERR = SINT - DS
C
          IF(I.LT.N) DS = S(I+1) - S(I)
C
          S(I) = S(I-1) + SQRT(SINT)
        ENDDO
C
        SERR = SERR / (S(N) - S(1))
        WRITE(*,*) IPASS, SERR
C
C------ re-spline X(S) and Y(S)
        CALL SEGSPL(X,XS,S,N)
        CALL SEGSPL(Y,YS,S,N)
C
        IF(ABS(SERR) .LT. 1.0E-7) RETURN
C
 100  CONTINUE
C
      RETURN
      END ! SPLNXY

      SUBROUTINE SEGSPL(X,XS,S,N)
C-----------------------------------------------
C     Splines X(S) array just like SPLINE,      |
C     but allows derivative discontinuities     |
C     at segment joints.  Segment joints are    |
C     defined by identical successive S values. |
C-----------------------------------------------
      DIMENSION X(N), XS(N), S(N)
C
      IF(S(1).EQ.S(2)  ) STOP 'SEGSPL:  First input point duplicated'
      IF(S(N).EQ.S(N-1)) STOP 'SEGSPL:  Last  input point duplicated'
C
      ISEG0 = 1
      DO 10 ISEG=2, N-2
        IF(S(ISEG).EQ.S(ISEG+1)) THEN
         NSEG = ISEG - ISEG0 + 1
         CALL SPLIND(X(ISEG0),XS(ISEG0),S(ISEG0),NSEG,-999.0,-999.0)
         ISEG0 = ISEG+1
        ENDIF
   10 CONTINUE
C
      NSEG = N - ISEG0 + 1
      CALL SPLIND(X(ISEG0),XS(ISEG0),S(ISEG0),NSEG,-999.0,-999.0)
C
      RETURN
      END ! SEGSPL

      SUBROUTINE SEGSPLD(X,XS,S,N,XS1,XS2)
C-----------------------------------------------
C     Splines X(S) array just like SPLIND,      |
C     but allows derivative discontinuities     |
C     at segment joints.  Segment joints are    |
C     defined by identical successive S values. |
C-----------------------------------------------
      DIMENSION X(N), XS(N), S(N)
C
      IF(S(1).EQ.S(2)  ) STOP 'SEGSPL:  First input point duplicated'
      IF(S(N).EQ.S(N-1)) STOP 'SEGSPL:  Last  input point duplicated'
C
      ISEG0 = 1
      DO 10 ISEG=2, N-2
        IF(S(ISEG).EQ.S(ISEG+1)) THEN
         NSEG = ISEG - ISEG0 + 1
         CALL SPLIND(X(ISEG0),XS(ISEG0),S(ISEG0),NSEG,XS1,XS2)
         ISEG0 = ISEG+1
        ENDIF
   10 CONTINUE
C
      NSEG = N - ISEG0 + 1
      CALL SPLIND(X(ISEG0),XS(ISEG0),S(ISEG0),NSEG,XS1,XS2)
C
      RETURN
      END ! SEGSPL

ext: f


text:
program cholspeed
  ! Tests the speed of modified Cholesky factorization

  use flips

  implicit none

  integer :: n,i
  character(len=1) :: type

  real(4), dimension(:,:), allocatable :: samat,scmat
  real(8), dimension(:,:), allocatable :: damat,dcmat
  complex(4), dimension(:,:), allocatable :: camat,ccmat
  complex(8), dimension(:,:), allocatable :: zamat,zcmat

  real(8) :: stime, etime, ttime

  write(*,*) 'Give size:'
  read(*,*) n
  write(*,*) 'Give type:'
  read(*,*) type

  select case (type)

  case ('s')

     allocate(samat(n,n),scmat(n,n))

     ! We do not really have to care about the symmetry
     call random_number(samat)

     ! Just to make sure that amat is positive definite
     do i = 1,n
        samat(i,i) = n
     end do

write(*,*) 'Start!'

     call cpu_time(stime)
     call modcholmat_s(n,samat,scmat)
     call cpu_time(etime)
     ttime = etime - stime

     deallocate(samat,scmat)

  case ('d')

     allocate(damat(n,n),dcmat(n,n))

     ! We do not really have to care about the symmetry
     call random_number(damat)

     ! Just to make sure that amat is positive definite
     do i = 1,n
        damat(i,i) = n
     end do

write(*,*) 'Start!'

     call cpu_time(stime)
     call modcholmat_d(n,damat,dcmat)
     call cpu_time(etime)
     ttime = etime - stime

     deallocate(damat,dcmat)

  case ('c')

     allocate(camat(n,n),ccmat(n,n),samat(n,n),scmat(n,n))

     ! We do not really have to care about the symmetry
     call random_number(samat)
     call random_number(scmat)

     camat = samat + (0.0,1.0) * scmat

     ! Just to make sure that amat is positive definite
     do i = 1,n
        camat(i,i) = n
     end do

write(*,*) 'Start!'

     call cpu_time(stime)
     call modcholmat_c(n,camat,ccmat)
     call cpu_time(etime)
     ttime = etime - stime

     deallocate(camat,ccmat,samat,scmat)

  case ('z')

     allocate(zamat(n,n),zcmat(n,n),damat(n,n),dcmat(n,n))

     ! We do not really have to care about the symmetry
     call random_number(damat)
     call random_number(dcmat)

     zamat = damat + (0.0,1.0) * dcmat

     ! Just to make sure that amat is positive definite
     do i = 1,n
        zamat(i,i) = n
     end do

write(*,*) 'Start!'

     call cpu_time(stime)
     call modcholmat_z(n,zamat,zcmat)
     call cpu_time(etime)
     ttime = etime - stime

     deallocate(zamat,zcmat,damat,dcmat)

  case default

     write(*,*) 'Unknown type. Exiting'
     stop

  end select

  write(*,*) 'Size:',n
  write(*,*) 'Type:',type
  write(*,*) 'CPU Time:',ttime

end program cholspeed

max_stars_repo_name: morispaa/flips

max_issues_repo_name: morispaa/flips

max_forks_repo_name: morispaa/flips


text:
c ------------------------------------------------------
      subroutine count_in_bins(histofill)
c ------------------------------------------------------

      implicit none

#include "local.inc"
#include "constant.inc"

	integer histofill
      real datawt, bakwt, mcwt

      datawt = 0.
	bakwt = 0.
	mcwt = 0.

	if (acceptance.ne.0) then
	   datawt = 1. / acceptance
         bakwt = (baknorm / acceptance)
	else
	   datawt = 1.
	   bakwt = baknorm
	end if

	mcwt = mcwtq2 * f2weight * wt_vtx * bpcor

c ---------------------------
c Count data events
c ---------------------------
      if (histofill.eq.1) then
         bindat = bindat + (datawt)
	   datraw = datraw + 1

c ---------------------------
c Count background events
c ---------------------------
      else if (histofill.eq.2) then
         binbgd = binbgd + (bakwt)
	   bakraw = bakraw + 1

c ---------------------------
c Count MC events
c ---------------------------
      else if (histofill.eq.3) then
         binmc = binmc + (mcwt)
	   mcraw = mcraw + 1

      end if

      return
      end

ext: fpp


text:
PROGRAM xcodeml_sample
  IMPLICIT NONE

  INTEGER :: my_integer
  REAL(KIND=8), DIMENSION(10) :: my_double_precision_real_array

  DO my_integer = 1, 10
    my_double_precision_real_array(my_integer) = my_integer ** 2
  END DO
END PROGRAM xcodeml_sample

max_stars_repo_name: peclatj/claw-compiler

max_issues_repo_name: peclatj/claw-compiler

max_forks_repo_name: peclatj/claw-compiler


text:
c -----------------------------------------------------------------------------
*dk midpt
      subroutine midpt(x,x1,x2,nn)

c...  This routine finds the midpoint x along the shorter great circle
c...  arc between points x1 and x2 on the unit sphere.

      integer, intent(in) :: nn
      real, intent(inout) :: x(nn,3)
      real, intent(in)    :: x1(nn,3), x2(nn,3)

      do j=1,3
         x(1,j) = x1(1,j) + x2(1,j)
      end do

      xnorm = 1./sqrt(x(1,1)**2 + x(1,2)**2 + x(1,3)**2)

      do j=1,3
         x(1,j) = xnorm*x(1,j)
      end do

      end

c -----------------------------------------------------------------------------
*dk grdgen
      subroutine grdgen(xn,nt)
      implicit none
c...  This routine generates the nodal coordinates xn for an
c...  icosahedral grid on the unit sphere.  The grid resolution
c...  corresponds to a subdivision of the edges of the original
c...  icosahedral triangles into nt equal parts.

c...  Note: nt is input as mt

      integer, intent(inout) :: nt
      real, intent(inout) :: xn((nt+1),(nt+1),10,3)

      real     :: fifthpi, w, cosw, sinw, lvt, nn, phi
      integer  :: id, k, i1, i2, j1, j2, l, l2, m, sgn

      fifthpi = 0.4*asin(1.)
      w       = 2.0*acos(1./(2.*sin(fifthpi)))
      cosw    = cos(w)
      sinw    = sin(w)
      lvt     = 1.45*log(real(nt))
      nn      = (nt+1)**2*10

 900  format (I4,I4,I4,I4,F12.8)

      do id=1,10

         sgn = 1.
         if(id .ge. 6) sgn = -1.
         phi = (2*mod(id, 5) - 3 + (id - 1)/5)*fifthpi

         xn(   1,   1,id,1) =  0.
         xn(   1,   1,id,2) =  0.
         xn(   1,   1,id,3) =  sgn
         xn(nt+1,   1,id,1) =  sinw*cos(phi)
         xn(nt+1,   1,id,2) =  sinw*sin(phi)
         xn(nt+1,   1,id,3) =  cosw*sgn
         xn(   1,nt+1,id,1) =  sinw*cos(phi + fifthpi + fifthpi)
         xn(   1,nt+1,id,2) =  sinw*sin(phi + fifthpi + fifthpi)
         xn(   1,nt+1,id,3) =  cosw*sgn
         xn(nt+1,nt+1,id,1) =  sinw*cos(phi + fifthpi)
         xn(nt+1,nt+1,id,2) =  sinw*sin(phi + fifthpi)
         xn(nt+1,nt+1,id,3) = -cosw*sgn

c         do k=0,nint(lvt-1)
c         do k=0,int(lvt-1)
         do k=0,lvt-1

            m  = 2**k
            l  = nt/m
            l2 = l/2

c...        rows of diamond--

            do j1=1,m+1
               do j2=1,m
                     i1 = (j1-1)*l + 1
                     i2 = (j2-1)*l + l2 + 1
                     call midpt(xn(i1,i2,id,1),xn(i1,i2-l2,id,1),
     &                          xn(i1,i2+l2,id,1),nn)
               end do
            end do

c...        columns of diamond--

            do j1=1,m+1
               do j2=1,m
                     i1 = (j2-1)*l + l2 + 1
                     i2 = (j1-1)*l + 1
                     call midpt(xn(i1,i2,id,1),xn(i1-l2,i2,id,1),
     &                          xn(i1+l2,i2,id,1),nn)
               end do
            end do

c...        diagonals of diamond--

            do j1=1,m
               do j2=1,m
                     i1 = (j1-1)*l + l2 + 1
                     i2 = (j2-1)*l + l2 + 1
                     call midpt(xn(i1,i2,id,1),xn(i1-l2,i2+l2,id,1),
     &                          xn(i1+l2,i2-l2,id,1),nn)
               end do
            end do

         end do

      end do

c...

      do k=1,3
         do id=1,10
            do j1=1,nt+1
               do i1=1,nt+1
                  write(6, 900) i1,j1,id,k,xn(i1,j1,id,k)
               enddo
            enddo
         enddo
      enddo

      end

c -----------------------------------------------------------------------------
*dk grdgen
      subroutine grdgentest(nt)
      implicit none
c...  This routine generates the nodal coordinates xn for an
c...  icosahedral grid on the unit sphere.  The grid resolution
c...  corresponds to a subdivision of the edges of the original
c...  icosahedral triangles into nt equal parts.

c...  Note: nt is input as mt

      integer, intent(inout) :: nt
      real                   :: xn((nt+1),(nt+1),10,3)

      real     :: fifthpi, w, cosw, sinw, lvt, nn, phi
      integer  :: id, k, i1, i2, j1, j2, l, l2, m, sgn

      fifthpi = 0.4*asin(1.)
      w       = 2.0*acos(1./(2.*sin(fifthpi)))
      cosw    = cos(w)
      sinw    = sin(w)
      lvt     = 1.45*log(real(nt))
      nn      = (nt+1)**2*10

 900  format (I4,I4,I4,I4,F12.8)
 901  format (F16.8,F16.8,F16.8)
 902  format (I4,I4,F16.8,F16.8,F16.8)

      do id=1,1

         sgn = 1.
         if(id .ge. 6) sgn = -1.
         phi = (2*mod(id, 5) - 3 + (id - 1)/5)*fifthpi

         xn(   1,   1,id,1) =  0.
         xn(   1,   1,id,2) =  0.
         xn(   1,   1,id,3) =  sgn
         write(6, 901) xn(1,1,id,1),xn(1,1,id,2),xn(1,1,id,3)
         xn(nt+1,   1,id,1) =  sinw*cos(phi)
         xn(nt+1,   1,id,2) =  sinw*sin(phi)
         xn(nt+1,   1,id,3) =  cosw*sgn
         write(6, 901) xn(nt+1,1,id,1),xn(nt+1,1,id,2),xn(nt+1,1,id,3)
         xn(   1,nt+1,id,1) =  sinw*cos(phi + fifthpi + fifthpi)
         xn(   1,nt+1,id,2) =  sinw*sin(phi + fifthpi + fifthpi)
         xn(   1,nt+1,id,3) =  cosw*sgn
         write(6, 901) xn(1,nt+1,id,1),xn(1,nt+1,id,2),xn(1,nt+1,id,3)
         xn(nt+1,nt+1,id,1) =  sinw*cos(phi + fifthpi)
         xn(nt+1,nt+1,id,2) =  sinw*sin(phi + fifthpi)
         xn(nt+1,nt+1,id,3) = -cosw*sgn
         write(6, 901) xn(nt+1,nt+1,id,1),xn(nt+1,nt+1,id,3),
     &        xn(nt+1,1,id,2)

c         do k=0,nint(lvt-1)
c         do k=0,int(lvt-1)
         do k=0,lvt-1

            m  = 2**k
            l  = nt/m
            l2 = l/2

c...        rows of diamond--

            do j1=1,m+1
               do j2=1,m
                     i1 = (j1-1)*l + 1
                     i2 = (j2-1)*l + l2 + 1
                     call midpt(xn(i1,i2,id,1),xn(i1,i2-l2,id,1),
     &                          xn(i1,i2+l2,id,1),nn)
                     write(6, 902) i1,i2,xn(i1,i2,id,1),xn(i1,i2,id,2),
     &                    xn(i1,i2,id,3)
               end do
            end do

c...        columns of diamond--

            do j1=1,m+1
               do j2=1,m
                     i1 = (j2-1)*l + l2 + 1
                     i2 = (j1-1)*l + 1
                     call midpt(xn(i1,i2,id,1),xn(i1-l2,i2,id,1),
     &                          xn(i1+l2,i2,id,1),nn)
                     write(6, 901) xn(i1,i2,id,1),xn(i1,i2,id,2),
     &                    xn(i1,i2,id,3)
               end do
            end do

c...        diagonals of diamond--

            do j1=1,m
               do j2=1,m
                     i1 = (j1-1)*l + l2 + 1
                     i2 = (j2-1)*l + l2 + 1
                     call midpt(xn(i1,i2,id,1),xn(i1-l2,i2+l2,id,1),
     &                          xn(i1+l2,i2-l2,id,1),nn)
                     write(6, 901) xn(i1,i2,id,1),xn(i1,i2,id,2),
     &                    xn(i1,i2,id,3)
               end do
            end do

         end do

      end do

      end


text:
C
C
C       $Id: terp1.f,v 1.1 1994-08-11 23:17:36 boote Exp $
C
C PACKAGE CUBSPL         NOTE--DOCUMENTATION FOR INDIVIDUAL ROUTINES
C                              FOLLOWS THE GENERAL PACKAGE INFORMATION
C
C LATEST REVISION        JANUARY 1985
C
C PURPOSE                TO PERFORM ONE AND TWO-DIMENSIONAL CUBIC SPLINE
C                        INTERPOLATION WITH CHOICE OF BOUNDARY
C                        CONDITIONS.  THE FUNCTION AND SELECTED
C                        DERIVATIVES MAY BE EVALUATED AT ANY POINT WHERE
C                        INTERPOLATION IS REQUIRED.  IN THE
C                        TWO-DIMENSIONAL CASE, THE GIVEN DATA POINTS
C                        MUST BE ON A RECTANGULAR GRID, WHICH NEED NOT
C                        BE EQUALLY SPACED.  THE PACKAGE CUBSPL CONTAINS
C                        SEVEN ROUTINES.
C
C                        SUBROUTINE COEFF1
C                          COMPUTES THE COEFFICIENTS FOR ONE-DIMENSIONAL
C                          CUBIC SPLINE INTERPOLATION USING ONE OF
C                          THE FOLLOWING BOUNDARY CONDITIONS AT
C                          EACH END OF THE RANGE.
C                            . SECOND DERIVATIVE GIVEN AT BOUNDARY.
C                            . FIRST DERIVATIVE GIVEN AT BOUNDARY.
C                            . PERIODIC BOUNDARY CONDITION.
C                            . FIRST DERIVATIVE DETERMINED BY FITTING A
C                              CUBIC TO THE FOUR POINTS NEAREST TO THE
C                              BOUNDARY.
C
C                        SUBROUTINE TERP1
C                          USING THE COEFFICIENTS COMPUTED BY COEFF1,
C                          THIS ROUTINE EVALUATES THE FUNCTION AND/OR
C                          FIRST AND SECOND DERIVATIVES AT ANY POINT
C                          WHERE INTERPOLATION IS REQUIRED.
C
C                        SUBROUTINE COEFF2
C                          COMPUTES THE COEFFICIENTS FOR TWO-DIMENSIONAL
C                          BICUBIC SPLINE INTERPOLATION WITH THE SAME
C                          CHOICE OF BOUNDARY CONDITIONS AS FOR COEFF1.
C
C                        FUNCTION TERP2
C                          USING THE COEFFICIENTS PRODUCED BY COEFF2,
C                          THIS SUBROUTINE EVALUATES THE FUNCTION OR A
C                          SELECTED DERIVATIVE AT ANY POINT WHERE
C                          TWO-DIMENSIONAL INTERPOLATION IS REQUIRED.
C
C                        SUBROUTINE TRIP
C                          A SIMPLE PERIODIC, TRIDIAGONAL LINEAR
C                          EQUATION SOLVER USED BY COEFF1.
C
C                        SUBROUTINE SEARCH
C                          PERFORMS A BINARY SEARCH IN A ONE-DIMENSIONAL
C                          FLOATING POINT TABLE ARRANGED IN ASCENDING
C                          ORDER.  THIS ROUTINE IS CALLED BY TERP1 AND
C                          TERP2.
C
C                        SUBROUTINE INTERP
C                          GIVEN COEFFICIENTS PROVIDED BY COEFF1 AND THE
C                          POSITION OF THE INTERPOLATION POINT IN THE
C                          INDEPENDENT VARIABLE TABLE, THIS ROUTINE
C                          PERFORMS ONE-DIMENSIONAL INTERPOLATION FOR
C                          THE FUNCTION VALUE, FIRST AND SECOND
C                          DERIVATIVE, AS DESIRED.  THIS ROUTINE IS
C                          CALLED BY TERP1 AND TERP2.
C
C USAGE                  FOR ONE-DIMENSIONAL CUBIC SPLINE INTERPOLATION,
C                        THE USER FIRST CALLS COEFF1 BY
C
C                          CALL COEFF1 (N,X,F,W,IOP,INT,WK)
C
C                        THIS SUBROUTINE RETURNS THE COEFFICIENTS
C                        NEEDED FOR THE SUBSEQUENT INTERPOLATION
C                        IN THE ARRAY W.  THE USER THEN CALLS
C                        SUBROUTINE TERP1 BY
C
C                          CALL TERP1 (N,X,F,W,Y,INT,TAB,ITAB)
C
C                        IN ORDER TO COMPUTE THE VALUE OF THE
C                        FUNCTION AND/OR ITS DERIVATIVES.  THE USER
C                        SPECIFIES Y, THE VALUE OF THE INDEPENDENT
C                        VARIABLE WHERE THE INTERPOLATION IS TO BE
C                        PERFORMED.  THE INTERPOLATED VALUES ARE
C                        RETURNED IN TAB.  THE PARAMETERS
C                        N,X,F,W, AND INT MUST NOT BE CHANGED
C                        BETWEEN THE SUBROUTINE CALLS.
C
C                        FOR TWO-DIMENSIONAL CUBIC SPLINE INTERPOLATION
C                        THE USER FIRST CALLS COEFF2 BY
C
C                          CALL COEFF2 (NX,X,NY,Y,F,FXX,FYY,FXXYY,
C                                       IDM,IBD,WK)
C
C                        THIS SUBROUTINE RETURNS THE COEFFICIENTS
C                        NEEDED FOR THE SUBSEQUENT INTERPOLATION IN
C                        THE ARRAYS FXX, FYY, FXXYY.  THE USER THEN
C                        CALLS THE ROUTINE TERP2 BY
C
C                          R = TERP2 (XB,YB,NX,X,NY,Y,F,FXX,FYY,FXXYY,
C                                     IDM,IXD,IYD)
C
C                        IN ORDER TO PERFORM THE INTERPOLATION AT THE
C                        POINT SPECIFIED BY THE VALUES OF XB AND YB.
C                        DEPENDING ON THE VALUES INPUT IN IXD AND IYD,
C                        THE ROUTINE RETURNS AN INTERPOLATED VALUE
C                        FOR THE FUNCTION OR ONE OF ITS PARTIAL
C                        DERIVATIVES.  THE PARAMETERS NX,X,NY,Y,F,FXX,
C                        FYY,FXXYY, AND IDM MUST NOT BE CHANGED
C                        BETWEEN THE CALLS TO COEFF2 AND TERP2.
C
C SPECIAL CONDITIONS     TABLES OF INDEPENDENT VARIABLES MUST BE
C                        ARRANGED IN ASCENDING ORDER.  FOR
C                        TWO-DIMENSIONAL INTERPOLATION, THE DATA POINTS
C                        MUST LIE ON A RECTANGULAR MESH.
C
C I/O                    NONE
C
C PRECISION              SINGLE
C
C REQUIRED LIBRARY       NONE
C FILES
C
C LANGUAGE               FORTRAN
C
C HISTORY                THIS PACKAGE IS BASED ON THE ROUTINES
C                            LA E 102A, SPL1D1
C                            LA E 103A, SPL1D2
C                            LA E 104A, SPL2D1
C                            LA E 105A, SPL2D2
C                        OF THE LOS ALAMOS CUBIC SPLINE PACKAGE WRITTEN
C                        BY THOMAS J. JORDAN AND BERTHA FAGAN, 1968.
C                        THE ROUTINES HAVE BEEN STREAMLINED AND
C                        STANDARDIZED.  THE ALGORITHM FOR
C                        TWO-DIMENSIONAL INTERPOLATION IS CONSIDERABLY
C                        MODIFIED.
C
C ALGORITHM              FOR ONE-DIMENSIONAL INTERPOLATION, THE CUBIC
C                        SPLINE IS EXPRESSED IN TERMS OF THE FUNCTION
C                        VALUES AND SECOND DERIVATIVES AT THE DATA
C                        POINTS.  THE SECOND DERIVATIVES ARE
C                        DETERMINED FROM A TRIDIAGONAL LINEAR SYSTEM
C                        WHICH DESCRIBES THE CONTINUITY OF THE FIRST
C                        DERIVATIVE AND INCORPORATES THE GIVEN
C                        BOUNDARY CONDITIONS.  COEFF1  SETS UP THIS
C                        SYSTEM AND CALLS SUBROUTINE TRIP TO SOLVE IT.
C
C                        THE CUBIC SEGMENT BETWEEN TWO ADJACENT
C                        TABULAR POINTS IS CONSTRUCTED FROM THE
C                        FUNCTION VALUES AND SECOND DERIVATIVES AT
C                        THESE POINTS.  THESE PROVIDE THE FOUR
C                        CONSTANTS NEEDED TO DEFINE THE CUBIC
C                        UNIQUELY.  FROM THIS CUBIC, VALUES OF THE
C                        FUNCTION AND ITS FIRST AND SECOND
C                        DERIVATIVES ARE READILY DETERMINED AT ANY
C                        INTERMEDIATE POINT.  ONE-DIMENSIONAL
C                        INTERPOLATION IS PERFORMED BY THE ROUTINE
C                        TERP1.  FOR TWO-DIMENSIONAL INTERPOLATION,
C                        THE BICUBIC SPLINE IS DESCRIBED IN TERMS OF
C                        VALUES OF F,FXX,FYY, AND FXXYY  AT EACH
C                        POINT ON THE GIVEN TWO-DIMENSIONAL
C                        RECTANGULAR GRID OF DATA POINTS.  HERE F
C                        IS THE FUNCTION VALUE,
C
C                          FXX = (D/DX)**2*F
C
C                        AND SO ON.  THE COEFFICIENTS ARE DETERMINED
C                        BY COEFF2, WHICH USES SUCCESSIVE
C                        APPLICATIONS OF COEFF1.
C
C                        1.  THE ARRAY FXX IS DETERMINED BY APPLYING
C                            COEFF1 TO F ALONG EACH LINE IN THE
C                            X-DIRECTION.
C
C                        2.  THE ARRAY FYY IS DETERMINED BY APPLYING
C                            COEFF1 TO F ALONG EACH LINE IN THE
C                            Y-DIRECTION.
C
C                        3.  FXXYY IS DETERMINED ON THE CONSTANT Y
C                            BOUNDARIES BY APPLYING COEFF1 TO FYY
C                            ALONG THESE BOUNDARIES.
C
C                        4.  THE REMAINDER OF THE ARRAY FXXYY IS
C                            DETERMINED BY APPLYING COEFF1 TO FXX
C                            ALONG EACH LINE IN THE Y-DIRECTION.
C
C                        THE BICUBIC WITHIN ANY RECTANGULAR ELEMENT
C                        OF THE GRID IS CONSTRUCTED FROM THE VALUES
C                        OF F,FXX,FYY,FXXYY AT THE FOUR CORNERS.
C                        THESE PROVIDE THE 16 CONSTANTS NECESSARY
C                        TO DEFINE THE BICUBIC UNIQUELY.  TO FIND
C                        THE VALUE OF F CORRESPONDING TO A POINT
C                        (XB,YB) WITHIN THE ELEMENTARY RECTANGLE,
C                        (X(I),Y(J)),(X(I+1),Y(J)),(X(I),Y(J+1)),
C                        (X(I+1),Y(J+1)), FIVE ONE DIMENSIONAL
C                        INTERPOLATIONS ARE PERFORMED.
C
C                        1.  F AT (XB,Y(J)) AND (XB,Y(J+1)) ARE
C                            FOUND BY INTERPOLATING F IN THE
C                            X-DIRECTION USING FXX. (TWO INTERPOLATIONS)
C
C                        2.  FYY AT (XB,Y(J)) AND (XB,Y(J+1)) ARE
C                            FOUND BY INTERPOLATING FYY IN THE
C                            X-DIRECTION USING FXXYY. (TWO
C                            INTERPOLATIONS.)
C
C                        3.  FINALLY F AT (XB,YB) IS FOUND BY
C                        INTERPOLATING BETWEEN F(XB,Y(J)) AND
C                        F(XB,Y(J+1)) IN THE Y-DIRECTION USING
C                        VALUES OF FYY(XB,Y(J)) AND FYY(XB,Y(J+1))
C                        OBTAINED ABOVE. (ONE INTERPOLATION).
C
C                        TWO-DIMENSIONAL INTERPOLATION IS PERFORMED
C                        IN TERP2.
C
C REFERENCES             FOR GREATER DETAIL, SEE J.L.WALSH,
C                        J.H.AHLBERG, E.N.NILSEN, BEST APPROXIMATION
C                        PROPERTIES OF THE SPLINE FIT, JOURNAL OF
C                        MATHEMATICS AND MECHANICS, VOL.II(1962),
C                        225-234.
C
C                        T.L. JORDAN, SMOOTHING AND MULTIVARIABLE
C                        INTERPOLATION WITH SPLINES, LOS ALAMOS
C                        REPORT, LA-3137, 1965.
C
C ACCURACY               NEAR MACHINE ACCURACY WAS OBTAINED WHEN
C                        INTERPOLATING A CUBIC IN ONE DIMENSION
C                        OR A BICUBIC IN TWO DIMENSIONS.
C
C PORTABILITY            FULLY PORTABLE WITH RESPECT TO FORTRAN 66.
C***********************************************************************
C
C SUBROUTINE COEFF1 (N,X,F,W,IOP,INT,WK)
C
C DIMENSION OF           X(N),F(INT*(N-1)+1),W(INT*(N-1)+1),IOP(2),
C ARGUMENTS              WK(3*N+1)
C
C PURPOSE                SUBROUTINE COEFF1 COMPUTES THE COEFFICIENTS
C                        FOR ONE-DIMENSIONAL CUBIC SPLINE
C                        INTERPOLATION USING ONE OF THE FOLLOWING
C                        BOUNDARY CONDITIONS AT EACH END OF THE
C                        RANGE
C
C                        .  SECOND DERIVATIVES GIVEN AT BOUNDARY
C                        .  FIRST DERIVATIVE GIVEN AT BOUNDARY
C                        .  PERIODIC BOUNDARY CONDITIONS
C                        .  FIRST DERIVATIVE CALCULATED BY FITTING A
C                           CUBIC TO THE FOUR POINTS NEAREST TO THE
C                           BOUNDARY
C
C                        NOTE THAT TERP1 MUST BE CALLED TO PERFORM
C                        THE INTERPOLATION.
C
C USAGE                  CALL COEFF1 (N,X,F,W,IOP,INT,WK)
C
C ARGUMENTS
C
C ON INPUT               N
C                          THE NUMBER OF DATA POINTS.  N MUST BE AT
C                          LEAST 4.
C
C                        X
C                          TABLE OF N INDEPENDENT VARIABLE VALUES IN
C                          ASCENDING ORDER.  DIMENSION OF X IN CALLING
C                          PROGRAM MUST BE AT LEAST N.
C
C                        F
C                          TABLE OF N CORRESPONDING DEPENDENT VARIABLE
C                          VALUES.  THE VALUES ARE SEPARATED BY INTERVAL
C                          INT.  THIS IS USUALLY UNITY FOR
C                          ONE-DIMENSIONAL INTERPOLATION.  OTHER VALUES
C                          ARE USEFUL WHEN COEFF1 IS INCORPORATED IN A
C                          TWO-DIMENSIONAL INTERPOLATION SCHEME (SEE
C                          COEFF2).  DIMENSION OF F IN THE CALLING
C                          PROGRAM MUST BE AT LEAST (INT*(N-1)+1).
C
C                        IOP
C                          TWO ELEMENT INTEGER ARRAY DEFINING BOUNDARY
C                          CONDITIONS AT X(1) AND X(N) ACCORDING TO THE
C                          FOLLOWING CODE.
C
C                          FOR IOP(1)
C                          = 1  SECOND DERIVATIVE GIVEN AT X(1).  PLACE
C                               VALUE OF SECOND DERIVATIVE IN W(1)
C                               BEFORE CALL TO COEFF1.
C                          = 2  FIRST DERIVATIVE GIVEN AT X(1).  PLACE
C                               VALUE OF FIRST DERIVATIVE IN W(1) BEFORE
C                               CALL.
C                          = 3  PERIODIC BOUNDARY CONDITION.  X(1) AND
C                               X(N) ARE EQUIVALENT POINTS.  F(1) AND
C                               F(INT*(N-1)+1) ARE EQUAL.
C                          = 4  THE FIRST DERIVATIVE AT X(1) IS
C                               CALCULATED BY FITTING A CUBIC TO POINTS
C                               X(1) THROUGH X(4).
C                          SIMILARLY, IOP(2) DEFINES THE BOUNDARY
C                          CONDITION AT X(N).  WHEN IOP(2) = 1 (OR 2),
C                          THE VALUE OF THE SECOND (OR FIRST) DERIVATIVE
C                          MUST BE PLACED IN W(INT*(N-1)+1).  NOTE THAT
C                          IF IOP(1) = 3, CONSISTENCY DEMANDS THAT
C                          IOP(2) = 3 ALSO.
C
C                        INT
C                          SPACING IN F AND W TABLES.  FOR
C                          ONE-DIMENSIONAL INTERPOLATION THIS WILL
C                          USUALLY BE UNITY.
C
C                        WK
C                          WORK AREA OF DIMENSION AT LEAST (3*N+1).
C
C ON OUTPUT              W
C                          TABLE OF SECOND DERIVATIVES CORRESPONDING TO
C                          GIVEN X AND F VALUES.  THE SEPARATION OF
C                          TABULAR ENTRIES IS INT (SEE ABOVE).
C                          DIMENSION OF W IN THE CALLING PROGRAM MUST BE
C                          AT LEAST (INT*(N-1)+1).
C
C                          THE ARRAYS X, F, W ARE USED AS INPUT FOR THE
C                          ROUTINE TERP1, WHICH PERFORMS INTERPOLATION
C                          AT A GIVEN VALUE OF THE INDEPENDENT VARIABLE.
C
C TIMING                 THE TIMING IS LINEARLY PROPORTIONAL TO N, THE
C                        NUMBER OF DATA POINTS.
C***********************************************************************
C
C     SUBROUTINE TERP1 (N,X,F,W,Y,INT,TAB,ITAB)
C
C
C DIMENSION OF           X(N),F(INT*(N-1)+1),W(INT*(N-1)+1),TAB(3),
C ARGUMENTS              ITAB(3)
C
C PURPOSE                USING THE COEFFICIENTS COMPUTED BY COEFF1,
C                        THIS ROUTINE EVALUATES THE FUNCTION AND/OR
C                        FIRST AND SECOND DERIVATIVES AT ANY POINT.
C
C USAGE                  CALL TERP1 (N,X,F,W,Y,INT,TAB,ITAB)
C
C ARGUMENTS
C
C ON INPUT               N
C                          THE NUMBER OF DATA POINTS.  N MUST BE AT
C                          LEAST 4.
C
C                        X
C                          TABLE OF N INDEPENDENT VARIABLE VALUES IN
C                          ASCENDING ORDER.  DIMENSION OF X IN THE
C                          CALLING PROGRAM MUST BE AT LEAST N.
C
C                        F
C                          TABLE OF N CORRESPONDING DEPENDENT VARIABLE
C                          VALUES SEPARATED BY INTERVAL INT, USUALLY
C                          UNITY FOR ONE-DIMENSIONAL INTERPOLATION.
C                          DIMENSION OF F IN THE CALLING PROGRAM MUST BE
C                          AT LEAST (INT*(N-1)+1).
C
C                        W
C                          TABLE OF SECOND DERIVATIVES COMPUTED BY
C                          COEFF1.  THE SEPARATION OF TABULAR ENTRIES IS
C                          INT.  DIMENSION OF W IN THE CALLING PROGRAM
C                          MUST BE AT LEAST (INT*(N-1)+1).
C
C                        Y
C                          VALUE OF THE INDEPENDENT VARIABLE AT WHICH
C                          INTERPOLATION IS REQUIRED.  IF Y LIES OUTSIDE
C                          THE RANGE OF THE TABLE, EXTRAPOLATION TAKES
C                          PLACE.
C
C                        INT
C                          SPACING OF TABULAR ENTRIES IN F AND W ARRAYS.
C                          THIS IS USUALLY UNITY FOR ONE-DIMENSIONAL
C                          INTERPOLATION.
C
C                        ITAB
C                          THREE ELEMENT INTEGER ARRAY DEFINING
C                          INTERPOLATION TO BE PERFORMED AT Y.
C                            IF ITAB(1) = 1, THE FUNCTION VALUE IS
C                                            RETURNED IN TAB(1).
C                            IF ITAB(2) = 1, THE FIRST DERIVATIVE IS
C                                            RETURNED IN TAB(2).
C                            IF ITAB(3) = 1, THE SECOND DERIVATIVE IS
C                                            RETURNED IN TAB(3).
C                          IF ITAB(I) = 0 FOR I = 1, 2 OR 3, THE
C                          CORRESPONDING FUNCTION VALUE OR DERIVATIVE IS
C                          NOT COMPUTED AND TAB(I) IS NOT REFERENCED.
C
C ON OUTPUT              TAB
C                          THREE ELEMENT ARRAY IN WHICH INTERPOLATED
C                          FUNCTION VALUE, FIRST AND SECOND DERIVATIVES
C                          ARE RETURNED AS DICTATED BY ITAB (SEE ABOVE).
C
C TIMING                 THIS PROCEDURE IS FAST.  THE MAXIMUM TIME FOR
C                        THE BINARY SEARCH IS PROPORTIONAL TO ALOG(N).
C                        THE TIME FOR FUNCTION AND DERIVATIVE EVALUATION
C                        IS INDEPENDENT OF N.
C***********************************************************************
C
C     SUBROUTINE COEFF2 (NX,X,NY,Y,F,FXX,FYY,FXXYY,IDM,IBD,WK)
C
C
C DIMENSION OF           X(NX),Y(NY),F(IDM,NY),FXX(IDM,NY),FYY(IDM,NY),
C ARGUMENTS              FXXYY(IDM,NY),IBD(4),WK(3*MAX0(NX,NY)+1)
C                        (IDM MUST BE .GE. NX)
C
C PURPOSE                SUBROUTINE COEFF2 COMPUTES THE COEFFICIENTS
C                        FOR TWO-DIMENSIONAL BICUBIC SPLINE
C                        INTERPOLATION WITH THE SAME CHOICE OF
C                        BOUNDARY CONDITIONS AS FOR COEFF1.  TERP2
C                        IS CALLED TO PERFORM INTERPOLATION.
C
C USAGE                  CALL COEFF2 (NX,X,NY,Y,F,FXX,FYY,FXXYY,IDM,IBD,
C                                     WK)
C
C ARGUMENTS
C
C ON INPUT               NX
C                          NUMBER OF GRID POINTS IN THE X-DIRECTION.  NX
C                          MUST BE AT LEAST 4.
C
C                        X
C                          TABLE OF NX VALUES OF THE FIRST INDEPENDENT
C                          VARIABLE ARRANGED IN ASCENDING ORDER.
C                          DIMENSION OF X IN THE CALLING PROGRAM MUST BE
C                          AT LEAST NX.
C
C                        NY
C                          NUMBER OF GRID POINTS IN THE Y-DIRECTION.  NY
C                          MUST BE AT LEAST 4.
C
C                        Y
C                          TABLE OF NY VALUES OF THE SECOND INDEPENDENT
C                          VARIABLE ARRANGED IN ASCENDING ORDER.
C                          DIMENSION OF Y IN THE CALLING PROGRAM MUST BE
C                          AT LEAST NY.
C
C                        F
C                          TWO DIMENSIONAL ARRAY OF FUNCTION VALUES AT
C                          THE GRID POINTS DEFINED BY THE ARRAYS X AND
C                          Y.  DIMENSION OF F IN THE CALLING PROGRAM IS
C                          (IDM, NYY) WHERE
C                              IDM .GE. NX
C                              NYY .GE. NY
C
C                        IDM
C                          FIRST DIMENSION IN THE CALLING PROGRAM OF
C                          ARRAYS F, FXX, FYY, FXXYY.  IDM MUST BE AT
C                          LEAST NX.
C
C                        IBD
C                          FOUR ELEMENT INTEGER ARRAY DEFINING BOUNDARY
C                          CONDITIONS ACCORDING TO THE FOLLOWING CODE.
C                          FOR IBD(1)
C                          = 1  THE SECOND DERIVATIVE OF F WITH RESPECT
C                               TO X IS GIVEN AT (X(1),Y(J)) FOR
C                               J = 1,NY,1.  VALUES OF THIS SECOND
C                               DERIVATIVE MUST BE PLACED IN FXX(1,J)
C                               FOR J = 1,NY,1, BEFORE CALLING COEFF2.
C                          = 2  THE FIRST DERIVATIVE OF F WITH RESPECT
C                               TO X IS GIVEN AT (X(1),Y(J)) FOR
C                               J = 1,NY,1.  VALUES OF THE DERIVATIVE
C                               MUST BE PLACED IN FXX(1,J) FOR
C                               J = 1,NY,1 BEFORE CALLING COEFF2.
C                          = 3  PERIODIC BOUNDARY CONDITION IN THE
C                               X-DIRECTION.  (X(1),Y(J)) AND
C                               AND (X(NX),Y(J)) ARE EQUIVALENT POINTS
C                               FOR J = 1,NY,1.  F(1,J) AND F(NX,J) ARE
C                               EQUAL.
C                          = 4  THE FIRST DERIVATIVE OF F WITH RESPECT
C                               TO X AT (X(1),Y(J)) IS COMPUTED BY
C                               FITTING A CUBIC TO F(1,J) THROUGH F(4,J)
C                               FOR J = 1,NY,1.
C
C                          SIMILARLY, IBD(2) DEFINES THE BOUNDARY
C                          CONDITION AT (X(NX),Y(J)) FOR J = 1,NY,1.
C                          WHEN IBD(2) = 1 (OR 2) THE VALUES OF THE
C                          SECOND (OR FIRST) DERIVATIVE OF F WITH
C                          RESPECT TO X ARE PLACED IN FXX(NX,J) FOR
C                          J = 1,NY,1.
C                            NOTE THAT IF IBD(1) = 3, CONSISTENCY
C                            REQUIRES THAT IBD(2) = 3 ALSO.
C                          FOR IBD(3)
C                          = 1  THE SECOND DERIVATIVE OF F WITH RESPECT
C                               TO Y IS GIVEN AT (X(I),Y(1)).  PLACE
C                               VALUES OF THE DERIVATIVE IN FYY(I,1) FOR
C                               I = 1,NX,1 BEFORE CALLING COEFF2.
C                          = 2  THE FIRST DERIVATIVE OF F WITH RESPECT
C                               TO Y IS GIVEN AT (X(I),Y(1)).  VALUES OF
C                               THIS DERIVATIVE MUST BE PLACED IN
C                               FYY(I,1) FOR I = 1,NX,1 BEFORE CALLING
C                               COEFF2.
C                          = 3  PERIODIC BOUNDARY CONDITION IN THE
C                               Y-DIRECTION.  (X(I),Y(1)) AND
C                               (X(I),Y(NY)) ARE EQUIVALENT POINTS.
C                               F(I,1) AND F(I,NY) ARE EQUAL.
C                          = 4  THE FIRST DERIVATIVE OF F WITH RESPECT
C                               TO Y AT (X(I),Y(1)) IS COMPUTED BY
C                               FITTING A CUBIC TO F(I,1) THROUGH F(I,4)
C                               FOR I = 1,NX,1.
C
C                          SIMILARY, IBD(4) DEFINES THE BOUNDARY
C                          CONDITION AT (X(I),Y(NY)) FOR I = 1,NX,1 AND
C                          GIVEN DERIVATIVE VALUES ARE PLACED IN
C                          FYY(I,NY).
C                            NOTE THAT CONSISTENCY DEMANDS THAT IF
C                            IBD(3) = 3, THEN IBD(4) = 3 ALSO.
C
C                        WK
C                          WORK AREA OF DIMENSION AT LEAST
C                          (3*MAX0(NX,NY)+1)
C
C ON OUTPUT              FXX
C                          ARRAY OF SECOND DERIVATIVES OF F WITH RESPECT
C                          TO X COMPUTED BY COEFF2.  FXX(I,J) IS
C                          DERIVATIVE AT (X(I),Y(J)).  AS FOR F,
C                          DIMENSION OF FXX IN THE CALLING PROGRAM IS
C                          (IDM,NYY).
C
C                        FYY
C                          ARRAY OF SECOND DERIVATIVES OF F WITH RESPECT
C                          TO Y COMPUTED BY COEFF2.  DIMENSION OF FYY IN
C                          THE CALLING PROGRAM IS (IDM,NYY).
C
C                        FXXYY
C                          ARRAY OF FOURTH DERIVATIVES
C                          (D/DX)**2*(D/DY)**2*F, COMPUTED BY COEFF2.
C                          DIMENSION OF FXXYY IN THE CALLING PROGRAM IS
C                          (IDM,NYY).
C
C                        THE ARRAYS X, Y, F, FXX, FYY, FXXYY ARE USED AS
C                        INPUT FOR THE ROUTINE TERP2 WHICH PERFORMS
C                        INTERPOLATION AT REQUIRED VALUES OF THE TWO
C                        INDEPENDENT VARIABLES.
C
C TIMING                 THE TIMING IS PROPORTIONAL TO NX*NY.
C***********************************************************************
C
C     FUNCTION TERP2 (XB,YB,NX,X,NY,Y,F,FXX,FYY,FXXYY,IDM,IXD,IYD)
C
C
C DIMENSION OF           X(NX),Y(NY),F(IDM,NY),FXX(IDM,NY),FYY(IDM,NY),
C ARGUMENTS              FXXYY(IDM,NY))
C                        (IDM MUST BE .GE. NX)
C
C PURPOSE                USING THE COEFFICIENTS PRODUCED BY COEFF2,
C                        THIS ROUTINE EVALUATES THE FUNCTION ON A
C                        SELECTED DERIVATIVE OF ANY POINT WHERE
C                        TWO-DIMENSIONAL INTERPOLATION IS REQUIRED.
C
C USAGE                  R = TERP2 (XB,YB,NX,X,NY,Y,F,FXX,FYY,FXXYY,IDM,
C                                   IXD,IYD)
C
C ARGUMENTS
C
C ON INPUT               XB, YB
C                          VALUES OF THE INDEPENDENT VARIABLES, X AND Y,
C                          AT WHICH INTERPOLATION IS REQUIRED.
C
C                        NX
C                          NUMBER OF GRID POINTS IN THE X-DIRECTION.  NX
C                          MUST BE AT LEAST 4.
C
C                        X
C                          TABLE OF NX VALUES OF THE INDEPENDENT
C                          VARIABLE, X, ARRANGED IN ASCENDING ORDER.
C                          DIMENSION OF X IN THE CALLING PROGRAM MUST BE
C                          AT LEAST NX.
C
C                        NY
C                          NUMBER OF GRID POINTS IN THE Y-DIRECTION.  NY
C                          MUST BE AT LEAST 4.
C
C                        Y
C                          TABLE OF NY VALUES OF THE INDEPENDENT
C                          VARIABLE, Y, ARRANGED IN ASCENDING ORDER.
C                          DIMENSION OF Y IN THE CALLING PROGRAM MUST BE
C                          AT LEAST NY.
C
C                        F
C                          TWO-DIMENSIONAL ARRAY OF FUNCTION VALUES AT
C                          GRID POINTS DEFINED BY THE ARRAYS X AND Y.
C                          DIMENSION OF F IN THE CALLING PROGRAM IS
C                          (IDM,NYY), WHERE
C                              IDM .GE. NX
C                              NYY .GE. NY
C
C                        FXX
C                          ARRAY OF SECOND DERIVATIVES OF F WITH RESPECT
C                          TO X COMPUTED BY COEFF2.  DIMENSION OF FXX IN
C                          THE CALLING PROGRAM IS (IDM,NYY).  SEE UNDER
C                          F ABOVE.
C
C                        FYY
C                          ARRAY OF SECOND DERIVATIVES OF F WITH RESPECT
C                          TO Y COMPUTED BY COEFF2.  DIMENSION OF FYY IN
C                          THE CALLING PROGRAM IS (IDM,NYY).
C
C                        FXXYY
C                          ARRAY OF FOURTH DERIVATIVES,
C                          (D/DX)**2*(D/DY)**2*F, COMPUTED BY COEFF2.
C                          DIMENSION OF FXXYY IN THE CALLING PROGRAM IS
C                          (IDM,NYY).
C
C                        IDM
C                          FIRST DIMENSION IN THE CALLING PROGRAM OF
C                          ARRAYS F, FXX, FYY AND FXXYY,
C                              IDM .GE. NX
C
C                        IXD, IYD
C                          DEFINE DERIVATIVE TO BE RETURNED BY THE
C                          FUNCTION TERP2.  IXD, IYD MAY EACH TAKE THE
C                          THE VALUES 0, 1, 2.  THE DERIVATIVE RETURNED
C                          IS (D/DX)**IXD*(D/DY)**IYD*F.
C                            NOTE THAT IF IXD = IYD = 0, THE FUNCTION
C                            VALUE ITSELF IS RETURNED.
C
C TIMING                 THIS PROCEDURE IS FAST.  THE MAXIMUM
C                        TIME FOR THE BINARY SEARCH IS PROPORTIONAL TO
C                        ALOG(NX*NY).  THE TIME FOR FUNCTION EVALUATION
C                        IS INDEPENDENT OF N.
C***********************************************************************
      SUBROUTINE TERP1 (N,X,F,W,Y,INT,TAB,ITAB)
C
      DIMENSION       X(1)       ,F(1)       ,W(1)       ,TAB(3)     ,
     1                ITAB(3)
      SAVE
C
C PERFORM SEARCH
C
      CALL SEARCH (Y,X,N,I)
C
C CARRY OUT INTERPOLATION (OR EXTRAPOLATION)
C
      CALL INTERP (N,X,F,W,Y,I,INT,TAB,ITAB)
      RETURN
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


program test_get_value

  use bmif_1_2, only: BMI_SUCCESS, BMI_FAILURE
  use bmisnowf
  use fixtures, only: status, print_array

  implicit none

  character (len=*), parameter :: config_file = "test.cfg"
  type (bmi_snow) :: m
  integer :: retcode

  retcode = test1()
  if (retcode.ne.BMI_SUCCESS) then
     stop BMI_FAILURE
  end if

  retcode = test2()
  if (retcode.ne.BMI_SUCCESS) then
     stop BMI_FAILURE
  end if

  retcode = test3()
  if (retcode.ne.BMI_SUCCESS) then
     stop BMI_FAILURE
  end if

contains

  ! Test getting plate_surface__temperature.
  function test1() result(code)
    character (len=*), parameter :: &
         var_name = "land_surface_air__temperature"
    integer, parameter :: size = 1
    real, parameter :: expected = 2.5
    real :: tval(size)
    integer :: i, code

    status = m%initialize(config_file)
    status = m%get_value(var_name, tval)
    status = m%finalize()

    print *, abs(tval(1) - expected)

    ! Visual inspection.
    write(*,*) "Test 1"

    code = BMI_SUCCESS
    do i = 1, size
       if (abs(tval(i) - expected) .gt. 1E-5) then
          code = BMI_FAILURE
          exit
       end if
    end do
  end function test1

  ! Test getting precipitation.
  function test2() result(code)
    character (len=*), parameter :: &
         var_name = "precipitation_mass_flux"
    integer, parameter :: size = 1
    real, parameter :: expected = 0
    real :: tval(size)
    integer :: i, code

    status = m%initialize(config_file)
    status = m%get_value(var_name, tval)
    status = m%finalize()

    print *, abs(tval(1) - expected)

    ! Visual inspection.
    write(*,*) "Test 2"

    code = BMI_SUCCESS
    do i = 1, size
       if (abs(tval(i) - expected) .gt. 1E-5) then
          code = BMI_FAILURE
          exit
       end if
    end do
  end function test2

  ! Test getting snow depth.
  function test3() result(code)
    character (len=*), parameter :: &
         var_name = "snowpack__depth"
    integer, parameter :: size = 1
    real, parameter :: expected = 0
    real :: tval(size)
    integer :: i, code

    status = m%initialize(config_file)
    status = m%get_value(var_name, tval)
    status = m%finalize()

    print *, abs(tval(1) - expected)

    ! Visual inspection.
    write(*,*) "Test 3"

    code = BMI_SUCCESS
    do i = 1, size
       if (abs(tval(i) - expected) .gt. 1E-5) then
          code = BMI_FAILURE
          exit
       end if
    end do
  end function test3

end program test_get_value


text:
      subroutine qzvec(nm,n,a,b,alfr,alfi,beta,z)
c
      integer i,j,k,m,n,en,ii,jj,na,nm,nn,isw,enm2
      double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n),z(nm,n)
      double precision d,q,r,s,t,w,x,y,di,dr,ra,rr,sa,ti,tr,t1,t2,w1,x1,
     x       zz,z1,alfm,almi,almr,betm,epsb
c
c     this subroutine is the optional fourth step of the qz algorithm
c     for solving generalized matrix eigenvalue problems,
c     siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
c
c     this subroutine accepts a pair of real matrices, one of them in
c     quasi-triangular form (in which each 2-by-2 block corresponds to
c     a pair of complex eigenvalues) and the other in upper triangular
c     form.  it computes the eigenvectors of the triangular problem and
c     transforms the results back to the original coordinate system.
c     it is usually preceded by  qzhes,  qzit, and  qzval.
c
c     on input
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c
c        n is the order of the matrices.
c
c        a contains a real upper quasi-triangular matrix.
c
c        b contains a real upper triangular matrix.  in addition,
c          location b(n,1) contains the tolerance quantity (epsb)
c          computed and saved in  qzit.
c
c        alfr, alfi, and beta  are vectors with components whose
c          ratios ((alfr+i*alfi)/beta) are the generalized
c          eigenvalues.  they are usually obtained from  qzval.
c
c        z contains the transformation matrix produced in the
c          reductions by  qzhes,  qzit, and  qzval, if performed.
c          if the eigenvectors of the triangular problem are
c          desired, z must contain the identity matrix.
c
c     on output
c
c        a is unaltered.  its subdiagonal elements provide information
c           about the storage of the complex eigenvectors.
c
c        b has been destroyed.
c
c        alfr, alfi, and beta are unaltered.
c
c        z contains the real and imaginary parts of the eigenvectors.
c          if alfi(i) .eq. 0.0, the i-th eigenvalue is real and
c            the i-th column of z contains its eigenvector.
c          if alfi(i) .ne. 0.0, the i-th eigenvalue is complex.
c            if alfi(i) .gt. 0.0, the eigenvalue is the first of
c              a complex pair and the i-th and (i+1)-th columns
c              of z contain its eigenvector.
c            if alfi(i) .lt. 0.0, the eigenvalue is the second of
c              a complex pair and the (i-1)-th and i-th columns
c              of z contain the conjugate of its eigenvector.
c          each eigenvector is normalized so that the modulus
c          of its largest component is 1.0 .
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      epsb = b(n,1)
      isw = 1
c     .......... for en=n step -1 until 1 do -- ..........
      do 800 nn = 1, n
         en = n + 1 - nn
         na = en - 1
         if (isw .eq. 2) go to 795
         if (alfi(en) .ne. 0.0d0) go to 710
c     .......... real vector ..........
         m = en
         b(en,en) = 1.0d0
         if (na .eq. 0) go to 800
         alfm = alfr(m)
         betm = beta(m)
c     .......... for i=en-1 step -1 until 1 do -- ..........
         do 700 ii = 1, na
            i = en - ii
            w = betm * a(i,i) - alfm * b(i,i)
            r = 0.0d0
c
            do 610 j = m, en
  610       r = r + (betm * a(i,j) - alfm * b(i,j)) * b(j,en)
c
            if (i .eq. 1 .or. isw .eq. 2) go to 630
            if (betm * a(i,i-1) .eq. 0.0d0) go to 630
            zz = w
            s = r
            go to 690
  630       m = i
            if (isw .eq. 2) go to 640
c     .......... real 1-by-1 block ..........
            t = w
            if (w .eq. 0.0d0) t = epsb
            b(i,en) = -r / t
            go to 700
c     .......... real 2-by-2 block ..........
  640       x = betm * a(i,i+1) - alfm * b(i,i+1)
            y = betm * a(i+1,i)
            q = w * zz - x * y
            t = (x * s - zz * r) / q
            b(i,en) = t
            if (dabs(x) .le. dabs(zz)) go to 650
            b(i+1,en) = (-r - w * t) / x
            go to 690
  650       b(i+1,en) = (-s - y * t) / zz
  690       isw = 3 - isw
  700    continue
c     .......... end real vector ..........
         go to 800
c     .......... complex vector ..........
  710    m = na
         almr = alfr(m)
         almi = alfi(m)
         betm = beta(m)
c     .......... last vector component chosen imaginary so that
c                eigenvector matrix is triangular ..........
         y = betm * a(en,na)
         b(na,na) = -almi * b(en,en) / y
         b(na,en) = (almr * b(en,en) - betm * a(en,en)) / y
         b(en,na) = 0.0d0
         b(en,en) = 1.0d0
         enm2 = na - 1
         if (enm2 .eq. 0) go to 795
c     .......... for i=en-2 step -1 until 1 do -- ..........
         do 790 ii = 1, enm2
            i = na - ii
            w = betm * a(i,i) - almr * b(i,i)
            w1 = -almi * b(i,i)
            ra = 0.0d0
            sa = 0.0d0
c
            do 760 j = m, en
               x = betm * a(i,j) - almr * b(i,j)
               x1 = -almi * b(i,j)
               ra = ra + x * b(j,na) - x1 * b(j,en)
               sa = sa + x * b(j,en) + x1 * b(j,na)
  760       continue
c
            if (i .eq. 1 .or. isw .eq. 2) go to 770
            if (betm * a(i,i-1) .eq. 0.0d0) go to 770
            zz = w
            z1 = w1
            r = ra
            s = sa
            isw = 2
            go to 790
  770       m = i
            if (isw .eq. 2) go to 780
c     .......... complex 1-by-1 block ..........
            tr = -ra
            ti = -sa
  773       dr = w
            di = w1
c     .......... complex divide (t1,t2) = (tr,ti) / (dr,di) ..........
  775       if (dabs(di) .gt. dabs(dr)) go to 777
            rr = di / dr
            d = dr + di * rr
            t1 = (tr + ti * rr) / d
            t2 = (ti - tr * rr) / d
            go to (787,782), isw
  777       rr = dr / di
            d = dr * rr + di
            t1 = (tr * rr + ti) / d
            t2 = (ti * rr - tr) / d
            go to (787,782), isw
c     .......... complex 2-by-2 block ..........
  780       x = betm * a(i,i+1) - almr * b(i,i+1)
            x1 = -almi * b(i,i+1)
            y = betm * a(i+1,i)
            tr = y * ra - w * r + w1 * s
            ti = y * sa - w * s - w1 * r
            dr = w * zz - w1 * z1 - x * y
            di = w * z1 + w1 * zz - x1 * y
            if (dr .eq. 0.0d0 .and. di .eq. 0.0d0) dr = epsb
            go to 775
  782       b(i+1,na) = t1
            b(i+1,en) = t2
            isw = 1
            if (dabs(y) .gt. dabs(w) + dabs(w1)) go to 785
            tr = -ra - x * b(i+1,na) + x1 * b(i+1,en)
            ti = -sa - x * b(i+1,en) - x1 * b(i+1,na)
            go to 773
  785       t1 = (-r - zz * b(i+1,na) + z1 * b(i+1,en)) / y
            t2 = (-s - zz * b(i+1,en) - z1 * b(i+1,na)) / y
  787       b(i,na) = t1
            b(i,en) = t2
  790    continue
c     .......... end complex vector ..........
  795    isw = 3 - isw
  800 continue
c     .......... end back substitution.
c                transform to original coordinate system.
c                for j=n step -1 until 1 do -- ..........
      do 880 jj = 1, n
         j = n + 1 - jj
c
         do 880 i = 1, n
            zz = 0.0d0
c
            do 860 k = 1, j
  860       zz = zz + z(i,k) * b(k,j)
c
            z(i,j) = zz
  880 continue
c     .......... normalize so that modulus of largest
c                component of each vector is 1.
c                (isw is 1 initially from before) ..........
      do 950 j = 1, n
         d = 0.0d0
         if (isw .eq. 2) go to 920
         if (alfi(j) .ne. 0.0d0) go to 945
c
         do 890 i = 1, n
            if (dabs(z(i,j)) .gt. d) d = dabs(z(i,j))
  890    continue
c
         do 900 i = 1, n
  900    z(i,j) = z(i,j) / d
c
         go to 950
c
  920    do 930 i = 1, n
            r = dabs(z(i,j-1)) + dabs(z(i,j))
            if (r .ne. 0.0d0) r = r * dsqrt((z(i,j-1)/r)**2
     x                                     +(z(i,j)/r)**2)
            if (r .gt. d) d = r
  930    continue
c
         do 940 i = 1, n
            z(i,j-1) = z(i,j-1) / d
            z(i,j) = z(i,j) / d
  940    continue
c
  945    isw = 3 - isw
  950 continue
c
      return
      end


text:
C
C**************************  VNOD3D ************************************
C
C  average Darcy velocitities on each tetrahedra to obtain the Darcy
C  velocities at each node
C
C***********************************************************************
C
      SUBROUTINE VNOD3D(N,NT,TP,TETRA,UU,VV,WW,UNOD,VNOD,WNOD)
C
      IMPLICIT  NONE
      INTEGER   I,K,II,INOD
      INTEGER   N,NT
      INTEGER   TP(*),TETRA(5,*)
      REAL*8    UU(*),VV(*),WW(*),UNOD(*),VNOD(*),WNOD(*)
C
      CALL INIT0R(N,UNOD)
      CALL INIT0R(N,VNOD)
      CALL INIT0R(N,WNOD)
      DO K=1,NT
         DO II=1,4
            INOD=TETRA(II,K)
            UNOD(INOD)=UNOD(INOD) + UU(K)
            VNOD(INOD)=VNOD(INOD) + VV(K)
            WNOD(INOD)=WNOD(INOD) + WW(K)
         END DO
      END DO
      DO I=1,N
         UNOD(I)=UNOD(I)/TP(I)
         VNOD(I)=VNOD(I)/TP(I)
         WNOD(I)=WNOD(I)/TP(I)
      END DO
C
      RETURN
      END


! { dg-do run }
! { dg-additional-options "-fcoarray=single -fdump-tree-original" }
!
! Check that {L,U}{,CO}BOUND intrinsics are properly simplified.
!
  type :: t
    integer :: c
  end type t

  type(t) :: d(3:8) = t(7)
  type(t) :: e[5:9,-1:*]

  if (lbound(d,      1) /= 3) call abort
  if (lbound(d(3:5), 1) /= 1) call abort
  if (lbound(d%c,    1) /= 1) call abort
  if (ubound(d,      1) /= 8) call abort
  if (ubound(d(3:5), 1) /= 3) call abort
  if (ubound(d%c,    1) /= 6) call abort

  if (lcobound(e,   1) /=  5) call abort
  if (lcobound(e%c, 1) /=  5) call abort
  if (lcobound(e,   2) /= -1) call abort
  if (lcobound(e%c, 2) /= -1) call abort
  if (ucobound(e,   1) /=  9) call abort
  if (ucobound(e%c, 1) /=  9) call abort
  ! no simplification for ucobound(e{,%c}, dim=2)
end
! { dg-final { scan-tree-dump-not "bound" "original" } }
! { dg-final { scan-tree-dump-not "abort" "original" } }


text:
           program demo_iostat
             implicit none
             integer :: stat, i
             open(88, file='test.dat')
             read(88, *, iostat=stat) i
             if(is_iostat_end(stat)) stop 'end of file'
           end program demo_iostat

max_stars_repo_name: urbanjost/fortran-intrinsic-manpages

max_issues_repo_name: urbanjost/fortran-intrinsic-manpages

max_forks_repo_name: urbanjost/fortran-intrinsic-manpages


text:
      SUBROUTINE INTINH (ASFTAB, START, END, PERM)

C  INTINH creates the structure network for testing interior ASFs.

C  Input parameters:
C    ASFTAB    : table containing ASF values
C    START,END : first,last column in the table (corresponding to aspect
C                identifier) from which values are taken.
C  Output parameters:
C    PERM      : permutation of locations used for structure network

C aspect identifier
      INTEGER    PLN,      PEDFG,    PEDCI
      PARAMETER (PLN  = 0, PEDFG=14, PEDCI=17)

C aspect source
C                bundled     individual
      INTEGER    PBUNDL,     PINDIV
      PARAMETER (PBUNDL = 0, PINDIV = 1)

C composition type
C                preconcatenate  postconcatenate  replace
      INTEGER    PCPRE,          PCPOST,          PCREPL
      PARAMETER (PCPRE = 0,      PCPOST = 1,      PCREPL = 2)

C off/on switch for edge flag and error handling mode
      INTEGER    POFF,     PON
      PARAMETER (POFF = 0, PON = 1)

      INTEGER    START,END, ASFTAB(5, PLN:PEDCI), PERM(14), ERRIND
      INTEGER    IX, EXPDX(14), SIZ

      REAL       XACT(4),XEXP(4),XFORM(3,3), YLOCEL

C  x-location of actual/expected interiors
      CALL SETRVS ('0.2,0.2,0.5,0.5', XACT, SIZ)
      CALL SETRVS ('0.6,0.6,0.9,0.9', XEXP, SIZ)

C  randomize order of interiors
      CALL RN1SHF (14, PERM)

C  set edge flag = OFF
      CALL PSIASF (PEDFG, PINDIV)
      CALL PSEDFG (POFF)

C  set up CSS as described above
      CALL LOCINT (PERM(1), XACT)
      CALL PEXST  (102)
      CALL LOCINT (PERM(9), XACT)
      CALL ASFSET (2, ASFTAB, START,END)
      CALL LOCINT (PERM(10), XACT)
      CALL XFINH  (PERM)

C  execute 104
      CALL PEXST (104)
C  now, cancel out transformation ...
      CALL IDMAT (3, XFORM)
      CALL PSLMT (XFORM, PCREPL)

      CALL LOCINT (PERM(13), XACT)
      CALL LOCINT (PERM(14), XACT)
      CALL PEXST (105)
      CALL PCLST

C  structure #102
      CALL POPST (102)
      CALL LOCINT (PERM(2), XACT)
      CALL ASFSET (3, ASFTAB, START,END)
      CALL PEXST (103)
      CALL LOCINT (PERM(8), XACT)
      CALL PCLST

C  structure #103
      CALL POPST (103)
      CALL LOCINT (PERM(3), XACT)
      CALL ASFSET (4, ASFTAB, START,END)
      CALL LOCINT (PERM(4), XACT)
      CALL PEXST (104)
      CALL LOCINT (PERM(7), XACT)
      CALL PCLST

C  structure #104
      CALL POPST (104)
C  polymarker 5 / 11
      CALL LOCINT (PERM(5), XACT)
      CALL ASFSET (5, ASFTAB, START,END)
C  polymarker 6 / 12
      CALL LOCINT (PERM(6), XACT)
      CALL PCLST

C  Expected attributes: structure #105
      CALL POPST (105)
C  values for expected ASFs:
      CALL SETVAL ('1,1,3,4,4,5,4,3,1,2,2,5,2,4', EXPDX)
      DO 400 IX =1,14
         CALL ASFSET (EXPDX(IX), ASFTAB, START,END)
         CALL LOCINT (PERM(IX), XEXP)
400   CONTINUE

C  draw labels
      CALL NUMLAB (14, 0.15, YLOCEL(1), 1.0/15)
      CALL PCLST

      END