text:
*> \brief \b DGET52
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       SUBROUTINE DGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
*                          ALPHAI, BETA, WORK, RESULT )
*
*       .. Scalar Arguments ..
*       LOGICAL            LEFT
*       INTEGER            LDA, LDB, LDE, N
*       ..
*       .. Array Arguments ..
*       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
*      $                   B( LDB, * ), BETA( * ), E( LDE, * ),
*      $                   RESULT( 2 ), WORK( * )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> DGET52  does an eigenvector check for the generalized eigenvalue
*> problem.
*>
*> The basic test for right eigenvectors is:
*>
*>                           | b(j) A E(j) -  a(j) B E(j) |
*>         RESULT(1) = max   -------------------------------
*>                      j    n ulp max( |b(j) A|, |a(j) B| )
*>
*> using the 1-norm.  Here, a(j)/b(j) = w is the j-th generalized
*> eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th
*> generalized eigenvalue of m A - B.
*>
*> For real eigenvalues, the test is straightforward.  For complex
*> eigenvalues, E(j) and a(j) are complex, represented by
*> Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that
*> eigenvector becomes
*>
*>                 max( |Wr|, |Wi| )
*>     --------------------------------------------
*>     n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| )
*>
*> where
*>
*>     Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j)
*>
*>     Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j)
*>
*>                         T   T  _
*> For left eigenvectors, A , B , a, and b  are used.
*>
*> DGET52 also tests the normalization of E.  Each eigenvector is
*> supposed to be normalized so that the maximum "absolute value"
*> of its elements is 1, where in this case, "absolute value"
*> of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
*> maximum "absolute value" norm of a vector v  M(v).
*> if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate
*> vector.  The normalization test is:
*>
*>         RESULT(2) =      max       | M(v(j)) - 1 | / ( n ulp )
*>                    eigenvectors v(j)
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] LEFT
*> \verbatim
*>          LEFT is LOGICAL
*>          =.TRUE.:  The eigenvectors in the columns of E are assumed
*>                    to be *left* eigenvectors.
*>          =.FALSE.: The eigenvectors in the columns of E are assumed
*>                    to be *right* eigenvectors.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The size of the matrices.  If it is zero, DGET52 does
*>          nothing.  It must be at least zero.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is DOUBLE PRECISION array, dimension (LDA, N)
*>          The matrix A.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of A.  It must be at least 1
*>          and at least N.
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is DOUBLE PRECISION array, dimension (LDB, N)
*>          The matrix B.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of B.  It must be at least 1
*>          and at least N.
*> \endverbatim
*>
*> \param[in] E
*> \verbatim
*>          E is DOUBLE PRECISION array, dimension (LDE, N)
*>          The matrix of eigenvectors.  It must be O( 1 ).  Complex
*>          eigenvalues and eigenvectors always come in pairs, the
*>          eigenvalue and its conjugate being stored in adjacent
*>          elements of ALPHAR, ALPHAI, and BETA.  Thus, if a(j)/b(j)
*>          and a(j+1)/b(j+1) are a complex conjugate pair of
*>          generalized eigenvalues, then E(,j) contains the real part
*>          of the eigenvector and E(,j+1) contains the imaginary part.
*>          Note that whether E(,j) is a real eigenvector or part of a
*>          complex one is specified by whether ALPHAI(j) is zero or not.
*> \endverbatim
*>
*> \param[in] LDE
*> \verbatim
*>          LDE is INTEGER
*>          The leading dimension of E.  It must be at least 1 and at
*>          least N.
*> \endverbatim
*>
*> \param[in] ALPHAR
*> \verbatim
*>          ALPHAR is DOUBLE PRECISION array, dimension (N)
*>          The real parts of the values a(j) as described above, which,
*>          along with b(j), define the generalized eigenvalues.
*>          Complex eigenvalues always come in complex conjugate pairs
*>          a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent
*>          elements in ALPHAR, ALPHAI, and BETA.  Thus, if the j-th
*>          and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1)
*>          is assumed to be equal to ALPHAR(j)/BETA(j).
*> \endverbatim
*>
*> \param[in] ALPHAI
*> \verbatim
*>          ALPHAI is DOUBLE PRECISION array, dimension (N)
*>          The imaginary parts of the values a(j) as described above,
*>          which, along with b(j), define the generalized eigenvalues.
*>          If ALPHAI(j)=0, then the eigenvalue is real, otherwise it
*>          is part of a complex conjugate pair.  Complex eigenvalues
*>          always come in complex conjugate pairs a(j)/b(j) and
*>          a(j+1)/b(j+1), which are stored in adjacent elements in
*>          ALPHAR, ALPHAI, and BETA.  Thus, if the j-th and (j+1)-st
*>          eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to
*>          be equal to  -ALPHAI(j)/BETA(j).  Also, nonzero values in
*>          ALPHAI are assumed to always come in adjacent pairs.
*> \endverbatim
*>
*> \param[in] BETA
*> \verbatim
*>          BETA is DOUBLE PRECISION array, dimension (N)
*>          The values b(j) as described above, which, along with a(j),
*>          define the generalized eigenvalues.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is DOUBLE PRECISION array, dimension (N**2+N)
*> \endverbatim
*>
*> \param[out] RESULT
*> \verbatim
*>          RESULT is DOUBLE PRECISION array, dimension (2)
*>          The values computed by the test described above.  If A E or
*>          B E is likely to overflow, then RESULT(1:2) is set to
*>          10 / ulp.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup double_eig
*
*  =====================================================================
      SUBROUTINE DGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
     $                   ALPHAI, BETA, WORK, RESULT )
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      LOGICAL            LEFT
      INTEGER            LDA, LDB, LDE, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
     $                   B( LDB, * ), BETA( * ), E( LDE, * ),
     $                   RESULT( 2 ), WORK( * )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, TEN
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            ILCPLX
      CHARACTER          NORMAB, TRANS
      INTEGER            J, JVEC
      DOUBLE PRECISION   ABMAX, ACOEF, ALFMAX, ANORM, BCOEFI, BCOEFR,
     $                   BETMAX, BNORM, ENORM, ENRMER, ERRNRM, SAFMAX,
     $                   SAFMIN, SALFI, SALFR, SBETA, SCALE, TEMP1, ULP
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, DLANGE
      EXTERNAL           DLAMCH, DLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           DGEMV
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          ABS, DBLE, MAX
*     ..
*     .. Executable Statements ..
*
      RESULT( 1 ) = ZERO
      RESULT( 2 ) = ZERO
      IF( N.LE.0 )
     $   RETURN
*
      SAFMIN = DLAMCH( 'Safe minimum' )
      SAFMAX = ONE / SAFMIN
      ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
*
      IF( LEFT ) THEN
         TRANS = 'T'
         NORMAB = 'I'
      ELSE
         TRANS = 'N'
         NORMAB = 'O'
      END IF
*
*     Norm of A, B, and E:
*
      ANORM = MAX( DLANGE( NORMAB, N, N, A, LDA, WORK ), SAFMIN )
      BNORM = MAX( DLANGE( NORMAB, N, N, B, LDB, WORK ), SAFMIN )
      ENORM = MAX( DLANGE( 'O', N, N, E, LDE, WORK ), ULP )
      ALFMAX = SAFMAX / MAX( ONE, BNORM )
      BETMAX = SAFMAX / MAX( ONE, ANORM )
*
*     Compute error matrix.
*     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
*
      ILCPLX = .FALSE.
      DO 10 JVEC = 1, N
         IF( ILCPLX ) THEN
*
*           2nd Eigenvalue/-vector of pair -- do nothing
*
            ILCPLX = .FALSE.
         ELSE
            SALFR = ALPHAR( JVEC )
            SALFI = ALPHAI( JVEC )
            SBETA = BETA( JVEC )
            IF( SALFI.EQ.ZERO ) THEN
*
*              Real eigenvalue and -vector
*
               ABMAX = MAX( ABS( SALFR ), ABS( SBETA ) )
               IF( ABS( SALFR ).GT.ALFMAX .OR. ABS( SBETA ).GT.
     $             BETMAX .OR. ABMAX.LT.ONE ) THEN
                  SCALE = ONE / MAX( ABMAX, SAFMIN )
                  SALFR = SCALE*SALFR
                  SBETA = SCALE*SBETA
               END IF
               SCALE = ONE / MAX( ABS( SALFR )*BNORM,
     $                 ABS( SBETA )*ANORM, SAFMIN )
               ACOEF = SCALE*SBETA
               BCOEFR = SCALE*SALFR
               CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1,
     $                     ZERO, WORK( N*( JVEC-1 )+1 ), 1 )
               CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ),
     $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
            ELSE
*
*              Complex conjugate pair
*
               ILCPLX = .TRUE.
               IF( JVEC.EQ.N ) THEN
                  RESULT( 1 ) = TEN / ULP
                  RETURN
               END IF
               ABMAX = MAX( ABS( SALFR )+ABS( SALFI ), ABS( SBETA ) )
               IF( ABS( SALFR )+ABS( SALFI ).GT.ALFMAX .OR.
     $             ABS( SBETA ).GT.BETMAX .OR. ABMAX.LT.ONE ) THEN
                  SCALE = ONE / MAX( ABMAX, SAFMIN )
                  SALFR = SCALE*SALFR
                  SALFI = SCALE*SALFI
                  SBETA = SCALE*SBETA
               END IF
               SCALE = ONE / MAX( ( ABS( SALFR )+ABS( SALFI ) )*BNORM,
     $                 ABS( SBETA )*ANORM, SAFMIN )
               ACOEF = SCALE*SBETA
               BCOEFR = SCALE*SALFR
               BCOEFI = SCALE*SALFI
               IF( LEFT ) THEN
                  BCOEFI = -BCOEFI
               END IF
*
               CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1,
     $                     ZERO, WORK( N*( JVEC-1 )+1 ), 1 )
               CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ),
     $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
               CALL DGEMV( TRANS, N, N, BCOEFI, B, LDA, E( 1, JVEC+1 ),
     $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
*
               CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC+1 ),
     $                     1, ZERO, WORK( N*JVEC+1 ), 1 )
               CALL DGEMV( TRANS, N, N, -BCOEFI, B, LDA, E( 1, JVEC ),
     $                     1, ONE, WORK( N*JVEC+1 ), 1 )
               CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC+1 ),
     $                     1, ONE, WORK( N*JVEC+1 ), 1 )
            END IF
         END IF
   10 CONTINUE
*
      ERRNRM = DLANGE( 'One', N, N, WORK, N, WORK( N**2+1 ) ) / ENORM
*
*     Compute RESULT(1)
*
      RESULT( 1 ) = ERRNRM / ULP
*
*     Normalization of E:
*
      ENRMER = ZERO
      ILCPLX = .FALSE.
      DO 40 JVEC = 1, N
         IF( ILCPLX ) THEN
            ILCPLX = .FALSE.
         ELSE
            TEMP1 = ZERO
            IF( ALPHAI( JVEC ).EQ.ZERO ) THEN
               DO 20 J = 1, N
                  TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
   20          CONTINUE
               ENRMER = MAX( ENRMER, TEMP1-ONE )
            ELSE
               ILCPLX = .TRUE.
               DO 30 J = 1, N
                  TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
     $                    ABS( E( J, JVEC+1 ) ) )
   30          CONTINUE
               ENRMER = MAX( ENRMER, TEMP1-ONE )
            END IF
         END IF
   40 CONTINUE
*
*     Compute RESULT(2) : the normalization error in E.
*
      RESULT( 2 ) = ENRMER / ( DBLE( N )*ULP )
*
      RETURN
*
*     End of DGET52
*
      END

max_stars_repo_name: haampie/libflame

max_issues_repo_name: haampie/libflame

max_forks_repo_name: haampie/libflame


