text:
module linalg
    implicit none

    integer,parameter   :: dp = kind(1.d0)

    interface GAssembly
        procedure   GAssembly1
        procedure   GAssembly2
        procedure   GAssembly3
    end interface

    contains

    !-------- Matrix inversion --------!
    function inv(A) result(Ainv)
        implicit none

        complex(dp),intent(in)       :: A(:,:)
        complex(dp)                  :: Ainv(size(A,1),size(A,2))

        integer                      :: s(2), nrows, ncols
        integer                      :: lda, lwork, info

        integer                      :: ipiv(size(A,1))
        complex(dp)                  :: work(size(A,1))

        external zgetrf
        external zgetri

        s = shape(A)
        nrows = s(1)
        ncols = s(2)
        lda = max(1,nrows)
        lwork = max(1,nrows)

        Ainv = A
        call zgetrf(nrows,ncols,Ainv,lda,ipiv,info)

        if (info .ne. 0) then
            print *, "Error in the LU factorization!"
        end if

        call zgetri(nrows,Ainv,lda,ipiv,work,lwork,info)

        if (info .ne. 0) then
            print *,"Error in the matrix inversion!"
        end if

    end function

    !------- Solve Linear System -------!
    function solve(A,B) result(X)
        implicit none

        complex(dp),intent(in)   :: A(:,:), B(:,:)
        complex(dp)              :: X(size(B,1),size(B,2))
        complex(dp)              :: Y(size(A,1),size(A,2))

        integer                  :: s(2), nra, nrb, nca, ncb
        integer                  :: lda, ldb, info

        integer                  :: ipiv(size(A,1))

        external zgesv

        s = shape(A)
        nra = s(1)
        nca = s(2)
        s = shape(B)
        nrb = s(1)
        ncb = s(2)
        lda = max(1,nra)
        ldb = max(1,nrb)

        X = B
        Y = A
        call zgesv(nra,ncb,Y,lda,ipiv,X,ldb,info)

        if (info .ne. 0) then
            print *,"Error in the linear system solver!"
        end if
    end function

    !---------------------------------------------------!
    ! Triple Matrix Product                             !
    !                                                   !
    ! D = A.B.C                                         !
    !                                                   !
    ! A: Complex                                        !
    ! B: Real                                           !
    ! C: Complex                                        !
    !---------------------------------------------------!
    function tri(A,B,C) result(D)
        implicit none

        complex(dp),intent(in)     :: A(:,:)
        complex(dp),intent(in)     :: B(:,:)
        complex(dp),intent(in)     :: C(:,:)
        complex(dp)                :: D(size(A,1),size(C,2))

        ! Hu & Shing matrix chain multiplication
        ! [m] ---B--- [k]
        !  |           |
        !  A           C
        !  |           |
        ! [n]---ABC-- [l]

        integer :: n,m,k,l
        integer :: s(2)

        s = shape(A)
        n = s(1)
        m = s(2)
        s = shape(C)
        k = s(1)
        l = s(2)

        if (n*m*k + n*k*l < m*k*l + n*m*l) then
            D = matmul(matmul(A,B),C)
        else
            D = matmul(A,matmul(B,C))
        end if
    end function

    !---------------------------------------------------!
    ! For a diagonal matrix B:                          !
    ! C = B.A.B+ if invert = false                      !
    ! C = B+.A.B if invert = true                       !
    !---------------------------------------------------!
    function d2tri(A,B,invert) result(C)
        use omp_lib
        implicit none

        complex(dp),intent(in)  :: A(:,:)
        complex(dp),intent(in)  :: B(:,:)
        logical,intent(in)      :: invert
        complex(dp)             :: C(size(A,1),size(A,2))

        integer     :: i,j
        complex(dp) :: p

        !$OMP PARALLEL DO
        do j = 1,size(B,1)
            p = B(j,j)
            if (invert .eqv. .false.) then
                p = conjg(p)
            end if
            do i = 1,size(A,1)
                C(i,j) = p*A(i,j)
                if (invert .eqv. .false.) then
                    C(i,j) = C(i,j)*B(i,i)
                else
                    C(i,j) = C(i,j)*conjg(B(i,i))
                end if
            end do
        end do
        !$OMP END PARALLEL DO
    end function

    !---------------------------------------------------!
    ! G = [(E+i.eta).I - H - S1 - S2]^(-1)                    !
    !---------------------------------------------------!
    function GAssembly3(E,H,S1,S2,eta) result(G)
        implicit none

        real,intent(in)         :: E
        complex(dp),intent(in)  :: H(:,:)
        complex(dp),intent(in)  :: S1(:,:)
        complex(dp),intent(in)  :: S2(:,:)
        real(dp),intent(in)     :: eta
        complex(dp)             :: G(size(H,1),size(H,2))

        integer             :: i

        ! G = (E+i.eta).I - H - S1 - S2
        G = -(H + S1 + S2)
        do i = 1,size(H,2)
            G(i,i) = G(i,i) + E + cmplx(0,1)*eta
        end do

        G = inv(G)
    end function

    !---------------------------------------------------!
    ! G = [(E+i.eta).I - H - S]^(-1)                    !
    !---------------------------------------------------!
    function GAssembly2(E,H,S,eta) result(G)
        implicit none

        real,intent(in)         :: E
        complex(dp),intent(in)  :: H(:,:)
        complex(dp),intent(in)  :: S(:,:)
        real(dp),intent(in)     :: eta
        complex(dp)             :: G(size(H,1),size(H,2))

        integer             :: i

        ! G = (E+i.eta).I - H - S
        G = -(H + S)
        do i = 1,size(H,2)
            G(i,i) = G(i,i) + E + cmplx(0,1)*eta
        end do

        G = inv(G)
    end function

    !---------------------------------------------------!
    ! G = [(E+i.eta).I - H]^(-1)                        !
    !---------------------------------------------------!
    function GAssembly1(E,H,eta) result(G)
        implicit none

        real,intent(in)         :: E
        complex(dp),intent(in)  :: H(:,:)
        real(dp),intent(in)     :: eta
        complex(dp)             :: G(size(H,1),size(H,2))

        integer             :: i

        ! G = (E+i.eta).I - H
        G = -H
        do i = 1,size(H,2)
            G(i,i) = G(i,i) + E + cmplx(0,1)*eta
        end do

        G = inv(G)
    end function

end module

max_stars_repo_name: StxGuy/GreenCheetah

max_issues_repo_name: StxGuy/GreenCheetah

max_forks_repo_name: StxGuy/GreenCheetah


text:
C$Procedure ADDCOM ( Add comments to DSK file )

      SUBROUTINE ADDCOM ( HANDLE, CMDFIL, INPFN, OUTFN, CMTFIL, APPFLG )
      IMPLICIT NONE

