!  VERSION 1.0  (NLSMC version)   2/5/99
!----------------------------------------------------------------------

!                         ====================
!                          subroutine MATRXO
!                         ====================

! *** MATRIX STORAGE SCHEME FOR CG ALGORITHM ***

!    This routine calculates the matrix elements off-diagonal subspace
! (pS=+/- 1, qS=0).  Nonsecular terms in the Hamiltonian is ignored
! due to the high field approximation.  This is done by not allowing
! the coupling between the spaces with different pS indices.  The only
! source of the coupling between diagonal and off-diagonal subspaces is
! pulse propagator.

!    M-symmetrization is included.  In the absence of nuclear Zeeman,
! the matrix becomes block-diagonal form and we need to keep only the
! block with jM=1 since the starting vector has non-zero elements in
! jM=1 and pulse propagator does not mix jM=1 and jM=-1 block.  In the
! presence of nuclear Zeeman, the block with jM=-1 is included in the
! calculation.

!    This routine is based upon EPRLL routines ver.1.3., whose features are
!      1. improved matrix storage algorithm (store only half of the matrix)
!      2. support nonaxial diffusion tensor R (in case of Brownian motion)
!      3. uses pruning scheme by eprbl and tnll.  eprbl and tnll generates
!         an <fileid>.ind file which keeps only the basis elements with
!         significant projection on the starting vector.

!    The phenomenological relaxation processes are also included as
!    T2ed, T2ef for the electronic relaxation.

!      T2ed : orientationally invariant electronic T2 relaxation
!             This defines the homogeneous line-width of ESR transition.
!      T2ef : angular variation of T2e values.  This angular variation
!             is naturally defined by the angle between B0 and the
!             molecular z-axis.  It is defined in the angle between
!             director and the molecular z-axis for convenience initially.
!             And then it is transformed to B0 is necessary.

!   *Note that the sign of these relaxation constants are reversed
!    with respect to the time domain formalism of the relaxation processes
!    since we are calculating the frequency spectrum directly.

!      Original program by G. Moro with modifications by D. Budil
!      Further modification by Sanghyuk Lee to calculate the slow-motional
!              2D-FT spectra.

!       flags
!       -----

!             fld   : True if |lr-lc| <=2.
!                     Used as flag for bandwidth of Liouville operator
!                     matrix elements. The limit on the L-bandwidth is due
!                     to the restriction of spin operator tensors to
!                     ranks 0 and 2.

!             fldmd : True if |lr-lc| <=2 or mr=mc, false otherwise.
!                     Used as flag for possible non-zero matrix
!                     elements. The L-bandwidth restriction is
!                     explained above and the diffusion operator
!                     matrix elements can be found only when mr=mc.

!             flk0  : True if lr=lc and kr=kc.
!                     Used as a flag for determing the existence of
!                     non-zero Liouville (A-tensor) and diffusion
!                     (Heisenberg spin exchange) matrix elements.

!	      fkm0  : True if kr=kc and mr=mc

!       Includes :
!               nlsdim.inc
!               eprprm.inc
!               indexf.inc
!               eprmat.inc
!               maxl.inc
!               rndoff.inc
!               physcn.inc

!       Uses:
!               cd2km.f
!               anxlk.f
!               ipar.f
!               w3j.f

!----------------------------------------------------------------------
!    subroutine matrxo(simparams_double, simparams_int, lstarts_offdiag, &
!                      jqe1,pi1,qi1,l1,jk1,k1,jm1,m1,ndimo,zmat,zdiag,izmat,jzmat,kzmat)
  subroutine matrxo(simparams_double, simparams_int, fad, fgd,lstarts_offdiag, &
                      jqe1,pi1,qi1,l1,jk1,k1,jm1,m1,ndimo,zmat,zdiag,izmat,&
                      jzmat,kzmat,nelim,nelre,ipt,cpot,lptmx,kptmx)
!    use eprmat
!    use basis
    implicit none

    include 'limits.inc' !now limits.inc has mxlval also
!    include 'simparm.inc'
!    include 'basis.inc'
!    include 'eprmat.inc'
    include 'rndoff.inc'
    include 'physcn.inc'
!   include 'mxlval.inc'

    logical :: fnz,fld,flk0,diftrm,sumtrm,fldmd
    logical :: newlr,newjkr,newkr,newjmr,newmr,newpnr, &
    newlc,newjkc,newkc,newjmc,newmc,newpnc
    integer :: ioldlr,ioldjkr,ioldkr,ioldjmr,ioldmr,ioldpnr, &
    ioldlc,ioldjkc,ioldkc,ioldjmc,ioldmc,ioldpnc

    double precision :: c,c1,c2,c3,c4,cdfd,cdfd1,cdfd2, &
    cdff,cdffu,cga0,cga2,cgam, &
    cliou,cliou0,cliou1,cliou2,clioua,clioua1, &
    clioua2,clioug,clioug1,clioug2, &
    cnorm,cnk,cnl,cnp,cpmkr,cplkc,cplmc,cnpr,cnpc,ct, &
    d1,d2,dsq2,dsq12,dsq18,dsq23, &
    fii,fki,ga,gg,ossc, &
    ra,ra1,ra2,rg,rg1,rg2,rpl,rmi,rz,rj,rp, &
    sa1,sg1,wliou1,wliou2,z, &
    t2ed,t2ef,ct2efi

    double precision :: zero,one
    parameter (zero=0.0d0,one=1.0d0)

    double precision :: d2km(2,5,5)
    integer :: i,j,l,li,ierr,iper, &
    ipnr,ipnc,ipnd,ipns,ipndab,ipnsab, &
    iqnr,iqnc,iqnd,iqndab, &
    lr,lc,ld,ldabs,lrprod,jkr,jkc,jkd,jmr,jmc,jmd, &
    kr,kc,kd,ks,kdabs,ksabs,kdi,ksi,kip, &
    mr,mc,md,ms,mdabs,msabs, &
    nrow,ncol,nel,nelr,neli, neltot

    integer :: ipar
    double precision :: w3j
    external ipar,w3j

    integer :: limcol

    double precision :: zmat(MXEL), zdiag(2,MXDIM)
    integer :: izmat(MXEL), jzmat(MXDIM+1), kzmat(MXDIM+1)
    double precision :: zmat_tmp(MXEL), zdiag_tmp(2,MXDIM)
    integer :: izmat_tmp(MXEL), jzmat_tmp(MXDIM+1), kzmat_tmp(MXDIM+1)

!    integer, dimension(:), allocatable :: jqe1,pi1,qi1,l1,jk1,k1,jm1,m1
    integer :: jqe1(MXDIM),pi1(MXDIM),qi1(MXDIM), &
               l1(MXDIM),jk1(MXDIM),k1(MXDIM),jm1(MXDIM), &
               m1(MXDIM)

    integer lstarts_offdiag(MXLVAL+1)

    double precision simparams_double(38)
    integer simparams_int(9)

    double precision :: gxx, gyy, gzz, axx, ayy, azz, dx, dy, dz, pml, pmxy, pmzz, djf, djfprp, oss, psi, &
                        ald, bed, gad, alm, bem, gam, c20, c22, c40, c42, &
                        c44, t2edi, t2ndi, t2efi, t1edi, t1ndi, b0, a0, g0, pl,pkxy, pkzz

    integer :: in2, ipdf, ist, igflg, iaflg, irflg, ml, mxy, mzz, lemx, lomx, kmx, mmx, ipnmx, jkmn, jmmn, &
               sptst, ipt, itd, itm, ipsi0, lband, kband, ldelta, kdelta, lptmx, kptmx, &
               nelre, nelim, ndimo, ndimo_in
    double precision :: faa(5), fgm(5), fam(2,5), fgd(2,5), fad(2,5), cpot(5,5), xlk(5,5), zeen

    integer :: inttmp1, inttmp2, cnt
    double precision :: tmp1, tmp2

    gxx=simparams_double(1); gyy=simparams_double(2); gzz=simparams_double(3);
    axx=simparams_double(4); ayy=simparams_double(5); azz=simparams_double(6);
    dx=simparams_double(7); dy=simparams_double(8); dz=simparams_double(9);
    pml=simparams_double(10); pmxy=simparams_double(11); pmzz=simparams_double(12);
    djf=simparams_double(13); djfprp=simparams_double(14);
    oss=simparams_double(15); psi=simparams_double(16);
    ald=simparams_double(17); bed=simparams_double(18); gad=simparams_double(19);
    alm=simparams_double(20); bem=simparams_double(21); gam=simparams_double(22);
    c20=simparams_double(23); c22=simparams_double(24);
    c40=simparams_double(25); c42=simparams_double(26); c44=simparams_double(27);
    t2edi=simparams_double(28);t2ndi=simparams_double(29); t2efi=simparams_double(30);
    t1edi=simparams_double(31); t1ndi=simparams_double(32);
    b0=simparams_double(33); a0=simparams_double(34); g0=simparams_double(35);
    pl=simparams_double(36); pkxy=simparams_double(37);
    pkzz=simparams_double(38);

    in2=simparams_int(1); ipdf=simparams_int(2); ist=simparams_int(3); lemx=simparams_int(4);
    lomx=simparams_int(5); kmx=simparams_int(6); mmx=simparams_int(7); ipnmx=simparams_int(8);
    ndimo_in=simparams_int(9)
    !write(*,*) 'w3j(8,2,7,2,0,-2)=',w3j(8,2,7,2,0,-2)
    kband=kptmx*2
    if (dabs(dx-dy) > rndoff) kband=kband+2
    if (lptmx >= 2) then
        lband=lptmx*2
    else
        lband=2
    end if

!######################################################################
    zeen = zero !deliberate, ignoring nuclear Zeeman
    fnz = .FALSE.
    if(abs(zeen) > rndoff) fnz= .TRUE.

!----------------------------------------------------------------------
!     define constants
!----------------------------------------------------------------------

    dsq2=dsqrt(2.0D0)
    dsq12=one/dsq2
    dsq18=one/dsqrt(8.0D0)
    dsq23=dsqrt(2.0D0/3.0D0)

!----------------------------------------------------------------------
!     Scale diffusion constants
!----------------------------------------------------------------------

    ct=g0*betae/hbar
    rpl=0.5d0*(dx+dy)/ct
    rmi=0.25d0*(dx-dy)/ct
    rz=dz/ct
    rj=djf/ct
    rp=djfprp/ct
    ossc=oss/ct

    t2ed=t2edi/ct
    t2ef=t2efi/ct

    fii=0.25d0*dble(in2*(in2+2))

!    Non-Brownian motional parameters

    if(ipdf /= 1) then
        pml=zero
        pmzz=zero
        pmxy=zero
    end if

!   Correction term for anisotropic viscosity

    if((ipdf == 2) .AND. (ist == 0)) then
        rpl=rpl+rp
        rz=rz+rp
    end if

!----------------------------------------------------------------------
!                                 2
!     Get Wigner rotation matrix D  (0,psi,0) for director tilt.
!                                 KM

!     Note: resultant matrix will be real-valued and d2km(2,i,j)=zero

!                              2
!           d2km(1,i+3,j+3) = d  (psi)
!                              ij
!----------------------------------------------------------------------

    call cd2km(d2km,zero,psi,zero)

!----------------------------------------------------------------------
!                               L
!     Calculate the quantities X  used in the potential-dependent
!                               K
!     portion of the diffusion operator
!----------------------------------------------------------------------

    !call anxlk(rpl,rmi,rz)
    call anxlk(ipt,lptmx,kptmx,cpot,rpl,rmi,rz,xlk)
!    do 567 i=1,5
!      write(*,*) (xlk(i,j),j=1,5)
!567 end do

!----------------------------------------------------------------------
!     Initialize counters
!----------------------------------------------------------------------

    nelr=0
    neli=0
    nel=0

    iper=1
    ioldlr=-1
    ioldjkr=0
    ioldkr=1000
    ioldjmr=0
    ioldmr=1000
    ioldpnr=-in2-1

