text:
C   18/01/93 303061342  MEMBER NAME  UNDULATOR *.FORT     M  E2FORT
        subroutine undulator(np,x,px,y,py,z,g,dv,sx,sy,sz,ulist)
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                 c
c  Undulator tracking subroutine.                                 c
c  K.Ohmi and E.Forest        KEK Report 92-14                    c
c                                                                 c
c                           written by K.Ohmi                     c
c                                                                 c
c                                        Dec  1994                c
c                                                                 c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      use tmacro
      use ffs_flag, only:calpol
      implicit real*8 (a-h,o-z)
c      implicit none
      parameter (nsli=1,npole=2,len=3,i_F0=4,i_G0=5,i_dph=6,
     &   i_Kz=7,i_Kx=8,i_Ky=9,i_Qx=10,i_Qy=11,lambda=12,i_ds=13)
      real*8 x(np),px(np),y(np),py(np),z(np),g(np),dv(np),
     $     sx(np),sy(np),sz(np)
      real*8 ulist(*)
      real*8 sinkz,sinkzp
      real*8 kx1,kx2,kx3,kx4,ff,gg,fx,fy,gx,gy,dpds
      real*8 det1a,fgx,fgy,tmp
      integer i,nslice
c
      nslice=ulist(1)+0.01
c
c   Lorentz factor gamma is stored in h0.
c   Factor gammabeta is stored in p0.
c     gamma=h0
c     gammabet=p0
c
c
c     Initialize dx/dy/theta(cost/sint) for inc/T(ENT|EXIT).inc
      dx=0.d0
      dy=0.d0
      theta=0.d0
c
      include 'inc/TENT.inc'
      do 10 i=1,np
c
c   Get Canonical variables
c       g(i)=pn-1.d0   Not good for precision
c       g(i)=g(i)*(g(i)+2.d0)
c   pn=|p|/p0=1+delta
       pn=1.d0+g(i)
       px(i)=px(i)*pn
       py(i)=py(i)*pn
 10    continue

      do 11 i=1,np
       do 21 j=1,nslice
         kx1=ulist(i_Kx)*x(i)
         kx2=ulist(i_Ky)*y(i)
         kx3=ulist(i_Qx)*x(i)
         kx4=ulist(i_Qy)*y(i)
         sinkz=sin(ulist(i_Kz)*j*ulist(i_ds))
         sinkzp=sin(ulist(i_Kz)*j*ulist(i_ds)+ulist(i_dph))
         ff=ulist(i_F0)/ulist(i_Kz)*cos(kx1)*cosh(kx2)*sinkz
     &       -ulist(i_Qy)*ulist(i_G0)/ulist(i_Qx)/ulist(i_Kz)*
     &        sinh(kx3)*sin(kx4)*sinkzp
         gg=-ulist(i_G0)/ulist(i_Kz)*cos(kx4)*cosh(kx3)*sinkzp
     &       +ulist(i_Kx)*ulist(i_F0)/ulist(i_Ky)/ulist(i_Kz)*
     &        sinh(kx2)*sin(kx1)*sinkz
c
         fx=-ulist(i_Kx)*ulist(i_F0)/
     &        ulist(i_Kz)*sin(kx1)*cosh(kx2)*sinkz
     &       -ulist(i_Qy)*ulist(i_G0)/
     &        ulist(i_Kz)*cosh(kx3)*sin(kx4)*sinkzp
         fy=ulist(i_Ky)*ulist(i_F0)/
     &        ulist(i_Kz)*cos(kx1)*sinh(kx2)*sinkz
     &       -ulist(i_Qy)*ulist(i_Qy)*ulist(i_G0)/ulist(i_Qx)/
     &        ulist(i_Kz)*sinh(kx3)*cos(kx4)*sinkzp
c
         gx=-ulist(i_G0)*ulist(i_Qx)/
     &        ulist(i_Kz)*cos(kx4)*sinh(kx3)*sinkzp
     &       +ulist(i_Kx)*ulist(i_Kx)*ulist(i_F0)/ulist(i_Ky)/
     &        ulist(i_Kz)*sinh(kx2)*cos(kx1)*sinkz
         gy=ulist(i_G0)*ulist(i_Qy)/
     &        ulist(i_Kz)*sin(kx4)*cosh(kx3)*sinkzp
     &       +ulist(i_Kx)*ulist(i_F0)/
     &        ulist(i_Kz)*cosh(kx2)*sin(kx1)*sinkz
         dpds=ulist(i_ds)/(1.+g(i))
c
         fx=dpds*fx
         gx=dpds*gx
         fy=dpds*fy
         gy=dpds*gy
         det1a=1.-fx-gy+fx*gy-fy*gx
c
         fgx=dpds*(ff*fx+gg*gx)
         fgy=dpds*(ff*fy+gg*gy)
c
         tmp=fx
         fx=(1.-gy)/det1a
         gx=gx/det1a
         fy=fy/det1a
         gy=(1.-tmp)/det1a
c
         tmp=px(i)-fgx
         px(i)=fx*tmp+gx*(py(i)-fgy)
         py(i)=fy*tmp+gy*(py(i)-fgy)
c
         fgx=(px(i)-ff)
         fgy=(py(i)-gg)
         x(i)=x(i)+dpds*fgx
         y(i)=y(i)+dpds*fgy
         z(i)=z(i)-0.5*(fgx*fgx+fgy*fgy)*(dpds*dpds)/ulist(i_ds)
 21    continue
 11    continue
c
c   Return to SAD variables
c
c   pn=|p|/p0=1+delta
      do 19 i=1,np
      pn=1.d0+g(i)
      px(i)=px(i)/pn
      py(i)=py(i)/pn
c   When you change momentum, change dv(i) as
      h1=sqrt(1.d0+(pn*p0)**2)
      dv(i)=-g(i)*(1.d0+pn)/h1/(h1+pn*h0)+dvfs
c   g(i)=sqrt(pn)-1.d0  Not good for precision
c      g(i)=g(i)/(1.d0+sqrt(pn))
 19   continue

      include 'inc/TEXIT.inc'
      return
      end

      subroutine undinit(p_in,ulist)
      use tmacro
      implicit real*8 (a-h,o-z)
      parameter (nsli=1,npole=2,len=3,i_F0=4,i_G0=5,i_dph=6,
     &   i_Kz=7,i_Kx=8,i_Ky=9,i_Qx=10,i_Qy=11,lambda=12,i_ds=13)
      parameter (nulist=20)
      real*8 p_in(70),ulist(nulist)
c
c  p_in(1) length , (2) Bx     , (3) By     , (4) Kx , (5) Qy ,
c      (6) dphase , (7) slice  , (8) Pole
c
      ulist(nsli)=p_in(7)
      ulist(npole)=p_in(8)
      ulist(len)=p_in(1)
      ulist(i_F0)=p_in(2)/pgev*c
      ulist(i_G0)=p_in(3)/pgev*c
      ulist(i_dph)=p_in(6)*pi/180.
c
      ulist(i_Kz)=pi2*ulist(npole)/ulist(len)
      ulist(i_Kx)=p_in(4)
      ulist(i_Ky)=sqrt(ulist(i_Kz)*ulist(i_Kz)
     &            +ulist(i_Kx)*ulist(i_Kx))
      ulist(i_Qy)=p_in(5)
      ulist(i_Qx)=sqrt(ulist(i_Kz)*ulist(i_Kz)
     &            +ulist(i_Qy)*ulist(i_Qy))
      ulist(lambda)=ulist(len)/ulist(npole)
c
      ulist(i_ds)=ulist(len)/ulist(nsli)
c
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
Frank Dogg is a cool guy, hes really into the hiphop and dj scene here in davis.

He runs a recording studio, but i see him do more house parties. Get to know him, and you might get to be in his tv show.

Ive seen him spin at little prague and Uncle Vitos Slice of N.Y. Uncle Vitos and I always see him in traffic in his frankdoggmobile van.

he also has a tv show on youtube at

http://www.youtube.com/frankdoggtv

