c $debug
c writem writes matrices
c objfun object function loglikelihood
      subroutine writem(itel)
      implicit double precision (a-h,o-z),integer (i-n)
      character name(6)*2,fmt1*13,fmt2*13
c
      parameter (maxg=2,maxn=10000,maxv=12,maxc=1,maxe=12,mpar=25)
      parameter (ncnln=0,nrowj=1)
c
      common / arra2 / cov(maxg,maxc,maxv,maxv),xm(maxg,maxc,maxv)
      common / wrkv /  wvec1(maxv),wvec2(maxv),wvec3(maxv)
      common / wrkm /  wmat1(maxv,maxv),wmat2(maxv,maxv),
     &                 wmat3(maxv,maxv),wmat4(maxv,maxv)
      common / xlis  / xmat(maxg,maxc,6,maxv,maxv),p(maxg,maxc)
      common / ilis  / imat(maxg,maxc,6,maxv,maxv),
     &                 idim(maxg,6,2),ip(maxg,maxc),ipat(maxg,6)
      common / scala / np,ng,npg(maxg),npc(0:maxg),nv(maxg),nc(maxg),
     &                 ne(maxg),itqn
      common / bits  / pi2,half,shalf,one,zero,two
      common / names / name
      data fmt1 /'(6(1x,f9.4))'/ 
      data fmt2 /'(6(1x,i10))'/
c itel = 0  write parameter matrices
c        1  write pattern matrices
c        2  write fitted cov and means
      if (itel.ge.2) write(16,888)
      do 1 ig=1,ng
      do 1 ic=1,nc(ig)
      write(16,988) ig,ic,nc(ig)
c not used
        do 2 im=1,6
           if (ipat(ig,im).eq.0.and.itel.eq.1) then
             if (im.gt.4) write(16,770) name(im)
             if (im.le.4) write(16,777) name(im)
             goto 2
         endif
        if (ipat(ig,im).gt.0.and.itel.eq.0) then
c fixed 
          if (im.le.4) write(16,776) name(im)                   
          if (im.gt.4) write(16,771) name(im)                   
c write fixed matrix 1 to 4
          if (im.lt.5) then
             do 30 i=1,idim(ig,im,1)
             indi=idim(ig,im,2)
             if (im.eq.3.or.im.eq.4) indi=i
             write(16,fmt1) (xmat(ig,ic,im,i,j),j=1,indi)
 30          continue
c mean vectors write fixed vector 5 to 6
               else
             write(16,fmt1) (xmat(ig,ic,im,i,1),i=1,idim(ig,im,1))
          endif
       endif
c pattern
      if (itel.eq.1.and.ipat(ig,im).eq.2) then
          write(16,775) name(im)
            if (im.lt.5) then
              do 31 i=1,idim(ig,im,1)
              indi=idim(ig,im,2)
              if (im.eq.3.or.im.eq.4) indi=i
              write(16,fmt2) (imat(ig,ic,im,i,j),j=1,indi)
 31           continue
               else
             write(16,fmt2) (imat(ig,ic,im,i,1),i=1,idim(ig,im,1))
           endif
      endif
  2   continue
c fitted
      if (itel.ge.2) then
         write(16,989) p(ig,ic),p(ig,ic)*dble(npg(ig))
         do 39 i=1,nv(ig)
         write(16,fmt1) (cov(ig,ic,i,j),j=1,i)
         do 39 j=1,i
         wmat1(i,j)=cov(ig,ic,i,j)
         wmat1(j,i)=cov(ig,ic,i,j)
 39      continue
c        call repdet(wmat1,maxv,nv(ig),wvec1)
         if (nv(ig).gt.1) then
          write(16,*) 'MATRIX SIGMA STANDARDIZED  '
          do 60 i=1,nv(ig)
          do 61 j=1,i
 61       wvec1(j)=cov(ig,ic,i,j)/dsqrt(cov(ig,ic,i,i)*cov(ig,ic,j,j))
 60       write(16,fmt1) (wvec1(j),j=1,i)
        endif
      endif
c moments
      if (itel.ge.2) then
        write(16,*) 'VECTOR of MEANS '
        write(16,fmt1) (xm(ig,ic,i),i=1,nv(ig))
        write(16,*) 'STANDARD DEVIATIONS '
        write(16,fmt1) (dsqrt(cov(ig,ic,i,i)),i=1,nv(ig))
      endif
      if (ic.eq.nc(ig)) goto 1
c
      if (itel.eq.0.or.itel.eq.3) then
        write(16,991) ic,p(ig,ic)
      endif
      if (itel.eq.1) then 
        if (ip(ig,ic).ne.0) then
        write(16,*) 'PROPORTION FREE PARAMETER (RE)NUMBERED  '
        write(16,fmt2) ip(ig,ic)
        endif
      endif
   1  continue
c
 775  format(a2,' FREE PARAMETERS (RE)NUMBERED  ')
 776  format('MATRIX ',a2)
 771  format('VECTOR ',a2)
 777  format('MATRIX ',a2,' NOT USED ')
 770  format('VECTOR ',a2,' NOT USED ')
 988  format(30('-'),'GROUP ',i2,' COMPONENT ',i2,' OF ',i2)
 888  format(/,'FITTED COVARIANCE MATRIX AND MEANS ')
 989  format('MATRIX SIGMA (P=',f5.3,' N=',f6.1,')')
 991  format('PROPORTION OF COMPONENT ',i1,' = ',f8.3)
      return 
      end
c
      subroutine objfun(iflag,npar,x,objf,g,nstate)      
c objfun for loglikelihood
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
c
      parameter (maxg=2,maxn=10000,maxv=12,maxc=1,maxe=12,mpar=25)
      parameter (ncnln=0,nrowj=1)
c
      parameter (mvt=maxc+maxv+1)
      double precision x(npar),g(npar)
      double precision objf,tmp,f,crit
c
      double precision binv(maxc,maxe,maxe)
      double precision bpsb(maxc,maxe,maxe)
      double precision binva(maxc,maxe)
      double precision binvl(maxc,maxv,maxe)
      double precision dxm(maxv),dcov(maxv,maxv),wcov(maxc,maxv,maxv)
      character mess*5,name(6)*2
c
c     common / summ / sigma(maxg,maxv,maxv),xmean(maxg,maxv)
      common / arra1 / obs(maxn,mvt),xmis,wmis(maxv,maxv)
      common / arra2 / cov(maxg,maxc,maxv,maxv),xm(maxg,maxc,maxv)
      common / iarra1 / misindi(maxn),misgrp(maxg)
c
      common / wrkm / wmat1(maxv,maxv),wmat2(maxv,maxv),
     &                wmat3(maxv,maxv),wmat4(maxv,maxv)
      common / wrkv / wvec1(maxv),wvec2(maxv),wvec3(maxv)
c
      common / xlis / xmat(maxg,maxc,6,maxv,maxv),p(maxg,maxc)
      common / ilis / imat(maxg,maxc,6,maxv,maxv),
     &                idim(maxg,6,2),ip(maxg,maxc),ipat(maxg,6)
      common / scala / np,ng,npg(maxg),npc(0:maxg),nv(maxg),
     &                 nc(maxg),ne(maxg),itqn
      common / bits  / pi2,half,shalf,one,zero,two
      common / names / name
c set to zero
      ifreq=25
      objf=zero
      do 203 j=1,ng
      do 203 i=npc(j-1)+1,npc(j)
 203  obs(i,nv(j)+nc(j)+1)=zero
      objf=zero
      if (itqn.ge.0) itqn=itqn+1
      do 1 i=1,npar
 1    g(i)=zero
c
c put parameters in place (proportions)
c
      do 1000 ig=1,ng
c
      p(ig,nc(ig))=one
c
      do 25 ic=1,nc(ig)-1
      indi=ip(ig,ic)
      if (indi.eq.0) goto 26
        p(ig,ic)=x(indi)
 26   p(ig,nc(ig))=p(ig,nc(ig))-p(ig,ic)
 25   continue
c
c stupid
      if (p(ig,nc(ig)).lt.zero) then
        mess='prop!'
        ifail=1
