text:
      subroutine termin
!***********************************************************************
!  Copyright, 1993, 2004,  The  Regents of the University of California.
!  This program was prepared by the Regents of the University of
!  California at Los Alamos National Laboratory (the University) under
!  contract No. W-7405-ENG-36 with the U.S. Department of Energy (DOE).
!  All rights in the program are reserved by the DOE and the University.
!  Permission is granted to the public to copy and use this software
!  without charge, provided that this Notice and any statement of
!  authorship are reproduced on all copies. Neither the U.S. Government
!  nor the University makes any warranty, express or implied, or
!  assumes any liability or responsibility for the use of this software.
C***********************************************************************
CD1
CD1 PURPOSE
CD1
CD1 Read I/O file names from terminal, set unit numbers and open files.
CD1
C***********************************************************************
CD2
CD2 REVISION HISTORY
CD2
CD2 Revision                    ECD
CD2 Date         Programmer     Number  Comments
CD2
CD2 22-NOV-93    Z. Dash        22      Initial implementation.
CD2
CD2 $Log:   /pvcs.config/fehm90/src/termin.f_a  $
!D2
!D2    Rev 2.5   06 Jan 2004 10:44:18   pvcs
!D2 FEHM Version 2.21, STN 10086-2.21-00, Qualified October 2003
!D2
!D2    Rev 2.4   29 Jan 2003 09:20:56   pvcs
!D2 FEHM Version 2.20, STN 10086-2.20-00
!D2
!D2    Rev 2.3   14 Nov 2001 13:15:58   pvcs
!D2 FEHM Version 2.12, STN 10086-2.12-00
!D2
!D2    Rev 2.2   06 Jun 2001 13:28:02   pvcs
!D2 FEHM Version 2.11, STN 10086-2.11-00
!D2
!D2    Rev 2.1   30 Nov 2000 12:12:28   pvcs
!D2 FEHM Version 2.10, STN 10086-2.10-00
!D2
!D2    Rev 2.0   Fri May 07 14:47:22 1999   pvcs
!D2 FEHM Version 2.0, SC-194 (Fortran 90)
CD2
CD2    Rev 1.7   Fri Feb 02 12:45:20 1996   hend
CD2 Updated Requirements Traceability
CD2
CD2    Rev 1.6   06/20/94 11:09:10   zvd
CD2
CD2
CD2    Rev 1.5   03/18/94 15:55:58   gaz
CD2 Added solve_new and cleaned up memory management.
CD2
CD2    Rev 1.4   02/28/94 11:59:44   zvd
CD2 Edited Rev 1.3 comment to reflect actual update.
CD2
CD2    Rev 1.3   02/28/94 11:52:38   zvd
CD2 Corrected problem of unit assignment for coordinate and zone
CD2 input files.
CD2
CD2    Rev 1.2   02/17/94 10:17:20   zvd
CD2 Fixed bug when opening alternate zone or coordinate data files.
CD2
CD2    Rev 1.1   02/14/94 11:37:20   zvd
CD2 Modified so terminal I/O is compatible with control file input.
CD2 Added single character response acceptance for terminal queries.
CD2
CD2    Rev 1.0   01/20/94 10:28:32   pvcs
CD2 original version in process of being certified
CD2
C***********************************************************************
CD3
CD3 INTERFACES
CD3
CD3 Formal Calling Parameters
CD3
CD3   None
CD3
CD3 Interface Tables
CD3
CD3   None
CD3
CD3 Files
CD3
CD3   See name descriptors under Global Variables
CD3
C***********************************************************************
CD4
CD4 GLOBAL OBJECTS
CD4
CD4 Global Constants
CD4
CD4   None
CD4
CD4 Global Types
CD4
CD4   None
CD4
CD4 Global Variables
CD4
CD4                            COMMON
CD4   Identifier      Type     Block  Description
CD4
CD4   cform           CHAR ARY faax   I/O file formats
CD4   cstats          CHAR ARY faax   I/O file status
CD4   ex              LOGICAL  faay   Logical for existence check (T/F)
CD4   iowork          CHAR ARY faax   File use/type designator
CD4   itert           INT      faai   Intermediate iteration counter
CD4   nmfil           CHAR ARY faax   I/O file names:
CD4                                   1 - Control file name
CD4                                   2 - Main input file name
CD4                                   3 - Coordinate input file name
CD4                                   4 - Zone input file name
CD4                                   5 - Main output file name
CD4                                   6 - Restart input file name
CD4                                   7 - Restart output file name
CD4                                   8 - Simulation history output file name
CD4                                   9 - Solute history output file name
CD4                                   10 - Contour plot output file name
CD4                                   11 - Dual porosity or dpdp contour plot
CD4                                          output file name
CD4                                   12 - Coefficient storage output file name
CD4                                   13 - Input check output file name
CD4                                   14 - Error output file name
CD4   nmfily          CHAR ARY faax   File messages
CD4   nmmax           INT      faay   Number of I/O files
CD4   nufilb          CHAR ARY faax   File unit numbers
CD4   suffix          CHAR ARY faax   Default file suffixes
CD4
CD4 Global Subprograms
CD4
CD4   Identifier      Type     Description
CD4
CD4   prefix                   Determine prefix of input file name
CD4   setunits                 Set file unit numbers
CD4
C***********************************************************************
CD5
CD5 LOCAL IDENTIFIERS
CD5
CD5 Local Constants
CD5
CD5   None
CD5
CD5 Local Types
CD5
CD5   None
CD5
CD5 Local variables
CD5
CD5   Identifier      Type     Description
CD5
CD5   i               INT      Loop index
CD5   namfil          CHAR     Input file name
CD5   numfil          INT ARY  Unit numbers for I/O files
CD5   yesno           CHAR     String used for input yes or no response
CD5
CD5 Local Subprograms
CD5
CD5   None
CD5
C***********************************************************************
CD6
CD6 FUNCTIONAL DESCRIPTION
CD6
C***********************************************************************
CD7
CD7 ASSUMPTIONS AND LIMITATIONS
CD7
CD7 None
CD7
C***********************************************************************
CD8
CD8 SPECIAL COMMENTS
CD8
CD8  Requirements from SDN: 10086-RD-2.20-00
CD8    SOFTWARE REQUIREMENTS DOCUMENT (RD) for the
CD8    FEHM Application Version 2.20
CD8
C***********************************************************************
CD9
CD9 REQUIREMENTS TRACEABILITY
CD9
CD9 2.6 Provide Input/Output Data Files
CD9 3.0 INPUT AND OUTPUT REQUIREMENTS
CD9
C***********************************************************************
CDA
CDA REFERENCES
CDA
CDA None
CDA
C***********************************************************************
CPS
CPS PSEUDOCODE
CPS
CPS BEGIN termin
CPS
CPS   initialize yesno flag to 'n'
CPS
CPS   FOR each I/O file (except error file)
CPS
CPS       set file unit number
CPS
CPS       IF this isn't the control or input file and the file name
CPS        prefix is different from the input file name prefix
CPS          IF this is the coordinate or zone file
CPS             set default file name to input file name
CPS          ELSE
CPS             set default file name using on input file name prefix and
CPS              default suffix
CPS          END IF
CPS       END IF
CPS
CPS       IF the file name should be set (flag is 'n')
CPS
CPS          initialize input file string to blank
CPS          determine file name or if using default or if not using
CPS
CPS          IF default input is to be used
CPS             set input file string to default
CPS          END IF
CPS
CPS          IF not using file
CPS
CPS             set input file string to 'not using'
CPS             set file unit number to 0
CPS
CPS          ELSE IF file status is old
CPS
CPS             determine if file exists
CPS             IF the file doesn't exist
CPS                IF this is the coefficient storage file
CPS                   set file status to new
CPS                ELSE
CPS                   write message to screen and request another name
CPS                   return to file name query
CPS                END IF
CPS             END IF
CPS
CPS             IF this is the control file
CPS
CPS                set control file unit number and name
CPS                return to calling routine to invoke control file I/O
CPS
CPS             ELSE IF this is the input file
CPS
CPS                set input unit number and name
CPS
CPS                call file_prefix to determine the input file name prefix
CPS                IF prefix is longer than 94 characters
CPS                   warn user that prefix was truncated to 94
CPS                    charcaters and set prefix length tho 94
CPS                END IF
CPS
CPS                determine if all file names should use the input file prefix
CPS                IF flag isn't set to 'n'
CPS                   set flag to 'y'
CPS                END IF
CPS
CPS             END IF
CPS
CPS          END IF

CPS       ELSE
CPS
CPS          IF this is the coordinate or zone input file
CPS             set file unit and name equal to the input file
CPS          ELSE
CPS             use the default file name
CPS          END IF
CPS
CPS          IF this is the initialization file
CPS             determine if file existss
CPS             IF the file does not exist
CPS                set file unit to 0
CPS             END IF
CPS          ELSE IF this is the coefficient storage file
CPS             determine if file existss
CPS             IF the file does not exist
CPS                set file status to new
CPS             END IF
CPS          END IF
CPS
CPS       END IF
CPS
CPS       IF the file unit number is 0
CPS
CPS          IF this is the input file
CPS             set unit to tty (5) and name to 'terminal'
CPS          ELSE if this is the coordiante or zone input file
CPS             set file unit and name equal to the input file
CPS          ELSE if this is the output file
CPS             set unit to tty (6) and name to 'terminal'
CPS          ELSE
CPS             set input file name string to 'not using'
CPS          END IF
CPS
CPS       ELSE
CPS
CPS          IF this isn't the save file or the coordinate and zone files
CPS           and they aren't the same as the input file
CPS             open and rewind the file
CPS          END IF
CPS
CPS       END IF
CPS
CPS       set file name to input file string
CPS
CPS   END FOR
CPS
CPS   call setunits to assign file unit numbers
CPS
CPS END termin
CPS
C***********************************************************************

      use comai
      use comxi
      implicit none

      integer i, numfil(13)
      character*1 yesno
      character*100 namfil

      yesno = 'n'
C     Make sure coefficient storage file starts out with status of old
      cstats(12) = 'old    '

      do i = 1, nmmax - 1