http://farm4.static.flickr.com/3430/3278589244_fd9387008e.jpg

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      INTEGER FUNCTION MB03ND( N, THETA, Q2, E2, PIVMIN, INFO )
C
C     SLICOT RELEASE 5.5.
C
C     Copyright (c) 2002-2012 NICONET e.V.
C
C     PURPOSE
C
C     To find the number of singular values of the bidiagonal matrix
C
C              |q(1) e(1)  .    ...    0   |
C              | 0   q(2) e(2)         .   |
C          J = | .                     .   |
C              | .                   e(N-1)|
C              | 0   ...     ...   0  q(N) |
C
C     which are less than or equal to a given bound THETA.
C
C     This routine is intended to be called only by other SLICOT
C     routines.
C
C     ARGUMENTS
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the bidiagonal matrix J.  N >= 0.
C
C     THETA   (input) DOUBLE PRECISION
C             Given bound.
C             Note: If THETA < 0.0 on entry, then MB03ND is set to 0
C                   as the singular values of J are non-negative.
C
C     Q2      (input) DOUBLE PRECISION array, dimension (N)
C             This array must contain the squares of the diagonal
C             elements q(1),q(2),...,q(N) of the bidiagonal matrix J.
C             That is, Q2(i) = J(i,i)**2 for i = 1,2,...,N.
C
C     E2      (input) DOUBLE PRECISION array, dimension (N-1)
C             This array must contain the squares of the superdiagonal
C             elements e(1),e(2),...,e(N-1) of the bidiagonal matrix J.
C             That is, E2(k) = J(k,k+1)**2 for k = 1,2,...,N-1.
C
C     PIVMIN  (input) DOUBLE PRECISION
C             The minimum absolute value of a "pivot" in the Sturm
C             sequence loop.
C             PIVMIN >= max( max( |q(i)|, |e(k)| )**2*sf_min, sf_min ),
C             where i = 1,2,...,N, k = 1,2,...,N-1, and sf_min is at
C             least the smallest number that can divide one without
C             overflow (see LAPACK Library routine DLAMCH).
C             Note that this condition is not checked by the routine.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value.
C
C     METHOD
C
C     The computation of the number of singular values s(i) of J which
C     are less than or equal to THETA is based on applying Sylvester's
C     Law of Inertia, or equivalently, Sturm sequences [1,p.52] to the
C     unreduced symmetric tridiagonal matrices associated with J as
C     follows. Let T be the following 2N-by-2N symmetric matrix
C     associated with J:
C
C               | 0   J'|
C          T =  |       |.
C               | J   0 |
C
C     (The eigenvalues of T are given by s(1),s(2),...,s(N),-s(1),-s(2),
C     ...,-s(N)). Then, by permuting the rows and columns of T into the
C     order 1, N+1, 2, N+2, ..., N, 2N it follows that T is orthogonally
C     similar to the tridiagonal matrix T" with zeros on its diagonal
C     and q(1), e(1), q(2), e(2), ..., e(N-1), q(N) on its offdiagonals
C     [3,4]. If q(1),q(2),...,q(N) and e(1),e(2),...,e(N-1) are nonzero,
C     Sylvester's Law of Inertia may be applied directly to T".
C     Otherwise, T" is block diagonal and each diagonal block (which is
C     then unreduced) must be analysed separately by applying
C     Sylvester's Law of Inertia.
C
C     REFERENCES
C
C     [1] Parlett, B.N.
C         The Symmetric Eigenvalue Problem.
C         Prentice Hall, Englewood Cliffs, New Jersey, 1980.
C
C     [2] Demmel, J. and Kahan, W.
C         Computing Small Singular Values of Bidiagonal Matrices with
C         Guaranteed High Relative Accuracy.
C         Technical Report, Courant Inst., New York, March 1988.
C
C     [3] Van Huffel, S. and Vandewalle, J.
C         The Partial Total Least-Squares Algorithm.
C         J. Comput. and Appl. Math., 21, pp. 333-341, 1988.
C
C     [4] Golub, G.H. and Kahan, W.
C         Calculating the Singular Values and Pseudo-inverse of a
C         Matrix.
C         SIAM J. Numer. Anal., Ser. B, 2, pp. 205-224, 1965.
C
C     [5] Demmel, J.W., Dhillon, I. and Ren, H.
C         On the Correctness of Parallel Bisection in Floating Point.
C         Computer Science Division Technical Report UCB//CSD-94-805,
C         University of California, Berkeley, CA 94720, March 1994.
C
C     NUMERICAL ASPECTS
C
C     The singular values s(i) could also be obtained with the use of
C     the symmetric tridiagonal matrix T = J'J, whose eigenvalues are
C     the squared singular values of J [4,p.213]. However, the method
C     actually used by the routine is more accurate and equally
C     efficient (see [2]).
C
C     To avoid overflow, matrix J should be scaled so that its largest
C     element is no greater than  overflow**(1/2) * underflow**(1/4)
C     in absolute value (and not much smaller than that, for maximal
C     accuracy).
C
C     With respect to accuracy the following condition holds (see [2]):
C
C     If the established value is denoted by p, then at least p
C     singular values of J are less than or equal to
C     THETA/(1 - (3 x N - 1.5) x EPS) and no more than p singular values
C     are less than or equal to
C     THETA x (1 - (6 x N-2) x EPS)/(1 - (3 x N - 1.5) x EPS).
C
C     CONTRIBUTOR
C
C     Release 3.0: V. Sima, Katholieke Univ. Leuven, Belgium, Apr. 1997.
C     Supersedes Release 2.0 routine MB03BD by S. Van Huffel, Katholieke
C     University, Leuven, Belgium.
C
C     REVISIONS
C
C     July 10, 1997.
C
C     KEYWORDS
C
C     Bidiagonal matrix, singular values.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO
      PARAMETER         ( ZERO = 0.0D0 )
C     .. Scalar Arguments ..
      INTEGER           INFO, N
      DOUBLE PRECISION  PIVMIN, THETA
C     .. Array Arguments ..
      DOUBLE PRECISION  E2(*), Q2(*)
C     .. Local Scalars ..
      INTEGER           J, NUMEIG
      DOUBLE PRECISION  R, T
C     .. External Subroutines ..
      EXTERNAL          XERBLA
C     .. Intrinsic Functions ..
      INTRINSIC         ABS
C     .. Executable Statements ..
C
C     Test the input scalar arguments.  PIVMIN is not checked.
C
      INFO = 0
      IF( N.LT.0 ) THEN
         INFO = -1
