c $debug
c cholinv 	invert symm pos def matrix
c cholesk 	cholesky decomposition
c fdh 		finite difference approx of hessian
c getdm		derivatives
c mks 		calculate model covariance matrices and mean vectors
c ppip 		put parameters in place
c getl 		calculate loglikelihood
c getdd1	derivatives....
c getpp   	calculates apost. assignment of subjects to components
c reader        reads part of input file (first 4 lines)
c repdet        gets det and writes det
c teltime       tells time.......
c check1        checks things to see if things are ok.....
      SUBROUTINE CHOLINV(A,IA,N,DET,IFAIL,P)
C invert pos def symm matrix by means of cholesky decomp
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DOUBLE PRECISION A(IA,N),DET,P(N),ZERO
      INTEGER          IFAIL,I,J,K
C
      DATA ZERO/ 0.D0 /
      IFAIL=0
      DET=ZERO
c do the cholesky decomp
      CALL CHOLESK(N,A,IA,P,DET,IFAIL)
      IF (IFAIL.NE.0) RETURN
      DO 10 I=1,N
      DO 11 J=1,I
  11  A(I,J)=ZERO
  10  A(I,I)=P(I)
      DO 15 J=1,N-1
      DO 15 I=J+1,N
      A(I,J)=ZERO
      DO 15 K=J,I-1
      A(I,J)=A(I,J)-P(I)*(A(K,I)*A(K,J))
  15  CONTINUE
C................LOWER TRIANGLE A CONTAINS THE INVERSE OF CHOLESKY INPUT
      DO 17 J=1,N-1
      DO 17 I=J+1,N
      A(J,I)=ZERO
      DO 18 K=I,N
  18  A(J,I)=A(J,I)+A(K,I)*A(K,J)
  17  CONTINUE
      DO 20 I=1,N
      P(I)=ZERO
      DO 21 K=I,N
  21  P(I)=P(I)+A(K,I)*A(K,I)
  20  A(I,I)=P(I)
C.......................UPPER TRIANGLE A CONTAINS THE INVERSE OF INPUT A
      DO 22 I=1,N
      DO 22 J=1,I
 22   A(I,J)=A(J,I)
      RETURN
      END
C-----------------------------------------------------------------------
      SUBROUTINE CHOLESK(N,A,IA,P,DET,IFAIL)
C cholesky decomp of pos def symm matrix
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      IMPLICIT INTEGER (I-N)
      DOUBLE PRECISION A(IA,N),DET,P(N),ZERO,ONE,TMP
      INTEGER          IFAIL,I,J,K
C
      DATA ONE /1.D0/, ZERO/ 0.D0 /
      IFAIL=1
C................................................*** 1) ELEMENT 1,1
      IF (A(1,1).LE.ZERO) RETURN
      P(1)=ONE/DSQRT(A(1,1))
C................................................*** 2) FIRST COLUMN
      DO 20 I=2,N
 20   A(1,I)=P(1)*A(1,I)
C................................................*** 3) NOW THE REST
      DO 30 I=2,N
      TMP=ZERO
      DO 35 K=1,I-1
 35   TMP=TMP+A(K,I)*A(K,I)
C......................................... ***  CHECK FOR SIGULARITY
      IF (TMP.GE.A(I,I)) RETURN
C
      P(I)=ONE/DSQRT(A(I,I)-TMP)
      DO 40 J=I,N
      IF (I.EQ.J) GOTO 40
      A(I,J)=P(I)
      TMP=ZERO
      DO 45 K=1,I-1
 45   TMP=TMP+A(K,J)*A(K,I)
      A(I,J)=A(I,J)*(A(J,I)-TMP)
 40   CONTINUE
 30   CONTINUE
C.................................................GET DETERMINANT OF A
      DET=ONE
      DO 50 I=1,N
      TMP=ONE/P(I)
 50   DET=DET*TMP*TMP
      IFAIL=0
      RETURN
      END
c
      subroutine fdh(x,g,hh,w,n,nn)
      implicit double precision (a-h,o-z),integer (i-n)  
c approximate hessian using exact gradients =  my central differences
      double precision sqrteta,stepsj,tempj,eps
      double precision hh(nn,n),w(n),x(n),g(n),f
      integer i,j
      eps=2.220446049250313d-16
      sqrteta=dsqrt(eps)
c     call objfun(1,n,x,f,g,nstate)      
      do 410 j=1,n
      stepsj=sqrteta*x(j)
      tempj=x(j)
      x(j)=tempj+stepsj
      stepsj=x(j)-tempj
      if (stepsj.eq.0.d0) stepsj=sqrteta
      call objfun(1,n,x,f,w,nstate)      
      x(j)=tempj-stepsj
      call objfun(1,n,x,f,g,nstate)      
      do 420 i=1,n
      hh(i,j)=((w(i)-g(i))/(2.d0*stepsj))
 420  continue
      x(j)=tempj
 410  continue
      return
      end
      subroutine getdm2(ig,ic,nv,nvobs,
     & ne,maxg,maxc,maxv,maxe,xmat,imat,
     & ipat,wmis,wmat1,wmat2,binvl2,xl2,xmisi,wvec2,dcov,dxm,
     & binvl,binva,bpsb,npar,g,weight)
c
c derivatives: d(SIGMA)/dM and d(MEANS)/dM where
c              M=Lambda, Psy, Theta, I-Beta, alfa or nu
c
      implicit double precision (a-h,o-z),integer (i-n)
      double precision bpsb(maxc,maxe,maxe),g(npar)
      double precision binva(maxc,maxe)
      double precision binvl(maxc,maxv,maxe)
      double precision dxm(maxv),dcov(maxv,maxv)
      double precision wmat1(maxv,maxv),wmat2(maxv,maxv),wvec2(maxv)
      double precision xmat(maxg,maxc,6,maxv,maxv)
      double precision weight,wmis(maxv,maxv)
c actually wmat3 wamt4 wvec3
      double precision binvl2(maxv,maxv),xmisi(maxv),xl2(maxv,maxv)
      integer ipat(maxg,6),imat(maxg,maxc,6,maxv,maxv),i,j,indi,iv,jv,k
      character name(6)*2
      common / bits  / pi2,half,shalf,one,zero,two
      common / names / name
c
      do 10 i=1,nv
      xmisi(i)=0.d0
      do 10 j=1,nvobs
      xmisi(i)=xmisi(i)+wmis(j,i)
 10   continue
      do 1 i=1,nvobs
      do 1 j=1,ne
      binvl2(i,j)=0.d0
      xl2(i,j)=0.d0
      do 1 k=1,nv
      xl2(i,j)=xl2(i,j)+wmis(i,k)*xmat(ig,ic,1,k,j) 
 1    binvl2(i,j)=binvl2(i,j)+wmis(i,k)*binvl(ic,k,j)
c  *********************** dM/d(B)
      if (ipat(ig,2).eq.2) then
          do 884 i=1,nvobs
          do 884 j=1,ne
 884      wmat2(i,j)=dxm(i)*xmat(ig,ic,6,j,1)
          do 881 i=1,ne
          do 881 j=1,ne
          wmat1(i,j)=zero
          do 881 k=1,nvobs
 881      wmat1(i,j)=wmat1(i,j)+binvl2(k,i)*wmat2(k,j)
          do 883 i=1,ne
          do 883 j=1,ne
          indi=imat(ig,ic,2,i,j)
          if (indi.eq.0) goto 883
          g(indi)=g(indi)-wmat1(i,j)*weight
 883      continue
        endif
c  *********************** dM/d(Lambda)
      if (ipat(ig,1).eq.2) then        
          do 151 i=1,nvobs
          do 151 j=1,ne
 151       wmat1(i,j)=dxm(i)*binva(ic,j)
          ii=0
          do 154 i=1,nv
          if (xmisi(i).eq.0.d0) goto 154
          ii=ii+1
          do 153 j=1,ne
          indi=imat(ig,ic,1,i,j)
          if (indi.eq.0) goto 153
          g(indi)=g(indi)+wmat1(ii,j)*weight
 153       continue 
 154       continue
      endif
c **************************  dM/d(nu)
      if (ipat(ig,5).eq.2) then
        iiv=0
        do 160 iv=1,nv
        if (xmisi(iv).eq.0.d0) goto 160
          iiv=iiv+1
          indi=imat(ig,ic,5,iv,1)
          if (indi.ne.0) then
           g(indi)=g(indi)+dxm(iiv)*weight
          endif
 160    continue
      endif
c ****************************  dM/d(alfa)
      if (ipat(ig,6).eq.2) then
          do 162 i=1,ne
          wvec2(i)=zero
          do 162 j=1,nvobs
  162     wvec2(i)=wvec2(i)+binvl2(j,i)*dxm(j)
          do 163 i=1,ne
          indi=imat(ig,ic,6,i,1)
          if (indi.eq.0) goto 163
          g(indi)=g(indi)+wvec2(i)*weight
 163      continue
      endif
c *******************************  d(cov)/d(te)
      if (ipat(ig,4).eq.2) then
         iiv=0
         do 170 iv=1,nv
         if (xmisi(iv).eq.0.d0) goto 170
         iiv=iiv+1
         jjv=0
         do 171 jv=1,nv 
         if (xmisi(jv).eq.0.d0) goto 171
         jjv=jjv+1
         indi=imat(ig,ic,4,iv,jv)
            if (indi.ne.0) then
             g(indi)=g(indi)+dcov(iiv,jjv)*weight
            endif
 171     continue
 170     continue
      endif
c ******************************* d(cov)/d(ps)
      if (ipat(ig,3).eq.2) then
               do 175 i=1,ne
               do 175 j=1,nvobs
               wmat2(i,j)=zero
               do 175 k=1,nvobs
 175           wmat2(i,j)=wmat2(i,j)+binvl2(k,i)*dcov(k,j)
               do 174 i=1,ne
               do 174 j=1,ne
               wmat1(i,j)=zero
               do 174 k=1,nvobs
 174           wmat1(i,j)=wmat1(i,j)+wmat2(i,k)*binvl2(k,j)
               do 176 i=1,ne
               do 176 j=1,ne
               indi=imat(ig,ic,3,i,j)
               if (indi.eq.0) goto 176            
               g(indi)=g(indi)+wmat1(i,j)*weight
 176           continue
      endif
