      program mpi_ver
      implicit none
      include 'mpif.h'
      integer(kind=kind(MPI_VERSION)), parameter :: zero = ichar('0')
      character, dimension(17), parameter :: mpiver_str =&
      (/ 'I', 'N', 'F', 'O', ':', 'M', 'P', 'I', '-', 'V', 'E', 'R', '[', &
        char(zero + MPI_VERSION), &
        '.', &
        char(zero + MPI_SUBVERSION), ']' /)
      print *, mpiver_str
      end program mpi_ver


text:
! -----------------------------------------------------------------
! Programmer(s): Cody J. Balos @ LLNL
! -----------------------------------------------------------------
! SUNDIALS Copyright Start
! Copyright (c) 2002-2019, Lawrence Livermore National Security
! and Southern Methodist University.
! All rights reserved.
!
! See the top-level LICENSE and NOTICE files for details.
!
! SPDX-License-Identifier: BSD-3-Clause
! SUNDIALS Copyright End
! -----------------------------------------------------------------
! This file contains a Fortran module for interfacing directly with
! the SUNDIALS SPGMR linear solver using the ISO_C_BINDING module.
! -----------------------------------------------------------------

module fsunlinsol_spfgmr_mod

  !======= Interfaces =========
  interface

     ! =================================================================
     ! Constructors
     ! =================================================================

     type(c_ptr) function FSUNLinSol_SPFGMR(y, pretype, maxl) &
          bind(C,name='SUNLinSol_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: y
       integer(c_int), value :: pretype
       integer(c_int), value :: maxl
     end function FSUNLinSol_SPFGMR

     ! =================================================================
     ! Destructors
     ! =================================================================

     subroutine FSUNLinSolFree_SPFGMR(LS) &
          bind(C,name='SUNLinSolFree_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
     end subroutine FSUNLinSolFree_SPFGMR

     ! =================================================================
     ! Setters
     ! =================================================================

     integer(c_int) function FSUNLinSol_SPFGMRSetPrecType(LS, pretype) &
          bind(C,name='SUNLinSol_SPFGMRSetPrecType')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: LS
       integer(c_int), value :: pretype
     end function FSUNLinSol_SPFGMRSetPrecType

     integer(c_int) function FSUNLinSol_SPFGMRSetGSType(LS, gstype) &
          bind(C,name='SUNLinSol_SPFGMRSetGSType')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: LS
       integer(c_int), value :: gstype
     end function FSUNLinSol_SPFGMRSetGSType

     integer(c_int) function FSUNLinSol_SPFGMRSetMaxRestarts(LS, maxrs) &
          bind(C,name='SUNLinSol_SPFGMRSetMaxRestarts')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: LS
       integer(c_int), value :: maxrs
     end function FSUNLinSol_SPFGMRSetMaxRestarts

     ! =================================================================
     ! Operations
     ! =================================================================

     integer(c_int) function FSUNLinSolGetType_SPFGMR(LS) &
          bind(C,name='SUNLinSolGetType_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
     end function FSUNLinSolGetType_SPFGMR

     integer(c_int) function FSUNLinSolInitialize_SPFGMR(LS) &
          bind(C,name='SUNLinSolInitialize_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
     end function FSUNLinSolInitialize_SPFGMR

     integer(c_int) function FSUNLinSolSetATimes_SPFGMR(LS, A_data, ATimes) &
          bind(C,name='SUNLinSolSetATimes_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: LS
       type(c_ptr),    value :: A_data
       type(c_funptr), value :: ATimes
     end function FSUNLinSolSetATimes_SPFGMR

     integer(c_int) function FSUNLinSolSetPreconditioner_SPFGMR(LS,     &
                                                               P_data, &
                                                               Pset,   &
                                                               Psol)   &
          bind(C,name='SUNLinSolSetPreconditioner_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: LS
       type(c_ptr),    value :: P_data
       type(c_funptr), value :: Pset
       type(c_funptr), value :: Psol
     end function FSUNLinSolSetPreconditioner_SPFGMR

     integer(c_int) function FSUNLinSolSetScalingVectors_SPFGMR(LS, s1, s2) &
          bind(C,name='SUNLinSolSetScalingVectors_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: LS
       type(c_ptr),    value :: s1
       type(c_ptr),    value :: s2
     end function FSUNLinSolSetScalingVectors_SPFGMR

     integer(c_int) function FSUNLinSolSetup_SPFGMR(LS, A) &
          bind(C,name='SUNLinSolSetup_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
       type(c_ptr), value :: A
     end function FSUNLinSolSetup_SPFGMR

     integer(c_int) function FSUNLinSolSolve_SPFGMR(LS, A, x, b, tol) &
          bind(C,name='SUNLinSolSolve_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr),    value :: LS
       type(c_ptr),    value :: A
       type(c_ptr),    value :: x
       type(c_ptr),    value :: b
       real(c_double), value :: tol
     end function FSUNLinSolSolve_SPFGMR

     integer(c_int) function FSUNLinSolNumIters_SPFGMR(LS) &
          bind(C,name='SUNLinSolNumIters_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
     end function FSUNLinSolNumIters_SPFGMR

     real(c_double) function FSUNLinSolResNorm_SPFGMR(LS) &
          bind(C,name='SUNLinSolResNorm_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
     end function FSUNLinSolResNorm_SPFGMR

     type(c_ptr) function FSUNLinSolResid_SPFGMR(LS) &
          bind(C,name='SUNLinSolResid_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
     end function FSUNLinSolResid_SPFGMR

     integer(c_long) function FSUNLinSolLastFlag_SPFGMR(LS) &
          bind(C,name='SUNLinSolLastFlag_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
     end function FSUNLinSolLastFlag_SPFGMR

     integer(c_int) function FSUNLinSolSpace_SPFGMR(LS, lenrwLS, leniwLS) &
          bind(C,name='SUNLinSolSpace_SPFGMR')
       use, intrinsic :: iso_c_binding
       implicit none
       type(c_ptr), value :: LS
       integer(c_long)    :: lenrwLS
       integer(c_long)    :: leniwLS
     end function FSUNLinSolSpace_SPFGMR

  end interface

end module fsunlinsol_spfgmr_mod

max_stars_repo_name: jrmie/math

max_issues_repo_name: jrmie/math

max_forks_repo_name: jrmie/math


text:
C     ( Last modified on 10 Aug 2005 at 15:54:33 )
      PROGRAM CLSSFY
C----------------------------------------------------------------------------
C
C  The purpose of this program is to add the classification characteristics
C  of the problem (expressed in SIF format and containing an explicit
C  classification to a database file (CLASSF.DB).  This classification must
C  occur in the file before the GROUPS or VARIABLES section of the SDIF,  and
*  must be of the form
C
C  * classification XXXn-XX-m-m
C
C  or
C
C  * CLASSIFICATION XXXn-XX-m-m
C
C  where n and m are integers containing from 1 to 9 digits.
C  The problem is read on standard input  using default device 4.
C  Other input is read on standard input using default device 5
C  and error messages are written on the error output, default device 0.
C
C  The database file is read (the default name is CLASSF.DB but this can be
C  changed by the user), the new problem classification is then inserted in
C  an updated copy of the database (the default name is CLASSF.UDB but this
C  can also be changed by the user).
C
C  Programming:  A.R. Conn and Ph. Toint, Aug 1992, for CGT Productions.
C
C--------- THE FOLLOWING SPECIFICATIONS MAY BE MODIFIED BY THE USER ----------
C
C  The SIF file is read on standard input using default device 4.
C  Standard input default definition is device 5.  Standard output default
C  definition is device 6. Standard error default definition is device 6.
C  Change them to whatever values are appropriate on your system.
C
      INTEGER      NSTDIN,     STDIN,     STDOUT,     STDERR
C                            (Units may not be 0 in Salford FORTRAN.)
C                            Change by Kristjan Jonasson, March 1993.
      PARAMETER  ( NSTDIN = 4, STDIN = 5, STDOUT = 6, STDERR = 6)

C
C  Classification database file name
C
      CHARACTER*32 DBNAME
      PARAMETER  ( DBNAME = 'CLASSF.DB' )
C
C  Input string for the classification database
C
      INTEGER      DBDVC
      PARAMETER  ( DBDVC = 55 )
C
C  Updated classification database file name
C
      CHARACTER*32 UDBNAM
      PARAMETER  ( UDBNAM = 'CLASSF.UDB' )
C
C  Output string for the updated database
C
      INTEGER      UDBDVC
      PARAMETER  ( UDBDVC = 56 )
C
C  Input string for mode of operation and problem name
C
      INTEGER      MPDVC
      PARAMETER  ( MPDVC = 57 )
      CHARACTER*32 MPNAM
      PARAMETER  ( MPNAM = 'CLSF.DAT' )
C
C---------------- END OF THE USER MODIFIABLE SPECIFICATION ------------------
C
C  Variable definitions
C
      CHARACTER*80 LINE
      CHARACTER*32 IFNAME
      CHARACTER*36 CLS, CLSDB
      CHARACTER*8  NAME, NAMEDB
      CHARACTER*1  CHOICE, OBJ, CON, SMTH, INTRST, INTVAR, UPPER,
     1             VARN, VARM
      LOGICAL      FOUND, WRITTN, AUTO, OVERW
      INTEGER      H, ILEN, IDER, K, NCON, NEWH, NVAR, START, TSIZE
C
C  General initializations
C
      FOUND  = .FALSE.
      AUTO   = .FALSE.
      OVERW  = .FALSE.
      NAME   = '?'
      OBJ    = '?'
      CON    = '?'
      SMTH   = '?'
      INTRST = '?'
      INTVAR = '?'
      VARN   = '?'
      VARM   = '?'
      IDER   = 0
      NVAR   = 0
      NCON   = 0
C
C  Open the mode and problem file
C
      OPEN ( UNIT = MPDVC, FILE = MPNAM, STATUS = 'UNKNOWN' )
C
C  Check for automatic execution
C
      READ  ( MPDVC, '( A1 )' ) CHOICE
      IF ( UPPER( CHOICE ) .EQ. 'Y' ) AUTO = .TRUE.
C
C  Input the initial name of the sif file with the classification you want
C  to add
C
 2    CONTINUE
      READ  ( MPDVC , FMT = '( A )', ERR = 4 ) IFNAME
      CLOSE ( MPDVC )
C                            Inspired by Kristjan Jonasson, March 1993.
C Since sed (or basename) does not exist in DOS, the following
C 9 lines allow that IFNAME be specified with or without the .SIF
C extension. This means that CLASSIFY.BAT and CLASSALL.BAT work.
      TSIZE=INDEX(IFNAME,' ')
      IF (TSIZE .GT. 4 .AND. IFNAME(TSIZE-4:TSIZE-1) .EQ. '.SIF') THEN
        TSIZE=TSIZE-4
      ELSE
        IFNAME(TSIZE:TSIZE) = '.'
        IFNAME(TSIZE+1:TSIZE+1) = 'S'
        IFNAME(TSIZE+2:TSIZE+2) = 'I'
        IFNAME(TSIZE+3:TSIZE+3) = 'F'
      ENDIF
      WRITE ( STDOUT, 6000 ) IFNAME(1:TSIZE+3)
      GO TO 100
 4    CONTINUE
      WRITE ( STDOUT, 1201 )
      GO TO 2
C
C  read the next line if not already at the end of the GROUPS section
C
100   CONTINUE
      OPEN ( UNIT = NSTDIN, FILE = IFNAME, STATUS = 'OLD' )
101   CONTINUE
      IF ( FOUND ) GO TO 200
      READ( NSTDIN, '(A80)', END = 99) LINE
C
C  exit when the GROUPS or VARIABLES section is encountered
C
      IF ( LINE(1:6) .EQ. 'GROUPS'     .OR.
     1     LINE(1:7) .EQ. 'COLUMNS'    .OR.
     1     LINE(1:9) .EQ. 'VARIABLES'  .OR.
     1     LINE(1:4) .EQ. 'ROWS'            ) GO TO 99
C
C  remove trailing blanks and compute line length
C
      DO 10 ILEN = 80, 1, -1
        IF ( LINE(ILEN:ILEN) .NE. ' ' ) GO TO 20
10    CONTINUE
      ILEN = 0
20    CONTINUE
C
C  find the problem's name
C
      IF ( LINE(1:4) .EQ. 'NAME' ) NAME = LINE(15:24)
C
C  skip non-comment lines
C
      IF ( LINE(1:1) .EQ. '*' ) THEN
C
C  find the position of the first word in the comment line
C
       DO 30 K = 2, ILEN
          IF ( LINE(K:K) .NE. ' ' ) GO TO 40
 30    CONTINUE
        K = 1
 40     CONTINUE
C
C  is it the classification ?
C
        IF ( K .GT. 1  .AND. LINE(K:K+13) .EQ. 'CLASSIFICATION' .OR.
     1       LINE(K:K+13) .EQ. 'classification' ) THEN
          FOUND = .TRUE.
C
C  find the classification string
C
          DO 50 START = K + 15, ILEN
             IF ( LINE(START:START) .NE. ' ' ) GO TO 60
 50       CONTINUE
          GO TO 300
 60       CONTINUE
C
C  objective classification
C
          OBJ  = LINE(START:START)
          IF ( OBJ .NE. 'N' .AND.
     1         OBJ .NE. 'C' .AND.
     1         OBJ .NE. 'L' .AND.
     1         OBJ .NE. 'Q' .AND.
     1         OBJ .NE. 'S' .AND.
     1         OBJ. NE. 'O'      ) THEN
            WRITE ( STDERR, 1001 ) NAME
          ENDIF
C
C  constraints classification
C
          CON  = LINE(START+1:START+1)
          IF ( CON .NE. 'U' .AND.
     1         CON .NE. 'X' .AND.
     1         CON .NE. 'B' .AND.
     1         CON .NE. 'N' .AND.
     1         CON .NE. 'L' .AND.
     1         CON .NE. 'Q' .AND.
     1         CON. NE. 'O'      ) THEN
            WRITE ( STDERR, 1002 ) NAME
          ENDIF
C
C  smoothness classification
C
          SMTH = LINE(START+2:START+2)
          IF ( SMTH .NE. 'R' .AND. SMTH .NE. 'I' ) THEN
            WRITE ( STDERR, 1003 ) NAME
          ENDIF
C
C  derivatives degree
C
          CALL GETINT( STDERR, LINE, START + 3, H, IDER )
          IF ( IDER. NE. 0 .AND. IDER .NE. 1 .AND. IDER .NE. 2 ) THEN
            WRITE ( STDERR, 1004 ) NAME
          ENDIF
          IF ( H .EQ. 80 ) THEN
            GO TO 300
          ENDIF
C
C  interest of the problem
C
          INTRST = LINE(H+1:H+1)
          IF ( INTRST .NE. 'A' .AND.
     1         INTRST .NE. 'M' .AND.
     1         INTRST .NE. 'R'      ) THEN
            WRITE ( STDERR, 1005 ) NAME
          ENDIF
C
C  internal variables?
C
          INTVAR = LINE(H+2:H+2)
          IF ( INTVAR .NE. 'Y' .AND. INTVAR .NE. 'N' ) THEN
            WRITE ( STDERR, 1006 ) NAME
          ENDIF
C
C  number of variables
C
          VARN = LINE(H+4:H+4)
          IF ( VARN .NE. 'V' ) THEN
            CALL GETINT( STDERR, LINE, H + 4, NEWH, NVAR )
            H = NEWH
            IF ( NVAR .LE. 0 ) THEN
              WRITE ( STDERR, 1007 ) NAME
              NVAR = 0
            ENDIF
            VARN = 'F'
            IF ( NEWH .EQ. 80 ) THEN
              GO TO 300
            ENDIF
          ELSE
            NVAR = 0
            H = H + 5
          ENDIF
C
C  number of constraints
C
          VARM = LINE(H+1:H+1)
          IF ( VARM .NE. 'V' ) THEN
            CALL GETINT( STDERR, LINE, H + 1, NEWH, NCON )
            IF ( NCON .LT. 0 ) THEN
              WRITE ( STDERR, 1008 ) NAME
              NCON = 0
            ENDIF
            VARM = 'F'
          ELSE
            NCON = 0
          ENDIF
        ENDIF
      ENDIF
C
C  loop for next line
C
      GO TO 101
C
C  incomplete classification
C
 300  CONTINUE
      WRITE ( STDERR, 1009) NAME
C
C  end of file or next SIF section encountered: missing classification
C
  99  CONTINUE
      WRITE ( STDERR, 1010 ) NAME
      STOP
C
C  format the new classification
C
 200  CONTINUE
      CLS = '                            '
      CLS(1:8) = NAME
      CLS(10:10) = OBJ
      CLS(11:11) = CON
      CLS(12:12) = SMTH
      WRITE ( CLS(13:13), '( I1 )' ) IDER
      CLS(14:14) = '-'
      CLS(15:15) = INTRST
      CLS(16:16) = INTVAR
      CLS(17:17) = '-'
      IF ( VARN .EQ. 'V' ) THEN
        CLS(26:26) = VARN
      ELSE
        WRITE ( CLS(18:26), '( I9 )' ) NVAR
      ENDIF
      CLS(27:27) = '-'
      IF ( VARM .EQ. 'V' ) THEN
        CLS(36:36) = VARM
      ELSE
        WRITE ( CLS(28:36), '( I9 )' ) NCON
      ENDIF
C
C  open the classification database
C
      OPEN ( UNIT = DBDVC, FILE = DBNAME, STATUS = 'UNKNOWN' )
C
C  open the updated classification database, according to the
C  conventions of the local operating system
C
      OPEN ( UNIT = UDBDVC, FILE = UDBNAM, STATUS = 'UNKNOWN' )
C
C  read the current line in the existing database
C
      WRITTN = .FALSE.
 400  CONTINUE
      READ ( DBDVC, '( A36 )', END = 900 ) CLSDB
      NAMEDB = CLSDB(1:8)
C
C  verify that the name of the problem to be classified is not already
C  present in the database
C
      IF ( NAMEDB .EQ. NAME ) THEN
        WRITE( STDOUT, 2000 ) NAMEDB
        IF ( .NOT. AUTO ) THEN
          IF ( CLSDB .EQ. CLS ) THEN
             WRITE( STDOUT, 2003 ) DBNAME
          ELSE
             WRITE( STDOUT, 2001 )
             READ  ( STDIN, '( A1 )' ) CHOICE
             IF ( UPPER( CHOICE ) .EQ. 'Y' ) THEN
               OVERW = .TRUE.
             ENDIF
          ENDIF
        ENDIF
      ENDIF
C
C  write the new problem classification, if it is not already present
C  in the existing database or if overwrite is specified
C
      IF ( WRITTN .OR. NAMEDB .LT. NAME ) THEN
        WRITE ( UDBDVC, '( A36 )' ) CLSDB
      ELSEIF ( NAMEDB .EQ. NAME ) THEN
        IF ( CLSDB .GT. CLS .AND. .NOT. OVERW) THEN
          WRITE ( UDBDVC, '( A36 )' ) CLS
          WRITTN = .TRUE.
          WRITE ( UDBDVC, '( A36 )' ) CLSDB
        ELSEIF ( CLSDB .LT. CLS .AND. .NOT. OVERW) THEN
          WRITE ( UDBDVC, '( A36 )' ) CLSDB
          WRITE ( UDBDVC, '( A36 )' ) CLS
          WRITTN = .TRUE.
        ELSE
          WRITE ( UDBDVC, '( A36 )' ) CLS
          WRITTN = .TRUE.
        ENDIF
      ELSE
        WRITE ( UDBDVC, '( A36 )' ) CLS
        WRITE ( UDBDVC, '( A36 )' ) CLSDB
        WRITTN = .TRUE.
      ENDIF
C
C  read next line of the existing database
C
      GO TO 400
C
C  close the database files
C
 900  CONTINUE
      CLOSE ( DBDVC )
      IF ( .NOT. WRITTN ) WRITE ( UDBDVC, '( A36 )' ) CLS
      CLOSE ( UDBDVC )
      STOP
C
C  non excutable statements
C
 1001 FORMAT( ' Problem ', A8,
     1        ' has impossible objective classification!' )
 1002 FORMAT( ' Problem ', A8,
     1        ' has impossible constraints classification!' )
 1003 FORMAT( ' Problem ', A8,
     1        ' has impossible smoothnes classification' )
 1004 FORMAT( ' Problem ', A8,
     1        ' has impossible degree of available derivatives!' )
 1005 FORMAT( ' Problem ', A8,
     1        ' has impossible problem interest specification!' )
 1006 FORMAT( ' Problem ', A8,
     1        ' has  internal variable specification!' )
 1007 FORMAT( ' Problem ', A8,
     1        ' has impossible number of variables!' )
 1008 FORMAT( ' Problem ', A8,
     1        ' has impossible number of constraints!' )
 1009 FORMAT( ' Incomplete classification for problem ', A8 )
 1010 FORMAT( ' Classification not found for Problem ', A8 )
 1011 FORMAT( ' ', A8, ' ', 60A1 )
 1201 FORMAT( '   *** YOU USED MORE THAN 28 CHARACTERS.  Please choose',
     1             ' again.'
     1        / )
 2000 FORMAT( ' Duplicate name: ', A8 )
 2001 FORMAT( ' Do you wish to overwrite the existing classification?',
     1        '  [<CR> = N] ? (N/Y)' )
 2003 FORMAT( ' Identical classification record already appears in ',
     1                A32 )
 6000 FORMAT( /' Your current SIF filename file is : ',
     1                A )
C
      END
C
C
C
      SUBROUTINE GETINT ( STDERR, LINE, START, HYPHEN, I )
C
C-----------------------------------------------------------------------------
C
C  The purpose of this subroutine is to read, in the string LINE, an
C  integer number, from the position START to the first following hyphen
C  or EOL. The position of the hyphen or EOL is returned in HYPHEN, while
C  the integer is returned in I.
C
C  Programming: Ph. Toint, Dec 1991, for CGT Productions.
C
C-----------------------------------------------------------------------------
C
      CHARACTER*80 LINE
      INTEGER      START, HYPHEN, I, STDERR, END
C
C   find the position of the hyphen
C
      END = START - 1
      DO 10 HYPHEN = START, 80
        IF ( LINE( HYPHEN:HYPHEN ) .EQ. '-' ) GO TO 20
        IF ( LINE( HYPHEN:HYPHEN ) .NE. ' ' ) END = END + 1
 10   CONTINUE
 20   CONTINUE
C
C   read the integer
C
      IF ( END .EQ. START ) THEN
        READ ( LINE( START:END ), FMT= '( I1 )', ERR = 30 ) I
      ELSE IF ( END .EQ. START + 1 ) THEN
        READ ( LINE( START:END ), FMT= '( I2 )', ERR = 30 ) I
      ELSE IF ( END .EQ. START + 2 ) THEN
        READ ( LINE( START:END ), FMT= '( I3 )', ERR = 30 ) I
      ELSE IF ( END .EQ. START + 3 ) THEN
        READ ( LINE( START:END ), FMT= '( I4 )', ERR = 30 ) I
      ELSE IF ( END .EQ. START + 4 ) THEN
        READ ( LINE( START:END ), FMT= '( I5 )', ERR = 30 ) I
      ELSE IF ( END .EQ. START + 5 ) THEN
        READ ( LINE( START:END ), FMT= '( I6 )', ERR = 30 ) I
      ELSE IF ( END .EQ. START + 6 ) THEN
        READ ( LINE( START:END ), FMT= '( I7 )', ERR = 30 ) I
      ELSE IF ( END .EQ. START + 7 ) THEN
        READ ( LINE( START:END ), FMT= '( I8 )', ERR = 30 ) I
      ELSE
        READ ( LINE( START:END ), FMT= '( I9 )', ERR = 30 ) I
      ENDIF
      RETURN
C
C   error in reading the integer
C
 30   CONTINUE
      WRITE ( STDERR, 1000 )
      RETURN
C
C   non executable statements
C
 1000 FORMAT(' Error in reading an integer!' )
      END
C
C-----------------------------------------------------------------------------
C
      CHARACTER*1 FUNCTION UPPER( CH )
C
C-----------------------------------------------------------------------------
C
C  The purpose of this function is to transform the character CH
C  to upper case, if it is not already.
C
C  Programming: Ph. Toint, Dec 1991, for CGT Productions.
C
C-----------------------------------------------------------------------------
C
      CHARACTER*1  CH
C
      INTRINSIC    ICHAR, CHAR
C
      INTEGER      ICH, LSTART
C
      LSTART = ICHAR( 'a' )
      ICH    = ICHAR( CH )
      IF ( ICH .GE. LSTART .AND. ICH .LE. ICHAR( 'z' ) ) THEN
        UPPER = CHAR( ICHAR( 'A' ) + ICH - LSTART )
      ELSE
        UPPER = CH
      ENDIF
      RETURN
      END

C
C-----------------------------------------------------------------------------
C

ext: f

max_forks_repo_name: dynamics-of-stellar-systems/dynamite


!-----------------------------------------------------------------------
      program test
      implicit none
      include 'ceedf.h'

      integer ceed,err
      integer u,v
      integer qf
      integer q,i
      parameter(q=8)
      real*8 uu(q)
      real*8 vv(q)
      character arg*32
      integer*8 uoffset,voffset

      call getarg(1,arg)
      call ceedinit(trim(arg)//char(0),ceed,err)

      call ceedqfunctioncreateidentity(ceed,1,ceed_eval_interp,&
     & ceed_eval_interp,qf,err)

      do i=0,q-1
        uu(i+1)=i*i
      enddo

      call ceedvectorcreate(ceed,q,u,err)
      uoffset=0
      call ceedvectorsetarray(u,ceed_mem_host,ceed_use_pointer,uu,uoffset,err)
      call ceedvectorcreate(ceed,q,v,err)
      call ceedvectorsetvalue(v,0.d0,err)

      call ceedqfunctionapply(qf,q,u,ceed_null,ceed_null,ceed_null,&
             &ceed_null,ceed_null,ceed_null,ceed_null,ceed_null,ceed_null,&
             &ceed_null,ceed_null,ceed_null,ceed_null,ceed_null,ceed_null,&
             &v,ceed_null,ceed_null,ceed_null,ceed_null,ceed_null,&
             &ceed_null,ceed_null,ceed_null,ceed_null,ceed_null,ceed_null,&
             &ceed_null,ceed_null,ceed_null,ceed_null,err)

      call ceedvectorgetarrayread(v,ceed_mem_host,vv,voffset,err)
      do i=1,q
        if (abs(vv(i+voffset)-(i-1)*(i-1)) > 1.0D-14) then
! LCOV_EXCL_START
          write(*,*) 'v(i)=',vv(i+voffset),', u(i)=',(i-1)*(i-1)
! LCOV_EXCL_STOP
        endif
      enddo
      call ceedvectorrestorearrayread(v,vv,voffset,err)

      call ceedvectordestroy(u,err)
      call ceedvectordestroy(v,err)
      call ceedqfunctiondestroy(qf,err)
      call ceeddestroy(ceed,err)
      end
!-----------------------------------------------------------------------


text:
      ASSIGN X TO 2
100   CONTINUE
      PRINT X ^3
      X = X + 1
      GOTO 100
      END

max_stars_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler

max_issues_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler

max_forks_repo_name: ElliotAlexander/EDSAC-FORTRAN-Compiler


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1995  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ############################################################
c     ##                                                        ##
c     ##  subroutine getmol2  --  get a Sybyl MOL2 format file  ##
c     ##                                                        ##
c     ############################################################
c
c
c     "getmol2" asks for a Sybyl MOL2 molecule file name,
c     then reads the coordinates from the file
c
c
      subroutine getmol2
      implicit none
      include 'files.i'
      include 'iounit.i'
      integer isyb
      integer freeunit
      logical exist
      character*120 sybylfile
c
c
c     try to get a filename from the command line arguments
c
      call nextarg (sybylfile,exist)
      if (exist) then
         call basefile (sybylfile)
         call suffix (sybylfile,'mol2','old')
         inquire (file=sybylfile,exist=exist)
      end if
c
c     ask for the user specified input structure filename
c
      do while (.not. exist)
         write (iout,10)
   10    format (/,' Enter a Sybyl MOL2 File Name :  ',$)
         read (input,20)  sybylfile
   20    format (a120)
         call basefile (sybylfile)
         call suffix (sybylfile,'mol2','old')
         inquire (file=sybylfile,exist=exist)
      end do
c
c     first open and then read the Sybyl MOL2 coordinates file
c
      isyb = freeunit ()
      open (unit=isyb,file=sybylfile,status='old')
      rewind (unit=isyb)
      call readmol2 (isyb)
      close (unit=isyb)
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
      subroutine gr2 (neq, t, y, ng, groot)
      integer neq, ng
      double precision t, y, groot
      dimension y(2), groot(1)
      dimension neq(*)
      groot(1) = y(1)
      return
      end

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
      subroutine DUSETN ( Nin, M1in, M2in)
c  Base name of file: DUCOMP
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 2001-06-18 DUCOMP Krogh  Replaced "M1 .eq.. 0" with "M1 .eq. 0"
C>> 1996-04-30 DUCOMP Krogh  Removed redundant save statement.
C>> 1994-11-02 DUCOMP Krogh  Changes to use M77CON
C>> 1994-08-04 DUCOMP CLL New subrs: [D|S]USETN & [D|S]UGETN.
C>> 1993-11-11 CLL Changed Entries to separate Subroutines.
C>> 1990-01-23 CLL Replace M with MPWR in call to IERM1
C>> 1987-12-07 DUCOMP Lawson  Initial code.
c
c  The file DUCOMP contains a set of program units to perform
c  computation on arrays representing the value of a quantity
c  and (optionally) its first and second partial derivatives
c  with respect to a set of N independent variables.
c
c  The lowest and highest order derivs desired are specified by
c  the integers M1 and M2 which must satisfy
c               0 .le. M1 .le. M2 .le. 2
c  The integers N, M1, and M2 must be set by the user by a call
c  to [D|S]USETN, which will put these values in the COMMON block UCOM1.
c
c  Each array containing a U-variable for which second partial
c  derivs are to be carried must be dimensioned at least
C  ((N+2)*(N+1))/2    EXAMPLES..
C              N = 1  2  3  4  5  6  7  8  9 10 11 12  13  14  15  16
C      DIMENSION = 3  6 10 15 21 28 36 45 55 66 78 91 105 120 136 153
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c        Subroutines
c
c    subroutine DUSETN ( N, M1, M2)
c    subroutine DUGETN ( N, M1, M2, L1, L2)
c    subroutine DUSET  ( VAL, KEY, Y9)
c
c    subroutine DUPRO  ( U,V,Y)
c    subroutine DUQUO  ( U1,V1,Y1)
c    subroutine DUSUM  ( U3,V3,Y3)
c    subroutine DUDIF  ( U4,V4,Y4)
c    subroutine DUSUM1 ( C5,V5,Y5)
c    subroutine DUDIF1 ( C6,V6,Y6)
c    subroutine DUPRO1 ( C7,V7,Y7)
c    subroutine DUQUO1 ( C8,V8,Y8)
c
c    subroutine DUSQRT (U,Y)
c    subroutine DUEXP  (U1,Y1)
c    subroutine DULOG  (U2,Y2)
c    subroutine DUPWRI (MPWR,U3,Y3)
c
c    subroutine DUSIN  (U,Y)
c    subroutine DUCOS  (U ,Y )
c    subroutine DUSINH (U ,Y )
c    subroutine DUCOSH (U ,Y )
c    subroutine DUATAN (U1,Y1)
c    subroutine DUATN2 (V2,U2,Y2)
c    subroutine DUASIN (U, Y)
c    subroutine DUACOS (U, Y)
c    subroutine DUTAN  (U, Y)
c    subroutine DUTANH (U, Y)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c--D replaces "?": ?UCOMP, ?UACOS, ?UASIN, ?UATAN, ?UATN2, ?UCOS
c--   & ?UCOSH, ?UDIF, ?UDIF1, ?UEXP, ?UGETN, ?ULOG, ?UPRO, ?UPRO1
c--   & ?UPWRI, ?UQUO, ?UQUO1, ?USET, ?USETN, ?USIN, ?USINH, ?USQRT
c--   & ?USUM, ?USUM1, ?UTAN, ?UTANH, ?UACS
c     Subprograms referecnced by both versions: ERMSG
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
c     C. L. Lawson, JPL, 1969 Dec 4
C     References:
c     1.   Wengert, R. E., A simple automatic derivative evaluation
c          program, Comm. ACM, 1, Aug 1964, 463-464.
c     2.   C. L. Lawson, Computing Derivatives using W-arithmetic and
c          U-arithmetic., JPL Appl Math TM 289, Sept 1971.
c     Revised by CLL for Fortran 77 in Jan 1987, Sept 1987.
c     ==================================================================
c     subroutine DUSETN ( Nin, M1in, M2in)
c
c  Nin [in]  Specifies the number of independent variables.
c     Require Nin .ge. 1.
c  M1in, M2in [in]  These specify the lowest and highest order derivs
c     desired.  These must satisfy 0 .le. M1in .le. M2in .le. 2.
c
c  This subr copies Nin, M1in, and M2in into common variables N, M1, and
c  M2.
c  It also computes L1 and L2 based on N, M1, and M2, and stores
c  L1 and L2 into common.  L1 and L2 specify the first and last
c  locations in a U-variable array that will be operated on when
c  dealing with derivs of orders M1 through M2 and
c  N independent variables.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer Nin, M1in, M2in
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      integer L1, L2
      common/UCOM2/L1,L2
      save /UCOM1/,/UCOM2/
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      N = Nin
      M1 = M1in
      M2 = M2in
C                                      SET L1 AND L2
      if(M1 .eq. 0) then
         L1 = 1
      elseif(M1 .eq. 1) then
         L1 = 2
      else
         L1 = N+2
      endif

      if(M2 .eq. 0) then
         L2 = 1
      elseif(M2 .eq. 1) then
         L2 = N+1
      else
         L2=1+N+((N*(N+1))/2)
      endif
      end
c     ==================================================================
      subroutine DUGETN ( Nout, M1out, M2out, L1out, L2out)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer Nout, M1out, M2out, L1out, L2out
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      integer L1, L2
      common/UCOM2/L1,L2
      save /UCOM1/,/UCOM2/
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      Nout = N
      M1out = M1
      M2out = M2
      L1out = L1
      L2out = L2
      end
c     ==================================================================
      subroutine DUSET( VAL, KEY, Y9)
c
c     A convenience routine for assigning an initial value to the
c     U-variable, Y9().  The parts of Y9() assigned will
c     depend on M1 and M2.
c
c     If M1 .le. 0 .le. M2, Y9() will be assigned the value, VAL.
c
c     If M1 .le. 1 .le. M2, the first partial derivs of Y9() will be
c     assigned depending on KEY.  KEY must satisfy 0 .le. KEY .le. N.
c     If KEY .eq.. 0, the U-variable is to be constant relative to the N
c     independent variables.  Thus all first partial derivs will be set
c     to zero.
c     If KEY is in [1,N], Y9() is being set as th KEYth independent
c     variable.  Thus the KEYth first partial deriv will be set to 1.0
c     and all other first partial derivs will be set to zero.
c
c     If M1 .le. 2 .le. M2, the second partial derivs of Y9() will be
c     set to zero.
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer KEY
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision VAL, Y9(*)
      double precision       ONE, ZERO
      parameter( ONE = 1.0D0, ZERO = 0.0D0)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if(M1-1)230,240,250
C                                      VALUE
  230 Y9(1)=VAL
      if(M2 .eq. 0) return
C                                      1ST PARTIALS
  240 continue
          do 245 I=1,N
  245     Y9(I+1) = ZERO
      if(KEY .ge. 1 .and. KEY .le. N) Y9(KEY+1) = ONE
      if(M2  .eq. 1) return
C                                      2ND PARTIALS
  250 continue
      K=N+2
          do 254 I=1,N
               do 252 J=1,I
               Y9(K)= ZERO
  252          K=K+1
  254     continue
      return
      end
c     ==================================================================
      subroutine DUPRO( U, V, Y)
c     Computes Y = U * V
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision       U(*),V(*),Y(*)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if(M2-1)60,40,10
C                                      2ND PARTIALS
   10 continue
      K=N+2
          do 30 I=1,N
               do 20 J=1,I
               Y(K)=U(1)*V(K)+U(K)*V(1)+U(I+1)*V(J+1)+U(J+1)*V(I+1)
               K=K+1
   20          continue
   30     continue
      if(M1 .eq. 2) return
C                                      1ST PARTIALS
   40 continue
          do 50 I=1,N
   50     Y(I+1)=U(1)*V(I+1) + U(I+1)*V(1)
      if(M1 .eq. 1) return
C                                      VALUE
   60 continue
      Y(1)=U(1)*V(1)
      return
      end
C
C
      subroutine DUQUO ( U1,V1,Y1)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision       U1(*),V1(*),Y1(*)
      double precision       VINV
      double precision       ONE, ZERO
      parameter( ONE = 1.0D0, ZERO = 0.0D0)
C
C     COMPUTE..        Y=U/V
C
      if(V1(1) .eq. ZERO) THEN
         call ERMSG('DUQUO',5,0,'Denominator is zero.','.')
         return
      end if
      if(M1-1)80,90,100
C                                      VALUE
   80 continue
      Y1(1)=U1(1)/V1(1)
      if(M2 .eq. 0) return
C                                      1ST PARTIAL DERIVS
   90 continue
      VINV=ONE/V1(1)
          do 95 I=1,N
   95     Y1(I+1)=(U1(I+1)-Y1(1)*V1(I+1))*VINV
      if(M2 .eq.1) return
      go to 110
C                                      2ND PARTIAL DERIVS
  100 VINV=ONE/V1(1)
  110 continue
      K=N+2
          do 114 I=1,N
               do 112 J=1,I
               Y1(K)= (U1(K)- Y1(1)*V1(K)-Y1(I+1)*V1(J+1)
     *                                   -Y1(J+1)*V1(I+1))*VINV
               K=K+1
  112          continue
  114     continue
      return
      end
C
C
C
      subroutine DUSUM  ( U3,V3,Y3)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer K
      integer L1, L2
      common/UCOM2/L1, L2
      save /UCOM2/
      double precision       U3(*),V3(*),Y3(*)
C
C     COMPUTE          Y=U+V      AND PARTIAL DERIVS
C
      do 148 K = L1,L2
  148     Y3(K)=U3(K)+V3(K)
      return
      end
C
C
      subroutine DUDIF  ( U4,V4,Y4)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer K
      integer L1, L2
      common/UCOM2/L1, L2
      save /UCOM2/
      double precision       U4(*),V4(*),Y4(*)
C
C     COMPUTE          Y=U-V      AND PARTIAL DERIVS
C
      do 154 K=L1,L2
  154     Y4(K)=U4(K)-V4(K)
      return
      end
C
C
      subroutine DUSUM1 ( C5,V5,Y5)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer K
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      integer L1, L2
      common/UCOM2/L1, L2
      save /UCOM1/,/UCOM2/
      double precision       C5,V5(*),Y5(*)
C
C     COMPUTE         Y=C+V   WHERE C IS A CONSTANT
C
      do 160 K=L1,L2
  160    Y5(K)=   V5(K)
      if(M1.eq.0) Y5(1)=C5+Y5(1)
      return
      end
C
C
C
      subroutine DUDIF1 ( C6,V6,Y6)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer K
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      integer L1, L2
      common/UCOM2/L1, L2
      save /UCOM1/,/UCOM2/
      double precision       C6,V6(*),Y6(*)
C
C     COMPUTE          Y=C-V  WHERE C IS A CONSTANT
C
      do 168 K=L1,L2
  168     Y6(K)=  -V6(K)
c
c         Using + in following statement because have just set
c         Y6(1) = -V6(1).
c
      IF (M1.eq.0) Y6(1)=C6+Y6(1)
      end
C
C
C
      subroutine DUPRO1( C7,V7,Y7)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer K
      integer L1, L2
      common/UCOM2/L1, L2
      save /UCOM2/
      double precision       C7,V7(*),Y7(*)
C
C     COMPUTE          Y=C*V  WHERE C IS A CONSTANT
C
      do 176 K=L1,L2
  176     Y7(K)=C7*V7(K)
      return
      end
C     ==================================================================
      subroutine DUQUO1(C8,V8,Y8)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision       C8, FAC, V8(*),Y8(*)
      double precision       ONE, ZERO
      parameter( ONE = 1.0D0, ZERO = 0.0D0)
C
C     COMPUTE          Y=C/V  WHERE C IS A CONSTANT
C
      if(V8(1) .eq. ZERO) then
         call ERMSG('DUQUO1',1,0,'Denominator is zero.','.')
         return
      endif
      if(M1-1)204,208,214
C                                      VALUE
  204 Y8(1)=C8/V8(1)
      if(M2 .eq. 0) return
C                                      1ST PARTIALS
  208 continue
      FAC = -Y8(1)/V8(1)
          do 212 I=1,N
  212     Y8(I+1) = V8(I+1)*FAC
      if(M2  .eq. 1) return
C                                      2ND PARTIALS
  214 FAC=-ONE/V8(1)
      K=N+2
          do 224 I=1,N
               do 220 J=1,I
               Y8(K)=(Y8(1)*V8(K)+Y8(I+1)*V8(J+1)+Y8(J+1)*V8(I+1))*FAC
  220          K=K+1
  224     continue
      return
      end
c     ------------------------------------------------------------------
      subroutine DUSQRT (U,Y)
C
C     COMPUTE..                                  Y=SQRT(U)
C
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer I, J, K
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision U(*),Y(*)
      double precision D1, FAC, FAC2
      double precision ZERO, ONE, HALF
      parameter(ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if(M1-1) 4, 6,16
    4 Y(1)= SQRT(U(1))
      if(M2 .eq. 0) return
C
    6 if(Y(1) .eq. ZERO) then
         call ERMSG('DUSQRT',2,0,
     *      'Deriv of sqrt(x) is infinite at x = 0','.')
         return
      endif
c
      D1 = HALF/Y(1)
          do 12 I=1,N
   12     Y(I+1)=U(I+1)*D1
      if(M2 .eq. 1)return
      go to 18
C
   16 continue
      D1 = HALF/Y(1)
   18 K=N+1
      FAC=ONE/Y(1)
          do 24 I=1,N
            FAC2 = FAC * Y(I+1)
               do 20 J=1,I
               K=K+1
   20          Y(K) = D1*U(K) - Y(J+1)*FAC2
   24     continue
      return
      end
C
C
C
      subroutine DUEXP  (U1,Y1)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision U1(*),Y1(*)
C                                                Y=EXP(U)
      if(M1-1)30,31,34
   30 Y1(1)= EXP(U1(1))
      if(M2 .eq. 0)return
C
   31 continue
          do 32 I=1,N
   32     Y1(I+1)=Y1(1)*U1(I+1)
      if(M2 .eq. 1)return
C
   34 continue
      K=N+1
          do 36 I=1,N
               do 36 J=1,I
               K=K+1
   36          Y1(K)=Y1(1) * (U1(K) + U1(I+1)*U1(J+1))
      return
      end
C
C
C
      subroutine DULOG  (U2,Y2)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision D1, U2(*),Y2(*)
      double precision ONE
      parameter( ONE = 1.0D0)
C                                      Y=LOG(U)
      if(M1-1)40,42,48
   40 Y2(1) = LOG(U2(1))
      if(M2 .eq. 0) return
C
   42 D1 = ONE/U2(1)
          do 44 I=1,N
   44     Y2(I+1)=U2(I+1) * D1
      if(M2 .eq. 1)return
      go to 49
C
   48 D1 = ONE/U2(1)
   49 continue
      K=N+1
          do 50 I=1,N
               do 50 J=1,I
               K=K+1
   50          Y2(K) = D1 * U2(K) - Y2(I+1) * Y2(J+1)
      return
      end
C
C
C
      subroutine DUPWRI (MPWR,U3,Y3)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, K, K1, K2
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      integer MPWR
      double precision FAC, FAC2, FM, U3(*),Y3(*)
      double precision ZERO, ONE, TWO
      parameter(ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0)
C                                      Y=U**MPWR    MPWR IS A CONSTANT
C     MPWR MAY BE POS.,ZERO, OR NEG.  IF MPWR = 0 THEN Y=1. INDEPENDENT
C     OF U.  IF MPWR IS NEG. THEN ERROR STOP IF U IS ZERO.
C
      if(MPWR  .eq. 0   ) go to 56
      if(U3(1) .NE. ZERO) go to 100
C                                      U = 0.  OR  MPWR = 0
      if(MPWR .lt. 0) then
         call IERM1('DUPWRI',4,0,
     *   'U**M is infinite when U = 0. and M < 0','M',MPWR,'.')
         return
      endif
C                            U=0. AND MPWR .GE. 0
   56 if(M1-1)64,68,84
   64 Y3(1)=ZERO
      if(MPWR .eq. 0) Y3(1)=ONE
      if(M2 .eq. 0) return
C
   68 if(MPWR .eq. 1) go to 74
          do 70 I=1,N
   70     Y3(I+1)=ZERO
      go to 80
   74     do 78 I=1,N
   78     Y3(I+1)=U3(I+1)
   80 if(M2 .eq. 1) return
C
   84 continue
      if(MPWR .eq. 2) then
         K=N+2
          do 94 I=1,N
               do 94 J=1,I
               Y3(K)=TWO*U3(I+1)*U3(J+1)
   94          K=K+1
         return
      endif

      K1 = N+1
      K2=K1-1+(N*(N+1))/2
      if(MPWR .eq. 1) then
         do 98 K=K1,K2
   98    Y3(K)=U3(K)
      else
         do 86 K=K1,K2
   86    Y3(K)=ZERO
      endif
      return
C                                 END..  U .eq. 0. .OR. MPWR .eq. 0
c
C                                 BEGIN..  U .NE. 0. .AND. MPWR .NE.0
  100 continue
      if(M1-1)104,108,112
  104 Y3(1) = U3(1)**MPWR
      if(M2 .eq. 0) return
C
  108 FM=MPWR
      FAC=FM*(Y3(1)/U3(1))
          do 110 I=1,N
  110     Y3(I+1) = FAC * U3(I+1)
      if(M2 .eq. 1) return
      go to 114
  112 FM = MPWR
      FAC= FM*(Y3(1)/U3(1))
  114 FAC2 = (FM-ONE)/U3(1)
      K=N+2
          do 116 I=1,N
               do 116 J=1,I
               Y3(K) = FAC*(FAC2*U3(I+1)*U3(J+1) + U3(K))
  116          K=K+1
      return
C                                  END..  U .NE. 0. .AND. MPWR .NE. 0
      END
c     ------------------------------------------------------------------
      subroutine DUSIN  (U,Y)
C                                       computes   Y=SIN(U)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer I, J, K, N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  D1, D2, FAC, U(*),Y(*)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if(M1 .LE. 0) THEN
         Y(1)= SIN(U(1))
         if(M2 .eq. 0) return
      end if
      D2 = -Y(1)
      D1 = COS(U(1))
      if(M1 .le. 1) go to 16
      go to 24
C
C                                      1ST PARTIALS
   16     do 20 I=1,N
   20     Y(I+1) = D1*U(I+1)
          if(M2 .eq. 1) return
C                                      2ND PARTIALS
   24 K=N+2
          do 28 I=1,N
            FAC = D2*U(I+1)
               do 28 J=1,I
               Y(K)= FAC * U(J+1) + D1*U(K)
   28          K=K+1
      return
      end
C
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine DUCOS  (U ,Y )
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  D1, D2, FAC, U(*),Y(*)
C
C     COMPUTE..                        Y=COS(U)
C
      if(M1 .LE. 0) THEN
         Y(1) = COS(U(1))
         if(M2 .eq. 0) return
      end if
      D2 = -Y(1)
      D1 = -SIN(U(1))
      if(M1 .le. 1) go to 16
      go to 24
C                                      1ST PARTIALS
   16     do 20 I=1,N
   20     Y(I+1) = D1*U(I+1)
          if(M2 .eq. 1) return
C                                      2ND PARTIALS
   24 K=N+2
          do 28 I=1,N
            FAC = D2*U(I+1)
               do 28 J=1,I
               Y(K)= FAC * U(J+1) + D1*U(K)
   28          K=K+1
      return
      end
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine DUSINH  (U ,Y )
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  FAC, D1, D2, U(*),Y(*)
C
C     COMPUTE..                        Y=SINH(U)
C
      if(M1 .LE. 0) THEN
         Y(1) = SINH(U(1))
         if(M2 .eq. 0) return
      end if
      D2 = Y(1)
      D1 = COSH(U(1))
      if(M1 .le. 1) go to 16
      go to 24
C                                      1ST PARTIALS
   16     do 20 I=1,N
   20     Y(I+1) = D1*U(I+1)
          if(M2 .eq. 1) return
C                                      2ND PARTIALS
   24 K=N+2
          do 28 I=1,N
            FAC = D2*U(I+1)
               do 28 J=1,I
               Y(K)= FAC * U(J+1) + D1*U(K)
   28          K=K+1
      return
      end
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine DUCOSH  (U ,Y )
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  FAC, D1, D2, U(*),Y(*)
C
C     COMPUTE..                        Y=COSH(U)
C
      if(M1 .LE. 0) THEN
         Y(1) = COSH(U(1))
         if(M2 .eq. 0) return
      end if
      D2 = Y(1)
      D1 = SINH(U(1))
      if(M1 .le. 1) go to 16
      go to 24
C                                      1ST PARTIALS
   16     do 20 I=1,N
   20     Y(I+1) = D1*U(I+1)
          if(M2 .eq. 1) return
C                                      2ND PARTIALS
   24 K=N+2
          do 28 I=1,N
            FAC = D2*U(I+1)
               do 28 J=1,I
               Y(K)= FAC * U(J+1) + D1*U(K)
   28          K=K+1
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine DUATAN (U1,Y1)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  D1, FAC, FAC2, U1(*),Y1(*)
      double precision       ONE
      parameter(ONE = 1.0D0)
C
C     COMPUTE..                        Y= ATAN(U)
C
      if(M1-1)30,32,36
   30 Y1(1) = ATAN(U1(1))
      if(M2 .eq. 0) return
C
   32 D1 = ONE/(ONE+U1(1)**2)
          do 34 I=1,N
   34     Y1(I+1)=D1*U1(I+1)
      if(M2 .eq. 1) return
      go to 38
C
   36 D1 = ONE/(ONE+U1(1)**2)
   38 FAC = -(U1(1)+U1(1))
      K=N+2
          do 40 I=1,N
            FAC2 = FAC * Y1(I+1)
               do 40 J=1,I
               Y1(K)= D1*U1(K) + FAC2 * Y1(J+1)
   40          K=K+1
      return
      end
C
C
C
      subroutine DUATN2(V2,U2,Y2)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  FAC1, FAC2, R, U2(*),V2(*),Y2(*)
      double precision       ONE
      parameter(ONE = 1.0D0)
C
C     COMPUTE   Y= ATAN2(V,U) = ATAN(V/U)  with 4-quadrant resolution.
C
      if(M1 .eq. 0) then
         Y2(1)= ATAN2(V2(1),U2(1))
         if(M2 .eq. 0) return
      endif
C
      if(abs(V2(1)) .ge. abs(U2(1))) then
         R=U2(1)/V2(1)
         FAC2=ONE/(V2(1)+R*U2(1))
         FAC1=R*FAC2
      else
         R=V2(1)/U2(1)
         FAC1=ONE/(V2(1)*R + U2(1))
         FAC2=R*FAC1
      endif
      if(M1 .eq. 2) go to 54
C                                      FIRST PARTIALS
          do 52 I=1,N
   52     Y2(I+1)=FAC1*V2(I+1)-FAC2*U2(I+1)
      if(M2 .eq. 1) return
C                                      SECOND PARTIALS
   54 K=N+2
          do 56 I=1,N
               do 56 J=1,I
               Y2(K)=FAC1*(V2(K) - U2(I+1)*Y2(J+1) - U2(J+1)*Y2(I+1) )
     *              -FAC2*(U2(K) + V2(I+1)*Y2(J+1) + V2(J+1)*Y2(I+1) )
   56          K=K+1
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine DUASIN (U1,Y1)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  U1(*),Y1(*)
C
C     COMPUTE..                        Y= ASIN(U)
C
      call DUACS (.false., U1,Y1)
      return
      end
C     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine DUACOS (U1,Y1)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  U1(*),Y1(*)
C
C     COMPUTE..                        Y= ACOS(U)
C
      call DUACS (.true., U1,Y1)
      return
      end
c     ==================================================================
      subroutine DUACS (ACOSIN, U1,Y1)
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      logical ACOSIN
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  FAC, FAC2, S1, U1(*),Y1(*)
      double precision       ONE, ZERO
      parameter(ONE = 1.0D0, ZERO = 0.0D0)
C
C     COMPUTE..                    Y= ACOS(U) if ACOSIN .eq. .true. and
C                                  Y= ASIN(U) if ACOSIN .eq. .false.
C
      if(M1 .eq. 0)then
         if(ACOSIN) then
            Y1(1) = ACOS(U1(1))
         else
            Y1(1) = ASIN(U1(1))
         endif
      endif
      if(M2 .eq. 0) return
C
      S1  = ONE - U1(1)**2
      if(S1 .eq. ZERO) then
c                                     Error condition.
         if(ACOSIN) then
            call ERMSG('DUACOS',1,0,
     *      'Deriv of ACOS(x) is infinite at x = -1 or +1','.')
         else
            call ERMSG('DUASIN',1,0,
     *      'Deriv of ASIN(x) is infinite at x = -1 or +1','.')
         endif
         return
      endif
      FAC = ONE/ sqrt(S1)
      if(ACOSIN) FAC = -FAC
      if(M1 .le. 1 .and. 1 .le. M2) then
c                                             First partials
          do 130 I=1,N
  130     Y1(I+1)=FAC*U1(I+1)
      endif
      if(M2 .eq. 1) return
C                                             Second partials
      K=N+2
          do 150 I=1,N
            FAC2 = U1(1) * Y1(I+1)
               do 140 J=1,I
               Y1(K)= FAC*( U1(K) + FAC2*Y1(J+1) )
               K=K+1
  140          continue
  150     continue
      return
      end
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine DUTAN  (U ,Y )
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  D1, D2, FAC2, U(*),Y(*)
      double precision       ONE, TWO
      parameter(ONE = 1.0D0, TWO = 2.0D0)
C
C     COMPUTE..                        Y=TAN(U)
C
      if(M1 .LE. 0) THEN
         Y(1) = TAN(U(1))
         if(M2 .eq. 0) return
      end if
      D1 = (ONE/COS(U(1)))**2
      D2 = TWO * Y(1) * D1
      if(M1 .le. 1) go to 216
      go to 224
C                                      1ST PARTIALS
  216     do 220 I=1,N
  220     Y(I+1) = D1*U(I+1)
          if(M2 .eq. 1) return
C                                      2ND PARTIALS
  224 K=N+2
         do 228 I=1,N
            FAC2 = D2 * U(I+1)
            do 228 J=1,I
               Y(K)= FAC2 * U(J+1) + D1 * U(K)
  228          K=K+1
      return
      end
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      subroutine DUTANH  (U ,Y )
c     -- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      integer i, j, k
      integer N, M1, M2
      common/UCOM1/N,M1,M2
      save /UCOM1/
      double precision  D1, D2, FAC2, U(*),Y(*)
      double precision       ONE, TWO
      parameter(ONE = 1.0D0, TWO = 2.0D0)
C
C     COMPUTE..                        Y=TANH(U)
C
      if(M1 .LE. 0) THEN
         Y(1) = TANH(U(1))
         if(M2 .eq. 0) return
      end if
      D1 = (ONE/COSH(U(1)))**2
      D2 = -TWO * Y(1) * D1
      if(M1 .le. 1) go to 216
      go to 224
C                                      1ST PARTIALS
  216     do 220 I=1,N
  220     Y(I+1) = D1*U(I+1)
          if(M2 .eq. 1) return
C                                      2ND PARTIALS
  224 K=N+2
         do 228 I=1,N
            FAC2 = D2 * U(I+1)
            do 228 J=1,I
               Y(K)= FAC2 * U(J+1) + D1 * U(K)
  228          K=K+1
      return
      end


subroutine compare_to_sac_file(file, y, n, b, dt)
  implicit none
  character(len=*) file
  real*4 y(*)
  integer n
  real*4 b, dt
  integer nerr
  real*4 y1(10000)
  integer n1
  real*4 dt1, b1
  integer nmax
  integer i

  integer isclosef

  nmax = 10000
  call rsac1(file, y1, n1, b1, dt1, nmax, nerr)
  if (nerr .ne. 0) then
     write(*,*) "Error reading file: ", file(1:len_trim(file)), "nerr: ", nerr
     return
  endif
  if(n .ne. n1) then
     write(*,*) "n: ", n, n1
  endif
  if(abs(b - b1) > 1e-7) then
     write(*,*) "b: ", b, b1
  endif
  if(abs(dt - dt1) > 1e-7) then
     write(*,*) "dt: ", dt, dt1
  endif
  do i = 1, n
     if(isclosef(y(i),y1(i)) .ne. 1) then
        write(*,*) "y[", i, "]: ", y(i),y1(i)
        call exit(-1)
     endif
  enddo
end subroutine compare_to_sac_file

subroutine test_revalresp()
  implicit none
  include "sacf.h"
  real*4 y(10000), b, dt
  integer n, nerr
  integer nmax
  real*8 limits(4)

  integer remove_evalresp_simple

  nmax = 10000
  limits(1) = 0.002
  limits(2) = 0.005
  limits(3) = 12.0
  limits(4) = 20.0
  ! Read in Raw Sac file

  call rsac1("raw.sac", y, n, b, dt, nmax, nerr)
  if(nerr .ne. 0) then
     write(*,*) "Error reading file: raw.sac nerr: ", nerr
     return
  endif

  if(remove_evalresp_simple(y, n, dt, limits) .ne. 0) then
     write(*,*) "Error removing instrument with evalresp"
     return
  endif

  ! write the deconvolved seismogram back to disk
  call wsac0("deconvolved_evr.sac", y, y, nerr)

  ! Check input/output
  call compare_to_sac_file("deconvolved_from_evr.sac",y,n,b,dt)

end subroutine test_revalresp

subroutine test_rpolezero()
  implicit none
  include 'sacf.h'

  real*4 y(10000), b, dt
  integer n, nerr
  integer nmax
  real*8 limits(4)

  integer remove_polezero_simple

  nmax = 10000
  limits(1) = 0.002
  limits(2) = 0.005
  limits(3) = 12.0
  limits(4) = 20.0

  ! Read in Raw Sac file
  call rsac1("raw.sac", y, n, b, dt, nmax, nerr)
  if(nerr .ne. 0) then
     write(*,*) "Error reading file: raw.sac nerr: ", nerr
     return
  endif
  if(remove_polezero_simple(y, n, dt, limits) .ne. 0) then
     write(*,*) "Error removing instrument with polezero"
     return
  endif
  ! write the deconvolved seismogram back to disk
  call wsac0("deconvolved_pz.sac", y, y, nerr)

  ! Check input/output
  call compare_to_sac_file("deconvolved_from_pz.sac", y,n,b,dt)

end subroutine test_rpolezero

program main
  implicit none
  write(*,*) "test evalresp"
  call test_revalresp()
  write(*,*) "test polezero"
  call test_rpolezero()
end program main


text:
program state_main

    use state_module, only: person_t
    implicit none

    type(person_t) :: person
    call person%hungry_state%no_hungry
    call person%work
    call person%work
    call person%eat
    call person%eat

end program state_main

!> Results shall be:

!  Ok, let us do work..
!  I am hungry, no work!!
!  Eatting..
!  Already baole!!

max_stars_repo_name: zoziha/Fortran-Design-Patterns

max_issues_repo_name: zoziha/Fortran-Design-Patterns

max_forks_repo_name: zoziha/Fortran-Design-Patterns


text:
        program main
        integer a,b,c
        external fct
        a = 1
        b = 2
        write (*, *), fct(b)
        c = fct(b) + 100
        write (*,*), a, b, c
        end

        function fct(arg)
        integer arg,fct
        arg = arg * 10
        fct = arg
        return
        end

max_stars_repo_name: dacmot/RevEngTools

max_issues_repo_name: dacmot/RevEngTools

max_forks_repo_name: dacmot/RevEngTools


text:
!-----------------------------------------------------------------------
!     Subroutine:    chisq_bprobes
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          05/26/2012
!     Description:   Calculates fit to B-Probe magnetic diagnostics.
!-----------------------------------------------------------------------
      SUBROUTINE chisq_bprobes(target,sigma,niter,iflag)
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stellopt_runtime
      USE stellopt_targets
      USE diagno_input_mod, ONLY: read_diagno_input
      USE safe_open_mod
!-----------------------------------------------------------------------
!     Input/Output Variables
!
!-----------------------------------------------------------------------
      IMPLICIT NONE
      REAL(rprec), INTENT(in)    ::  target(nprobes)
      REAL(rprec), INTENT(in)    ::  sigma(nprobes)
      INTEGER,     INTENT(in)    ::  niter
      INTEGER,     INTENT(inout) ::  iflag

!-----------------------------------------------------------------------
!     Local Variables
!
!-----------------------------------------------------------------------
      INTEGER :: iunit, dex, itmp, ik,ier
      REAL(rprec), ALLOCATABLE :: xp(:),yp(:),zp(:),modb(:),flux(:)
!----------------------------------------------------------------------
!     BEGIN SUBROUTINE
!----------------------------------------------------------------------
      IF (iflag < 0) RETURN
      IF (niter >= 0) THEN
         ! Read the bprobes file
         iunit = 12
         CALL safe_open(iunit,ier,'diagno_bth.'//TRIM(proc_string_old),'old','formatted')
         IF (ier < 0) THEN; iflag=ier; RETURN; END IF
         dex = 0
         DO
            READ(iunit,*,iostat=ik)
            IF (ik/=0) EXIT
            dex = dex + 1
         END DO
         REWIND(iunit)
         dex = dex - 1 ! Avoid the # line
         ALLOCATE(xp(dex),yp(dex),zp(dex),modb(dex),flux(dex))
         DO ik = 1, dex
            READ(iunit,*) itmp,xp(ik),yp(ik),zp(ik),modb(ik),flux(ik)
         END DO
         CLOSE(iunit)
         IF (iflag == 1) WRITE(iunit_out,'(A,2(2X,I3.3))') 'B_PROBES ',dex,7
         IF (iflag == 1) WRITE(iunit_out,'(A)') 'X  Y  Z  |B|  TARGET  SIGMA  FLUX'
         DO ik = 1, nprobes
            IF (sigma(ik) < bigno) THEN
               mtargets = mtargets + 1
               targets(mtargets) = target(ik)
               sigmas(mtargets)  = sigma(ik)
               vals(mtargets)    = flux(ik)
            END IF
            IF (iflag == 1 .and. ik <= dex) WRITE(iunit_out,'(7ES22.12E3)') xp(ik),yp(ik),zp(ik),modb(ik),target(ik),sigma(ik),flux(ik)
         END DO
         DEALLOCATE(xp,yp,zp,modb,flux)
      ELSE
         DO ik =1, nprobes
            IF (sigma(ik) < bigno) THEN
               lneed_magdiag = .true.
               mtargets = mtargets + 1
               IF (niter == -2) target_dex(mtargets)=jtarget_bprobe
            END IF
         END DO
         CALL read_diagno_input(TRIM(id_string),iflag)
      END IF
      RETURN
!----------------------------------------------------------------------
!     END SUBROUTINE
!----------------------------------------------------------------------
      END SUBROUTINE chisq_bprobes

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


! (from http://patorjk.com/software/taag/ ANSI shadow)

!   ██████╗ ██████╗ ██████╗ ██████╗
!   ██╔══██╗██╔══██╗╚════██╗██╔══██╗
!   ██████╔╝██████╔╝ █████╔╝██║  ██║
!   ██╔═══╝ ██╔══██╗ ╚═══██╗██║  ██║
!   ██║     ██████╔╝██████╔╝██████╔╝
!   ╚═╝     ╚═════╝ ╚═════╝ ╚═════╝

!------------------------------------------------------------------------------!
!>  Peeling Ballooning in 3D
!------------------------------------------------------------------------------!
!>  \author
!!  Toon Weyens,
!!  Contact: weyenst@gmail.com
!------------------------------------------------------------------------------!
!>  \version    2.41
!!  \date       2012-2020
!!  Copyright (C) 2020Toon Weyens
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License as published by
!!  the Free Software Foundation, either version 3 of the License, or
!!  (at your option) any later version.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the License along with this program.
!!  If not, contact weyenst@gmail.com
!------------------------------------------------------------------------------!
!>  \see
!!  References:
!!  \cite weyens2014theory
!!  \cite Weyens2017PB3D
!------------------------------------------------------------------------------!
#define CHCKERR if(ierr.ne.0) then; call sudden_stop(ierr); end if
program PB3D
    use num_vars, only: prog_name, prog_style, rank, rich_restart_lvl, eq_style
    use str_utilities, only: r2str, i2str
    use messages
    use grid_vars, only: grid_type
    use eq_vars, only: eq_1_type, eq_2_type
    use X_vars, only: X_1_type, X_2_type
    use vac_vars, only: vac_type
    use sol_vars, only: sol_type
    use HDF5_vars, only: init_HDF5
    use driver_eq, only: run_driver_eq
    use driver_X, only: run_driver_X
    use driver_sol, only: run_driver_sol
    use files_ops, only: open_input, open_output, parse_args, init_files, &
        &close_output
    use input_ops, only: read_input_opts, read_input_eq, print_output_in
    use input_utilities, only: dealloc_in
    use MPI_ops, only: start_MPI, stop_MPI, broadcast_input_opts, sudden_stop
    use eq_ops, only: calc_normalization_const, normalize_input
    use eq_utilities, only: do_eq, eq_info
    use rich_ops, only: init_rich, term_rich, do_rich, start_rich_lvl, &
        &stop_rich_lvl
    use rich_vars, only: rich_info, &
        &rich_lvl
#if ldebug
    use num_vars, only: ltest
    use test, only: generic_tests
#endif

    implicit none

    ! local variables
    integer :: ierr                                                             ! error
    type(grid_type), target :: grid_eq                                          ! equilibrium grid
    type(grid_type), pointer :: grid_eq_B => null()                             ! field-aligned equilibrium grid
    type(grid_type), target :: grid_X                                           ! perturbation grid
    type(grid_type), pointer :: grid_X_B => null()                              ! field-aligned perturbation grid
    type(grid_type) :: grid_sol                                                 ! solution grid
    type(eq_1_type) :: eq_1                                                     ! flux equilibrium variables
    type(eq_2_type) :: eq_2                                                     ! metric equilibrium variables
    type(X_1_type) :: X_1                                                       ! vectorial perturbation variables
    type(X_2_type) :: X_2                                                       ! integrated tensorial perturbation variables
    type(vac_type) :: vac                                                       ! vacuum variables
    type(sol_type) :: sol                                                       ! solution variables

    !------------------------------!
    !   Initialize some routines   !
    !------------------------------!
    ierr = start_MPI()                                                          ! start MPI
    CHCKERR
    prog_name = 'PB3D'                                                          ! program name
    prog_style = 1                                                              ! main part
    call print_hello()                                                          ! print message with time, etc
    call init_output()                                                          ! initialize output utilities
    call init_files()                                                           ! initialize file operations
    call init_time()                                                            ! initialize time
    call init_HDF5()                                                            ! initialize HDF5

    !-----------------------------------------------------------------------!
    !   Read the user-provided input file and the equilibrium code output   !
    !-----------------------------------------------------------------------!
    call start_time
    call writo('Initialization')
    call lvl_ud(1)
    if (rank.eq.0) then                                                         ! only master
        ierr = parse_args()                                                     ! parse argument (options are used in open_input)
        CHCKERR
        ierr = open_input()                                                     ! open the input files
        CHCKERR
        ierr = read_input_opts()                                                ! read input options file
        CHCKERR
        if (rich_restart_lvl.eq.1) then                                         ! start at first level
            ierr = read_input_eq()                                              ! read input equilibrium file
            CHCKERR
            ierr = calc_normalization_const()                                   ! set up normalization constants
            CHCKERR
            call normalize_input()                                              ! normalize the input
            ierr = open_output()                                                ! open output file
            CHCKERR
            ierr = print_output_in('in')                                        ! print input outputs
            CHCKERR
            call dealloc_in()                                                   ! clean up input from equilibrium codes
        else                                                                    ! restart at higher level
            ierr = open_output()                                                ! open output file
            CHCKERR
        end if
    end if
    ierr = broadcast_input_opts()                                               ! broadcast input options to other processors
    CHCKERR
    call writo('')
    call passed_time
    call writo('')
    call lvl_ud(-1)

#if ldebug
    !-------------------!
    !   Do some tests   !
    !-------------------!
    if (ltest) then
        call start_time
        call writo('Generic Tests')
        call lvl_ud(1)
        ierr = generic_tests()                                                  ! run generic test
        CHCKERR
        call writo('')
        call passed_time
        call writo('')
        call lvl_ud(-1)
    end if
#endif

    !----------------------------------------------!
    !   Initialize Richardson Extrapolation Loop   !
    !----------------------------------------------!
    call start_time
    ierr = init_rich()
    CHCKERR
    call stop_time

    RICH: do while(do_rich())
        !----------------------------!
        !   Start Richardson level   !
        !----------------------------!
        call start_time
        call writo('Richardson level '//trim(i2str(rich_lvl)))
        call lvl_ud(1)
        ierr = start_rich_lvl()                                                 ! start Richardson level, setting n_r_sol and other variables
        CHCKERR
        call stop_time
        call writo('')
        call lvl_ud(-1)

        PAR: do while(do_eq())
            !-----------------------------------!
            !   Main Driver: Equilibrium part   !
            !-----------------------------------!
            call start_time
            call writo('Equilibrium driver'//trim(rich_info())//&
                &trim(eq_info()))
            call lvl_ud(1)
            ierr = run_driver_eq(grid_eq,grid_eq_B,eq_1,eq_2,vac)               ! equilibrium driver
            CHCKERR
            call writo('')
            call passed_time
            call writo('')
            call lvl_ud(-1)

            !------------------------------------!
            !   Main driver: Perturbation part   !
            !------------------------------------!
            call start_time
            call writo('Perturbation driver'//trim(rich_info())//&
                &trim(eq_info()))
            call lvl_ud(1)
            ierr = run_driver_X(grid_eq,grid_eq_B,grid_X,grid_X_B,eq_1,eq_2,&
                &X_1,X_2)                                                       ! perturbation driver
            CHCKERR
            call writo('')
            call passed_time
            call writo('')
            call lvl_ud(-1)
        end do PAR

        !--------------------------------!
        !   Main driver: Solution part   !
        !--------------------------------!
        call start_time
        call writo('Solution driver'//trim(rich_info()))
        call lvl_ud(1)
        ierr = run_driver_sol(grid_eq,grid_X,grid_sol,X_2,vac,sol)              ! solution driver
        CHCKERR
        call writo('')
        call passed_time
        call writo('')
        call lvl_ud(-1)

        call start_time
        call stop_rich_lvl()
        call stop_time
    end do RICH

    !---------------------------------------!
    !   Stop Richarson Extrapolation Loop   !
    !---------------------------------------!
    call start_time
    call term_rich()                                                            ! stop Richardson loop
    call stop_time

    !--------------!
    !   Clean up   !
    !--------------!
    call writo('Clean up')
    call lvl_ud(1)
    if (eq_style.eq.2) call eq_2%dealloc()
    ierr = stop_MPI(grid_eq=grid_eq,&
        &grid_eq_B=grid_eq_B,&
        &grid_X=grid_X,&
        &grid_X_B=grid_X_B,&
        &grid_sol=grid_sol,&
        &eq_1=eq_1,&
        &eq_2=eq_2,&
        &X_1=X_1,&
        &X_2=X_2,&
        &vac=vac,&
        &sol=sol)
    CHCKERR
    call close_output
    call lvl_ud(-1)

    call print_goodbye
end program PB3D


text:
program kaprekar
    use, intrinsic :: iso_fortran_env, only : error_unit
    implicit none
    integer, parameter :: nr_digits = 4
    integer, dimension(:), allocatable :: digits
    integer :: number, prev_number, small, large, istat

    number = get_argument()
    allocate(digits(nr_digits), stat=istat)
    if (istat /= 0) then
        write (error_unit, '(A)') 'Error: can not allocate array'
        stop 5
    end if
    prev_number = 0
    do while (number /= prev_number)
        call number2digits(number, digits)
        call sort_digits(digits)
        small = digits2number(digits)
        large = reverse_digits2number(digits)
        prev_number = number
        number = large - small
    end do
    print '(I0)', number

contains

    integer function get_argument() result(number)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        character(len=80) :: buffer
        integer :: istat, first_digit, n, i
        logical :: okay

        if (command_argument_count() /= 1) then
            write (error_unit, '(A)') &
                'Error: expecting a single number as argument'
            stop 1
        end if
        call get_command_argument(1, buffer)
        read (buffer, '(I10)', iostat=istat) number
        if (istat /= 0) then
            write (error_unit, '(A)') &
                'Error: expecting as integer as argument'
            stop 2
        end if
        if (number < 1 .or.  number > 9999) then
            write (error_unit, '(A)') &
                'Error: argument should be between 1 and 9999'
            stop 3
        endif
        okay = .false.
        first_digit = mod(number, 10)
        n = number/10
        do i = 2, 4
            if (first_digit /= mod(n, 10)) then
                okay = .true.
            exit
            end if
            n = n/10
        end do
        if (.not.  okay) then
            write (error_unit, '(A)') &
                'Error: number must have at least two distinct digits'
            stop 4
        end if

    end function get_argument

    subroutine number2digits(number, digits)
        implicit none
        integer, value :: number
        integer, dimension(:), intent(inout) :: digits
        integer :: i

        do i = 1, size(digits)
            digits(i) = mod(number, 10)
            number = number/10
        end do

    end subroutine number2digits

    integer function digits2number(digits) result(number)
        implicit none
        integer, dimension(:), intent(in) :: digits
        integer :: factor, i

        number = 0
        factor = 1
        do i = 1, size(digits)
            number = number + factor*digits(i)
            factor = factor*10
        end do

    end function digits2number

    integer function reverse_digits2number(digits) result(number)
        implicit none
        integer, dimension(:), intent(in) :: digits
        integer :: factor, i

        number = 0
        factor = 1
        do i = size(digits), 1, -1
            number = number + factor*digits(i)
            factor = factor*10
        end do

    end function reverse_digits2number

    subroutine sort_digits(digits)
        implicit none
        integer, dimension(:), intent(inout) :: digits
        integer :: i, j, index, tmp

        do i = 1, size(digits)
            index = i
            do j = i + 1, size(digits)
                if (digits(index) < digits(j)) &
                    index = j
            end do
            tmp = digits(index)
            digits(index) = digits(i)
            digits(i) = tmp
        end do

    end subroutine sort_digits

end program kaprekar

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC