text:
!***********************************************************************
!
!     file  lusol8b.f
!
!     lu8adc   lu8adr   lu8dlc   lu8dlr   lu8mod   lu8rpr
!
!     These routines call lu7asv and lu7bak in lusol7b.f.
!
! 08 Jun 2004: lusol8b.f is essentially lu8b.for from VMS days.
!              integer*4 changed to integer  .
! 15 Sep 2004: Test nout. gt. 0 to protect write statements.
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu8adc( mode, m, n, v, w,
     $                   lena, luparm, parmlu,
     $                   a, indc, indr, ip, iq,
     $                   lenc, lenr, locc, locr,
     $                   inform, diag, vnorm )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena), v(m), w(n)
      integer            indc(lena), indr(lena), ip(m), iq(n)
      integer            lenc(n), lenr(m)
      integer            locc(n), locr(m)

!     ------------------------------------------------------------------
!     lu8adc  updates the LU factorization  A = L*U  when the vector
!     a(new)  is added to become the  n-th  column of  A.
!     The dimension of  A  is assumed to increase from  n-1  to  n
!     (not from  n  to  n+1).
!
!     If  mode = 1,  v(*)  must contain  a(new).
!     If  mode = 2,  v(*)  must satisfy  L*v = a(new).
!     On exit,  L*v = a(new)  in both cases.
!
!     The array  w(*)  is not used or altered.
!
!     On entry, all elements of  locc  are assumed to be zero.
!     On a successful exit (inform ne 7), this will again be true.
!
!     On exit:
!     inform =  0  if the rank of U stayed the same.
!     inform =  1  if the rank of U increased by 1.
!     inform =  7  if the update was not completed (lack of storage).
!
!     -- Feb 1985: Last  F66 version.
!     12 May 1988: First F77 version.  Now uses lu8rpc.
!     ------------------------------------------------------------------

      nout   = luparm(1)
      lprint = luparm(2)
      nrank  = luparm(16)

      ! Set  locc(n)  in case the user forgot.
      ! Then let lu8rpc do the job (but stop it from printing messages).

      locc(n)   =  0
      iq(n)     =  n
      luparm(2) = -1

      call lu8rpc( 0, mode, m, n, n, v, w,
     $             lena, luparm, parmlu,
     $             a, indc, indr, ip, iq,
     $             lenc, lenr, locc, locr,
     $             inform, diag, vnorm )

      if (inform .le. 0) then
         if (nrank .eq. n - 1) then
            if (nout. gt. 0  .and.  lprint .ge. 0)
     &           write(nout, 1100) m, n, diag, vnorm
         end if
      else if (inform .eq. 7) then
         if (nout. gt. 0  .and.  lprint .ge. 0)
     &        write(nout, 1700) lena
      end if

      luparm(2) = lprint
      return

 1100 format(/ ' lu8adc  warning.  Rank did not increase',
     $         ' after adding a column.'
     $       / ' m =', I8, '    n =', i8,
     $         '    diag =', 1p, e12.2, '    vnorm =', e12.2)
 1700 format(/ ' lu8adc  error...  Insufficient storage.',
     $         '    lena =', I8)

      end ! subroutine lu8adc

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu8adr( m, n, w,
     $                   lena, luparm, parmlu,
     $                   a, indc, indr, ip, iq,
     $                   lenc, lenr, locc, locr,
     $                   inform, diag )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena), w(n)
      integer            indc(lena), indr(lena), ip(m), iq(n)
      integer            lenc(n), lenr(m)
      integer            locc(n), locr(m)

!     ------------------------------------------------------------------
!     lu8adr  updates the LU factorization  A = L*U  when the vector
!     w(*)  is added to become the  m-th  row of  A.
!     The dimension of  A  is assumed to increase from  m-1  to  m
!     (not from  m  to  m+1).
!
!     w(*)  is not altered.
!
!     On entry, all elements of  locc  are assumed to be zero.
!     On a successful exit (inform ne 7), this will again be true.
!
!     On exit:
!     inform =  0  if the rank of U stayed the same.
!     inform =  1  if the rank of U increased by 1.
!     inform =  7  if the update was not completed (lack of storage).
!
!     -- Feb 1985: Last  F66 version.
!     17 May 1988: First F77 version.
!     ------------------------------------------------------------------

      parameter        ( zero = 0.0d+0 )

      nout   = luparm(1)
      lprint = luparm(2)
      nrank  = luparm(16)
      lenl   = luparm(23)
      lenu   = luparm(24)
      lrow   = luparm(25)
      small  = parmlu(3)
      diag   = zero
      nrank1 = nrank + 1

      ! Compress row file if necessary.

      minfre = n
      nfree  = lena - lenl - lrow
      if (nfree .ge. minfre) go to 100
      call lu1rec( m, .true., luparm, lrow, ilast,
     $             lena, a, indr, lenr, locr )
      nfree  = lena - lenl - lrow
      if (nfree .lt. minfre) go to 970

      ! Pack the nonzeros of  w  at the end of the row file.
      ! Go backwards in order to set  kfirst.

  100 locr(m) = lrow + 1
      ip(m)   = m
      lenr(m) = 0

      do 120 k = n, 1, -1
         j          = iq(k)
         if (abs( w(j) ) .le. small) go to 120
         kfirst     = k
         lrow       = lrow + 1
         a(lrow)    = w(j)
         indr(lrow) = j
  120 continue

      lenw    = lrow + 1 - locr(m)
      lenu    = lenu + lenw
      lenr(m) = lenw

!     ------------------------------------------------------------------
!     Triangularize the new row.
!     ------------------------------------------------------------------

      if (lenw .gt. 0) then

         ! Swap the new row into position nrank1.

         ip(m)      = ip(nrank1)
         ip(nrank1) = m

         ! Perform a forward sweep to eliminate subdiagonal elements.

         if (kfirst .le. nrank1) then
            call lu7for( m, n, kfirst, nrank1,
     $                   lena, luparm, parmlu,
     $                   lenl, lenu, lrow,
     $                   a, indc, indr, ip, iq, lenr, locc, locr,
     $                   inform, diag )
            if (inform .eq. 7) go to 970
         end if

         ! See if the rank increases.

         if (nrank .lt. n) then
            nrank  = nrank1
            call lu7rnk( m, n, 0,
     $                   lena, luparm, parmlu,
     $                   lenl, lenu, lrow, nrank,
     $                   a, indc, indr, ip, iq, lenr, locc, locr,
     $                   inform, diag )
         end if
      end if