c *******************************  d(cov)/d(ly)
      if (ipat(ig,1).gt.0) then
c note ipat(1)>0 because part of the next derivative (d(cov)/d(be))
c                                                is calculated here
         do 191 i=1,nvobs
         do 191 j=1,ne
         wmat1(i,j)=zero
         do 191 k=1,nvobs
 191     wmat1(i,j)=wmat1(i,j)+dcov(i,k)*xl2(k,j)
         do 192 i=1,nvobs
         do 192 j=1,ne
         wmat2(i,j)=zero
         do 192 k=1,ne
 192     wmat2(i,j)=wmat2(i,j)+wmat1(i,k)*bpsb(ic,k,j)
         ii=0
         do 193 i=1,nv
         if (xmisi(i).eq.0.d0) goto 193
         ii=ii+1
         do 194 j=1,ne
         indi=imat(ig,ic,1,i,j)
            if (indi.eq.0) goto 194
         g(indi)=g(indi)+wmat2(ii,j)*two*weight
 194     continue
 193     continue
      endif
c  ******************************* d(cov)/d(be)
      if (ipat(ig,2).eq.2) then
       do 481 i=1,ne
       do 481 j=1,ne
       wmat1(i,j)=zero
       do 481 k=1,nvobs
 481   wmat1(i,j)=wmat1(i,j)+binvl2(k,i)*wmat2(k,j)
       do 482 i=1,ne
       do 482 j=1,ne
        indi=imat(ig,ic,2,i,j)
        if (indi.eq.0) goto 482
        g(indi)=g(indi)-wmat1(i,j)*two*weight
 482   continue
      endif
      return 
      end

      subroutine getdm(ig,ic,nv,ne,maxg,maxc,maxv,maxe,xmat,imat,
     & ipat,wmat1,wmat2,wvec2,dcov,dxm,binvl,binva,bpsb,npar,g,weight)
c derivatives: d(SIGMA)/dM and d(MEANS)/dM where
c              M=Lambda, Psy, Theta, I-Beta, alfa or nu
      implicit double precision (a-h,o-z),integer (i-n)
      double precision bpsb(maxc,maxe,maxe),g(npar)
      double precision binva(maxc,maxe)
      double precision binvl(maxc,maxv,maxe)
      double precision dxm(maxv),dcov(maxv,maxv)
      double precision wmat1(maxv,maxv),wmat2(maxv,maxv),wvec2(maxv)
      double precision xmat(maxg,maxc,6,maxv,maxv)
      double precision weight
      integer ipat(maxg,6),imat(maxg,maxc,6,maxv,maxv),i,j,indi,iv,jv,k
      character name(6)*2
      common / bits  / pi2,half,shalf,one,zero,two
      common / names / name
c dM/dB
      if (ipat(ig,2).eq.2) then
          do 884 i=1,nv
          do 884 j=1,ne
 884      wmat2(i,j)=dxm(i)*xmat(ig,ic,6,j,1)
          do 881 i=1,ne
          do 881 j=1,ne
          wmat1(i,j)=zero
          do 881 k=1,nv
 881      wmat1(i,j)=wmat1(i,j)+binvl(ic,k,i)*wmat2(k,j)
          do 883 i=1,ne
          do 883 j=1,ne
          indi=imat(ig,ic,2,i,j)
          if (indi.eq.0) goto 883
          g(indi)=g(indi)-wmat1(i,j)*weight
 883      continue
        endif
c  *********************** dM/d(Lambda)
      if (ipat(ig,1).eq.2) then        
          do 151 i=1,nv
          do 151 j=1,ne
151       wmat1(i,j)=dxm(i)*binva(ic,j)
          do 153 i=1,nv
          do 153 j=1,ne
          indi=imat(ig,ic,1,i,j)
          if (indi.eq.0) goto 153
          g(indi)=g(indi)+wmat1(i,j)*weight
153       continue 
      endif
c **************************  dM/d(nu)
      if (ipat(ig,5).eq.2) then
        do 160 iv=1,nv
        indi=imat(ig,ic,5,iv,1)
         if (indi.ne.0) then
          g(indi)=g(indi)+dxm(iv)*weight
         endif
 160    continue
      endif
c ****************************  dM/d(alfa)
      if (ipat(ig,6).eq.2) then
          do 162 i=1,ne
          wvec2(i)=zero
          do 162 j=1,nv
  162     wvec2(i)=wvec2(i)+binvl(ic,j,i)*dxm(j)
          do 163 i=1,ne
          indi=imat(ig,ic,6,i,1)
          if (indi.eq.0) goto 163
          g(indi)=g(indi)+wvec2(i)*weight
 163      continue
      endif