C$ Abstract
C
C     Add comments to the output DSK file created by MKDSK.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     DSK
C
C$ Keywords
C
C     FILES
C     TOPOGRAPHY
C
C$ Declarations

      INCLUDE 'mkdsk.inc'

      CHARACTER*(*)         CMDFIL
      INTEGER               HANDLE
      CHARACTER*(*)         CMTFIL
      CHARACTER*(*)         INPFN
      CHARACTER*(*)         OUTFN
      LOGICAL               APPFLG

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     CMDFIL     I   Name of setup file.
C     HANDLE     I   Handle of DSK file.
C     CMTFIL     I   Name of comment file.
C     INPFN      I   Input shape file name.
C     OUTFN      I   Output DSK file name.
C     APPFLG     I   Append flag.
C
C$ Detailed_Input
C
C     CMDFIL         is the name of the MKDSK setup file.
C
C     HANDLE         is the handle associated with the output
C                    DSK file.
C
C     CMTFIL         is the name of a text file containing comments
C                    to be added to the DSK file.  If CMTFIL is
C                    blank, this is interpreted to mean there is no
C                    comment file.
C
C     INPFN          is the name of the input shape file to be
C                    converted to DSK format.
C
C     OUTFN          is the name of the output DSK file resulting
C                    from conversion of the shape file.
C
C     APPFLG         is a logical flag which is .TRUE. if the output
C                    data are to be appended to an existing DSK file
C                    and .FALSE. if the output DSK file is new.
C
C$ Detailed_Output
C
C     None.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If an error occurs while reading the comment file, the
C        error will be diagnosed by routines in the call tree of
C        this routine.
C
C     2) If an error occurs while attempting to add comments to the
C        output DSK file, the error will be diagnosed by routines in
C        the call tree of this routine.
C
C$ Files
C
C     Normally, the mapping implemented by this routine is defined
C     by a kernel variable introduced into the kernel pool by loading
C     a SPICE text kernel.
C
C$ Particulars
C
C     This routine adds comments to the comment area of the output
C     DSK file. These are:
C
C        - The contents of a comment file, if any, specified in
C          the setup file.
C
C        - The run time and date, the names of the setup, input,
C          and output files, and an indication of whether the
C          output file was new or appended to.
C
C        - The contents of the setup file.
C
C$ Examples
C
C     None.
C
C$ Restrictions
C
C     1) This routine is intended for use only within the MKDSK
C        program.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N.J. Bachman    (JPL)
C
C$ Version
C
C-    MKDSK Version 2.0.0, 07-FEB-2017 (NJB)
C
C        Now writes the MKDSK version string to the comment area
C        of the output file.
C
C-    MKDSK Version 1.0.0, 15-APR-2010 (NJB)
C
C-&

C
C     SPICELIB functions
C
      INTEGER               FRSTNP
      INTEGER               RTRIM

      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local parameters
C
      INTEGER               CLINSZ
      PARAMETER           ( CLINSZ = 80 )

      INTEGER               MXCMNT
      PARAMETER           ( MXCMNT = 10000 )

C
C     Local variables
C

      CHARACTER*(SHRTLN)    ASTRLN
      CHARACTER*(FILSIZ)    CMNBUF ( LBCELL : MXCMNT )
      CHARACTER*(CLINSZ)    LINE
      CHARACTER*(SHRTLN)    TSTAMP

      DOUBLE PRECISION      TVEC   ( 6 )

      INTEGER               CMNUNT
      INTEGER               L
      INTEGER               M

      LOGICAL               EOF

C
C     Saved variables
C

C
C     Save large buffers to avoid stack problems in
C     some C environments.
C
      SAVE

      IF ( RETURN() ) THEN
         RETURN
      END IF

      CALL CHKIN ( 'ADDCOM' )

C
C     Set the maximum size of the comment line buffer.
C
      CALL SSIZEC ( MXCMNT, CMNBUF )
      L = 0

C
C     Comment area separator line.
C
      ASTRLN = '****************************************' //
     .         '****************************************'

C
C     If the comment file was provided we copy its content to the
C     comment area.
C
      IF ( CMTFIL .NE. ' ' ) THEN
C
C        We open the comment file, copy text from it to the comment
C        buffer line by line, clean non-printing characters from the
C        lines on the fly and dump the buffer to the comment area
C        when it's full. We repeat until all comments have been copied.
C
         CALL TXTOPR ( CMTFIL, CMNUNT )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'ADDCOM' )
            RETURN
         END IF

C
C        Insert top comment separator line.
C
         CMNBUF ( 1 ) = ASTRLN
         CMNBUF ( 2 ) = ' '
         L            = 2

C
C        Get next comment line.
C
         CALL READLN ( CMNUNT, LINE, EOF )

         DO WHILE (  .NOT.  ( EOF .OR. FAILED() )  )
C
C           Replace non-printing characters with spaces.
C
            DO WHILE ( FRSTNP ( LINE ) .NE. 0 )
               M          = FRSTNP( LINE )
               LINE (M:M) = ' '
            END DO

            IF ( L .LT. MXCMNT ) THEN
C
C              Store line in the buffer.
C
               L            = L + 1
               CMNBUF ( L ) = LINE ( : RTRIM(LINE) )

            ELSE
C
C              Buffer is full. Set the cardinality of the comment
C              buffer and write it to DSK comment area. Reset
C              counter.
C
               CALL SCARDC ( L,         CMNBUF    )
               CALL DASAC  ( HANDLE, L, CMNBUF(1) )

               L = 0

            END IF

C
C           Get next comment line.
C
            CALL READLN ( CMNUNT, LINE, EOF )

         END DO

C
C        Dump the rest of the buffer into the comment area.
C
         IF ( L .NE. 0 ) THEN

            CALL SCARDC ( L,         CMNBUF    )
            CALL DASAC  ( HANDLE, L, CMNBUF(1) )

            L = 0

         END IF

C
C        Close comment file.
C
         CLOSE ( CMNUNT )

      END IF

C
C     Add a header preceding contents of the setup file and containing
C     setup file name and current CPU time.
C
      CALL CPUTIM ( TVEC )
      TSTAMP = 'YYYY-MM-DDTHR:MN:SC'

      CALL DPFMT ( TVEC(1), '0YYY',  TSTAMP(1:4)   )
      CALL DPFMT ( TVEC(2), '0M',    TSTAMP(6:7)   )
      CALL DPFMT ( TVEC(3), '0D',    TSTAMP(9:10)  )
      CALL DPFMT ( TVEC(4), '0h',    TSTAMP(12:13) )
      CALL DPFMT ( TVEC(5), '0m',    TSTAMP(15:16) )
      CALL DPFMT ( TVEC(6), '0s',    TSTAMP(18:19) )

      CMNBUF ( 1 ) = ' '
      CMNBUF ( 2 ) =  ASTRLN
      CMNBUF ( 3 ) = 'MKDSK VERSION:       ' // VER
      CMNBUF ( 4 ) = 'MKDSK RUN DATE/TIME: ' // TSTAMP(:RTRIM(TSTAMP))
      CMNBUF ( 5 ) = 'MKDSK SETUP FILE:    ' // CMDFIL(:RTRIM(CMDFIL))
      CMNBUF ( 6 ) = 'MKDSK INPUT FILE:    ' // INPFN (:RTRIM(INPFN))
      CMNBUF ( 7 ) = 'MKDSK OUTPUT FILE:   ' // OUTFN (:RTRIM(OUTFN))

      IF ( APPFLG ) THEN
         CMNBUF ( 8 ) = 'OUTPUT FILE STATUS:    EXISTING FILE'
      ELSE
         CMNBUF ( 8 ) = 'OUTPUT FILE STATUS:    NEW FILE'
      END IF

      CMNBUF ( 9  ) = ASTRLN
      CMNBUF ( 10 ) = ' '

      L             = 10

