/RecoHI/HiJetAlgos/plugins/bpmpd-2.11.f

Large files files are truncated, but you can click here to view the full file

C Version for CMSSW, modified to suppress gfortran warnings
C -------------------------------------------------------------------
c Primal-dual method with supernodal cholesky factorization
c               Version 2.11 (1996 December)
c  Written by Cs. Meszaros, MTA SzTAKI, Budapest, Hungary
c        Questions, remarks to the e-mail address:
c               meszaros@lutra.sztaki.hu
c
c  All rights reserved ! Free for academic and research use only !
c  Commercial users are required to purchase a software license.
c
c Related publications:
c
c    Meszaros, Cs.: Fast Cholesky Factorization for Interior Point Methods
c       of Linear Programming. Computers & Mathematics with Applications,
c       Vol. 31. No.4/5 (1996) pp. 49-51.
c
c    Meszaros, Cs.: The "inexact" minimum local fill-in ordering algorithm.
c       Working Paper WP 95-7, Computer and Automation Institute, Hungarian
c       Academy of Sciences
c
c    Maros I., Meszaros Cs.: The Role of the Augmented System in Interior
c       Point Methods.  European Journal of Operations Researches
c       (submitted)
c
c ===========================================================================
c
c  Callable interface
c
c  Standard form: ax-s=b    u>=x,s>=l
c
c  remarks:
c          EQ  rows    0  >= s >=   0
c          GT  rows  +inf >= s >=   0
c          LT  rows    0  >= s >= -inf
c          FR  rows  +inf >= s >= -inf
c
c  input:   obj           objective function (to be minimize)       (n)
c           rhs           right-hand side                           (m)
c           lbound        lower bounds                              (m+n)
c           ubound        upper bounds                              (m+n)
c           colpnt        pointer to the columns                    (n+1)
c           rowidx        row indices                               (nz)
c           nonzeros      nonzero values                            (nz)
c           big           practical +inf
c
c  output: code           termination code
c          xs             primal values
c          dv             dual values
c          dspr           dual resuduals
c
c Input arrays will be destroyed !
c
c ===========================================================================
c
      subroutine solver(
     x obj,rhs,lbound,ubound,diag,odiag,xs,dxs,dxsn,up,dspr,ddspr,
     x ddsprn,dsup,ddsup,ddsupn,dv,ddv,ddvn,prinf,upinf,duinf,scale,
     x nonzeros,
     x vartyp,slktyp,colpnt,ecolpnt,count,vcstat,pivots,invprm,
     x snhead,nodtyp,inta1,prehis,rowidx,rindex,
     x code,opt,iter,corect,fixn,dropn,fnzmax,fnzmin,addobj,
     x bigbou,big,ft)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      common/initv/ prmin,upmax,dumin,stamet,safmet,premet,regul
      real*8        prmin,upmax,dumin
      integer*4     stamet,safmet,premet,regul
c
      integer*4 fixn,dropn,code,iter,corect,fnzmin,fnzmax,ft
      real*8  addobj,opt,big,
     x obj(n),rhs(m),lbound(mn),ubound(mn),scale(mn),diag(mn),odiag(mn),
     x xs(mn),dxs(mn),dxsn(mn),up(mn),dspr(mn),ddspr(mn),ddsprn(mn),
     x dsup(mn),ddsup(mn),ddsupn(mn),dv(m),ddv(m),ddvn(m),
     x nonzeros(cfree),prinf(m),upinf(mn),duinf(mn),bigbou
      integer*4 vartyp(n),slktyp(m),colpnt(n1),ecolpnt(mn),
     x count(mn),vcstat(mn),pivots(mn),invprm(mn),snhead(mn),
     x nodtyp(mn),inta1(mn),prehis(mn),rowidx(cfree),rindex(rfree)
c
      common/numer/ tplus,tzer
      real*8        tplus,tzer
      common/ascal/ objnor,rhsnor,scdiff,scpass,scalmet
      real*8        objnor,rhsnor,scdiff
      integer*4     scpass,scalmet
c ---------------------------------------------------------------------------
      integer*4 i,j,k,active,pnt1,pnt2,prelen,freen
      real*8 scobj,scrhs,sol,lbig
      character*99 buff
C CMSSW: Temporary integer array needed to avoid reusing REAL*8 for
C integer storage
      integer*4 pmbig(m),ppbig(m),dmbig(n),dpbig(n)
      integer*4 iwork1(mn+mn),iwork2(mn+mn),iwork3(mn+mn),iwork4(mn+mn),
     &     iwork5(mn+mn)
c ---------------------------------------------------------------------------
c
c inicializalas
c
      if(cfree.le.(nz+1)*2)then
        write(buff,'(1x,a)')'Not enough memory, realmem < nz !'
        call mprnt(buff)
        code=-2
        goto 50
      endif
      if(rfree.le.nz)then
        write(buff,'(1x,a)')'Not enough memory, intmem < nz !'
        call mprnt(buff)
        code=-2
        goto 50
      endif
      iter=0
      corect=0
      prelen=0
      fnzmin=cfree
      fnzmax=-1
      scobj=1.0d+0
      scrhs=1.0d+0
      code=0
      lbig=0.9d+0*big
      if(bigbou.gt.lbig)then
        lbig=bigbou
        big=lbig/0.9d+0
      endif
      do i=1,mn
        scale(i)=1.0d+0
      enddo
c
c Remove fix variables and free rows
c
      do i=1,n
        vartyp(i)=0
        if(abs(ubound(i)-lbound(i)).le.tplus*(abs(lbound(i)+1.0d0)))then
          vartyp(i)= 1
          vcstat(i)=-2-1
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            rhs(rowidx(j))=rhs(rowidx(j))-ubound(i)*nonzeros(j)
          enddo
          addobj=addobj+obj(i)*lbound(i)
        else
          vcstat(i)=0
        endif
      enddo
      do i=1,m
        slktyp(i)=0
        j=i+n
        if((ubound(j).gt.lbig).and.(lbound(j).lt.-lbig))then
          vcstat(j)=-2-1
        else
          vcstat(j)=0
        endif
      enddo
c
c   p r e s o l v e r
c
      call timer(k)
      if(premet.gt.0)then
        write(buff,'(1x)')
        call mprnt(buff)
        write(buff,'(1x,a)')'Process: presolv'
        call mprnt(buff)
        call presol(colpnt,rowidx,nonzeros,rindex,nonzeros(nz+1),
     x  snhead,snhead(n1),nodtyp,nodtyp(n1),vcstat,vcstat(n1),
     x  ecolpnt,count,ecolpnt(n1),count(n1),
C CMSSW: Prevent REAL*8 reusage warning
C Was:  vartyp,dxsn(n1),dxs(n1),diag(n1),odiag(n1),
     x  vartyp,dxsn(n1),dxs(n1),pmbig,ppbig,
     x  ubound,lbound,ubound(n1),lbound(n1),rhs,obj,prehis,prelen,
C CMSSW: Prevent REAL*8 reusage warning
C Was:  addobj,big,pivots,invprm,dv,ddv,dxsn,dxs,diag,odiag,premet,code)
     x  addobj,big,pivots,invprm,dv,ddv,dxsn,dxs,dmbig,dpbig,premet,
     x  code)
        write(buff,'(1x,a)')'Presolv done...'
        call mprnt(buff)
        if(code.ne.0)goto 45
      endif
c
c Remove lower bounds
c
      call stndrd(ubound,lbound,rhs,obj,nonzeros,
     x vartyp,slktyp,vcstat,colpnt,rowidx,addobj,tplus,tzer,lbig,big)
c
c Scaling before aggregator
c
      i=iand(scalmet,255)
      j=iand(scpass,255)
      if(i.gt.0)call mscale(colpnt,rowidx,nonzeros,obj,rhs,ubound,
     x vcstat,scale,upinf,i,j,scdiff,ddsup,dxsn,dxs,snhead)
c
c Aggregator
c
      if(premet.gt.127)then
        write(buff,'(1x)')
        call mprnt(buff)
        write(buff,'(1x,a)')'Process: aggregator'
        call mprnt(buff)
        call aggreg(colpnt,rowidx,nonzeros,rindex,
     x  vcstat,vcstat(n1),ecolpnt,count,ecolpnt(n1),count(n1),
     x  rhs,obj,prehis,prelen,pivots,vartyp,slktyp,invprm,snhead,
     x  nodtyp,inta1,inta1(n1),dv,addobj,premet,code)
        write(buff,'(1x,a)')'Aggregator done...'
        call mprnt(buff)
        if(code.ne.0)goto 55
      endif
c
c Scaling after aggregator
c
      i=scalmet/256
      j=scpass/256
      if(i.gt.0)call mscale(colpnt,rowidx,nonzeros,obj,rhs,
     x ubound,vcstat,scale,upinf,i,j,scdiff,ddsup,dxsn,dxs,snhead)
c
      call timer(j)
      write(buff,'(1x)')
      call mprnt(buff)
      write(buff,'(1x,a,f8.2,a)')
     x 'Time for presolv, scaling and aggregator: ',0.01*(j-k),' sec.'
      call mprnt(buff)
c
c cleaning
c
      do i=1,mn
        xs(i)=0.0d+0
        dspr(i)=0.0d+0
        dsup(i)=0.0d+0
        up(i)=0.0d+0
      enddo
      do i=1,m
        dv(i)=0.0d+0
      enddo
c
c Is the problem solved ?
c
      fixn=0
      dropn=0
      freen=0
      do i=1,n
        if(vcstat(i).le.-2)then
          fixn=fixn+1
        else if(vartyp(i).eq.0) then
          freen=freen+1
        endif
      enddo
      do i=1,m
        if(vcstat(i+n).le.-2)dropn=dropn+1
      enddo
      active=mn-fixn-dropn
      if(active.eq.0)code=2
      if(code.gt.0)then
        opt=addobj
        write(buff,'(1x,a)')'Problem is solved by the pre-solver'
        call mprnt(buff)
        if(code.gt.0)goto 55
        goto 50
      endif
c
c Presolve statistics
c
      if(premet.gt.0)then
        i=0
        j=0
        do k=1,n
          if(vcstat(k).gt.-2)then
            i=i+count(k)-ecolpnt(k)+1
            if(j.lt.count(k)-ecolpnt(k)+1)j=count(k)-ecolpnt(k)+1
          endif
        enddo
        write(buff,'(1x,a22,i8)')'Number of rows       :',(m-dropn)
        call mprnt(buff)
        write(buff,'(1x,a22,i8)')'Number of columns    :',(n-fixn)
        call mprnt(buff)
        write(buff,'(1x,a22,i8)')'Free variables       :',freen
        call mprnt(buff)
        write(buff,'(1x,a22,i8)')'No. of nonzeros      :',i
        call mprnt(buff)
        write(buff,'(1x,a22,i8)')'Longest column count :',j
        call mprnt(buff)
      endif
c
c Incrase rowidx by n
c
      j=colpnt(1)
      k=colpnt(n+1)-1
      do i=j,k
        rowidx(i)=rowidx(i)+n
      enddo
      active=mn-fixn-dropn
c
c Normalize obj and rhs
c
      if(objnor.gt.tzer)then
        call scalobj(obj,scobj,vcstat,objnor)
      endif
      if(rhsnor.gt.tzer)then
        call scalrhs(rhs,scrhs,vcstat,rhsnor,ubound,xs,up)
      endif
c
c Calling phas12
c
      sol=scobj*scrhs
      i=mn+mn
      call timer(k)
      call phas12(
     x obj,rhs,ubound,diag,odiag,xs,dxs,dxsn,up,dspr,ddspr,
     x ddsprn,dsup,ddsup,ddsupn,dv,ddv,ddvn,nonzeros,prinf,upinf,duinf,
     x vartyp,slktyp,colpnt,ecolpnt,count,vcstat,pivots,invprm,
     x snhead,nodtyp,inta1,rowidx,rindex,
C CMSSW: Prevent REAL*8 reusage warning
C Was: dxs,dxsn,ddspr,ddsprn,ddsup,ddsupn,
     x dxs,iwork1,iwork2,iwork3,iwork4,iwork5,
     x code,opt,iter,corect,fixn,dropn,active,fnzmax,fnzmin,addobj,
     x sol,ft,i)
      call timer(j)
      write(buff,'(1x,a,f11.2,a)')'Solver time ',0.01*(j-k),' sec.'
      call mprnt(buff)
c
c Decrease rowidx by n
c
      j=colpnt(1)
      k=colpnt(n+1)-1
      do i=j,k
        rowidx(i)=rowidx(i)-n
      enddo
c
c Rescaling
c
  55  do i=1,m
        rhs(i)=rhs(i)*scrhs*scale(i+n)
        ubound(i+n)=ubound(i+n)*scrhs*scale(i+n)
        xs(i+n)=xs(i+n)*scrhs*scale(i+n)
        up(i+n)=up(i+n)*scrhs*scale(i+n)
        dv(i)=dv(i)*scobj/scale(i+n)
        dspr(i+n)=dspr(i+n)/scale(i+n)*scobj
        dsup(i+n)=dsup(i+n)/scale(i+n)*scobj
      enddo
c
      do i=1,n
        obj(i)=obj(i)*scobj*scale(i)
        ubound(i)=ubound(i)*scrhs/scale(i)
        pnt1=colpnt(i)
        pnt2=colpnt(i+1)-1
        do j=pnt1,pnt2
          nonzeros(j)=nonzeros(j)*scale(i)*scale(rowidx(j)+n)
        enddo
c
        xs(i)=xs(i)/scale(i)*scrhs
        up(i)=up(i)/scale(i)*scrhs
        dspr(i)=dspr(i)*scale(i)*scobj
        dsup(i)=dsup(i)*scale(i)*scobj
      enddo
c
c Postprocessing
c
  45  call pstsol(colpnt,rowidx,nonzeros,vcstat,vcstat(n1),
     x vartyp,slktyp,ubound,lbound,ubound(n1),lbound(n1),rhs,obj,xs,
     x inta1,ddvn,prehis,prelen,big)
c
  50  return
      end
c
c ===========================================================================
c
      subroutine stndrd(ubound,lbound,rhs,obj,nonzeros,
     x vartyp,slktyp,vcstat,colpnt,rowidx,addobj,tplus,tzer,lbig,big)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      integer*4 vartyp(n),slktyp(m),vcstat(mn),colpnt(n1),rowidx(nz)
      real*8 ubound(mn),lbound(mn),rhs(m),obj(n),nonzeros(nz),
     x addobj,tplus,tzer,lbig,big
c
      integer*4 i,j,k,pnt1,pnt2
c
c generate standard form, row modification
c
      k=0
      do 150 i=1,m
        j=i+n
        if(vcstat(j).gt.-2)then
          if(abs(ubound(j)-lbound(j)).le.tplus*(abs(lbound(j))+1d0))then
            slktyp(i)=0
            ubound(j)=0.0d+00
            rhs(i)=rhs(i)+lbound(j)
            goto 150
          endif