c *******************************  d(cov)/d(te)
      if (ipat(ig,4).eq.2) then
         do 170 iv=1,nv
         do 170 jv=1,nv 
         indi=imat(ig,ic,4,iv,jv)
            if (indi.ne.0) then
             g(indi)=g(indi)+dcov(iv,jv)*weight
            endif
 170     continue
      endif
c ******************************* d(cov)/d(ps)
      if (ipat(ig,3).eq.2) then
               do 175 i=1,ne
               do 175 j=1,nv
               wmat2(i,j)=zero
               do 175 k=1,nv
 175           wmat2(i,j)=wmat2(i,j)+binvl(ic,k,i)*dcov(k,j)
               do 174 i=1,ne
               do 174 j=1,ne
               wmat1(i,j)=zero
               do 174 k=1,nv
 174           wmat1(i,j)=wmat1(i,j)+wmat2(i,k)*binvl(ic,k,j)
               do 176 i=1,ne
               do 176 j=1,ne
               indi=imat(ig,ic,3,i,j)
               if (indi.eq.0) goto 176            
               g(indi)=g(indi)+wmat1(i,j)*weight
 176           continue
      endif
c *******************************  d(cov)/d(ly)
      if (ipat(ig,1).gt.0) then
c note ipat(1)>0 because part of the next derivative (d(cov)/d(be))
c                                                is calculated here
         do 191 i=1,nv
         do 191 j=1,ne
         wmat1(i,j)=zero
         do 191 k=1,nv
 191     wmat1(i,j)=wmat1(i,j)+dcov(i,k)*xmat(ig,ic,1,k,j)
         do 192 i=1,nv
         do 192 j=1,ne
         wmat2(i,j)=zero
         do 192 k=1,ne
 192     wmat2(i,j)=wmat2(i,j)+wmat1(i,k)*bpsb(ic,k,j)
         do 193 i=1,nv
         do 193 j=1,ne
         indi=imat(ig,ic,1,i,j)
            if (indi.eq.0) goto 193
         g(indi)=g(indi)+wmat2(i,j)*two*weight
 193     continue
      endif
c  ******************************* d(cov)/d(be)
      if (ipat(ig,2).eq.2) then
       do 481 i=1,ne
       do 481 j=1,ne
       wmat1(i,j)=zero
       do 481 k=1,nv
 481   wmat1(i,j)=wmat1(i,j)+binvl(ic,k,i)*wmat2(k,j)
       do 482 i=1,ne
       do 482 j=1,ne
        indi=imat(ig,ic,2,i,j)
        if (indi.eq.0) goto 482
        g(indi)=g(indi)-wmat1(i,j)*two*weight
 482   continue
      endif
      return 
      end
c
      subroutine mks(ig,maxg,nc,maxv,maxc,maxe,nv,ne,wmat1,
     &               xmat,binv,bpsb,binvl,cov,binva,xm)
c calculates (SIGMA) and (MEANS) given lisrel model
      implicit double precision (a-h,o-z)
      implicit integer (i-n)
      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 xmat(maxg,maxc,6,maxv,maxv)
      double precision cov(maxg,maxc,maxv,maxv),xm(maxg,maxc,maxv)
      double precision wmat1(maxv,maxv)
      integer nc(maxg)
      character name(6)*2
      common / bits  / pi2,half,shalf,one,zero,two
      common / names / name
c
      do 1 ic=1,nc(ig)
      do 700 i=1,ne           
      do 700 j=1,ne         
 700  wmat1(i,j)=xmat(ig,ic,2,i,j) 
      call invert(wmat1,ne,tmp,maxv)
      if (tmp.eq.zero) stop ' strange: (I-B)~ is singular '
      do 701 i=1,ne
      do 701 j=1,ne
 701  binv(ic,i,j)=wmat1(i,j)
c get (I-B)~*PS*(I-B)~`
      do 702 i=1,ne
      do 702 j=1,ne
      wmat1(i,j)=zero
      do 702 k=1,ne
 702  wmat1(i,j)=wmat1(i,j)+binv(ic,i,k)*xmat(ig,ic,3,k,j)
      do 703 i=1,ne
      do 703 j=1,ne
      bpsb(ic,i,j)=zero
      do 703 k=1,ne
  703 bpsb(ic,i,j)=bpsb(ic,i,j)+wmat1(i,k)*binv(ic,j,k) 
c get LY(I-B)~*PS*(I-B)~`
      do 130 i=1,nv
      do 130 j=1,ne
      wmat1(i,j)=zero
      binvl(ic,i,j)=zero
      do 130 k=1,ne
      binvl(ic,i,j)=binvl(ic,i,j)+xmat(ig,ic,1,i,k)*binv(ic,k,j) 
 130  wmat1(i,j)=wmat1(i,j)+xmat(ig,ic,1,i,k)*bpsb(ic,k,j)
c get LY (I-B)~*PS*(I-B)~`LY`
      do 131 i=1,nv
      do 131 j=1,nv
      cov(ig,ic,i,j)=zero
      do 131 k=1,ne
 131  cov(ig,ic,i,j)=cov(ig,ic,i,j)+wmat1(i,k)*xmat(ig,ic,1,j,k)