!     ------------------------------------------------------------------
!     Set inform for normal exit.
!     ------------------------------------------------------------------
      if (nrank .ne. nrank1) then
         inform = 0
         if (nrank .eq. m - 1) then
            if (nout. gt. 0  .and.  lprint .ge. 0)
     &           write(nout, 1100) m, n, diag
         end if
      else
         inform = 1
      end if
      go to 990

      ! Not enough storage.

  970 inform = 7
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1700) lena

      ! Exit.

  990 luparm(10) = inform
      luparm(15) = luparm(15) + 1
      luparm(16) = nrank
      luparm(23) = lenl
      luparm(24) = lenu
      luparm(25) = lrow
      return

 1100 format(/ ' lu8adr  warning.  Rank did not increase',
     $         ' after adding a row.'
     $       / ' m =', i8, '    n =', i8,
     $         '    diag =', 1p, e12.2)
 1700 format(/ ' lu8adr  error...  Insufficient storage.',
     $         '    lena =', i8)

      end ! subroutine lu8adr

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu8dlc( m, n, jdel,
     $                   lena, luparm, parmlu,
     $                   a, indc, indr, ip, iq,
     $                   lenc, lenr, locc, locr,
     $                   inform )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena)
      integer            indc(lena), indr(lena), ip(m), iq(n)
      integer            lenc(n), lenr(m)
      integer            locc(n), locr(m)

!     ------------------------------------------------------------------
!     lu8dlc  updates the LU factorization  A = L*U  when column  jdel
!     is deleted.  The dimension of  A  decreases from  n  to  n-1.
!     The column permutation  iq(*)  and the column indices in  U
!     are altered accordingly.
!
!     NOTE:  The calling program must change  n  to  n-1.
!
!     In some cases it may be more efficient to use  lu8rpc  to
!     replace column  jdel  by zero, leaving it in the current position.
!
!     On entry, all elements of  locc  are assumed to be zero.
!     On a successful exit (inform ne 7), this will again be true.
!
!     On exit:
!     inform = -1  if the rank of U decreased by 1.
!     inform =  0  if the rank of U stayed the same.
!     inform =  7  if the update was not completed (lack of storage).
!
!     -- Feb 1985: Last  F66 version.
!     17 May 1988: First F77 version.
!     ------------------------------------------------------------------

      nout   = luparm(1)
      lprint = luparm(2)
      nrank  = luparm(16)
      lenl   = luparm(23)
      lenu   = luparm(24)
      lrow   = luparm(25)
      n1     = n - 1
      if (jdel .lt. 1) go to 980
      if (jdel .gt. n) go to 980

      ! Remove column jdel from U, and set kdel so that iq(kdel) = jdel.

      call lu7zap( m, n, jdel, kdel,
     $             lena, lenu, lrow, nrank,
     $             a, indr, ip, iq, lenr, locr )

      ! Renumber columns of  U  that are to the right of column  jdel.
      ! In effect, those columns are shifted one place to the left.

      if (jdel .lt. n) then
         do k = 1, nrank
            i      = ip(k)
            lr1    = locr(i)
            lr2    = lr1 + lenr(i) - 1

            do l = lr1, lr2
               j     = indr(l)
               if (j .gt. jdel) indr(l) = j - 1
            end do
         end do

         do k = 1, n
            j     = iq(k)
            if (j .gt. jdel) iq(k) = j - 1
         end do
      end if

      ! Perform cyclic permutations to move column kdel to the end
      ! and the corresponding row to position nrank.
      ! Then eliminate the resulting row spike.

      call lu7cyc( kdel, nrank, ip )
      call lu7cyc( kdel, n    , iq )

      call lu7for( m, n1, kdel, nrank,
     $             lena, luparm, parmlu,
     $             lenl, lenu, lrow,
     $             a, indc, indr, ip, iq, lenr, locc, locr,
     $             inform, diag )
      if (inform .eq. 7) go to 970

      ! See if the rank decreased.

      CALL LU7RNK( M, N1, 0,
     $             LENA, LUPARM, PARMLU,
     $             LENL, LENU, LROW, NRANK,
     $             A, INDC, INDR, IP, IQ, LENR, LOCC, LOCR,
     $             INFORM, DIAG )
      go to 990

      ! Not enough storage.

  970 inform = 7
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1700) lena
      go to 990

      ! jdel  is out of range.

  980 inform = 8
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1800) m, n, jdel

      ! Exit.

  990 luparm(10) = inform
      luparm(15) = luparm(15) + 1
      luparm(16) = nrank
      luparm(23) = lenl
      luparm(24) = lenu
      luparm(25) = lrow
      return

 1700 format(/ ' lu8dlc  error...  Insufficient storage.',
     $         '    lena =', i8)
 1800 format(/ ' lu8dlc  error...  jdel  is out of range.',
     $         '    m =', i8, '    n =', i8, '    jdel =', i8)

      end ! subroutine lu8dlc

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu8mod( mode, m, n, beta, v, w,
     $                   lena, luparm, parmlu,
     $                   a, indc, indr, ip, iq,
     $                   lenc, lenr, locc, locr,
     $                   lenlv, inform )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena), v(m), w(n)
      integer            indc(lena), indr(lena), ip(m), iq(n)
      integer            lenc(n), lenr(m)
      integer            locc(n), locr(m)
      integer            lenlv(3)