C
C        Error return.
C
         CALL XERBLA( 'MB03ND', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      IF ( N.EQ.0 .OR. THETA.LT.ZERO ) THEN
         MB03ND = 0
         RETURN
      END IF
C
      NUMEIG = N
      T = -THETA
      R = T
      IF ( ABS( R ).LT.PIVMIN ) R = -PIVMIN
C
      DO 20 J = 1, N - 1
         R = T - Q2(J)/R
         IF ( ABS( R ).LT.PIVMIN ) R = -PIVMIN
         IF ( R.GT.ZERO ) NUMEIG = NUMEIG - 1
         R = T - E2(J)/R
         IF ( ABS( R ).LT.PIVMIN ) R = -PIVMIN
         IF ( R.GT.ZERO ) NUMEIG = NUMEIG - 1
   20 CONTINUE
C
      R = T - Q2(N)/R
      IF ( ABS( R ).LT.PIVMIN ) R = -PIVMIN
      IF ( R.GT.ZERO ) NUMEIG = NUMEIG - 1
      MB03ND = NUMEIG
C
      RETURN
C *** Last line of MB03ND ***
      END

max_stars_repo_name: bgin/MissileSimulation

max_issues_repo_name: bgin/MissileSimulation

max_forks_repo_name: bgin/MissileSimulation


program main

use test, only: test_value
use io, only: stdout

use map

implicit none

! utmgeo
! stereo

write(stdout,*) 'map_module unit test passed'
end


module mod_sort
    use omp_lib
    use m_mrgrnk

    implicit none

    integer, parameter :: max_simple_sort_size = 20
    private
    public :: parallel_sort
contains

    subroutine parallel_sort(A, order)
        real, intent(in), dimension(:) :: A
        integer, intent(out), dimension(size(A)) :: order

        integer :: len, from, middle, to, nthreads, thread, chunk, chunk2, i

        len = size(A)
        nthreads = omp_get_max_threads()
        chunk = len/nthreads

        !----------------------------------------
        ! Initialize order
        !----------------------------------------
        do i = 1, len
            order(i) = i
        end do

        !----------------------------------------
        ! Sort each chunk
        !----------------------------------------
        !$OMP parallel do default(none)          &
        !$OMP firstprivate(chunk, len, nthreads) &
        !$OMP shared(A, order) private(from, to)
        do thread = 0, nthreads
            from = thread*chunk + 1
            to = min((thread + 1)*chunk, len)

            call mrgrnk(A(from:to), order(from:to))
            order(from:to) = order(from:to) + from - 1

        end do
        !$OMP end parallel do

        !----------------------------------------
        ! Merge pieces together
        !----------------------------------------
        i = 1
        chunk2 = chunk
        do while (chunk2 < size(A))

            !$OMP parallel do default(none)  &
            !$OMP shared(chunk2, A, order)   &
            !$OMP private(from, middle, to)
            do thread = 0, ceiling(.5*size(A)/chunk2)
                from = thread*2*chunk2 + 1
                middle = (thread*2 + 1)*chunk2
                to = (thread*2 + 2)*chunk2

                middle = min(middle, size(A))
                to = min(to, size(A))
                if (from < to) then
                    call merge(A, order, from, middle, to)
                end if
            end do
            !$OMP end parallel do
            chunk2 = chunk2*2
            i = i + 1
        end do
    end subroutine parallel_sort

    !> Merge two parts of A, ordered by order from left to right
  !! around middle.
    subroutine merge(A, order, left, middle, right)
        real, intent(in), dimension(:) :: A
        integer, intent(inout), dimension(size(A)) :: order

        integer, intent(in) :: left, middle, right

        integer :: leftA, rightA, leftB, rightB
        integer :: iA, iB, i
        integer :: lenA, lenB

        integer, dimension(left:middle) :: orderA
        integer, dimension(middle + 1:right) :: orderB

        ! copy order
        orderA = order(left:middle)
        orderB = order(middle + 1:right)

        ! more explicit variables
        leftA = left
        rightA = middle
        leftB = middle + 1
        rightB = right

        ! initialize iA, iB to their leftmost position
        iA = leftA
        iB = leftB

        i = leftA

        do while ((iA <= rightA) .and. (iB <= rightB))
            if (A(orderA(iA)) <= A(orderB(iB))) then
                order(i) = orderA(iA)
                iA = iA + 1
            else
                order(i) = orderB(iB)
                iB = iB + 1
            end if

            i = i + 1
        end do

        ! either A or B still have elements, append them to the new order
        do while (iA <= rightA)
            order(i) = orderA(iA)
            iA = iA + 1

            i = i + 1

        end do
        do while (iB <= rightB)
            order(i) = orderB(iB)
            iB = iB + 1

            i = i + 1
        end do

    end subroutine merge

end module mod_sort


text:
subroutine MPI_Main
  implicit none
  include 'mpif.h'

  integer :: myrank, ierr, numranks, req
  double precision :: inval, outval, expect

  call MPI_Init(ierr)
  call MPI_Comm_rank(MPI_COMM_WORLD, myrank, ierr)
  call MPI_Comm_size(MPI_COMM_WORLD, numranks, ierr)

  inval = myrank + 1
  expect = (numranks*(numranks+1))/2

  call MPI_Allreduce(inval, outval, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
                     MPI_COMM_WORLD, ierr)

  if (myrank .eq. 0) then
    if (outval .eq. expect) then
      print *, 'MPI_Allreduce test passes'
    else
      print *, 'MPI_Allreduce test failed'
    end if
  end if

  call MPI_Iallreduce(inval, outval, 1, MPI_DOUBLE_PRECISION, MPI_SUM, &
                      MPI_COMM_WORLD, req, ierr)

  if (myrank .eq. 0) then
    call MPI_Wait(req, MPI_STATUS_IGNORE, ierr)
    if (outval .eq. expect) then
      print *, 'MPI_Iallreduce test passes'
    else
      print *, 'MPI_Iallreduce test failed'
    end if
  end if

  call MPI_Finalize(ierr)

end subroutine

max_stars_repo_name: scottkwarren/config-db

max_issues_repo_name: scottkwarren/config-db

max_forks_repo_name: scottkwarren/config-db


module LATIN_indicator

  use Structure
  use Matrix_op
  use global_stage
  use local_stage
  !use local_stage_vec

  implicit none

  contains

  function indicator1(G,L)  result(ind)

    !type(Input), intent(inout) :: I
    type(global), intent(inout) :: G
    type(local), intent(inout) :: L

    double precision :: ind
    ind =  100*( norm2( L%d_Ep - G%d_Ep )/norm2( L%d_Ep ) + norm2( L%sigma - G%sigma )/norm2( L%sigma ) )

  end function indicator1

end module LATIN_indicator


text:
! { dg-do compile }
! Tests the fix for PR32827, in which IMPORT :: my_type put the
! symbol into the interface namespace, thereby generating an error
! when the declaration of 'x' is compiled.
!
! Contributed by Douglas Wells <sysmaint@contek.com>
!
subroutine func1(param)
  type :: my_type
    integer :: data
  end type my_type
  type(my_type) :: param
  param%data = 99
end subroutine func1

subroutine func2(param)
  type :: my_type
    integer :: data
  end type my_type
  type(my_type) :: param
  param%data = 21
end subroutine func2

  type :: my_type
    integer :: data
  end type my_type

  interface
    subroutine func1(param)
      import :: my_type
      type(my_type) :: param
    end subroutine func1
  end interface
  interface
    subroutine func2(param)
      import
      type(my_type) :: param
    end subroutine func2
  end interface

  type(my_type) :: x
  call func1(x)
  print *, x%data
  call func2(x)
  print *, x%data
end

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


!> @file
!! \brief The driver program to solve a linear system with default options.
!!
!! <pre>
!! -- Distributed SuperLU routine (version 3.2) --
!! Lawrence Berkeley National Lab, Univ. of California Berkeley.
!! October, 2012
!! </pre>
!
      program f_pzdrive
!
! Purpose
! =======
!
! The driver program F_PDDRIVE.
!
! This example illustrates how to use F_PDGSSVX with the full
! (default) options to solve a linear system.
!
! Seven basic steps are required:
!   1. Create C structures used in SuperLU_DIST
!   2. Initialize the MPI environment and the SuperLU process grid
!   3. Set up the input matrix and the right-hand side
!   4. Set the options argument
!   5. Call f_pdgssvx
!   6. Release the process grid and terminate the MPI environment
!   7. Release all structures
!
!
      use superlu_mod
!      implicit none
      include 'mpif.h'
      integer maxn, maxnz, maxnrhs
      parameter ( maxn = 10000, maxnz = 100000, maxnrhs = 10 )
      integer rowind(maxnz), colptr(maxn)
      double complex  values(maxnz), b(maxn), berr(maxnrhs)
      integer n, m, nnz, nprow, npcol, ldb, init
      integer*4 iam, info, i, ierr, ldb4, nrhs

      integer(superlu_ptr) :: grid
      integer(superlu_ptr) :: options
      integer(superlu_ptr) :: ScalePermstruct
      integer(superlu_ptr) :: LUstruct
      integer(superlu_ptr) :: SOLVEstruct
      integer(superlu_ptr) :: A
      integer(superlu_ptr) :: stat

! Initialize MPI environment
      call mpi_init(ierr)

! Check malloc
!      call f_check_malloc(iam)

! Create Fortran handles for the C structures used in SuperLU_DIST
      call f_create_gridinfo_handle(grid)
      call f_create_options_handle(options)
      call f_create_ScalePerm_handle(ScalePermstruct)
      call f_create_LUstruct_handle(LUstruct)
      call f_create_SOLVEstruct_handle(SOLVEstruct)
      call f_create_SuperMatrix_handle(A)
      call f_create_SuperLUStat_handle(stat)

! Initialize the SuperLU_DIST process grid
      nprow = 2
      npcol = 2
      call f_superlu_gridinit(MPI_COMM_WORLD, nprow, npcol, grid)

! Bail out if I do not belong in the grid.
      call get_GridInfo(grid, iam=iam)
      if ( iam >= nprow * npcol ) then
         go to 100
      endif
      if ( iam == 0 ) then
         write(*,*) ' Process grid ', nprow, ' X ', npcol
      endif

! Read Harwell-Boeing matrix, and adjust the pointers and indices
! to 0-based indexing, as required by C routines.
      if ( iam == 0 ) then
         open(file = "../EXAMPLE/cg20.cua", status = "old", unit = 5)
         call zhbcode1(m, n, nnz, values, rowind, colptr)
         close(unit = 5)
!
         do i = 1, n+1
            colptr(i) = colptr(i) - 1
         enddo
         do i = 1, nnz
            rowind(i) = rowind(i) - 1
         enddo
      endif

! Distribute the matrix to the process gird
      call  f_zcreate_dist_matrix(A, m, n, nnz, values, rowind, colptr, grid)

! Setup the right hand side
      call get_CompRowLoc_Matrix(A, nrow_loc=ldb)
      do i = 1, ldb
         b(i) = 1.0
      enddo
      nrhs = 1
      ldb4 = ldb

! Set the default input options
      call f_set_default_options(options)

! Change one or more options
!      call set_superlu_options(options,Fact=FACTORED)
!      call set_superlu_options(options,ParSymbFact=YES)

! Initialize ScalePermstruct and LUstruct
      call get_SuperMatrix(A, nrow=m, ncol=n)
      call f_ScalePermstructInit(m, n, ScalePermstruct)
      call f_LUstructInit(m, n, LUstruct)

! Initialize the statistics variables
      call f_PStatInit(stat)

! Call the linear equation solver
      call f_pzgssvx(options, A, ScalePermstruct, b, ldb4, nrhs, &
                     grid, LUstruct, SOLVEstruct, berr, stat, info)

      if (info == 0) then
         write (*,*) 'Backward error: ', (berr(i), i = 1, nrhs)
      else
         write(*,*) 'INFO from f_pdgssvx = ', info
      endif

! Deallocate the storage allocated by SuperLU_DIST
      call f_PStatFree(stat)
      call f_Destroy_CompRowLoc_Mat_dist(A)
      call f_ScalePermstructFree(ScalePermstruct)
      call f_Destroy_LU(n, grid, LUstruct)
      call f_LUstructFree(LUstruct)
      call get_superlu_options(options, SolveInitialized=init)
      if (init == YES) then
         call f_zSolveFinalize(options, SOLVEstruct)
      endif

! Release the SuperLU process grid
100   call f_superlu_gridexit(grid)

! Deallocate the C structures pointed to by the Fortran handles
      call f_destroy_gridinfo_handle(grid)
      call f_destroy_options_handle(options)
      call f_destroy_ScalePerm_handle(ScalePermstruct)
      call f_destroy_LUstruct_handle(LUstruct)
      call f_destroy_SOLVEstruct_handle(SOLVEstruct)
      call f_destroy_SuperMatrix_handle(A)
      call f_destroy_SuperLUStat_handle(stat)

! Check malloc
!      call f_check_malloc(iam)

! Terminate the MPI execution environment
      call mpi_finalize(ierr)

      stop
      end


text:
c !DEC$ IF DEFINED (NEED_BLAS)
      SUBROUTINE sswap (N, SX, INCX, SY, INCY)
      USE LAPREC, ONLY: WP => SP
      IMPLICIT NONE
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER, INTENT(IN) :: N, INCX, INCY
      REAL(WP), INTENT(INOUT) :: SX(*), SY(*)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER :: IX, IY
      REAL(WP), DIMENSION(:), ALLOCATABLE :: STEMP
C-----------------------------------------------
c
c     interchanges two vectors.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
c
      IF (n .le. 0) RETURN

      ALLOCATE (stemp(n), stat=ix)
      IF (ix .ne. 0) STOP 'Allocation error in swap'

      IF (incx.ne.1 .or. incy.ne.1) THEN
c
c       code for unequal increments or equal increments not equal
c         to 1
c
         ix = 1
         iy = 1
         IF (incx .lt. 0) ix = ((-n) + 1)*incx + 1
         IF (incy .lt. 0) iy = ((-n) + 1)*incy + 1
         stemp(:n) = sx(ix:(n-1)*incx+ix:incx)
         sx(ix:(n-1)*incx+ix:incx) = sy(iy:(n-1)*incy+iy:incy)
         sy(iy:(n-1)*incy+iy:incy) = stemp(:n)

      ELSE
c
c     code for both increments equal to 1
c
         stemp(:n) = sx(:n)
         sx(:n) = sy(:n)
         sy(:n) = stemp(:n)
      END IF

      DEALLOCATE (stemp)

      END SUBROUTINE sswap

      SUBROUTINE dswap (N, SX, INCX, SY, INCY)
      USE LAPREC, ONLY: WP => DP
      IMPLICIT NONE
C-----------------------------------------------
C   D U M M Y   A R G U M E N T S
C-----------------------------------------------
      INTEGER, INTENT(IN) :: N, INCX, INCY
      REAL(WP), INTENT(INOUT) :: SX(*), SY(*)
C-----------------------------------------------
C   L O C A L   V A R I A B L E S
C-----------------------------------------------
      INTEGER :: IX, IY
      REAL(WP), DIMENSION(:), ALLOCATABLE :: STEMP
C-----------------------------------------------
c
c     interchanges two vectors.
c     jack dongarra, linpack, 3/11/78.
c     modified 12/3/93, array(1) declarations changed to array(*)
c
c
      IF (n .le. 0) RETURN

      ALLOCATE (stemp(n), stat=ix)
      IF (ix .ne. 0) STOP 'Allocation error in swap'

      IF (incx.ne.1 .or. incy.ne.1) THEN
c
c       code for unequal increments or equal increments not equal
c         to 1
c
         ix = 1
         iy = 1
         IF (incx .lt. 0) ix = ((-n) + 1)*incx + 1
         IF (incy .lt. 0) iy = ((-n) + 1)*incy + 1
         stemp(:n) = sx(ix:(n-1)*incx+ix:incx)
         sx(ix:(n-1)*incx+ix:incx) = sy(iy:(n-1)*incy+iy:incy)
         sy(iy:(n-1)*incy+iy:incy) = stemp(:n)

      ELSE
c
c     code for both increments equal to 1
c
         stemp(:n) = sx(:n)
         sx(:n) = sy(:n)
         sy(:n) = stemp(:n)
      END IF

      DEALLOCATE (stemp)

      END SUBROUTINE dswap
c !DEC$ ELSE
c       SUBROUTINE bla_stub
c       END SUBROUTINE bla_stub
c !DEC$ ENDIF


text:
program test_convert_ubt_and_wbv
  use mod_orrb
  implicit none
  !
  ! Note: This does not cover large rank cases (e.g. ubw=n-1).
  !
  character(len=40) :: test_name
  character(len=*), parameter :: fmt="(A40, 'Time: ',ES8.2,', lbw: ',I3,', ubw: ', I3,', error: ',ES8.2, ', ', A10)"

  real(kind=dp) :: t0, t1
  integer(kind=int32) :: na, ubwa, lbwa
  type(error_info) :: error
  real(kind=dp), parameter :: tol=3e-15, c=2.0
  !
  real(kind=dp), dimension(:,:), allocatable :: a0_d, a1_d
  complex(kind=dp), dimension(:,:), allocatable :: a0_z, a1_z

  type(d_wbv), allocatable :: wbv_d
  type(z_wbv), allocatable :: wbv_z
  type(d_ubt), allocatable :: ubt_d
  type(z_ubt), allocatable :: ubt_z

  call initialize_errors

  print *
  print *, "--------------------------------"
  print *
  print *, "Real WBV and UBT Conversion Tests:"
  print *
  na=40
  lbwa=5; ubwa=3
  wbv_d=d_random_wbv(na,lbwa,ubwa,error=error)
  a0_d=general(wbv_d,error)
  call cpu_time(t0)
  ubt_d=ubt_of(wbv_d,error)
  call cpu_time(t1)
  a1_d=general(ubt_d,error)
  test_name = "Real WBV to UBT, n=40;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,ubt_d%lbw, &
       ubwa, ubt_d%ubw, t0,t1,c*tol,error)

  na=1
  lbwa=0; ubwa=0
  wbv_d=d_random_wbv(na,lbwa,ubwa,error=error)
  a0_d=general(wbv_d,error)
  call cpu_time(t0)
  ubt_d=ubt_of(wbv_d,error)
  call cpu_time(t1)
  a1_d=general(ubt_d,error)
  test_name = "Real WBV to UBT, n=1;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,ubt_d%lbw, &
       ubwa, ubt_d%ubw, t0,t1,c*tol,error)

  na=2
  lbwa=1; ubwa=1
  wbv_d=d_random_wbv(na,lbwa,ubwa,error=error)
  a0_d=general(wbv_d,error)
  call cpu_time(t0)
  ubt_d=ubt_of(wbv_d,error)
  call cpu_time(t1)
  a1_d=general(ubt_d,error)
  test_name = "Real WBV to UBT, n=2;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,ubt_d%lbw, &
       ubwa, ubt_d%ubw, t0,t1,c*tol,error)

  na=3
  lbwa=1; ubwa=1
  wbv_d=d_random_wbv(na,lbwa,ubwa,error=error)
  a0_d=general(wbv_d,error)
  call cpu_time(t0)
  ubt_d=ubt_of(wbv_d,error)
  call cpu_time(t1)
  a1_d=general(ubt_d,error)
  test_name = "Real WBV to UBT, n=3;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,ubt_d%lbw, &
       ubwa, ubt_d%ubw, t0,t1,c*tol,error)

  na=4
  lbwa=2; ubwa=2
  wbv_d=d_random_wbv(na,lbwa,ubwa,error=error)
  a0_d=general(wbv_d,error)
  call cpu_time(t0)
  ubt_d=ubt_of(wbv_d,error)
  call cpu_time(t1)
  a1_d=general(ubt_d,error)
  test_name = "Real WBV to UBT, n=4;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,ubt_d%lbw, &
       ubwa, ubt_d%ubw, t0,t1,c*tol,error)
  print *

  na=40
  lbwa=5; ubwa=3
  ubt_d=d_random_ubt(na,lbwa,ubwa,error=error)
  a0_d=general(ubt_d,error)
  call cpu_time(t0)
  wbv_d=wbv_of(ubt_d,error)
  call cpu_time(t1)
  a1_d=general(wbv_d,error)
  test_name = "Real UBT to WBV, n=40;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,wbv_d%lbw, &
       ubwa, wbv_d%ubw, t0,t1,c*tol,error)

  na=1
  lbwa=0; ubwa=0
  ubt_d=d_random_ubt(na,lbwa,ubwa,error=error)
  a0_d=general(ubt_d,error)
  call cpu_time(t0)
  wbv_d=wbv_of(ubt_d,error)
  call cpu_time(t1)
  a1_d=general(wbv_d,error)
  test_name = "Real UBT to WBV, n=1;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,wbv_d%lbw, &
       ubwa, wbv_d%ubw, t0,t1,c*tol,error)

  na=2
  lbwa=1; ubwa=1
  ubt_d=d_random_ubt(na,lbwa,ubwa,error=error)
  a0_d=general(ubt_d,error)
  call cpu_time(t0)
  wbv_d=wbv_of(ubt_d,error)
  call cpu_time(t1)
  a1_d=general(wbv_d,error)
  test_name = "Real UBT to WBV, n=2;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,wbv_d%lbw, &
       ubwa, wbv_d%ubw, t0,t1,c*tol,error)

  na=3
  lbwa=1; ubwa=1
  ubt_d=d_random_ubt(na,lbwa,ubwa,error=error)
  a0_d=general(ubt_d,error)
  call cpu_time(t0)
  wbv_d=wbv_of(ubt_d,error)
  call cpu_time(t1)
  a1_d=general(wbv_d,error)
  test_name = "Real UBT to WBV, n=3;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,wbv_d%lbw, &
       ubwa, wbv_d%ubw, t0,t1,c*tol,error)

  na=4
  lbwa=2; ubwa=2
  ubt_d=d_random_ubt(na,lbwa,ubwa,error=error)
  a0_d=general(ubt_d,error)
  call cpu_time(t0)
  wbv_d=wbv_of(ubt_d,error)
  call cpu_time(t1)
  a1_d=general(wbv_d,error)
  test_name = "Real UBT to WBV, n=4;"
  call d_output_result(test_name,a0_d,a1_d,lbwa,wbv_d%lbw, &
       ubwa, wbv_d%ubw, t0,t1,c*tol,error)

  print *
  print *, "--------------------------------"
  print *
  print *, "Complex WBV and UBT Conversion Tests:"
  print *

  na=40
  lbwa=5; ubwa=3
  wbv_z=z_random_wbv(na,lbwa,ubwa,error=error)
  a0_z=general(wbv_z,error)
  call cpu_time(t0)
  ubt_z=ubt_of(wbv_z,error)
  call cpu_time(t1)
  a1_z=general(ubt_z,error)
  test_name = "Complex WBV to UBT, n=40;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,ubt_z%lbw, &
       ubwa, ubt_z%ubw, t0,t1,c*tol,error)

  na=1
  lbwa=0; ubwa=0
  wbv_z=z_random_wbv(na,lbwa,ubwa,error=error)
  a0_z=general(wbv_z,error)
  call cpu_time(t0)
  ubt_z=ubt_of(wbv_z,error)
  call cpu_time(t1)
  a1_z=general(ubt_z,error)
  test_name = "Complex WBV to UBT, n=1;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,ubt_z%lbw, &
       ubwa, ubt_z%ubw, t0,t1,c*tol,error)

  na=2
  lbwa=1; ubwa=1
  wbv_z=z_random_wbv(na,lbwa,ubwa,error=error)
  a0_z=general(wbv_z,error)
  call cpu_time(t0)
  ubt_z=ubt_of(wbv_z,error)
  call cpu_time(t1)
  a1_z=general(ubt_z,error)
  test_name = "Complex WBV to UBT, n=2;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,ubt_z%lbw, &
       ubwa, ubt_z%ubw, t0,t1,c*tol,error)

  na=3
  lbwa=1; ubwa=1
  wbv_z=z_random_wbv(na,lbwa,ubwa,error=error)
  a0_z=general(wbv_z,error)
  call cpu_time(t0)
  ubt_z=ubt_of(wbv_z,error)
  call cpu_time(t1)
  a1_z=general(ubt_z,error)
  test_name = "Complex WBV to UBT, n=3;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,ubt_z%lbw, &
       ubwa, ubt_z%ubw, t0,t1,c*tol,error)

  na=4
  lbwa=2; ubwa=2
  wbv_z=z_random_wbv(na,lbwa,ubwa,error=error)
  a0_z=general(wbv_z,error)
  call cpu_time(t0)
  ubt_z=ubt_of(wbv_z,error)
  call cpu_time(t1)
  a1_z=general(ubt_z,error)
  test_name = "Complex WBV to UBT, n=4;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,ubt_z%lbw, &
       ubwa, ubt_z%ubw, t0,t1,c*tol,error)
  print *

  na=40
  lbwa=5; ubwa=3
  ubt_z=z_random_ubt(na,lbwa,ubwa,error=error)
  a0_z=general(ubt_z,error)
  call cpu_time(t0)
  wbv_z=wbv_of(ubt_z,error)
  call cpu_time(t1)
  a1_z=general(wbv_z,error)
  test_name = "Complex UBT to WBV, n=40;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,wbv_z%lbw, &
       ubwa,wbv_z%ubw,t0,t1,c*tol,error)

  na=1
  lbwa=0; ubwa=0
  ubt_z=z_random_ubt(na,lbwa,ubwa,error=error)
  a0_z=general(ubt_z,error)
  call cpu_time(t0)
  wbv_z=wbv_of(ubt_z,error)
  call cpu_time(t1)
  a1_z=general(wbv_z,error)
  test_name = "Complex UBT to WBV, n=1;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,wbv_z%lbw, &
       ubwa,wbv_z%ubw,t0,t1,c*tol,error)

  na=2
  lbwa=1; ubwa=1
  ubt_z=z_random_ubt(na,lbwa,ubwa,error=error)
  a0_z=general(ubt_z,error)
  call cpu_time(t0)
  wbv_z=wbv_of(ubt_z,error)
  call cpu_time(t1)
  a1_z=general(wbv_z,error)
  test_name = "Complex UBT to WBV, n=2;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,wbv_z%lbw, &
       ubwa,wbv_z%ubw,t0,t1,c*tol,error)

  na=3
  lbwa=1; ubwa=1
  ubt_z=z_random_ubt(na,lbwa,ubwa,error=error)
  a0_z=general(ubt_z,error)
  call cpu_time(t0)
  wbv_z=wbv_of(ubt_z,error)
  call cpu_time(t1)
  a1_z=general(wbv_z,error)
  test_name = "Complex UBT to WBV, n=3;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,wbv_z%lbw, &
       ubwa,wbv_z%ubw,t0,t1,c*tol,error)

  na=4
  lbwa=2; ubwa=2
  ubt_z=z_random_ubt(na,lbwa,ubwa,error=error)
  a0_z=general(ubt_z,error)
  call cpu_time(t0)
  wbv_z=wbv_of(ubt_z,error)
  call cpu_time(t1)
  a1_z=general(wbv_z,error)
  test_name = "Complex UBT to WBV, n=4;"
  call z_output_result(test_name,a0_z,a1_z,lbwa,wbv_z%lbw, &
       ubwa,wbv_z%ubw,t0,t1,c*tol,error)