c get LY (I-B)~*PS*(I-B)~`LY`+TE
      do 132 i=1,nv
      do 132 j=1,nv
 132  cov(ig,ic,i,j)=cov(ig,ic,i,j)+xmat(ig,ic,4,i,j)
c get (I-B)~*AL
      do 133 i=1,ne
      binva(ic,i)=zero
      do 133 j=1,ne
 133  binva(ic,i)=binva(ic,i)+binv(ic,i,j)*xmat(ig,ic,6,j,1)
c get LY (I-B)~*AL
      do 735 i=1,nv
      xm(ig,ic,i)=zero
      do 735 j=1,ne
 735  xm(ig,ic,i)=xm(ig,ic,i)+xmat(ig,ic,1,i,j)*binva(ic,j)
c get NU + LY (I-B)~*AL
      do 134 i=1,nv
 134  xm(ig,ic,i)=xm(ig,ic,i)+xmat(ig,ic,5,i,1)
  1   continue
      return
      end
      subroutine ppip
     &(ig,nc,maxg,maxc,maxv,imat,ipat,idim,xmat,x,npar)
c put parameters in x into model matrices 
      implicit double precision (a-h,o-z), integer (i-n)
      integer ipat(maxg,6),imat(maxg,maxc,6,maxv,maxv),idim(maxg,6,2)
      integer ii,indi,i,j,nc(maxg)
      double precision x(npar),xmat(maxg,maxc,6,maxv,maxv)
      do 1 ic=1,nc(ig)
      do 10 ii=1,6
       if (ipat(ig,ii).eq.2) then
         do 11 i=1,idim(ig,ii,1)
         do 11 j=1,idim(ig,ii,2)
         indi=imat(ig,ic,ii,i,j)
          if (indi.ne.0) xmat(ig,ic,ii,i,j)=x(indi)
  11     continue
        endif
  10  continue
   1  continue
      return 
      end
c
      subroutine getl(ig,maxg,nc,maxc,maxv,maxe,wmat1,wmat2,cov,wmis,
     & wvec1,wvec2,p,wcov,obs,xmis,misindi,
     & npg,npc,nv,maxn,mvt,np,xm,objf,ifail,mess)
c calculate loglikelihood
      implicit double precision (a-h,o-z), integer (i-n)
      double precision cov(maxg,maxc,maxv,maxv),wcov(maxc,maxv,maxv)
      double precision wmat1(maxv,maxv),wvec2(maxv),p(maxg,maxc)
      double precision wmat2(maxv,maxv),wvec1(maxv)
      double precision obs(maxn,mvt),xm(maxg,maxc,maxv)
      double precision wmis(maxv,maxv)
      double precision det,f,tmp,tmp1
      integer npg(maxg),npc(0:maxg),nc(maxg),misindi(maxn)
      character mess*5,name(6)*2
      common / bits  / pi2,half,shalf,one,zero,two
      common / names / name
c
      do 30 ic=1,nc(ig)
c
      do 32 i=1,nv
      do 32 j=1,nv
  32  wmat1(i,j)=cov(ig,ic,i,j)
      ifail=0
c invert covariance matrix
      call CHOLINV(wmat1,maxv,nv,DET,IFAIL,wvec2)
c next can very well occur
      if (ifail.ne.0) then 
       mess='SSing'
       return
      endif
      do 33 i=1,nv
      do 33 j=1,nv
 33   wcov(ic,i,j)=wmat1(i,j)
      tmp=one/dsqrt(det)
      tmp1=pi2**(-dble(nv)/two)
c
      do 200 icase=npc(ig-1)+1,npc(ig)
c
      if (misindi(icase).eq.0) then
c
      do 43 iv=1,nv
 43   wvec2(iv)=obs(icase,iv)-xm(ig,ic,iv)
      call mvn(wvec2,wmat1,nv,det,tmp,f,maxv,tmp1)
      obs(icase,nv+ic)=f
      obs(icase,nv+nc(ig)+1)=obs(icase,nv+nc(ig)+1)+p(ig,ic)*f
      endif
 200  continue
c
      do 201 icase=npc(ig-1)+1,npc(ig)
c
      if (misindi(icase).gt.0) then
c
      do 119 iv=1,nv
      do 119 jv=1,nv
 119  wmis(iv,jv)=0.d0
      nvobs=0
      do 120 iv=1,nv
      if (obs(icase,iv).ne.xmis) then      
      nvobs=nvobs+1
      wvec2(nvobs)=obs(icase,iv)-xm(ig,ic,iv)
      wmis(nvobs,iv)=1.d0
      endif
 120  continue
c
      do 2 i=1,nvobs
      do 2 j=1,nv
      wmat2(i,j)=0.d0 
      do 2 k=1,nv
      wmat2(i,j)=wmat2(i,j)+wmis(i,k)*cov(ig,ic,k,j)
 2    continue
      do 3 i=1,nvobs
      do 3 j=1,nvobs
      wmat1(i,j)=0.d0
      do 3 k=1,nv
      wmat1(i,j)=wmat1(i,j)+wmat2(i,k)*wmis(j,k)
 3    continue
      ifail=0