!----------------------------------------------------------------------
!     **** Loop over rows ***
!----------------------------------------------------------------------
    !write(*,*) 'rpl,rmi,dx,dy,rz',rpl,rmi,dx,dy,rz
    do 200 nrow=1,ndimo
        lr=l1(nrow)
        jkr=jk1(nrow)
        kr=k1(nrow)
        jmr=jm1(nrow)
        mr=m1(nrow)
        ipnr=pi1(nrow)
        iqnr=qi1(nrow)

        newlr=lr.ne.ioldlr
        newjkr=newlr.or.(jkr.ne.ioldjkr)
        newkr=newjkr.or.(kr.ne.ioldkr)
        newjmr=newkr.or.(jmr.ne.ioldjmr)
        newmr=newjmr.or.(mr.ne.ioldmr)
        newpnr=newmr.or.(ipnr.ne.ioldpnr)

        if(newlr) then
            ioldlr=lr
            lrprod=lr*(lr+1)
        end if

        if(newjkr) ioldjkr=jkr

        if(newkr) then
            ioldkr=kr

        !         ** Calculate isotropic part of diffusion operator

            c1=dble(lrprod)
            if(ipdf /= 0) then
                c2=one+pml*c1
                if(dabs(pl) > rndoff) c2=c2**pl
                c1=c1/c2
            end if
            cdfd=rpl*c1

            c1=dble(kr*kr)
            c2=rz
            c3=rpl
            if(ipdf /= 0) then
                c4=one+pmzz*c1
                if(dabs(pkzz) > rndoff) c4=c4**pkzz
                c2=c2/c4
                c4=one+pmxy*c1
                if(dabs(pkxy) > rndoff) c4=c4**pkxy
                c3=c3/c4
            end if
            cdfd=cdfd+c1*(c2-c3)

        !         ** Discrete jump motion

            if(ist /= 0) then
                i=kr/ist
                if((i*ist) /= kr) cdfd=cdfd+rj
            end if

        !         Rotationally invariant line width : T2 electronic relaxation

            if(iper /= 0) cdfd=cdfd+t2ed

            cdfd1=cdfd

        end if

        if(newjmr) ioldjmr=jmr

        if(newmr) then
            ioldmr=mr
            cpmkr=ipar(mr+kr)

        !   Anisotropic viscosity term in diffusion operator

            if(ipdf == 2) cdfd=cdfd1+dble(mr*mr)*(rj-rp)

        end if

        if(newpnr) then
            ioldpnr=ipnr
            if((ipnr == 0) .AND. (mr == 0)) then
                cnpr=dsq2
            else
                cnpr=1.0D0
            end if
        end if

    !----------------------------------------------------------------------
    !   **** Loop over columns (note only upper diagonal) ****
    !----------------------------------------------------------------------

        ioldlc=-1
        ioldjkc=0
        ioldkc=1000
        ioldjmc=0
        ioldmc=1000
        ioldpnc=-in2-1

        jzmat(nrow)=neli+1
        kzmat(nrow)=nelr+1
!        write(*,*) jzmat(nrow), kzmat(nrow)

!        do 300 ncol=nrow,ndimo
        if(lr.le.(lemx-lband-1)) then
            limcol=lstarts_offdiag(lr+lband+2)
        else
            limcol=ndimo
        end if
        do 300 ncol=nrow,limcol
            lc=l1(ncol)
            jkc=jk1(ncol)
            kc=k1(ncol)
            jmc=jm1(ncol)
            mc=m1(ncol)
            ipnc=pi1(ncol)
            iqnc=qi1(ncol)

            newlc=lc.ne.ioldlc
            newjkc=newlc.or.(jkc.ne.ioldjkc)
            newkc=newjkc.or.(kc.ne.ioldkc)
            newjmc=newkc.or.(jmc.ne.ioldjmc)
            newmc=newjmc.or.(mc.ne.ioldmc)
            newpnc=newmc.or.(ipnc.ne.ioldpnc)

            if(newlc) then
                ioldlc=lc
                ld=lr-lc
                ldabs=iabs(ld)
                fld=ldabs.le.2
                cnl=dsqrt((2.D0*lr+1.D0)*(2.D0*lc+1.D0))

            !           Angular variation of the line width : T2 electronic relaxation

                ct2efi=w3j(lr,2,lc,mr,0,-mr)*w3j(lr,2,lc,kr,0,-kr)
                ct2efi=cnl*cpmkr*t2ef*ct2efi*d2km(1,3,3)

            end if

            if(newjkc) then
                ioldjkc=jkc
                jkd=jkr-jkc
            end if

            if(newkc) then
                ioldkc=kc
                kd=kr-kc
                ks=kr+kc
                kdabs=iabs(kd)
                ksabs=iabs(ks)
                cplkc=ipar(lc+kc)

                flk0=(ld.eq.0).and.(kd.eq.0)

                if((kr == 0) .AND. (kc == 0)) then
                    cnk=0.5D0
                else if((kr /= 0) .AND. (kc /= 0)) then
                    cnk=one
                else
                    cnk=dsq12
                end if

            !---------------------------------------------------------------------
            !       Calculate R(mu,l) as in Meirovitch, Igner, Igner, Moro, & Freed
            !       J. Phys. Chem. 77 (1982) p. 3935, Eqs. A42 & A44
            !---------------------------------------------------------------------

                ra=zero
                rg=zero
                if(fld) then
                    ra1=zero
                    rg1=zero
                    if(kdabs <= 2) then
                        if(jkd == 0) then
                            ga=fad(1,kd+3)
                            gg=fgd(1,kd+3)
                        else
                            ga=fad(2,kd+3)*jkr
                            gg=fgd(2,kd+3)*jkr
                        end if
                        z=w3j(lr,2,lc,kr,-kd,-kc)
                        ra1=ga*z
                        rg1=gg*z
                    end if

                    ra2=zero
                    rg2=zero
                    if(ksabs <= 2) then
                        if(jkd == 0) then
                            ga=fad(1,ks+3)
                            gg=fgd(1,ks+3)
                        else
                            ga=fad(2,ks+3)*jkr
                            gg=fgd(2,ks+3)*jkr
                        end if
                        z=w3j(lr,2,lc,kr,-ks,kc)
                        ra2=ga*z
                        rg2=gg*z
                    end if
                !                                      --- for g tensor
                    rg=rg1+cplkc*jkc*rg2

                !                                      --- for A tensor
                    if(in2 /= 0) ra=ra1+cplkc*jkc*ra2
                end if

            !----------------------------------------------------------------------
            !       Calculate off-diagonal terms of the diffusion operator,
            !       including the potential-dependent portion and the rhombic
            !       component of the isotropic diffusion operator.

            !       See Meirovitch, Igner, Igner, Moro,and Freed,
            !       J. Chem. Phys. 77 (1982) pp.3933-3934, and especially
            !       Eqns. A22-24 and A40 for more information.
            !----------------------------------------------------------------------

            !           -------- Rhombic part of isotropic diffusion operator

                if((ld == 0) .AND. (kr+2 == kc)) then
                    cdff=rmi*dsqrt(dble((lr+kc-1)*(lr+kc)* &
                    (lr-kc+1)*(lr-kc+2)))/cnk
                else
                    cdff=zero
                end if

            !           -------- Potential-dependent part of diffusion operator

                if((ipt /= 0) .AND. (ldabs <= lband) .AND. (ipar(ks) == 1) &
                 .AND. (jkd == 0) .AND. ((kdabs <= kband) .OR. &
                (ksabs <= kband))) then

                    kdi=1+kdabs/2
                    ksi=1+ksabs/2
                    c1=cpmkr*cnl*cnk

                !             Loop over L and K indices in sum over terms in potential

                    do 444 l=0,lband,2
                        cdffu=0.0D0
                        li=1+l/2
                        if (ksi <= li .AND. abs(xlk(li,ksi)) >= rndoff) &
                        cdffu=cplkc*jkc*xlk(li,ksi)*w3j(lr,l,lc,kr,-ks,kc)
                        if (kdi <= li .AND. abs(xlk(li,kdi)) >= rndoff) &
                        cdffu=cdffu+xlk(li,kdi)*w3j(lr,l,lc,kr,-kd,-kc)

                        if (abs(cdffu) > rndoff) cdff=cdff+ &
                        w3j(lr,l,lc,mr,0,-mr)*c1*cdffu