text:
!     This software is Copyright (c) 2012-2015 The Regents of the
!     University of California. All Rights Reserved.
!
!     Permission to copy and modify this software and its documentation
!     for educational, research and non-profit purposes, without fee,
!     and without a written agreement is hereby granted, provided that
!     the above copyright notice, this paragraph and the following three
!     paragraphs appear in all copies.
!
!     Permission to make commercial use of this software may be obtained
!     by contacting:
!
!     Technology Transfer Office
!     9500 Gilman Drive, Mail Code 0910
!     University of California
!     La Jolla, CA 92093-0910
!     (858) 534-5815
!     invent@ucsd.edu
!
!     This software program and documentation are copyrighted by The
!     Regents of the University of California. The software program and
!     documentation are supplied "as is", without any accompanying
!     services from The Regents. The Regents does not warrant that the
!     operation of the program will be uninterrupted or error-free. The
!     end-user understands that the program was developed for research
!     purposes and is advised not to rely exclusively on the program for
!     any reason.
!
!     IN NO EVENT SHALL THE UNIVERSITY OF CALIFORNIA BE LIABLE TO ANY
!     PARTY FOR DIRECT, INDIRECT, SPECIAL, INCIDENTAL, OR CONSEQUENTIAL
!     DAMAGES, INCLUDING LOST PROFITS, ARISING OUT OF THE USE OF THIS
!     SOFTWARE AND ITS DOCUMENTATION, EVEN IF THE UNIVERSITY OF
!     CALIFORNIA HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!     THE UNIVERSITY OF CALIFORNIA SPECIFICALLY DISCLAIMS ANY
!     WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
!     OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE
!     SOFTWARE PROVIDED HEREUNDER IS ON AN "AS IS" BASIS, AND THE
!     UNIVERSITY OF CALIFORNIA HAS NO OBLIGATIONS TO PROVIDE
!     MAINTENANCE, SUPPORT, UPDATES, ENHANCEMENTS, OR MODIFICATIONS.

      SUBROUTINE GE (N, A, B)

      IMPLICIT NONE

      INTEGER, INTENT(IN) :: N
      REAL*8, INTENT(IN) :: A(N,N)
      REAL*8, INTENT(INOUT) :: B(N)

      INTEGER m, ipv, i, j
      REAL*8 pivot, saveEl
      REAL*8, ALLOCATABLE :: C(:,:)

      IF (N .EQ. 2) THEN
         pivot  = A(1,1)*A(2,2) - A(2,1)*A(1,2)
         saveEl = (B(1)*A(2,2) - B(2)*A(1,2))/pivot
         B(2)   = (B(2)*A(1,1) - B(1)*A(2,1))/pivot
         B(1)   = saveEl
         RETURN
      END IF
      ALLOCATE(C(N,N+1))

      C(:N,:N) = A
      C(:,N+1) = B

      DO m=1,N-1
         ipv = m
         pivot = ABS(C(m,m))
         DO i=m+1,N
            IF (ABS(C(i,m)) .GT. pivot) THEN
               ipv = i
               pivot = ABS(C(i,m))
            END IF
         END DO
         IF (pivot .LT. 2D0*EPSILON(pivot)) THEN
            PRINT *, 'Singular matrix'
            STOP
         END IF
         IF (ipv .NE. m) THEN
            DO j=m, N+1
               saveEl = C(m,j)
               C(m,j) = C(ipv,j)
               C(ipv,j) = saveEl
            END DO
            DO j=1, m-1
            END DO
         END IF

         DO i=m+1,N
            saveEl = C(i,m)/C(m,m)
            C(i,m) = 0D0
            DO j=m+1,N+1
               C(i,j) = C(i,j) - saveEl*C(m,j)
            END DO
         END DO
      END DO

      DO j=N,1,-1
         DO i=j+1,N
            C(j,N+1) = C(j,N+1) - C(j,i)*C(i,N+1)
         END DO
         C(j,N+1) = C(j,N+1)/C(j,j)
      END DO

      B = C(:,N+1)

      DEALLOCATE(C)

      RETURN
      END SUBROUTINE GE

ext: f


text:
At the west end of the Arboretum there is a large wooden Address(gazebo, 38.529283, 121.762326) surrounded by gardens. It was built in 1966.

The gazebo is in the middle of the Carolee Shields White Flower Garden near the Shields Oak Grove.  The garden surrounding the Gazebo is a moon garden, meant to be viewed at night. All of the plants in that garden are whiteflowered varieties.

This gazebo and the grassy area surrounding the gazebo are great places to relax, study and enjoy the outdoors on a nice day. The setting is absolutely beautiful and it is probably a great location for a wedding  in fact, many people do host their wedding ceremonies out there. Info about reserving the Gazebo for a ceremony: http://events.ucdavis.edu/events/ Because this area is so far from the central part of campus, I dont think there is much traffic here, which is another reason why it is a great place to study.

Makes a pretty good Makeout Spots Makeout Spot if you can get it all to yourselves.  It is especially romantic early in the day when the Sun is just coming up since most of the flowers are particularly fragrant.  There is also a fairly low amount of traffic in the area that early.

The arboretum gazebo also makes an excellent place to chill with a guitar or harmonica. Recommended for acoustic rock and small, semiprivate shows.  The wooden structure allows for amazing acoustic sound, as well as allows for some privacy for our shyer musicians.

Map of it all located http://arboretum.ucdavis.edu/visitor_map.aspx on the official site.  As you can see, most of the Arboretum (and what most people whove seen it consider the nicer part) is the part furthest away from Campus campus and the Borders end.  Nice gardens, tons of flowers, many paths and lots of benches...its really nice. Just park on the street right by the Vet Med Teaching Hospital VMTHs Llamas llamas and enjoy. Users/EdwinS ES

There are often many burned out candles left on the Gazebo floor, likely left behind after some sort of ceremony.  This happens many times throughout the year. If you reserve it through Conference and Event Services for an event they will make sure it is powerwashed before you arrive.

 Strange things that can happen at the Arboretum Gazebo

The gazebo was once the site of a birthday party for http://www.livejournal.com/users/littlenick/ Nick Arora that turned into a brawl between several East Davis High School kids hopped up on scooby snacks and Users/KrisFricke Kris Frickes head. The offenders were vanquished by a group of attendees unofficially called the Vigilantes of the Central Valley. They fled; Kris Frickes head (along with the rest of him) survived in an unsurpassed feat of resilience and the masses scattered. When the officials arrived to attend to the wounded, the party ended. The keg had been emptied already anyway. Users/GautamJoshi

When I was a freshman in Tercero B, this is the place kids from the dorm came to drink most of the time.  Its a somewhat longish walk, but not too bad.  I remember hanging out here one night, and after people started to get tipsy, I walked off, circled around to the street, flipped on a big Mag Light, and walked loudly toward the gazebo, before saying in what I thought was a pretty good cop voice, OK, kids, Im going to need to see some ID.  Everyone but one guy goes dead silent.  That one guy?  Holy crap, I dont think Ive ever seen anyone run that fast in my life.  By the time everyone cracked up he was halfway back to the dorms, and when we got back, he started laughing at people who werent as fast as him and got (he thought) stuck with tickets. Users/TomGarberson

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE SOR(A,B,C,D,E,F,U,JMAX,RJAC)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(JMAX,JMAX),B(JMAX,JMAX),C(JMAX,JMAX),
     *    D(JMAX,JMAX),E(JMAX,JMAX),F(JMAX,JMAX),U(JMAX,JMAX)
      PARAMETER(MAXITS=1000,EPS=1.D-5,ZERO=0.D0,HALF=.5D0,QTR=.25D0,ONE=
     *1.D0)
      ANORMF=ZERO
      DO 12 J=2,JMAX-1
        DO 11 L=2,JMAX-1
          ANORMF=ANORMF+ABS(F(J,L))
11      CONTINUE
12    CONTINUE
      OMEGA=ONE
      DO 15 N=1,MAXITS
        ANORM=ZERO
        DO 14 J=2,JMAX-1
          DO 13 L=2,JMAX-1
            IF(MOD(J+L,2).EQ.MOD(N,2))THEN
              RESID=A(J,L)*U(J+1,L)+B(J,L)*U(J-1,L)+
     *            C(J,L)*U(J,L+1)+D(J,L)*U(J,L-1)+
     *            E(J,L)*U(J,L)-F(J,L)
              ANORM=ANORM+ABS(RESID)
              U(J,L)=U(J,L)-OMEGA*RESID/E(J,L)
            ENDIF
13        CONTINUE
14      CONTINUE
        IF(N.EQ.1) THEN
          OMEGA=ONE/(ONE-HALF*RJAC**2)
        ELSE
          OMEGA=ONE/(ONE-QTR*RJAC**2*OMEGA)
        ENDIF
        IF((N.GT.1).AND.(ANORM.LT.EPS*ANORMF))RETURN
15    CONTINUE
      PAUSE 'MAXITS exceeded'
      END

ext: for


text:
       program ph2dt

c Version 1.3 - 08/2010  same as 1.2, but now compiles with gfortran
c                         (rcs removed)
c Version 1.2 - 07/2010
c Version 1.1 - 10/2004
c Author: Felix Waldhauser, felixw@ldeo.columbia.edu
c
c started 03/1999
c 01-03/2001  clean up & bug fixes by Bruce Julian, Fred Klein, Keith
c             Richards-Dinger, Felix Waldhauser
c 10/2004     Version 1.1: fixed errors listed in BugList to V 1.0.
c 05/2005     accomodate format for neg magnitudes
c 07/2010     version 1.2
c 08/2010     now compiles with gfortran  (rcs removed and mod problem fixed)
c 08/2010     version 1.3
c
c Purpose:
c Reads and filters absolute travel-time data from network catalogs
c to form travel-time data for pairs of earthquakes. Control file or
c interactive input of parameters to optimize linkage between events
c (i.e. optimize quality and minimize number of links between events).
c See hypoDD user guide for information on how to use ph2dt.

c For a user guide to ph2dt see USGS open-file report:
c    Waldhauser, F., HypoDD: A computer program to compute double-difference
c    earthquake locations,  U.S. Geol. Surv. open-file report , 01-113,
c    Menlo Park, California, 2001.
c

c--Reads the input file specified om the command line (ph2dt ph2dt.inp)
c  to get the input station and phase file names, and the numerical parameters.
c--Writes the files (input files to hypoDD):
c  dt.ct	the travel times for the common stations for close
c               earthquake pairs
c  event.dat	earthquake list in dd format
c  event.sel	selected earthquake list in dd format (should be the same
c               as event.dat)
c  ph2dt.log	log file
c--Most earthquake selection (dates, magnitudes, etc) happens before ph2dt
c  processes phase to differential time data.
c--Use the files events.select to select events by ID (one per line).
c--FWK 2/9/01

c The standard input format (free format) is:
c#  yr  mn dy hr mn  sc      lat      lon      dep  mag  eh  ez  res  id
c sta  t-t  wght pha
c e.g.:
c# 1997  6 13 13 44 37.33  49.2424 -123.6192   3.42 3.43 0.0 0.0 0.0  1
c NAB      4.871   1    P
c BIB      5.043   0.5  P
c SHB      7.043   1    S
c WPB      8.934   1    P

c See hypoDD user guide for a description of the parameters.

	implicit none

	include 'ph2dt.inc'

	real		a_dist(MOBS)
	real		alat
	real		alon
	real		aoffs(MEV)
	real		a_time1(MOBS)
	real		a_time2(MOBS)
	doubleprecision	atoangle	! ASCII-to-angle function
	real		avoff
	real		avoff_str
	real		maxoff_str
	real		a_wtr(MOBS)
	real		az
	real		blat
	real		blon
	real		b_time1(MOBS)
	real		b_time2(MOBS)
	character	buf1*20		! Input buffer
	character	buf2*20		! Input buffer
	real		b_wtr(MOBS)
	integer		cuspid(MEV)
	integer		date(MEV)
	character	dattim*25	! Current time
	real		del
	real		depth(MEV)
	real		dist
	real		dlat
	real		dlon
	integer		dy
	logical		ex
	real		herr
	integer		hr(MEV)
	integer		i
	integer		iargc
	integer		icusp(MEV)
	integer		idata
	integer		iformat
	integer		ii
	integer		iimp
	integer		iindx(MOBS)
	integer		inb
	integer		indx(MEV)
	integer		iobs
	integer		iobsP
	integer		ipair
	integer		ipair_str
	integer		ista
	integer		itake
	integer		itmp
	integer		j
	integer		k
	integer		kk
        real            KMPERDEG
	integer		l
	real		lat(MEV)
	integer		limobs_pair
	integer		log
	real		lon(MEV)
	integer		m
	real		mag
	real		maxdist
	integer		maxobs_pair
	real		maxoffset
	integer		minute(MEV)
	integer		minobs_pair
	real		minwght
	integer		mnb
	integer		mo
	integer		n1
	integer		n2
	integer		n3
	integer		n4
	integer		n5
	integer		n6
	integer		n7
	integer		n8
	integer		narguments
	integer		ncusp
	integer		nerr
	integer		nev
	integer		nobs
	integer		nobs_ct(MEV)
	integer		npair
	integer		npha
	integer		nsta
	integer		ok
	real		PI
	real		p_time(MEV,MOBS)
	real		p_wghtr(MEV,MOBS)
	real		res
	real		rtime
	real		sec(MEV)
	real		s_lat(MSTA)
	real		s_lon(MSTA)
	integer		sscanf3		! String-reading function
	integer		trimlen
	real		vel
	real		verr
	real		wtr
	integer		yr
	character	a_lab(MOBS)*7
	character	a_pha(MOBS)*1
	character	b_lab(MOBS)*7
	character	b_pha(MOBS)*1
	character	fn0*80
	logical		fn0ex
	character	fn2*80
	logical		fn2ex
	character	fn9*80
	character	line*180
	character	p_pha(MEV,MOBS)*1
	character	p_sta(MEV,MSTA)*7
	character	s_lab(MSTA)*7
	character	str1*1
	character	str30*30
	character	str*40
	character	take(MEV,MEV)*1

	parameter	(PI=3.141593)
        parameter       (KMPERDEG=111.1949266)