c invert covariance matrix
      call CHOLINV(wmat1,maxv,nvobs,DET,IFAIL,wvec1)
c next can very well occur
      if (ifail.ne.0) then 
       mess='SSing'
       return
      endif
      tmp=one/dsqrt(det)
      tmp1=pi2**(-dble(nvobs)/two)
      call mvn(wvec2,wmat1,nvobs,det,tmp,f,maxv,tmp1)
      obs(icase,nv+ic)=f
      obs(icase,nv+nc(ig)+1)=obs(icase,nv+nc(ig)+1)+p(ig,ic)*f
c
      endif
 201  continue
 30   continue
c
      do 205 i=npc(ig-1)+1,npc(ig)
c
      if (obs(i,nv+nc(ig)+1).le.zero) then 
       ifail=1
       mess='obs=0'
       return
      endif
c
      objf=objf+dlog(obs(i,nv+nc(ig)+1))
 205  continue
c     objf=-objf
      return 
      end
c
      subroutine getdd2(ig,ic,is,maxn,mvt,maxv,obs,xmis,
     &      wmis,wvec1,wvec2,wvec3,wcov,xm,wmat1,wmat2,cov,
     &      dxm,dcov,p,nv,nvobs,np,nc,npc,maxc,maxg)
c calculate dL\d(SIGMA) and dL=d(MEAN) where L is loglikelihood 
      implicit double precision (a-h,o-z), integer (i-n)
      double precision wcov(maxc,maxv,maxv),dcov(maxv,maxv)
      double precision wmat1(maxv,maxv),wvec1(maxv),p(maxg,maxc)
      double precision obs(maxn,mvt),xm(maxg,maxc,maxv),dxm(maxv)
      double precision xtmp,wmis(maxv,maxv),wmat2(maxv,maxv)
      double precision wvec3(maxv),xmis
      double precision cov(maxg,maxc,maxv,maxv)
      integer npc(0:maxg),nc(maxg)
      integer i,j,iv,jv,ic
      character name(6)*2
      common / bits  / pi2,half,shalf,one,zero,two
      common / names / name
c     do 444 i=1,nv
c     dxm(i)=zero
c     do 444 j=1,nv
c444  dcov(i,j)=zero
c
c     do 220 is=npc(ig-1)+1,npc(ig)
c     if (misindi(is).eq.1) goto 220
      xtmp=(p(ig,ic)*obs(is,nv+ic))/obs(is,nv+nc(ig)+1)
      do 119 iv=1,nv
      do 119 jv=1,nv
 119  wmis(iv,jv)=0.d0
      nvobs=0
      do 120 iv=1,nv
      if (obs(is,iv).ne.xmis) then      
      nvobs=nvobs+1
      wvec3(nvobs)=obs(is,iv)-xm(ig,ic,iv)
      wmis(nvobs,iv)=1.d0
      endif
 120  continue
      do 2 i=1,nvobs
      do 2 j=1,nv
      wmat1(i,j)=0.d0 
      do 2 k=1,nv
      wmat1(i,j)=wmat1(i,j)+wmis(i,k)*cov(ig,ic,k,j)
 2    continue
      do 3 i=1,nvobs
      do 3 j=1,nvobs
      wmat2(i,j)=0.d0
      do 3 k=1,nv
      wmat2(i,j)=wmat2(i,j)+wmat1(i,k)*wmis(j,k)
 3    continue
      ifail=0
c invert covariance matrix
      call CHOLINV(wmat2,maxv,nvobs,DET,IFAIL,wvec2)
      do 54 iv=1,nvobs
      wvec1(iv)=zero
      do 54 jv=1,nvobs
      wvec1(iv)=wvec1(iv)+wmat2(iv,jv)*wvec3(jv)
 54   continue
      do 55 iv=1,nvobs
 55   dxm(iv)=xtmp*wvec1(iv)
c **************************  cov dL/d(cov)
      do 60 iv=1,nvobs
      do 60 jv=1,nvobs
 60   wmat1(iv,jv)=wvec1(iv)*wvec1(jv)
      do 62 iv=1,nvobs
      do 62 jv=1,nvobs
 62   wmat1(iv,jv)=(wmat1(iv,jv)-wmat2(iv,jv))*xtmp*half 
      do 221 i=1,nvobs
      do 221 j=1,nvobs
 221  dcov(i,j)=wmat1(i,j)
c
c220  continue
      return
      end
      subroutine getdd1(ig,ic,is,maxn,mvt,maxv,obs,
     &                  wvec1,wcov,xm,wmat1,
     &                  dxm,dcov,p,nv,np,nc,npc,maxc,maxg)