!                        if(cdff.gt.rndoff) write(*,*) 'cdffu=',cdffu
                    444 END DO

                end if
            end if

            if(newjmc) then
                ioldjmc=jmc
                jmd=jmr-jmc
            end if

            if(newmc) then
                ioldmc=mc
                md=mr-mc
                ms=mr+mc
                mdabs=iabs(md)
                msabs=iabs(ms)
                cplmc=ipar(lc+mc)

            !----------------------------------------------------------------------
            !     set flags and zero out matrix elements
            !     if no non-zero elements are possible
            !----------------------------------------------------------------------

                if(fld .OR. (md == 0)) then
                    fldmd=.true.
                else
                    cliou=0.0D0
                    cgam=0.0D0
                    fldmd=.false.
                end if

            !----------------------------------------------------------------------
            !     get appropriate 3-J symbols if non-zero Liouville
            !     operator matrix elements are possible
            !----------------------------------------------------------------------

                wliou1=zero
                wliou2=zero
                if(fld) then
                    wliou1=w3j(lr,2,lc,mr,-md,-mc)
                    wliou2=w3j(lr,2,lc,mr,-ms,mc)
                end if

            end if

            if(newpnc) then
                ioldpnc=ipnc
                ipnd=ipnr-ipnc
                ipndab=iabs(ipnd)
                ipns=ipnr+ipnc
                ipnsab=iabs(ipns)

                if((ipnc == 0) .AND. (mc == 0)) then
                    cnpc=dsq2
                else
                    cnpc=1.0D0
                end if

                cnp=1.0D0/(cnpr*cnpc)
                cnorm=cnl*cnk*cnp*cpmkr

            !----------------------------------------------------------------------
            !     get director tilt dependence of Liouville
            !     operator matrix elements
            !----------------------------------------------------------------------

                d1=zero
                d2=zero
                if((ipndab <= 2) .AND. (mdabs <= 2)) &
                d1=d2km(1,ipnd+3,md+3)*wliou1

                if((ipnsab <= 2) .AND. (msabs <= 2)) &
                d2=d2km(1,ipns+3,ms+3)*wliou2

            end if

            if(fldmd) then
                iqnd=iqnr-iqnc
                iqndab=iabs(iqnd)

            !----------------------------------------------------------------------
            !    Liouville operator

            !    Note: see Meirovitch, Igner,Igner,Moro, and Freed
            !    J. Chem. Phys. 77 (1982) Appendix B p.3936-3937 for details.

            !          Sa and Sg in this program also contain the
            !          appropriate Clebsch-Gordon coefficient, and fki is
            !          Ki in text.

            !          The manner in which the isotropic contributions from the
            !          spin Hamiltonian are included in this program is not clear
            !          from the text, and  deserves some comments.
            !          When Lc=Lr, both the (script) l=0 and l=2 ISTO's must
            !          be included, whereas only l=2 terms appear off the diagonal.
            !          The Clebsch-Gordan coefficients for these terms are
            !          c(1,0,1,0;0,0) and c(1,+/-1,1,-/+1;0,0), or -/+ 1/sqrt(3),
            !          respectively. Here, cga0 is set to +/- 1 instead of -/+
            !          1/sqrt(3) because the extra factors are already included in g0
            !          and a0 (note sign change). The isotropic terms are not
            !          normalized by any factor.
            !	   The cnk factor is not needed because the l=0 ISTO components
            !	   have not been transformed by the K-symmetrization.
            !          The factor cnl*cpmkr is canceled by the product of
            !          w3j(lr,mr,0,0,lr,-mr) and w3j(lr,kr,0,0,lr,-kr), which
            !          therefore do not need to be calculated.

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

                cliou=zero

                if(fld) then
                    clioua=zero
                    clioug=zero
                    cliou0=zero
                    cliou1=zero
                    cliou2=zero

                !        **   Hyperfine interaction   **

                    if((jmd == 0) .AND. (in2 /= 0)) then
                    !                                         *** difference term ***
                        clioua1=zero
                        cliou1=zero
                        if((ipndab == iqndab) .AND. (ipndab <= 1) .AND. (mdabs.le.2))  then

                            cga0=zero
                            cga2=zero
                            if(ipnd /= 0) then
                            !                                         *** pseudosecular term ***
                            !					      (delta pI .ne. 0)
                                kip=iqnr*iqnd+ipnr*ipnd
                                kip=kip*(kip-2)
                                fki=dsqrt(fii-0.25D0*kip)
                                sa1=-iper*ipnd*fki*dsq18
                                cga2=dsq12
                            else
                            !                                         *** secular term ***
                            !					      (delta pS,pI .eq. 0)
                                sa1=iper*iqnr*0.5D0
                                cga0=one
                                cga2=dsq23
                            end if
                            clioua1=sa1*d1*ra*cga2
                            if(flk0 .AND. (jkd == 0) .AND. (md == 0)) &
                            cliou1=sa1*cga0*a0
                        end if

                    !                                         *** sum term ***
                        clioua2=zero
                        cliou2=zero
                        if((ipnsab == iqndab) .AND. (ipnsab <= 1) .AND. (msabs.le.2))  then

                            cga0=zero
                            cga2=zero
                            if(ipns /= 0) then
                            !                                         *** pseudosecular term ***
                            !					      (delta pI .ne. 0)
                                kip=iqnr*iqnd+ipnr*ipns
                                kip=kip*(kip-2)
                                fki=dsqrt(fii-0.25D0*kip)
                                sa1=-iper*ipns*fki*dsq18
                                cga2=dsq12
                            else
                            !                                         *** secular term ***
                            !					      (delta pS,pI .eq. 0)
                                sa1=iper*iqnr*0.5D0
                                cga0=one
                                cga2=dsq23
                            end if
                            clioua2=sa1*d2*ra*cga2
                            if(flk0 .AND. (jkd == 0) .AND. (ms == 0)) &
                            cliou2=sa1*cga0*a0
                        end if

                        clioua=clioua1+jmc*cplmc*clioua2
                        cliou0=cliou1+cliou2

                    end if

                !        **   Electronic Zeeman interaction   **

                    cliou1=zero
                    cliou2=zero
                    clioug1=zero
                    clioug2=zero
                    if((jmd == 0) .AND. (iqnd == 0)) then

                        if(flk0 .AND. (jkd == 0) .AND. (md == 0) .AND. &
                        (ipnd == 0))            cliou1=iper*b0

                        if((ipnd == 0) .AND. (abs(rg) > rndoff)) then
                            sg1=iper*dsq23
                            clioug1=sg1*d1*rg
                        end if

                        if(flk0 .AND. (jkd == 0) .AND. (ms == 0) .AND. &
                        (ipns == 0))            cliou2=iper*b0

                        if((ipns == 0) .AND. (abs(rg) > rndoff)) then
                            sg1=iper*dsq23
                            clioug2=sg1*d2*rg
                        end if

                        clioug=clioug1+jmc*cplmc*clioug2
                        cliou0=(cliou0+cliou1+cliou2)*cnp

                    end if

                !        **   Nuclear Zeeman interaction   **

                    cliou1=zero
                    cliou2=zero
                    if(fnz .AND. (jmd /= 0) .AND. flk0 .AND. (jkd == 0).and.(iqnd == 0)) then
                        if((ipnd == 0) .AND. (md == 0)) cliou1=ipnr*zeen
                        if((ipns == 0) .AND. (ms == 0)) cliou2=ipnr*zeen
                        cliou0=cliou0+cnp*(cliou1+cliou2)
                    end if

                    cliou=cnorm*(clioua+clioug)+cliou0
                    !write(*,*) 'cliou=', cliou
                end if

            !----------------------------------------------------------------------
            !           Diffusion operator
            !----------------------------------------------------------------------

                cgam=zero

            !      **   Heisenberg exchange   **

                if((dabs(ossc) > rndoff) .AND. (ipnd == 0) .AND. flk0 &
                 .AND. (md == 0) .AND. (jkd == 0) .AND. (jmd == 0)) then
                    c=zero
                    if(iqnd == 0) c=one
                    if((ipnr == 0) .AND. (lr == 0)) c=c-one/dble(in2+1)
                    cgam=c*ossc
                end if

            !      **   Potential-dependent terms   **

                if((ipnd == 0) .AND. (iqnd == 0) .AND. (md == 0) .AND. &
                (jkd == 0) .AND. (jmd == 0)) then
                    cgam=cgam+cdff
                    if (nrow == 1 .and. ncol == 316) write(*,*) 'pot-dep term (1,316) cgam=', cgam
                    if(fld .AND. (kd == 0))   cgam=cgam+ct2efi
                end if

            end if

        !----------------------------------------------------------------------
        !         Store matrix elements
        !----------------------------------------------------------------------
!            write(*,*) cgam
!            if(nrow==91 .and. ncol==541) write(*,*) '91,541, cgam, cliou=',cgam, cliou
            if(nrow == ncol) then
!                write(*,*) 'cdfd=',cdfd
                cgam=cgam+cdfd
                zdiag(1,nrow)=cgam
                zdiag(2,nrow)=-cliou
            else

                if(abs(cliou) > rndoff) then
                    nel=nel+1
                    if(nel > MXEL) then
                        ierr=1
                        write(*,*) 'Critical error; nel > MXEL!!!'
                        return
                    else
                        neli=neli+1
                        !write(*,*) '-cliou', -cliou
                        zmat(neli)=-cliou
                        !write(*,*) 'neli 1+, ncol=', ncol
                        izmat(neli)=ncol
                    end if
                end if

                if(abs(cgam) > rndoff) then
                    nel=nel+1
                    if(nel > MXEL) then
                        ierr=1
                        write(*,*) 'Critical error; nel > MXEL!!!'
                        return
                    else
                        nelr=nelr+1
                        !write(*,*) cgam
                        zmat(MXEL-nelr+1)=cgam
                        !write(*,*) 'nelr 1+, ncol=', ncol
                        izmat(MXEL-nelr+1)=ncol
                    end if
                end if
            end if

        !----------------------------------------------------------------------
        !     Increment column counter
        !----------------------------------------------------------------------

        !  ****test*****
        !         if (nrow.eq.1) then
        !         if (ncol.eq.1)       write(6,332)
        !         write(6,333) ncol,iqec,lc,jkc,kc,mc,ipnc,iqnc
        ! 332     format(2x,' ncol  iqec   lc    jkc   mc   ipnc  iqnc')
        ! 333     format(2x,8i6)
        !         end if

        !----------------------------------------------------------------------
        !         end loop over columns
        !----------------------------------------------------------------------
        300 END DO
    !----------------------------------------------------------------------
    !       end loop over rows
    !----------------------------------------------------------------------
    200 END DO

    ierr=0
    neltot=nel
    nelre=nelr
    nelim=neli

    jzmat(ndimo+1)=neli+1
    kzmat(ndimo+1)=nelr+1
    go to 1500
     write(*,*) nelre, nelim, ndimo
     open(unit=9,file='matrx.mtxx',status='unknown',access='sequential')
     read(9,*) (zdiag_tmp(1,i),zdiag_tmp(2,i),i=1,ndimo)
     read(9,*) (jzmat_tmp(i),i=1,ndimo+1)
     read(9,*) (izmat_tmp(i),zmat_tmp(i),i=1,nelim)
     read(9,*) (kzmat_tmp(i),i=1,ndimo+1)
     read(9,*) (izmat_tmp(MXEL-i+1),zmat_tmp(MXEL-i+1),i=1,nelre)
     do 33 i=1,ndimo
       zdiag(1,i) = zdiag(1,i) - zdiag_tmp(1,i)
       zdiag(2,i) = zdiag(2,i) - zdiag_tmp(2,i)
33   end do
     do 34 i=1,ndimo+1
       jzmat(i)=jzmat(i)-jzmat_tmp(i) !(jzmat(i),i=1,ndimo+1)
34   end do
     do 36 i=1,nelim
       izmat(i)=izmat(i)-izmat_tmp(i)
       zmat(i)=zmat(i)-zmat_tmp(i)
36   end do
     do 37 i=1,ndimo+1
       kzmat(i)=kzmat(i)-kzmat_tmp(i)
37   end do
     do 38 i=1,nelre   !where is neltot????
       izmat(MXEL-i+1)=izmat(MXEL-i+1)-izmat_tmp(MXEL-i+1)
       zmat(MXEL-i+1)=zmat(MXEL-i+1)-zmat_tmp(MXEL-i+1)
38   end do
     close(unit=9)
1500 continue

!   open (unit=9,file='matlab_real_new.mtxx',status='unknown', access='sequential')
!   cnt=0
!   do 44 i=1,ndimo
!     write (9,*) i,i,zdiag(1,i) ! get the diagonal element out first
!     do 45 j=kzmat(i),(kzmat(i+1)-1)
!       cnt=cnt+1 ! about to print an element in izmat, so move 1 step down
!       write (9,*) i,izmat(MXEL+1-cnt),zmat(MXEL+1-cnt)
!       write (9,*) izmat(MXEL+1-cnt),i,zmat(MXEL+1-cnt)
!45   continue
!44 continue
!   close(unit=9)

!    open (unit=9,file='matlab_imag_new.mtxx',status='unknown', access='sequential')
!    cnt=0
!    do 441 i=1,ndimo
!     write (9,*) i,i,zdiag(2,i) ! get the diagonal element out first
!     do 451 j=jzmat(i),(jzmat(i+1)-1)
!       cnt=cnt+1 ! about to print an element in izmat, so move 1 step down
!       write (9,*) i,izmat(cnt),zmat(cnt)
!       write (9,*) izmat(cnt),i,zmat(cnt)
!451   continue
!441 continue
!    close(unit=9)

!    write(*,*)  gxx, gyy, gzz, axx, ayy, azz, dx, dy, dz, pml, pmxy, pmzz, djf, djfprp, oss, psi, &
!                        ald, bed, gad, alm, bem, gam, c20, c22, c40, c42, &
!                        c44, t2edi, t2ndi, t2efi, t1edi, t1ndi, b0, a0, g0, pl,pkxy, pkzz

    return
    end subroutine matrxo


text:
function flin ( n, jsearch, l, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc )

!*****************************************************************************80
!
!! FLIN is the function of one variable to be minimized by MINNY.
!
!  Discussion:
!
!    F(X) is a scalar function of a vector argument X.
!
!    A minimizer of F(X) is sought along a line or parabola.
!
!    This function has been modified, by removing the occurrence of a
!    common block, so that it looks more like a "normal" function that does
!    not rely so strongly on peculiarities of FORTRAN.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    30 July 2016
!
!  Author:
!
!    Original FORTRAN77 version by Richard Brent.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization with Derivatives,
!    Prentice Hall, 1973,
!    Reprinted by Dover, 2002.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of variables.
!
!    Input, integer ( kind = 4 ) JSEARCH, indicates the kind of search.
!    If J is a legal column index, linear search in direction of V(*,JSEARCH).
!    Otherwise, then the search is parabolic, based on X, Q0 and Q1.
!
!    Input, real ( kind = 8 ) L, is the parameter determining the particular
!    point at which F is to be evaluated.
!    For a linear search, L is the step size.
!    For a quadratic search, L is a parameter which specifies
!    a point in the plane of X, Q0 and Q1.
!
!    Input, external F, is the name of the function to be minimized.
!    The function should have the form
!      function f(x,n)
!      integer ( kind = 4 ) n
!      real ( kind = 8 ) f
!      real ( kind = 8 ) x(n)
!    and accepts X and N as input, returning in F the function value.
!
!    Input, real ( kind = 8 ) X(N), the base point of the search.
!
!    Input/output, integer ( kind = 4 ) NF, the function evaluation counter.
!
!    Input, real ( kind = 8 ) V(N,N), a matrix whose columns constitute
!    search directions.
!
!    Input, real ( kind = 8 ) Q0(N), Q1(N), two auxiliary points used to
!    determine the plane when a quadratic search is performed.
!
!    Input, real ( kind = 8 ) QD0, QD1, values needed to compute the
!    coefficients QA, QB, QC.
!
!    Output, real ( kind = 8 ) QA, QB, QC, coefficients used to combine
!    Q0, X, and A1 if a quadratic search is used.
!
!    Output, real ( kind = 8 ) FLIN, the value of the function at the
!    minimizing point.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ), external :: f
  real ( kind = 8 ) flin
  integer ( kind = 4 ) i
  integer ( kind = 4 ) jsearch
  real ( kind = 8 ) l
  integer ( kind = 4 ) nf
  real ( kind = 8 ) q0(n)
  real ( kind = 8 ) q1(n)
  real ( kind = 8 ) qa
  real ( kind = 8 ) qb
  real ( kind = 8 ) qc
  real ( kind = 8 ) qd0
  real ( kind = 8 ) qd1
  real ( kind = 8 ) t(n)
  real ( kind = 8 ) v(n,n)
  real ( kind = 8 ) x(n)