c setpar:
c--Only standard format is now supported. Conversion from other formats occurs externally
      iformat= 0        ! 0=standard; 1=NCSN (non-y2k/y2k);
                        ! 2=pkf_hrsn (Bill's S-plus); 3= arcvel (stefano)
c file with cuspids to select for
      fn9= 'events.select'

      log= 20
      open(log,file='ph2dt.log',status='unknown')
      str= 'starting ph2dt (v1.3 - 08/2010)...'
      call datetime(dattim)
      write(6,'(a40,a)') str, dattim
      write(log,'(a40,a)') str, dattim

c--- get input parameter file name:
      narguments = iargc()
      if(narguments.lt.1) goto 2
      call getarg(1,str30)
      inquire(FILE= str30,exist=ex)
      if(.not. ex) stop' >>> ERROR OPENING INPUT PARAMETER FILE.'

c open/read input parameter file:
      open(1,status='unknown',file=str30(1:trimlen(str30)))
      idata= 2
      l=1
1     read (1,'(a)',end=111) line
      if(line(1:1).eq.'*' .or. line(2:2).eq.'*') goto 1
      if(l.eq.1) read (line,'(a)',err=112) fn2
      if(l.eq.2) read (line,'(a)',err=112) fn0

      if(l.eq.3) then
        read (line,*,err=112) minwght,maxdist,
     &         maxoffset,mnb,limobs_pair,minobs_pair,maxobs_pair
        close (1)
        goto 34
      endif
      l=l+1
      goto 1

c--Error statements
111   stop '** Premature end of command file'
112   write (*,*) '** Error reading data from input line:'
      write (*,'(a)') line
      stop

c--- Get modified (high-resolution) station file name:
2     write (6,'(a)') 'STATION FILE [<ret> station.dat]: '
      read (5,'(a)') fn2
      if (trimlen(fn2).le.1) then
         fn2= 'station.dat' !default input file name
      else
         fn2= fn2(1:trimlen(fn2))
      endif
      inquire (FILE= fn2,exist=fn2ex)
      if (.not. fn2ex) then
          stop' --> ERROR OPENING STATION FILE.'
      endif

c--- Get phase input file name:
      write (6,'(a)') 'PHASE INPUT FILE [<ret> phase.dat]: '
      read (5,'(a)') fn0
      if(trimlen(fn0).le.1) then
         fn0= 'phase.dat'       !default input file name
      else
         fn0= fn0(1:trimlen(fn0))
      endif
      inquire(FILE= fn0,exist=fn0ex)
      if(.not. fn0ex) then
          stop' --> ERROR OPENING PHASE FILE.'
      endif

c--- Set phase format: only standard format now supported
      iformat= 0

c--- get max pick quality
      write(6,'(a)') 'Min. pick weight/quality [<ret> 0]: '
      read(5,'(a)') str
      read(str,*,err=6)minwght
      goto 7
6     minwght= 0
7     continue

c--- Get max event-station distance
      write(6,'(a)') 'Max. distance between event pair and '//
     & 'station [<ret> 200 km]: '
      read(5,'(a)') str
      read(str,*,err=8)maxdist
      goto 9
8     maxdist= 200
9     continue

c--- Get max interevent offset (km) for which dtimes are calculated
      write(6,'(a)') 'Max. epicentral separation [<ret> 10 km]: '
      read(5,'(a)') str
      read(str,*,err=10)maxoffset
      goto 11
10    maxoffset= 10
11    continue

c--- Get number of nearest neighbours
      write(6,'(a)') 'Max. number of nearest neighbours per '//
     & 'event [<ret> 10]: '
      read(5,'(a)') str
      read(str,*,err=20)mnb
      goto 21
20    mnb= 10
21    continue

c--- Get number of dt-obs per pair so it counts as neighbour
      write(6,'(a)') 'Min. number of links necessary to '//
     & 'define a neighbor [<ret> 8]: '
      read(5,'(a)') str
      read(str,*,err=25)limobs_pair
      goto 26
25    limobs_pair= 8
26    continue

c--- Get min number of dt-obs per pair
      write(6,'(a)') 'Min. number of links saved per pair [<ret> 8]: '
      read(5,'(a)') str
      read(str,*,err=27)minobs_pair
      goto 28
27    minobs_pair=  8
28    continue

c--- Get max number of dt-obs per pair
      write(6,'(a)') 'Max. number of links saved per pair [<ret> 50]: '
      read(5,'(a)') str
      read(str,*,err=29)maxobs_pair
      goto 30
29    maxobs_pair= 50
30    continue

34    continue
      idata=2

c--- open files:
      open(12,file='dt.ct',status='unknown')
      open(14,file='event.dat',status='unknown')
      open(15,file='event.sel',status='unknown')

c--- read icusp's
      ncusp= 0
      inquire(FILE= fn9,exist=ex)
      if(ex) then
         open (1,file=fn9,status='unknown')
         i=1
35       read (1,'(a)',end=36,err=37) line
         read (line,*,err=37) icusp(i)
         i=i+1
         goto 35

c--Get out of read loop on error or end of file
37       write (6,*) '** Error reading id number:'
         write (6,*) line
         stop
36       ncusp= i-1
         close(1)
      endif

      write (*,'(/,"reading data ...")')
      write (log,'("reading data ...")')
c--- read stations
      open (2,file=fn2, status='unknown')
      i=1
40    read(2,'(a)',end=48)line

c       Split into fields separated by white space
	if (sscanf3(line, "%s%s%s", s_lab(i), buf1, buf2) .ne. 3) then
	  write (6,*) line
	  stop '** Bad station line'
	endif
	call rpad(s_lab(i))

c       Convert strings to numbers, interpreting colons, if any.
        s_lat(i) = atoangle(buf1)
        s_lon(i) = atoangle(buf2)

42      if (i.gt.MSTA) stop'>>> Increase MSTA in ph2dt.inc.'
        i= i+1
      goto 40
48    continue

      nsta= i-1
      write (*,*) '> stations = ', nsta
      write (log,*) '> stations = ', nsta
      close(2)

c--- check for double station:
      do i=1,nsta-1
        do j=i+1,nsta
           if (s_lab(i).eq.s_lab(j)) then
               write (*,*) s_lab(i)
               stop'>>> This station is listed twice in station file!'
           endif
        enddo
      enddo

      if(idata.eq.2) goto 666
666   continue	!no cross data

c--- read absolute network travel times:
      open (2,file=fn0, status='unknown')
      i= 1
      ii= 1
      npha= 0
100   continue
c--Inoperable code for reading other formats removed by FWK
      goto 131   !standard format

c--Fatal error on read messages, added by FWK
1290  write (6,*) line
      write (*,'(a)') line
      stop '** Bad earthquake line'
1292  write (6,*) line
      write (*,'(a)') line
      stop '** Bad phase line'

c--- start standard data format
c read header:
130   read (2,'(a)',end=159) line  		! read header line
      if (line(1:1).eq.'#') goto 160 		! store previous event
131   if (ii.eq.1) read(2,'(a)',end=200) line  	! read header line
      if (line(1:1).eq.'#') then
         read (line,*,err=1290) str1,yr,mo,dy,hr(i),minute(i),sec(i),
     &   lat(i),lon(i),depth(i),mag,herr,verr,res,cuspid(i)
         date(i)= yr*10000 + mo*100 + dy
         rtime= hr(i)*1000000 + minute(i)*10000 + sec(i)*100

         ii= ii+1
         k= 1  !phase counter
      else
c read phase data lines
         read (line,*,err=1292) p_sta(i,k),p_time(i,k),p_wghtr(i,k),
     &   p_pha(i,k)
         if (p_wghtr(i,k).ge.0.and.p_wghtr(i,k).lt.minwght) goto 130
         k= k+1
         if (k.gt.MOBS) stop'>>> Increase MOBS ph2dt.inc!'
      endif
      goto 130
c--- end standard format

c--- processing for all formats starts here:
159   if(k.eq.1) goto 200
160   nobs_ct(i)=k-1
      itake= 1
      if (ncusp.gt.0) then
         itake= 0
         do k=1,ncusp
            if (cuspid(i).eq.icusp(k)) itake= 1
         enddo
      endif
c--- event selection
c write header to total event list file:
      write(14,612)date(i),int(rtime),lat(i),
     &            lon(i),depth(i),mag,herr,verr,res,
     &            cuspid(i)

c--Keep an event only if it has more than minobs_pair observations,
c  and is on ID list if one was given
c--Faulty selection code removed by FWK
      if (itake.eq.1 .and. nobs_ct(i).ge.minobs_pair) then

c write event to selected event list file:
         write (15,612) date(i),int(rtime),lat(i),
     &            lon(i),depth(i),mag,herr,verr,res,
     &            cuspid(i)
         npha= npha+nobs_ct(i)
         i= i+1
         if (i.gt.MEV) stop'>>> Increase MEV in ph2dt.inc!'
      endif

cfw612   format (i8,2x,i8,2x,f8.4,2x,f9.4,2x,
cfw     &       f9.3,2x,f3.1,2x,f6.2,2x,f6.2,2x,f5.2,1x,i10)
612   format (i8,2x,i8,2x,f8.4,2x,f9.4,2x,
     &       f9.3,2x,f4.1,2x,f6.2,2x,f6.2,2x,f5.2,1x,i10)

      if (iformat.eq.0.and.line(1:1).ne.'#') goto 200 	!no ne event
      goto 100

200   nev= i-1
      write(*,*)'> events total = ',ii-1
      write(*,*)'> events selected = ',nev
      write(*,*)'> phases = ',npha
      write(log,*)'> events total = ',ii-1
      write(log,*)'> events selected = ',nev
      write(log,*)'> phases = ',npha

c--- form dtimes:
      write(*,'("forming dtimes...")')
      write(log,'("forming dtimes...")')
      write(log,*)'Reporting missing stations (STA) and '
      write(log,*)'   outliers (STA,ID1,ID2,OFFSET (km),T1,T2,T1-T2):'
      n1= 0
      n2= 0
      n3= 0
      n4= 0
      n5= 0
      n6= 0
      n7= 0
      n8= 0
      nerr= 0

      do i=1,nev
         do j=1,nev
           take(i,j)= '1'
         enddo
      enddo

      ipair= 1
      ipair_str= 1
      avoff= 0
      avoff_str= 0
      maxoff_str= 0
      do i=1,nev

c find nearest neighb:
        do j=1,nev
          dlat= lat(i) - lat(j)
          dlon= lon(i) - lon(j)
          aoffs(j)= sqrt( (dlat*KMPERDEG)**2 +
     &                    (dlon*(cos(lat(i)*PI/180)*KMPERDEG))**2 +
     &                 (depth(i)-depth(j))**2)
        enddo
        aoffs(i)= 99999
        call INDEXX(nev,aoffs,INDX)

        inb= 0
        nobs= 0
        do m=1,nev-1   ! same event last
          if(inb.ge.mnb) goto 400    ! next event

          k= indx(m)	! get nearest first
          if(take(k,i).eq.'0') then	! already selec as strong neighb
             inb= inb+1
             goto 350
          elseif(take(k,i).eq.'9') then	! already selec as weak neighb
             goto 350
          endif

          n1=n1+1

c check for max interevent offset:
          if(aoffs(indx(m)).gt.maxoffset) goto 400	! next event

c search for common stations/phases:
          iobs= 0
          iobsP= 0
          iimp= 0   !obs needs to be included regardless weight or dist
          do j=1,nobs_ct(i)
            do l=1,nobs_ct(k)
               if(p_sta(i,j).eq.p_sta(k,l).and.
     &            p_pha(i,j).eq.p_pha(k,l)) then

                  if(p_pha(i,j).eq.'P') n3= n3+1
                  if(p_pha(i,j).eq.'S') n6= n6+1
c check for station label in station file:
                  do ii= 1,nsta
                    ok= 0
                    if(iformat.eq.1) then   !moved station labels..
                      if( (p_sta(i,j)(6:6).eq.' '.and.
     &                  p_sta(i,j)(1:5).eq.s_lab(ii)(1:5)) .or.
     &                 (p_sta(i,j)(6:6).ne.' '.and.
     &                  p_sta(i,j)(1:6).eq.s_lab(ii)(1:6))) ok= 1
                    else
                      if(p_sta(i,j).eq.s_lab(ii)) ok= 1
                    endif
c select for station - pair centroid  distance
                    if(ok.eq.1) then
                         alat= s_lat(ii)
                         alon= s_lon(ii)
                         blat= (lat(i)+lat(k))/2
                         blon= (lon(i)+lon(k))/2
	                 call delaz(alat,alon,blat,blon,del,dist,az)

c delete far away stations
                         if(dist.gt.maxdist) then
                             n5= n5+1
                             goto 300
                         endif
                         goto 250
                    endif
                  enddo
                  n4= n4+1
                  write(log,*)'Station not in station file: ',p_sta(i,j)
                  goto 300
250               continue
                  ista= ii

c get average weight: 8/16/00
                   wtr = (abs(p_wghtr(i,j)) + abs(p_wghtr(k,l)))/2

c remove outliers above the separation-delaytime line:
                  if(p_pha(i,j).eq.'P') vel= 4.
                  if(p_pha(i,j).eq.'S') vel= 2.3
                  if(abs(p_time(i,j)-p_time(k,l)).gt.
     &               aoffs(indx(m))/vel + 0.5) then
                     write(log,'(a,a7,2i9,4f9.3)')'Outlier: ',
     & p_sta(i,j),cuspid(i),cuspid(k),aoffs(indx(m)),p_time(i,j),
     & p_time(k,l),p_time(i,j)-p_time(k,l)
                     nerr= nerr+1
                     goto 300
                  endif

                  iobs=iobs+1
                  if(p_pha(i,j).eq.'P') iobsP= iobsP+1
                  nobs=nobs+1
                  a_lab(iobs)= s_lab(ista)
                  a_time1(iobs)= p_time(i,j)
                  a_time2(iobs)= p_time(k,l)
                  a_wtr(iobs)= wtr
                  a_dist(iobs)= dist     !distance to station
                  a_pha(iobs)= p_pha(i,j)

                  if(p_wghtr(i,j).lt.0 .or. p_wghtr(k,l).lt.0) then
                     a_dist(iobs)= 0  	! set to 0 so it will be selected first
                     iimp= iimp+1
                  endif
                  goto 300
               endif
            enddo
300         continue   ! next station / observation
          enddo

          itmp= iobs
          if(iobs.gt.maxobs_pair) itmp= min(maxobs_pair+iimp,iobs)   !add important ones
          if(iobs.ge.minobs_pair) then  !pair selected

c sort obs by distance:
              if(iobs.gt.1) then
                call INDEXX(iobs,a_dist,iindx)
                do kk=1,iobs
                    b_lab(kk)=   a_lab(iindx(kk))
                    b_time1(kk)= a_time1(iindx(kk))
                    b_time2(kk)= a_time2(iindx(kk))
                    b_wtr(kk)=   a_wtr(iindx(kk))
                    b_pha(kk)=   a_pha(iindx(kk))
                enddo
              else
                b_lab(1)=   a_lab(1)
                b_time1(1)= a_time1(1)
                b_time2(1)= a_time2(1)
                b_wtr(1)=   a_wtr(1)
                b_pha(1)=   a_pha(1)
              endif

c write out delay times:
              write(12,'(a1,1x,i9,1x,i9)')'#',cuspid(i),cuspid(k)
              do kk= 1,itmp
                 write(12,'(a7,1x,f7.3,1x,f7.3,1x,f6.4,1x,a1)')
     &           b_lab(kk)(1:7),b_time1(kk),b_time2(kk), b_wtr(kk),
     &           b_pha(kk)
                 if(b_pha(kk).eq.'P') n7= n7+1
                 if(b_pha(kk).eq.'S') n8= n8+1
              enddo
              avoff= avoff + aoffs(k)
              ipair= ipair+1
          endif

          if(iobs.ge.limobs_pair) then	! select as strong neighbor
             take(i,k)= '0'
             inb= inb+1
             ipair_str= ipair_str+1
             avoff_str= avoff_str + aoffs(k)
             if(aoffs(k).gt.maxoff_str)maxoff_str= aoffs(k)
          else				! weak neighbor
             take(i,k)= '9'
          endif

350       continue  ! next neighbour
        enddo
400     continue  ! next event
        if(inb.lt.mnb) n2=n2+1
      enddo
      npair= ipair-1
      avoff= avoff/npair
      avoff_str= avoff_str/(ipair_str-1)

      write(*,*)'> P-phase pairs total = ',n3
      write(*,*)'> S-phase pairs total = ',n6
      write(*,*)'> outliers = ',nerr,' (',nerr*100/(n3+n6),'%)'
      write(*,*)'> phases at stations not in station list = ',n4
      write(*,*)'> phases at distances larger than MAXDIST = ',n5
      if(n3.gt.0)
     & write(*,*)'> P-phase pairs selected = ',n7,' (',n7*100/n3,'%)'
      if(n6.gt.0)
     & write(*,*)'> S-phase pairs selected = ',n8,' (',n8*100/n6,'%)'
      write(*,*)'> weakly linked events = ',n2,' (',n2*100/nev,'%)'
      write(*,*)'> linked event pairs = ',ipair
      write(*,*)'> average links per pair = ',(n7+n8)/ipair
      write(*,*)'> average offset (km) betw. linked events = ',
     & avoff
      write(*,*)'> average offset (km) betw. strongly linked events = ',
     & avoff_str
      write(*,*)'> maximum offset (km) betw. strongly linked events = ',
     & maxoff_str
      write(log,*)'> P-phase pairs total = ',n3
      write(log,*)'> S-phase pairs total = ',n6
      write(log,*)'> outliers = ',nerr,' (',nerr*100/(n3+n6),'%)'
      write(log,*)'> phases at stations not in station list = ',n4
      write(log,*)'> phases at distances larger than MAXDIST = ',n5
      if(n3.gt.0)
     & write(log,*)'> P-phase pairs selected = ',n7,' (',n7*100/n3,'%)'
      if(n6.gt.0)
     & write(log,*)'> S-phase pairs selected = ',n8,' (',n8*100/n6,'%)'
      write(log,*)'> weakly linked events = ',n2,' (',n2*100/nev,'%)'
      write(log,*)'> linked event pairs = ',ipair
      write(log,*)'> average links per pair = ',(n7+n8)/ipair
      write(log,*)'> average offset (km) betw. linked events = ',
     & avoff
      write(log,*)'> avg. offset (km) betw. strongly linked events = ',
     & avoff_str
      write(log,*)'> max. offset (km) betw. strongly linked events = ',
     & maxoff_str

      close(2)
      close(12)
      close(14)
      close(15)

	call datetime(dattim)

	write (log,'(/,"Done.  ",a)') dattim
	close (log)

	write (*,'(/,"Done.  ",a)') dattim

      write(*,'(/,a)')'Output files: dt.ct; event.dat; '//
     & 'event.sel; ph2dt.log'
      write(*,'(a)')'ph2dt parameters were: '
      write(*,*)'(minwght,maxdist,maxsep,maxngh,minlnk,'//
     & 'minobs,maxobs)'
      write(*,*)minwght,maxdist,maxoffset,mnb,limobs_pair,minobs_pair,
     & maxobs_pair

      end


text:
MODULE ex3_mod
CONTAINS
ATTRIBUTES(GLOBAL) SUBROUTINE sumation(a,b,d)
IMPLICIT NONE
INTEGER,INTENT(in) :: a(:,:)
INTEGER,INTENT(in) :: b(:,:)
INTEGER,INTENT(out):: d(:,:)
INTEGER :: i,j,nx,ny

i = blockDim%x*(blockIdx%x-1)+threadIdx%x
j = blockDim%y*(blockIdx%y-1)+threadIdx%y
nx=size(a,1) ; ny=size(a,2)
IF ((i<=nx).and.(j<=ny)) THEN
	d(i,j) = a(i,j)+ b(i,j)
END IF
END SUBROUTINE sumation
END MODULE ex3_mod

!!==========================================
PROGRAM ex3_presentation
USE cudafor
USE ex3_mod
IMPLICIT NONE
INTEGER, PARAMETER :: nx=30000000,ny=4
INTEGER :: a(nx,ny),b(nx,ny),d(nx,ny)
INTEGER, device :: a_d(nx,ny),b_d(nx,ny),d_d(nx,ny)
TYPE(dim3) :: gridDim,blockDim
a=10
b=20
blockDim=dim3(512,2,1)
gridDim=dim3(ceiling(real(nx)/blockDim%x),ceiling(real(ny)/blockDim%y),1)

a_d=a ; b_d=b
CALL sumation <<<gridDim,blockDim>>>(a_d,b_d,d_d)
d=d_d
IF(all(d==30)) THEN
      WRITE(* ,*)'**** OK ****'
ELSE
      WRITE(* ,*)'Program Failed'
END IF
END PROGRAM ex3_presentation

max_stars_repo_name: danhtaihoang/cudafortran

max_issues_repo_name: danhtaihoang/cudafortran

max_forks_repo_name: danhtaihoang/cudafortran


text:
!
!     ------------------------------------------------------------------
!     B N D I N V
!     ------------------------------------------------------------------
!
      SUBROUTINE BNDINV(A, EL, N, DETERM, EPSIL, ITEST, NSIZE)
!
!       DOUBLE PRECISION MATRIX INVERSION SUBROUTINE
!       FROM "DLYTAP".
!
!*      DOUBLE PRECISION E,F
!*      DOUBLE PRECISION A,EL,D,DSQRT,C,S,DETERP
!-----------------------------------------------
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:08:49   1/ 6/07
!...Modified by Charlotte Froese Fischer
!                     Gediminas Gaigalas  10/05/17
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: N
      INTEGER , INTENT(OUT) :: ITEST
      INTEGER , INTENT(IN) :: NSIZE
      REAL(DOUBLE) , INTENT(OUT) :: DETERM
      REAL(DOUBLE) , INTENT(IN) :: EPSIL
      REAL(DOUBLE) , INTENT(INOUT) :: A(NSIZE,NSIZE)
      REAL(DOUBLE) , INTENT(INOUT) :: EL(NSIZE,NSIZE)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: ISL2, K000FX, INDSNL, I, J, N1, M, K, J1, I1, KS
      REAL(DOUBLE) :: D, C, S, DETERP, F, E, EPSILP, RAT
!-----------------------------------------------
      IF (N < 2) GO TO 140
      ISL2 = 0
      K000FX = 2
      IF (ISL2 == 0) INDSNL = 2
      IF (ISL2 == 1) INDSNL = 1
!       CALL SLITET(2,INDSNL)
!       CALL OVERFL(K000FX)
!       CALL DVCHK(K000FX)
!
!       SET EL = IDENTITY MATRIX
      DO I = 1, N
         EL(I,:N) = 0.0D0
         EL(I,I) = 1.0D0
      END DO
!
!       TRIANGULARIZE A, FORM EL
!
      N1 = N - 1
      M = 2
      DO J = 1, N1
         DO I = M, N
            IF (A(I,J) == 0.0D0) CYCLE
            D = DSQRT(A(J,J)*A(J,J)+A(I,J)*A(I,J))
            C = A(J,J)/D
            S = A(I,J)/D
            DO K = J, N
               D = C*A(J,K) + S*A(I,K)
               A(I,K) = C*A(I,K) - S*A(J,K)
               A(J,K) = D
            END DO
            DO K = 1, N
               D = C*EL(J,K) + S*EL(I,K)
               EL(I,K) = C*EL(I,K) - S*EL(J,K)
               EL(J,K) = D
            END DO
         END DO
         M = M + 1
      END DO
!       CALL OVERFL(K000FX)
!       GO TO (140,51),K000FX
!
!       CALCULATE THE DETERMINANT
      DETERP = A(1,1)
      DO I = 2, N
         DETERP = DETERP*A(I,I)
      END DO
      DETERM = DETERP
!       CALL OVERFL(K000FX)
!       GO TO (140,520,520),K000FX
!
!       IS MATRIX SINGULAR
      F = A(1,1)
      E = A(1,1)
      DO I = 2, N
         IF (DABS(F) < DABS(A(I,I))) F = A(I,I)
         IF (DABS(E) <= DABS(A(I,I))) CYCLE
         E = A(I,I)
      END DO
      EPSILP = EPSIL
      IF (EPSILP <= 0) EPSILP = 1.0E-8
      RAT = E/F
      IF (ABS(RAT) < EPSILP) GO TO 130
!
!       INVERT TRIANGULAR MATRIX
      J = N
      DO J1 = 1, N
!       CALL SLITE(2)
         I = J
         ISL2 = 1
         DO I1 = 1, J
!       CALL SLITET(2,K000FX)
            IF (ISL2 == 0) K000FX = 2
            IF (ISL2 == 1) THEN
               K000FX = 1
               ISL2 = 0
            ENDIF
            SELECT CASE (K000FX)
            CASE DEFAULT
               A(I,J) = 1.0D0/A(I,I)
            CASE (2)
               KS = I + 1
               D = 0.0D0
               D = SUM(A(I,KS:J)*A(KS:J,J))
               A(I,J) = -D/A(I,I)
            END SELECT
 1003       CONTINUE
            I = I - 1
         END DO
         J = J - 1
      END DO
!       CALL OVERFL(K000FX)
!       GO TO (140,103,103),K000FX

!103    CALL DVCHK(K000FX)
!       GO TO (140,105),K000FX
!
!       PREMULTIPLY EL BY INVERTED TRIANGULAR MATRIX
      M = 1
      DO I = 1, N
         DO J = 1, N
            D = 0.0D0
            D = SUM(A(I,M:N)*EL(M:N,J))
            EL(I,J) = D
         END DO
         M = M + 1
      END DO
!       CALL OVERFL(K000FX)
!       GO TO (140,123,123),K000FX
!
!       RECOPY EL TO A
      A(:N,:N) = EL(:N,:N)
      ITEST = 0
!126    IF(INDSNL.EQ.1)CALL SLITE(2)
  126 CONTINUE
      IF (INDSNL == 1) ISL2 = 1
      RETURN
!
  130 CONTINUE
      ITEST = 1
      GO TO 126
  140 CONTINUE
      ITEST = -1
      GO TO 126
      END SUBROUTINE BNDINV

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


!*==zggev3.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief  ZGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download ZGGEV3 + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE ZGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
!                          VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          JOBVL, JOBVR
!       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION   RWORK( * )
!       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
!      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
!      $                   WORK( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> ZGGEV3 computes for a pair of N-by-N complex nonsymmetric matrices
!> (A,B), the generalized eigenvalues, and optionally, the left and/or
!> right generalized eigenvectors.
!>
!> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
!> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
!> singular. It is usually represented as the pair (alpha,beta), as
!> there is a reasonable interpretation for beta=0, and even for both
!> being zero.
!>
!> The right generalized eigenvector v(j) corresponding to the
!> generalized eigenvalue lambda(j) of (A,B) satisfies
!>
!>              A * v(j) = lambda(j) * B * v(j).
!>
!> The left generalized eigenvector u(j) corresponding to the
!> generalized eigenvalues lambda(j) of (A,B) satisfies
!>
!>              u(j)**H * A = lambda(j) * u(j)**H * B
!>
!> where u(j)**H is the conjugate-transpose of u(j).
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] JOBVL
!> \verbatim
!>          JOBVL is CHARACTER*1
!>          = 'N':  do not compute the left generalized eigenvectors;
!>          = 'V':  compute the left generalized eigenvectors.
!> \endverbatim
!>
!> \param[in] JOBVR
!> \verbatim
!>          JOBVR is CHARACTER*1
!>          = 'N':  do not compute the right generalized eigenvectors;
!>          = 'V':  compute the right generalized eigenvectors.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrices A, B, VL, and VR.  N >= 0.
!> \endverbatim
!>
!> \param[in,out] A
!> \verbatim
!>          A is COMPLEX*16 array, dimension (LDA, N)
!>          On entry, the matrix A in the pair (A,B).
!>          On exit, A has been overwritten.
!> \endverbatim
!>
!> \param[in] LDA
!> \verbatim
!>          LDA is INTEGER
!>          The leading dimension of A.  LDA >= max(1,N).
!> \endverbatim
!>
!> \param[in,out] B
!> \verbatim
!>          B is COMPLEX*16 array, dimension (LDB, N)
!>          On entry, the matrix B in the pair (A,B).
!>          On exit, B has been overwritten.
!> \endverbatim
!>
!> \param[in] LDB
!> \verbatim
!>          LDB is INTEGER
!>          The leading dimension of B.  LDB >= max(1,N).
!> \endverbatim
!>
!> \param[out] ALPHA
!> \verbatim
!>          ALPHA is COMPLEX*16 array, dimension (N)
!> \endverbatim
!>
!> \param[out] BETA
!> \verbatim
!>          BETA is COMPLEX*16 array, dimension (N)
!>          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
!>          generalized eigenvalues.
!>
!>          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
!>          underflow, and BETA(j) may even be zero.  Thus, the user
!>          should avoid naively computing the ratio alpha/beta.
!>          However, ALPHA will be always less than and usually
!>          comparable with norm(A) in magnitude, and BETA always less
!>          than and usually comparable with norm(B).
!> \endverbatim
!>
!> \param[out] VL
!> \verbatim
!>          VL is COMPLEX*16 array, dimension (LDVL,N)
!>          If JOBVL = 'V', the left generalized eigenvectors u(j) are
!>          stored one after another in the columns of VL, in the same
!>          order as their eigenvalues.
!>          Each eigenvector is scaled so the largest component has
!>          abs(real part) + abs(imag. part) = 1.
!>          Not referenced if JOBVL = 'N'.
!> \endverbatim
!>
!> \param[in] LDVL
!> \verbatim
!>          LDVL is INTEGER
!>          The leading dimension of the matrix VL. LDVL >= 1, and
!>          if JOBVL = 'V', LDVL >= N.
!> \endverbatim
!>
!> \param[out] VR
!> \verbatim
!>          VR is COMPLEX*16 array, dimension (LDVR,N)
!>          If JOBVR = 'V', the right generalized eigenvectors v(j) are
!>          stored one after another in the columns of VR, in the same
!>          order as their eigenvalues.
!>          Each eigenvector is scaled so the largest component has
!>          abs(real part) + abs(imag. part) = 1.
!>          Not referenced if JOBVR = 'N'.
!> \endverbatim
!>
!> \param[in] LDVR
!> \verbatim
!>          LDVR is INTEGER
!>          The leading dimension of the matrix VR. LDVR >= 1, and
!>          if JOBVR = 'V', LDVR >= N.
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
!>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
!> \endverbatim
!>
!> \param[in] LWORK
!> \verbatim
!>          LWORK is INTEGER
!>          The dimension of the array WORK.
!>
!>          If LWORK = -1, then a workspace query is assumed; the routine
!>          only calculates the optimal size of the WORK array, returns
!>          this value as the first entry of the WORK array, and no error
!>          message related to LWORK is issued by XERBLA.
!> \endverbatim
!>
!> \param[out] RWORK
!> \verbatim
!>          RWORK is DOUBLE PRECISION array, dimension (8*N)
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          =1,...,N:
!>                The QZ iteration failed.  No eigenvectors have been
!>                calculated, but ALPHA(j) and BETA(j) should be
!>                correct for j=INFO+1,...,N.
!>          > N:  =N+1: other then QZ iteration failed in DHGEQZ,
!>                =N+2: error return from DTGEVC.
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date January 2015
!
!> \ingroup complex16GEeigen
!
!  =====================================================================
      SUBROUTINE ZGGEV3(Jobvl,Jobvr,N,A,Lda,B,Ldb,Alpha,Beta,Vl,Ldvl,Vr,&
     &                  Ldvr,Work,Lwork,Rwork,Info)
      IMPLICIT NONE
!*--ZGGEV3220
!
!  -- LAPACK driver routine (version 3.6.1) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     January 2015
!
!     .. Scalar Arguments ..
      CHARACTER Jobvl , Jobvr
      INTEGER Info , Lda , Ldb , Ldvl , Ldvr , Lwork , N
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION Rwork(*)
      COMPLEX*16 A(Lda,*) , Alpha(*) , B(Ldb,*) , Beta(*) , Vl(Ldvl,*) ,&
     &           Vr(Ldvr,*) , Work(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ZERO , ONE
      PARAMETER (ZERO=0.0D0,ONE=1.0D0)
      COMPLEX*16 CZERO , CONE
      PARAMETER (CZERO=(0.0D0,0.0D0),CONE=(1.0D0,0.0D0))
!     ..
!     .. Local Scalars ..
      LOGICAL ilascl , ilbscl , ilv , ilvl , ilvr , lquery
      CHARACTER chtemp
      INTEGER icols , ierr , ihi , ijobvl , ijobvr , ileft , ilo , in , &
     &        iright , irows , irwrk , itau , iwrk , jc , jr , lwkopt
      DOUBLE PRECISION anrm , anrmto , bignum , bnrm , bnrmto , eps ,   &
     &                 smlnum , temp
      COMPLEX*16 x
!     ..
!     .. Local Arrays ..
      LOGICAL ldumma(1)
!     ..
!     .. External Subroutines ..
      EXTERNAL DLABAD , XERBLA , ZGEQRF , ZGGBAK , ZGGBAL , ZGGHD3 ,    &
     &         ZHGEQZ , ZLACPY , ZLASCL , ZLASET , ZTGEVC , ZUNGQR ,    &
     &         ZUNMQR
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      DOUBLE PRECISION DLAMCH , ZLANGE
      EXTERNAL LSAME , DLAMCH , ZLANGE
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , DBLE , DIMAG , MAX , SQRT
!     ..
!     .. Statement Functions ..
      DOUBLE PRECISION ABS1
!     ..
!     .. Statement Function definitions ..
      ABS1(x) = ABS(DBLE(x)) + ABS(DIMAG(x))
!     ..
!     .. Executable Statements ..
!
!     Decode the input arguments
!
      IF ( LSAME(Jobvl,'N') ) THEN
         ijobvl = 1
         ilvl = .FALSE.
      ELSEIF ( LSAME(Jobvl,'V') ) THEN
         ijobvl = 2
         ilvl = .TRUE.
      ELSE
         ijobvl = -1
         ilvl = .FALSE.
      ENDIF
!
      IF ( LSAME(Jobvr,'N') ) THEN
         ijobvr = 1
         ilvr = .FALSE.
      ELSEIF ( LSAME(Jobvr,'V') ) THEN
         ijobvr = 2
         ilvr = .TRUE.
      ELSE
         ijobvr = -1
         ilvr = .FALSE.
      ENDIF
      ilv = ilvl .OR. ilvr
!
!     Test the input arguments
!
      Info = 0
      lquery = (Lwork==-1)
      IF ( ijobvl<=0 ) THEN
         Info = -1
      ELSEIF ( ijobvr<=0 ) THEN
         Info = -2
      ELSEIF ( N<0 ) THEN
         Info = -3
      ELSEIF ( Lda<MAX(1,N) ) THEN
         Info = -5
      ELSEIF ( Ldb<MAX(1,N) ) THEN
         Info = -7
      ELSEIF ( Ldvl<1 .OR. (ilvl .AND. Ldvl<N) ) THEN
         Info = -11
      ELSEIF ( Ldvr<1 .OR. (ilvr .AND. Ldvr<N) ) THEN
         Info = -13
      ELSEIF ( Lwork<MAX(1,2*N) .AND. .NOT.lquery ) THEN
         Info = -15
      ENDIF
!
!     Compute workspace
!
      IF ( Info==0 ) THEN
         CALL ZGEQRF(N,N,B,Ldb,Work,Work,-1,ierr)
         lwkopt = MAX(1,N+INT(Work(1)))
         CALL ZUNMQR('L','C',N,N,N,B,Ldb,Work,A,Lda,Work,-1,ierr)
         lwkopt = MAX(lwkopt,N+INT(Work(1)))
         IF ( ilvl ) THEN
            CALL ZUNGQR(N,N,N,Vl,Ldvl,Work,Work,-1,ierr)
            lwkopt = MAX(lwkopt,N+INT(Work(1)))
         ENDIF
         IF ( ilv ) THEN
            CALL ZGGHD3(Jobvl,Jobvr,N,1,N,A,Lda,B,Ldb,Vl,Ldvl,Vr,Ldvr,  &
     &                  Work,-1,ierr)
            lwkopt = MAX(lwkopt,N+INT(Work(1)))
            CALL ZHGEQZ('S',Jobvl,Jobvr,N,1,N,A,Lda,B,Ldb,Alpha,Beta,Vl,&
     &                  Ldvl,Vr,Ldvr,Work,-1,Rwork,ierr)
            lwkopt = MAX(lwkopt,N+INT(Work(1)))
         ELSE
            CALL ZGGHD3(Jobvl,Jobvr,N,1,N,A,Lda,B,Ldb,Vl,Ldvl,Vr,Ldvr,  &
     &                  Work,-1,ierr)
            lwkopt = MAX(lwkopt,N+INT(Work(1)))
            CALL ZHGEQZ('E',Jobvl,Jobvr,N,1,N,A,Lda,B,Ldb,Alpha,Beta,Vl,&
     &                  Ldvl,Vr,Ldvr,Work,-1,Rwork,ierr)
            lwkopt = MAX(lwkopt,N+INT(Work(1)))
         ENDIF
         Work(1) = DCMPLX(lwkopt)
      ENDIF
!
      IF ( Info/=0 ) THEN
         CALL XERBLA('ZGGEV3 ',-Info)
         RETURN
      ELSEIF ( lquery ) THEN
         RETURN
      ENDIF
!
!     Quick return if possible
!
      IF ( N==0 ) RETURN
!
!     Get machine constants
!
      eps = DLAMCH('E')*DLAMCH('B')
      smlnum = DLAMCH('S')
      bignum = ONE/smlnum
      CALL DLABAD(smlnum,bignum)
      smlnum = SQRT(smlnum)/eps
      bignum = ONE/smlnum
!
!     Scale A if max element outside range [SMLNUM,BIGNUM]
!
      anrm = ZLANGE('M',N,N,A,Lda,Rwork)
      ilascl = .FALSE.
      IF ( anrm>ZERO .AND. anrm<smlnum ) THEN
         anrmto = smlnum
         ilascl = .TRUE.
      ELSEIF ( anrm>bignum ) THEN
         anrmto = bignum
         ilascl = .TRUE.
      ENDIF
      IF ( ilascl ) CALL ZLASCL('G',0,0,anrm,anrmto,N,N,A,Lda,ierr)
!
!     Scale B if max element outside range [SMLNUM,BIGNUM]
!
      bnrm = ZLANGE('M',N,N,B,Ldb,Rwork)
      ilbscl = .FALSE.
      IF ( bnrm>ZERO .AND. bnrm<smlnum ) THEN
         bnrmto = smlnum
         ilbscl = .TRUE.
      ELSEIF ( bnrm>bignum ) THEN
         bnrmto = bignum
         ilbscl = .TRUE.
      ENDIF
      IF ( ilbscl ) CALL ZLASCL('G',0,0,bnrm,bnrmto,N,N,B,Ldb,ierr)
!
!     Permute the matrices A, B to isolate eigenvalues if possible
!
      ileft = 1
      iright = N + 1
      irwrk = iright + N
      CALL ZGGBAL('P',N,A,Lda,B,Ldb,ilo,ihi,Rwork(ileft),Rwork(iright), &
     &            Rwork(irwrk),ierr)
!
!     Reduce B to triangular form (QR decomposition of B)
!
      irows = ihi + 1 - ilo
      IF ( ilv ) THEN
         icols = N + 1 - ilo
      ELSE
         icols = irows
      ENDIF
      itau = 1
      iwrk = itau + irows
      CALL ZGEQRF(irows,icols,B(ilo,ilo),Ldb,Work(itau),Work(iwrk),     &
     &            Lwork+1-iwrk,ierr)
!
!     Apply the orthogonal transformation to matrix A
!
      CALL ZUNMQR('L','C',irows,icols,irows,B(ilo,ilo),Ldb,Work(itau),  &
     &            A(ilo,ilo),Lda,Work(iwrk),Lwork+1-iwrk,ierr)
!
!     Initialize VL
!
      IF ( ilvl ) THEN
         CALL ZLASET('Full',N,N,CZERO,CONE,Vl,Ldvl)
         IF ( irows>1 ) CALL ZLACPY('L',irows-1,irows-1,B(ilo+1,ilo),   &
     &                              Ldb,Vl(ilo+1,ilo),Ldvl)
         CALL ZUNGQR(irows,irows,irows,Vl(ilo,ilo),Ldvl,Work(itau),     &
     &               Work(iwrk),Lwork+1-iwrk,ierr)
      ENDIF
!
!     Initialize VR
!
      IF ( ilvr ) CALL ZLASET('Full',N,N,CZERO,CONE,Vr,Ldvr)
!
!     Reduce to generalized Hessenberg form
!
      IF ( ilv ) THEN
!
!        Eigenvectors requested -- work on whole matrix.
!
         CALL ZGGHD3(Jobvl,Jobvr,N,ilo,ihi,A,Lda,B,Ldb,Vl,Ldvl,Vr,Ldvr, &
     &               Work(iwrk),Lwork+1-iwrk,ierr)
      ELSE
         CALL ZGGHD3('N','N',irows,1,irows,A(ilo,ilo),Lda,B(ilo,ilo),   &
     &               Ldb,Vl,Ldvl,Vr,Ldvr,Work(iwrk),Lwork+1-iwrk,ierr)
      ENDIF
!
!     Perform QZ algorithm (Compute eigenvalues, and optionally, the
!     Schur form and Schur vectors)
!
      iwrk = itau
      IF ( ilv ) THEN
         chtemp = 'S'
      ELSE
         chtemp = 'E'
      ENDIF
      CALL ZHGEQZ(chtemp,Jobvl,Jobvr,N,ilo,ihi,A,Lda,B,Ldb,Alpha,Beta,  &
     &            Vl,Ldvl,Vr,Ldvr,Work(iwrk),Lwork+1-iwrk,Rwork(irwrk), &
     &            ierr)
      IF ( ierr/=0 ) THEN
         IF ( ierr>0 .AND. ierr<=N ) THEN
            Info = ierr
         ELSEIF ( ierr>N .AND. ierr<=2*N ) THEN
            Info = ierr - N
         ELSE
            Info = N + 1
         ENDIF
         GOTO 100
      ENDIF
!
!     Compute Eigenvectors
!
      IF ( ilv ) THEN
         IF ( .NOT.(ilvl) ) THEN
            chtemp = 'R'
         ELSEIF ( ilvr ) THEN
            chtemp = 'B'
         ELSE
            chtemp = 'L'
         ENDIF
!
         CALL ZTGEVC(chtemp,'B',ldumma,N,A,Lda,B,Ldb,Vl,Ldvl,Vr,Ldvr,N, &
     &               in,Work(iwrk),Rwork(irwrk),ierr)
         IF ( ierr/=0 ) THEN
            Info = N + 2
            GOTO 100
         ENDIF
!
!        Undo balancing on VL and VR and normalization
!
         IF ( ilvl ) THEN
            CALL ZGGBAK('P','L',N,ilo,ihi,Rwork(ileft),Rwork(iright),N, &
     &                  Vl,Ldvl,ierr)
            DO jc = 1 , N
               temp = ZERO
               DO jr = 1 , N
                  temp = MAX(temp,ABS1(Vl(jr,jc)))
               ENDDO
               IF ( temp>=smlnum ) THEN
                  temp = ONE/temp
                  DO jr = 1 , N
                     Vl(jr,jc) = Vl(jr,jc)*temp
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
         IF ( ilvr ) THEN
            CALL ZGGBAK('P','R',N,ilo,ihi,Rwork(ileft),Rwork(iright),N, &
     &                  Vr,Ldvr,ierr)
            DO jc = 1 , N
               temp = ZERO
               DO jr = 1 , N
                  temp = MAX(temp,ABS1(Vr(jr,jc)))
               ENDDO
               IF ( temp>=smlnum ) THEN
                  temp = ONE/temp
                  DO jr = 1 , N
                     Vr(jr,jc) = Vr(jr,jc)*temp
                  ENDDO
               ENDIF
            ENDDO
         ENDIF
      ENDIF
!
!     Undo scaling if necessary
!
!
 100  IF ( ilascl ) CALL ZLASCL('G',0,0,anrmto,anrm,N,1,Alpha,N,ierr)
!
      IF ( ilbscl ) CALL ZLASCL('G',0,0,bnrmto,bnrm,N,1,Beta,N,ierr)
!
      Work(1) = DCMPLX(lwkopt)
!
!     End of ZGGEV3
!
      END SUBROUTINE ZGGEV3


text:
      subroutine sfvbak(nsfmax,nnsfmx,nsf,nsfnum,nsfnod,
     1                  sfv,sfvnum,sfvnod)
c
      implicit none
      integer  i,j
      integer  nsfmax,nnsfmx,nsf
      integer  nsfnum(nsfmax),nsfnod(nsfmax,nnsfmx)
      integer  sfv(2),sfvnum(2,nsfmax),sfvnod(2,nsfmax,nnsfmx)
c
      nsf = sfv(1)
      do i=1,nsf
         nsfnum(i) = sfvnum(1,i)
         do j=1,nsfnum(i)
            nsfnod(i,j) = sfvnod(1,i,j)
         end do
      end do
c
      return
      end


text:
program factorial_function
    use, intrinsic :: iso_fortran_env, only : i8 => INT64
    implicit none
    integer(kind=i8), parameter :: n_min = 0_i8, n_max = 10_i8
    integer(kind=i8) :: n

    do n = n_min, n_max
        print *, n, factorial(n)
    end do

contains

    integer(kind=i8) function factorial(n)
        use, intrinsic :: iso_fortran_env, only : i8 => INT64
        implicit none
        integer(kind=i8), intent(in) :: n
        integer(kind=i8) :: factor
        factorial  = 1_i8
        do factor = 2_i8, n
            factorial = factorial*factor
        end do
    end function factorial

end program factorial_function

max_stars_repo_name: gjbex/FortranForProgrammers

max_issues_repo_name: gjbex/FortranForProgrammers

max_forks_repo_name: gjbex/FortranForProgrammers


module compute_n_steady_module

  use ml_layout_module
  use multifab_physbc_module
  use define_bc_module
  use bc_module
  use bndry_reg_module
  use ml_solve_module
  use probin_reactdiff_module, only: nspecies, mg_verbose, cg_verbose, &
       implicit_diffusion_rel_eps, implicit_diffusion_abs_eps, D_Fick

  implicit none

  private

  public :: compute_n_steady

contains

  subroutine compute_n_steady(mla,n_steady,dx,dt,the_bc_tower)

    ! compute div D_k grad n_steady_k = 0 with inhomogeneous boundary conditions

    type(ml_layout), intent(in   ) :: mla
    type(multifab) , intent(inout) :: n_steady(:)
    real(kind=dp_t), intent(in   ) :: dx(:,:),dt
    type(bc_tower) , intent(in   ) :: the_bc_tower

    ! for multigrid solver; (alpha - div beta grad) phi = rhs
    type(multifab) :: alpha(mla%nlevel)
    type(multifab) :: rhs(mla%nlevel)
    type(multifab) :: phi(mla%nlevel)
    type(multifab) :: beta(mla%nlevel,mla%dim)

    ! for diffusion multigrid - not used but needs to be passed in
    type(bndry_reg) :: fine_flx(mla%nlevel)

    integer :: i,dm,n,nlevs,spec

    type(bl_prof_timer),save :: bpt

    nlevs = mla%nlevel
    dm = mla%dim

    ! if we are periodic in all directions, set n_steady=0 and return
    if (all(mla%pmask(1:dm))) then
       do n=1,nlevs
          call multifab_setval(n_steady(n),0.d0,all=.true.)
       end do
       return
    end if

    call build(bpt,"compute_n_steady")

    ! for multigrid solver; (alpha - div beta grad) phi = rhs
    do n=1,nlevs
       call multifab_build(alpha(n),mla%la(n),1,0)
       call multifab_build(rhs(n)  ,mla%la(n),1,0)
       call multifab_build(phi(n)  ,mla%la(n),1,1)
       do i=1,dm
          call multifab_build_edge(beta(n,i),mla%la(n),1,0,i)
       end do
    end do

    ! stores beta*grad phi/dx_fine on coarse-fine interfaces
    ! this gets computed inside of ml_cc_solve
    ! we pass it back out because some algorithms (like projection methods)
    ! use this information
    do n = 1,nlevs
       call bndry_reg_build(fine_flx(n),mla%la(n),ml_layout_get_pd(mla,n))
    end do

    do n=1,nlevs
       call multifab_setval(alpha(n),0.d0,all=.true.)
       call multifab_setval(rhs(n),0.d0,all=.true.)
    end do

    do spec=1,nspecies

       ! beta = D_k
       do n=1,nlevs
          do i=1,dm
             call multifab_setval(beta(n,i),D_Fick(spec),all=.true.)
          end do
       end do

       ! initial guess for phi
       ! zero on the interior, ghost cells filled with Dirichlet value
       do n=1,nlevs
          call multifab_setval(phi(n),0.d0,all=.true.)
          call multifab_physbc(phi(n),1,scal_bc_comp+spec-1,1, &
                               the_bc_tower%bc_tower_array(n),dx_in=dx(n,:))
       end do

       ! solve the Poisson equation
       call ml_cc_solve(mla,rhs,phi,fine_flx,alpha,beta,dx, &
                        the_bc_tower,scal_bc_comp+spec-1, &
                        stencil_order=1, &
                        verbose=mg_verbose, &
                        cg_verbose=cg_verbose, &
                        eps=implicit_diffusion_rel_eps, &
                        abs_eps=implicit_diffusion_abs_eps)

       ! copy solution into n_new
       do n=1,nlevs
          call multifab_copy_c(n_steady(n),spec,phi(n),1,1,0)
       end do

    end do

    do n=1,nlevs
       call multifab_destroy(alpha(n))
       call multifab_destroy(rhs(n))
       call multifab_destroy(phi(n))
       do i=1,dm
          call multifab_destroy(beta(n,i))
       end do
    end do

    do n = 1,nlevs
       call bndry_reg_destroy(fine_flx(n))
    end do

    call destroy(bpt)

  end subroutine compute_n_steady

end module compute_n_steady_module


subroutine sub_read
	use mod_param
	implicit none

   !Guess for the value function:
    open(15, file='V0.txt', status='old')
	read(15,*) V0_0
	close(15)
	V0 = reshape(V0_0,[nagrid, negrid, ntgrid])

   !Used for decompositions:
    if (fix_occpol == 1) then
        open(15, file='occpol_fixed_note.txt', status='old')
	    read(15,*) occpol_fixed_0
	    close(15)
	    occpol_fixed = reshape(occpol_fixed_0,[nagrid, negrid, ntgrid])
    endif

    if (fix_kpol == 1) then

        open(15, file='policycap_fixed_note.txt', status='old')
	    read(15,*) policycap_fixed_0
	    close(15)
	    policycap_fixed = reshape(policycap_fixed_0,[nagrid, negrid, ntgrid])

        open(15, file='Icap_fixed_note.txt', status='old')
	    read(15,*) Icap_fixed_0
	    close(15)
	    Icap_fixed = reshape(Icap_fixed_0,[nagrid, negrid, ntgrid])

    endif

    if (fix_apol_nd == 1) then

        open(15, file='apolse0_fixed_note.txt', status='old')
	    read(15,*) apolse0_fixed_note_0
	    close(15)
	    apolse0_fixed_note = reshape(apolse0_fixed_note_0,[nagrid, negrid, ntgrid])

    endif

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

   	!Reading in the 3 grids for assets:
    open(15, file='agrid.txt', status='old')
	read(15,*) agrid
	close(15)

    open(15, file='agrid_fine.txt', status='old')
	read(15,*) agrid_fine
	close(15)

    open(15, file='agrid_dist.txt', status='old')
	read(15,*) agrid_dist
	close(15)

    !Reading in the grid for epsilon:
    open(15, file='eps_grid.txt', status='old')
	read(15,*) eps_grid
	close(15)

    !Reading in the grid for theta:
    open(15, file='theta.txt', status='old')
	read(15,*) theta
	close(15)

    !Reading in the grid for phi (fraction of income taxes evaded):
    open(15, file='phigrid.txt', status='old')
	read(15,*) phigrid
	close(15)

    !Reading in the grid for kapital (business capital):
    open(15, file='kgrid.txt', status='old')
	read(15,*) kgrid
	close(15)

	!Reading in transition matrix for epsilon:
	open(15, file='Peps.txt', status='old')
	read(15,*) Peps0
	close(15)
    Peps = reshape(Peps0,[negrid, negrid])

	!Reading in the transition matrix for theta:
	open(15, file='Ptheta.txt', status='old')
	read(15,*) Ptheta0
	close(15)
    Ptheta = reshape(Ptheta0,[ntgrid, ntgrid])

    !Reading in the grid for labor supply
	open(15, file='lgrid.txt', status='old')
	read(15,*) lgrid
	close(15)

    !Reading in the grid for labor hirings:
	open(15, file='ngrid.txt', status='old')
	read(15,*) ngrid
	close(15)

    !Reading in the grid for labor supply entre:
	open(15, file='legrid.txt', status='old')
	read(15,*) legrid
	close(15)

end subroutine sub_read


text:
C This file is part of PREPRO.
C
C    Author: Dermott (Red) Cullen
C Copyright: (C) International Atomic Energy Agency
C
C PREPRO is free software; you can redistribute it and/or modify it
C under the terms of the MIT License; see LICENSE file for more details.

c=======================================================================
c
c     Program ENDF2C
c     ==============
c     Convert ENDF Data to Standard FORTRAN, C and C++ Format.
c
c     Version 2014-1 Feb. 2014 * Initial version.
C             2014-2 Oct. 2014 * Changed from D to E exponential form
C                                to improve compatibility between
C                                computer languages.
C             2015-1 Jan. 2015 * General updates for release with
C                                PREPRO2015.
C                              * Changed ENDF data filenames from ENDF2C
C                                to ENDFB, to agree with PREPRO default
C                                definitions.
C                              * Added code name (to be compatible
C                                with PREPRO output), but NOT TIME (to
C                                keep this code as computer independent
C                                as possible).
C             2017-1 May  2017 * Updated based on user feedbsck
C             2018-1 Jan. 2018 * Added on-line output for ALL ENDERROR
C             2019-1 June 2019 * Added /UNITS/ to allow correct output
C                                at end = output either o.k. or error.
C             2020-1 Feb. 2020 * Identical to 2019-1.
C             2021-1 Jan. 2021 * Updated for FOTRAN 2018
C
c     Purpose
c     ==================================================================
c     This code is designed for,
c     1) ENDF Data in any ENDF format = ENDF-1 through ENDF-6.
c     2) On any type of computer = 32 or 64 bit system/compiler
c
c     This code tries to keep things as simple as possible
c     1) There are NO INPUT PARAMETERS.
c     2) It reads  an ENDF formatted file named ENDFB.IN
c     3) It writes an ENDF formatted file named ENDFB.OUT
c     4) It writes a report file named          ENDF2C.LST
c
c     Author's Message
c     ----------------
c     I consider insuring that ENDF data is in a standard, officially
c     approved format for FORTRAN, C and C++ is SO IMPORTANT this code
c     does only one thing - and only one thing - and it does it in the
c     simplest possible manner - efficiency is NOT a consideration -
c     ONLY accuracy and general utility of the ENDF data is considered.
c
c     Method
c     ------
c     Other codes that attempt to do the same thing - including codes
c     written be me decades ago - are very complicated, and therefore
c     ERROR PRONE because they try to deal with each and every variant
c     in which data can be coded in the ENDF format. Needless to say
c     this means that every time the ENDF formats and procedures change
c     these codes MUSE also be changed.
c
c     In contrast, ENDF2C uses my almost 50 years of experience dealing
c     with the ENDF format to realize that except for the comments at
c     the beginning for each evaluation (MF/MT=1/451), every line of
c     ENDF data is IDENTICAL - in every version of the ENDF format, from
c     the original ENDF to today's ENDF-6. So to translate ENDF data
c     into an official format I do not have to consider differences in
c     each section (MF/MT) of data.
c
c     Every line of ENDF is divided into 6 fields, each 11 columns wide.
c     Each of the 6 fields is either, blank, integer or floating point.
c     Floating point fields ALL include a decimal point (.). So that ALL
c     this code does is convert every floating point field to standard
c     format.
c
c     In order to insure that this PRESERVES the accuracy of the data
c     this is done by reading and writing each ENDF line as characters.
c     Blank and integer fields are copied exactly as read. ALL floating
c     point number that are read are converted internally from character
c     to floating point - they are then converted back into characters
c     in a standard, officially approved format, for output.
c
c     As a last step to insure the accuracy of results the characters
c     to be output are again converted from characters to floating
c     point, and the numerical value that is output is compared to the
c     numerical value originally read, and if there is ANY DIFFERENCE
c     the characters strings read and written are listed in the output:
c     the characters strings read and written as well as the difference
c     is listed in the output report (ENDF2C.LST) and on the screen.
c
c     Running Time
c     ------------
c     It takes only seconds to translate an ENDF formatted evaluation,
c     so running time need not be a consideration. Concentrate on
c     keeping it simple and reliable - that should be your focus.
c
c     Documentation
c     -------------
c     ALL of my codes that process ENDF data and change it in ANY WAY
c     document what they have done by adding comment lines at the end
c     of the comment section (MF/MT=1/451) of each evaluation. This
c     allows data users to determine the pedigree of the data they are
c     using, by reading these comments. This code documents what is has
c     done by adding the following 2 comment lines.
c
c ***************** Program ENDF2C (Version 2021-1) ***************
c      Convert ENDF Data to Standard FORTRAN, C and C++ Format
c
c     WARNING - This documentation is IMPORTANT to data users and it
c     should not be deleted.
c
c     Written by
c     ------------------------------------
c     Dermott E. Cullen
c     University of California (retired)
c-----Present Home Address----------------------------------------------
c     Dermott E. Cullen
c     1466 Hudson Way
c     Livermore, CA 94550
c     U.S.A.
c     Telephone  925-443-1911
c     E. Mail    RedCullen1@Comcast.net
c     Website    RedCullen1.net/HOMEPAGE.NEW
c
c=======================================================================
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save
      character*8 codename
      character*1 LINE66
      integer*4 OUTP,OTAPE
      common/endf1/LINE66(11,6)
      common/endf2/MAT,MF,MT,NOSEQ,IVERSE
      common/ENDFIO/INP,OUTP,ITAPE,OTAPE
      common/namecode/codename