!     ------------------------------------------------------------------
!     lu8mod  updates the LU factorization  A = L*U  when the
!     m by n matrix  A  is subjected to a rank-one modification
!     to become
!                    A  +  beta * v * w(transpose)
!     for the given scalar  beta  and vectors  v, w.
!
!     If mode = 1, v(*) must contain the vector  v.
!                  On exit, v(*) will contain  y  satisfying  L*y = v.
!
!     If MODE = 2, v(*) must already contain   y  satisfying  L*y = v.
!
!     In both cases, w(*) is unaltered.
!
!     On entry, the elements of  locc  are assumed to be zero.
!     On a successful exit (inform ne 7), this will again be true.
!
!     On exit:
!     inform = -1  if the rank of U decreased by 1.
!     inform =  0  if the rank of U stayed the same.
!     inform =  1  if the rank of U increased by 1.
!     inform =  7  if the update was not completed (lack of storage).
!
!     -- Mar 1985: Last  F66 version.
!     18 May 1988: First F77 version.
!     ------------------------------------------------------------------

      logical            singlr

      nout   = luparm(1)
      lprint = luparm(2)
      nrank  = luparm(16)
      lenl   = luparm(23)
      lenu   = luparm(24)
      lrow   = luparm(25)
      small  = parmlu(3)
      utol1  = parmlu(4)
      nrank0 = nrank
      lenl0  = lenl
      singlr = .FALSE.
      ! If necessary, solve  L*v(new) = v.

      if (mode .eq. 1) then
         call lu6sol( 1, m, n, v, w, lena, luparm, parmlu,
     $                a, indc, indr, ip, iq, lenc, lenr, locc, locr,
     $                inform )
      end if

!     ------------------------------------------------------------------
!     Find the first nonzero in  w  (in pivotal column order).
!     ------------------------------------------------------------------
      do k = 1, n
         kfirst = k
         j      = iq(k)
         if (abs( w(j) ) .gt. small) go to 120
      end do
      go to 900

!     ------------------------------------------------------------------
!     Eliminate any nonzeros in  v  below the trapezoid.
!     --- This shouldn't ever happen. REMEMBER TO ERASE THIS ------
!     ------------------------------------------------------------------
  120 if (nrank .lt. m) then
         jelm   = 0

         call lu7elm( m, n, jelm, v,
     $                lena, luparm, parmlu,
     $                lenl, lenu, lrow, nrank,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )

         if (inform .eq. 1) nrank = nrank+1;
         if (inform .eq. 7) go to 970

      end if

!     ------------------------------------------------------------------
!     Find the last nonzero in  v  (in pivotal row order).
!     ------------------------------------------------------------------
      do k = nrank, 1, -1
         klast  = k
         i      = ip(k)
         if (abs( v(i) ) .gt. small) go to 220
      end do
      go to 900

!     ------------------------------------------------------------------
!     Perform a backward sweep of eliminations to reduce part of  v
!     to a multiple of the unit vector  e(iw),  where  iw = ip(klast).
!     Elements  ip(kfirst+1),  ip(kfirst+2),  ...,  ip(klast)  of  v
!     are involved.
!     L, U  and  ip  are updated accordingly.
!     U  will then be trapezoidal except for row  iw = ip(klast).
!     ------------------------------------------------------------------
  220 if (kfirst+1 .lt.  klast) then
         call lu7bak( m, n, kfirst, klast, v,
     $                lena, luparm, parmlu,
     $                lenl, lenu, lrow,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform )
         if (inform .ne. 0) go to 970
      end if
      lenlv(1) = lenl - lenl0
      lenl0 = lenl
!     ------------------------------------------------------------------
!     Pack the nonzeros of  w  in pivotal order in front of  L.
!     (We will treat the packed  w  much like a normal row of  U.)
!     Set markers on  w  and initialize the corresponding
!     elements of  indc(*),  which are later used by  lu7asv.
!     ------------------------------------------------------------------

      lfree  = lena - lenl
      minfre = n + 1 - kfirst
      nfree  = lfree - lrow
      if (nfree .ge. minfre) go to 310
      call lu1rec( m, .true., luparm, lrow, ilast,
     $             lena, a, indr, lenr, locr )
      nfree  = lfree - lrow
      if (nfree .lt. minfre) go to 970

  310 lw     = lfree + 1

      do 320   k  = n, kfirst, -1
         j        = iq(k)
         if (abs( w(j) ) .le. small) go to 320
         lw       = lw - 1
         a(lw)    = w(j)
         indr(lw) = j
         indc(lw) = 0
         locc(j)  = lw
  320 continue

      lw1    = lw
      lw2    = lfree
      lenw   = lw2 + 1 - lw1
      lfree  = lfree - lenw

!     ------------------------------------------------------------------
!     Add multiples of  w  to the first  kfirst  rows of  U.
!     (This does not alter the trapezoidal form of  U.)
!     ------------------------------------------------------------------

      do 450 k  = 1, kfirst
         iv     = ip(k)
         cv     = v(iv)
         if (abs( cv ) .le. small) go to 450

!        ===============================================================
!        Compress storage if necessary, so there will be room if
!        row  iv  has to be moved to the end.
!        ===============================================================
         minfre = n
         nfree  = lfree - lrow
         if (nfree .ge. minfre) go to 420
         call lu1rec(m, .true., luparm, lrow,
     $               ilast, lena, a,indr,lenr, locr)
         nfree  = lfree - lrow
         if (nfree .lt. minfre) go to 970

!        ===============================================================
!        Set  v  =  v  +  wmult * w.
!        ===============================================================
  420    wmult  = beta * cv
         call lu7asv( m, n, iv, lenw, lw1, lw2, k, wmult,
     $                lena, luparm, parmlu,
     $                lenu, lrow,
     $                a, indc, indr, lenr, locc, locr )
  450 continue

!     ------------------------------------------------------------------
!     Add a multiple of  w  to row  iw  of  U.
!     ------------------------------------------------------------------

      if (kfirst .lt. klast) then
         minfre = n
         nfree  = lfree - lrow
         if (nfree .ge. minfre) go to 500
         call lu1rec( m, .true., luparm, lrow, ilast,
     $                lena, a,indr,lenr,locr )
         nfree  = lfree - lrow
         if (nfree .lt. minfre) go to 970

  500    iw     = ip(klast)
         marker = m + 1
         wmult  = beta * v(iw)
         call lu7asv( m, n, iw, lenw, lw1, lw2, marker, wmult,
     $                lena, luparm, parmlu,
     $                lenu, lrow,
     $                a, indc, indr, lenr, locc, locr )
      end if