!
!  The search is linear.
!
  if ( 1 <= jsearch ) then

    t(1:n) = x(1:n) + l * v(1:n,jsearch)
!
!  The search is along a parabolic space curve.
!
  else

    qa =                 l * ( l - qd1 ) /       ( qd0 + qd1 ) / qd0
    qb = - ( l + qd0 ) *     ( l - qd1 ) / qd1                 / qd0
    qc =   ( l + qd0 ) * l               / qd1 / ( qd0 + qd1 )

    t(1:n) = qa * q0(1:n) + qb * x(1:n) + qc * q1(1:n)

  end if
!
!  The function evaluation counter NF is incremented.
!
  nf = nf + 1
!
!  Evaluate the function.
!
  flin = f ( t, n )

  return
end
subroutine minfit ( n, tol, a, q )

!*****************************************************************************80
!
!! MINFIT computes the singular value decomposition of an N by N array.
!
!  Discussion:
!
!    This is an improved version of the EISPACK routine MINFIT
!    restricted to the case M = N and P = 0.
!
!    The singular values of the array A are returned in Q.  A is
!    overwritten with the orthogonal matrix V such that U * diag(Q) = A * V,
!    where U is another orthogonal matrix.
!
!    Thanks to Andreas Zuend for pointing out a potential for overflow in
!    the computation z = sqrt ( f*f + 1 ), 22 March 2012.
!
!    Thanks to Martin Horvat for correcting the initial assignment of S,
!    and for moving the line "L=I" back within the loop, 25 June 2013.
!
!    Thanks to Martin Horvat for moving the assignment "X=MAX(X,Y)"
!    into the loop, 29 October 2013.
!
!    Further modifications to this function were made, to remove the use
!    of statement labels, so that the function looks more like a "normal"
!    function that relies much less on the peculiarities of FORTRAN,
!    27 July 2016.
!
!    Removed argument M.  A must be dimensions (N,N).
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    30 July 2016
!
!  Author:
!
!    Original FORTRAN77 version by Richard Brent.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization with Derivatives,
!    Prentice Hall, 1973,
!    Reprinted by Dover, 2002.
!
!    James Wilkinson, Christian Reinsch,
!    Handbook for Automatic Computation,
!    Volume II, Linear Algebra, Part 2,
!    Springer Verlag, 1971.
!
!    Brian Smith, James Boyle, Jack Dongarra, Burton Garbow, Yasuhiko Ikebe,
!    Virginia Klema, Cleve Moler,
!    Matrix Eigensystem Routines, EISPACK Guide,
!    Lecture Notes in Computer Science, Volume 6,
!    Springer Verlag, 1976,
!    ISBN13: 978-3540075462,
!    LC: QA193.M37.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the order of the matrix A.
!
!    Input, real ( kind = 8 ) TOL, a tolerance which determines when a vector
!    (a column or part of a column of the matrix) may be considered
!    "essentially" equal to zero.
!
!    Input/output, real ( kind = 8 ) A(N,N).  On input, an N by N array whose
!    singular value decomposition is desired.  On output, the
!    SVD orthogonal matrix factor V.
!
!    Input/output, real ( kind = 8 ) Q(N), the singular values.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) a(n,n)
  real ( kind = 8 ) c
  real ( kind = 8 ) e(n)
  real ( kind = 8 ) eps
  real ( kind = 8 ) f
  real ( kind = 8 ) g
  real ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ii
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) kk
  integer ( kind = 4 ) kt
  integer ( kind = 4 ), parameter :: kt_max = 30
  integer ( kind = 4 ) l
  integer ( kind = 4 ) l2
  integer ( kind = 4 ) ll2
  real ( kind = 8 ) q(n)
  real ( kind = 8 ) r8_hypot
  real ( kind = 8 ) s
  logical skip
  real ( kind = 8 ) temp
  real ( kind = 8 ) tol
  real ( kind = 8 ) x
  real ( kind = 8 ) y
  real ( kind = 8 ) z
!
!  Householder's reduction to bidiagonal form.
!
  if ( n == 1 ) then
    q(1) = a(1,1)
    a(1,1) = 1.0D+00
    return
  end if

  eps = epsilon ( eps )
  g = 0.0D+00
  x = 0.0D+00

  do i = 1, n

    e(i) = g
    l = i + 1

    s = sum ( a(i:n,i) ** 2 )

    g = 0.0D+00

    if ( tol <= s ) then

      f = a(i,i)

      g = sqrt ( s )
      if ( 0.0D+00 <= f ) then
        g = - g
      end if

      h = f * g - s
      a(i,i) = f - g

      do j = l, n

        f = dot_product ( a(i:n,i), a(i:n,j) ) / h

        a(i:n,j) = a(i:n,j) + f * a(i:n,i)

      end do

    end if

    q(i) = g

    s = sum ( a(i,l:n) ** 2 )

    g = 0.0D+00

    if ( tol <= s ) then

      if ( i /= n ) then
        f = a(i,i+1)
      end if

      g = sqrt ( s )
      if ( 0.0D+00 <= f ) then
        g = - g
      end if

      h = f * g - s

      if ( i /= n ) then

        a(i,i+1) = f - g
        e(l:n) = a(i,l:n) / h

        do j = l, n

          s = dot_product ( a(j,l:n), a(i,l:n) )

          a(j,l:n) = a(j,l:n) + s * e(l:n)

        end do

      end if

    end if

    y = abs ( q(i) ) + abs ( e(i) )

    x = max ( x, y )

  end do
!
!  Accumulation of right-hand transformations.
!
  a(n,n) = 1.0D+00
  g = e(n)
  l = n

  do i = n - 1, 1, -1

    if ( g /= 0.0D+00 ) then

      h = a(i,i+1) * g

      a(l:n,i) = a(i,l:n) / h

      do j = l, n

        s = dot_product ( a(i,l:n), a(l:n,j) )

        a(l:n,j) = a(l:n,j) + s * a(l:n,i)

      end do

    end if

    a(i,l:n) = 0.0D+00
    a(l:n,i) = 0.0D+00
    a(i,i) = 1.0D+00

    g = e(i)

    l = i

  end do
!
!  Diagonalization of the bidiagonal form.
!
  eps = eps * x

  do k = n, 1, -1

    kt = 0

    do

      kt = kt + 1

      if ( kt_max < kt ) then
        e(k) = 0.0D+00
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'MINFIT - Fatal error!'
        write ( *, '(a)' ) '  The QR algorithm failed to converge.'
        stop 1
      end if

      skip = .false.

      do l2 = k, 1, -1

        l = l2

        if ( abs ( e(l) ) <= eps ) then
          skip = .true.
          exit
        end if

        if ( l /= 1 ) then
          if ( abs ( q(l-1) ) <= eps ) then
            exit
          end if
        end if

      end do
!
!  Cancellation of E(L) if 1 < L.
!
      if ( .not. skip ) then

        c = 0.0D+00
        s = 1.0D+00

        do i = l, k

          f = s * e(i)
          e(i) = c * e(i)
          if ( abs ( f ) <= eps ) then
            exit
          end if
          g = q(i)
!
!  q(i) = h = sqrt(g*g + f*f).
!
          h = r8_hypot ( f, g )

          q(i) = h

          if ( h == 0.0D+00 ) then
            g = 1.0D+00
            h = 1.0D+00
          end if

          c =   g / h
          s = - f / h

        end do

      end if
!
!  Test for convergence for this index K.
!
      z = q(k)

      if ( l == k ) then
        if ( z < 0.0D+00 ) then
          q(k) = - z
          a(1:n,k) = - a(1:n,k)
        end if
        exit
      end if
!
!  Shift from bottom 2*2 minor.
!
      x = q(l)
      y = q(k-1)
      g = e(k-1)
      h = e(k)
      f = ( ( y - z ) * ( y + z ) + ( g - h ) * ( g + h ) ) &
        / ( 2.0D+00 * h * y )

      g = r8_hypot ( f, 1.0D+00 )

      if ( f < 0.0D+00 ) then
        temp = f - g
      else
        temp = f + g
      end if

      f = ( ( x - z ) * ( x + z ) + h * ( y / temp - h ) ) / x
!
!  Next QR transformation.
!
      c = 1.0D+00
      s = 1.0D+00

      do i = l + 1, k

        g = e(i)
        y = q(i)
        h = s * g
        g = g * c

        z = r8_hypot ( f, h )

        e(i-1) = z

        if ( z == 0.0D+00 ) then
          f = 1.0D+00
          z = 1.0D+00
        end if

        c = f / z
        s = h / z
        f =   x * c + g * s
        g = - x * s + g * c
        h = y * s
        y = y * c

        do j = 1, n
          x = a(j,i-1)
          z = a(j,i)
          a(j,i-1) = x * c + z * s
          a(j,i) = - x * s + z * c
        end do

        z = r8_hypot ( f, h )

        q(i-1) = z

        if ( z == 0.0D+00 ) then
          f = 1.0D+00
          z = 1.0D+00
        end if

        c = f / z
        s = h / z
        f =   c * g + s * y
        x = - s * g + c * y

      end do

      e(l) = 0.0D+00
      e(k) = f
      q(k) = x

    end do

  end do

  return
end
subroutine minny ( n, jsearch, nits, d2, x1, f1, fk, f, x, t, h, v, q0, q1, &
  nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 )