c calculate dL\d(SIGMA) and dL=d(MEAN) where L is loglikelihood 
      implicit double precision (a-h,o-z), integer (i-n)
      double precision wcov(maxc,maxv,maxv),dcov(maxv,maxv)
      double precision wmat1(maxv,maxv),wvec1(maxv),p(maxg,maxc)
      double precision obs(maxn,mvt),xm(maxg,maxc,maxv),dxm(maxv)
      double precision xtmp
      integer npc(0:maxg),nc(maxg)
      integer i,j,iv,jv,ic
      character name(6)*2
      common / bits  / pi2,half,shalf,one,zero,two
      common / names / name
c     do 444 i=1,nv
c     dxm(i)=zero
c     do 444 j=1,nv
c444  dcov(i,j)=zero
c
c     do 220 is=npc(ig-1)+1,npc(ig)
c     if (misindi(is).eq.1) goto 220
      xtmp=(p(ig,ic)*obs(is,nv+ic))/obs(is,nv+nc(ig)+1)
      do 54 iv=1,nv
      wvec1(iv)=zero
      do 54 jv=1,nv
 54   wvec1(iv)=wvec1(iv)+wcov(ic,iv,jv)*(obs(is,jv)-xm(ig,ic,jv))
      do 55 iv=1,nv
 55   dxm(iv)=dxm(iv)+xtmp*wvec1(iv)
c **************************  cov dL/d(cov)
      do 60 iv=1,nv
      do 60 jv=1,nv
 60   wmat1(iv,jv)=wvec1(iv)*wvec1(jv)
      do 62 iv=1,nv
      do 62 jv=1,nv
 62   wmat1(iv,jv)=(wmat1(iv,jv)-wcov(ic,iv,jv))*xtmp*half 
      do 221 i=1,nv
      do 221 j=1,nv
 221  dcov(i,j)=dcov(i,j)+wmat1(i,j)
c220  continue
      return
      end
      subroutine 
     &getpp(nv,npg,p,obs,xmis,nc,maxn,mvt,maxc,iunit,ent,entr,bu,maxg,
     &ng,npc,ig)
      implicit double precision (a-h,o-z),integer (i-n)
c calculates entropy measure and apost assigment to categories
      double precision obs(maxn,mvt),p(maxg,maxc),tmp,entr,ent
      double precision bu(maxc)
      integer ic,i,np,nc(maxg),j,iunit,ind,npc(0:maxg),npg(maxg)
      ent=0.d0
      entr=0.d0
      do 400 i=npc(ig-1)+1,npc(ig)
      ic=0
      tmp=0.d0
      do 401 j=1,nc(ig)
      bu(j)=(p(ig,j)*obs(i,nv+j))/obs(i,nv+nc(ig)+1)
c
       if (bu(j).gt.0.d0) then
        entr=entr-bu(j)*dlog(bu(j))
c      else
c       ind=ind+1
       endif
         if (bu(j).gt.tmp) then
          tmp=bu(j)
          ic=j
         endif
 401  continue
c unit 4 = postp.dat
      write(iunit,912) ig,i,ic,(bu(j),j=1,nc(ig))
 400  continue
      ent=entr
      entr=1.d0-(entr/(dble(npg(ig))*dble(nc(ig))))
c     if (ind.eq.2) then
       write(iunit,913)  ig,ent,entr       
c     else
c      write(iunit,*) ' postp. probs. zero: no entropy '
c      entr=dble(ind)
c     endif
 913  format(' group ',i2,' en(tau) ',f8.3,' entropy ',f7.3)
 912  format(' group ',i2,' case ',i4,' component ',i2,1x,6(f6.3,1x))
      return 
      end
      subroutine reader(ig,inf2,np,nv,nc,ne,ng,iw1,maxqn,ipat,xmis,maxg)
c reads "user friendly" input 
      implicit double precision (a-h,o-z),integer (i-n)
c
      parameter (nkeyw=15)
      character xl1*60,xl2*60,xl3*60,bl3*3,inf2*12
      character txt(nkeyw)*3,bit*79,nn*13,bl*1,xll*210
      integer ipat(maxg,6),np,nv,ne,nc,iw1,maxqn
      double precision xmis
c
      data txt/ 'FI=','NP=','NV=','NC=','NE=','QC=',
     &'IT=','LY=','BE=','PS=','TE=','NU=','AL=','NG=','MI='/
      data bl3,bl/ '   ',' '/
c
c ME IT=500 QC=10 
c DA FI=LY10 NP=400 NG=1
c MO NV=5 NC=2 NE=1 LY=2 BE=0 PS=1 NU=2 TE=2 AL=0
c
 99   format(a79)
      open(7,file='scratch')
      rewind 7
c read title
 555  read(1,'(a60)') bit