!     ------------------------------------------------------------------
!     Cancel the markers on  w.
!     ------------------------------------------------------------------
      do lw  = lw1, lw2
         jw       = indr(lw)
         locc(jw) = 0
      end do
!     ------------------------------------------------------------------
!     Apply a forward sweep to eliminate the nonzeros in row  iw.
!     ------------------------------------------------------------------
      if (kfirst .gt. klast) then
         if (klast .lt. nrank) go to 900
      else
         call lu7for( m, n, kfirst, klast,
     $                lena, luparm, parmlu,
     $                lenl, lenu, lrow,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )
         if (inform .eq. 7) go to 970
         lenlv(2) = lenl - lenl0
         lenl0 = lenl
      end if

!     ------------------------------------------------------------------
!     Test for singularity in column klast (if klast .le. nrank).
!     The code is similar to part of lu8rpc with klast in place of krep.
!     ------------------------------------------------------------------

      if (klast .le. nrank) then
         diag   = zero
         iw     = ip(klast)
         singlr = lenr(iw) .eq. 0

         if (.not. singlr) then
            l1     = locr(iw)
            j1     = indr(l1)
            singlr = j1 .ne. jrep

            if (.not. singlr) then
               diag   = a(l1)
               singlr = abs( diag ) .le. utol1
            end if
         end if

         if ( singlr  .and.  klast .lt. nrank ) then

            ! Perform cyclic permutations to move column klast
            ! to the end and the corresponding row to position nrank.
            ! Then eliminate the resulting row spike.

            call lu7cyc( klast, nrank, ip )
            call lu7cyc( klast, n    , iq )

            call lu7for( m, n, klast, nrank,
     $                   lena, luparm, parmlu,
     $                   lenl, lenu, lrow,
     $                   a, indc, indr, ip, iq, lenr, locc, locr,
     $                   inform, diag )
            lenlv(3) = lenl - lenl0
            if (inform .eq. 7) go to 970
         end if
         ! Find the best column to be in position nrank.
         ! If nothing satisfactory exists, nrank will be decreased.

         jsing  = 0
         call lu7rnk( m, n, jsing,
     $                lena, luparm, parmlu,
     $                lenl, lenu, lrow, nrank,
     $                a, indc, indr, ip, iq, lenr, locc, locr,
     $                inform, diag )
      end if
!     ------------------------------------------------------------------
!     Set inform for exit.
!     ------------------------------------------------------------------

  900 if (nrank .eq. nrank0) then
         inform =  0
      else if (nrank .lt. nrank0) then
         inform = -1
      else
         inform =  1
      end if
      go to 990

      ! Not enough storage.

  970 inform = 7
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1700) lena

      ! Exit.

  990 luparm(10) = inform
      luparm(15) = luparm(15) + 1
      luparm(16) = nrank
      luparm(23) = lenl
      luparm(24) = lenu
      luparm(25) = lrow
      return

 1700 format(/ ' lu8mod  error...  Insufficient storage.',
     $         '    lena =', i8)

      end ! subroutine lu8mod

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine lu8rpr( mode1, mode2, m, n, irep, v, w, wnew,
     $                   lena, luparm, parmlu,
     $                   a, indc, indr, ip, iq,
     $                   lenc, lenr, locc, locr,
     $                   lenlv, li, lj, lv,
     $                   inform )

      implicit           double precision (a-h,o-z)
      integer            luparm(30)
      double precision   parmlu(30), a(lena), v(m), w(n), wnew(n)
      integer            indc(lena), indr(lena), ip(m), iq(n)
      integer            lenc(n), lenr(m)
      integer            locc(n), locr(m)
      integer            lenlv(3),li(3*m), lj(3*m)
      double precision   lv(3*m)

!     ------------------------------------------------------------------
!     lu8rpr  updates the LU factorization  A = L*U  when row  irep
!     (a vector  w(old)) is replaced by some vector  w(new).
!     The update is implemented as the rank-one modification
!
!           A(new)  =  A  -  e(irep) * ( w(old) - w(new) )'
!
!     with variations determined by  mode1  and  mode2  as follows.
!
!     mode1
!     -----
!       0     w(old)  is assumed to be zero.
!             w(*)    need not be set on entry, but will be altered.
!
!       1     w(old)  is assumed to be unknown.  it will be computed
!                     from the LU factors of  A.
!             w(*)    need not be set on entry, but will be altered.
!
!       2     w(*)    must contain  w(old).
!                     On exit, it will contain  w(old) - w(new).
!
!       3     w(*)    must contain  w(old) - w(new).  It is not altered.
!                     wnew(*) is not used.
!
!       4     w(*)    must contain  w(new) - w(old).  It is not altered.
!                     wnew(*) is not used.
!
!     If  mode1 = 3 or 4,  mode2  is not used.  It may be set to 0 or 1.
!     Otherwise,  mode2  is used as follows.
!
!     mode2
!     -----
!       0     w(new)  is assumed to be zero.   wnew(*)  is not used.
!
!       1     wnew(*) must contain the new row.
!
!     v(*)    is a work array of length  m.
!     On entry, all elements of  locc  are assumed to be zero.
!     On a successful exit (inform ne 7), this will again be true.
!
!
!     On exit:
!     inform = -1  if the rank of U decreased by 1.
!     inform =  0  if the rank of U stayed the same.
!     inform =  1  if the rank of U increased by 1.
!     inform =  7  if the update was not completed (lack of storage).
!     inform =  8  if jrep is not between 1 and n.
!
!     -- Mar 1985: Last  F66 version.
!     18 May 1988: First F77 version.
!     ------------------------------------------------------------------

      parameter        ( zero = 0.0d+0,  one = 1.0d+0 )

      lenL0 = luparm(23)
      lenlv(1:3) = 0
      nout   = luparm(1)
      lprint = luparm(2)
      if (irep .lt. 1) go to 980
      if (irep .gt. m) go to 980

      if (mode1 .eq. 0) then

         ! The old row  irep  is zero.

         do j = 1, n
            w(j) = zero
         end do
      else if (mode1 .eq. 1) then

         ! Compute row irep as the vector  w  =  A(transpose) * e(irep).

         do i = 1, m
            v(i)  = zero
         end do
         v(irep) = one
         call lu6mul( 6, m, n, v, w, lena, luparm, parmlu,
     $                a, indc, indr, ip, iq, lenc, lenr, locc, locr )
      end if

      if (mode1 .le. 2) then

         ! Compute the difference except if  wnew = 0.

         if (mode2 .gt. 0) then
            do j = 1, n
               w(j)  = w(j) - wnew(j)
            end do
         end if
      end if