ccc          if((ubound(j).gt.lbig).and.(lbound(j).lt.-lbig))then
ccc            vcstat(j)=-2
ccc            slktyp(i)=0
ccc            goto 150
ccc          endif
          if(lbound(j).lt.-lbig)then
            slktyp(i)=2
            lbound(j)=-ubound(j)
            ubound(j)=big
            rhs(i)=-rhs(i)
            k=k+1
          else
            slktyp(i)=1
          endif
          rhs(i)=rhs(i)+lbound(j)
          ubound(j)=ubound(j)-lbound(j)
          if(ubound(j).lt.lbig)slktyp(i)=-slktyp(i)
        else
          slktyp(i)=0
        endif
 150  continue
c
c negate reverse rows
c
      if(k.gt.0)then
        do i=1,n
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            if(abs(slktyp(rowidx(j))).ge.2)nonzeros(j)=-nonzeros(j)
          enddo
        enddo
      endif
c
c column modification
c
      do 155 i=1,n
        if(vcstat(i).gt.-2)then
ccc          if(abs(ubound(i)-lbound(i)).le.tplus*(abs(lbound(i))+1d0))then
ccc            vcstat(i)=-2
ccc            vartyp(i)= 1
ccc            do j=colpnt(i),colpnt(i+1)-1
ccc              rhs(rowidx(j))=rhs(rowidx(j))-nonzeros(j)*lbound(i)
ccc            enddo
ccc            addobj=addobj+obj(i)*lbound(i)
ccc            goto 155
ccc          endif
          if((ubound(i).gt.lbig).and.(lbound(i).lt.-lbig))then
            vartyp(i)=0
            goto 155
          endif
          if(lbound(i).lt.-lbig)then
            vartyp(i)=2
            lbound(i)=-ubound(i)
            ubound(i)=big
            obj(i)=-obj(i)
            do j=colpnt(i),colpnt(i+1)-1
              nonzeros(j)=-nonzeros(j)
            enddo
          else
            vartyp(i)=1
          endif
          if(abs(lbound(i)).gt.tzer)then
            if(ubound(i).lt.lbig)ubound(i)=ubound(i)-lbound(i)
            do j=colpnt(i),colpnt(i+1)-1
              rhs(rowidx(j))=rhs(rowidx(j))-nonzeros(j)*lbound(i)
            enddo
            addobj=addobj+obj(i)*lbound(i)
          endif
          if(ubound(i).lt.lbig)vartyp(i)=-vartyp(i)
        endif
 155  continue
      return
      end
c
c ===========================================================================
c Primal-dual method with supernodal cholesky factorization
c               Version 2.11 (1996 December)
c  Written by Cs. Meszaros, MTA SzTAKI, Budapest, Hungary
c           e-mail: meszaros@lutra.sztaki.hu
c                      see "bpmain.f"
c
c code=-2 General memory limit (no solution)
c code=-1 Memory limit during iterations
c code= 0
c code= 1 No optimum
c code= 2 Otimal solution
c code= 3 Primal Infeasible
c code= 4 Dual Infeasible
c
c ===========================================================================
c
      subroutine phas12(
     x obj,rhs,bounds,diag,odiag,xs,dxs,dxsn,up,dspr,ddspr,
     x ddsprn,dsup,ddsup,ddsupn,dv,ddv,ddvn,nonzeros,prinf,upinf,duinf,
     x vartyp,slktyp,colpnt,ecolpnt,count,vcstat,pivots,invprm,
     x snhead,nodtyp,inta1,rowidx,rindex,
     x rwork1,iwork1,iwork2,iwork3,iwork4,iwork5,
     x code,opt,iter,corect,fixn,dropn,active,fnzmax,fnzmin,addobj,
     x scobj,factim,mn2)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      common/mscal/ varadd,slkadd,scfree
      real*8        varadd,slkadd,scfree
c
      common/numer/ tplus,tzer
      real*8        tplus,tzer
c
      common/param/ palpha,dalpha
      real*8        palpha,dalpha
c
      common/factor/ tpiv1,tpiv2,tabs,trabs,lam,tfind,order,supdens
      real*8         tpiv1,tpiv2,tabs,trabs,lam,tfind,order,supdens
c
      common/toler/ tsdir,topt1,topt2,tfeas1,tfeas2,feas1,feas2,
     x              pinfs,dinfs,inftol,maxiter
      real*8        tsdir,topt1,topt2,tfeas1,tfeas2,feas1,feas2,
     x              pinfs,dinfs,inftol
      integer*4     maxiter
c
      common/initv/ prmin,upmax,dumin,stamet,safmet,premet,regul
      real*8        prmin,upmax,dumin
      integer*4     stamet,safmet,premet,regul
c
      integer*4 fixn,dropn,active,code,iter,corect,fnzmin,fnzmax,mn2
      real*8  addobj,scobj,opt
c
      common/predp/ ccstop,barset,bargrw,barmin,mincor,maxcor,inibar
      real*8        ccstop,barset,bargrw,barmin
      integer*4     mincor,maxcor,inibar
c
      common/predc/ target,tsmall,tlarge,center,corstp,mincc,maxcc
      real*8        target,tsmall,tlarge,center,corstp
      integer*4     mincc,maxcc

      common/itref/ tresx,tresy,maxref
      real*8        tresx,tresy
      integer*4     maxref
c
      real*8 obj(n),rhs(m),bounds(mn),diag(mn),odiag(mn),xs(mn),
     x dxs(mn),dxsn(mn),up(mn),dspr(mn),ddspr(mn),ddsprn(mn),dsup(mn),
     x ddsup(mn),ddsupn(mn),dv(m),ddv(m),ddvn(m),nonzeros(cfree),
     x prinf(m),upinf(mn),duinf(mn),rwork1(mn)

      integer*4 vartyp(n),slktyp(m),colpnt(n1),ecolpnt(mn),count(mn),
     x vcstat(mn),pivots(mn),invprm(mn),snhead(mn),nodtyp(mn),
     x inta1(mn),rowidx(cfree),rindex(rfree),factim,
     x iwork1(mn2),iwork2(mn2),iwork3(mn2),iwork4(mn2),iwork5(mn2)
c
c ---------------------------------------------------------------------------
c
      integer*4 i,j,err,factyp,pphase,dphase,t1,t2,opphas,odphas
      real*8 pinf,dinf,uinf,prelinf,drelinf,popt,dopt,cgap,
     x prstpl,dustpl,barpar,oper,maxstp,pinfrd,dinfrd,objerr,nonopt,
     x oprelinf,odrelinf,opinf,odinf,ocgap
      integer*4 corr,corrc,barn,fxp,fxd,fxu,nropt
      character*99 buff,sbuff
      character*1 wmark
c
c to save parameters
c
      integer*4 maxcco,mxrefo
      real*8 lamo,spdeno,bargro,topto
C CMSSW: Temporary integer array needed to avoid reusing REAL*8 for
C integer storage
      integer*4 inta12(mn)
c
c --------------------------------------------------------------------------
c
 101  format(1x,' ')
 102  format(1x,'It-PC   P.Inf   D.Inf  U.Inf   Actions           ',
     x 'P.Obj           D.Obj  Barpar')
 103  format(1x,'------------------------------------------------',
     x '------------------------------')
 104  format(1x,I2,a1,I1,I1,' ',1PD7.1,' ',1PD7.1,' ',1PD6.0,
     x ' ',I2,' ',I3,' ',I3,' ',1PD15.8,' ',1PD15.8,' ',1PD6.0)
c
c Saving parameters
c
      maxcco=maxcc
      mxrefo=maxref
      lamo=lam
      spdeno=supdens
      bargro=bargrw
      topto=topt1
c
c Include dummy ranges if requested
c
      if(regul.gt.0)then
        do i=1,m
          if(slktyp(i).eq.0)then
            slktyp(i)=-1
            bounds(i+n)=0.0d+0
          endif
        enddo
      endif
c
c Other initialization
c
      nropt=0
      factim=0
      wmark='-'
      fxp=0
      fxd=0
      fxu=0
c
      call stlamb(colpnt,vcstat,rowidx,inta1,fixn,dropn,factyp)
      call timer(t1)
      j=0
      do i=1,n
        if((vcstat(i).gt.-2).and.(vartyp(i).eq.0))j=j+1
      enddo
      if((j.gt.0).and.(scfree.lt.tzer))factyp=1
c
c Initial scaling matrix (diagonal)
c
      call fscale (vcstat,diag,odiag,vartyp,slktyp)
      do i=1,m
        dv(i)=0.0d+0
      enddo

ccc      i=2*rfree
ccc      j=400
ccc      call paintmat(m,n,nz,i,rowidx,colpnt,rindex,j,'matrix01.pic')

c
c Initial factorization
c
      fnzmax=0
      if(factyp.eq.1)then
        call ffactor(ecolpnt,vcstat,colpnt,rowidx,
     x  iwork4,pivots,count,nonzeros,diag,
     x  iwork1,iwork1(mn+1),iwork2,iwork2(mn+1),inta1,iwork5,
     x  iwork5(mn+1),iwork3,iwork3(mn+1),iwork4(mn+1),rindex,
     x  rwork1,fixn,dropn,fnzmax,fnzmin,active,oper,xs,slktyp,code)
        if(code.ne.0)goto 999
        call supnode(ecolpnt,count,rowidx,vcstat,pivots,snhead,
     x  invprm,nodtyp)
      else
c
c minimum local fill-in ordering
c
        i=int(tfind)
        if(order.lt.1.5)i=0
        if(order.lt.0.5)i=-1
        call symmfo(inta1,pivots,ecolpnt,vcstat,
     x  colpnt,rowidx,nodtyp,rindex,iwork3,invprm,
     x  count,snhead,iwork1,iwork1(mn+1),iwork2,iwork2(mn+1),
     x  iwork4,iwork4(mn+1),iwork3(mn+1),iwork5,iwork5(mn+1),
C CMSSW: Prevent REAL*8 reusage warning
C Was:  nonzeros,fnzmax,oper,i,rwork1,code
     x  nonzeros,fnzmax,oper,i,inta12,code)
        if(code.ne.0)goto 999
        call supnode(ecolpnt,count,rowidx,vcstat,pivots,snhead,
     x  invprm,nodtyp)
        popt=trabs
        trabs=tabs
        call nfactor(ecolpnt,vcstat,rowidx,pivots,count,nonzeros,
     x  diag,err,rwork1,iwork2,iwork2(mn+1),dropn,slktyp,
     x  snhead,iwork3,invprm,nodtyp,dv,odiag)
        trabs=popt
      endif
      fnzmin=fnzmax
c
c Compute centrality and iterative refinement power
c
      if(fnzmin.eq.0)fnzmin=1
      cgap=oper/fnzmin/10.0d+0
      j=0
  78  if(cgap.ge.1.0d+0)then
        cgap=cgap/2
        j=j+1
        goto 78
      endif
      if(j.eq.0)j=1
      if(maxcc.le.0d+0)then
        maxcc=-maxcc
      else
        if(j.le.maxcc)maxcc=j
      endif
      if(mincc.gt.maxcc)maxcc=mincc
      cgap=log(1.0d+0+oper/fnzmin/5.0d+0)/log(2.0d+00)
      if(maxref.le.0)then
        maxref=-maxref
      else
        maxref=int(cgap*maxref)
      endif
      if(maxref.le.0)maxref=0
      write(buff,'(1x,a,i2)')'Centrality correction Power:',maxcc
      call mprnt(buff)
      write(buff,'(1x,a,i2)')'Iterative refinement  Power:',maxref
      call mprnt(buff)
c
c Starting point
c
      call initsol(xs,up,dv,dspr,dsup,rhs,obj,bounds,vartyp,slktyp,
     x vcstat,colpnt,ecolpnt,pivots,rowidx,nonzeros,diag,rwork1,
     x count)
      call timer(t2)
c
      write(buff,'(1x,a,f12.2,a)')'FIRSTFACTOR TIME :',
     x (dble(t2-t1)*0.01d+0),' sec'
      call mprnt(buff)
c
      maxstp=1.0d+0
      iter=0
      corect=0
      corr=0
      corrc=0
      barn=0
      cgap=0.0d+0
      do i=1,mn
        if(vcstat(i).gt.-2)then
          if(i.le.n)then
            j=vartyp(i)
          else
            j=slktyp(i-n)
          endif
          if(j.ne.0)then
            cgap=cgap+xs(i)*dspr(i)
            barn=barn+1
          endif
          if(j.lt.0)then
            cgap=cgap+up(i)*dsup(i)
            barn=barn+1
          endif
        endif
      enddo
      if(barn.lt.1)barn=1

ccc      i=2*rfree
ccc      j=350
ccc      call paintaat(mn,nz,pivotn,i,rowidx,ecolpnt,count,rindex,
ccc     x j,pivots,iwork1,iwork1(mn+1),iwork2,iwork2(mn+1),iwork3,
ccc     x iwork3(mn+1),'normal01.pic')

ccc      i=2*rfree
ccc      j=400
ccc      call paintata(mn,nz,pivotn,i,rowidx,ecolpnt,count,rindex,
ccc     x j,pivots,iwork1,iwork1(mn+1),iwork2,iwork2(mn+1),iwork3,
ccc     x 'atapat01.pic')


ccc      i=2*rfree
ccc      j=350
ccc      err=nz
ccc      call paintfct(mn,cfree,pivotn,i,rowidx,ecolpnt,count,rindex,
ccc     x j,pivots,iwork2,err,'factor01.pic')
c
c Initialize for the iteration loop
c
      do i=1,n
        if((vcstat(i).gt.-2).and.(vartyp(i).ne.0))then
          if(xs(i).gt.dspr(i))then
            vcstat(i)=1
          else
            vcstat(i)=0
          endif
        endif
      enddo
      do i=1,m
        if((vcstat(i+n).gt.-2).and.(slktyp(i).ne.0))then
          if(xs(i+n).gt.dspr(i+n))then
            vcstat(i+n)=1
          else
            vcstat(i+n)=0
          endif
        endif
      enddo
      opphas=0
      odphas=0
      pinfrd=1.0d+0
      dinfrd=1.0d+0
      barpar=0.0d+0
c
c main iteration loop
c
  10  if(mod(iter,20).eq.0)then
        write(buff,101)
        call mprnt(buff)
        write(buff,102)
        call mprnt(buff)
        write(buff,103)
        call mprnt(buff)
      endif
c
c Infeasibilities
c
      call cprinf(xs,prinf,slktyp,colpnt,rowidx,nonzeros,
     x rhs,vcstat,pinf)
      call cduinf(dv,dspr,dsup,duinf,vartyp,slktyp,colpnt,rowidx,
     x nonzeros,obj,vcstat,dinf)
      call cupinf(xs,up,upinf,bounds,vartyp,slktyp,vcstat,
     x uinf)
c
c Objectives
c
      call cpdobj(popt,dopt,obj,rhs,bounds,xs,dv,dsup,
     x vcstat,vartyp,slktyp)
      popt=scobj*popt+addobj
      dopt=scobj*dopt+addobj
c
c Stopping criteria
c
      call stpcrt(prelinf,drelinf,popt,dopt,cgap,iter,
     x code,pphase,dphase,maxstp,pinf,uinf,dinf,
     x prinf,upinf,duinf,nonopt,pinfrd,dinfrd,
     x prstpl,dustpl,obj,rhs,bounds,xs,dxs,dspr,ddspr,dsup,ddsup,dv,ddv,
     x up,addobj,scobj,vcstat,vartyp,slktyp,
     x oprelinf,odrelinf,opinf,odinf,ocgap,opphas,odphas,sbuff)
c
      write(buff,104)iter,wmark,corr,corrc,pinf,dinf,uinf,fxp,fxd,fxu,
     x popt,dopt,barpar
      call mprnt(buff)
      if(code.ne.0)then
        write(buff,'(1x)')
        call mprnt(buff)
        call mprnt(sbuff)
        goto 90
      endif
c
c P-D solution modification
c
      call pdmodi(xs,dspr,vcstat,vartyp,slktyp,cgap,popt,
     x dopt,prinf,duinf,upinf,colpnt,rowidx,nonzeros,pinf,uinf,dinf)
c
c Fixing variables / dropping rows / handling dual slacks
c
      i=fixn
      call varfix(vartyp,slktyp,rhs,colpnt,rowidx,nonzeros,
     x xs,up,dspr,dsup,vcstat,fixn,dropn,addobj,scobj,obj,bounds,
     x duinf,dinf,fxp,fxd,fxu)
      if(fixn.ne.i)then
        call supupd(pivots,invprm,snhead,nodtyp,vcstat,ecolpnt)
        call cprinf(xs,prinf,slktyp,colpnt,rowidx,nonzeros,
     x  rhs,vcstat,pinf)
        call cupinf(xs,up,upinf,bounds,vartyp,slktyp,vcstat,
     x  uinf)
      endif
c
c Compute gap
c
      cgap=0.0d+0
      do i=1,mn
        if(vcstat(i).gt.-2)then
          if(i.le.n)then
            j=vartyp(i)
          else
            j=slktyp(i-n)
          endif
          if(j.ne.0)then
            cgap=cgap+xs(i)*dspr(i)
            if(j.lt.0)then
              cgap=cgap+up(i)*dsup(i)
            endif
          endif
        endif
      enddo
c
c Computation of the scaling matrix
c
      objerr=abs(dopt-popt)/(abs(popt)+1.0d+0)
      call cdiag(xs,up,dspr,dsup,vartyp,slktyp,vcstat,diag,odiag)
      pinfrd=pinf
      dinfrd=dinf
c
c The actual factorization
c
  50  err=0
      call timer(t1)
      if (factyp.eq.1) then
        call mfactor(ecolpnt,vcstat,colpnt,rowidx,pivots,
     x  count,iwork4,nonzeros,diag,err,rwork1,iwork2,iwork2(mn+1),
     x  dropn,slktyp,snhead,iwork3,invprm,nodtyp,dv,odiag)
      else
        call nfactor(ecolpnt,vcstat,rowidx,pivots,count,nonzeros,
     x  diag,err,rwork1,iwork2,iwork2(mn+1),dropn,slktyp,
     x  snhead,iwork3,invprm,nodtyp,dv,odiag)
      endif
      call timer(t2)
      if(err.gt.0)then
        do i=1,mn
          diag(i)=odiag(i)
        enddo
        call newsmf(colpnt,pivots,rowidx,nonzeros,ecolpnt,count,
     x  vcstat,invprm,snhead,nodtyp,iwork1,rwork1,iwork2,iwork3,
     x  iwork4,code)
        if(code.lt.0)then
          write(buff,'(1x)')
          call mprnt(buff)
          goto 90
        endif
        goto 50
      endif
      factim=factim+t2-t1
c
c We are in the finish ?
c
      wmark(1:1)='-'
      if(objerr.gt.1.0d+0)objerr=1.0d+0
      if(objerr.lt.topt1)objerr=topt1
      if((objerr.le.topt1*10.0d+0).and.(pphase+dphase.eq.4))then
         if(bargrw.gt.0.1d+0)bargrw=0.1d+0
         nropt=nropt+1
         if(nropt.eq.5)then
           nropt=0
           topt1=topt1*sqrt(10.d+0)
           write(buff,'(1x,a)')'Near otptimal but slow convergence.'
           call mprnt(buff)
         endif
         wmark(1:1)='+'
      endif
c
c primal-dual predictor-corrector direction
c
      call  cpdpcd(xs,up,dspr,dsup,prinf,duinf,upinf,
     x dxsn,ddvn,ddsprn,ddsupn,dxs,ddv,ddspr,ddsup,bounds,
     x ecolpnt,count,pivots,vcstat,diag,odiag,rowidx,nonzeros,
     x colpnt,vartyp,slktyp,barpar,corr,prstpl,dustpl,barn,cgap)
      corect=corect+corr
c
c primal-dual centality-correction
c
      call  cpdccd(xs,up,dspr,dsup,upinf,
     x dxsn,ddvn,ddsprn,ddsupn,dxs,ddv,ddspr,ddsup,bounds,
     x ecolpnt,count,pivots,vcstat,diag,odiag,rowidx,nonzeros,
     x colpnt,vartyp,slktyp,barpar,corrc,prstpl,dustpl)
      corect=corect+corrc
c
c compute steplengths
c
      iter=iter+1
      prstpl=prstpl*palpha
      dustpl=dustpl*dalpha
c
c compute the new primal-dual solution
c
      call cnewpd(prstpl,xs,dxs,up,upinf,dustpl,dv,ddv,dspr,
     x ddspr,dsup,ddsup,vartyp,slktyp,vcstat,maxstp)
c
c End main loop
c
      goto 10
c
 90   opt=(dopt-popt)/(abs(popt)+1.0d+0)
      write(buff,'(1x,a,1PD11.4,a,1PD18.10)')
     x 'ABSOLUTE infeas.   Primal  :',pinf,   '    Dual         :',dinf
      call mprnt(buff)
      write(buff,'(1x,a,1PD11.4,a,1PD18.10)')
     x 'PRIMAL :  Relative infeas. :',prelinf,'    Objective    :',popt
      call mprnt(buff)
      write(buff,'(1x,a,1PD11.4,a,1PD18.10)')
     x 'DUAL   :  Relative infeas. :',drelinf,'    Objective    :',dopt
      call mprnt(buff)
      write(buff,'(1x,a,1PD11.4,a,1PD18.10)')
     x 'Complementarity gap        :',cgap,'    Duality gap  :',opt
      call mprnt(buff)
      opt=popt
c
c Restoring parameters
c
 999  maxcc=maxcco
      maxref=mxrefo
      lam=lamo
      supdens=spdeno
      bargrw=bargro
      topt1=topto
      return
      end
c
c ===========================================================================
c ===========================================================================
c
      subroutine mscale(colpnt,rowidx,nonzeros,
     x obj,rhs,ubound,vcstat,scale,scalen,scalmet,scpass,scdiff,
     x ddsup,ddsupn,dxs,snhead)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      integer*4 colpnt(n1),rowidx(nz),vcstat(mn),
     x scalmet,scpass,snhead(mn)
      real*8    nonzeros(cfree),obj(n),rhs(m),ubound(mn),scale(mn),
     x scalen(mn),scdiff,ddsup(mn),ddsupn(mn),dxs(mn)
c
      integer*4 i
      character*99 buff
c
      write(buff,'(1x)')
      call mprnt(buff)
      write(buff,'(1x,a)')'Process: scaling'
      call mprnt(buff)
c
      do i=1,mn
        scalen(i)=1.0d+0
      enddo
c
      if((scalmet.eq.2).or.(scalmet.eq.4))then
        call scale1(ubound,nonzeros,colpnt,obj,scalen,vcstat,
     x  rowidx,rhs,ddsup,scpass,scdiff,snhead,nonzeros(nz+1))
      endif
      if((scalmet.eq.3).or.(scalmet.eq.5))then
        call scale2(ubound,nonzeros,colpnt,obj,scalen,vcstat,
     x  rowidx,rhs,scpass,scdiff,ddsup,ddsupn,dxs,snhead)
      endif
      if((scalmet.gt.0).and.(scalmet.le.3))then
        call sccol2(ubound,nonzeros,colpnt,obj,scalen,
     x  vcstat,rowidx)
        call scrow2(rhs,ubound,nonzeros,rowidx,colpnt,ddsup,
     x  scalen,vcstat)
      endif
c
      do i=1,mn
        scale(i)=scale(i)*scalen(i)
      enddo
c
      write(buff,'(1x,a)')'Scaling done...'
      call mprnt(buff)
      return
      end
c
c ============================================================================
c
      subroutine scale1(bounds,rownzs,colpnt,obj,scale,
     x vcstat,rowidx,rhs,work1,scpass,scdif,veclen,
     x lognz)

      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      common/numer/ tplus,tzer
      real*8        tplus,tzer

      real*8 bounds(mn),rownzs(cfree),obj(n),scale(mn),
     x  rhs(m),work1(mn),scdif,lognz(nz)
      integer*4 rowidx(cfree),colpnt(n1),vcstat(mn),scpass,veclen(mn)
c
      real*8 defic,odefic
      integer*4 pass,i,j,pnt1,pnt2,nonz
      character*99 buff
c
      pass=0
      nonz=0
      defic= 1.0d+0
      odefic=0.0d+0
      do i=1,mn
        veclen(i)=0
      enddo
      do i=1,n
        if(vcstat(i).gt.-2)then
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            if((abs(rownzs(j)).gt.tzer).and.
     x      (vcstat(rowidx(j)+n).gt.-2))then
              lognz(j)=log(abs(rownzs(j)))
              veclen(i)=veclen(i)+1
              veclen(rowidx(j)+n)=veclen(rowidx(j)+n)+1
              nonz=nonz+1
              odefic=odefic+abs(lognz(j))
            else
              lognz(j)=0.0d+0
            endif
          enddo
        endif
      enddo
      do i=1,mn
        if(veclen(i).eq.0)veclen(i)=1
        scale(i)=0.0d+0
      enddo
      if(nonz.eq.0)goto 999
      odefic=exp(odefic/dble(nonz))
      if(odefic.le.scdif)goto 999
  10  write(buff,'(1x,a,i2,a,d12.6)')'Pass',pass,'. Average def.',odefic
      call mprnt(buff)
      call sccol1(colpnt,scale,
     x vcstat,rowidx,veclen,lognz)
      pass=pass+1
      call scrow1(rowidx,colpnt,work1,scale,vcstat,defic,veclen,lognz)
      defic=exp(defic/dble(nonz))
      if(defic.le.scdif)goto 999
      if(pass.ge.scpass)goto 999
      if(odefic.le.defic)goto 999
      odefic=defic
      goto 10
 999  write(buff,'(1x,a,i2,a,d12.6)')'Pass',pass,'. Average def.',defic
      call mprnt(buff)
c
c Scaling
c
      do i=1,mn
        scale(i)=exp(scale(i))
      enddo
      do i=1,n
        pnt1=colpnt(i)
        pnt2=colpnt(i+1)-1
        do j=pnt1,pnt2
          rownzs(j)=rownzs(j)/scale(i)/scale(rowidx(j)+n)
        enddo
        obj(i)=obj(i)/scale(i)
        bounds(i)=bounds(i)*scale(i)
      enddo
      do i=1,m
        rhs(i)=rhs(i)/scale(i+n)
        bounds(i+n)=bounds(i+n)/scale(i+n)
      enddo
      return
      end
c
c ============================================================================
c
      subroutine scrow1(rowidx,colpnt,
     x maxi,scale,excld,ss,veclen,lognz)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      real*8 lognz(nz),maxi(mn),scale(mn),ss
      integer*4 rowidx(cfree),colpnt(n1),excld(mn),veclen(mn)
      common/numer/ tplus,tzer
      real*8        tplus,tzer
c ---------------------------------------------------------------------------
      integer*4 i,j,pnt1,pnt2
      real*8 sol
c ---------------------------------------------------------------------------
      ss=0
      do i=1,m
        maxi(i)=0.0d+0
      enddo
      do i=1,n
        if(excld(i).gt.-2)then
          sol=scale(i)
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            if(excld(rowidx(j)+n).gt.-2)then
              maxi(rowidx(j))=maxi(rowidx(j))+lognz(j)-sol
              ss=ss+abs(lognz(j)-sol-scale(rowidx(j)+n))
            endif
          enddo
        endif
      enddo
      do i=1,m
        scale(n+i)=maxi(i)/veclen(i+n)
      enddo
      return
      end
c
c ===========================================================================
c
      subroutine sccol1(colpnt,scale,
     x excld,rowidx,veclen,lognz)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      real*8 scale(mn),lognz(nz)
      integer*4 colpnt(n1),excld(mn),rowidx(cfree),veclen(mn)
      common/numer/ tplus,tzer
      real*8        tplus,tzer
c ---------------------------------------------------------------------------
      integer*4 i,j,pnt1,pnt2
      real*8 ma
c ---------------------------------------------------------------------------
      do i=1,n
        ma=0.0d+0
        if(excld(i).gt.-2)then
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            ma=ma+lognz(j)-scale(rowidx(j)+n)
          enddo
          scale(i)=ma/veclen(i)
        endif
      enddo
      return
      end
c
c ===========================================================================
c
      subroutine scrow2(rhs,bounds,rownzs,rowidx,
     x colpnt,maxi,scale,excld)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      common/numer/ tplus,tzer
      real*8        tplus,tzer
c
      real*8 rownzs(cfree),bounds(mn),rhs(m),maxi(m),scale(mn)
      integer*4 rowidx(cfree),colpnt(n1),excld(mn)
c ---------------------------------------------------------------------------
      integer*4 i,j,pnt1,pnt2,k
      real*8 sol
c ---------------------------------------------------------------------------
      do i=1,m
        maxi(i)=0
      enddo
      do i=1,n
        if(excld(i).gt.-2)then
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            k=rowidx(j)
            sol=abs(rownzs(j))
            if (maxi(k).lt.sol)maxi(k)=sol
          enddo
        endif
      enddo
      do i=1,m
        if(maxi(i).le.tzer)maxi(i)=1.0d+0
        scale(n+i)=maxi(i)*scale(n+i)
        rhs(i)=rhs(i)/maxi(i)
        bounds(i+n)=bounds(i+n)/maxi(i)
      enddo
      do i=1,n
        pnt1=colpnt(i)
        pnt2=colpnt(i+1)-1
        do j=pnt1,pnt2
          k=rowidx(j)
          rownzs(j)=rownzs(j)/maxi(k)
        enddo
      enddo
      return
      end
c
c ===========================================================================

c
      subroutine sccol2(bounds,rownzs,colpnt,obj,scale,
     x excld,rowidx)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      real*8 rownzs(cfree),bounds(mn),obj(n),scale(mn)
      integer*4 colpnt(n1),excld(mn),rowidx(cfree)
      common/numer/ tplus,tzer
      real*8        tplus,tzer
c ---------------------------------------------------------------------------
      integer*4 i,j,pnt1,pnt2
      real*8 sol,ma
c ---------------------------------------------------------------------------
      do i=1,n
        if(excld(i).gt.-2)then
          ma=0
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            if(excld(rowidx(j)+n).gt.-2)then
              sol=abs(rownzs(j))
              if (ma.lt.sol)ma=sol
            endif
          enddo
          if (ma.le.tzer)ma=1.0d+0
          scale(i)=ma*scale(i)
          do j=pnt1,pnt2
            rownzs(j)=rownzs(j)/ma
          enddo
          obj(i)=obj(i)/ma
          bounds(i)=bounds(i)*ma
        endif
      enddo
      return
      end
c
c ===========================================================================
c
      subroutine scalobj(obj,scobj,excld,objnor)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      real*8 obj(n),scobj,objnor
      integer*4 excld(n),i
      character*99 buff
c ---------------------------------------------------------------------------
      scobj=0.0d+0
      do i=1,n
        if(excld(i).gt.-2)then
          if (abs(obj(i)).gt.scobj)scobj=abs(obj(i))
        endif
      enddo
      scobj=scobj/objnor
      if(scobj.lt.1.0d-08)scobj=1.0d-08
      write(buff,'(1x,a,d8.2)')'Obj. scaled ',scobj
      call mprnt(buff)
      do i=1,n
        obj(i)=obj(i)/scobj
      enddo
      return
      end
c
c ===========================================================================
c
      subroutine scalrhs(rhs,scrhs,excld,rhsnor,bounds,xs,up )
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      real*8 rhs(m),scrhs,rhsnor,bounds(mn),xs(mn),up(mn)
      integer*4 excld(mn),i
      character*99 buff
c ---------------------------------------------------------------------------
      scrhs=0.0d+0
      do i=1,m
        if(excld(i+n).gt.-2)then
          if(abs(rhs(i)).gt.scrhs)scrhs=abs(rhs(i))
        endif
      enddo
      scrhs=scrhs/rhsnor
      if(scrhs.lt.1.0d-08)scrhs=1.0d-08
      write(buff,'(1x,a,d8.2)')'Rhs. scaled ',scrhs
      call mprnt(buff)
      do i=1,m
        rhs(i)=rhs(i)/scrhs
      enddo
      do i=1,mn
        bounds(i)=bounds(i)/scrhs
        xs(i)=xs(i)/scrhs
        up(i)=up(i)/scrhs
      enddo
      return
      end
c
c ============================================================================
c Curtis-Reid Scaling algorithm
c ============================================================================
c
      subroutine scale2(bounds,rownzs,colpnt,obj,sc,
     x  vcstat,rowidx,rhs,scpass,scdif,scm1,rk,logsum,count)

      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      common/numer/ tplus,tzer
      real*8        tplus,tzer
c
      real*8 bounds(mn),rownzs(cfree),obj(n),sc(mn),
     x  rhs(m),scdif,scm1(mn),rk(mn),logsum(mn)
      integer*4 rowidx(cfree),colpnt(n1),vcstat(mn),scpass,count(mn)
c
      integer*4 i,j,in,pnt1,pnt2,pass
      real*8 logdef,s,qk,qkm1,ek,ekm1,ekm2,sk,skm1
      character*99 buff
c
      pass=0
      do i=1,mn
       count(i)=0
       logsum(i)=0.0d+0
      enddo
      logdef=0.0d+0
      in=0
      do i=1,n
        if(vcstat(i).gt.-2)then
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            if(vcstat(rowidx(j)+n).gt.-2)then
              if(abs(rownzs(j)).gt.tzer)then
                s=log(abs(rownzs(j)))
                count(rowidx(j)+n)=count(rowidx(j)+n)+1
                count(i)=count(i)+1
                logsum(i)=logsum(i)+s
                logsum(rowidx(j)+n)=logsum(rowidx(j)+n)+s
                logdef=logdef+s*s
                in=in+1
              endif
            endif
          enddo
        endif
      enddo
      do i=1,mn
       if((vcstat(i).le.-2).or.(count(i).eq.0))count(i)=1
      enddo
      logdef=sqrt(logdef)/dble(in)
      logdef=exp(logdef)
      write(buff,'(1x,a,i2,a,d12.6)')'Pass',pass,'. Average def.',logdef
      call mprnt(buff)
      if(logdef.le.scdif)then
        do i=1,mn
          sc(i)=1.0d+0
        enddo
        goto 999
      endif
c
c Initialize
c
      do i=1,m
        sc(i+n)=logsum(i+n)/count(i+n)
        rk(i+n)=0
      enddo
      sk=0
      do i=1,n
        if(vcstat(i).gt.-2)then
          s=logsum(i)
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          do j=pnt1,pnt2
            s=s-logsum(rowidx(j)+n)/count(rowidx(j)+n)
          enddo
        else
          s=0
        endif
        rk(i)=s
        sk=sk+s*s/count(i)
        sc(i)=0.0d+0
      enddo
      do i=1,mn
        scm1(i)=sc(i)
      enddo
      ekm1=0
      ek=0
      qk=1.0d+0
c
c Curtis-Reid scaling
c
  10  pass=pass+1
        do i=1,m
          rk(i+n)=ek*rk(i+n)
        enddo
        do i=1,n
          if(vcstat(i).gt.-2)then
            pnt1=colpnt(i)
            pnt2=colpnt(i+1)-1
            s=rk(i)/count(i)
            do j=pnt1,pnt2
              if(vcstat(rowidx(j)+n).gt.-2)
     x        rk(rowidx(j)+n)=rk(rowidx(j)+n)+s
            enddo
          endif
        enddo
        skm1=sk
        sk=0.0d+0
        do i=1,m
          rk(i+n)=-rk(i+n)/qk
          sk=sk+rk(i+n)*rk(i+n)/count(i+n)
        enddo
        ekm2=ekm1
        ekm1=ek
        ek=qk*sk/skm1
        qkm1=qk
        qk=1-ek
        if(pass.gt.scpass)goto 20
c
c Update Column-scale factors
c
        do i=1,n
          if(vcstat(i).gt.-2)then
            s=sc(i)
            sc(i)=s+(rk(i)/count(i)+ekm1*ekm2*(s-scm1(i)))/qk/qkm1
            scm1(i)=s
          endif
        enddo
c
c even pass
c
        do i=1,n
          if(vcstat(i).gt.-2)then
            s=ek*rk(i)
            pnt1=colpnt(i)
            pnt2=colpnt(i+1)-1
            do j=pnt1,pnt2
              if(vcstat(rowidx(j)+n).gt.-2)
     x        s=s+rk(rowidx(j)+n)/count(rowidx(j)+n)
            enddo
            s=-s/qk
          else
            s=0
          endif
          rk(i)=s
        enddo
        skm1=sk
        sk=0.0d+0
        do i=1,n
          sk=sk+rk(i)*rk(i)/count(i)
        enddo
        ekm2=ekm1
        ekm1=ek
        ek=qk*sk/skm1
        qkm1=qk
        qk=1-ek
c
c Update Row-scale factors
c
        do i=1,m
          j=i+n
          if(vcstat(j).gt.-2)then
            s=sc(j)
            sc(j)=s+(rk(j)/count(j)+ekm1*ekm2*(s-scm1(j)))/qk/qkm1
            scm1(j)=s
          endif
        enddo
      goto 10
c
c Syncronize Column factors
c
  20  do i=1,n
        if(vcstat(i).gt.-2)then
          sc(i)=sc(i)+(rk(i)/count(i)+ekm1*ekm2*(sc(i)-scm1(i)))/qkm1
        endif
      enddo
c
c Scaling
c
      logdef=0
      do i=1,mn
        if(vcstat(i).gt.-2)then
          sc(i)=exp(sc(i))
        else
          sc(i)=1.0d+0
        endif
      enddo
      do i=1,n
        pnt1=colpnt(i)
        pnt2=colpnt(i+1)-1
        do j=pnt1,pnt2
          rownzs(j)=rownzs(j)/sc(i)/sc(rowidx(j)+n)
          if((vcstat(rowidx(j)+n).gt.-2).and.
     x      (abs(rownzs(j)).gt.tzer))then
            s=log(abs(rownzs(j)))
            logdef=logdef+s*s
          endif
        enddo
        obj(i)=obj(i)/sc(i)
        bounds(i)=bounds(i)*sc(i)
      enddo
      do i=1,m
        rhs(i)=rhs(i)/sc(i+n)
        bounds(i+n)=bounds(i+n)/sc(i+n)
      enddo
      logdef=sqrt(logdef)/dble(in)
      logdef=exp(logdef)
      pass=pass-1
      write(buff,'(1x,a,i2,a,d12.6)')'Pass',pass,'. Average def.',logdef
      call mprnt(buff)
 999  return
      end
c
c ============================================================================
c ===========================================================================
c
      subroutine stlamb(colpnt,vcstat,rowidx,cnt,fixn,dropn,p)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      integer*4 colpnt(n1),vcstat(mn),rowidx(nz),cnt(mn),
     x          fixn,dropn,p
c
      common/factor/ tpiv1,tpiv2,tabs,trabs,lam,tfind,order,supdens
      real*8         tpiv1,tpiv2,tabs,trabs,lam,tfind,order,supdens
c
      common/setden/ maxdense,densgap,setlam,denslen
      real*8         maxdense,densgap
      integer*4      setlam,denslen
c
      integer*4 i,j,pnt1,pnt2,cn,lcn,lcd,ndn,z,maxcn
      real*8    la
      character*99 buff
c
c ---------------------------------------------------------------------------
c

C CMSSW: Explicit initialization needed
      ndn=0

      write(buff,'(1X)')
      call mprnt(buff)
      do i=1,m
        cnt(i)=0
      enddo
      if((m-dropn).ge.(n-fixn))then
        cnt(1)=m-dropn-n+fixn
      endif
      maxcn=0
      do i=1,n
        if(vcstat(i).gt.-2)then
          pnt1=colpnt(i)
          pnt2=colpnt(i+1)-1
          cn=0
          do j=pnt1,pnt2
            if(vcstat(rowidx(j)).gt.-2)cn=cn+1
          enddo
          if(cn.gt.0)cnt(cn)=cnt(cn)+1
          vcstat(i)=cn
          if(maxcn.lt.cn)maxcn=cn
        endif
      enddo
      if(setlam.lt.0)goto 70
c
      cn =maxcn
      lcd=maxcn
      lcn=maxcn
      z=0
C CMSSW: Explicit integer conversion needed
      pnt1=int((n-fixn+m-dropn)*maxdense)
      pnt2=0
      if((m-dropn).ge.1.5*(n-fixn))then
        maxdense=1.0
      endif
      if((m-dropn).ge.2.5*(n-fixn))then
        lcn=1
        lcd=2
        goto 60
      endif
c
      do while ((pnt2.le.pnt1).and.(cn.gt.0))
        if(cnt(cn).eq.0)then
          z=z+1
        else
          if(z.gt.0)then
            if((densgap*cn*cn).le.(cn+z+1)*(cn+z+1))then
              lcd=cn+z+1
              lcn=cn
              ndn=pnt2
            endif
            z=0
          endif
          pnt2=pnt2+cnt(cn)
        endif
        cn=cn-1
      enddo
c
  60  write(buff,'(1X,A,I6)')'Largest sparse column length :',lcn
      call mprnt(buff)
      if((maxcn.le.denslen).or.(lcn.eq.maxcn))then
        write(buff,'(1X,A)')'Problem has no dense columns'
        call mprnt(buff)
        lcn=maxcn
      else
        write(buff,'(1X,A,I6)')'Smallest dense column length :',lcd
        call mprnt(buff)
        write(buff,'(1X,A,I6)')'Number of dense columns      :',ndn
        call mprnt(buff)
      endif
      la=lcn+0.5
      la=la/m
      write(buff,'(1X,A,F7.4)')'Computed density parameter   : ',la
      call mprnt(buff)
      if(la.gt.lam)then
        lam=la
      else
        write(buff,'(1X,A,F7.4)') 'Parameter reset to value    : ',lam
        call mprnt(buff)
      endif
  70  lam=lam*m
      p=1
      if((lam.ge.maxcn).and.(setlam.le.0))p=2
      if(supdens.le.lam)supdens=lam
c
      write(buff,'(1X)')
      call mprnt(buff)
      return
      end
c
c ===========================================================================
c ===========================================================================
c
      subroutine symmfo(inta1,pivots,ecolpnt,vcstat,
     x colpnt,rowidx,rowpnt,colindex,perm,invperm,
     x count,inta2,inta3,inta4,inta5,inta6,inta7,inta8,inta9,
     x inta10,inta11,nonzeros,l,oper,tfind,inta12,code)
c
      common/dims/ n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
      integer*4    n,n1,m,mn,nz,cfree,pivotn,denwin,rfree
c
      integer*4 inta1(mn),ecolpnt(mn),pivots(mn),vcstat(mn),
     x colpnt(n1),rowidx(cfree),rowpnt(mn),colindex(rfree),
     x perm(mn),invperm(mn),count(mn),inta2(mn),inta3(mn),inta4(mn),
     x inta5(mn),inta6(mn),inta7(mn),inta8(mn),inta9(mn),inta10(mn),
     x inta11(mn),tfind,inta12(mn),l,code
c
      real*8   nonzeros(cfree),oper
c
      integer*4 i,j,k,t1,tt1,t2,p1,p2,pnt,pnt1,pnt2,aatnz
      character*99 buff
c
c ---------------------------------------------------------------------------
c
    1 format(1x,'Building aat                 time:',f9.2,' sec')
    2 format(1x,'Building ordering    list    time:',f9.2,' sec')
    4 format(1x,'Symbolic factorisation       time:',f9.2,' sec')
    5 format(1x,'Total symbolic phase         time:',f9.2,' sec')
    6 format(1x,'Sub-diagonal nonzeros in aat     :',i9)
    7 format(1x,'Sub-diagonal nonzeros in L       :',i9)
    8 format(1x,'NONZEROS         :',i12)
    9 format(1x,'OPERATIONS       :',f13.0)
   10 format(1x,'Minimum Local Fill-in Ordering with Power:',i3)
   11 format(1x,'Minimum Degree Ordering (Power=0)')
   12 format(1x,'Without Ordering')
c
      call timer (tt1)
      if(tfind.lt.0)then
        write(buff,12)
      else if(tfind.eq.0)then
        write(buff,11)
      else
        write(buff,10)tfind
      endif
      oper=0.0d+0
      call mprnt(buff)
      do i=1,nz
        rowidx(i)=rowidx(i)-n
      enddo
      if(rfree.lt.nz)then
        write(buff,'(1x,a)')'Not enough integer memory'
        call mprnt(buff)
        code=-2
        goto 999
      endif
      if(cfree.lt.2*nz)then
        write(buff,'(1x,a)')'Not enough real memory'
        call mprnt(buff)
        code=-2
        goto 999
      endif
c
c If no ordering...
c
      if(tfind.lt.0)then
        t2=tt1
        do i=1,m
          perm(i)=i
        enddo
        goto 50
      endif
c
c Otherwise...
c
      do i=1,n
        inta2(i)=i
      enddo
      call transps(n,m,nz,colpnt,rowidx,nonzeros,
     x rowpnt,colindex,nonzeros(nz+1),inta2)
      k=1
      l=m
      do i=1,m
        pivots(i)=0
        if(vcstat(i+n).le.-2)then
          invperm(l)=i
          l=l-1
        else
          invperm(k)=i
          k=k+1
        endif
      enddo
      call transps(m,n,nz,rowpnt,colindex,nonzeros(nz+1),
     x colpnt,rowidx,nonzeros,invperm)
      do i=1,n
        p1=colpnt(i)
        if(vcstat(i).le.-2)then
          p2=colpnt(i)-1
        else
          p2=colpnt(i+1)-1
…