c-----------------------------------------------------------------------
c
c     Define Code name for output report and initialize time
c
c-----------------------------------------------------------------------
      codename = 'ENDF2C  '
c-----2019/1/8 - not used unit numbers, defined for PREPRO compatible
      INP   = 2
      OUTP  = 3
      ITAPE = 10
      OTAPE = 11
c-----Initialize timer.
      call TIMER
c-----------------------------------------------------------------------
c
c     Define I/O Files
c
c-----------------------------------------------------------------------
      open( 3,file='ENDF2C.LST')
      open(10,file='ENDFB.IN')
      open(11,file='ENDFB.OUT')
c-----------------------------------------------------------------------
c
c     Identify Code.
c
c-----------------------------------------------------------------------
      write(*,10)
      write(3,10)
   10 format(//
     1 ' ',13('*'),'**** Program ENDF2C (Version 2021-1) ***',
     2 12('*')/5(' '),
     3 ' Convert ENDF Data to Standard FORTRAN, C and C++ Format',
     4 5(' ')/1x,71('-')/' ENDF Tape Label (TPID)'/1x,71('-'))
c-----Initialize previous MAT read and ENDF format version #
      MATLAST = -99999
      IVERSE  = 0
      NOSEQ   = 0
c-----------------------------------------------------------------------
c
c     copy TPID - and identify it in output listing
c
c-----------------------------------------------------------------------
      read (10,20) LINE66,MAT,MF,MT
      write( *,20) LINE66,MAT,MF,MT,NOSEQ
      write( 3,20) LINE66,MAT,MF,MT,NOSEQ
   20 format(66A1,i4,i2,i3,i5)