!*****************************************************************************80
!
!! MINNY minimizes a scalar function of N variables along a line.
!
!  Discussion:
!
!    MINNY minimizes F along the line from X in the direction V(*,J) unless
!    J is less than 1, when a quadratic search is made in the plane
!    defined by Q0, Q1 and X.
!
!    If FK = true, then F1 is FLIN(X1).  Otherwise X1 and F1 are ignored
!    on entry unless final FX is greater than F1.
!
!    This function was modified by removing the common blocks
!    and the use of labeled statements, 28 July 2016.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    30 July 2016
!
!  Author:
!
!    Original FORTRAN77 version by Richard Brent.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization with Derivatives,
!    Prentice Hall, 1973,
!    Reprinted by Dover, 2002.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of variables.
!
!    Input, integer ( kind = 4 ) JSEARCH, indicates the kind of search.
!    If JSEARCH is a legal column index, linear search in direction of V(*,J).
!    Otherwise, the search is parabolic, based on X, Q0 and Q1.
!
!    Input, integer ( kind = 4 ) NITS, the maximum number of times the interval
!    may be halved to retry the calculation.
!
!    Input/output, real ( kind = 8 ) D2, is either zero, or an approximation to
!    the value of (1/2) times the second derivative of F.
!
!    Input/output, real ( kind = 8 ) X1, on entry, an estimate of the
!    distance from X to the minimum along V(*,J), or, if J = 0, a curve.
!    On output, the distance between X and the minimizer that was found.
!
!    Input/output, real ( kind = 8 ) F1, ?
!
!    Input, logical FK; if FK is TRUE, then on input F1 contains
!    the value FLIN(X1).
!
!    Input, external real ( kind = 8 ) F, is the name of the function to
!    be minimized.  The function should have the form
!      function f(x,n)
!      integer ( kind = 4 ) n
!      real ( kind = 8 ) f
!      real ( kind = 8 ) x(n)
!    and accepts X and N as input, returning in F the function value.
!
!    Input/output, real ( kind = 8 ) X(N), ?
!
!    Input, real ( kind = 8 ) T, ?
!
!    Input, real ( kind = 8 ) H, ?
!
!    Input, real ( kind = 8 ) V(N,N), a matrix whose columns are direction
!    vectors along which the function may be minimized.
!
!    Input, real ( kind = 8 ) Q0(N), an auxiliary point used to define
!    a curve through X.
!
!    Input, real ( kind = 8 ) Q1(N), an auxiliary point used to define
!    a curve through X.
!
!    Input/output, integer ( kind = 4 ) NL, the number of linear searches.
!
!    Input/output, integer ( kind = 4 ) NF, the number of function evaluations.
!
!    Input, real ( kind = 8 ) DMIN, an estimate for the smallest eigenvalue.
!
!    Input, real ( kind = 8 ) LDT, the length of the step.
!
!    Input/output, real ( kind = 8 ) FX, the value of F(X,N).
!
!    Input/output, real ( kind = 8 ) QA, QB, QC, ?
!
!    Input, real ( kind = 8 ) QD0, QD1, ?.
!
  implicit none

  real ( kind = 8 ) d1
  real ( kind = 8 ) d2
  real ( kind = 8 ) dmin
  logical dz
  real ( kind = 8 ), external :: f
  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) f2
  logical fk
  real ( kind = 8 ) flin
  real ( kind = 8 ) fm
  real ( kind = 8 ) fx
  real ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ) jsearch
  integer ( kind = 4 ) k
  real ( kind = 8 ) ldt
  real ( kind = 8 ) m2
  real ( kind = 8 ) m4
  real ( kind = 8 ) machep
  integer ( kind = 4 ) n
  integer ( kind = 4 ) nf
  integer ( kind = 4 ) nits
  integer ( kind = 4 ) nl
  logical ok
  real ( kind = 8 ) q0(n)
  real ( kind = 8 ) q1(n)
  real ( kind = 8 ) qa
  real ( kind = 8 ) qb
  real ( kind = 8 ) qc
  real ( kind = 8 ) qd0
  real ( kind = 8 ) qd1
  real ( kind = 8 ) s
  real ( kind = 8 ) sf1
  real ( kind = 8 ) small
  real ( kind = 8 ) sx1
  real ( kind = 8 ) t
  real ( kind = 8 ) t2
  real ( kind = 8 ) temp
  real ( kind = 8 ) v(n,n)
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) x1
  real ( kind = 8 ) x2
  real ( kind = 8 ) xm

  machep = epsilon ( machep )
  small = machep ** 2
  m2 = sqrt ( machep )
  m4 = sqrt ( m2 )
  sf1 = f1
  sx1 = x1
  k = 0
  xm = 0.0D+00
  fm = fx
  f0 = fx
  dz = ( d2 < machep )
!
!  Find the step size.
!
  s = sqrt ( sum ( x(1:n) ** 2 ) )

  if ( dz ) then
    temp = dmin
  else
    temp = d2
  end if

  t2 = m4 * sqrt ( abs ( fx ) / temp + s * ldt ) + m2 * ldt
  s = m4 * s + t
  if ( dz .and. s < t2 ) then
    t2 = s
  end if

  t2 = max ( t2, small )
  t2 = min ( t2, 0.01D+00 * h )

  if ( fk .and. f1 <= fm ) then
    xm = x1
    fm = f1
  end if

  if ( .not. fk .or. abs ( x1 ) < t2 ) then

    if ( 0.0D+00 <= x1 ) then
      temp = 1.0D+00
    else
      temp = - 1.0D+00
    end if

    x1 = temp * t2
    f1 = flin ( n, jsearch, x1, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc )

  end if

  if ( f1 <= fm ) then
    xm = x1
    fm = f1
  end if
!
!  Evaluate FLIN at another point and estimate the second derivative.
!
  do

    if ( dz ) then

      if ( f1 <= f0 ) then
        x2 = 2.0D+00 * x1
      else
        x2 = - x1
      end if

      f2 = flin ( n, jsearch, x2, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc )

      if ( f2 <= fm ) then
        xm = x2
        fm = f2
      end if

      d2 = ( x2 * ( f1 - f0 ) - x1 * ( f2 - f0 ) ) &
        / ( ( x1 * x2 ) * ( x1 - x2 ) )

    end if
!
!  Estimate the first derivative at 0.
!
    d1 = ( f1 - f0 ) / x1 - x1 * d2
    dz = .true.
!
!  Predict the minimum.
!
    if ( d2 <= small ) then

      if ( 0.0D+00 <= d1 ) then
        x2 = - h
      else
        x2 = h
      end if

    else

      x2 = ( - 0.5D+00 * d1 ) / d2

    end if

    if ( h < abs ( x2 ) ) then

      if ( x2 <= 0.0D+00 ) then
        x2 = - h
      else
        x2 = h
      end if

    end if
!
!  Evaluate F at the predicted minimum.
!
    ok = .true.

    do

      f2 = flin ( n, jsearch, x2, f, x, nf, v, q0, q1, qd0, qd1, qa, qb, qc )

      if ( nits <= k .or. f2 <= f0 ) then
        exit
      end if

      k = k + 1

      if ( f0 < f1 .and. 0.0D+00 < x1 * x2 ) then
        ok = .false.
        exit
      end if

      x2 = 0.5D+00 * x2

    end do

    if ( ok ) then
      exit
    end if

  end do
!
!  Increment the one-dimensional search counter.
!
  nl = nl + 1

  if ( fm < f2 ) then
    x2 = xm
  else
    fm = f2
  end if
!
!  Get a new estimate of the second derivative.
!
  if ( small < abs ( x2 * ( x2 - x1 ) ) ) then
    d2 = ( x2 * ( f1 - f0 ) - x1 * ( fm - f0 ) ) / ( ( x1 * x2 ) * ( x1 - x2 ) )
  else
    if ( 0 < k ) then
      d2 = 0.0D+00
    end if
  end if

  d2 = max ( d2, small )

  x1 = x2
  fx = fm

  if ( sf1 < fx ) then
    fx = sf1
    x1 = sx1
  end if
!
!  Update X for linear but not parabolic search.
!
  if ( 1 <= jsearch ) then

    x(1:n) = x(1:n) + x1 * v(1:n,jsearch)

  end if

  return
end
function praxis ( t0, h0, n, prin, x, f )

!*****************************************************************************80
!
!! PRAXIS seeks an N-dimensional minimizer X of a scalar function F(X).
!
!  Discussion:
!
!    PRAXIS returns the minimum of the function F(X,N) of N variables
!    using the principal axis method.  The gradient of the function is
!    not required.
!
!    The approximating quadratic form is
!
!      Q(x') = F(x,n) + (1/2) * (x'-x)' * A * (x'-x)
!
!    where X is the best estimate of the minimum and
!
!      A = inverse(V') * D * inverse(V)
!
!    V(*,*) is the matrix of search directions;
!    D(*) is the array of second differences.
!
!    If F(X) has continuous second derivatives near X0, then A will tend
!    to the hessian of F at X0 as X approaches X0.
!
!    Thanks to Andreas Zuend for pointing out an error in the form of the
!    call to the routine r8mat_print (), 22 March 2012.
!
!    This function was modified by eliminating the use of labeled statements,
!    and removing the common blocks, 28 July 2016.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    28 July 2016
!
!  Author:
!
!    Original FORTRAN77 version by Richard Brent.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization with Derivatives,
!    Prentice Hall, 1973,
!    Reprinted by Dover, 2002.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) T0, is a tolerance.  PRAXIS attempts to return
!    praxis = f(x) such that if X0 is the true local minimum near X, then
!    norm ( x - x0 ) < T0 + sqrt ( EPSILON ( X ) ) * norm ( X ),
!    where EPSILON ( X ) is the machine precision for X.
!
!    Input, real ( kind = 8 ) H0, is the maximum step size.  H0 should be
!    set to about the maximum distance from the initial guess to the minimum.
!    If H0 is set too large or too small, the initial rate of
!    convergence may be slow.
!
!    Input, integer ( kind = 4 ) N, the number of variables.
!
!    Input, integer ( kind = 4 ) PRIN, controls printing intermediate results.
!    0, nothing is printed.
!    1, F is printed after every n+1 or n+2 linear minimizations.
!       final X is printed, but intermediate X is printed only
!       if N is at most 4.
!    2, the scale factors and the principal values of the approximating
!       quadratic form are also printed.
!    3, X is also printed after every few linear minimizations.
!    4, the principal vectors of the approximating quadratic form are
!       also printed.
!
!    Input/output, real ( kind = 8 ) X(N), is an array containing on entry a
!    guess of the point of minimum, on return the estimated point of minimum.
!
!    Input, external real ( kind = 8 ) F, is the name of the function to be
!    minimized.  The function should have the form
!      function f(x,n)
!      integer ( kind = 4 ) n
!      real ( kind = 8 ) f
!      real ( kind = 8 ) x(n)
!    and accepts X and N as input, returning in F the function value.
!
!    Output, real ( kind = 8 ) PRAXIS, the function value at the minimizer.
!
!  Local parameters:
!
!    Local, real ( kind = 8 ) DMIN, an estimate for the smallest eigenvalue.
!
!    Local, real ( kind = 8 ) FX, the value of F(X,N).
!
!    Local, logical ILLC, is TRUE if the system is ill-conditioned.
!
!    Local, real ( kind = 8 ) LDT, the length of the step.
!
!    Local, integer ( kind = 4 ) NF, the number of function evaluations.
!
!    Local, integer ( kind = 4 ) NL, the number of linear searches.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) d(n)
  real ( kind = 8 ) d2
  real ( kind = 8 ) df
  real ( kind = 8 ) dmin
  real ( kind = 8 ) dn
  real ( kind = 8 ) dni
  real ( kind = 8 ), external :: f
  real ( kind = 8 ) f1
  logical fk
  real ( kind = 8 ) fx
  real ( kind = 8 ) h
  real ( kind = 8 ) h0
  integer ( kind = 4 ) i
  integer ( kind = 4 ) ii
  logical illc
  integer ( kind = 4 ) j
  integer ( kind = 4 ) jsearch
  integer ( kind = 4 ) k
  integer ( kind = 4 ) k2
  integer ( kind = 4 ) kl
  integer ( kind = 4 ) kt
  integer ( kind = 4 ) ktm
  real ( kind = 8 ) large
  real ( kind = 8 ) ldfac
  real ( kind = 8 ) lds
  real ( kind = 8 ) ldt
  real ( kind = 8 ) m2
  real ( kind = 8 ) m4
  real ( kind = 8 ) machep
  integer ( kind = 4 ) nits
  integer ( kind = 4 ) nl
  integer ( kind = 4 ) nf
  real ( kind = 8 ) praxis
  integer ( kind = 4 ) prin
  real ( kind = 8 ) q0(n)
  real ( kind = 8 ) q1(n)
  real ( kind = 8 ) qa
  real ( kind = 8 ) qb
  real ( kind = 8 ) qc
  real ( kind = 8 ) qd0
  real ( kind = 8 ) qd1
  real ( kind = 8 ) qf1
  real ( kind = 8 ) r
  real ( kind = 8 ) r8_uniform_01
  real ( kind = 8 ) s
  real ( kind = 8 ) scbd
  integer ( kind = 4 ) seed
  real ( kind = 8 ) sf
  real ( kind = 8 ) sl
  real ( kind = 8 ) small
  real ( kind = 8 ) t
  real ( kind = 8 ) t0
  real ( kind = 8 ) t2
  real ( kind = 8 ) v(n,n)
  real ( kind = 8 ) value
  real ( kind = 8 ) vlarge
  real ( kind = 8 ) vsmall
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) y(n)
  real ( kind = 8 ) z(n)