!     ------------------------------------------------------------------
!     Set  v = a unit vector  and do the rank-one modification
!     (but don't let lu8mod print anything).
!     ------------------------------------------------------------------

      do i = 1, m
         v(i)  = zero
      end do

      v(irep)   =   one
      beta      = - one
      if (mode1 .eq. 4) beta = one
      luparm(2) = - 1

      call lu8mod( 1, m, n, beta, v, w,
     $             lena, luparm, parmlu,
     $             a, indc, indr, ip, iq, lenc, lenr, locc, locr,
     $             lenlv, inform )
      if (inform .eq. 7) go to 970

      lenlv_total = lenlv(1) + lenlv(2) + lenlv(3)
      do j =1, lenlv_total
         li(j) = indc(lena-lenL0-j+1)
         lj(j) = indr(lena-lenL0-j+1)
         lv(j) = a(lena-lenL0-j+1)
      end do
      go to 990

      ! Not enough storage.

  970 inform = 7
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1700) lena
      go to 990

      ! irep  is out of range.

  980 inform = 8
      if (nout. gt. 0  .and.  lprint .ge. 0)
     &     write(nout, 1800) m, n, irep

      ! Exit.

  990 luparm(2)  = lprint
      luparm(10) = inform
      return

 1700 format(/ ' lu8rpr  error...  Insufficient storage.',
     $         '    lena =', i8)
 1800 format(/ ' lu8rpr  error...  irep  is out of range.',
     $         '    m =', i8, '    n =', i8, '    irep =', i8)

      end ! subroutine lu8rpr

max_stars_repo_name: Rioghasarig/sr-cur

max_issues_repo_name: Rioghasarig/sr-cur

max_forks_repo_name: Rioghasarig/sr-cur


text:
program crm

!       Main module.

use vars
use hbuffer
use microphysics
use sgs
use tracers
use movies, only: init_movies
use held_suarez_mod
use python_caller
implicit none

integer k, icyc, nn, nstatsteps
double precision cputime, oldtime, init_time, elapsed_time !bloss wallclocktime
double precision usrtime, systime
integer itmp1(1), oldstep
!-------------------------------------------------------------------
! determine the rank of the current task and of the neighbour's ranks

call task_init()
!------------------------------------------------------------------
! print time, version, etc

if(masterproc) call header()
!------------------------------------------------------------------
! Initialize timing library.  2nd arg 0 means disable, 1 means enable

   call t_setoptionf (1, 0)
   call t_initializef ()

   call t_startf ('total')
   call t_startf ('initialize')
!------------------------------------------------------------------
! Get initial time of job

call t_stampf(init_time,usrtime,systime)
!------------------------------------------------------------------

call init()                    ! initialize some statistics arrays
call setparm()                 ! set all parameters and constants
if (dopython) call initialize_python_caller()

!------------------------------------------------------------------
! Initialize or restart from the save-dataset:

if(nrestart.eq.0) then
   day=day0
   call setgrid() ! initialize vertical grid structure
   call setdata() ! initialize all variables
elseif(nrestart.eq.1) then
   call read_all()
   call setgrid() ! initialize vertical grid structure
   call diagnose()
   call sgs_init()
   call micro_init()  !initialize microphysics
elseif(nrestart.eq.2) then  ! branch run
   call read_all()
   call setgrid() ! initialize vertical grid structure
   call diagnose()
   call setparm() ! overwrite the parameters
   call sgs_init()
   call micro_init()  !initialize microphysics
   nstep = 0
   day0 = day
else
   print *,'Error: confused by value of NRESTART'
   call task_abort()
endif

call init_movies()
call stat_2Dinit(1) ! argument of 1 means storage terms in stats are reset
call tracers_init() ! initialize tracers
call setforcing()

if(masterproc) call printout()
!------------------------------------------------------------------
!  Initialize statistics buffer:

call hbuf_init()

!------------------------------------------------------------------
nstatis = nstat/nstatfrq
nstat = nstatis * nstatfrq
nstatsteps = 0
call t_stopf ('initialize')
!------------------------------------------------------------------
!   Main time loop
!------------------------------------------------------------------
call t_stampf(oldtime,usrtime,systime)
oldstep = nstep

do while(nstep.lt.nstop.and.nelapse.gt.0)
  nstep = nstep + 1
  time = time + dt
  day = day0 + nstep*dt/86400.
  nelapse = nelapse - 1
!------------------------------------------------------------------
!  Check if the dynamical time step should be decreased
!  to handle the cases when the flow being locally linearly unstable
!------------------------------------------------------------------

  ncycle = 1

  call kurant()

  total_water_before = total_water()
  total_water_evap = 0.
  total_water_prec = 0.
  total_water_ls = 0.

  do icyc=1,ncycle

     icycle = icyc
     dtn = dt/ncycle
     dt3(na) = dtn
     dtfactor = dtn/dt

     if(mod(nstep,nstatis).eq.0.and.icycle.eq.ncycle) then
        nstatsteps = nstatsteps + 1
        dostatis = .true.
        if(masterproc) print *,'Collecting statistics...'
     else
        dostatis = .false.
     endif

     !bloss:make special statistics flag for radiation,since it's only updated at icycle==1.
     dostatisrad = .false.
     if(mod(nstep,nstatis).eq.0.and.icycle.eq.1) dostatisrad = .true.

!---------------------------------------------
!  	the Adams-Bashforth scheme in time

     call abcoefs()

!---------------------------------------------
!  	initialize stuff:

     call zero()

      if (doheldsuarez) call held_suarez_driver()
!-----------------------------------------------------------
!       Buoyancy term:

     call buoyancy()

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

     total_water_ls =  total_water_ls - total_water()

!------------------------------------------------------------
!       Large-scale and surface forcing:

     call forcing()

!----------------------------------------------------------
!       Nadging to sounding:

     call nudging()