contains

  subroutine d_output_result(name,a0,a1,lbw0,lbw1,ubw0,ubw1,t0,t1,bnd,error)
    character(len=*) :: name
    real(kind=dp), dimension(:,:) :: a0, a1
    real(kind=dp) :: bnd, t0, t1
    integer(kind=int32) :: lbw0, lbw1, ubw0, ubw1
    type(error_info) :: error

    real(kind=dp) :: berr
    character(len=10) :: test_result

    if (error%code > 0) then
       print *, "Calling error in test: ", name
    else
       berr = maxabs(a1-a0)
       if (lbw0==lbw1 .and. ubw0==ubw1 .and. berr < bnd) then
          test_result="PASSED"
       else
          test_result="    FAILED"
       end if
       write (*,fmt) name, t1-t0, lbw1, ubw1, berr, test_result
    end if
  end subroutine d_output_result

  subroutine z_output_result(name,a0,a1,lbw0,lbw1,ubw0,ubw1,t0,t1,bnd,error)
    character(len=*) :: name
    complex(kind=dp), dimension(:,:) :: a0, a1
    real(kind=dp) :: bnd, t0, t1
    integer(kind=int32) :: lbw0, lbw1, ubw0, ubw1
    type(error_info) :: error

    real(kind=dp) :: berr
    character(len=10) :: test_result

    if (error%code > 0) then
       print *, "Calling error in test: ", name
    else
       berr = maxabs(a1-a0)
       if (lbw0==lbw1 .and. ubw0==ubw1 .and. berr < bnd) then
          test_result="PASSED"
       else
          test_result="    FAILED"
       end if
       write (*,fmt) name, t1-t0, lbw1, ubw1, berr, test_result
    end if
  end subroutine z_output_result