C     **** Append file suffix to prefix from input file name ****
         if (i .gt. 2 .and.
     *        nmfil(2)(1:itert) .ne. nmfil(i)(1:itert)) then
            nmfil(i) = blank
            if (i .eq. 3 .or. i .eq. 4) then
               nmfil(i) = nmfil(2)
            else
               nmfil(i)(1:itert) = nmfil(2)(1:itert)
               nmfil(i)(itert+1:itert+5) = suffix(i)
            end if
         endif

         numfil(i) = nufilb(i)

         if (yesno .eq. 'n') then
            namfil = blank
 10         write(6, 6000)  iowork(i), nmfil(i)
 6000       format(/, 1x, 'Enter name for ', a6, 1x,
     *           '-- default file name: ', a100, /, 1x,
     *           '[(name/na or not using), RETURN = DEFAULT]')

            read (5, '(a100)')  namfil
            if (namfil(1:4) .eq. '    ') then
               namfil = nmfil(i)
            end if
            if (namfil(1:4) .eq. 'na  ' .or.
     *           namfil(1:4) .eq. 'not ') then
               namfil = nmfily(3)
               numfil(i) = 0
            else if (cstats(i) .eq. 'old    ')  then
               inquire (file = namfil, exist = ex)
               if (.not. ex) then
                  if (iowork(i) .eq. 'isstor') then
                     cstats(i) = 'new    '
                  else
                     write(6,*) iowork(i), ' Input file does not exist,'
     *                    , ' enter another'
                     go to 10
                  endif
               endif
               if (iowork(i) .eq. 'iocntl') then
                  iocntl = numfil(i)
                  nmfil(i) = namfil
                  return
               else if (iowork(i) .eq. 'inpt  ') then
                  inpt = numfil(i)
                  nmfil(i) = namfil
                  call file_prefix(namfil, itert)
                  if (itert .gt. 94) then
                     itert = 94
                     write (6, *) ' **** Input file name prefix was',
     *                    ' truncated to 94 characters ****'
                  end if
                  write (6, *)
                  write(6, *) ' Do you want all file names of the form '
     *                 , namfil(1:itert), '.*',' ? [(y/n), RETURN = y]',
     *                 ' *** Note: If "y" incoor and inzone will equal',
     *                 ' inpt ***'

                  read(5, '(a1)') yesno(1:1)
                  if (yesno(1:1) .ne. 'n') yesno = 'y'
               end if
            end if
         else
C     **** Set incoor and inzone to inpt ****
            if (iowork(i) .eq. 'incoor' .or.
     *           iowork(i) .eq. 'inzone') then
               namfil = nmfil(2)
               nmfil(i) = nmfil(2)
               numfil(i) = numfil(2)
            else
               namfil = nmfil(i)
            end if
C     **** Set iread to 'not using' if file not found ****
            if (iowork(i) .eq. 'iread ') then
               inquire (file = namfil, exist = ex)
               if (.not. ex) then
                  numfil(i) = 0
               end if
C     **** Set isstor to status new if file not found ****
            else if (iowork(i) .eq. 'isstor') then
               inquire (file = namfil, exist = ex)
               if (.not. ex) then
                  cstats(i) = 'new    '
               end if
            end if
         end if

         if (numfil(i) .eq. 0)  then
            if (iowork(i) .eq. 'inpt  ')  then
               numfil(i) = 5
               namfil = nmfily(1)
            else if (iowork(i) .eq. 'incoor' .or.
     *              iowork(i) .eq. 'inzone')  then
               numfil(i) = numfil(2)
               namfil = nmfil(2)
            else if (iowork(i) .eq. 'iout  ')  then
               numfil(i) = 6
               namfil = nmfily(2)
            else
               namfil = nmfily(3)
            end if
         else  if (numfil(i) .eq. 5)  then
            namfil = nmfily(1)
         else if (numfil(i) .eq. 6)  then
            namfil = nmfily(2)
         else
            if (iowork(i) .ne. 'isave ' .and.
     *           .not. ((iowork(i) .eq. 'incoor' .or.
     *           iowork(i) .eq. 'inzone') .and.
     *           namfil .eq. nmfil(2))) then
               open (numfil(i), file = namfil, status = cstats(i),
     *              form = cform (i))
               rewind  numfil(i)
            else if ((iowork(i) .eq. 'incoor' .or.
     *              iowork(i) .eq. 'inzone') .and.
     *              nmfil(i) .eq. nmfil(2)) then
               numfil(i) = numfil(2)
            end if
         end if
         nmfil(i) = namfil
      end do

      call setunits (numfil)

      return
      end

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
Angelo Moreno has been director of the DHS Orchestra Davis High School Orchestra program since 2000 and current leads the Symphony Orchestra, the Chamber Orchestra, and the Baroque Ensemble.  He received both his bachelors and masters degrees in music from UC Davis, and received the Harmony in Our Lives Award in 2011.  He is also recipient of the 2013 http://calmusiced.com/index.php/cmea/99stateawardrecepients2 Richard Levin State Orchestra Educator Award and the Sacramento News & Reviews 2009 Music Educators Outstanding Achievement Award.  His achievements include developing the DHS Orchestra program from an enrollment of about 35 students when he first arrived to more than four times that size in recent years.  Other DHS Orchestra DHS Orchestra Program achievements under his direction can be found at their Davis Wiki DHS Orchestra article.  Moreno is also one of the directors of the https://sacramentoyouthsymphony.org/ Sacramento Youth Symphony.

Media

  Hudson, Jeff; 2011 Nov. 9, http://www.davisenterprise.com/arts/music/morenotobehonoredatharmonyinourlivesconcertssunday/ Moreno to be honored at Harmony in Our Lives concerts Sunday, Davis Enterprise, pg. A1.  Retrieved 25 January 2013.
  Kincaid, Jennifer; 2009 Feb. 19, http://www.newsreview.com/sacramento/musicalinspiration/content?oid914385 Musical Inspiration: 2009 Sacramento Music Educators Outstanding Achievement Award recipient Angelo Moreno, Sacramento News & Review. Retrieved 22 Jan. 2013.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c=======================================================================