!----------------------------------------------------------
!   	spange-layer damping near the upper boundary:

     if(dodamping) call damping()

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

     total_water_ls =  total_water_ls + total_water()

!---------------------------------------------------------
!   Ice fall-out

      if(docloud) then
          call ice_fall()
      end if

!----------------------------------------------------------
!     Update scalar boundaries after large-scale processes:

     call boundaries(3)

!---------------------------------------------------------
!     Update boundaries for velocities:

      call boundaries(0)

!-----------------------------------------------
!     surface fluxes:

     if(dosurface) call surface()
     call micro_flux() !bloss: moved here for consistency of isotopic fluxes

!-----------------------------------------------------------
!  SGS physics:

     if (dosgs) call sgs_proc()

!----------------------------------------------------------
!     Fill boundaries for SGS diagnostic fields:

     call boundaries(4)
!-----------------------------------------------
!       advection of momentum:

     call advect_mom()

!----------------------------------------------------------
!	SGS effects on momentum:
     if (dosgs) call boundaries(5)
     if(dosgs) call sgs_mom()

!-----------------------------------------------------------
!       Coriolis force:

     if(docoriolis) call coriolis()

!---------------------------------------------------------
!       compute rhs of the Poisson equation and solve it for pressure.

     call pressure()

!---------------------------------------------------------
!       find velocity field at n+1/2 timestep needed for advection of scalars:
!  Note that at the end of the call, the velocities are in nondimensional form.

     call adams()

!----------------------------------------------------------
!     Update boundaries for all prognostic scalar fields for advection:

     call boundaries(2)

!---------------------------------------------------------
!      advection of scalars :

     call advect_all_scalars()

!-----------------------------------------------------------
!    Convert velocity back from nondimensional form:

      call uvw()

!----------------------------------------------------------
!     Update boundaries for scalars to prepare for SGS effects:

     call boundaries(3)

!---------------------------------------------------------
!      SGS effects on scalars :

     if (dosgs) call sgs_scalars()

!-----------------------------------------------------------
!       Handle upper boundary for scalars

     if(doupperbound) call upperbound()

!-----------------------------------------------------------
!       Cloud condensation/evaporation and precipitation processes:

      if (docloud.or.dosmoke) call micro_proc()

!----------------------------------------------------------
!  Tracers' physics:

      call tracers_physics()

!-----------------------------------------------------------
!	Radiation

      if(dolongwave.or.doshortwave) then
         ! run on first step to initialize solin
         call radiation()
      end if

!--------------
! Neural network

      if (dopython .and. (mod(nstep-1, npython) == 0)) call state_to_python(dtn)
      if (usepython) then
         call apply_nn_forcings(dtn)
      end if

!-----------------------------------------------------------
!    Compute diagnostic fields:

      call diagnose()

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

! Rotate the dynamic tendency arrays for Adams-bashforth scheme:

      nn=na
      na=nc
      nc=nb
      nb=nn
      ! call python

   end do ! icycle

  total_water_after = total_water()

!----------------------------------------------------------
!  collect statistics, write save-file, etc.

   call stepout(nstatsteps)

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

   ! bloss:  only stop when stat and 2D fields are output.
   if ( (nstep.lt.nstop).and.(mod(nstep,nstat).eq.0) .AND. &
        ( (.NOT.save2Davg) .OR. (mod(nstep,nsave2D).eq.0) ) ) then

     if(masterproc) then
       call t_stampf(cputime,usrtime,systime)
       write(*,999) cputime-oldtime, float(nstep-oldstep)*dt/3600.
999    format('CPU TIME = ',f12.4, ' OVER ', f8.2,' MODEL HOURS')

       elapsed_time = cputime-init_time
       if ((elapsed_time+2.*(cputime-oldtime)).gt.60.*float(nelapsemin)) then
         nelapse=0 ! Job will stop when nelapse=0
         write(*,*) 'Job stopping -- nelapsemin reached'
       end if

       oldtime = cputime
       oldstep = nstep
     end if

      if(dompi) then
        if(masterproc) itmp1 = nelapse
        call task_bcast_integer(0,itmp1,1)
        nelapse = itmp1(1)
      end if
   end if

end do ! main loop

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

   call t_stopf('total')
   if(masterproc) call t_prf(rank)

if(masterproc) write(*,*) 'Finished with SAM, exiting...'
call task_stop()

end program crm

max_stars_repo_name: kbren/uwnet

max_issues_repo_name: kbren/uwnet

max_forks_repo_name: kbren/uwnet


text:
*
*  Program enddo - demonstrate how the END DO FORTRAN extension
*  works.
*
*  Written by: Charles Severance - 10Mar92
*
      integer I
*
*     Using the ENDDO it is not necesssary to specify a statement label
*     on the DO statement.  This makes code look much cleaner.
*
      PRINT *,'Starting a DO LOOP...'
      DO I=1,10
       PRINT *,'I=',I
      ENDDO
      PRINT *,'Done with the DO LOOP'
*
*     It is also possible to nest DO loops using the END DO statement.
*     FORTRAN matches the DO and END DO from the inside out.
*
      PRINT *,'Starting the second DO LOOP'
      DO I = 1,5,2
       PRINT *,'I=',I
       DO J=50,70,10
         PRINT *,I,J
       ENDDO
      ENDDO
      PRINT *,'Done with the second DO LOOP'
      END

ext: f


text:
        subroutine idd_poweroftwo(m,l,n)
c
c       computes l = floor(log_2(m)) and n = 2**l.
c
c       input:
c       m -- integer whose log_2 is to be taken
c
c       output:
c       l -- floor(log_2(m))
c       n -- 2**l
c
        implicit none
        integer l,m,n
c
c
        l = 0
        n = 1
c
 1000   continue
          l = l+1
          n = n*2
        if(n .le. m) goto 1000
c
        l = l-1
        n = n/2