!
!  Initialization.
!
  machep = epsilon ( machep )
  small = machep * machep
  vsmall = small * small
  large = 1.0D+00 / small
  vlarge = 1.0D+00 / vsmall
  m2 = sqrt ( machep )
  m4 = sqrt ( m2 )
  seed = 123456789
!
!  Heuristic numbers:
!
!  If the axes may be badly scaled (which is to be avoided if
!  possible), then set SCBD = 10.  Otherwise set SCBD = 1.
!
!  If the problem is known to be ill-conditioned, initialize ILLC = true.
!
!  KTM is the number of iterations without improvement before the
!  algorithm terminates.  KTM = 4 is very cautious; usually KTM = 1
!  is satisfactory.
!
  scbd = 1.0D+00
  illc = .false.
  ktm = 1

  if ( illc ) then
    ldfac = 0.1D+00
  else
    ldfac = 0.01D+00
  end if

  kt = 0
  nl = 0
  nf = 1
  fx = f ( x, n )
  qf1 = fx
  t = small + abs ( t0 )
  t2 = t
  dmin = small
  h = h0
  h = max ( h, 100.0D+00 * t )
  ldt = h
!
!  The initial set of search directions V is the identity matrix.
!
  v(1:n,1:n) = 0.0D+00
  do i = 1, n
    v(i,i) = 1.0D+00
  end do

  d(1:n) = 0.0D+00
  qa = 0.0D+00
  qb = 0.0D+00
  qc = 0.0D+00
  qd0 = 0.0D+00
  qd1 = 0.0D+00
  q0(1:n) = x(1:n)
  q1(1:n) = x(1:n)

  if ( 0 < prin ) then
    call print2 ( n, x, prin, fx, nf, nl )
  end if
!
!  The main loop starts here.
!
  do

    sf = d(1)
    d(1) = 0.0D+00
!
!  Minimize along the first direction V(*,1).
!
    jsearch = 1
    nits = 2
    d2 = d(1)
    s = 0.0D+00
    value = fx
    fk = .false.

    call minny ( n, jsearch, nits, d2, s, value, fk, f, x, t, &
      h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 )

    d(1) = d2

    if ( s <= 0.0D+00 ) then
      v(1:n,1) = - v(1:n,1)
    end if

    if ( sf <= 0.9D+00 * d(1) .or. d(1) <= 0.9D+00 * sf ) then
      d(2:n) = 0.0D+00
    end if
!
!  The inner loop starts here.
!
    do k = 2, n

      y(1:n) = x(1:n)

      sf = fx

      if ( 0 < kt ) then
        illc = .true.
      end if

      do

        kl = k
        df = 0.0D+00
!
!  A random step follows, to avoid resolution valleys.
!
        if ( illc ) then

          do i = 1, n
            r = r8_uniform_01 ( seed )
!           call random_number ( harvest = r )
            s = ( 0.1D+00 * ldt + t2 * 10.0D+00 ** kt ) * ( r - 0.5D+00 )
            z(i) = s
            x(1:n) = x(1:n) + s * v(1:n,i)
          end do

          fx = f ( x, n )
          nf = nf + 1

        end if
!
!  Minimize along the "non-conjugate" directions V(*,K),...,V(*,N).
!
        do k2 = k, n

          sl = fx

          jsearch = k2
          nits = 2
          d2 = d(k2)
          s = 0.0D+00
          value = fx
          fk = .false.

          call minny ( n, jsearch, nits, d2, s, value, fk, f, x, t, &
            h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 )

          d(k2) = d2

          if ( illc ) then
            s = d(k2) * ( ( s + z(k2) ) ** 2 )
          else
            s = sl - fx
          end if

          if ( df <= s ) then
            df = s
            kl = k2
          end if

        end do
!
!  If there was not much improvement on the first try, set
!  ILLC = true and start the inner loop again.
!
        if ( illc ) then
          exit
        end if

        if ( abs ( 100.0D+00 * machep * fx ) <= df ) then
          exit
        end if

        illc = .true.

      end do

      if ( k == 2 .and. 1 < prin ) then
        call r8vec_print ( n, d, '  The second difference array:' )
      end if
!
!  Minimize along the "conjugate" directions V(*,1),...,V(*,K-1).
!
      do k2 = 1, k - 1

        jsearch = k2
        nits = 2
        d2 = d(k2)
        s = 0.0D+00
        value = fx
        fk = .false.

        call minny ( n, jsearch, nits, d2, s, value, fk, f, x, t, &
          h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 )

        d(k2) = d2

      end do

      f1 = fx
      fx = sf
      lds = 0.0D+00

      do i = 1, n
        sl = x(i)
        x(i) = y(i)
        sl = sl - y(i)
        y(i) = sl
        lds = lds + sl ** 2
      end do

      lds = sqrt ( lds )
!
!  Discard direction V(*,kl).
!
!  If no random step was taken, V(*,KL) is the "non-conjugate"
!  direction along which the greatest improvement was made.
!
      if ( small < lds ) then

        do j = kl - 1, k, -1
          v(1:n,j+1) = v(1:n,j)
          d(j+1) = d(j)
        end do

        d(k) = 0.0D+00

        v(1:n,k) = y(1:n) / lds
!
!  Minimize along the new "conjugate" direction V(*,k), which is
!  the normalized vector:  (new x) - (old x).
!
        jsearch = k
        nits = 4
        d2 = d(k)
        value = f1
        fk = .true.

        call minny ( n, jsearch, nits, d2, lds, value, fk, f, x, t, &
          h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 )

        d(k) = d2

        if ( lds <= 0.0D+00 ) then
          lds = - lds
          v(1:n,k) = - v(1:n,k)
        end if

      end if

      ldt = ldfac * ldt
      ldt = max ( ldt, lds )

      if ( 0 < prin ) then
        call print2 ( n, x, prin, fx, nf, nl )
      end if

      t2 = m2 * sqrt ( sum ( x(1:n) ** 2 ) ) + t
!
!  See whether the length of the step taken since starting the
!  inner loop exceeds half the tolerance.
!
      if ( 0.5D+00 * t2 < ldt ) then
        kt = - 1
      end if

      kt = kt + 1

      if ( ktm < kt ) then

        if ( 0 < prin ) then
          call r8vec_print ( n, x, '  X:' )
        end if

        praxis = fx

        return

      end if

    end do
!
!  The inner loop ends here.
!
!  Try quadratic extrapolation in case we are in a curved valley.
!
    call quad ( n, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, fx, qf1, &
      qa, qb, qc, qd0, qd1 )

    d(1:n) = 1.0D+00 / sqrt ( d(1:n) )

    dn = maxval ( d(1:n) )

    if ( 3 < prin ) then
      call r8mat_print ( n, n, v, '  The new direction vectors:' )
    end if

    do j = 1, n
      v(1:n,j) = ( d(j) / dn ) * v(1:n,j)
    end do
!
!  Scale the axes to try to reduce the condition number.
!
    if ( 1.0D+00 < scbd ) then

      do i = 1, n
        z(i) = max ( m4, sqrt ( sum ( v(i,1:n) ** 2 ) ) )
      end do

      s = minval ( z(1:n) )

      do i = 1, n

        sl = s / z(i)
        z(i) = 1.0D+00 / sl

        if ( scbd < z(i) ) then
          sl = 1.0D+00 / scbd
          z(i) = scbd
        end if

        v(i,1:n) = sl * v(i,1:n)

      end do

    end if
!
!  Calculate a new set of orthogonal directions before repeating
!  the main loop.
!
!  Transpose V for MINFIT:
!
    v(1:n,1:n) = transpose ( v(1:n,1:n) )
!
!  Call MINFIT to find the singular value decomposition of V.
!
!  This gives the principal values and principal directions of the
!  approximating quadratic form without squaring the condition number.
!
    call minfit ( n, vsmall, v, d )
!
!  Unscale the axes.
!
    if ( 1.0D+00 < scbd ) then

      do i = 1, n
        v(i,1:n) = z(i) * v(i,1:n)
      end do

      do i = 1, n

        s = sqrt ( sum ( v(1:n,i) ** 2 ) )
        d(i) = s * d(i)
        v(1:n,i) = v(1:n,i) / s

      end do

    end if

    do i = 1, n

      dni = dn * d(i)

      if ( large < dni ) then
        d(i) = vsmall
      else if ( dni < small ) then
        d(i) = vlarge
      else
        d(i) = 1.0D+00 / dni ** 2
      end if

    end do
!
!  Sort the singular values and singular vectors.
!
    call svsort ( n, d, v )
!
!  Determine the smallest eigenvalue.
!
    dmin = max ( d(n), small )
!
!  The ratio of the smallest to largest eigenvalue determines whether
!  the system is ill conditioned.
!
    if ( dmin < m2 * d(1) ) then
      illc = .true.
    else
      illc = .false.
    end if

    if ( 1 < prin ) then

      if ( 1.0D+00 < scbd ) then
        call r8vec_print ( n, z, '  The scale factors:' )
      end if

      call r8vec_print ( n, d, '  Principal values of the quadratic form:' )

    end if

    if ( 3 < prin ) then
      call r8mat_print ( n, n, v, '  The principal axes:' )
    end if
!
!  The main loop ends here.
!
  end do

  if ( 0 < prin ) then
    call r8vec_print ( n, x, '  X:' )
  end if

  praxis = fx

  return
end
subroutine print2 ( n, x, prin, fx, nf, nl )