c-----------------------------------------------------------------------
c
c     Check for expected TPID format = file may not have TPID.
c
c-----------------------------------------------------------------------
      if(MF.eq.0.and.MT.eq.0) then
      NOSEQ = 0
      write(11,20) LINE66,MAT,MF,MT,NOSEQ               ! Assume TPID
      write( *,30)
      write( 3,30)
   30 format(1x,71('-')/' Input       Output         Fractional Diff.',
     5 '    MAT MF  MT  Seq# Format'/
     6       1x,71('-'))
      else
      write( *,40)                                      ! NOT TPID
      write( 3,40)
   40 format(1x,71('-')/' WARNING - This DOES NOT look like a TPID.',
     1       ' Expect MF = MT = 0 on TPID.'/
     2       '           Will Assume file has NO TPID',
     3       ' and TRY to Translate this line.'/
     4 1x,71('-')/' Input       Output         Fractional Diff.',
     5 '    MAT MF  MT  Seq# Format'/
     6       1x,71('-'))
      go to 60                                            ! Translate
      endif
c-----------------------------------------------------------------------
c
c     Loop over ALL ENDF lines to TEND line
c
c-----------------------------------------------------------------------
   50 read (10,20,end=80) LINE66,MAT,MF,MT
c-----Initialize ENDF Format vesion # for each new MAT
   60 if(MAT.ne.MATLAST) then
      MATLAST = MAT
      IVERSE  = 0
      NOSEQ   = 0
      endif