c
c
        return
        end


      subroutine tag_denerror(tag,tagl1,tagl2,tagl3,tagh1,tagh2,tagh3, &
                              set,clear, &
                              den,denl1,denl2,denl3,denh1,denh2,denh3, &
                              lo,hi,nd,domlo,domhi, &
                              delta,xlo,problo,time,level)
      use amrex_fort_module, only : rt => amrex_real
      use probdata_module
      use turbforce_module, only : forcing_time_scale_max, stop_forcing
      implicit none

      integer set, clear, nd, level
      integer tagl1,tagl2,tagl3,tagh1,tagh2,tagh3
      integer denl1,denl2,denl3,denh1,denh2,denh3
      integer lo(3), hi(3), domlo(3), domhi(3)
      integer tag(tagl1:tagh1,tagl2:tagh2,tagl3:tagh3)
      real(rt) den(denl1:denh1,denl2:denh2,denl3:denh3,nd)
      real(rt) delta(3), xlo(3), problo(3), time

      real(rt) :: ax,ay,az
      integer          :: n_lo, n_hi
      integer          :: i, j, k, domhalf
      domhalf = (domhi(1)-domlo(1))/2
      n_lo = domlo(1) + domhalf - domhalf/(2**(level+1))
      n_hi = domlo(1) + domhalf + domhalf/(2**(level+1))

      ! use this to start at t = 0.25
      !if (time .ge. 0.125*(1 + 0.5*level)*stop_forcing*forcing_time_scale_max) then
      ! use this to test with a = 1.0
      if (time .ge. (1 + 0.5*level)*stop_forcing*forcing_time_scale_max) then
      ! use this to test with a = 0.5
      !if (time .ge. 0.5*(1 + 0.5*level)*stop_forcing*forcing_time_scale_max) then
         !$OMP parallel do private(k,j,i)
         do k = lo(3), hi(3)
            do j = lo(2), hi(2)
               do i = lo(1), hi(1)
                  if (i.ge.n_lo .and. i.lt.n_hi .and. &
                      j.ge.n_lo .and. j.lt.n_hi .and. &
                      k.ge.n_lo .and. k.lt.n_hi) then
                     tag(i,j,k) = set
                  endif
               enddo
            enddo
         enddo
      endif

!     Tag on regions of high density
!     if (level .lt. max_denerr_lev) then
!        do k = lo(3), hi(3)
!           do j = lo(2), hi(2)
!              do i = lo(1), hi(1)
!                if (den(i,j,k,1) .ge. denerr) then
!                    tag(i,j,k) = set
!                 endif
!              enddo
!           enddo
!        enddo
!     endif

!     Tag on regions of high density gradient
!     if (level .lt. max_dengrad_lev) then
!        do k = lo(3), hi(3)
!           do j = lo(2), hi(2)
!              do i = lo(1), hi(1)
!                 ax = ABS(den(i+1,j,k,1) - den(i,j,k,1))
!                 ay = ABS(den(i,j+1,k,1) - den(i,j,k,1))
!                 az = ABS(den(i,j,k+1,1) - den(i,j,k,1))
!                 ax = MAX(ax,ABS(den(i,j,k,1) - den(i-1,j,k,1)))
!                 ay = MAX(ay,ABS(den(i,j,k,1) - den(i,j-1,k,1)))
!                 az = MAX(az,ABS(den(i,j,k,1) - den(i,j,k-1,1)))
!                 if ( MAX(ax,ay,az) .ge. dengrad) then
!                    tag(i,j,k) = set
!                 endif
!              enddo
!           enddo
!        enddo
!     endif

      end subroutine tag_denerror


text:
      SUBROUTINE ZLARFY( UPLO, N, V, INCV, TAU, C, LDC, WORK )
*
*  -- LAPACK auxiliary test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      CHARACTER          UPLO
      INTEGER            INCV, LDC, N
      COMPLEX*16         TAU
*     ..
*     .. Array Arguments ..
      COMPLEX*16         C( LDC, * ), V( * ), WORK( * )
*     ..
*
*  Purpose
*  =======
*
*  ZLARFY applies an elementary reflector, or Householder matrix, H,
*  to an n x n Hermitian matrix C, from both the left and the right.
*
*  H is represented in the form
*
*     H = I - tau * v * v'
*
*  where  tau  is a scalar and  v  is a vector.
*
*  If  tau  is  zero, then  H  is taken to be the unit matrix.
*
*  Arguments
*  =========
*
*  UPLO    (input) CHARACTER*1
*          Specifies whether the upper or lower triangular part of the
*          Hermitian matrix C is stored.
*          = 'U':  Upper triangle
*          = 'L':  Lower triangle
*
*  N       (input) INTEGER
*          The number of rows and columns of the matrix C.  N >= 0.
*
*  V       (input) COMPLEX*16 array, dimension
*                  (1 + (N-1)*abs(INCV))
*          The vector v as described above.
*
*  INCV    (input) INTEGER
*          The increment between successive elements of v.  INCV must
*          not be zero.
*
*  TAU     (input) COMPLEX*16
*          The value tau as described above.
*
*  C       (input/output) COMPLEX*16 array, dimension (LDC, N)
*          On entry, the matrix C.
*          On exit, C is overwritten by H * C * H'.
*
*  LDC     (input) INTEGER
*          The leading dimension of the array C.  LDC >= max( 1, N ).
*
*  WORK    (workspace) COMPLEX*16 array, dimension (N)
*
*  =====================================================================
*
*     .. Parameters ..
      COMPLEX*16         ONE, ZERO, HALF
      PARAMETER          ( ONE = ( 1.0D+0, 0.0D+0 ),
     $                   ZERO = ( 0.0D+0, 0.0D+0 ),
     $                   HALF = ( 0.5D+0, 0.0D+0 ) )
*     ..
*     .. Local Scalars ..
      COMPLEX*16         ALPHA, tmp
      integer            i, ix, iy
*     ..
*     .. External Subroutines ..
      EXTERNAL           ZAXPY, ZHEMV, ZHER2
*     ..
*     .. External Functions ..
      COMPLEX*16         ZDOTC
      EXTERNAL           ZDOTC
*     ..
*     .. Executable Statements ..
*
      IF( TAU.EQ.ZERO )
     $   RETURN
*
*     Form  w:= C * v
*
      CALL ZHEMV( UPLO, N, ONE, C, LDC, V, INCV, ZERO, WORK, 1 )