end program test_convert_ubt_and_wbv

max_stars_repo_name: m-a-stewart/BandGivensWeight

max_issues_repo_name: m-a-stewart/BandGivensWeight

max_forks_repo_name: m-a-stewart/BandGivensWeight


! { dg-do run }
! { dg-options "-fcheck=all" }
!
! PR 50225: [OOP] The allocation status for polymorphic allocatable function results is not set properly
!
! Contributed by Arjen Markus <arjen.markus895@gmail.com>

module points2d

  implicit none

  type point2d
      real :: x, y
  end type

contains

 subroutine print( point )
   class(point2d) :: point
   write(*,'(2f10.4)') point%x, point%y
 end subroutine

 subroutine random_vector( point )
   class(point2d) :: point
   call random_number( point%x )
   call random_number( point%y )
   point%x = 2.0 * (point%x - 0.5)
   point%y = 2.0 * (point%y - 0.5)
 end subroutine

 function add_vector( point, vector )
   class(point2d), intent(in)  :: point, vector
   class(point2d), allocatable :: add_vector
   allocate( add_vector )
   add_vector%x = point%x + vector%x
   add_vector%y = point%y + vector%y
 end function

end module points2d

program random_walk

  use points2d
  implicit none

  type(point2d), target   :: point_2d, vector_2d
  class(point2d), pointer :: point, vector
  integer :: i

  point  => point_2d
  vector => vector_2d

  do i=1,2
    call random_vector(point)
    call random_vector(vector)
    call print(add_vector(point, vector))
  end do

end program random_walk