!*****************************************************************************80
!
!! PRINT2 prints certain data about the progress of the iteration.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    22 May 2006
!
!  Author:
!
!    Original FORTRAN77 version by Richard Brent.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization with Derivatives,
!    Prentice Hall, 1973,
!    Reprinted by Dover, 2002.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of variables.
!
!    Input, real ( kind = 8 ) X(N), the current estimate of the minimizer.
!
!    Input, integer ( kind = 4 ) PRIN, the user-specifed print level.
!    0, nothing is printed.
!    1, F is printed after every n+1 or n+2 linear minimizations.
!       final X is printed, but intermediate X is printed only
!       if N is at most 4.
!    2, the scale factors and the principal values of the approximating
!       quadratic form are also printed.
!    3, X is also printed after every few linear minimizations.
!    4, the principal vectors of the approximating quadratic form are
!       also printed.
!
!    Input, real ( kind = 8 ) FX, the smallest value of F(X) found so far.
!
!    Input, integer ( kind = 4 ) NF, the number of function evaluations.
!
!    Input, integer ( kind = 4 ) NL, the number of linear searches.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) fx
  integer ( kind = 4 ) nf
  integer ( kind = 4 ) nl
  integer ( kind = 4 ) prin
  real ( kind = 8 ) x(n)

  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) '  Linear searches      ', nl
  write ( *, '(a,i8)' ) '  Function evaluations ', nf

  if ( n <= 4 .or. 2 < prin ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'X:'
    write ( *, '(5g14.6)' ) x(1:n)
  end if

  write ( *, '(a)' ) ' '
  write ( *, '(a,g14.6)' ) '  The function value FX: ', fx

  return
end
subroutine quad ( n, f, x, t, h, v, q0, q1, nl, nf, dmin, ldt, fx, qf1, &
  qa, qb, qc, qd0, qd1 )

!*****************************************************************************80
!
!! QUAD seeks to minimize the scalar function F along a particular curve.
!
!  Discussion:
!
!    The minimizer to be sought is required to lie on a curve defined
!    by Q0, Q1 and X.
!
!    This function was modified by removing the common blocks,
!    28 July 2016.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    30 July 2016
!
!  Author:
!
!    Original FORTRAN77 version by Richard Brent.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization with Derivatives,
!    Prentice Hall, 1973,
!    Reprinted by Dover, 2002.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of variables.
!
!    Input, external real ( kind = 8 ) F, is the name of the function to
!    be minimized.  The function should have the form
!      function f(x,n)
!      integer ( kind = 4 ) n
!      real ( kind = 8 ) f
!      real ( kind = 8 ) x(n)
!    and accepts X and N as input, returning in F the function value.
!
!    Input/output, real ( kind = 8 ) X(N), ?
!
!    Input, real ( kind = 8 ) T, ?
!
!    Input, real ( kind = 8 ) H, ?
!
!    Input, real ( kind = 8 ) V(N,N), the matrix of search directions.
!
!    Input/output, real ( kind = 8 ) Q0(N), Q1(N), auxiliary points used to \
!    define a curve through X.
!
!    Input/output, integer ( kind = 4 ) NL, the number of linear searches.
!
!    Input/output, integer ( kind = 4 ) NF, the number of function evaluations.
!
!    Input, real ( kind = 8 ) DMIN, an estimate for the smallest eigenvalue.
!
!    Input, real ( kind = 8 ) LDT, the length of the step.
!
!    Input/output, real ( kind = 8 ) FX, the value of F(X,N).
!
!    Input/output, real ( kind = 8 ) QF1, QA, QB, QC, QD0, QD1 ?
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) dmin
  real ( kind = 8 ), external :: f
  logical fk
  real ( kind = 8 ) fx
  real ( kind = 8 ) h
  integer ( kind = 4 ) i
  integer ( kind = 4 ) jsearch
  real ( kind = 8 ) l
  real ( kind = 8 ) ldt
  integer ( kind = 4 ) nf
  integer ( kind = 4 ) nits
  integer ( kind = 4 ) nl
  real ( kind = 8 ) q0(n)
  real ( kind = 8 ) q1(n)
  real ( kind = 8 ) qa
  real ( kind = 8 ) qb
  real ( kind = 8 ) qc
  real ( kind = 8 ) qd0
  real ( kind = 8 ) qd1
  real ( kind = 8 ) qf1
  real ( kind = 8 ) s
  real ( kind = 8 ) t
  real ( kind = 8 ) temp
  real ( kind = 8 ) v(n,n)
  real ( kind = 8 ) value
  real ( kind = 8 ) x(n)

  temp = fx
  fx   = qf1
  qf1  = temp

  call r8vec_swap ( n, x, q1 )

  qd1 = sqrt ( sum ( ( x(1:n) - q1(1:n) ) ** 2 ) )

  l = qd1
  s = 0.0D+00

  if ( qd0 <= 0.0D+00 .or. qd1 <= 0.0D+00 .or. nl < 3 * n * n ) then

    fx = qf1
    qa = 0.0D+00
    qb = 0.0D+00
    qc = 1.0D+00

  else

    jsearch = -1
    nits = 2
    value = qf1
    fk = .true.

    call minny ( n, jsearch, nits, s, l, value, fk, f, x, t, &
      h, v, q0, q1, nl, nf, dmin, ldt, fx, qa, qb, qc, qd0, qd1 )

    qa =                 l * ( l - qd1 )       / ( qd0 + qd1 ) / qd0
    qb = - ( l + qd0 )     * ( l - qd1 ) / qd1                 / qd0
    qc =   ( l + qd0 ) * l               / qd1 / ( qd0 + qd1 )

  end if

  qd0 = qd1

  do i = 1, n
    s = q0(i)
    q0(i) = x(i)
    x(i) = qa * s + qb * x(i) + qc * q1(i)
  end do

  return
end
function r8_hypot ( x, y )

!*****************************************************************************80
!
!! R8_HYPOT returns the value of sqrt ( X^2 + Y^2).
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    22 March 2012
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X, Y, the arguments.
!
!    Output, real ( kind = 8 ) R8_HYPOT, the value of sqrt ( X^2 + Y^2 ).
!
  implicit none

  real ( kind = 8 ) a
  real ( kind = 8 ) b
  real ( kind = 8 ) c
  real ( kind = 8 ) r8_hypot
  real ( kind = 8 ) x
  real ( kind = 8 ) y

  if ( abs ( x ) < abs ( y ) ) then
    a = abs ( y )
    b = abs ( x )
  else
    a = abs ( x )
    b = abs ( y )
  end if
!
!  A contains the larger value.
!
  if ( a == 0.0D+00 ) then
    c = 0.0D+00
  else
    c = a * sqrt ( 1.0D+00 + ( b / a ) ** 2 )
  end if

  r8_hypot = c

  return
end
function r8_uniform_01 ( seed )

!*****************************************************************************80
!
!! R8_UNIFORM_01 returns a unit pseudorandom R8.
!
!  Discussion:
!
!    An R8 is a real ( kind = 8 ) value.
!
!    For now, the input quantity SEED is an integer variable.
!
!    This routine implements the recursion
!
!      seed = 16807 * seed mod ( 2^31 - 1 )
!      r8_uniform_01 = seed / ( 2^31 - 1 )
!
!    The integer arithmetic never requires more than 32 bits,
!    including a sign bit.
!
!    If the initial seed is 12345, then the first three computations are
!
!      Input     Output      R8_UNIFORM_01
!      SEED      SEED
!
!         12345   207482415  0.096616
!     207482415  1790989824  0.833995
!    1790989824  2035175616  0.947702
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    05 July 2006
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Paul Bratley, Bennett Fox, Linus Schrage,
!    A Guide to Simulation,
!    Springer Verlag, pages 201-202, 1983.
!
!    Pierre L'Ecuyer,
!    Random Number Generation,
!    in Handbook of Simulation,
!    edited by Jerry Banks,
!    Wiley Interscience, page 95, 1998.
!
!    Bennett Fox,
!    Algorithm 647:
!    Implementation and Relative Efficiency of Quasirandom
!    Sequence Generators,
!    ACM Transactions on Mathematical Software,
!    Volume 12, Number 4, pages 362-376, 1986.
!
!    Peter Lewis, Allen Goodman, James Miller
!    A Pseudo-Random Number Generator for the System/360,
!    IBM Systems Journal,
!    Volume 8, pages 136-143, 1969.
!
!  Parameters:
!
!    Input/output, integer ( kind = 4 ) SEED, the "seed" value, which should
!    NOT be 0. On output, SEED has been updated.
!
!    Output, real ( kind = 8 ) R8_UNIFORM_01, a new pseudorandom variate,
!    strictly between 0 and 1.
!
  implicit none

  integer ( kind = 4 ), parameter :: i4_huge = 2147483647
  integer ( kind = 4 ) k
  real ( kind = 8 ) r8_uniform_01
  integer ( kind = 4 ) seed

  if ( seed == 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'R8_UNIFORM_01 - Fatal error!'
    write ( *, '(a)' ) '  Input value of SEED = 0.'
    stop 1
  end if

  k = seed / 127773

  seed = 16807 * ( seed - k * 127773 ) - k * 2836

  if ( seed < 0 ) then
    seed = seed + i4_huge
  end if

  r8_uniform_01 = real ( seed, kind = 8 ) * 4.656612875D-10

  return
end
subroutine r8mat_print ( m, n, a, title )

!*****************************************************************************80
!
!! R8MAT_PRINT prints an R8MAT.
!
!  Discussion:
!
!    An R8MAT is an array of R8 values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    12 September 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, the number of rows in A.
!
!    Input, integer ( kind = 4 ) N, the number of columns in A.
!
!    Input, real ( kind = 8 ) A(M,N), the matrix.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none

  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) a(m,n)
  character ( len = * ) title

  call r8mat_print_some ( m, n, a, 1, 1, m, n, title )

  return
end
subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

!*****************************************************************************80
!
!! R8MAT_PRINT_SOME prints some of an R8MAT.
!
!  Discussion:
!
!    An R8MAT is an array of R8 values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    26 March 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) M, N, the number of rows and columns.
!
!    Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed.
!
!    Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print.
!
!    Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none

  integer ( kind = 4 ), parameter :: incx = 5
  integer ( kind = 4 ) m
  integer ( kind = 4 ) n

  real ( kind = 8 ) a(m,n)
  character ( len = 14 ) ctemp(incx)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) i2hi
  integer ( kind = 4 ) i2lo
  integer ( kind = 4 ) ihi
  integer ( kind = 4 ) ilo
  integer ( kind = 4 ) inc
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j2
  integer ( kind = 4 ) j2hi
  integer ( kind = 4 ) j2lo
  integer ( kind = 4 ) jhi
  integer ( kind = 4 ) jlo
  character ( len = * ) title

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( title )

  do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx

    j2hi = j2lo + incx - 1
    j2hi = min ( j2hi, n )
    j2hi = min ( j2hi, jhi )

    inc = j2hi + 1 - j2lo

    write ( *, '(a)' ) ' '

    do j = j2lo, j2hi
      j2 = j + 1 - j2lo
      write ( ctemp(j2), '(i8,6x)' ) j
    end do

    write ( *, '(''  Col   '',5a14)' ) ctemp(1:inc)
    write ( *, '(a)' ) '  Row'
    write ( *, '(a)' ) ' '

    i2lo = max ( ilo, 1 )
    i2hi = min ( ihi, m )

    do i = i2lo, i2hi

      do j2 = 1, inc

        j = j2lo - 1 + j2

        if ( a(i,j) == real ( int ( a(i,j) ), kind = 8 ) ) then
          write ( ctemp(j2), '(f8.0,6x)' ) a(i,j)
        else
          write ( ctemp(j2), '(g14.6)' ) a(i,j)
        end if

      end do

      write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j), j = 1, inc )

    end do

  end do

  return
end
subroutine r8vec_print ( n, a, title )

!*****************************************************************************80
!
!! R8VEC_PRINT prints an R8VEC.
!
!  Discussion:
!
!    An R8VEC is a vector of R8 values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    22 August 2000
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of components of the vector.
!
!    Input, real ( kind = 8 ) A(N), the vector to be printed.
!
!    Input, character ( len = * ) TITLE, a title.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) a(n)
  integer ( kind = 4 ) i
  character ( len = * ) title

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( title )
  write ( *, '(a)' ) ' '
  do i = 1, n
    write ( *, '(2x,i8,2x,g16.8)' ) i, a(i)
  end do

  return
end
subroutine r8vec_swap ( n, a1, a2 )

!*****************************************************************************80
!
!! R8VEC_SWAP swaps the entries of two R8VECs.
!
!  Discussion:
!
!    An R8VEC is a vector of R8 values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    04 December 2004
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the number of entries in the arrays.
!
!    Input/output, real ( kind = 8 ) A1(N), A2(N), the vectors to swap.
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) a1(n)
  real ( kind = 8 ) a2(n)
  real ( kind = 8 ) a3(n)

  a3(1:n) = a1(1:n)
  a1(1:n) = a2(1:n)
  a2(1:n) = a3(1:n)

  return
end
subroutine svsort ( n, d, v )