ce    write(16,'(a60)') bit
      if (bit(1:6).eq.bl3//bl3) goto 555
      if (bit(1:2).ne.'TI') stop ' line 1 should start `TI` '
c read three next lines (DA, ME, MO) strip DA ME and MO
c safe assumption: keywords occupy no more that 70 spaces 
c SAFE assumption: space between keywords no more than 9 spaces 
      read(1,99) bit
      i1=index(bit(1:79),bl3//bl3//bl3)-1
      xl1(1:i1-3)=bit(4:i1)
ce    write(16,99) bit
      read(1,99) bit
      i2=index(bit(1:79),bl3//bl3//bl3)-1
      xl2(1:i2-3)=bit(4:i2)
ce    write(16,99) bit
      read(1,99) bit
      i3=index(bit(1:79),bl3//bl3//bl3)-1
      xl3(1:i3-3)=bit(4:i3)
ce    write(16,99) bit
c glue together lines 
      il=i1+i2+i3-6
      xll(1:il)=xl1(1:i1-3)//bl//xl2(1:i2-3)//bl//xl3(1:i3-3)//bl

c seek keywords and, of these, the assigned values
c
      do 100 ii=1,nkeyw
      i=index(xll(1:il),txt(ii))
      if (i.eq.0) then 
       write(16,999) txt(ii)
       stop ' fatal see output  '  
      endif
 999  format(' keyword ',a3,' missing or +9 spaces between keywords!')
      j=index(xll(i:il),bl)
      nn=xll(i+3:j+i-1)
c write info as characters
      write(7,'(a13)') nn
 100  continue
ce    WRITE(6,*) 
      close(7,status='keep')
      open(7,file='scratch',status='old')
      rewind 7
c and read info in appropriate type
      read(7,'(a12)',err=500) inf2
      read(7,*,err=502) np
      read(7,*,err=503) nv
      read(7,*,err=504) nc
      read(7,*,err=505) ne
      read(7,*,err=508) iw1
      read(7,*,err=509) maxqn
      read(7,*,err=511) ipat(ig,1)
      read(7,*,err=512) ipat(ig,2)
      read(7,*,err=513) ipat(ig,3)
      read(7,*,err=514) ipat(ig,4)
      read(7,*,err=515) ipat(ig,5)
      read(7,*,err=516) ipat(ig,6)
      read(7,*,err=518) ng         
      read(7,*,err=519) xmis
c     close(7,status='delete')
      return
c error messages 
 500  stop ' error in reading data file name '
 502  stop ' error in NP= (number of cases)'
 503  stop ' error in NV= (number of variables)'
 504  stop ' error in NC= (number of components)'
 505  stop ' error in NE= (number of eta variables)'
 508  stop ' error in QC= (convergence crit QN)'
 509  stop ' error in MQ= (max iterations QN)'
 511  stop ' error in LY= '
 512  stop ' error in BE= '
 513  stop ' error in PS= '
 514  stop ' error in TE= '
 515  stop ' error in NU= '
 516  stop ' error in AL= '
 518  stop ' error in NG= '
 519  stop ' error in MI= '
      end
c
      subroutine repdet(wmat,maxv,m,wvec)
c claculates and reports determinant of symm. matrix
      implicit double precision (a-h,o-z),integer (i-n)
      double precision det,wmat(maxv,maxv),wvec(maxv)
      integer m,ifail
      external cholesk
      ifail=0
      call CHOLESK(m,wmat,maxv,wvec,DET,IFAIL)
      if (ifail.ne.0) then
      write(16,*) ' -> NOT POSITIVE DEFINITE ! '
      else
      write(16,999) det
      endif
 999  format(' DETERMINANT: ',E12.5)
      return
      end
c
      subroutine teltime(txt)
c the time
      character txt*10
      integer i1,i2,i3,i4
      call gettim(i1,i2,i3,i4)
      write(16,917) txt,i1,i2,i3,i4
 917  format
     & (53x,a10,' at ',i2,':',i2,':',i2,':',i2)
      return
      end
c
      subroutine check1(maxqn,iwk,ne,maxe,nc,
     &       maxn,np,maxc,nv,maxv,ipat,name,ng,maxg)
c checks on input
      implicit integer (i-n), double precision (a-h,o-z)
      character name(6)*2
      integer ipat(maxg,6),iwk,maxqn,ne,maxe,nc,maxc
      integer nv,maxv
      if (maxqn.le.0) stop ' strange value IT= '
      if (iwk.le.0)   stop ' strange value QC= ' 
      if (ne.gt.maxe) stop ' ne .gt. maxe (setpar.inc) '
      if (nc.gt.maxc) stop ' nc .gt. maxc (setpar.inc) '
      if (np.gt.maxn) stop ' np .gt. maxn (setpar.inc) '
      if (nv.gt.maxv) stop ' nv .gt. maxv (setpar.inc) '
      if (ng.gt.maxg) stop ' ng .gt. maxg (setpar.inc) '
      if (maxe.gt.maxv) stop ' maxe.gt.maxv...gt not allowed (do eq!) '
      do 42 j=1,ng
      do 42 i=1,6
       if (ipat(j,i).lt.0.or.ipat(j,i).gt.2) then 
       write(16,916) name(i),ipat(j,i)
       stop ' error see output '
      endif
  42  continue
 916  format(' parameter matrix ',a2,' cannot be assigned ',i1,/,
     &' should be 0 (not used) 1 (fixed) or 2 (free) ')
      return
      end