text:
       function symstore_2s(T)
        implicit none

        integer :: w
        type(Tensor2s), intent(in) :: T
        type(Tensor2s) :: symstore_2s

        symstore_2s%a6 = T%a6
       end function symstore_2s

       function symstore_2sa(T)
        implicit none

        integer :: w
        real(kind=8), dimension(6), intent(in) :: T
        type(Tensor2s) :: symstore_2sa

        symstore_2sa%a6 = T

       end function symstore_2sa

       function symstore_4s(T)
        implicit none

        integer :: w1,w2
        type(Tensor4s), intent(in) :: T
        type(Tensor4s) :: symstore_4s

        symstore_4s%a6b6 = T%a6b6

       end function symstore_4s

       function symstore_4sa(T)
        implicit none

        integer :: w1,w2
        real(kind=8), dimension(6,6), intent(in) :: T
        type(Tensor4s) :: symstore_4sa

        symstore_4sa%a6b6 = T

       end function symstore_4sa

       function symstore_2(T)
        implicit none

        integer :: w
        type(Tensor2), intent(in) :: T
        type(Tensor2s) :: symstore_2
        integer :: i

        w = 1

        symstore_2%a6 = 0.d0
        do i=1,3
         symstore_2%a6(i) = T%ab(i,i)
        enddo
        symstore_2%a6(4) = T%ab(1,2)*dble(w)
        symstore_2%a6(5) = T%ab(2,3)*dble(w)
        symstore_2%a6(6) = T%ab(3,1)*dble(w)

       end function symstore_2

       function symstore_4(T)
        implicit none

        integer :: v1,v2
        type(Tensor4), intent(in) :: T
        type(Tensor4s) :: symstore_4
        real(kind=8) :: w
        integer :: i,j

        v1 = 1
        v2 = 1
        w = dble(v1)

        symstore_4%a6b6 = 0.d0
        do i=1,3
         do j=1,3
          symstore_4%a6b6(i,j) = T%abcd(i,i,j,j)
         enddo
        enddo

        symstore_4%a6b6(4,4) = T%abcd(1,2,1,2)!+T%abcd(2,1,1,2))/2.d0*w
        symstore_4%a6b6(5,5) = T%abcd(2,3,2,3)!+T%abcd(3,2,2,3))/2.d0*w
        symstore_4%a6b6(6,6) = T%abcd(3,1,3,1)!+T%abcd(1,3,3,1))/2.d0*w

        do i=1,3
         symstore_4%a6b6(i,4) = T%abcd(i,i,1,2)
         symstore_4%a6b6(i,5) = T%abcd(i,i,2,3)
         symstore_4%a6b6(i,6) = T%abcd(i,i,3,1)
         symstore_4%a6b6(4,i) = T%abcd(1,2,i,i)
         symstore_4%a6b6(5,i) = T%abcd(2,3,i,i)
         symstore_4%a6b6(6,i) = T%abcd(3,1,i,i)
        enddo

        symstore_4%a6b6(4,5) = T%abcd(1,2,2,3)!+T%abcd(2,1,2,3))/2.d0*w
        symstore_4%a6b6(4,6) = T%abcd(1,2,3,1)!+T%abcd(2,1,3,1))/2.d0*w

        symstore_4%a6b6(5,4) = T%abcd(2,3,1,2)!+T%abcd(3,2,1,2))/2.d0*w
        symstore_4%a6b6(5,6) = T%abcd(2,3,3,1)!+T%abcd(3,2,3,1))/2.d0*w

        symstore_4%a6b6(6,4) = T%abcd(3,1,1,2)!+T%abcd(1,3,1,2))/2.d0*w
        symstore_4%a6b6(6,5) = T%abcd(3,1,2,3)!+T%abcd(1,3,2,3))/2.d0*w

       end function symstore_4

max_stars_repo_name: kengwit/ttb

max_issues_repo_name: kengwit/ttb

max_forks_repo_name: kengwit/ttb


text:
!	Created by Kerry Meyer on 2/8/08.
!
!	Parameters, variables, and equations for Rayleigh optical depth calculations.
!	Constants and formulas from:
!		Bodhaine et al., 1999: On Rayleigh Optical Depth Calculations. J. Atmos. and Ocean Tech., 16. 1854-1861.
!
!	Parameters:
!		CO2_ppm		= concentration of CO2 in ppm
!		F_Ar		= depolarization factor for argon
!		F_CO2		= depolarization factor for CO2
!		ndens_da	= number density of dry air
!		av_num		= Avogadro's number
!		g0			= gravitational constant
!
!	Variables:
!		atm_press	= atmospheric pressure at each level from atmospheric profile
!		F_O2		= depolarization factor for O2 at each channel
!		F_N2		= depolarization factor for N2 at each channel
!		F_air		= depolarization factor for dry air at each channel
!		refin_air	= refractive index of dry air at each channel
!		ray_cross	= Rayleigh scattering cross-section at each channel
!		mw_air		= molecular weight of dry air
!		g1,g2		= gravity, corrected for altitude

module dis_rayleigh

use dis_global_pars

implicit none
public

real, parameter :: CO2_ppm = 360, &
				   F_Ar = 1.0, &
				   F_CO2 = 1.15, &
				   ndens_da = 2.546899E19, &
				   av_num = 6.0221367E23, &
				   g0 = 9.832

real, dimension(maxulv) :: atm_press, atm_height

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

contains

subroutine readprof (atmprof,nlyr,atm_press,atm_height)

use dis_global_pars

implicit none

!	Dummy variables
character(len=100), intent(in) :: atmprof

real, dimension(maxulv), intent(out) :: atm_press, atm_height

integer, intent(in) :: nlyr

!	Local variables
integer :: i
logical :: debug
debug = .false.
	open (50,file=atmprof)

	!print*,'nlyr: ',nlyr
	do i = 1, nlyr+1

		read (50,*) atm_height(i), atm_press(i)
        if (debug) print*,i,atm_height(i), atm_press(i)
	end do

	close (50)

end subroutine readprof

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

subroutine rayleigh_tau_calc (nlyr,ich,atm_press,atm_height,wvlo,wvhi,rayleigh_tau)

use dis_global_pars

implicit none

!	Dummy variables
real, dimension(:), intent(inout) :: rayleigh_tau
real, dimension(:), intent(in) :: atm_press, atm_height
real, intent(in) :: wvlo, wvhi

integer, intent(in) :: nlyr

!	Local variables
real :: wvln
real :: ray_cross_num, ray_cross_den, ray_cross
real :: F_O2, F_N2, F_air, refin_air
real :: mw_air, g1, g2

integer :: i, ich
logical :: debug
debug = .false.
!	Calculate molecular weight of air (corresponding to CO2 concentration)
	mw_air = 15.0556*(CO2_ppm/1000000.) + 28.9595

!	Calculate Rayleigh scattering cross-section
	wvln = (1. / ((wvlo + wvhi) / 2.)) * 10000.
	F_N2 = 1.034 + 3.17E-4*(1./(wvln**2))
	F_O2 = 1.096 + 1.385E-3*(1./(wvln**2)) + 1.448E-4*(1./(wvln**4))
	F_air = (78.084*F_N2 + 20.946*F_O2 + 0.934*F_Ar + CO2_ppm*F_CO2/10000.) / &
				(78.084 + 20.946 + 0.934 + CO2_ppm/10000.)
	refin_air = ((8060.51 + (2480990./(132.274-wvln**-2)) + (17455.7/(39.32957-wvln**-2)))*&
					(1. + 0.54*((CO2_ppm/1000000.)-0.0003)) / 10**8) + 1.
	ray_cross_num = 24. * (acos(-1.)**3) * ((refin_air**2 - 1.)**2)
	ray_cross_den = ((wvln*0.0001)**4) * ((refin_air**2 + 2.)**2) * ndens_da
	ray_cross_den = ray_cross_den * ndens_da
	ray_cross = (ray_cross_num / ray_cross_den) * F_air

	do i = 1, nlyr
		g1 = g0 * (1.-(2.*atm_height(i+1)/6.37E6))
		g2 = g0 * (1.-(2.*atm_height(i)/6.37E6))
		rayleigh_tau(i) = ((ray_cross*av_num) / mw_air) * (atm_press(i+1)/g1 - atm_press(i)/g2) * 10.
		if (debug) print*,i,rayleigh_tau(i)
	end do

end subroutine rayleigh_tau_calc

end module dis_rayleigh

max_stars_repo_name: GEOS-ESM/AeroApps

max_issues_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


text:
      MODULE el53_I
      INTERFACE
!
      SUBROUTINE EL53(JJA,JJB,JA,JB,JC,JD,IREZ,ICOLBREI)
      INTEGER, INTENT(IN) :: JJA,JJB,JA,JB,JC,JD,IREZ,ICOLBREI
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