c       goto 444
      endif
c put parameters in place....
      call ppip(ig,nc,maxg,maxc,maxv,imat,ipat,idim,xmat,x,npar)
c calculate cov and xm 
      call mks(ig,maxg,nc,maxv,maxc,maxe,nv(ig),ne(ig),wmat1,
     &               xmat,binv,bpsb,binvl,cov,binva,xm)
c
      if (iflag.eq.-200.or.iflag.eq.-300) goto 1000
c
c get loglikelihood......
c
c nc>1 group ig is a mixture
       call getl(ig,maxg,nc,maxc,maxv,maxe,wmat1,wmat2,cov,wmis,
     &   wvec1,wvec2,p,wcov,obs,xmis,misindi,
     &   npg,npc,nv(ig),maxn,mvt,np,xm,objf,ifail,mess)
c check ifail and penalize if necessary
 444        if (ifail.ne.0) then 
c error ! punish !
c              if (mess.eq.'obs=0') write(16,999)(x(i),i=1,npar)
c999           format(5(5(1x,f10.5))/)
               write(*,777) itqn,mess 
 777           format(61x,i3,1x,a5)
               objf=100000.d0
               do 666 i=1,npar
 666           g(i)=100000.d0
               return
             endif
c -600  ... to check gradients, to get logl
       if (iflag.eq.-600) goto 1000
c
c derivatives  derivatives chains chains chains 
c *** p:  dL/dp logl with respect to proportion
c
       do 50 ic=1,nc(ig)-1
       do 50 is=npc(ig-1)+1,npc(ig)
       indi=ip(ig,ic)
       if (indi.ne.0) g(indi)=g(indi)+
     &      (obs(is,nv(ig)+ic)-obs(is,nv(ig)+nc(ig)))/
     &       obs(is,nv(ig)+nc(ig)+1)
  50   continue
c group ig, 
       do 52 ic=1,nc(ig)
c getdd1 dxm and dcov for non missing
       do 55 i=1,nv(ig)
       dxm(i)=0.d0
       do 55 j=1,nv(ig)
 55    dcov(i,j)=0.d0
       isnomis=0
       do 53 is=npc(ig-1)+1,npc(ig)
       if (misindi(is).eq.0) then
       isnomis=1
       call getdd1(ig,ic,is,maxn,mvt,maxv,obs,wvec1,wcov,xm,wmat1,
     &             dxm,dcov,p,nv(ig),np,nc,npc,maxc,maxg)
       endif
  53   continue
       if (isnomis.eq.1) then
       call getdm(ig,ic,nv(ig),ne(ig),maxg,maxc,maxv,maxe,xmat,imat,
     & ipat,wmat1,wmat2,wvec2,
     & dcov,dxm,binvl,binva,bpsb,npar,g,one)
       endif
       if (misgrp(ig).eq.1) then
       do 58 is=npc(ig-1)+1,npc(ig)
       if (misindi(is).gt.0) then
       call getdd2(ig,ic,is,maxn,mvt,maxv,obs,xmis,
     &             wmis,wvec1,wvec2,wvec3,
     &             wcov,xm,wmat1,wmat2,cov,
     &             dxm,dcov,p,nv(ig),nvobs,np,nc,npc,maxc,maxg)
       call getdm2(ig,ic,nv(ig),nvobs,ne(ig),
     & maxg,maxc,maxv,maxe,xmat,imat,
     & ipat,wmis,wmat1,wmat2,wmat3,wmat4,wvec1,wvec2,
     & dcov,dxm,binvl,binva,bpsb,npar,g,one)
       endif
 58    continue
       endif
c
  52  continue
c next group
1000  continue
c oridnary ml
      objf=-objf
      gl=zero
      do 80 i=1,npar
      gl=gl+g(i)*g(i)
 80   g(i)=-g(i)
c     if (mod(itqn,ifreq).eq.0.and.itqn.gt.0) write(*,888) itqn,objf,gl
 888  format(' qn it ',i4,1x,' logl ',e14.8,' len g ',e14.8)
      return
      end