C
C     Now we will copy contents of the setup file to the comment area
C     using exactly the same procedure: open the setup file, copy
C     text from the file into the buffer line by line, clean
C     non-printing characters from the lines on the fly and dump the
C     buffer to the comment area when it's full. We repeat until all
C     setup lines have been copied.
C
      CALL TXTOPR ( CMDFIL, CMNUNT )
      EOF = .FALSE.

      DO WHILE ( .NOT. EOF )
C
C        Read next line.
C
         CALL READLN ( CMNUNT, LINE, EOF )

         IF ( .NOT. EOF ) THEN
C
C           Replace non-printing character with spaces.
C
            DO WHILE ( FRSTNP ( LINE ) .NE. 0 )
               M            = FRSTNP( LINE )
               LINE ( M:M ) = ' '
            END DO

            IF ( L .LT. MXCMNT ) THEN
C
C              Store line on buffer.
C
               L           = L + 1
               CMNBUF( L ) = LINE ( : RTRIM(LINE) )

            ELSE
C
C              Buffer is full. Set the cardinality of the comment
C              buffer and write it to DSK comment area. Reset counter
C              and store the last line that we have obtained in the
C              first line of the buffer.
C
               CALL SCARDC ( L,         CMNBUF    )
               CALL DASAC  ( HANDLE, L, CMNBUF(1) )

               L = 0

            END IF

         END IF

      END DO

      CLOSE ( CMNUNT )

C
C     Add "bottom of the comments" separator line.
C
      IF ( L .LE. MXCMNT - 2 ) THEN

         CMNBUF ( L + 1 ) = ' '
         CMNBUF ( L + 2 ) = ASTRLN
         L = L + 2

      ELSE

C
C        Dump current contents of the comment buffer, first. After
C        that stick separator at the top of the buffer.
C
         CALL SCARDC ( L, CMNBUF )
         CALL DASAC  ( HANDLE, L, CMNBUF(1) )

         CMNBUF ( 1 ) = ' '
         CMNBUF ( 2 ) = ASTRLN
         L            = 2

      END IF

C
C     Dump the buffer one more time, if required.
C
      IF ( L .GT. 0 ) THEN

         CALL SCARDC ( L,         CMNBUF    )
         CALL DASAC  ( HANDLE, L, CMNBUF(1) )

      END IF

      CALL CHKOUT ( 'ADDCOM' )
      RETURN
      END

ext: f

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
! { dg-do compile }
!
! PR 44565: [4.6 Regression] [OOP] ICE in gimplify_expr with array-valued generic TBP
!
! Contributed by Hans-Werner Boschmann <boschmann@tp1.physik.uni-siegen.de>

module ice6

  type :: t
   contains
     procedure :: get_array
     generic :: get_something => get_array
  end type

contains

  function get_array(this)
    class(t) :: this
    real,dimension(2) :: get_array
  end function get_array

  subroutine do_something(this)
    class(t) :: this
    print *,this%get_something()
  end subroutine do_something

end module ice6