text:
* $Id$
c-----------------------------------------------------------------------
C This files contains a set of routines used to convert pair
c quantities into quartet quantities. These routines used to
c be included in different files, usually where they are called
c from. They are called with ...BL(IADDRESS).. double precision
c parameters which then are of the INTEGER type in the routines.
c It causes problems of incompatibile types for some compilers.
c-----------------------------------------------------------------------
c List of subroutines with previous location :
c
c amshift.f:      subroutine convr3(bl,m,nbls,npij,npkl,idx1,idx2,
c
c derivat.f:      subroutine conv24x(nbls,npij,npkl,idx1,idx2 ,
c derivat.f:      subroutine conv24r(nbls,npij,idx1,xab,xabq)
c
c spec_calcint.f: subroutine conv1x_1(nbls1,mmax1,npij,lcij, idx1,indx,
c spec_calcint.f: subroutine conv1der_1(nbls1,npij,lci,idx1,indx, aa,
c spec_calcint.f: subroutine conv1der_2(nbls1,npij,lci,idx1,indx, aa,
c spec_calcint.f: subroutine conv1x_2(nbls1,mmax1,npij,lcij, idx1,indx,
c spec_calcint.f: subroutine conv2x(nbls1,nfumax1,npkl,lckl, idx2,indx,
c-----------------------------------------------------------------------
c=======================================================================
      subroutine convr3(bl,m,nbls,npij,npkl,idx1,idx2,
     *                   xab,xcd, ixabn,ixcdn)
      implicit real*8 (a-h,o-z)
      dimension bl(*)
      dimension idx1(*),idx2(*)
      dimension xab(npij,3),xcd(npkl,3)
c
      nbls1=nbls
      nbls2=nbls*2
      nbls3=nbls*3
      nbls1=nbls1*m
      nbls2=nbls2*m
      nbls3=nbls3*m
      call getmem(nbls3,ixabn)
      call getmem(nbls3,ixcdn)
c
       ixab1=ixabn-1
       ixcd1=ixcdn-1
c
      ijklnmr=0
      do 100 ijkl=1,nbls
      ijpar=idx1(ijkl)
      klpar=idx2(ijkl)
c
      xab1=xab(ijpar,1)
      xab2=xab(ijpar,2)
      xab3=xab(ijpar,3)
      xcd1=xcd(klpar,1)
      xcd2=xcd(klpar,2)
      xcd3=xcd(klpar,3)
c
        do 100 nmr=1,m
        ijklnmr=ijklnmr+1
        bl(ixab1+ijklnmr)      =xab1
        bl(ixab1+ijklnmr+nbls1)=xab2
        bl(ixab1+ijklnmr+nbls2)=xab3
c
        bl(ixcd1+ijklnmr)      =xcd1
        bl(ixcd1+ijklnmr+nbls1)=xcd2
        bl(ixcd1+ijklnmr+nbls2)=xcd3
c
  100 continue
      return
      end
c=======================================================================
      subroutine conv24x(nbls,npij,npkl,idx1,idx2 ,
     *                  xab ,xcd, xyab, xycd ,
     *                  xabq,xcdq,xyabq,xycdq )
      implicit real*8 (a-h,o-z)
c
      dimension idx1(nbls),idx2(nbls)
      dimension xab(npij,3) ,xcd(npkl,3) ,xyab(npij,3) ,xycd(npkl,3)
      dimension xabq(nbls,3),xcdq(nbls,3),xyabq(nbls,3),xycdq(nbls,3)
c
      do 100 ijkl=1,nbls
      ijpar=idx1(ijkl)
      klpar=idx2(ijkl)
        do 150 i=1,3
        xabq(ijkl,i)=xab(ijpar,i)
        xcdq(ijkl,i)=xcd(klpar,i)
        xyabq(ijkl,i)=xyab(ijpar,i)
        xycdq(ijkl,i)=xycd(klpar,i)
  150   continue
  100 continue
c
      end
c=======================================================================
      subroutine conv24r(nbls,npij,idx1,xab,xabq)
      implicit real*8 (a-h,o-z)
c
      dimension idx1(nbls)
      dimension xab(npij,3),xabq(nbls,3)
c
      do 100 ijkl=1,nbls
      ijpar=idx1(ijkl)
c     klpar=idx2(ijkl)
        do 150 i=1,3
        xabq(ijkl,i)=xab(ijpar,i)
c       xcdq(ijkl,i)=xcd(klpar,i)
  150   continue
  100 continue
      end
c=======================================================================
      subroutine conv1x_1(nbls1,mmax1,npij,lcij, idx1,indx,
     *                    abnia,xpn,abnix,xpnx )
c-------------------------------------------------------------------
c npij = number of uniqe pairs now
c-------------------------------------------------------------------
c
      implicit real*8 (a-h,o-z)
      dimension idx1(*),indx(*)
      dimension xpn(npij,3,*)
      dimension abnia(npij,mmax1,*)
c
      dimension xpnx(nbls1,3)
      dimension abnix(nbls1,mmax1)
c
      do 10 i=1,nbls1
      ijkl=indx(i)
      ijpar=idx1(ijkl)
        xpnx(i,1)=xpn(ijpar,1,lcij)
        xpnx(i,2)=xpn(ijpar,2,lcij)
        xpnx(i,3)=xpn(ijpar,3,lcij)
   10 continue
c
      do 20 m=1,mmax1
      do 20 i=1,nbls1
      ijkl=indx(i)
      ijpar=idx1(ijkl)
        abnix(i,m)=abnia(ijpar,m,lcij)
   20 continue
      end
c=======================================================================
      subroutine conv1x_2(nbls1,mmax1,npij,lcij, idx1,indx,xpn,xpnx )
c
c npij = number of uniqe pairs now
c
c
      implicit real*8 (a-h,o-z)
      dimension idx1(*),indx(*)
      dimension xpn(npij,3,*)
      dimension xpnx(nbls1,3)
c
      do 10 i=1,nbls1
      ijkl=indx(i)
      ijpar=idx1(ijkl)
        xpnx(i,1)=xpn(ijpar,1,lcij)
        xpnx(i,2)=xpn(ijpar,2,lcij)
        xpnx(i,3)=xpn(ijpar,3,lcij)
   10 continue
c
      end
c=======================================================================
      subroutine conv1der_1(nbls1,npij,lci,idx1,indx, aa, aax)
c
c npij = number of uniqe pairs now
c exponents are already rescaled by 2 in precal2a_1
c
      implicit real*8 (a-h,o-z)
      dimension idx1(*),indx(*)
      dimension aa(npij,*)
c output :
      dimension aax(nbls1)
c
      do 10 i=1,nbls1
      ijkl=indx(i)
      ijpar=idx1(ijkl)
cccc    aax(i)=aa(ijpar,lci)*2.0d0   ! already rescaled
        aax(i)=aa(ijpar,lci)
   10 continue
c
      end
c=======================================================================
      subroutine conv1der_2(nbls1,npij,lci,idx1,indx, aa, aax)
c
c npij = number of uniqe pairs now
c
      implicit real*8 (a-h,o-z)
      dimension idx1(*),indx(*)
      dimension aa(npij,*)
c output :
      dimension aax(nbls1)
c-------------------------------------------------------------------
c this is for iroute=2 only: all exponents in a block are the same:
c exponents are already rescaled by 2 in precal2a_2
c-------------------------------------------------------------------
      aax(1)=aa(1,lci)
c
c     ijkl1 =indx(1)
c     ijpar1=idx1(ijkl1)
c     aax(1)=aa(ijpar1,lci)*2.d0
c----------------------------------------------
c     do 10 i=1,nbls1
c     ijkl=indx(i)
c     ijpar=idx1(ijkl)
c       aax(i)=aa(ijpar,lci)*2.0d0
c  10 continue
c----------------------------------------------
c
      end
c=======================================================================
      subroutine conv2x(nbls1,nfumax1,npkl,lckl, idx2,indx,
     *                  habcd,nfumax, habcdx  )
c
c npkl = number of uniqe pairs now
c
      implicit real*8 (a-h,o-z)
      dimension idx2(*),indx(*)
      dimension habcd(npkl,3,nfumax,*)
      dimension habcdx(nbls1,3,nfumax)
c
         do 32 ifu=1,nfumax1
         do 32 i=1,nbls1
         ijkl=indx(i)
         klpar=idx2(ijkl)
           habcdx(i,1,ifu)=habcd(klpar,1,ifu,lckl)
           habcdx(i,2,ifu)=habcd(klpar,2,ifu,lckl)
           habcdx(i,3,ifu)=habcd(klpar,3,ifu,lckl)
   32    continue
      end
c=======================================================================

max_stars_repo_name: dinisAbranches/nwchem

max_issues_repo_name: dinisAbranches/nwchem

max_forks_repo_name: dinisAbranches/nwchem


text:
! a driver that sets up a circular MAC velocity field and initializes
! particles at various radii.  We then advect the particles around
! for one period.  We can check the accurary of the particle advection
! by comparing the initial radius to the final radius.

subroutine varden()

  use BoxLib
  use f2kcli
  use ml_boxarray_module
  use layout_module
  use multifab_module
  use ml_cc_restriction_module
  use define_bc_module
  use bl_mem_stat_module
  use bl_timer_module
  use box_util_module
  use bl_IO_module
  use fabio_module
  use variables, only: nscal, init_variables, rho_comp, spec_comp
  use geometry, only:  spherical, &
                       init_spherical, init_center, &
                       init_cutoff, initialize_dx
  use probin_module, only: prob_lo, prob_hi, pmask, &
                           test_set, cflfac, &
                           stop_time, vel_amp, &
                           max_levs
  use runtime_init_module
  use bl_constants_module
  use particle_module
  use test_particles_module, only: init_umac_2d

  implicit none

  integer, parameter :: NCOMPS = 1   !Number of components in state multifab
  integer, parameter :: NUMGHOST = 1 !Number of ghost cells
  integer, parameter :: MV_COMP = 1  !Index for velocity magnitude

  integer :: i,n,comp
  integer :: ng_s, ng_um
  integer :: nlevs, dm

  type(ml_layout) :: mla

  type(multifab), allocatable :: umac(:,:)
  type(multifab), allocatable :: s(:)

  real(kind=dp_t), pointer :: ump(:,:,:,:), vmp(:,:,:,:), wmp(:,:,:,:)
  real(kind=dp_t), pointer :: sp(:,:,:,:)

  real(kind=dp_t), pointer :: dx(:,:)

  integer, allocatable :: lo(:),hi(:)

  type(ml_boxarray) :: mba

  type(bc_tower) ::  the_bc_tower

  real(kind=dp_t) :: t, dt, maxR

  type(particle_container) :: particles
  real(kind=dp_t), allocatable :: point(:)

  integer :: nparticles
  real(kind=dp_t) :: dr_part, dtheta_part, r, theta

  ! general Maestro initializations
  call runtime_init()
  call init_spherical()
  call init_center()

  call init_variables()

  ! setup the grid
  call read_a_hgproj_grid(mba, test_set)

  call ml_layout_build(mla,mba,pmask)

  ! check for proper nesting
  if (.not. ml_boxarray_properly_nested(mla%mba, 3, pmask)) then
     call bl_error('fixed_grids not properly nested')
  end if

  ! initialize nlevs
  nlevs = mla%nlevel

  if (nlevs .ne. max_levs) then
     call bl_error('varden.f90: nlevs .ne. max_levs not supported yet')
  end if

  ! initialize dm
  dm = mla%dim

  ! initialize boundary conditions
  call initialize_bc(the_bc_tower,nlevs,pmask)
  do n = 1,nlevs
     call bc_tower_level_build(the_bc_tower,n,mla%la(n))
  end do

  ! allocate states
  allocate(umac(nlevs,dm))
  allocate(s(nlevs))

  if (spherical == 1) then
     call bl_error("ERROR: test_particles not defined for spherical = 1")
  endif

  ! maximum length from center
  if (dm == 2) then
     maxR = max(HALF*(prob_hi(1)-prob_lo(1)), &
                HALF*(prob_hi(2)-prob_lo(2)))
  else
     maxR = max(HALF*(prob_hi(1)-prob_lo(1)), &
                HALF*(prob_hi(2)-prob_lo(2)), &
                HALF*(prob_hi(3)-prob_lo(3)))
  endif

  ! build MAC velocity
  do n = 1,nlevs
     call multifab_build(s(n), mla%la(n), NCOMPS, NUMGHOST)

     do comp=1,dm
        call multifab_build_edge(umac(n,comp), mla%la(n),NCOMPS,NUMGHOST,comp)
     end do
  end do

  ! initialize_dx
  call initialize_dx(dx,mba,nlevs)

  ! other allocations
  allocate(lo(dm))
  allocate(hi(dm))

  ! initialize the particles
  call build(particles)

  allocate(point(dm))

  nparticles = 20
  dr_part = maxR/nparticles
  dtheta_part = TWO*M_PI/nparticles

  do i = 1, nparticles
     theta = dble(i-ONE)*dtheta_part
     r = dble(i-HALF)*dr_part

     point(1) = r*cos(theta) + HALF*(prob_lo(1) + prob_hi(1))
     point(2) = r*sin(theta) + HALF*(prob_lo(2) + prob_hi(2))

     call add(particles,point,mla,dx,prob_lo)
  enddo

  ! initialize the MAC velocity field  ! negative velocity
  ng_um = nghost(umac(1,1))
  ng_s  = nghost(s(1))

  do n = 1, nlevs
     do i = 1, nfabs(umac(n,1))

        ump  => dataptr(umac(n,1),i)
        vmp  => dataptr(umac(n,2),i)

        sp => dataptr(s(n), i)

        lo = lwb(get_box(umac(n,1), i))
        hi = upb(get_box(umac(n,1), i))

        select case (dm)
        case (2)
           call init_umac_2d(ump(:,:,1,MV_COMP), vmp(:,:,1,MV_COMP), ng_um, &
                              sp(:,:,1,MV_COMP), ng_s, &
                             lo, hi, dx(n,:))

        case (3)
           wmp  => dataptr(umac(n,3),i)
           call bl_error("ERROR: 3-d not implemented")
        end select
     end do
  enddo

  ! write out a plotfile and timestamp that contains the
  ! magnitude of the velocity field
  call timestamp(particles, 'timestamp', s, &
                (/MV_COMP/), (/"magvel"/), ZERO)
  call fabio_ml_multifab_write_d(s,mla%mba%rr(:,1), &
                                 "magvel_field")

  ! compute the initial timestep -- dt = dx / |U|.  The velocity field
  ! has the magnitude |U| = A r, where A is the amplitude and r is the
  ! radius.
  dt = cflfac*dx(nlevs,1)/(vel_amp*maxR)

  ! advance the particls using the MAC velocity field
  t = ZERO
  do while (t < stop_time)

     print *, 't = ', t, 'dt = ', dt

     ! advance density according to rho_t + (rho U)_x = 0
     call particle_container_move_advect(particles,mla,umac,dx,dt, &
                                         prob_lo,prob_hi)

     ! update the time
     t = t + dt

     ! write out a particle timestamp
     call timestamp(particles, 'timestamp', s, &
                    (/MV_COMP/), (/"magvel"/), t)

     ! adjust the timestep, if necessary
     if (t + dt > stop_time) then
        dt = stop_time - t
     endif

  end do

  print *, 'finished evolution, t = ', t

  ! clean-up
  do n = 1,nlevs
     call destroy(s(n))
     do comp=1,dm
        call destroy(umac(n,comp))
     end do
  end do

  call destroy(mla)
  call destroy(mba)

  deallocate(s)
  deallocate(umac)

  call bc_tower_destroy(the_bc_tower)

  call destroy(particles)

  call runtime_close()

end subroutine varden

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO


text:
C NCLFORTSTART
      subroutine evssm (nvar, lssm, cssm, neval, weval, wevec,info,iopt)
      implicit none
                                                  ! input
      integer          nvar, lssm, info, neval, iopt
      double precision cssm(lssm)
                                                  ! output
      double precision weval(neval), wevec(nvar,neval)
C NCLEND
      integer          il  , iu  , mout, ldz, n, nn
     +                ,iwork(5*nvar), ifail(nvar)
      double precision vl  , vu  , abstol, work(8*nvar)
      character*1      jobz, rang, uplo

c .   jobz    = "V"       ; compute eigen{values/vectors}
c .   rang    = "I"       ; compute user specified iL to iU eigenvalues
c .   uplo    = "U"       ; input is lower(L)/upper(U) triangle
c .   vl      = 1d20      ; not used if:  range= "A" or "I"
c .   vu      = 1d20      ; not used
c .   neval   = iu-il+1   ; number of eigenvalues and vectors
c .   abstol  = 0.0       ; this mean use the default tolerance
c .   mout    = 0         ; ouput number of eigenvalues

      jobz    = "V"
      rang    = "I"
      uplo    = "U"

      il      = max(nvar-neval+1,1)
      iu      = nvar

      mout    = -99
      ldz     = nvar

      vl      = 1d20
      vu      = 1d20
      abstol  = 0.0d0

      call dspevx(jobz, rang, uplo, nvar, cssm, vl, vu, il, iu
     +           ,abstol ,mout, weval, wevec, ldz, work, iwork
     +           ,ifail, info)

      if (iopt.eq.0) then

c reverse the order of arrays returned from "dspevx" so
c .   largest eigenvalues/vectors are first.

          do n=1,neval
             work(n) = weval(n)
          end do
          do n=1,neval
             weval(n) = work(neval+1-n)
          end do

          do n=1,nvar
              do nn=1,neval
                 work(nn) = wevec(n,nn)
              end do
              do nn=1,neval
                 wevec(n,nn) = work(neval+1-nn)
              end do
          end do
      end if

      return
      end

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
This has multiple meanings.  You are probably looking for information on one of the following:

  Department of Nutrition  The UC Davis Department of Nutrition
  Nutrition Major  An academic major offered by the aforementioned department.

This is a Disambiguation disambiguation page  a navigational aid which lists other pages that might otherwise share the same title. If an article link referred you here, you might want to go back and fix it to point directly to the intended page.


!----------------------------------------------------------------------
!----------------------------------------------------------------------
!   bvp :    A nonlinear ODE eigenvalue problem
!----------------------------------------------------------------------
!----------------------------------------------------------------------

SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
!--------- ----

  IMPLICIT NONE
  INTEGER, INTENT(IN) :: NDIM, IJAC, ICP(*)
  DOUBLE PRECISION, INTENT(IN) :: U(NDIM), PAR(*)
  DOUBLE PRECISION, INTENT(OUT) :: F(NDIM)
  DOUBLE PRECISION, INTENT(INOUT) :: DFDU(NDIM,*), DFDP(NDIM,*)

  DOUBLE PRECISION pi

  pi = 4.d0*ATAN(1.d0)

  F(1) = U(2)
  F(2) = -( PAR(1)*pi )**2 * U(1) + U(1)**2

END SUBROUTINE FUNC

SUBROUTINE STPNT(NDIM,U,PAR,T)
!--------- -----

  IMPLICIT NONE
  INTEGER, INTENT(IN) :: NDIM
  DOUBLE PRECISION, INTENT(INOUT) :: U(NDIM), PAR(*)
  DOUBLE PRECISION, INTENT(IN) :: T

  PAR(1:2) = 0.0d0

  U = 0.0d0

END SUBROUTINE STPNT

SUBROUTINE BCND(NDIM,PAR,ICP,NBC,U0,U1,FB,IJAC,DBC)
!--------- ----

  IMPLICIT NONE
  INTEGER, INTENT(IN) :: NDIM, ICP(*), NBC, IJAC
  DOUBLE PRECISION, INTENT(IN) :: PAR(*), U0(NDIM), U1(NDIM)
  DOUBLE PRECISION, INTENT(OUT) :: FB(NBC)
  DOUBLE PRECISION, INTENT(INOUT) :: DBC(NBC,*)

  FB = (/ U0(1), U1(1) /)

END SUBROUTINE BCND

SUBROUTINE ICND(NDIM,PAR,ICP,NINT,U,UOLD,UDOT,UPOLD,FI,IJAC,DINT)
!--------- ----

  IMPLICIT NONE
  INTEGER, INTENT(IN) :: NDIM, ICP(*), NINT, IJAC
  DOUBLE PRECISION, INTENT(IN) :: PAR(*)
  DOUBLE PRECISION, INTENT(IN) :: U(NDIM), UOLD(NDIM), UDOT(NDIM), UPOLD(NDIM)
  DOUBLE PRECISION, INTENT(OUT) :: FI(NINT)
  DOUBLE PRECISION, INTENT(INOUT) :: DINT(NINT,*)

  FI = (/ U(1)-PAR(2) /)

END SUBROUTINE ICND

SUBROUTINE FOPT
END SUBROUTINE FOPT

SUBROUTINE PVLS
END SUBROUTINE PVLS