!*****************************************************************************80
!
!! SVSORT descending sorts singular values D and adjusts V.
!
!  Discussion:
!
!    A simple bubble sort is used on D.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    25 February 2002
!
!  Author:
!
!    Original FORTRAN77 version by Richard Brent.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    Richard Brent,
!    Algorithms for Minimization with Derivatives,
!    Prentice Hall, 1973,
!    Reprinted by Dover, 2002.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) N, the length of D, and the order of V.
!
!    Input/output, real ( kind = 8 ) D(N), the vector to be sorted.
!    On output, the entries of D are in descending order.
!
!    Input/output, real ( kind = 8 ) V(N,N), an N by N array to be adjusted
!    as D is sorted.  In particular, if the value that was in D(I) on input is
!    moved to D(J) on output, then the input column V(*,I) is moved to
!    the output column V(*,J).
!
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) d(n)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) j
  integer ( kind = 4 ) j2
  integer ( kind = 4 ) j3
  real ( kind = 8 ) t
  real ( kind = 8 ) v(n,n)

  do j = 1, n - 1

    j3 = j;
    do j2 = j + 1, n
      if ( d(j3) < d(j2) ) then
        j3 = j2
      end if
    end do

    t     = d(j)
    d(j)  = d(j3)
    d(j3) = t

    do i = 1, n
      t       = v(i,j)
      v(i,j)  = v(i,j3)
      v(i,j3) = t
    end do

  end do

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    31 May 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    18 May 2013
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer ( kind = 4 ) d
  integer ( kind = 4 ) h
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) s
  integer ( kind = 4 ) values(8)
  integer ( kind = 4 ) y

  call date_and_time ( values = values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

  if ( h < 12 ) then
    ampm = 'AM'
  else if ( h == 12 ) then
    if ( n == 0 .and. s == 0 ) then
      ampm = 'Noon'
    else
      ampm = 'PM'
    end if
  else
    h = h - 12
    if ( h < 12 ) then
      ampm = 'PM'
    else if ( h == 12 ) then
      if ( n == 0 .and. s == 0 ) then
        ampm = 'Midnight'
      else
        ampm = 'AM'
      end if
    end if
  end if

  write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) &
    d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm )

  return
end

max_stars_repo_name: revitmt/ChemicalKinetics

max_issues_repo_name: revitmt/ChemicalKinetics

max_forks_repo_name: revitmt/ChemicalKinetics


text:
      SUBROUTINE E09407 (USEPRM)

C   E09407 tests the handling of error 407.

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /ERRINF/ ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR,
     1                USRERR, ERRSAV,      FUNSAV,      FILSAV,
     2                EFCNT, EFID
      INTEGER         ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR(10),
     1                USRERR, ERRSAV(200), FUNSAV(200), FILSAV(200),
     2                EFCNT, EFID(100)
      COMMON /ERRCHR/ CURCON,     ERRSRS,    ERRMRK,    ERFLNM,
     1                CONTAB
      CHARACTER       CURCON*200, ERRSRS*40, ERRMRK*20, ERFLNM*80,
     1                CONTAB(40)*150

      INTEGER     USEPRM, ISIZ, IDUM1,ARID,ARNM,NSTID,STLIST(10)
      LOGICAL     STREQ
      CHARACTER   MSG*100

      CURCON = 'some of the specified structures do not exist on '//
     1         'the archive file'
      MSG = 'it should ignore the structures that do not exist.'
      CALL SETVS ('407', EXPERR, EXPSIZ)
      ERRSRS = '8'
      CALL ESETUP (USEPRM)

      ARID = 11
      CALL AVARNM (ARNM)
      CALL POPARF (ARID, ARNM)
      CALL PDASAR (ARID)
      CALL CREST

      CALL PARAST (ARID)
      CALL SETVS ('100,102,101', STLIST, ISIZ)

      CALL RDSTAR (ARID, ISIZ, STLIST)

      CALL PDST (100)
      CALL PDST (101)
      CALL PRAST (ARID)
      CALL PQSID (0, ERRIND, NSTID, IDUM1)
      CALL CHKINQ ('pqsid', ERRIND)
      CALL TSTWAN (STREQ('O**O') .AND. NSTID.EQ.0, MSG)

      CALL CREST

      CALL PARAST (ARID)
      CALL SETVS ('100,102', STLIST, ISIZ)

      CALL RDSNAR (ARID, ISIZ, STLIST)

      CALL PCLARF (ARID)

      CALL ENDERR

      END

ext: f


text:
*dk,ssum
      function ssum(n,x,ix)
C
C ######################################################################
C
C        $Log: ssum.f,v $
C        Revision 2.00  2007/11/03 00:49:13  spchu
C        Import to CVS
C
CPVCS
CPVCS       Rev 1.21   02 Oct 2007 12:40:28   spchu
CPVCS    original version
C
C ######################################################################
C
       implicit real*8 (a-h,o-z)
      dimension x(n)
      ssum=0.0
      do i=1,n,ix
         ssum=ssum+x(i)
      enddo
      goto 9999
 9999 continue
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


module comm

use, intrinsic :: iso_c_binding, only: sp=>C_FLOAT, dp=>C_DOUBLE, int64=>C_LONG_LONG, sizeof=>c_sizeof

implicit none (type, external)
private
public :: rand_init, sp, dp, J, pi, debug

complex(dp),parameter :: J = (0._dp, 1._dp)
real(dp),parameter :: pi = 4._dp*atan(1._dp)
logical :: debug = .false.

contains

subroutine rand_init(repeatable, image_distinct)
logical, intent(in) :: repeatable, image_distinct

call random_init(repeatable, image_distinct)

end subroutine rand_init

end module comm


text:
c
c
c     ###################################################
c     ##  COPYRIGHT (C)  1997  by  Jay William Ponder  ##
c     ##              All Rights Reserved              ##
c     ###################################################
c
c     ###############################################################
c     ##                                                           ##
c     ##  function sigmoid  --  general sigmoidal functional form  ##
c     ##                                                           ##
c     ###############################################################
c
c
c     "sigmoid" implements a normalized sigmoidal function on the
c     interval [0,1]; the curves connect (0,0) to (1,1) and have
c     a cooperativity controlled by beta, they approach a straight
c     line as beta -> 0 and get more nonlinear as beta increases
c
c
      function sigmoid (beta,x)
      implicit none
      real*8 beta,x
      real*8 sigmoid
      real*8 expmax
      real*8 expmin
      real*8 expterm
c
c
c     compute the value of the normalized sigmoidal function
c
      if (beta .eq. 0.0d0) then
         sigmoid = x
      else
         expmax = 1.0d0 / (exp(-beta) + 1.0d0)
         expmin = 1.0d0 / (exp(beta) + 1.0d0)
         expterm = 1.0d0 / (exp(beta*(2.0d0*x-1.0d0)) + 1.0d0)
         sigmoid = (expmax - expterm) / (expmax - expmin)
      end if
      return
      end

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


text:
C    @(#)getxct.f	20.6 8/19/99
      subroutine getxct (jt, kt, mt, nt, vk, qk, ityp, dvdq)

C     Convert BX buses from state PV into state PQ_discrete.

      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha.inc'
      include 'ipfinc/alpha2.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/ecvar.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/slnopt.inc'
      include 'ipfinc/tbx.inc'
      include 'ipfinc/xdata.inc'
      include 'ipfinc/xblim.inc'

      userek = xdata(5,nt)
      usecap = xdata(6,nt)
      vsqr = vk ** 2
      qmin_gen = tbx(4,jt) + qloadu(kt)
      qmax_gen = tbx(3,jt) + qloadu(kt)
      qmax = -qloadu(kt) - amin1 (0.0, userek * vsqr / bmva)
      qmin = -qloadu(kt) - amax1 (0.0, usecap * vsqr / bmva)
      qneut = qmin - vsqr * userek / bmva

      totrek = userek
      totcap = usecap
      isw = 0

      q1 = qneut
      q0 = qneut
      if (qk .lt. qneut .and. userek .lt. 0.0) then
        q3 = tbx(3,jt)
        isw = 1
      else if (qk .gt. qneut .and. usecap .gt. 0.0) then
        q3 = tbx(4,jt)
        isw = 2
      else if (userek .lt. 0.0) then
        q3 = tbx(3,jt)
        isw = 1
      else
        q3 = tbx(4,jt)
        isw = 2
      endif

      do k = 7, 21, 2
        n = xdata(k,nt)
        if (n .eq. 0) go to 160
        if ((isw .eq. 1) .and. (xdata(k+1,nt) .gt. 0.0)) go to 150
        if ((isw .eq. 2) .and. (xdata(k+1,nt) .lt. 0.0)) go to 150
        do l = 1, n
          q2 = q1
          q1 = q1 + xdata(k+1,nt) * vsqr / bmva
          if (isw .eq. 2) then
            if (qk - qmax_gen .lt. q1) go to 170
          else
            if (qk - qmin_gen .gt. q1) go to 170
          endif
        enddo
  150   continue
      enddo

  160 if (isw .eq. 2 .and. abs(qk - qmax_gen - q1) .gt. option(7)) then
        totrek = 0.0
        go to 180
      else if (isw .eq. 1 .and. abs(qk - qmin_gen - q1) .gt. option(7))
     &   then
        totcap = 0.0
        go to 180
      endif

  170 continue
      qnew = bestep (kt, mt, nt, qk, q1, q2, dvdq)
      used = (qnew - q0) * bmva / vsqr

      if (isw .eq. 2) then
        totrek = 0
        totcap = used
      else
        totrek = used
        totcap = 0
      endif

  180 xdata(5,nt) = totrek
      xdata(6,nt) = totcap
      bkku(kt) = bkku(kt) + (totrek + totcap - userek - usecap) / bmva
      qnetu(kt) = dmin1 (tbx(3,jt), dmax1 (tbx(4,jt), dble(q3)))
c
c     Initialize xblim(1,*), xvlim(1,*) as largest B where Vk < Vmin
c     Initialize xblim(2,*), xvlim(2,*) as smallest B where Vk > Vmax
c
      bnew = totrek + totcap
      if (bnew .lt. xdata(4,nt)) then
        if (vk .lt. vlimn(kt)) then
          if (xvlim(1,nt) .eq. 0.0) then
            xblim(1,nt) = bnew
            xvlim(1,nt) = vk
          else if (vk .gt. xvlim(1,nt)) then
            xblim(1,nt) = bnew
            xvlim(1,nt) = vk
          endif
        else if (bnew+0.01 .lt. xblim(1,nt)) then
          if (xvlim(1,nt) .eq. 0.0) then
            xblim(1,nt) = bnew
            xvlim(1,nt) = vk
          else if (xvlim(1,nt) .gt. vlimn(kt) .and.
     &             vk .lt. xvlim(1,nt)) then
            xblim(1,nt) = bnew
            xvlim(1,nt) = vk
          else if (vk .gt. xvlim(1,nt)) then
            xblim(1,nt) = bnew
            xvlim(1,nt) = vk
          endif
        else if (abs(bnew-xblim(1,nt)) .lt. 0.01) then
          xvlim(1,nt) = vk
        endif
      endif
      if (bnew .gt. xdata(3,nt)) then
        if (vk .gt. vlimx(kt)) then
          if (xvlim(2,nt) .eq. 0.0) then
            xblim(2,nt) = bnew
            xvlim(2,nt) = vk
          else if (vk .lt. xvlim(2,nt)) then
            xblim(2,nt) = bnew
            xvlim(2,nt) = vk
          endif
        else if (bnew-0.01 .gt. xblim(2,nt)) then
          if (xvlim(2,nt) .eq. 0.0) then
            xblim(2,nt) = bnew
            xvlim(2,nt) = vk
          else if (xvlim(2,nt) .lt. vlimx(kt) .and.
     &             vk .gt. xvlim(2,nt)) then
            xblim(2,nt) = bnew
            xvlim(2,nt) = vk
          else if (vk .lt. xvlim(2,nt)) then
            xblim(2,nt) = bnew
            xvlim(2,nt) = vk
          endif
        else if (abs(bnew-xblim(2,nt)) .lt. 0.01) then
          xvlim(2,nt) = vk
        endif
      endif
C                                                                      *
      if (idswb .ne. 0) then
        write (dbug,200) jt, kt, userek, usecap,
     &                   totrek, totcap
  200   format (' GETXCT/ Type X ', 2i6, ' X_initial ', 2f10.1,
     &          ' X_final ', 2f10.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


text:
program blackjack

use, intrinsic :: iso_c_binding, only: c_int
use game
implicit none (external)

integer(c_int) :: win, cards(52)
character(2) :: argv
integer :: i

debug = .false.
call get_command_argument(1, argv, status=i)
if(i==0 .and. argv=='-d') debug = .true.

call mix(cards)
print '(52I3)', cards
win = hand(cards)

end program

max_stars_repo_name: fortran-gaming/NIST-Blackjack

max_issues_repo_name: fortran-gaming/NIST-Blackjack

max_forks_repo_name: fortran-gaming/NIST-Blackjack