c-----------------------------------------------------------------------
c
c     special treatment for SEND/FEND/MEND/TEND
c
c-----------------------------------------------------------------------
      if(MT.le.0) then               ! SEND/FEND/MEND/TEND?
      if(MF.gt.0) then
      NOSEQ = 99999                  ! SEND NOSEQ = 99999
      else
      NOSEQ = 0                      ! FEND/MEND/TEND NOSEQ = 0
      endif
      write(11,70) MAT,MF,MT,NOSEQ ! Output with blank fields
   70 format(66X,i4,i2,i3,i5)
      NOSEQ = 0                      ! Initialze before next line.
      if(MAT.lt.0) then
      go to 100                      ! TEND = End of ENDF data
      else
      go to 50                       ! Otherwise continue reading
      endif
      endif
c-----------------------------------------------------------------------
c
c     special treatment MF/MT=1/451 = Add Documentation
c
c-----------------------------------------------------------------------
      if(MF.eq.1.and.MT.eq.451) then
      call FILE1
      go to 50
      endif
c-----------------------------------------------------------------------
c
c     No other special cases - translate 6 data fields to standard form
c
c-----------------------------------------------------------------------
      NOSEQ = NXTSEQ(NOSEQ)
      call STANDARD(6)
c-----Output standard format
      write(11,20) LINE66,MAT,MF,MT,NOSEQ
      go to 50