! { dg-final { cleanup-modules "ice6" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


! { dg-do compile }
! PR fortran/65173
program min_obj
   implicit none
   integer, parameter :: a = 128
   type :: param_t
      integer :: n= 0
      real*8, dimension(256), allocatable :: x           ! { dg-error "must have a deferred shape" }
      real*8, dimension(2,256), allocatable :: bounds    ! { dg-error "must have a deferred shape" }
      character(a), dimension(256), allocatable :: names ! { dg-error "must have a deferred shape" }
   end type param_t
   contains
      subroutine extrace_params_from_section ( )
         character(*), dimension(), parameter :: &      ! { dg-error "expression in array specification" }
         & char_params = ['element', 'parametrization']
      end subroutine extrace_params_from_section
end program min_obj


!--------------------------------------------------
!..Homework 2
!..Team 31
!..AE305 - Numerical Methods
!
!..This program calculates the chamber pressure,
!..burn rate and specific impulse for given rocket
!..using RK4 method with constant time step and
!..outputs results to a file.
!--------------------------------------------------
Module rocket_params_vars
  Implicit none
  integer, parameter :: no_eqs = 2
  real*8, parameter :: p_a = 101325.d0, &
                     rho_p = 1140.d0, &
                         n = 0.305d0, &
                         a = 0.0000555d0, &
                       r_0 = 0.05d0, &
                       r_f = 0.15d0, &
                         L = 1.25d0, &
                       T_c = 2810.d0, &
                     R_cst = 365.d0, &
                      gama = 1.25d0, &
                      grav = 9.81d0, &
                        pi = 4.*atan(1.d0)
  real*8 :: p_c, p_c_dot, &
          r, r_dot, &
          p_c_old, r_old, &
          A_star

  contains

  Function m_n_dot(p_c_)
    ! calculate m_n dot to use in ODE
    real*8 :: m_n_dot, p_c_
    m_n_dot = p_c_*A_star*sqrt(gama/(R_cst*T_c))*((gama+1)/2)**(-(gama+1)/(2*(gama-1)))
    return
  End Function m_n_dot

  Function rho_c(p_c_)
    ! calculate rho_c to use in ODE
    real*8 :: rho_c, p_c_
    rho_c = p_c_/(R_cst*T_c)
    return
  End Function rho_c

  Function f_cor(rad)
   ! perimeter factor
   Implicit none
    real*8 :: f_cor, eta, rad
     eta = (r_f-rad)/(r_f-r_0)
     f_cor = 1.d0
     if( eta < 0 )then
         f_cor = 0
     else if (eta <= 0.15)then
         f_cor = 1-exp(-7*eta)
     end if
   return
  End Function f_cor

  Function I_sp(p_c_)
    ! specific impulse
    Implicit none
    real*8 :: I_sp, p_c_
    ! Here, absolute of value inside square root is taken because
    ! there is really small error caused by numerical computing makes
    ! value inside of the square root negative for some intervals
    ! which results in runtime error. This negative value is so small
    ! that its practically zero. Therefore, absolute can be used without
    ! effecting accuracy of result. Using abs solved the runtime error.
      I_sp = (1/grav)*sqrt(abs(((2*gama*R_cst*T_c)/(gama-1))*(1-((p_a/p_c_)**((gama-1)/gama)))))
    return
  End function I_sp

  Function ODEs( p_c_ ) result ( k )
    Implicit none
    real*8 :: p_c_
    real*8, dimension( no_eqs ) :: k
    r_dot = a*p_c**n
    k( 1 ) =  r_dot
    k( 2 ) =  R_cst*T_c*(((f_cor(r)*2* a * p_c ** n)/r)*(rho_p-rho_c(p_c_))-m_n_dot(p_c_)/(pi*L*r**2))
    return
  End function

  Subroutine RK4( t, del_t )
    Implicit none
    integer, parameter :: no_stages = 4
    real*8 :: t, del_t
    real*8, dimension( no_eqs )  ::  k( no_stages, no_eqs ), phi(no_eqs)
    p_c_old = p_c
    r_old = r
    k(1, : ) = ODEs( p_c ) !Slopes at the beginning of the interval are calculated

  !Slopes at p_1=0.5 fraction of the interval are calculated
    p_c = p_c_old + 0.5d0 * del_t * k(1, 2)
    r = r_old + 0.5d0 * del_t * k(1, 1)
    k(2, : ) = ODEs( p_c )

  !Slopes at p_2=0.5 fraction of the interval are calculated
    p_c = p_c_old + 0.5d0 * del_t * k(2, 2)
    r = r_old + 0.5d0 * del_t * k(2, 1)
    k(3, : ) = ODEs( p_c )

  !Slopes at p_3=1.0 fraction of the interval are calculated
    p_c = p_c_old + del_t * k(3, 2)
    r = r_old + del_t * k(3, 1)
    k(4, : ) = ODEs( p_c )

    phi(:) = (1.0/6)* (k(1, : )+2*k(2, : )+2*k(3, : )+k(4, : )) !weighted slopes are calculated
    p_c = p_c_old + phi(2) * del_t ! p_c is updated
    r = r_old + phi(1) * del_t ! r is updated

    t = t + del_t ! time is updated
  return
  End subroutine

End module

Program Rocket_Perf
  Use rocket_params_vars
  Implicit none
  character*40 :: fname
  real*8 ::  dt, time, isp
  integer :: nstep = 0
  real*8 :: th_radius

  write(*,'(a)') 'Enter throat radius in centimeters :'
  write(*,'(a)',advance='no') ':>'
  read*, th_radius

  A_star = pi * (th_radius/100)**2

  print*,'enter time step size, dt [s] : '
  read*, dt

!.. initial params
  p_c =  p_a    ! chamber pressure
  r = r_0      ! chamber radius
  isp = I_sp(p_c)

  write(*,'(a)',advance='no')' Enter the output file name [rocket.dat]:>'
  read(*,'(a)') fname
  if( fname .eq. ' ') fname = 'rocket.dat'
  open(1 ,file=fname, form='formatted')

  write(1,*)' "t [s]" "p_c [MPa]" "I_sp[s]" "r_dot[m/s]"'

  time = 0.d0
  nstep = 0

  do while( p_c >= p_a )
    nstep = nstep + 1

    ! update solution and time
    call RK4(time, dt)

    ! calculate isp
    isp = I_sp(p_c)

    ! print on screen and store soln
    if( nstep == 1 .or. mod(nstep,5)==0) &
      write(1,'(8(4x,e12.6))') time, p_c*1.d-6 ,isp, r_dot
      write(*,'(8(a,e9.3))')' t = ',time,' s,   p_c = ', p_c*1d-6,' MPa, I_sp = ', isp, ' s r_dot = ',r_dot,' m/s'

  enddo

  ! Print the interval number for bonus part
  write(*,*) 'Number of time intervals needed = ', nstep

  close(1)

  stop
End


text:
module mpi_f08

  use mpi, only: mpi_status_size

  implicit none
  save

  ! TODO: Make everything private by default

  private mpi_status_size

  logical, parameter :: mpi_async_protects_nonblocking = .false.
  logical, parameter :: mpi_subarrays_supported = .false.

  type :: mpi_comm
     integer :: comm
  end type mpi_comm

  type :: mpi_datatype
     integer :: datatype
  end type mpi_datatype

  type :: mpi_status
     integer :: status(mpi_status_size)
     integer :: mpi_source
     integer :: mpi_tag
     integer :: mpi_error
     integer :: padding
  end type mpi_status

  type(mpi_comm) :: mpi_comm_null
  type(mpi_comm) :: mpi_comm_self
  type(mpi_comm) :: mpi_comm_world

  type(mpi_datatype) :: mpi_int
  type(mpi_datatype) :: mpi_integer

  interface operator (.eq.)
     module procedure mpi_comm_eq
     module procedure mpi_datatype_eq
  end interface operator (.eq.)

  interface operator (.ne.)
     module procedure mpi_comm_ne
     module procedure mpi_datatype_ne
  end interface operator (.ne.)

  interface mpi_send
     module procedure mpi_send_impl
  end interface mpi_send

  interface mpi_recv
     module procedure mpi_recv_impl
  end interface mpi_recv

  interface mpi_get_count
     module procedure mpi_get_count_impl
  end interface mpi_get_count

  interface mpi_comm_size
     module procedure mpi_comm_size_impl
  end interface mpi_comm_size

  interface mpi_comm_rank
     module procedure mpi_comm_rank_impl
  end interface mpi_comm_rank

  interface mpi_init
     module procedure mpi_init_impl
  end interface mpi_init

  interface mpi_finalize
     module procedure mpi_finalize_impl
  end interface mpi_finalize

  interface mpi_abort
     module procedure mpi_abort_impl
  end interface mpi_abort

  interface mpi_wtime
     double precision function mpi_wtime()
       implicit none
     end function mpi_wtime
  end interface mpi_wtime

  interface mpi_wtick
     double precision function mpi_wtick()
       implicit none
     end function mpi_wtick
  end interface mpi_wtick

contains

  logical function mpi_comm_eq(comm1, comm2)
    type(mpi_comm), intent(in) :: comm1, comm2
    mpi_comm_eq = comm1%comm .eq. comm2%comm
  end function mpi_comm_eq

  logical function mpi_comm_ne(comm1, comm2)
    type(mpi_comm), intent(in) :: comm1, comm2
    mpi_comm_ne = comm1%comm .ne. comm2%comm
  end function mpi_comm_ne

  logical function mpi_datatype_eq(datatype1, datatype2)
    type(mpi_datatype), intent(in) :: datatype1, datatype2
    mpi_datatype_eq = datatype1%datatype .eq. datatype2%datatype
  end function mpi_datatype_eq

  logical function mpi_datatype_ne(datatype1, datatype2)
    type(mpi_datatype), intent(in) :: datatype1, datatype2
    mpi_datatype_ne = datatype1%datatype .ne. datatype2%datatype
  end function mpi_datatype_ne

  subroutine mpi_send_impl(buf, count, datatype, dest, tag, comm, ierror)
    use mpi, only: mpi_send
    !gcc$ attributes no_arg_check :: buf
    type(*), intent(in) :: buf
    integer, intent(in) :: count
    type(mpi_datatype), intent(in) :: datatype
    integer, intent(in) :: dest
    integer, intent(in) :: tag
    type(mpi_comm), intent(in) :: comm
    integer, optional, intent(out) :: ierror

    integer ierror1

    call mpi_send(buf, count, datatype%datatype, dest, tag, comm%comm, ierror1)

    if (present(ierror)) ierror = ierror1
  end subroutine mpi_send_impl

  subroutine mpi_recv_impl(buf, count, datatype, dest, tag, comm, status, ierror)
    use mpi, only: mpi_status_size, mpi_source, mpi_tag, mpi_error, mpi_recv
    !gcc$ attributes no_arg_check :: buf
    type(*) :: buf
    integer, intent(in) :: count
    type(mpi_datatype), intent(in) :: datatype
    integer, intent(in) :: dest
    integer, intent(in) :: tag
    type(mpi_comm), intent(in) :: comm
    type(mpi_status), intent(out) :: status
    integer, optional, intent(out) :: ierror

    integer status1(mpi_status_size)
    integer ierror1

    call mpi_recv(buf, count, datatype%datatype, dest, tag, comm%comm, status1, ierror1)

    status%status = status1(1:mpi_status_size)
    status%mpi_source = status1(mpi_source)
    status%mpi_tag = status1(mpi_tag)
    status%mpi_error = status1(mpi_error)

    if (present(ierror)) ierror = ierror1
  end subroutine mpi_recv_impl

  subroutine mpi_get_count_impl(status, datatype, count, ierror)
    use mpi, only: mpi_status_size, mpi_source, mpi_tag, mpi_error, mpi_get_count
    implicit none
    type(mpi_status), intent(in) :: status
    type(mpi_datatype), intent(in) :: datatype
    integer, intent(out) :: count
    integer, optional, intent(out) :: ierror

    integer status1(mpi_status_size)
    integer ierror1

    status1(1:mpi_status_size) = status%status
    status1(mpi_source) = status%mpi_source
    status1(mpi_tag) = status%mpi_tag
    status1(mpi_error) = status%mpi_error

    call mpi_get_count(status1, datatype%datatype, count, ierror1)

    if (present(ierror)) ierror = ierror1
  end subroutine mpi_get_count_impl

  subroutine mpi_comm_size_impl(comm, size, ierror)
    use mpi, only: mpi_comm_size
    implicit none
    type(mpi_comm), intent(in) :: comm
    integer, intent(out) :: size
    integer, optional, intent(out) :: ierror

    integer ierror1

    call mpi_comm_size(comm%comm, size, ierror1)

    if (present(ierror)) ierror = ierror1
  end subroutine mpi_comm_size_impl

  subroutine mpi_comm_rank_impl(comm, rank, ierror)
    use mpi, only: mpi_comm_rank
    implicit none
    type(mpi_comm), intent(in) :: comm
    integer, intent(out) :: rank
    integer, optional, intent(out) :: ierror

    integer ierror1

    call mpi_comm_rank(comm%comm, rank, ierror1)

    if (present(ierror)) ierror = ierror1
  end subroutine mpi_comm_rank_impl

  subroutine mpi_init_impl(ierror)
    use mpi, only: mpi_init
    implicit none
    integer, optional, intent(out) :: ierror

    integer ierror1

    call mpi_init(ierror1)

    if (present(ierror)) ierror = ierror1
  end subroutine mpi_init_impl

  subroutine mpi_finalize_impl(ierror)
    use mpi, only: mpi_finalize
    implicit none
    integer, optional, intent(out) :: ierror

    integer ierror1

    call mpi_finalize(ierror1)

    if (present(ierror)) ierror = ierror1
  end subroutine mpi_finalize_impl

  subroutine mpi_abort_impl(comm, errorcode, ierror)
    use mpi, only: mpi_abort
    implicit none
    type(mpi_comm), intent(in) :: comm
    integer, intent(in) :: errorcode
    integer, optional, intent(out) :: ierror

    integer ierror1

    call mpi_abort(comm%comm, errorcode, ierror1)

    if (present(ierror)) ierror = ierror1
  end subroutine mpi_abort_impl

end module mpi_f08

max_stars_repo_name: giordano/MPItrampoline

max_issues_repo_name: giordano/MPItrampoline

max_forks_repo_name: giordano/MPItrampoline


text:
      SUBROUTINE RHPOIL(NEED, F, N, NF, PT, R, RD, RHOI, YN, W)
      COMMON /FUDGE/ NFUDGE
      INTEGER NFUDGE
      INTEGER NEED(2), N, NF, RHOI(6)
      DOUBLE PRECISION F, PT(3), R(*), RD(*), W(N), YN(2,N)
C PT = PHI AND THETA (WHEN PS == P, I.E. RHOI(2) == RHOI(3))
C
      DOUBLE PRECISION INVCN, LPN, PNORMS, DR7MDC
      EXTERNAL INVCN, LPN, PNORMS, DR7MDC
      INTEGER ERRFLG, I, IM, WCOMP
      DOUBLE PRECISION CI, E, PHI, PHIRI, PHIMRI, PSI, PSI1, PSI2,
     1                 RI, T, T1, T2, THETA, YI
      DOUBLE PRECISION DATAN, DEXP, DLOG,  SQRT
      DOUBLE PRECISION CNN, EIGHT, EXPMAX, EXPMIN, FOUR, HALF, ONE, TWO,
     1                 TWOPI, ZERO
      DATA CNN/0.D+0/, EXPMAX/0.D+0/, EIGHT/8.D+0/, EXPMIN/0.D+0/,
     1     FOUR/4.0D+0/, HALF/0.5D+0/, ONE/1.D+0/, TWO/2.D+0/,
     2     TWOPI/0.D+0/, ZERO/0.D+0/
C
C *** BODY ***
C
      IM = RHOI(1)
      WCOMP = RHOI(6)
      IF (IM .LE. 0) GO TO 800
      IF (IM .GT. 13) GO TO 800
      IF (EXPMAX .GT. ZERO) GO TO 10
         EXPMAX = TWO * DLOG(DR7MDC(5))
         EXPMIN = TWO * DLOG(DR7MDC(2))
         TWOPI = EIGHT * DATAN(ONE)
 10   IF (NEED(1) .EQ. 2) GO TO 240
      F = ZERO
      GO TO (20,20,40,60,80,80,100,120,140,160,180,220,180), IM
C
C *** POISSON RHO (AND CONVENTIONAL IRLS) ***
C
 20   DO 30 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         F = F + YN(2,I)*RI - YN(1,I)*DLOG(RI)
 30      CONTINUE
      GO TO 999
C
C *** LOG LINEAR POISSON ***
C
 40   DO 50 I = 1, N
         E = ZERO
         RI = R(I)
         IF (RI .GT. EXPMAX) GO TO 800
         IF (RI .GT. EXPMIN) E = EXP(RI)
         F = F + YN(2,I)*E - YN(1,I)*RI
         R(I) = E
 50      CONTINUE
      GO TO 999
C
C *** SQUARE-ROOT LINEAR POISSON ***
C
 60   DO 70 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         F = F + YN(2,I)*RI**2 - TWO*YN(1,1)*DLOG(RI)
 70      CONTINUE
      GO TO 999
C
C *** BINOMIAL RHO (AND CONVENTIONAL IRLS) ***
C
 80   DO 90 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         IF (RI .GE. ONE) GO TO 800
         F = F - YN(1,I)*DLOG(RI) - (YN(2,I) - YN(1,I))*DLOG(ONE-RI)
 90      CONTINUE
      GO TO 999
C
C *** BINOMIAL LOGISTIC RHO ***
C
 100  DO 110 I = 1, N
         RI = R(I)
         IF (RI .GE. EXPMAX) GO TO 800
         E = ZERO
         IF (RI .GT. EXPMIN) E = DEXP(RI)
         F = F + YN(2,I)*DLOG(ONE + E) - YN(1,I)*RI
         R(I) = E
 110     CONTINUE
      GO TO 999
C
C *** PROBIT ***
C
 120  DO 130 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         F = F - YI*LPN(RI) - (YN(2,I)-YI)*LPN(-RI)
 130     CONTINUE
        IF (NFUDGE .GT. 0) WRITE(*,*) 'NFUDGE =', NFUDGE
        NFUDGE = 0
      GO TO 999
C
C *** WEIBULL ***
C
 140  DO 150 I = 1, N
         RI = R(I)
         IF (RI .GE. EXPMAX) GO TO 800
         E = ZERO
         IF (RI .GT. EXPMIN) E = DEXP(RI)
         R(I) = E
         T = ZERO
         IF (-E .GT. EXPMIN) T = DEXP(-E)
         F = F + (YN(2,I) - YN(1,I))*E - YN(1,I)*DLOG(ONE - T)
 150     CONTINUE
      GO TO 999
C
C  *** GAMMA ERRORS ***
C
 160  DO 170 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         F = F + YN(1,I)*RI - YN(2,I)*DLOG(RI)
 170     CONTINUE
      GO TO 999
C
C  ***  PREGIBON ERRORS ***
C
C      *** IN THIS CASE, YN(1,I) = Y(I), YN(2,I) = LOG(Y(I))
C      *** AND YN(I,J), J = N+1(1)2*N, I = 1 OR 2 = SCRATCH
C
 180  IF (NF .GT. 1) GO TO 190
      RHOI(4) = 0
      RHOI(5) = 0
 190  I = N + N + 3
C     *** THE YLOG ARRAY PASSED TO PREGRV MUST BE AT LEAST N+2 LONG
      IF (NEED(2) .NE. RHOI(4)) GO TO 200
         I = I + 3*N
         RHOI(5) = NF
         GO TO 210
 200  RHOI(4) = NF
 210  IF (IM .EQ. 11) THEN
        CALL PREGRH(0, F, N, NF, PT, R, RD, RHOI, YN(1,N+1), YN,
     1            YN(1,I))
      ELSE
        CALL PREGRV(0, F, N, NF, PT, R, RD, RHOI, YN(1,N+1), YN,
     1            YN(1,I))
        END IF
      GO TO 999
C
C *** LEAST-SQUARES ***
C
 220  DO 230 I = 1, N
        E = R(I) - YN(1,I)
        F = F + E*E
 230    CONTINUE
      F = HALF * F
      GO TO 999
C
 240  GO TO (250,270,310,350,400,420,460,500,570,620,660,780,660), IM
C
C *** IRLS POISSON DERIVATIVES ***
C
 250  DO 260 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         R(I) = YN(2,I) - YN(1,I) / RI
         RD(I) = YN(2,I) / RI
 260     CONTINUE
      GO TO 820
C
C *** POISSON DERIVATIVES ***
C
 270  DO 300 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         YI = YN(1,I)
         CI = YN(2,I)
         E = YI / RI
         R(I) = CI - E
         RD(I) = E / RI
         GO TO (300, 280, 280, 290), WCOMP
 280     W(I) = CI / RI
         GO TO 300
 290     IF (YI .LE. ZERO) THEN
             W(I) = HALF * CI / RI
         ELSE
            T1 = CI*RI + YI*(DLOG(E/CI) - ONE)
            IF (T1 .NE. ZERO) THEN
               T = R(I)
               W(I) = T*T / (T1+T1)
            ELSE
               W(I) = RD(I)
               END IF
            END IF
 300     CONTINUE
      GO TO 810
C
C *** LOG LINEAR POISSON ***
C
 310  DO 340 I = 1, N
         YI = YN(1,I)
         CI = YN(2,I)
         RI = CI*R(I)
         R(I) = RI - YI
         RD(I) = RI
         GO TO (340,340,320,330), WCOMP
 320     T = RI/YI
         IF (T .EQ. ONE) THEN
            W(I) = YI
         ELSE
            W(I) = YI * ((T - ONE) / DLOG(T))
            ENDIF
         GO TO 340
 330     T1 = RI + YI*(DLOG(YI/RI) - ONE)
         IF (T1 .NE. ZERO) THEN
            T = RI - YI
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 340     CONTINUE
      IF (WCOMP .LE. 2) GO TO 820
      GO TO 999
C
C *** SQUARE-ROOT LINEAR POISSON ***
C
 350  DO 390 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         YI = YN(1,I)
         CI = YN(2,I)
         E = YI / RI
         R(I) = TWO * (CI*RI - E)
         RD(I) = TWO * (CI + E/RI)
         GO TO (390, 360, 370, 380), WCOMP
 360     W(I) = FOUR * CI
         GO TO 390
 370     T1 = RI -  SQRT(YI/CI)
         IF (T1 .NE. ZERO) THEN
            T = CI*RI - YI/RI
            W(I) = (T+T) / T1
         ELSE
            W(I) = RD(I)
            END IF
         GO TO 390
 380     T1 = CI*RI*RI - YI + YI*DLOG(YI/(CI*RI*RI))
         IF (T1 .NE. ZERO) THEN
            T = CI*RI - YI/RI
            T = T / T1
            W(I) = T + T
         ELSE
            W(I) = RD(I)
            END IF
 390     CONTINUE
      GO TO 810
C
C *** IRLS BINOMIAL ***
C
 400  DO 410 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         IF (RI .GE. ONE) GO TO 800
         YI = YN(1,I)
         CI = YN(2,I)
         T = ONE / (ONE - RI)
         R(I) = (CI - YI) * T  -  YI / RI
         RD(I) = T * CI / RI
 410     CONTINUE
      GO TO 820
C
C *** BINOMIAL ***
C
 420  DO 450 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
         IF (RI .GE. ONE) GO TO 800
         YI = YN(1,I)
         T = ONE / (ONE - RI)
         CI = (YN(2,I) - YI) * T
         YI = YI / RI
         R(I) = CI - YI
         RD(I) = T*CI + YI/RI
         GO TO (450,430,430,440), WCOMP
 430     W(I) = T*YN(2,I) / RI
         GO TO 450
 440     YI = YN(1,I)
         CI = YN(2,I)
         T2 = YI / CI
         T1 = (YI - CI)*DLOG((ONE - RI)/(ONE - T2)) + YI*DLOG(T2/RI)
         IF (T1 .NE. ZERO) THEN
            T = (CI*RI - YI)/(RI * (ONE - RI))
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 450     CONTINUE
      GO TO 810
C
C *** BINOMIAL LOGISTIC ***
C
 460  DO 490 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         CI = YN(2,I)
         T = ONE / (ONE + RI)
         T1 = T * RI * CI
         R(I) = T1 - YI
         RD(I) = T * T1
         GO TO (490,490,470,480), WCOMP
 470     T1 = (ONE + RI)*DLOG(RI*(CI-YI)/YI)
         IF (T1 .NE. ZERO) THEN
            W(I) = ((CI - YI)*RI - YI) / T1
         ELSE
            W(I) = RD(I)
            END IF
         GO TO 490
 480     T1 = CI*DLOG((ONE+RI)*(ONE - YI/CI)) + YI*DLOG(YI/(RI*(CI-YI)))
         IF (T1 .NE. ZERO) THEN
            T = ((CI - YI)*RI - YI) / (ONE + RI)
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 490     CONTINUE
      IF (WCOMP .LE. 2) GO TO 820
      GO TO 999
C
C *** PROBIT ***
C
 500  IF (CNN .LE. ZERO) CNN = ONE /  SQRT(TWOPI)
      DO 560 I = 1, N
         RI = R(I)
         YI = YN(1,I)
         CI = YN(2,I) - YI
         E = ZERO
         T = -HALF * RI**2
         IF (T .GT. EXPMIN) E = CNN * DEXP(T)
         PHIRI = PNORMS(RI)
         IF (WCOMP .EQ. 2)
     1          W(I) = YN(2,I) * (E / PHIRI) * (E / (ONE - PHIRI))
         IF (PHIRI .LE. ZERO) GO TO 510
            PHIRI = ONE / PHIRI
            T1 = E*PHIRI*YI
            T2 = T1*(RI + PHIRI*E)
            T1 = -T1
            GO TO 520
 510     T1 = YI * (RI + ONE/RI)
         T2 = YI * (ONE - ONE/RI**2)
 520     PHIMRI = PNORMS(-RI)
         IF (PHIMRI .LE. ZERO) GO TO 530
            PHIMRI = ONE / PHIMRI
            T = E*CI*PHIMRI
            R(I) = T + T1
            RD(I) = T*(PHIMRI*E - RI) + T2
            GO TO (560,560,540,550), WCOMP
 530     R(I) = CI*(RI + ONE/RI) + T1
         RD(I) = CI*(ONE - ONE/RI**2) + T2
         GO TO (560,560,540,550), WCOMP
 540     T = RI - INVCN(YI/YN(2,I), ERRFLG)
         IF (ERRFLG .NE. 0) THEN
            WRITE(*,*) 'ERROR FROM INVCN: I, YI, YN(1,I), YN(2,I) ='
     1                  , I, YI, YN(1,I), YN(2,I)
             GO TO 800
             END IF
         IF (T .NE. ZERO) THEN
             W(I) = R(I) / T
         ELSE
             W(I) = RD(I)
             END IF
         GO TO 560
 550     T2 = CI
         CI = YN(2,I)
         T1 = T2*(DLOG(T2/CI) - LPN(-RI))
         IF (YI .GT. ZERO) T1 = T1 + YI*(DLOG(YI/CI) - LPN(RI))
         IF (T1 .NE. ZERO) THEN
             T = R(I)
             W(I) = T*T / (T1+T1)
         ELSE
             W(I) = RD(I)
             END IF
 560     CONTINUE
      GO TO 810
C
C *** WEIBULL ***
C
 570  DO 610 I = 1, N
         RI = R(I)
         E = ZERO
         IF (-RI .GT. EXPMIN) E = DEXP(-RI)
         T = RI / (ONE - E)
         CI = YN(2,I)*RI
         YI = YN(1,I)*T
         R(I) = CI - YI
         RD(I) = CI - YI*(ONE - E*T)
         GO TO (570,580,590,600), WCOMP
 580     W(I) = E*CI*RI / (ONE - E)
         GO TO 610
 590     T1 = DLOG(-RI / DLOG(ONE - YN(1,I)/YN(2,I)))
         IF (T1 .NE. ZERO) THEN
            W(I) = (CI - YI) / T1
         ELSE
            W(I) = RD(I)
            END IF
         GO TO 610
 600     YI = YN(1,I)
         CI = YN(2,I)
         T2 = YI / CI
         CI = CI - YI
         T1 = CI*(RI + DLOG(ONE - T2)) + YI*(DLOG(T2/(ONE - E)))
         IF (T1 .NE. ZERO) THEN
            T = CI - YI
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 610     CONTINUE
      GO TO 810
C
C  *** GAMMA ERRORS ***
C
 620  DO 650 I = 1, N
         RI = R(I)
         IF (RI .LE. ZERO) GO TO 800
C        F = F + YN(1,I)*RI - YN(2,I)*DLOG(RI)
         T = YN(2,I)/RI
         T1 = ONE
         R(I) = YN(1,I) - T
         RD(I) = T/RI
         GO TO (650,650,630,640), WCOMP
 630     W(I) = YN(1,I) / RI
         GO TO 650
 640     T2 = YN(1,I) * RI / YN(2,I)
         T1 = T2 - ONE
         T = T1*RD(I)*T1
         IF (T .GT. ZERO) THEN
            T2 = T1 - DLOG(T2)
            T = T / (T2+T2)
            END IF
         W(I) = T
 650     CONTINUE
      IF (WCOMP .LE. 2) GO TO 820
      GO TO 999
C
C ***  PREGIBON ERRORS ***
C
 660  IF (WCOMP .GE. 2) CALL DV7CPY(N, W, R)
      I = N + N + 3
      IF (RHOI(4) .EQ. NF) GO TO 670
         I = I + 3*N
         IF (RHOI(5) .EQ. NF) GO TO 670
         WRITE(6,*) 'HELP! NF =', NF, ' BUT RHOI =', RHOI
         GO TO 800
 670  IF (IM .EQ. 11) THEN
         CALL PREGRH(1, F, N, NF, PT, R, RD, RHOI, YN(1,N+1), YN,
     1            YN(1,I))
      ELSE
         CALL PREGRV(1, F, N, NF, PT, R, RD, RHOI, YN(1,N+1), YN,
     1            YN(1,I))
         END IF
      IF (NF .EQ. 0) GO TO 999
      GO TO (820,680,700,720), WCOMP
 680  PSI = PT(3)
      T = (TWO - PT(2))*PSI - TWO
      T1 = PSI*PSI
      DO 690 I = 1, N
 690        W(I) = YN(2,I) * T1 * W(I)**T
      GO TO 999
 700  T = ONE / PT(3)
      DO 710 I = 1, N
         T1 = W(I) - ONE
         IF (T1 .NE. ZERO) THEN
            YI = YN(1,I)
            W(I) = R(I) / (W(I) - YI**T)
         ELSE
            W(I) = RD(I)
            END IF
 710     CONTINUE
      GO TO 999
 720  PHI = PT(1)
      THETA = PT(2)
      PSI = PT(3)
      IF (THETA .EQ. ONE) GO TO 740
      IF (THETA .EQ. TWO) GO TO 760
      T1 = ONE - THETA
      T2 = TWO - THETA
      PSI1 = PSI * T1
      PSI2 = PSI * T2
      DO 730 I = 1, N
         RI = W(I)
         YI = YN(1,I)
         T = YI**T2
         E = YN(2,I)/PHI * ((T - YI*RI**PSI1)/T1 - (T - RI**PSI2)/T2)
         IF (E .NE. ZERO) THEN
            T = R(I)
            W(I) = T*T / (E+E)
         ELSE
            W(I) = RD(I)
            END IF
 730     CONTINUE
      GO TO 999
 740  DO 750 I = 1, N
         RI = W(I)
         YI = YN(1,I)
         T1 = YN(2,I)/PHI * (RI**PSI - YI + YI*(DLOG(YI)-PSI*DLOG(RI)))
         IF (T1 .NE. ZERO) THEN
            T = R(I)
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 750     CONTINUE
      GO TO 999
 760  DO 770 I = 1, N
         RI = W(I)
         YI = YN(1,I)
         T1 = YI*RI**(-PSI) - ONE + PSI*DLOG(RI) - DLOG(YI)
         IF (T1 .NE. ZERO) THEN
            T1 = T1 * YN(2,I) / PHI
            T = R(I)
            W(I) = T*T / (T1+T1)
         ELSE
            W(I) = RD(I)
            END IF
 770     CONTINUE
      GO TO 999
C
C *** LEAST SQUARES ***
C
 780  DO 790 I = 1, N
         R(I) = R(I) - YN(1,I)
         RD(I) = ONE
 790     CONTINUE
      GO TO 820
C
 800  NF = 0
      GO TO 999
C
 810  IF (WCOMP .GT. 1) GO TO 999
 820  CALL DV7CPY(N, W, RD)
C
 999  RETURN
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
There are two old men that hang out at the very center of the MU nearly every day. Right by the entrance to the section of the Coffee House that has a piano. One old man always has on a fishing hat and the other one is bald. Rumor has it from Users/GarySaylin that one of them is an emeritus art professor. Apparently they are two Town Characters that are worth talking to.

The guy in the fishing hat taught an art history freshman seminar in winter 2002. It was the craziest class I ever took. I serve him on the hot food line sometimes. He doesnt remember me, but tells me Im a wonderful gentleman. JoshuaHeller

20080930 11:40:42 nbsp That would be Professor Emeritus Seymour Howard...he practices Tai Chi on campus last I heard, on Sunday mornings over by the Surge buildings. Hes a great guy to talk with about history, art, politics, whatever. Users/ArthurFrane

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
module ControlParameterClass
    use, intrinsic :: iso_fortran_env
    use std
    implicit none
    type::ControlParameter_
        real(real64)   :: Tol
        integer(int32) :: SimMode
        integer(int32) :: ItrTol
        integer(int32) :: Timestep
    contains
        procedure,public :: init => initControlPara
        procedure,public :: remove => initControlPara
        procedure,public :: SetControlPara => SetControlPara
        procedure,public :: set => SetControlPara
        procedure,public :: export => exportControlPara
        procedure,public :: open => openControlPara
        procedure,public :: save => saveControlPara
    end type ControlParameter_

contains

! ########################################################
subroutine openControlPara(obj,path,name)
    class(ControlParameter_),intent(inout):: obj

    character(*),intent(in) :: path
    character(*),optional,intent(in) :: name
    type(IO_) :: f

    if(present(name) )then
        call f%open(trim(path)//"/"//trim(adjustl(name))//"/","ControlPara",".prop")
        read(f%fh,*) obj%Tol
        read(f%fh,*) obj%SimMode
        read(f%fh,*) obj%ItrTol
        read(f%fh,*) obj%Timestep
        call f%close()
    else
        call f%open(trim(path)//"/ControlPara/","ControlPara",".prop")
        read(f%fh,*) obj%Tol
        read(f%fh,*) obj%SimMode
        read(f%fh,*) obj%ItrTol
        read(f%fh,*) obj%Timestep
        call f%close()
    endif
end subroutine
! ########################################################

! ########################################################
subroutine saveControlPara(obj,path,name)
    class(ControlParameter_),intent(inout):: obj

    character(*),intent(in) :: path
    character(*),optional,intent(in) :: name
    type(IO_) :: f

    if(present(name) )then
        call execute_command_line("mkdir -p "//trim(path)//"/"//trim(adjustl(name)))
        call f%open(trim(path)//"/"//trim(adjustl(name))//"/","ControlPara",".prop")
        write(f%fh,*) obj%Tol
        write(f%fh,*) obj%SimMode
        write(f%fh,*) obj%ItrTol
        write(f%fh,*) obj%Timestep
        call f%close()
    else
        call execute_command_line("mkdir -p "//trim(path)//"/ControlPara")
        call f%open(trim(path)//"/ControlPara/","ControlPara",".prop")
        write(f%fh,*) obj%Tol
        write(f%fh,*) obj%SimMode
        write(f%fh,*) obj%ItrTol
        write(f%fh,*) obj%Timestep
        call f%close()
    endif
end subroutine
! ########################################################

! ########################################################
subroutine initControlPara(obj)
    class(ControlParameter_),intent(inout):: obj

    obj%Tol=1.0e-16
    obj%SimMode=1
    obj%ItrTol=1000
    obj%Timestep=1
end subroutine
! ########################################################

!######### Import Control Parameters #########
subroutine SetControlPara(obj,OptionalTol,OptionalItrTol,OptionalTimestep,OptionalSimMode)
    class(ControlParameter_),intent(out)::obj
    real(real64),optional,intent(in)::OptionalTol
    integer(int32),optional,intent(in)::OptionalSimMode,OptionalItrTol,OptionalTimestep

    if(present(OptionalSimMode) )then
        obj%SimMode = OptionalSimMode
    else
        obj%SimMode = 1
    endif

    if(present(OptionalTol) )then
        obj%Tol = OptionalTol
    else
        obj%Tol = 1.0e-15
    endif

    if(present(OptionalItrTol) )then
        obj%ItrTol = OptionalItrTol
    else
        obj%ItrTol = 100
    endif

    if(present(OptionalTimestep) )then
        obj%Timestep = OptionalTimestep
    else
        obj%Timestep = 100
    endif

end subroutine SetControlPara
!######### Import Control Parameters #########

subroutine exportControlPara(obj,restart,path)
    class(ControlParameter_),intent(inout) :: obj
    logical,optional,intent(in) :: restart
    character(*),intent(in) :: path
    type(IO_) :: f

    if(present(restart) )then
        call execute_command_line("mkdir -p "//trim(path)//"/ControlPara")
        call f%open(trim(path)//"/ControlPara/","ControlPara",".res")
        write(f%fh,*) obj%Tol
        write(f%fh,*) obj%SimMode
        write(f%fh,*) obj%ItrTol
        write(f%fh,*) obj%Timestep
        call f%close()
        return
    endif

end subroutine

end module ControlParameterClass

max_stars_repo_name: kazulagi/plantFEM

max_issues_repo_name: kazulagi/plantFEM

max_forks_repo_name: kazulagi/plantFEM


text:
! RUN: %flang_fc1 -fsyntax-only %s
! RUN: rm -rf %t && mkdir %t
! RUN: touch %t/empty.f90
! RUN: %flang_fc1 -fsyntax-only %t/empty.f90

max_stars_repo_name: mkinsner/llvm

max_issues_repo_name: mkinsner/llvm

max_forks_repo_name: mkinsner/llvm