c
c     subroutine ALLOCRB
c
c     Risk Budgeting allocation strategy
c     (Multi-Volatility-constrained)
c
c     Max[ w'*rho - k*w'*Q*w]
c      s.t.
c     sqrt[w'*cov_i*w] <= sigma(i)  (i-th risk budgeting constraint)
c     C*w <= b                      (general linear constraints)
c     Cinf <= w <= Csup             (lower/upper bounds)
c
c        w   : portfolio weights
c        Q   : covariance matrix
c        rho : assets performance
c
c-----------------------------------------------------------------------
      SUBROUTINE allocrb ( n, cov, Q, kappa, rho,
     &                     neq, nin, ccst, bcst,
     &                     cinf, csup, nbc, class, sigma,
     &                     iwork, dwork, wopt, lambda, info )
c-----------------------------------------------------------------------
c
c     INPUT
c       n      : portfolio size                                  integer
c       cov    : covariance matrix (n*n)                          double
c       Q      : matrix (n*n)                                     double
c       kappa  : risk aversion coefficient                        double
c       rho    : expected returns vector (n)                      double
c       neq    : number equality constraints                     integer
c       nin    : number inequality constraints                   integer
c       ccst   : matrix of constraints (n*(neq+nin))              double
c       bcst   : vector initial of constraints (neq+nin)          double
c       cinf   : lower bound (n)                                  double
c       csup   : upper bound (n)                                  double
c       nbc    : number of risk budgeting constraints (>=1)      integer
c       class  : block definition (n)                            integer
c       sigma  : volatilities budgets constraints (nbc)           double
c
c     WORKSPACE
c       iwork  : 2*nbc + 5*n + neq + 2*nin + 12                  integer
c       dwork  : nbc*(nbc+21)/2 +
c                n*(8*n+20+2*neq+2*nin)+4*neq+6*nin+14            double
c     OUTPUT
c       wopt   : optimal portfolio (n)                            double
c       lambda : optimal lambda (nbc)                             double
c       info   : diagnostic argument                             integer
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, info, neq, nin, nbc
      INTEGER class(*)
      DOUBLE PRECISION kappa
      DOUBLE PRECISION cov(*), Q(*), rho(*), cinf(*), csup(*),
     &                 ccst(n,*), bcst(*), wopt(*), sigma(*), lambda(*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER i, j, npk, totitr, totsim, ncst, infoBFGS, test,
     &        infot, izero
      INTEGER piwbf, pisim, pinfo, pin, pineq, pinin, piclass, piqp,
     &        pydual, pbinf, pbsup, pgrad, pwbf, pdsim,
     &        pdcov, pdrho, pdcinf, pdcsup, pdsigma, pdc, pdb,
     &        pdcov1, pdQ, pdkappa, pdqp, pdlagr, pdccs, pdbcs
      DOUBLE PRECISION funct, epsbfg, tmp, scal
      DOUBLE PRECISION INFINI, MINFINI, ZERO, EPS, ONE, COEF
      PARAMETER (INFINI=1.E+15,ZERO=0.E0,EPS=1.E-30,ONE=1E0,COEF=-0.5E0)
      PARAMETER (MINFINI = -1.E10)
c
c     external subroutines
      EXTERNAL simrb, simext, testsdp, OVTMCV, IVX, YV, bfgsmvol,
     &         checklinbox, initfeas, IVI
c
c     intrinsic functions
c
c-----------------------------------------------------------------------
c
c
c     initializations
      info = 0          ! diagnostic argument
      ncst = neq + nin  ! number total of constrains
      izero = 0
c
c     pointers for integer work space : iwork
c     ---------------------------------------
      piwbf = 1
c     piwbf : pointer for BFGSMVOL who needs (2*nbc + 1)
c     so (2*nbc + 1)
      pisim  = piwbf + (2*nbc + 1)
c     pisim : pointer for simulator SIMRB (5*n + neq + 2*nin + 11)
c     so (2*nbc + 5*n + neq + 2*nin + 12)
      pinfo = pisim
c     pinfo : pointer for info (errors of simul), so 1 more
      pin   = pinfo + ( 1 )
c     pin   : pointer for value of n, so 1 more
      pineq = pin + ( 1 )
c     pineq : pointer for value of neq, so 1 more
      pinin = pineq + ( 1 )
c     pinin : pointer for value of nin, so 1 more
      piclass = pinin + ( 1 )
c     piclass: pointer for value of class, so n more
      piqp = piclass + ( n )
c     piqp : pointer for SIMRB
c
c     Total size of dwork array = 2*nbc + 5*n + neq + 2*nin + 12
c
c     pointers for double precision work space  : dwork
c     -------------------------------------------------
      pydual = 1
c     pydual: pointer for initial value of ydual (nbc)
c     so (nbc)
      pbinf = pydual + ( nbc )
c     pbinf : pointer for inferior bounds binf (nbc)
c     so (2*nbc)
      pbsup = pbinf + ( nbc )
c     pbsup : pointer for superior bounds bsup (nbc)
c     so (3*nbc)
      pgrad = pbsup + ( nbc )
c     pgrad : pointer for gradient at solution (nbc)
c     so (4*nbc)
      pwbf = pgrad + ( nbc )
c     pwbf  : pointer for BFGSMVOL ( nbc*(nbc+11)/2 )
c     so ( nbc*(nbc+19)/2 )
      pdsim = pwbf + ( nbc*(nbc+11)/2 )
c     pdsim  : pointer for the simulator SIMRB
c             n*(7*n+20+2*neq+2*nin)+4*neq+6*nin+14+nbc
c     so
      pdcov   = pdsim
c     pdcov : pointer for cov. matrix (n*n)
      pdQ     = pdcov + ( n*n )
c     pdQ   : pointer for Q matrix (n*n)
      pdkappa = pdQ + ( n*n )
c     pdcov : pointer for kappa (1)
      pdrho   = pdkappa + ( 1 )
c     pdrho : pointer for expected returns (n)
      pdcinf  = pdrho + ( n )
c     pdcinf : pointer for lower bounds (n)
      pdcsup  = pdcinf + ( n )
c     pdcsup : pointer for upper bounds (n)
      pdsigma = pdcsup + ( n )
c     pdsigma : pointer for risk budget constraint (nbc)
      pdc     = pdsigma + ( nbc )
c     pdc : pointer for matrix constraint (n*(neq + nin))
      pdb     = pdc + ( n*(neq + nin) )
c     pdb : pointer for vector constraint (neq + nin)
      pdqp    = pdb + ( neq + nin )
c     pdqp : pointer for QP (SIMRB), (n*n + 6*n + 2*nin)
      pdlagr  = pdqp + ( n*n + 6*n + 2*nin )
c     pdlagr : pointer for QP parameter (2*n*n + 3*n + neq + nin)
      pdcov1  = pdlagr + ( 2*n*n + 3*n + neq + nin )
c     pdcov1  : temporary cov. matrix ( n*n )
      pdccs = pdcov1 + ( n*n )
c     needs n*(neq+nin+1)
      pdbcs = pdccs + n*(neq+nin+1)
c     needs neq+nin+1
c
c     Total size of dwork array = 4*nbc + ( nbc*(nbc+11)/2 )
c                               + n*(7*n+20+2*neq+2*nin)+4*neq+6*nin+14+nbc (SIMRB)
c                               + n*n + n*(neq+nin+1) + neq+nin+1
c           = nbc*(nbc+21)/2 + n*(8*n+21+3*neq+3*nin)+5*neq+7*nin+15
c
c-----------------------------------------------------------------------
c
c      TODO - test if classes are not all null
c
c-----------------------------------------------------
c     satisfiability conditions (well-defined problem)
c-----------------------------------------------------
c
c     matrices are SDP ?
      CALL testsdp (n, cov, iwork(pisim), dwork(pdsim), info)
      IF (info .NE. 0) THEN
         info = -108
         RETURN
      ENDIF

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'TESTSDP COV ', info
c      close(unit=1)

      CALL YV ( n*n, Q, dwork(pdQ) )
      IF (kappa .LT. EPS) CALL YV ( n*n, cov, dwork(pdQ) )
      CALL testsdp (n, dwork(pdQ), iwork(pisim), dwork(pdsim), info)
      IF (info .NE. 0) THEN
         info = -109
         RETURN
      ENDIF

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'TESTSDP Q ', info
c      close(unit=1)

c
c     problem initialization and feasibility check
      CALL initfeas ( n, cinf, csup, neq, nin, ccst, bcst,
     &                iwork(pisim), dwork(pdsim), wopt, info)
      IF (info .EQ. 1111) THEN
        info = -100
        RETURN
      ENDIF

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'INITFEAS 1 ', info
c      close(unit=1)

c
c------------------------------------------------------------
c     construction of work vectors (communication with simrb)
c------------------------------------------------------------
c
c     initialization
      CALL IVX ( nbc, dwork(pydual), ONE )   ! dual
      CALL IVX ( nbc, dwork(pbinf), ZERO )   ! lower bounds
      CALL IVX ( nbc, dwork(pbsup), INFINI ) ! upper bounds
c
c     construct of iwork (communication with simrb)
      iwork(pinfo) = 0       ! diagnostic argument
      iwork(pin)   = n       ! number of asset(s)
      iwork(pineq) = neq     ! number of eq. constraint(s)
      iwork(pinin) = nin     ! number of ineq. constraint(s)
      CALL YVI ( n, class, iwork(piclass) ) ! cov. matrix partition
c
c     construct dwork (communication with simrb)
      CALL YV(n*n, cov, dwork(pdcov))     ! cov. matrix
      CALL YV(n*n, Q, dwork(pdQ))         ! Q matrix
      CALL YV(n, rho, dwork(pdrho))       ! rho
      dwork(pdkappa) = kappa              ! kappa coef.
      CALL YV(n, cinf, dwork(pdcinf))     ! lower bounds
      CALL YV(n, csup, dwork(pdcsup))     ! upper bounds
      CALL YV(nbc, sigma, dwork(pdsigma)) ! risk budget(s)
c
c     case rho(1) = rho(2) = ... = rho(n)
      tmp = rho(1)
      test = 0
      DO i = 2,n
       IF (ABS(tmp - rho(i)) .GT. EPS) THEN
        test = 1
       ENDIF
      ENDDO
      IF (test .EQ. 0) THEN
        dwork(pdrho + 1) = dwork(pdrho + 1)
     &                     + 1.E-4*(MAX(dwork(pdrho + 1),1.E-6))
      ENDIF
c
c     case ncst <> 0 (neq=0, nin=0)
      IF (ncst .NE. 0) THEN
        CALL YM(n, ncst, ccst, dwork(pdc))  ! constraint matrix
        CALL YV(ncst, bcst, dwork(pdb))     ! constraint vector
      ENDIF
c
c-------------------------------------------------------------------
c     Non Linear optimization (BFGS)
c-------------------------------------------------------------------
      epsbfg = 1.E-10   ! precision stop test BFGS
      CALL bfgsmvol ( simrb, simext, nbc, dwork(pydual), epsbfg,
     &                dwork(pbinf), dwork(pbsup),
     &                iwork(piwbf), dwork(pwbf),
     &                funct, dwork(pgrad), totitr, totsim, info)

      CALL YV(nbc, dwork(pydual), lambda) ! ydual -> lambda
      infoBFGS = info ! BFGS diagnostic argument

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'BFGSMVOL ', infoBFGS
c      close(unit=1)

c
c-------------------------------------------------
c     Optimize with optimal dual solution (lambda)
c-------------------------------------------------
c
      CALL IMX ( n, n, dwork(pdcov), ZERO ) ! initialize pdcov
c
c     construct cov. matrix with optimal dual solution
      DO i = 1,nbc
         CALL IVI(n, iwork(piwbf), izero) ! initialize piwbf
c
c        dimension of block i
         npk = 0
         DO j = 1,n
            IF (class(j) .EQ. i) THEN
               npk = npk + 1
               iwork(piwbf + npk - 1) = j
            ENDIF
         ENDDO
         IF (npk .GT. 1) THEN
c
c           initialize temporary cov. matrix
            CALL IMX ( npk, npk, dwork(pdcov1), ZERO )
c
c           extract block i
            CALL YMCPI ( n, cov, npk, iwork(piwbf), dwork(pdcov1), info)
c
c           product lambda(i)*cov(i)
            CALL PMX (npk, npk, dwork(pdcov1), lambda(i), dwork(pdcov1))
c
c           sub-block -> cov. matrix
            CALL YMCPIR (npk, dwork(pdcov1), n, iwork(piwbf),
     &                   dwork(pdcov), info)
         ENDIF
      ENDDO
c
c     product kappa*Q -> pdcov2
      CALL PMX(n, n, Q, kappa, dwork(pdcov1))
c
c     kappa*Q + Sum(Cov_i)
      CALL SM(n, n, dwork(pdcov), dwork(pdcov1), dwork(pdcov))
c
c     linear part: l(i) = -0.5*rho(i)
      scal = -2.
      CALL PVX ( n, rho, COEF, dwork(pdrho) )
c
c     case rho(1) = rho(2) = ... = rho(n)
      tmp = rho(1)
      test = 0
      DO i = 2,n
       IF (ABS(tmp - rho(i)) .GT. EPS) THEN
        test = 1
       ENDIF
      ENDDO
      IF (test .EQ. 0) THEN
        dwork(pdrho + 1) = dwork(pdrho + 1)
     &                     + 1.E-4*(MAX(dwork(pdrho + 1),1.E-6))
      ENDIF
c
c     problem initialization and feasibility check
      CALL initfeas ( n, cinf, csup, neq, nin, ccst, bcst,
     &                iwork(piqp), dwork(pdqp), wopt, info)
      IF (info .EQ. 1111) THEN
        info = -100
        RETURN
      ENDIF

c     open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'INITFEAS 2', info
c      close(unit=1)

c
c     quadratic solver QP
      CALL qp ( n, dwork(pdcov), dwork(pdrho), neq, nin,
     &          dwork(pdc), dwork(pdb), cinf, csup,
     &          iwork(piqp), dwork(pdqp),
     &          dwork(pdlagr), wopt, info )

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'QUAPRO ', info
c      close(unit=1)

c
c----------------------------------------------------------------------
c     test if the ith-volatility budget constraint is too high/small
c     (budget constraint not attainable) -> info = 111/110
c----------------------------------------------------------------------
c
      IF (info .NE. 1001) THEN
        DO i = 1,nbc
            IF (ABS(lambda(i)) .LT. 1.E-15) THEN
                info = 111
            ENDIF
            IF (ABS(lambda(i)) .GT. 1.E+15) THEN
                info = 110
            ENDIF
        ENDDO
      ENDIF
c
c     checking non-emptiness of linear constraints
c      CALL checklinbox ( n, wopt,
c     &                   neq, nin, ccst, bcst, cinf, csup,
c     &                   dwork(pwbf), infot)
c      IF (infot .LT. 0) THEN
c        info = infot
c        RETURN
c      ENDIF

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'CHECKLINBOX ', info
c      close(unit=1)

c
c     construct cov. matrix with optimal dual solution
      DO i = 1,nbc
         DO j = 1,n
            iwork(piwbf + i - 1) = 0
         ENDDO
c
c        dimension of the i-th block
         npk = 0
         DO j = 1,n
            IF (class(j) .EQ. i) THEN
               npk = npk + 1
               iwork(piwbf + npk - 1) = j
            ENDIF
         ENDDO

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'TEST 1', info
c      close(unit=1)

         IF (npk .GT. 1) THEN
c
c           initialize temporary cov. matrix
            CALL IMX ( npk, npk, dwork(pdcov1), ZERO )
            CALL IMX ( n, n, dwork(pdcov), ZERO )

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'TEST 2', info
c      close(unit=1)

c
c           extracting the i-th block
            CALL YMCPI(n, cov, npk, iwork(piwbf), dwork(pdcov1), infot)

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'TEST 3', info
c      close(unit=1)

c
c           sub-block -> cov. matrix
            CALL YMCPIR (npk, dwork(pdcov1), n, iwork(piwbf),
     &                   dwork(pdcov), infot)

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'TEST 4', info
c      close(unit=1)

c
c           optimal portfolio volatility: sqrt[wopt'*Cov*wopt]
            CALL OVTMCV(n, dwork(pdcov), wopt, scal)
            scal = SQRT(scal)
            IF (scal .GT. (sigma(i)+1.E-6)) THEN
                info = -120
                RETURN
            ENDIF
         ENDIF
      ENDDO

c      open(unit=1,file='NUMXallocRiskBudgeting.txt',status='unknown')
c      write(1, *) 'TEST 5', info
c      close(unit=1)

c
c     BFGS diagnostic
      IF ((infoBFGS .NE. 0) .AND. (info .EQ. 0)) THEN
        info = infoBFGS
      ENDIF
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
Program Ejemplo15
!
!		Aqui veremos el Uso del lazo CYCLE
!
IMPLICIT NONE

INTEGER :: ii, jj, kk

ii = 10
jj = 10

 Loop2: Do ii = 1, 15

			Print *, ' ii : ', ii

			IF (ii .LT. 10) CYCLE Loop2

 Loop1:		DO kk = 1, 20

				Print *,'Loop 1'
					IF (kk.lt.12) CYCLE Loop1
					Print *,'kk : ', kk

			END DO Loop1

		END Do  Loop2

End Program Ejemplo15

max_stars_repo_name: EdgardoBonzi/Fortran-Examples

max_issues_repo_name: EdgardoBonzi/Fortran-Examples

max_forks_repo_name: EdgardoBonzi/Fortran-Examples


text:
!***********************************************************************
! Read and makes available climate data (tmin, tmax, precip, potential
! solar radiation, potential evapotranspieration) and/or transpiration
! on, by HRU from files pre-processed Data Files available for other
! PRMS modules
!***********************************************************************
      MODULE PRMS_CLIMATE_HRU
        USE PRMS_MODULE, ONLY: MAXFILE_LENGTH
        ! Local Variables
        INTEGER, SAVE :: Precip_unit, Tmax_unit, Tmin_unit, Et_unit, Swrad_unit, Transp_unit
        INTEGER, SAVE :: Humidity_unit, Windspeed_unit
        CHARACTER(LEN=11), SAVE :: MODNAME
        ! Control Parameters
        CHARACTER(LEN=MAXFILE_LENGTH), SAVE :: Tmin_day, Tmax_day, Precip_day, Potet_day, Swrad_day, Transp_day
        CHARACTER(LEN=MAXFILE_LENGTH), SAVE :: Humidity_day, Windspeed_day
        INTEGER, SAVE :: Cbh_check_flag, Cbh_binary_flag
        ! Declared Variables
        DOUBLE PRECISION, SAVE :: Basin_windspeed
        REAL, ALLOCATABLE :: Humidity_hru(:), Windspeed_hru(:)
        ! Declared Parameters
        REAL, SAVE, ALLOCATABLE :: Rain_cbh_adj(:, :), Snow_cbh_adj(:, :), Potet_cbh_adj(:, :)
        REAL, SAVE, ALLOCATABLE :: Tmax_cbh_adj(:, :), Tmin_cbh_adj(:, :)
      END MODULE PRMS_CLIMATE_HRU

      INTEGER FUNCTION climate_hru()
      USE PRMS_CLIMATE_HRU
      USE PRMS_MODULE, ONLY: Process, Nhru, Climate_transp_flag, Orad_flag, Model, &
     &    Climate_precip_flag, Climate_temp_flag, Climate_potet_flag, Climate_swrad_flag, &
     &    Start_year, Start_month, Start_day, Humidity_cbh_flag, Windspeed_cbh_flag
      USE PRMS_BASIN, ONLY: Active_hrus, Hru_route_order, Hru_area, Basin_area_inv, MM2INCH
      USE PRMS_CLIMATEVARS, ONLY: Solrad_tmax, Solrad_tmin, Basin_temp, &
     &    Basin_tmax, Basin_tmin, Tmaxf, Tminf, Tminc, Tmaxc, Tavgf, &
     &    Tavgc, Hru_ppt, Hru_rain, Hru_snow, Prmx, Pptmix, Newsnow, &
     &    Precip_units, Tmax_allrain_f, Adjmix_rain, &
     &    Basin_ppt, Basin_potet, Potet, Basin_snow, Basin_rain, &
     &    Basin_horad, Orad, Swrad, Basin_potsw, Basin_swrad, Basin_obs_ppt, &
     &    Transp_on, Basin_transp_on, Tmax_allsnow_f, Basin_humidity
      USE PRMS_SET_TIME, ONLY: Nowmonth, Jday
      USE PRMS_SOLTAB, ONLY: Soltab_basinpotsw, Hru_cossl, Soltab_potsw
      IMPLICIT NONE
! Functions
      INTRINSIC ABS, DBLE, SNGL
      INTEGER, EXTERNAL :: declparam, control_integer, getparam, control_string, declvar
      EXTERNAL :: read_error, precip_form, temp_set, find_header_end, find_current_time
      EXTERNAL :: read_cbh_date, check_cbh_value, check_cbh_intvalue, print_module
! Local Variables
      INTEGER :: yr, mo, dy, i, hr, mn, sec, jj, ierr, istop, missing, ios
      DOUBLE PRECISION :: sum_obs
      REAL :: tmax_hru, tmin_hru, ppt, harea
      CHARACTER(LEN=80), SAVE :: Version_climate_hru
!***********************************************************************
      climate_hru = 0
      ierr = 0
      IF ( Process(:3)=='run' ) THEN
        IF ( Climate_temp_flag==1 ) THEN
          IF ( Cbh_binary_flag==0 ) THEN
            READ ( Tmax_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Tmaxf(i), i=1,Nhru)
          ELSE
            READ ( Tmax_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Tmaxf(i), i=1,Nhru)
          ENDIF
          IF ( Cbh_check_flag==1 ) CALL read_cbh_date(yr, mo, dy, 'Tmaxf', ios, ierr)
          IF ( Cbh_binary_flag==0 ) THEN
            READ ( Tmin_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Tminf(i), i=1,Nhru)
          ELSE
            READ ( Tmin_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Tminf(i), i=1,Nhru)
          ENDIF
          IF ( Cbh_check_flag==1 ) CALL read_cbh_date(yr, mo, dy, 'Tminf', ios, ierr)
          Basin_tmax = 0.0D0
          Basin_tmin = 0.0D0
          Basin_temp = 0.0D0
        ENDIF

        IF ( Climate_precip_flag==1 ) THEN
          IF ( Cbh_binary_flag==0 ) THEN
            READ ( Precip_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Hru_ppt(i), i=1,Nhru)
          ELSE
            READ ( Precip_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Hru_ppt(i), i=1,Nhru)
          ENDIF
          IF ( Cbh_check_flag==1 ) CALL read_cbh_date(yr, mo, dy, 'Hru_ppt', ios, ierr)
          Basin_ppt = 0.0D0
          Basin_rain = 0.0D0
          Basin_snow = 0.0D0
          sum_obs = 0.0D0
        ENDIF

        IF ( Climate_potet_flag==1 ) THEN
          IF ( Cbh_binary_flag==0 ) THEN
            READ ( Et_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Potet(i), i=1,Nhru)
          ELSE
            READ ( Et_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Potet(i), i=1,Nhru)
          ENDIF
          IF ( Cbh_check_flag==1 ) CALL read_cbh_date(yr, mo, dy, 'Potet', ios, ierr)
          Basin_potet = 0.0D0
        ENDIF

        IF ( Climate_swrad_flag==1 ) THEN
          IF ( Orad_flag==0 ) THEN
            IF ( Cbh_binary_flag==0 ) THEN
              READ ( Swrad_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Swrad(i), i=1,Nhru)
            ELSE
              READ ( Swrad_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Swrad(i), i=1,Nhru)
            ENDIF
          ELSE
            IF ( Cbh_binary_flag==0 ) THEN
              READ ( Swrad_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Swrad(i), i=1,Nhru), Orad
            ELSE
              READ ( Swrad_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Swrad(i), i=1,Nhru), Orad
            ENDIF
          ENDIF
          IF ( Cbh_check_flag==1 ) CALL read_cbh_date(yr, mo, dy, 'Swrad', ios, ierr)
          Basin_swrad = 0.0D0
        ENDIF

        IF ( Climate_transp_flag==1 ) THEN
          IF ( Cbh_binary_flag==0 ) THEN
            READ ( Transp_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Transp_on(i), i=1,Nhru)
          ELSE
            READ ( Transp_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Transp_on(i), i=1,Nhru)
          ENDIF
          IF ( Cbh_check_flag==1 ) CALL read_cbh_date(yr, mo, dy, 'Transp_on', ios, ierr)
          Basin_transp_on = 0
        ENDIF

        IF ( Humidity_cbh_flag==1 ) THEN
          IF ( Cbh_binary_flag==0 ) THEN
            READ ( Humidity_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Humidity_hru(i), i=1,Nhru)
          ELSE
            READ ( Humidity_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Humidity_hru(i), i=1,Nhru)
          ENDIF
          IF ( Cbh_check_flag==1 ) CALL read_cbh_date(yr, mo, dy, 'Humidity_hru', ios, ierr)
          Basin_humidity = 0.0D0
        ENDIF

        IF ( Windspeed_cbh_flag==1 ) THEN
          IF ( Cbh_binary_flag==0 ) THEN
            READ ( Windspeed_unit, *, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Windspeed_hru(i), i=1,Nhru)
          ELSE
            READ ( Windspeed_unit, IOSTAT=ios ) yr, mo, dy, hr, mn, sec, (Windspeed_hru(i), i=1,Nhru)
          ENDIF
          IF ( Cbh_check_flag==1 ) CALL read_cbh_date(yr, mo, dy, 'Windspeed_hru', ios, ierr)
          Basin_windspeed = 0.0D0
        ENDIF

        IF ( ierr/=0 ) STOP

        missing = 0
        DO jj = 1, Active_hrus
          i = Hru_route_order(jj)
          harea = Hru_area(i)

          IF ( Climate_temp_flag==1 ) THEN
            IF ( Cbh_check_flag==1 ) THEN
              CALL check_cbh_value('Tmaxf', Tmaxf(i), -99.0, 150.0, missing)
              CALL check_cbh_value('Tminf', Tminf(i), -99.0, 150.0, missing)
            ENDIF
            tmax_hru = Tmaxf(i) + Tmax_cbh_adj(i, Nowmonth)
            tmin_hru = Tminf(i) + Tmin_cbh_adj(i, Nowmonth)
            CALL temp_set(i, tmax_hru, tmin_hru, Tmaxf(i), Tminf(i), &
     &                    Tavgf(i), Tmaxc(i), Tminc(i), Tavgc(i), harea)
          ENDIF

          IF ( Climate_potet_flag==1 ) THEN
            IF ( Cbh_check_flag==1 ) CALL check_cbh_value('Potet', Potet(i), 0.0, 50.0, missing)
            Potet(i) = Potet(i)*Potet_cbh_adj(i, Nowmonth)
            Basin_potet = Basin_potet + DBLE( Potet(i)*harea )
          ENDIF

          IF ( Climate_swrad_flag==1 ) THEN
            IF ( Cbh_check_flag==1 ) CALL check_cbh_value('Swrad', Swrad(i), 0.0, 1000.0, missing)
            Basin_swrad = Basin_swrad + DBLE( Swrad(i)*harea )
          ENDIF

          IF ( Climate_transp_flag==1 ) THEN
            IF ( Cbh_check_flag==1 ) CALL check_cbh_intvalue('Transp_on', Transp_on(i), 0, 1, missing)
            IF ( Transp_on(i)==1 ) Basin_transp_on = 1
          ENDIF

          IF ( Climate_precip_flag==1 ) THEN
            IF ( Cbh_check_flag==1 ) CALL check_cbh_value('Hru_ppt', Hru_ppt(i), 0.0, 30.0, missing)

!******Initialize HRU variables
            Pptmix(i) = 0
            Newsnow(i) = 0
            Prmx(i) = 0.0
            Hru_rain(i) = 0.0
            Hru_snow(i) = 0.0

            IF ( Hru_ppt(i)>0.0 ) THEN
              IF ( Precip_units==1 ) Hru_ppt(i) = Hru_ppt(i)*MM2INCH
              ppt = Hru_ppt(i)
              CALL precip_form(ppt, Hru_ppt(i), Hru_rain(i), Hru_snow(i), &
     &                         Tmaxf(i), Tminf(i), Pptmix(i), Newsnow(i), &
     &                         Prmx(i), Tmax_allrain_f(i,Nowmonth), &
     &                         Rain_cbh_adj(i,Nowmonth), Snow_cbh_adj(i,Nowmonth), &
     &                         Adjmix_rain(i,Nowmonth), harea, sum_obs, Tmax_allsnow_f(i,Nowmonth))
            ELSEIF ( Hru_ppt(i)<0.0 ) THEN
              PRINT *, 'ERROR, negative precipitation value entered in CBH File, HRU:', i
              CALL print_date(0)
              ierr = 1
!              Hru_ppt(i) = 0.0
            ENDIF
          ENDIF

          IF ( Humidity_cbh_flag==1 ) THEN
            IF ( Cbh_check_flag==1 ) CALL check_cbh_value('Humidity_hru', Humidity_hru(i), 0.0, 100.0, missing)
            IF ( missing==0 ) Basin_humidity = Basin_humidity + DBLE( Humidity_hru(i)*harea )
          ENDIF

          IF ( Windspeed_cbh_flag==1 ) THEN
            IF ( Cbh_check_flag==1 ) CALL check_cbh_value('Windspeed_hru', Windspeed_hru(i), 0.0, 400.0, missing)
            IF ( missing==0 ) Basin_windspeed = Basin_windspeed + DBLE( Windspeed_hru(i)*harea )
          ENDIF
        ENDDO

        IF ( missing==1 .OR. ierr==1 ) THEN
          CALL print_date(0)
          STOP
        ENDIF

        IF ( Climate_temp_flag==1 ) THEN
          Basin_tmax = Basin_tmax*Basin_area_inv
          Basin_tmin = Basin_tmin*Basin_area_inv
          Basin_temp = Basin_temp*Basin_area_inv
          Solrad_tmax = Basin_tmax
          Solrad_tmin = Basin_tmin
        ENDIF

        IF ( Climate_precip_flag==1 ) THEN
          Basin_ppt = Basin_ppt*Basin_area_inv
          Basin_obs_ppt = sum_obs*Basin_area_inv
          Basin_rain = Basin_rain*Basin_area_inv
          Basin_snow = Basin_snow*Basin_area_inv
        ENDIF
        IF ( Climate_potet_flag==1 ) Basin_potet = Basin_potet*Basin_area_inv
        IF ( Climate_swrad_flag==1 ) THEN
          Basin_horad = Soltab_basinpotsw(Jday)
          IF ( Orad_flag==0 ) Orad = SNGL( (DBLE(Swrad(1))*Hru_cossl(1)*Basin_horad)/Soltab_potsw(Jday,1) ) ! ??bad assumption using HRU 1
          Basin_swrad = Basin_swrad*Basin_area_inv
          Basin_potsw = Basin_swrad
        ENDIF
        IF ( Humidity_cbh_flag==1 ) Basin_humidity = Basin_humidity*Basin_area_inv
        IF ( Windspeed_cbh_flag==1 ) Basin_windspeed = Basin_windspeed*Basin_area_inv

      ELSEIF ( Process(:4)=='decl' ) THEN
        Version_climate_hru = 'climate_hru.f90 2017-10-23 15:57:00Z'
        MODNAME = 'climate_hru'

        IF ( control_integer(Cbh_check_flag, 'cbh_check_flag')/=0 ) Cbh_check_flag = 1
        IF ( control_integer(Cbh_binary_flag, 'cbh_binary_flag')/=0 ) Cbh_binary_flag = 0

        IF ( Climate_temp_flag==1 .OR. Model==99 ) CALL print_module(Version_climate_hru, 'Temperature Distribution    ', 90)
        IF ( Climate_precip_flag==1 .OR. Model==99 ) CALL print_module(Version_climate_hru, 'Precipitation Distribution  ', 90)
        IF ( Climate_swrad_flag==1 .OR. Model==99 ) CALL print_module(Version_climate_hru, 'Solar Radiation Distribution', 90)
        IF ( Climate_potet_flag==1 .OR. Model==99 ) CALL print_module(Version_climate_hru, 'Potential Evapotranspiration', 90)
        IF ( Climate_transp_flag==1 .OR. Model==99 ) CALL print_module(Version_climate_hru, 'Transpiration Distribution  ', 90)
        IF ( Humidity_cbh_flag==1 .OR. Model==99 ) THEN
          ALLOCATE ( Humidity_hru(Nhru) )
          IF ( declvar(MODNAME, 'humidity_hru', 'nhru', Nhru, 'real', &
     &         'Relative humidity of each HRU', &
     &         'percentage', Humidity_hru)/=0 ) CALL read_error(3, 'humidity_hru')
        ENDIF
        IF ( Windspeed_cbh_flag==1 .OR. Model==99 ) THEN
          IF ( declvar(MODNAME, 'basin_windspeed', 'one', 1, 'double', &
     &         'Basin area-weighted average wind speed', &
     &         'meters/second', Basin_windspeed)/=0 ) CALL read_error(3, 'basin_windspeed')
          ALLOCATE ( Windspeed_hru(Nhru) )
          IF ( declvar(MODNAME, 'windspeed_hru', 'nhru', Nhru, 'real', &
     &         'Area of HRU that is impervious', &
     &         'meters/second', Windspeed_hru)/=0 ) CALL read_error(3, 'windspeed_hru')
        ENDIF

!   Declared Parameters
        IF ( Climate_temp_flag==1 .OR. Model==99 ) THEN
          ALLOCATE ( Tmax_cbh_adj(Nhru,12) )
          IF ( declparam(MODNAME, 'tmax_cbh_adj', 'nhru,nmonths', 'real', &
     &         '0.0', '-10.0', '10.0', &
     &         'Monthly maximum temperature adjustment factor for each HRU', &
     &         'Monthly (January to December) adjustment factor to maximum air temperature for each HRU,'// &
     &         ' estimated on the basis of slope and aspect', &
     &         'temp_units')/=0 ) CALL read_error(1, 'tmax_cbh_adj')

          ALLOCATE ( Tmin_cbh_adj(Nhru,12) )
          IF ( declparam(MODNAME, 'tmin_cbh_adj', 'nhru,nmonths', 'real', &
     &         '0.0', '-10.0', '10.0', &
     &         'Monthly minimum temperature adjustment factor for each HRU', &
     &         'Monthly (January to December) adjustment factor to minimum air temperature for each HRU,'// &
     &         ' estimated on the basis of slope and aspect', &
     &         'temp_units')/=0 ) CALL read_error(1, 'tmin_cbh_adj')
        ENDIF

        IF ( Climate_precip_flag==1 .OR. Model==99 ) THEN
          ALLOCATE ( Rain_cbh_adj(Nhru,12) )
          IF ( declparam(MODNAME, 'rain_cbh_adj', 'nhru,nmonths', 'real', &
     &         '1.0', '0.5', '2.0', &
     &         'Rain adjustment factor, by month for each HRU', &
     &         'Monthly (January to December) adjustment factor to'// &
     &         ' measured precipitation determined to be rain on'// &
     &         ' each HRU to account for differences in elevation, and so forth', &
     &         'decimal fraction')/=0 ) CALL read_error(1, 'rain_cbh_adj')

          ALLOCATE ( Snow_cbh_adj(Nhru,12) )
          IF ( declparam(MODNAME, 'snow_cbh_adj', 'nhru,nmonths', 'real', &
     &         '1.0', '0.5', '2.0', &
     &         'Snow adjustment factor, by month for each HRU', &
     &         'Monthly (January to December) adjustment factor to'// &
     &         ' measured precipitation determined to be snow on'// &
     &         ' each HRU to account for differences in elevation, and so forth', &
     &         'decimal fraction')/=0 ) CALL read_error(1, 'snow_cbh_adj')
        ENDIF

        IF ( Climate_potet_flag==1 .OR. Model==99 ) THEN
          ALLOCATE ( Potet_cbh_adj(Nhru,12) )
          IF ( declparam(MODNAME, 'potet_cbh_adj', 'nhru,nmonths', 'real', &
     &         '1.0', '0.5', '1.5', &
     &         'Potential ET adjustment factor, by month for each HRU', &
     &         'Monthly (January to December) adjustment factor to'// &
     &         ' potential evapotranspiration specified in CBH Files for each HRU', &
     &         'decimal fraction')/=0 ) CALL read_error(1, 'potet_cbh_adj')
        ENDIF

      ELSEIF ( Process(:4)=='init' ) THEN
        Basin_humidity = 0.0D0
        Basin_windspeed = 0.0D0
        IF ( Humidity_cbh_flag==1 ) Humidity_hru = 0.0
        IF ( Windspeed_cbh_flag==1 ) Windspeed_hru = 0.0

        istop = 0
        ierr = 0

        IF ( Climate_precip_flag==1 ) THEN
          IF ( getparam(MODNAME, 'rain_cbh_adj', Nhru*12, 'real', Rain_cbh_adj)/=0 ) CALL read_error(2, 'rain_cbh_adj')
          IF ( getparam(MODNAME, 'snow_cbh_adj', Nhru*12, 'real', Snow_cbh_adj)/=0 ) CALL read_error(2, 'snow_cbh_adj')

          IF ( control_string(Precip_day, 'precip_day')/=0 ) CALL read_error(5, 'precip_day')
          CALL find_header_end(Precip_unit, Precip_day, 'precip_day', ierr, 1, Cbh_binary_flag)
          IF ( ierr==1 ) THEN
            istop = 1
          ELSE
            CALL find_current_time(Precip_unit, Start_year, Start_month, Start_day, ierr, Cbh_binary_flag)
            IF ( ierr==-1 ) THEN
              PRINT *, 'for first time step, CBH File: ', Precip_day
              istop = 1
            ENDIF
          ENDIF
        ENDIF

        IF ( Climate_temp_flag==1 ) THEN
          IF ( getparam(MODNAME, 'tmax_cbh_adj', Nhru*12, 'real', Tmax_cbh_adj)/=0 ) CALL read_error(2, 'tmax_cbh_adj')
          IF ( getparam(MODNAME, 'tmin_cbh_adj', Nhru*12, 'real', Tmin_cbh_adj)/=0 ) CALL read_error(2, 'tmin_cbh_adj')

          IF ( control_string(Tmax_day, 'tmax_day')/=0 ) CALL read_error(5, 'tmax_day')
          IF ( control_string(Tmin_day, 'tmin_day')/=0 ) CALL read_error(5, 'tmin_day')
          CALL find_header_end(Tmax_unit, Tmax_day, 'tmax_day', ierr, 1, Cbh_binary_flag)
          IF ( ierr==1 ) THEN
            istop = 1
          ELSE
            CALL find_current_time(Tmax_unit, Start_year, Start_month, Start_day, ierr, Cbh_binary_flag)
            IF ( ierr==-1 ) THEN
              PRINT *, 'for first time step, CBH File: ', Tmax_day
              istop = 1
            ENDIF
          ENDIF
          CALL find_header_end(Tmin_unit, Tmin_day, 'tmin_day', ierr, 1, Cbh_binary_flag)
          IF ( ierr==1 ) THEN
            istop = 1
          ELSE
            CALL find_current_time(Tmin_unit, Start_year, Start_month, Start_day, ierr, Cbh_binary_flag)
            IF ( ierr==-1 ) THEN
              PRINT *, 'for first time step, CBH File: ', Tmin_day
              istop = 1
            ENDIF
          ENDIF
        ENDIF

        IF ( Climate_potet_flag==1 ) THEN
          IF ( getparam(MODNAME, 'potet_cbh_adj', Nhru*12, 'real', Potet_cbh_adj)/=0 ) CALL read_error(2, 'potet_cbh_adj')
          IF ( control_string(Potet_day, 'potet_day')/=0 ) CALL read_error(5, 'potet_day')
          CALL find_header_end(Et_unit, Potet_day, 'potet_day', ierr, 1, Cbh_binary_flag)
          IF ( ierr==1 ) THEN
            istop = 1
          ELSE
            CALL find_current_time(Et_unit, Start_year, Start_month, Start_day, ierr, Cbh_binary_flag)
            IF ( ierr==-1 ) THEN
              PRINT *, 'for first time step, CBH File: ', Potet_day
              istop = 1
            ENDIF
          ENDIF
        ENDIF

        IF ( Climate_transp_flag==1 ) THEN
          IF ( control_string(Transp_day, 'transp_day')/=0 ) CALL read_error(5, 'transp_day')
          CALL find_header_end(Transp_unit, Transp_day, 'transp_day', ierr, 1, Cbh_binary_flag)
          IF ( ierr==1 ) THEN
            istop = 1
          ELSE
            CALL find_current_time(Transp_unit, Start_year, Start_month, Start_day, ierr, Cbh_binary_flag)
            IF ( ierr==-1 ) THEN
              PRINT *, 'for first time step, CBH File: ', Transp_day
              istop = 1
            ENDIF
          ENDIF
        ENDIF

        IF ( Climate_swrad_flag==1 ) THEN
          IF ( control_string(Swrad_day, 'swrad_day')/=0 ) CALL read_error(5, 'swrad_day')
          CALL find_header_end(Swrad_unit, Swrad_day, 'swrad_day', ierr, 1, Cbh_binary_flag)
          IF ( ierr==1 ) THEN
            istop = 1
          ELSE
            CALL find_current_time(Swrad_unit, Start_year, Start_month, Start_day, ierr, Cbh_binary_flag)
            IF ( ierr==-1 ) THEN
              PRINT *, 'for first time step, CBH File: ', Swrad_day
              istop = 1
            ENDIF
          ENDIF
        ENDIF

        IF ( Humidity_cbh_flag==1 ) THEN
          IF ( control_string(Humidity_day, 'humidity_day')/=0 ) CALL read_error(5, 'humidity_day')
          ierr = 2 ! signals routine to ignore CBH file requirement and use a parameter
          CALL find_header_end(Humidity_unit, Humidity_day, 'humidity_day', ierr, 1, Cbh_binary_flag)
          IF ( ierr==1 ) THEN
            istop = 1
          ELSEIF ( ierr==2 ) THEN
            Humidity_cbh_flag = 0
          ELSE
            CALL find_current_time(Humidity_unit, Start_year, Start_month, Start_day, ierr, Cbh_binary_flag)
            IF ( ierr==-1 ) THEN
              PRINT *, 'for first time step, CBH File: ', Humidity_day
              istop = 1
            ENDIF
          ENDIF
        ENDIF

        IF ( Windspeed_cbh_flag==1 ) THEN
          IF ( control_string(Windspeed_day, 'windspeed_day')/=0 ) CALL read_error(5, 'windspeed_day')
          CALL find_header_end(Windspeed_unit, Windspeed_day, 'windspeed_day', ierr, 1, Cbh_binary_flag)
          IF ( ierr==1 ) THEN
            istop = 1
          ELSE
            CALL find_current_time(Windspeed_unit, Start_year, Start_month, Start_day, ierr, Cbh_binary_flag)
            IF ( ierr==-1 ) THEN
              PRINT *, 'for first time step, CBH File: ', Windspeed_day
              istop = 1
            ENDIF
          ENDIF
        ENDIF

        IF ( istop==1 ) STOP 'ERROR in climate_hru'

      ENDIF

      END FUNCTION climate_hru

!***********************************************************************
!     Read a day in the CBH File
!***********************************************************************
      SUBROUTINE read_cbh_date(Year, Month, Day, Var, Ios, Iret)
      USE PRMS_SET_TIME, ONLY: Nowyear, Nowmonth, Nowday
! Argument
      INTEGER, INTENT(IN) :: Year, Month, Day, Ios
      CHARACTER(LEN=*), INTENT(IN) :: Var
      INTEGER, INTENT(INOUT) :: Iret
! Functions
      EXTERNAL :: print_date
! Local Variables
      INTEGER :: right_day
!***********************************************************************
      right_day = 1
      IF ( Year/=Nowyear .OR. Month/=Nowmonth .OR. Day/=Nowday ) right_day = 0
      IF ( Ios/=0 .OR. right_day==0 ) THEN
        PRINT *, 'ERROR, reading CBH File, variable: ', Var, ' IOSTAT=', Ios
        IF ( Ios==-1 ) THEN
          PRINT *, '       End-of-File found'
        ELSEIF ( right_day==0 ) THEN
          PRINT *, '       Wrong day found'
        ELSE
          PRINT *, '       Invalid data value found'
        ENDIF
        CALL print_date(0)
        Iret = 1
      ENDIF
      END SUBROUTINE read_cbh_date

!***********************************************************************
!     Check CBH value limits
!***********************************************************************
      SUBROUTINE check_cbh_value(Var, Var_value, Lower_val, Upper_val, Missing)
! Argument
      REAL, INTENT(IN) :: Lower_val, Upper_val
      REAL, INTENT(INOUT) :: Var_value
      CHARACTER(LEN=*), INTENT(IN) :: Var
      INTEGER, INTENT(INOUT) :: Missing
! Functions
      !INTRINSIC ISNAN
      EXTERNAL :: print_date
!***********************************************************************
      !IF ( ISNAN(Var_value) ) THEN
      !  PRINT *, 'ERROR, NaN value found for variable: ', Var
      !  Var_value = 0.0
      !  Missing = 1
      !  CALL print_date(0)
      !ELSEIF ( Var_value<Lower_val .OR. Var_value>Upper_val ) THEN
      IF ( Var_value<Lower_val .OR. Var_value>Upper_val ) THEN
        PRINT *, 'ERROR, bad value, variable: ', Var, ' Value:', Var_value
        PRINT *, '       lower bound:', Lower_val, ' upper bound:', Upper_val
        Missing = 1
        CALL print_date(0)
      ENDIF
      END SUBROUTINE check_cbh_value

!***********************************************************************
!     Check CBH integer value limits
!***********************************************************************
      SUBROUTINE check_cbh_intvalue(Var, Var_value, Lower_val, Upper_val, Missing)
! Argument
      INTEGER, INTENT(IN) :: Var_value, Lower_val, Upper_val
      CHARACTER(LEN=*), INTENT(IN) :: Var
      INTEGER, INTENT(INOUT) :: Missing
! Functions
      EXTERNAL :: print_date
!***********************************************************************
      IF ( Var_value<Lower_val .OR. Var_value>Upper_val ) THEN
        PRINT *, 'ERROR, bad value, variable: ', Var, ' Value:', Var_value
        PRINT *, '       lower bound:', Lower_val, ' upper bound:', Upper_val
        Missing = 1
        CALL print_date(0)
      ENDIF
      END SUBROUTINE check_cbh_intvalue

max_stars_repo_name: PJ-Thompson/GSFLOWurban

max_issues_repo_name: PJ-Thompson/GSFLOWurban

max_forks_repo_name: PJ-Thompson/GSFLOWurban


text:
program main
  use mod_monolis
  implicit none
  type(monolis_structure) :: mat !> 疎行列変数
  integer(kint) :: nnode, nelem, nbase_func, i
  real(kdouble) :: resid
  integer(kint), allocatable :: elem(:,:)
  real(kdouble), allocatable :: node(:,:), coef(:), B(:), X(:)
  character :: fname*128

  call monolis_set_debug(.true.)

  call monolis_global_initialize()

  fname = "node.dat"
  call monolis_input_mesh_node(fname, nnode, node)

  fname = "elem.dat"
  call monolis_input_mesh_elem(fname, nelem, nbase_func, elem)

  open(20, file = "coef.dat", status = "old")
    read(20,*) i
    if(i /= nelem) stop "ncoef /= nelem"
    allocate(coef(nelem), source = 0.0d0)
    do i = 1, nelem
      read(20,*) coef(i)
    enddo
  close(20)

  call monolis_initialize(mat, "./") !> 疎行列変数の初期化

  call monolis_get_nonzero_pattern(mat, nnode, nbase_func, 1, nelem, elem)

  do i = 1, nelem
    if(elem(1,i) == elem(2,i))then
      call monolis_add_scalar_to_sparse_matrix(mat, elem(1,i), elem(2,i), 1, 1, coef(i))
    else
      call monolis_add_scalar_to_sparse_matrix(mat, elem(1,i), elem(2,i), 1, 1, coef(i))
      call monolis_add_scalar_to_sparse_matrix(mat, elem(2,i), elem(1,i), 1, 1, coef(i))
    endif
  enddo

  allocate(B(nnode), source = 0.0d0)
  allocate(X(nnode), source = 0.0d0)

  X = 1.0d0
  call monolis_matvec_product(mat, X, B)
  X = 0.0d0

  call monolis_solve(mat, B, X)

  resid = 0.0d0
  do i = 1, nnode
    resid = resid + (X(i) - 1.0d0)**2
  enddo
  write(*,*)"resid: ", dsqrt(resid)

  call monolis_finalize(mat) !> 疎行列変数の解放

  call monolis_global_finalize()

end program main

max_stars_repo_name: nqomorita/monolis

max_issues_repo_name: nqomorita/monolis

max_forks_repo_name: nqomorita/monolis


text:
SUBROUTINE SIMULATION(ALPHA, BETA, GAMMA, DELTA, ARRAYS) BIND(C ,NAME='f_routine')
    USE ISO_C_BINDING
    IMPLICIT NONE
    INTEGER (C_LONG), VALUE :: alpha
    REAL (C_DOUBLE), INTENT(INOUT) :: beta
    REAL (C_DOUBLE), INTENT(OUT) :: gamma
    REAL (C_DOUBLE),DIMENSION(*),INTENT(IN) :: DELTA
    TYPE, BIND(C) :: PASS
        INTEGER (C_INT) :: LENC, LENF
        TYPE (C_PTR) :: C, F
    END TYPE PASS
    TYPE (PASS), INTENT(INOUT) :: ARRAYS
    REAL (C_FLOAT), ALLOCATABLE, TARGET, SAVE :: ETA(:)
    REAL (C_FLOAT), POINTER :: C_ARRAY(:)
    integer i,j
    write(*,'("In Fortran called from C alpha=",i4,&
              " beta=",f10.2," gamma=",f10.2)')&
            alpha,beta,gamma
    gamma=0.0
    do i=1,alpha
        gamma=gamma+beta*delta(i)
    enddo
    beta=1234.0

    !...
    write(*,*)"! Associate C_ARRAY with an array allocated in C"
    CALL C_F_POINTER (ARRAYS%C, C_ARRAY, (/ARRAYS%LENC/) )
    if(c_associated(ARRAYS%C, c_loc(C_ARRAY)))then
        write(*,*)'ARRAYS%C, C_ARRAY point to same target'
    else
        write(*,*)'ARRAYS%C, C_ARRAY do not point to same target'
        stop
    endif
    !...
    write(*,*)"! Allocate an array and make it available in C"
    ARRAYS%LENF = 100
    ALLOCATE (ETA(ARRAYS%LENF))
    ARRAYS%F = C_LOC(ETA)
    j=min(ARRAYS%lenc,ARRAYS%lenf)
    write(*,*)"Fortan fills the array for C"
    do i=1,j
        ETA(i)=C_ARRAY(i)*2
    enddo
END SUBROUTINE SIMULATION

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
C    @(#)reord2.f	20.4 11/11/97
      subroutine reord2(ktot, kdebug)

      include 'ipfinc/parametr.inc'

      include 'ipfinc/elim2.inc'
      include 'ipfinc/ikk.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/prt.inc'
      include 'ipfinc/red2.inc'

c     Note: This variable can be changed with a symbolic debugger
      common /term_debug/ iterm

      do i = 1, 800
        indexx(i) = 0
      enddo
      korder(0) = 0
      ich2(0) = 0
      ich1(0) = 0
      nopt = 0
      ncmax = 0
      ncmin = 0
      nzd = 0
      ix0 = 200
      do i = 1, ktot
        j =  -kownt(1, i)
        if (j .lt. 0) then
          kownt(1, i) = -j
          kownt(2, i) = 0
        else if (j .gt. 0) then
          nzd = nzd + 1
          iv = ivaln(i)
          iv = max0(iv, -ix0)
          iv = min0(iv, 800-ix0)
          kownt(1, i) =  -j
          kownt(2, i) = iv
          ncmin = min0(ncmin, iv)
          ncmax = max0(ncmax, iv)
          k = indexx(iv+ix0)
          ich1(i) = k
          if (k .gt. 0) ich2(k) = i
          indexx(iv+ix0) = i
          ich2(i) = 0
          ikk(4, i) = 1
        endif
      enddo

c     Define KBSKNT and KBRKNT to pertain to eliminated and envelope
c     network only.

      kbrknt = 0
      kbsknt = 0
      do i = 1, ktot
        if (kownt(1, i) .lt. 0) then

          kbsknt = kbsknt + 1
          l1 = loc(i)
          do while (l1 .ne. 0)
            kbrknt = kbrknt + 1
            l1 = korder(l1)
          enddo
        elseif (kownt(1, i) .ne. 0) then
          ic = 0
          l1 = loc(i)
          do while (l1 .ne. 0)
            k = kolum(l1)
            if (kownt(1, k) .lt. 0) ic = ic + 1
            l1 = korder(l1)
          enddo

          if (ic .gt. 0) ikk(3, i) = 1
        endif
      enddo
      do i = 1, ktot
        if (kownt(1, i) .gt. 0) then

          if (ikk(3, i) .ne. 0) then
            l1 = loc(i)
            kbsknt = kbsknt + 1
            do while (l1 .ne. 0)
              k = kolum(l1)
              if (kownt(1, k) .ge. 0) then
                if (ikk(3, k) .eq. 0) goto 100
              endif
              kbrknt = kbrknt + 1
  100         l1 = korder(l1)
            enddo
          endif
        endif
      enddo
      nmin = kbsknt + kbrknt
      if (kdebug .ne. 0) write (dbug, 10000)
10000 format ('1 NODE ELIMINATED  STEP  VALENCY  NODES REMAINING  BRAN'
     & , 'ches remaining  ratio  sparsity factor equivalent network ',
     & 'size '//)

c     NZD is the count of nodes in the eliminated system. All or part
c     may be eliminated. Only these nodes are candidates for elimination

      ksw = 0
      nelim = 0
      if (kdebug .eq. 3) call debug(ktot)
  110 continue
      do ncn = ncmin, ncmax
        i = indexx(ncn+ix0)
        if (i .ne. 0) goto 120
      enddo
      goto 140

  120 if (kdebug .eq. 3) call debug(ktot)
      ich2(0) = 0
      ich1(0) = 0

      ncmin = ncn

c     check validity of ICH1 <--> ICH2

      if (kdebug .ne. 0) then
        call rcheck(ncmin, ncmax, inum)
        j = nzd - nelim
        if (inum+nelim .ne. nzd) then
          write (dbug, 10010) nelim, inum, j, ncmin, ncmax
10010     format (' PARTIAL LOSS OF SUBSYSTEM AT ELIMINATION STEP ',
     &     5i6)
        endif
      endif
      inext = ich1(i)
      iv = ivaln(i)
      indexx(ncn+ix0) = inext
      ich2(inext) = 0
      ich1(i) = 0
      kownt(2, i) = 20000
      nelim = nelim + 1
      ikk(4, i) = 0
      ikk(5, i) = nelim
      if (nelim .ne. nzd) then
        if (ksw .ne. 1) then
          if (iv .ne. -ix0) then
            if (iv .lt. 800-ix0) goto 130
          endif
          iv = ivaln(i)
  130     continue
          kbsknt = kbsknt - 1
          kbrknt = kbrknt + iv
          ic = kbsknt + kbrknt
          if (kdebug .ne. 0) then
            write (dbug, 10020) i, nelim, iv, kbsknt, kbrknt, ic
10020       format (i13, i10, i7, i16, i18, 25x, i20)
          endif
          if (ic .le. nmin) then
            nmin = ic
            nopt = nelim
          endif

c         Build IX array of adjacent nodes

          itot = 0
          l1 = loc(i)
          lx = l1
          do while (l1 .ne. 0)
            itot = itot + 1
            ix(itot) = kolum(l1)
            kolum(l1) = 100000
            lx = l1
            l1 = korder(l1)
          enddo
          korder(last) = loc(i)
          loc(i) = 0
          last = lx
          if (kdebug .gt. 1) then
            write (dbug, 10030) i, itot, (ix(j), j = 1, itot)
10030       format (' ELIMINATING NODE ', i4, ' BRANCHES ', 11i6/(38x,
     &       10i6))
          endif
          if (last .ne. 0) then
            if (itot .eq. 0) goto 110

            call elimin(i, kerr)
c
c           Debug toplolgical symmetry check
c
            if (iterm .ne. 0) then
              write (dbug, 10032) i, nelim
              write (*, 10032) i, nelim
10032         format (' Eliminating node ', i4, ' at step ', i4)
              call ck_topol (ktot, kerr)
            endif

            if (kerr .eq. 0) then
              do j = 1, itot
                jbs = ix(j)
                j2 = kownt(1, jbs)
                iv = kownt(2, jbs)
                if (j2 .le. 0) then
                  ich2(0) = 0
                  ich1(0) = 0
                  j1 = ivaln(jbs)
                  j1 = max0(j1, -ix0)
                  j1 = min0(j1, 800-ix0)
                  if (j1 .ne. iv) then
                    if (ich2(jbs) .eq. 0) then
                      k = ich1(jbs)
                      indexx(iv+ix0) = k
                      ich2(k) = 0
                    else
                      l = ich2(jbs)
                      k = ich1(jbs)
                      if (l .gt. 0) ich1(l) = k
                      if (k .gt. 0) ich2(k) = l
                    endif
                    ich1(0) = 0
                    ich2(0) = 0
                    k = indexx(j1+ix0)
                    indexx(j1+ix0) = jbs
                    ich1(jbs) = k
                    if (k .gt. 0) ich2(k) = jbs
                    ich2(jbs) = 0
                    ncmin = min0(ncmin, j1)
                    ncmax = max0(ncmax, j1)
                  endif
                  kownt(1, jbs) = min0(0, j2-icon(j))
                  kownt(2, jbs) = j1
                endif
              enddo
              goto 110
            endif
          endif

c         Curtail futher elimination -- storage limit reached.

          write (errbuf(1), 10040) nelim
10040     format (
     &     ' AN EXTREMELY DENSE NETWORK WAS DEVELOPED AFTER ELIMINATINg
     &'
     &     , i4, ' nodes from system. ')
          call prterx('W', 1)
          ksw = 1
        endif
        goto 110
      endif
      if (kbsknt+kbrknt-1 .le. nmin) nopt = nelim
  140 if (nelim .gt. nopt) then
        write (outbuf, 10050) nopt, nelim
10050   format ('0 AN OPTIMALLY DETERMINED NETWORK HAS ELIMINATED ',
     &   i4, ' NOdes from a total of ', i4, ' candidates.')
        call prtout(1)
      endif
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


!*****************************************************************************************
!> author: Jacob Williams
!
!  Brent algorithms for minimization and root solving without derivatives.
!
!### See also
!
!  1. R. Brent, "Algorithms for Minimization Without Derivatives",
!     Prentice-Hall, Inc., 1973.

    module brent_module

    use kind_module,  only: wp
    use numbers_module

    implicit none

    private

    type,public :: brent_class
        !! the main class
        procedure(func),pointer :: f => null()  !! function to be minimized
    contains
        procedure :: set_function
        procedure :: minimize => fmin
        procedure :: find_zero => zeroin
    end type brent_class

    abstract interface
        function func(me,x) result(f)
            !! Interface to the function to be minimized.
            !! It should evaluate f(x) for any x in the interval (ax,bx)
            import :: wp,brent_class
            implicit none
            class(brent_class),intent(inout) :: me
            real(wp),intent(in) :: x
            real(wp) :: f
        end function func
    end interface

    !unit test routine:
    public :: brent_test

    contains
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 7/19/2014
!
!  Set the function to be minimized.

    subroutine set_function(me,f)

    implicit none

    class(brent_class),intent(inout) :: me
    procedure(func)                  :: f

    me%f => f

    end subroutine set_function
!*****************************************************************************************

!*****************************************************************************************
!>
!  An approximation x to the point where f attains a minimum on
!  the interval (ax,bx) is determined.
!
!  the method used is a combination of golden section search and
!  successive parabolic interpolation. convergence is never much slower
!  than that for a fibonacci search. if f has a continuous second
!  derivative which is positive at the minimum (which is not at ax or
!  bx), then convergence is superlinear, and usually of the order of
!  about 1.324.
!
!  the function f is never evaluated at two points closer together
!  than eps*abs(fmin) + (tol/3), where eps is approximately the square
!  root of the relative machine precision. if f is a unimodal
!  function and the computed values of f are always unimodal when
!  separated by at least eps*abs(x) + (tol/3), then fmin approximates
!  the abcissa of the global minimum of f on the interval ax,bx with
!  an error less than 3*eps*abs(fmin) + tol. if f is not unimodal,
!  then fmin may approximate a local, but perhaps non-global, minimum to
!  the same accuracy.
!
!  this function subprogram is a slightly modified version of the
!  algol 60 procedure localmin given in richard brent, algorithms for
!  minimization without derivatives, prentice - hall, inc. (1973).
!
!### See also
!  [1] http://www.netlib.org/fmm/fmin.f

    function fmin(me,ax,bx,tol) result(xmin)

    implicit none

    class(brent_class),intent(inout) :: me
    real(wp),intent(in) :: ax   !! left endpoint of initial interval
    real(wp),intent(in) :: bx   !! right endpoint of initial interval
    real(wp),intent(in) :: tol  !! desired length of the interval of
                                !! uncertainty of the final result (>=0)
    real(wp)            :: xmin !! abcissa approximating the point where
                                !! f attains a minimum

    real(wp) :: a,b,d,e,xm,p,q,r,tol1,tol2,u,v,w
    real(wp) :: fu,fv,fw,fx,x
    real(wp) :: abs,sqrt,sign

    real(wp),parameter :: c = (three-sqrt(five))/two  !! squared inverse of golden ratio
    real(wp),parameter :: half = 0.5_wp
    real(wp),parameter :: eps = sqrt(epsilon(one))

    ! initialization

    a = ax
    b = bx
    v = a + c*(b - a)
    w = v
    x = v
    e = zero
    fx = me%f(x)
    fv = fx
    fw = fx

    do    !  main loop starts here

        xm = half*(a + b)
        tol1 = eps*abs(x) + tol/three
        tol2 = two*tol1

        !  check stopping criterion

        if (abs(x - xm) <= (tol2 - half*(b - a))) exit

        ! is golden-section necessary

        if (abs(e) <= tol1) then

            !  a golden-section step

            if (x >= xm) then
                e = a - x
            else
                e = b - x
            end if
            d = c*e

        else

            !  fit parabola

            r = (x - w)*(fx - fv)
            q = (x - v)*(fx - fw)
            p = (x - v)*q - (x - w)*r
            q = two*(q - r)
            if (q > zero) p = -p
            q =  abs(q)
            r = e
            e = d

            !  is parabola acceptable

            if ((abs(p) >= abs(half*q*r)) .or. (p <= q*(a - x)) .or. (p >= q*(b - x))) then

                !  a golden-section step

                if (x >= xm) then
                    e = a - x
                else
                    e = b - x
                end if
                d = c*e

            else

                !  a parabolic interpolation step

                d = p/q
                u = x + d

                !  f must not be evaluated too close to ax or bx

                if (((u - a) < tol2) .or. ((b - u) < tol2)) d = sign(tol1, xm - x)

            end if

        end if

        !  f must not be evaluated too close to x

        if (abs(d) >= tol1) then
            u = x + d
        else
            u = x + sign(tol1, d)
        end if
        fu = me%f(u)

        !  update a, b, v, w, and x

        if (fu <= fx) then
            if (u >= x) a = x
            if (u < x) b = x
            v = w
            fv = fw
            w = x
            fw = fx
            x = u
            fx = fu
        else
            if (u < x) a = u
            if (u >= x) b = u
            if (fu <= fw .or. w == x) then
                v = w
                fv = fw
                w = u
                fw = fu
            else if (fu <= fv .or. v == x .or. v == w ) then
                v = u
                fv = fu
            end if
        end if

    end do    !  end of main loop

    xmin = x

    end function fmin
!*****************************************************************************************

!*****************************************************************************************
!>
!  Find a zero of the function \( f(x) \) in the given interval
!  \( [a_x,b_x] \) to within a tolerance \( 4 \epsilon |x| + tol \),
!  where \( \epsilon \) is the relative machine precision defined as
!  the smallest representable number such that \( 1.0 + \epsilon > 1.0 \).
!
!  It is assumed that \( f(a_x) \) and \( f(b_x) \) have opposite signs.
!
!### References
!  * R. P. Brent, "[An algorithm with guaranteed convergence for
!    finding a zero of a function](http://maths-people.anu.edu.au/~brent/pd/rpb005.pdf)",
!    The Computer Journal, Vol 14, No. 4., 1971.
!  * R. P. Brent, "[Algorithms for minimization without derivatives](http://maths-people.anu.edu.au/~brent/pub/pub011.html)",
!    Prentice-Hall, Inc., 1973.
!
!### See also
!  1. [zeroin.f](http://www.netlib.org/go/zeroin.f) from Netlib

    subroutine zeroin(me,ax,bx,tol,xzero,fzero,iflag,fax,fbx)

    use iso_fortran_env, only: error_unit

    implicit none

    class(brent_class),intent(inout) :: me
    real(wp),intent(in)              :: ax      !! left endpoint of initial interval
    real(wp),intent(in)              :: bx      !! right endpoint of initial interval
    real(wp),intent(in)              :: tol     !! desired length of the interval of uncertainty of the final result (>=0)
    real(wp),intent(out)             :: xzero   !! abscissa approximating a zero of `f` in the interval `ax`,`bx`
    real(wp),intent(out)             :: fzero   !! value of `f` at the root (`f(xzero)`)
    integer,intent(out)              :: iflag   !! status flag (`-1`=error, `0`=root found)
    real(wp),intent(in),optional     :: fax     !! if `f(ax)` is already known, it can be input here
    real(wp),intent(in),optional     :: fbx     !! if `f(ax)` is already known, it can be input here

    real(wp),parameter :: eps   = epsilon(one)  !! original code had d1mach(4)
    real(wp) :: a,b,c,d,e,fa,fb,fc,tol1,xm,p,q,r,s

    tol1 = eps+one

    a=ax
    b=bx

    if (present(fax)) then
        fa = fax
    else
        fa=me%f(a)
    end if
    if (present(fbx)) then
        fb = fbx
    else
        fb=me%f(b)
    end if

    !check trivial cases first:
    if (fa==zero) then

        iflag = 0
        xzero = a
        fzero = fa

    elseif (fb==zero) then

        iflag = 0
        xzero = b
        fzero = fb

    elseif (fa*(fb/abs(fb))<zero) then  ! check that f(ax) and f(bx) have different signs

        c=a
        fc=fa
        d=b-a
        e=d

        do

            if (abs(fc)<abs(fb)) then
                a=b
                b=c
                c=a
                fa=fb
                fb=fc
                fc=fa
            end if

            tol1=two*eps*abs(b)+0.5_wp*tol
            xm = 0.5_wp*(c-b)
            if ((abs(xm)<=tol1).or.(fb==zero)) exit

            ! see if a bisection is forced
            if ((abs(e)>=tol1).and.(abs(fa)>abs(fb))) then
                s=fb/fa
                if (a/=c) then
                    ! inverse quadratic interpolation
                    q=fa/fc
                    r=fb/fc
                    p=s*(two*xm*q*(q-r)-(b-a)*(r-one))
                    q=(q-one)*(r-one)*(s-one)
                else
                    ! linear interpolation
                    p=two*xm*s
                    q=one-s
                end if
                if (p<=zero) then
                    p=-p
                else
                    q=-q
                end if
                s=e
                e=d
                if (((two*p)>=(three*xm*q-abs(tol1*q))) .or. &
                    (p>=abs(0.5_wp*s*q))) then
                    d=xm
                    e=d
                else
                    d=p/q
                end if
            else
                d=xm
                e=d
            end if

            a=b
            fa=fb
            if (abs(d)<=tol1) then
                if (xm<=zero) then
                    b=b-tol1
                else
                    b=b+tol1
                end if
            else
                b=b+d
            end if
            fb=me%f(b)
            if ((fb*(fc/abs(fc)))>zero) then
                c=a
                fc=fa
                d=b-a
                e=d
            end if

        end do

        iflag = 0
        xzero = b
        fzero = fb

    else

        iflag = -1
        write(error_unit,'(A)')&
            'Error in zeroin: f(ax) and f(bx) do not have different signs.'

    end if

    end subroutine zeroin
!*****************************************************************************************

!*****************************************************************************************
!> author: Jacob Williams
!  date: 7/16/2014
!
!  Test of the [[fmin]] and [[zeroin]] functions.

    subroutine brent_test()

    implicit none

    real(wp) :: r,fzero
    integer :: iflag

    real(wp),parameter :: ax = zero
    real(wp),parameter :: bx = two*pi
    real(wp),parameter :: tol = 1.0e-6_wp

    type,extends(brent_class) :: myfunc_type
        integer :: i = 0    !! function counter
    end type myfunc_type
    type(myfunc_type) :: myfunc

    write(*,*) ''
    write(*,*) '---------------'
    write(*,*) ' brent_test'
    write(*,*) '---------------'
    write(*,*) ''

    call myfunc%set_function(sin_func)    !set the function

    !call fmin:
    ! [the minimum is at 270 deg]
    myfunc%i = 0
    r = myfunc%minimize(ax,bx,tol)
    write(*,*) 'minimum of sin(x) at: ', r*180.0_wp/pi,' deg'
    write(*,*) 'number of function calls: ', myfunc%i

    !call zeroin:
    ! [the root is at pi]
    myfunc%i = 0
    call myfunc%find_zero(ax+0.0001_wp,bx/two+0.0002,tol,r,fzero,iflag)
    write(*,*) 'root of sin(x) at: ', r*180.0_wp/pi,' deg'
    write(*,*) 'number of function calls: ', myfunc%i

    contains

        function sin_func(me,x) result(f)
        !! Example function to minimize: sin(x)

        implicit none

        class(brent_class),intent(inout) :: me
        real(wp),intent(in) :: x
        real(wp) :: f

        f = sin(x)

        select type (me)
        class is (myfunc_type)
            me%i = me%i + 1 !number of function calls
        end select

        end function sin_func

    end subroutine brent_test
 !*****************************************************************************************

!*****************************************************************************************
    end module brent_module
!*****************************************************************************************