*
**    zdotc crashes on some systems (MKL, MacOS); just inline it. -mgates
      tmp = ZERO
      ix = 1
      iy = 1
      do i = 1, n
         tmp = tmp + conjg(work(ix)) * v(iy)
         ix = ix + 1
         iy = iy + incv
      end do
      ALPHA = -HALF*TAU*tmp
**    ALPHA = -HALF*TAU*ZDOTC( N, WORK, 1, V, INCV )
      CALL ZAXPY( N, ALPHA, V, INCV, WORK, 1 )
*
*     C := C - v * w' - w * v'
*
      CALL ZHER2( UPLO, N, -TAU, V, INCV, WORK, 1, C, LDC )
*
      RETURN
*
*     End of ZLARFY
*
      END

max_stars_repo_name: maxhutch/magma

max_issues_repo_name: maxhutch/magma

max_forks_repo_name: maxhutch/magma


text:
module blocks

use cinter, only: mvaddch, endwin
use shapes, only: Piece
use fields, only: field

implicit none (type, external)
private
public :: game_over, draw_piece, freeze

contains

subroutine draw_piece(P)
class(piece), intent(in) :: P

integer :: i, j, x, y

! not concurrent since "mvaddch" remembers its position
do i = 1, P%Ny
  y = (i-1) + P%y
  if (y<0) cycle

  do j = 1, P%Nx
    x = (j-1) + P%x

    if (P%values(i, j) /= 0) then
      ! zero-indexed, so another -1
      call mvaddch(y=y-1, x=(P%x0-1) + x-1, &
                    ch=P%ch(P%values(i, j)))
    endif
  end do

end do
end subroutine draw_piece

subroutine freeze(F, P, NP)
! Called when a piece has hit another and freezes
class(piece), intent(inout) :: P, NP
class(field), intent(inout) :: F

integer :: i, y, ix, x, Nx, ixs
character(120):: buf

if(.not. P%landed) return

x = P%x
Nx = P%Nx

! not concurrent due to impure "game over"
ix = max(1, x)
ixs = min(F%W, ix + Nx - (ix-x) - 1)
do i = 1, P%Ny
  if (all(P%values(i, :) == 0)) cycle

  y = i-1 + P%y

  if (y <= 1)  then
    write(buf,'(A12,I3,A3,I3)') 'freeze: x=',x,'y=',y
    call game_over(F, P, buf)
  endif

  where(F%screen(y, ix:ixs) == 1 .or. P%values(i,ix-x+1:ixs-x+1) ==1)
    F%screen(y, ix:ixs) = 1
  endwhere
end do

call F%clear_lines()
call P%spawn_block(F, NP)
end subroutine freeze

subroutine game_over(F, P, msg)
class(field), intent(in) :: F
class(piece), intent(in), optional :: P
character(*), intent(in), optional :: msg
integer :: i

call endwin()

do i = 1,size(F%screen,1)
  print '(100I1)',F%screen(i,:)
enddo

if(present(P)) print *, P%why
if(present(msg)) print *, msg

print *,' '
print *, 'Level:', F%level
Print *, 'Score:', F%score
print *, 'Number of Blocks:', F%Nblock
print *, 'Number of Lines Cleared:', F%Ncleared
print *, 'Block Sequence: ', F%blockseq(:F%Nblock)

if (F%debug) then
  write(F%udbg,*) 'Block Sequence: ', F%blockseq(:F%Nblock)
  close(F%udbg)
endif

stop 'Goodbye from Tetran'

end subroutine game_over

end module blocks

max_stars_repo_name: fortran-gaming/tetran

max_issues_repo_name: fortran-gaming/blocktran

max_forks_repo_name: fortran-gaming/blocktran


text:
program verify_bounds
    implicit none
    integer, parameter :: NR_CASES = 10000
    integer :: i, start_val, end_val, step

    do i = 1, NR_CASES
        call generate_bounds(start_val, end_val, step)
        call test_do_loop(start_val, end_val, step)
    end do

contains

    subroutine test_do_loop(start_val, end_val, step)
        implicit none
        integer, value :: start_val, end_val, step
        integer :: i, expected_i, last_val, expected_last_val, &
                   nr_iters, expected_nr_iters

        nr_iters = 0
        last_val = start_val
        do i = start_val, end_val, step
            nr_iters = nr_iters + 1
            last_val = i
        end do
        expected_nr_iters = (end_val - start_val + step)/step
        if (nr_iters /= expected_nr_iters) then
            print '(5(A, I0))', 'start = ', start_val, &
                                ', end = ', end_val, &
                                ' step = ', step, &
                                ', expected nr. iters = ', expected_nr_iters, &
                                ', nr. iters = ', nr_iters
        end if
        if (nr_iters /= 0) then
            expected_last_val = start_val + (nr_iters - 1)*step
            if (last_val /= expected_last_val) then
                print '(5(A, I0))', 'start = ', start_val, &
                                    ', end = ', end_val, &
                                    ' step = ', step, &
                                    ', expected last = ', expected_last_val, &
                                    ', last = ', last_val
            end if
        end if
        expected_i = start_val + expected_nr_iters*step
        if (i /= expected_i) then
            print '(5(A, I0))', 'start = ', start_val, &
                                ', end = ', end_val, &
                                ' step = ', step, &
                                ', expected i = ', expected_i, &
                                ', i = ', i
        end if
    end subroutine test_do_loop

    subroutine generate_bounds(start_val, end_val, step)
        implicit none
        integer, intent(out) :: start_val, end_val, step

        step = 0
        do while (step == 0)
            step = random_int(-100, 100)
        end do
        if (step >= 0) then
            do while (.true.)
                start_val = random_int(-100, 100)
                end_val = random_int(-100, 100)
                if (start_val <= end_val) exit
            end do
        else
            do while (.true.)
                start_val = random_int(-100, 100)
                end_val = random_int(-100, 100)
                if (start_val >= end_val) exit
            end do
        end if
    end subroutine generate_bounds

    function random_int(from_val, to_val) result(r_int)
        implicit none
        integer ,value :: from_val, to_val
        integer ::  r_int
        real :: r

        call random_number(r)
        r_int = int(from_val + r*(to_val - from_val))
    end function random_int

end program verify_bounds

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC