c $debug
c nlaux2a
c no setpar.inc // commons scale & bits no dimensions
c.....list of subroutines....
c confun	dummy !! remove or rename when using nonl. constraints
c infofit       AIC BIC CAIC
c gethes        central diff approx to hessian
c readm 	reads input file after first 4 lines
c lastwo 	simple error routine
c sortn		sorts parameter numbers (numerical recipe routine)
c fdh2 		grad check... =obsolete
c invert	invert nonsymm matrix 
c mvn 		calculate multivariate density
c dista         calc Yung distance
c .............................................................
      subroutine confun(mode,ncnln,n,nrowj,needc,x,cstr,cjac,nstate)
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      integer mode,nstate,n,ncnln,nrowj,needc(*)
      double precision x(n),cstr(nrowj),cjac(nrowj,n)
      character name(6)*2
c     common / scala / np,nv,nc,ne,itqn,item,iem
c     common / bits  / pi2,half,shalf,one,zero,two
c     common / names / name
      if (ncnln.eq.0) return
      if (nstate.eq.1) then
       do 5 i=1,nrowj
       do 5 j=1,n
  5    cjac(i,j)=0.d0 
      endif
      return
      end
      subroutine infofit(objf,npar,np)
c AIC BAIC CAIC
      double precision aic,bic,caic,objf,one,two
      integer npar,np
      data one,two / 1.d0,2.d0 /
      caic=two*(objf)+dble(npar)*(dlog(dble(np))+one) 
      bic=two*(objf)+dble(npar)*dlog(dble(np))
      aic=two*(objf)+two*dble(npar)
c     write(16,928) aic,bic,caic
 928  format(/,'AIC ',f10.3,' BIC ',f10.3,' CAIC 'f10.3)
      return 
      end
      subroutine gethes(hess,x,g,bu,bl,cjac,cn,npar,mpar,nhmax,
     &      nrowj,ncnln,nctotl,ntot,ifail)
c approx hessian
      implicit double precision (a-h,o-z),integer (i-n)
      double precision hess(nhmax,nhmax),cn(nrowj),cjac(nrowj,mpar)
      double precision bu(nctotl),bl(nctotl),x(mpar),g(mpar)
      double precision zero,two,det,one
      external fdh,cholinv,confun,invert
      data zero,two,one/0.d0,2.d0,1.d0/
      do 422 i=1,npar
      do 422 j=1,npar
 422  hess(i,j)=zero
      call fdh(x,g,hess,bu,npar,nhmax)
      do 111 i=1,npar-1
      do 111 j=i+1,npar
 111  hess(i,j)=(hess(i,j)+hess(j,i))/two 
      do 112 i=1,npar
      do 112 j=1,i
 112  hess(i,j)=hess(j,i)
c
      write(17,*) ' hessian dimension = ',npar
      do 218 i=1,ntot
 218  write(17,'(6(1x,f12.6))') (hess(i,j),j=1,i)
c     if no nonlinear constraints 
      if (ncnln.eq.0) then
        ifail=0
        call cholinv(hess,mpar,npar,tmp,ifail,bu)
         if (ifail.ne.0) return
      else 
c correction for nonlinear constraints
         call confun(mode,ncnln,npar,nrowj,needc,x,cn,cjac,1)
         do 201 i=1,nrowj
         do 201 j=1,npar
         hess(i+npar,j)=cjac(i,j)
 201     hess(j,i+npar)=hess(i+npar,j)
         do 202 i=1,npar
         do 202 j=1,npar
         do 202 k=1,nrowj
 202     hess(i,j)=hess(i,j)+cjac(k,i)*cjac(k,j)
c
         write(17,*) ' hessian dimension = ',ntot
         do 208 i=1,ntot
 208     write(17,'(6(1x,f12.6))') (hess(i,j),j=1,i)
c
         call invert(hess,ntot,det,nhmax)
            if (det.eq.zero) then 
              ifail=1
              return
            endif
       endif
      do 92 i=1,npar
      bu(i)=dsqrt(hess(i,i))
 92   bl(i)=x(i)/bu(i)
      return
      end
      subroutine hessing(ncnln)
c a messenger
      character txt*8
      write(*,*) ' HESSIAN NOT POS DEF......'
      if (ncnln.eq.0) then
      write(16,*) ' HESSIAN NOT POS DEF (cholinv) '
       else
      write(16,*) ' HESSIAN NOT POS DEF (invert) '
       endif
      return
      end 
      subroutine readm(ig,nv,ne,nc,idim,ipat,xmat,imat,ng,maxg,
     &           maxc,maxv,p,ip,iwork,liwork,name,npar,mpar,ifail)
c read matrices from input file
      implicit double precision (a-h,o-z),integer (i-n)
c
      integer idim(maxg,6,2),ipat(maxg,6),imat(maxg,maxc,6,maxv,maxv),
     & ip(maxg,maxc),nc(maxg)
      double precision xmat(maxg,maxc,6,maxv,maxv),p(maxg,maxc)
      integer iwork(liwork)
      character name(6)*2,tag*25,sp*25
      external writem
      data sp/'                        '/
c proportions
      p(ig,nc(ig))=1.d0
c I-B   
      do 10 ic=1,nc(ig)
      do 11 i1=1,ne  
      do 12 i2=1,ne
 12   xmat(ig,ic,2,i1,i2)=0.d0
 11   xmat(ig,ic,2,i1,i1)=1.d0
 10   continue
c count no. matrices used
      imt=0
      do 150 i=1,6
 150  if (ipat(ig,i).ne.0) imt=imt+1
c ly -- fill in dimensions
      idim(ig,1,1)=nv
      idim(ig,1,2)=ne
c be
      idim(ig,2,1)=ne
      idim(ig,2,2)=ne
c ps
      idim(ig,3,1)=ne
      idim(ig,3,2)=ne
c te
      idim(ig,4,1)=nv
      idim(ig,4,2)=nv
c ty or nu
      idim(ig,5,1)=nv
      idim(ig,5,2)=1
c al
      idim(ig,6,1)=ne
      idim(ig,6,2)=1
c read matrices first tag=name...................................
      do 2 ic=1,nc(ig)
      iim=0
 777  read(1,'(a2)') tag
      if (tag.eq.sp) goto 777
      iim=iim+1
c name identifies number ly=1 be=2 ps=3 te=4 nu=5 al=6 (fixed !)
      itag=-1
      do 40 i=1,6
      if (tag.eq.name(i)) itag=i
  40  continue
      if (itag.eq.-1) then 
         write(16,929) ig,ic,tag
 929     format(' group ',i2,' component ',i1,1x,a2,' unknown name ')
         ifail=-9
         return
      endif
c presence of name implies ipat(ig,itag).ne.0
          iend=idim(ig,itag,1)
          jend=idim(ig,itag,2)
          if (itag.eq.3.or.itag.eq.4) then 
           read(1,*,err=900)((xmat(ig,ic,itag,i,j),j=1,i),i=1,iend)
            do 565 i=1,iend
            do 565 j=1,i    
 565        xmat(ig,ic,itag,j,i)=xmat(ig,ic,itag,i,j)
            else
            read(1,*,err=900)((xmat(ig,ic,itag,i,j),j=1,jend),i=1,iend)
             endif
         if (ipat(ig,itag).eq.2) then
          iend=idim(ig,itag,1)
          jend=idim(ig,itag,2)
          if (itag.eq.3.or.itag.eq.4) then  
           read(1,*,err=901)((imat(ig,ic,itag,i,j),j=1,i),i=1,iend)
            do 568 i=1,iend
            do 568 j=1,i
 568        imat(ig,ic,itag,j,i)=imat(ig,ic,itag,i,j)
           else
            read(1,*,err=901)((imat(ig,ic,itag,i,j),j=1,jend),i=1,iend)
            endif
          endif

      if (imt.gt.iim) goto 777
      if (ic.lt.nc(ig)) then
       read(1,*,err=902) p(ig,ic)
       read(1,*,err=903) ip(ig,ic)
      endif
      if (nc(ig).eq.1) p(ig,1)=1.d0
 2    continue
c     write(16,*) ' ORIGINAL NUMBERING OF PARAMETERS '
c     call writem(1)
c error trap: if no parameters present, read to end-of-file and goto 444
c...................................... renumber parameters
c count parameters 
c
      if (ig.lt.ng) return
      icount=0
      do 200 igg=1,ng
      do 200 ic=1,nc(igg)
      do 200 im=1,6
      do 200 i=1,idim(igg,im,1)
      do 200 j=1,idim(igg,im,2)
      indi=imat(igg,ic,im,i,j)
      if (indi.eq.0) goto 200
c take into account equalities
       do 300 icou=1,icount
       if (iwork(icou).eq.indi) goto 200
 300  continue
c icount counts parameter and iwork contains original parameter number
      icount=icount+1
      iwork(icount)=indi     
 200  continue
c do the same for proportions
      do 210 igg=1,ng
      do 210 ic=1,nc(igg)-1
      indi=ip(igg,ic)
      if (indi.eq.0) goto 210
      do 310 icou=1,icount
      if (iwork(icou).eq.indi) goto 210
 310  continue
      icount=icount+1
      iwork(icount)=indi
 210  continue
c
c take into account the cases of zero parameters npar=zero
      if (icount.eq.0) then
       npar=0
       return
      endif
c icount=number of parameters, iwork = original numbering
      if (icount.gt.1) call sort(icount,iwork) 
c renumber
      do 400 igg=1,ng
      do 400 ic=1,nc(igg)
      do 400 im=1,6
      do 400 i=1,idim(igg,im,1)
      do 400 j=1,idim(igg,im,2)
      indi=imat(igg,ic,im,i,j)
c
       do 410 ii=1,icount
       if (indi.eq.iwork(ii)) then 
         imat(igg,ic,im,i,j)=ii
       endif
c
 410   continue
 400   continue
       do 500 igg=1,ng
       do 500 ic=1,nc(igg)
       indi=ip(igg,ic)
       do 510 ii=1,icount
       if (indi.eq.iwork(ii)) then
        ip(igg,ic)=ii
       endif
 510   continue
 500   continue
c  icount=npar..
      npar=icount
      return
c error
 900     write(16,800) name(itag),ic
         call lastwo
 901     write(16,801) name(itag),ic
         call lastwo
 902     write(16,802) ic
         call lastwo
 903     write(16,803) ic
         call lastwo
 800     format(' error reading fixed ',a2,' of component ',i1)
 801     format(' error reading free ',a2,' of component ',i1)
 802     format(' error reading fixed proportion in component ',i1)
 803     format(' error reading free proportion in component ',i1)
      return
      end
c
         subroutine lastwo
c report last to lines in event of certain errors
         implicit double precision (a-h,o-z),integer (i-n)
         character x*80
         backspace 1
         read(1,920,err=910) x    
         write(16,*) ' LAST LINE READ -1 '
         write(16,920) x
         read(1,920,err=910) x    
 920     format(a79)
         write(16,*) ' LAST LINE READ   '  
         write(16,920) x
 910     stop ' FATAL ERROR READING FILE - SEE OUTFILE '
         end    
      SUBROUTINE SORT(N,idim)
c sort when n>1
      integer idim(N),l,n,ir,rra,j,i
      L=N/2+1
      IR=N
10    CONTINUE
        IF(L.GT.1)THEN
          L=L-1
          RRA=idim(L)
        ELSE
          RRA=idim(IR)
          idim(IR)=idim(1)
          IR=IR-1
          IF(IR.EQ.1)THEN
            idim(1)=RRA
            RETURN
          ENDIF
        ENDIF
        I=L
        J=L+L
20      IF(J.LE.IR)THEN
          IF(J.LT.IR)THEN
            IF(idim(J).LT.idim(J+1))J=J+1
          ENDIF
          IF(RRA.LT.idim(J))THEN
            idim(I)=idim(J)
            I=J
            J=J+J
          ELSE
            J=IR+1
          ENDIF
        GO TO 20
        ENDIF
        idim(I)=RRA
      GO TO 10
      END
c
      subroutine fdh2(x,g,w,n)
c routine to check loglikelihood derivatives....obsolete..can be deleted
      implicit double precision (a-h,o-z),integer (i-n)  
      double precision sqrteta,stepsj,tempj,eps
      double precision w(n),x(n),g(n),f,f2
      integer j
      data eps/ 2.220446049250313D-016 /
      call objfun(-600,n,x,f,g,nstate)      
      sqrteta=dsqrt(eps)
      do 410 j=1,n
      stepsj=sqrteta*x(j)
      tempj=x(j)
      x(j)=x(j)+stepsj
      stepsj=x(j)-tempj
      if (stepsj.eq.0.d0) stepsj=sqrteta
      call objfun(-600,n,x,f2,g,nstate)      
      w(j)=(f2-f)/stepsj
 420  continue
      x(j)=tempj
 410  continue
      call objfun(0,n,x,f,g,nstate)      
      return
      end
      subroutine invert(help,n,det,m)
c..........................Cooley and Lohnes (1966)
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      double precision  help(m,n)
      double precision  one,zero,t,det
      integer           j,k,l
      data zero /0.d0/,one/1.d0/
c
      det=one 
      do 10 j=1,n
      t=help(j,j)
c     det=det*t
      help(j,j)=one 
      do 20 k=1,n
c............................check for singularity
      if (t.eq.zero) then
      det=zero 
      return
      endif
 20   help(j,k)=help(j,k)/t
      do 10 k=1,n
      if (k.eq.j) goto 10
      t=help(k,j)
      help(k,j)=zero
      do 30 l=1,n
 30   help(k,l)=help(k,l)-help(j,l)*t
 10   continue
      return
      end
c
      subroutine mvn(x,s,nv,det,tmp,f,maxv,tmp1)
c calculate mvn density
      implicit double precision (a-h,o-z),integer (i-n)
      double precision s(maxv,nv),x(maxv),f,det,tmp,wmis(maxv,maxv)
      double precision tmp1
      character name(6)*2
      common / bits /  pi2,half,shalf,one,zero,two
      common / names / name
      f=zero
      do 1 i=1,nv
      do 1 j=1,nv
  1   f=f+x(i)*s(i,j)*x(j)
      f=tmp1*tmp*dexp(shalf*f)
      return 
      end
      subroutine dista(cov,xm,work,wmat1,wvec1,maxv,maxc,nc,nv,lwork,
     &                 maxg,ig)
c Yung's distances between components, like Mahalonobis distances
      implicit integer(i-n), double precision(a-h,o-z)
      double precision cov(maxg,maxc,maxv,maxv),xm(maxg,maxc,maxv)
      double precision wmat1(maxv,maxv),tmp,wvec1(maxv),work(lwork)
      integer nc,nv,ic,jc,i,j,k,ind
c
      if (nc.eq.1) return
      write(16,981) 
  981 format(/,23('- '))
      ind=0 
      do 1 ic=1,nc
      do 2 i=1,nv
      do 2 j=1,i
      wmat1(i,j)=cov(ig,ic,i,j)
      wmat1(j,i)=wmat1(i,j)
 2    continue 
      call CHOLINV(wmat1,maxv,Nv,tmp,I,wvec1)
c I is ifail always 0, otherwise error would have occurred earlier
      do 1 jc=1,nc
      ind=ind+1
      if (ic.eq.jc) then
         work(ind)=0.d0
         goto 1
      endif
c
      tmp=0.d0
      do 3 i=1,nv
      do 3 j=1,nv
 3    tmp=tmp+(xm(ig,ic,i)-
     &   xm(ig,jc,i))*wmat1(i,j)*(xm(ig,ic,j)-xm(ig,jc,j)) 
      work(ind)=dsqrt(tmp)
 1    continue
      do 4 i=2,nc
      do 4 j=1,i-1  
      i1=nc*(i-1)+j
      i2=nc*(j-1)+i
      if (work(i1).ge.work(i2)) then
       write(16,999) ig,i,j,work(i1)
        else
       write(16,999) ig,j,i,work(i2)
      endif
  4   continue
 999  format(' GROUP ',i2,' COMPONENTS ',i1,1x,i1,
     &       ' Yung`s DISTANCE ',f5.2)
      write(16,922)             
 922  format(46('-'))
      return
      end