c-----------------------------------------------------------------------
c
c     End of ENDF input
c
c-----------------------------------------------------------------------
c-----End of File BEFORE TEND
   80 write( *,90)
      write( 3,90)
   90 format(1x,71('-')/
     1 ' ERROR - ENDF file does not end with TEND Line.'/1x,71('-'))
      call ENDERROR
c-----TEND line read
  100 call endit
      end
      Subroutine FILE1
C=======================================================================
c
c     Add documentation to MF/MT=1/451
c
c     ENDF-4 Format: 2-nd line N1 > 0, N2 = 0
c         -5       : 3-rd line N1 = 0, N2 = 0
c         -6       : 4-th line N1 > 0, N2 > 0
c
C=======================================================================
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save
      character*1 LINE66,LINE44(44),LINE11(11),MTDICT(44)
      common/endf1/LINE66(11,6)
      common/endf2/MAT,MF,MT,NOSEQ,IVERSE
      equivalence (LINE66(1,1),LINE44(1)),(LINE66(1,6),LINE11)
c-----Initialize ENDF format version
      IVERSE = 0
c-----------------------------------------------------------------------
c
c     first line already read - translate C1 and C2 and output
c
c-----------------------------------------------------------------------
      NOSEQ = NXTSEQ(NOSEQ)
      call STANDARD(2)
      write(11,10) LINE66,MAT,MF,MT,NOSEQ
c-----------------------------------------------------------------------
c
c     Loop over UP TO next 3 lines = fewer for earlier ENDF versions.
c
c-----------------------------------------------------------------------
      do 60 iline=1,3
      read (10,10) LINE66,MAT,MF,MT
   10 format(66A1,i4,i2,i3,i5)
c-----Convert C1 and C2 fields
      NOSEQ = NXTSEQ(NOSEQ)
      call STANDARD(2)
c-----Define ENDF format version based on N1 and N2
      call in9(XN1,LINE66(1,5))
      call in9(XN2,LINE66(1,6))
      N1IN = XN1
      N2IN = XN2
c-----Based on which card read
      go to (20,30,40),iline
   20 N1A = N1IN
      N2A = N2IN
      if(N1A.le.0.or.N2A.ne.0) go to 50
      IVERSE = 4                              ! IV Format
      NWD    = N1A
      go to 70
   30 N1B = N1IN
      if(N2A.gt.0) go to 50
      IVERSE = 5                              ! V  Format
      NWD    = N1B
      go to 70
   40 N1C = N1IN
      IVERSE = 6                              ! VI Format
      NWD    = N1C
      go to 70
   50 write(11,10) LINE66,MAT,MF,MT,NOSEQ   ! Copy line = no change.
   60 CONTINUE
c-----Output updated comment count
   70 NWDOUT = NWD + 2
      write(11,80) LINE44,NWDOUT,LINE11,MAT,MF,MT,NOSEQ ! New count
   80 format(44A1,i11,11A1,i4,i2,i3,i5)
c-----------------------------------------------------------------------
c
c     Copy comments = insure first 22 coilumns are blank.
c
c-----------------------------------------------------------------------
      do iwd = 1,NWD
      read (10,10) LINE66,MAT,MF,MT
      NOSEQ = NXTSEQ(NOSEQ)
      write(11,10) LINE66,MAT,MF,MT,NOSEQ
      enddo
c-----------------------------------------------------------------------
c
c     Add comments + new sequence numbers
c
c-----------------------------------------------------------------------
      NOSEQ = NXTSEQ(NOSEQ)
      write(11,90) MAT,MF,MT,NOSEQ
      NOSEQ = NXTSEQ(NOSEQ)
      write(11,100) MAT,MF,MT,NOSEQ
   90 format(
     1 ' ',13('*'),'**** Program ENDF2C (Version 2021-1) ***',
     2 12('*'),i4,i2,i3,i5)
  100 format(5(' '),
     1 'Convert ENDF Data to Standard FORTRAN, C and C++ Format',
     2 6(' '),i4,i2,i3,i5)
c-----------------------------------------------------------------------
c
c     Copy to SEND
c
c-----------------------------------------------------------------------
  110 read (10,120) MTDICT,MAT,MF,MT
  120 format(22x,44A1,i4,i2,i3,i5)
      if(MT.gt.0) then
      NOSEQ = NXTSEQ(NOSEQ)
      write(11,120) MTDICT,MAT,MF,MT,NOSEQ
      go to 110
      endif
      NOSEQ = 99999                     ! SEND line
      write(11,130) MAT,MF,MT,NOSEQ
  130 format(66x,i4,i2,i3,i5)
      NOSEQ = 0                         ! Initialize for next line
      return
      end
      Subroutine STANDARD(NFIELD)
C=======================================================================
c
c     Scan ALL other lines
c
c     NFIELD = # of fields to scan = 2 or 6
c
C=======================================================================
      implicit real*8 (a-h,o-z)
      implicit integer*4 (i-n)
      save
      character*1 LINE66
      character*1 LINESAVE(11)
      common/endf1/LINE66(11,6)
      common/endf2/MAT,MF,MT,NOSEQ,IVERSE
c-----------------------------------------------------------------------
c
c     Scan fields for . No other test is necessary
c
c-----------------------------------------------------------------------
      do 30 i = 1,NFIELD
      do k = 1,11
      if(LINE66(k,i).eq.'.') go to 10
      enddo
c-----no . = nothing to change
      go to 30
c-----save original form
   10 do kk=1,11
      LINESAVE(kk) = LINE66(kk,i)
      enddo
c-----translate
      call IN9 (zin,LINE66(1,i))
c-----Try for 10 digit output - switch to 9 if negative.
      call OUT10(zin,LINE66(1,i))
c-----------------------------------------------------------------------
c
c     Compare original and final = Print if significantly different
c
c-----------------------------------------------------------------------
      call IN9 (zix,LINE66(1,i))
      if(DABS(zin).le.0.0d+0) then   ! Avoid divding by zero.
      diff = (zin - zix)
      else
      diff = (zin - zix)/zin
      endif
      if(dabs(diff).gt.1.0d-10) then
      write(*,20) LINESAVE,(LINE66(kk,i),kk=1,11),diff,
     1 MAT,MF,MT,NOSEQ,IVERSE
      write(3,20) LINESAVE,(LINE66(kk,i),kk=1,11),diff,
     1 MAT,MF,MT,NOSEQ,IVERSE
   20 format(1x,11a1,1x,11a1,2x,1pd20.12,i5,i3,i4,i6,i5)
      endif
c-----end of field loop
   30 continue
      return
      end

max_stars_repo_name: IAEA-NDS/PREPRO

max_issues_repo_name: IAEA-NDS/PREPRO

max_forks_repo_name: IAEA-NDS/PREPRO


text:
	SUBROUTINE GDU_MTCH ( nexp, gdattm, level, ivcord, parm,
     +			      levtyp, nlev, rlevel, icord, nparm,
     +			      prmlst, ntime, timfnd, match, iret )
C************************************************************************
C* GDU_MTCH								*
C*									*
C* This routine checks for a match between the user input grid info	*
C* and the grid header read from a file.				*
C*									*
C* GDU_MTCH ( NEXP, GDATTM, LEVEL, IVCORD, PARM, LEVTYP, NLEV, RLEVEL,  *
C*            ICORD, NPARM, PRMLST, NTIME, TIMFND, MATCH, IRET )	*
C*									*
C* Input parameters:							*
C*	NEXP		INTEGER		Max number of levels		*
C*	GDATTM (2)	CHAR*		Date/times from the grid header	*
C*	LEVEL (2)	INTEGER		Levels from the grid header	*
C*	IVCORD		INTEGER		Vert coord from the grid header	*
C*	PARM		CHAR*		Parameter from the grid header	*
C*	LEVTYP		INTEGER		Level type from the user	*
C*					  0 = no levels input		*
C*					  1 = list of levels		*
C*					  2 = range of levels		*
C*					  3 = all levels		*
C*	NLEV		INTEGER		Number of levels		*
C*	RLEVEL (NEXP,2)	REAL		Levels or range from the user	*
C*	ICORD		INTEGER		Vert coord from the user	*
C*	NPARM		INTEGER		Number of parameters		*
C*	PRMLST (NPARM)	CHAR*		Parameter list from the user	*
C*	NTIME		INTEGER		Number of times			*
C*	TIMFND (NTIME)	CHAR*		List of times from the user	*
C*									*
C* Output parameters:							*
C*	MATCH		LOGICAL		Match flag			*
C*	IRET		INTEGER		Return code			*
C*									*
C**									*
C* Log: 								*
C* S. Jacobs/NCEP	9/96						*
C* S. Maxwell/GSC	9/96		Clean up.			*
C************************************************************************
	CHARACTER*(*)	gdattm(2), parm, prmlst(*), timfnd(*)
	INTEGER		level(2)
	REAL		rlevel( NEXP, 2 )
	LOGICAL		match
C*
	CHARACTER	tgdtm*36
C------------------------------------------------------------------------
	iret  = 0
	match = .false.
C
C*	If the user entered ALL, then match is true.
C
	IF  ( timfnd(1) .eq. 'ALL' )  THEN
	    match = .true.
	  ELSE
C
C*	    Concatenate the times.
C
	    CALL ST_LSTR ( gdattm(1), len1, ier )
	    CALL ST_LSTR ( gdattm(2), len2, ier )
	    IF ( len1 .gt. 0 ) THEN
		IF ( len2 .gt. 0 ) THEN
		    tgdtm = gdattm(1)(:len1) // ':' // gdattm(2)(:len2)
		  ELSE
		    tgdtm = gdattm(1)
		END IF
	    END IF
C
C*	    Check all the user input times for a match.
C
	    DO  k = 1, ntime
		IF  ( tgdtm .eq. timfnd(k) )  match = .true.
	    END DO
	END IF
C
C*	Check the vertical coordinate.
C
	IF  ( match )  THEN
	    match = .false.
C
C*	    If the user entered ALL, then match is true.
C
	    IF  ( icord .eq. -2 )  THEN
		match = .true.
	      ELSE
C
C*		Otherwise, check for a match.
C
		IF  ( ivcord .eq. icord )  match = .true.
	    END IF
	END IF
C
C*	Check the levels.
C
	IF  ( match )  THEN
	    match = .false.
C
C*	    If the user entered ALL, then match is true.
C
	    IF  ( levtyp .eq. 3 )  THEN
		match = .true.
	      ELSE IF  ( levtyp .eq. 2 )  THEN
C
C*		Check the range of levels.
C
		IF  ( ( ( ( rlevel(1,1) .le. level(1) ) .and.
     +		          ( rlevel(2,1) .ge. level(1) ) ) .or.
     +		        ( ( rlevel(1,1) .ge. level(1) ) .and.
     +		          ( rlevel(2,1) .le. level(1) ) ) ) .and.
     +		      ( level(2) .eq. -1 ) )  match = .true.

	      ELSE
C
C*		Otherwise, check all the individual levels.
C
		DO  k = 1, nlev
		    IF  ( ( level(1) .eq. rlevel(k,1) ) .and.
     +			  ( level(2) .eq. rlevel(k,2) ) )
     +			match = .true.
		END DO
	    END IF
	END IF
C
C*	Check the list of parameters.
C
	IF  ( match )  THEN
	    match = .false.
C
C*	    If the user entered ALL, then match is true.
C
	    IF  ( prmlst(1) .eq. 'ALL' )  THEN
		match = .true.
	      ELSE
C
C*		Otherwise, check list of parameters for a match.
C
		DO  k = 1, nparm
		    IF  ( parm .eq. prmlst(k) )  match = .true.
		END DO
	    END IF
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
	options/extend_source
!------------------------------------------------------------------------------

	PROGRAM MISSION_SCAN
C
C	include		'wind_examples:wind_tm_routine_def.for'
C	include		'wind_examples:wind_return_code_def.for'
C
        COMMON /MF37OUT/ ACOUNT,XASV,YASV,ZASV,XDSV,YDSV,ZDSV,TESV,
     1          DENSV,ETASV,SWEPOTSV
	INTEGER*4 FHIST
	character*80	stream,FILENAME
	character*120	Junk
c	CHARACTER*80	YYYYMMDD
	real*8 		scet8,scet,scetstr,scetend
	INTEGER*4 	OK,GET_STREAM_NAME,TMCHANNEL
	INTEGER*4	W_CHANNEL_OPEN,W_CHANNEL_CLOSE
	integer*4	yyyy,mon,dd,hh,mm,ss,dds
        real*4          acount(5,10,10),xasv(10,10),yasv(10,10),zasv(10,10)
        real*4          xdsv(10,10),ydsv(10,10),zdsv(10,10),TESV(10,10)
        REAL*4          DENSV(10,10),SWEPOTSV(10,10),ETASV(10,10)
C
C	DATA SCET8/4697.5D00/		! NOV-11-1994 = TURN ON DAY
C	DATA SCET8/5313.5D00/		! JUL-19-1996 = TEST
C	DATA SCET8/5479.0D00/		! JAN-01-1997 = TEST
C	DATA SCET8/5708.0D00/		! AUG-18-1997 = 3B TEST
C	DATA SCET8/5245.5D00/		! MAY-12-1996 = START OF LONG ORB
C	DATA SCET8/5233.5D00/		! MAY-1-1996
C	DATA SCET8/5195.D00/		! MAR-24-1996
C	DATA SCET8/5217.D00/		! APR-15-1996
c	DATA SCET8/5256.D00/		! MAY-23-1996 =   " LOW BIT RATE
c	DATA SCET8/5260.D00/		! MAY-27-1996 =    test
C	DATA SCET8/4741.0D00/		! JAN-01-1995 = START OF LONG ORB
C
C	WRITE(88,1012)
C 1012	format(2X,'TIME OF CHANGE',4X,'DRV  RX  RY  RZ  VX  VY  VZ  DRV V' )
	READ(5,*) SCETSTR,SCETEND
	SCET8 = SCETSTR
	PRINT*,'DAYS',SCET8,SCETEND
C	OPEN(UNIT=56,FILE='MAKEFILE37.RESULTS',STATUS='OLD',ACCESS='APPEND')
C	OPEN(UNIT=57,FILE='MAKEFILE18.VOLTS',STATUS='OLD',ACCESS='APPEND')
 100	continue
	call w_ur8_to_ymd(scet8,yyyy,mon,dd,hh,mm,ss,ms)
	PRINT*,'YMD',YYYY,MON,DD,HH,MM
	WRITE(stream,30) YYYY,MON,DD
 30	FORMAT(I4.4,I2.2,I2.2)
	OK = GET_STREAM_NAME(stream)
C
C	call w_ur8_to_ydoy(scet8,yyyy,idoy,msec)
C
C	print*,'in main, stream=',stream
C
	ok = w_channel_open(tmchannel,stream)
C
	print*,'open, TMCHANNEL,ok=',TMCHANNEL,ok
	if (ok.ne.1) stop 'Cannot open t/m channel'
c
	scet = 0.
	call w_channel_position(tmchannel,scet)
	print*,'t/m file starts at scet',scet
	dds = scet
C
C	open SWE data file and advance through header
C
        WRITE(FILENAME,123) YYYY,MON,DD
C 123   FORMAT('MONTHLY:[KELLOGG.EFLUX]FLUX_',I4.4,I2.2,I2.2,'.ASC')
C 123    FORMAT('USER_A:[KELLOGG.WIND.PAPERS]FLUX_',I4.4,I2.2,I2.2,'.ASC')
 123    FORMAT('DATA_A:[KELLOGG.WIND]SWE_H2_',I4.4,I2.2,I2.2,'.DAT')
        print*,'recheck time',scet,yyyy,mon,dd
        print*,' name of swe file ',filename
c        OPEN(UNIT=85,FILE=FILENAME,STATUS='OLD',IOSTAT=I85ERR,READONLY)
c	PRINT*,'OPEN 85, IOSTAT= ',I85ERR
c	IF (I85ERR.NE.0) THEN
c	  SCET8 = SCET8 + 1.D00
c	  ERR = W_CHANNEL_CLOSE(TMCHANNEL)
c	  GO TO 100
c	ENDIF
c	DO N = 1,9
c	  READ(85,1083) JUNK
c	  PRINT*,JUNK
c	ENDDO
c 1083	FORMAT(A)
c	WRITE(86,*) SCET,YYYY,MON,DD
C
	call makefile(tmchannel)
C	CALL MAKEFILE
	print*,'closing unit 85'
	CLOSE(UNIT=83)
	CLOSE(UNIT=84)
	CLOSE(UNIT=85)
	print*,'closing tm channel, tmchannel=',tmchannel
	OK = W_CHANNEL_CLOSE(TMCHANNEL)
	scet8 = scet8 + 1.d00
c	if(scet8.lt.4720.d00) go to 100
C	if(scet8.lt.5330.d00) go to 100    ! APPROX END OF LONG ORB AUG 1996
C	if(scet8.lt.4964.d00) go to 100    ! APPROX END OF LONG ORB AUG 1995
C	if(scet8.lt.5511.d00) go to 100    ! FEB-01-1997
c	if(scet8.lt.5262.d00) go to 100    ! MAY 29-1997
C	if(scet8.lt.5253.d00) go to 100    ! MAY 20-1996
C	if(scet8.lt.5201.d00) go to 100    ! MAR 30-1996
C	if(scet8.lt.5223.d00) go to 100    ! MAR 30-1996
c	if(scet8.lt.5286.d00) go to 100    ! ONE MONTH LARGE-1997
C	if(scet8.lt.5708.5d00) go to 100    ! TEST
	if(scet8.lt.scetend) go to 100
C
C	WRITE(89,*) SCETSTR,SCETEND
C	write(89,1088) FHIST,NHISTST
C	print*,'fhist',fhist
 1088	FORMAT(14I5)
C
       type*,'write accumulated samples as function of angle'
C
 1087   FORMAT(10F7.2)

	print*, 'STOP IN MISSION_SCAN38'
	CLOSE(UNIT=86)
	STOP
	END
	options/extend_source
!------------------------------------------------------------------------------
	integer*4	function	get_stream_name(stream)
! This routine gets the user's TM stream type specification.
!
	implicit	none
	character*(*)	stream
	CHARACTER*80	YYYYMMDD
	integer*4	iq
	DATA IQ /1/
C
C10	  write(6,6)
C	  write(6,7)
C	  read(5,5,err=10,end=20) iq, stream
C	print*,'in get_, initial stream=  ',stream
c 	if (iq .lt. 1) then
c	   stream = 'offline'
c	else if (stream(1:1) .eq. 'o' .or. stream(1:1) .eq. 'O') then
c	   stream = 'offline'
c	else
c	   ! assume the user entered the TIME of an offline file
	   YYYYMMDD = STREAM(1:8)
	   WRITE(STREAM,30) YYYYMMDD
 30	   FORMAT('wi_lz_wav_',A8,'_v*.dat')
	   PRINT*,'in get_stream_name, file= ',STREAM
c	end if

	get_stream_name = 1

 20	return
c
  5	format(q,a)
  6	format(1x,'Enter TM stream type [O=offline (default), R=realtime ]: ',$)
  7	format(1x,'or type desired time as YYYYMMDD, e.g. 19961113  ',$)
c
	end