      program etaspa
c-----------------------------------------------------------------------
c  space-time version of the etas model.
c-----------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(npa=22222,nda=22222,mdl=43333,nr=020)
      parameter(npall=5*npa)
      common/data/ xx(nda),yy(nda),zz(nda),zmg(nda),
     &            nn,mm,tx,ty,tz,ind(nda)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /swgrad/ w1,w2,w4,w5,w7,ptdet,w00,w01
      common /prihe5/ adp1(npa),aup1(nr,npa),adp2(npa),aup2(nr,npa),
     & adp4(npa),aup4(nr,npa),adp5(npa),aup5(nr,npa),
     & adp7(npa),aup7(nr,npa)
      common /ccc/ isw,ipr
      common /xad/xd(npall),m
      common /prupd/ iupd
      common /abicc/abic,fmin
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      common /fixpar/ a1,a2,a3,aa4,a5,a6,a7
      common /hesapprox/nhesapp
      common /simpdiep/dist,eps
      dimension x(9)
      dimension xdd(npall)
      external func17,funcsb
c     data dist,eps /.1d1, 1.d-3/
c     include'mpif.h'
c
c     time1=dclock()
      myrank=0
c     nhesapp=0
c     nhesapp=1
      iupd=0
cyo
c     w00=0.1d4
      w00=0.0d0
c     w01=1.d0/0.25
c
c     w01=0.1d-1
c      call mpi_init(ierr)
c      call mpi_comm_size(mpi_comm_world,nprocs,ierr)
c      call mpi_comm_rank(mpi_comm_world,myrank,ierr)
c      write(*,*)myrank, ' in ', nprocs
c
      call redata(xx,yy,zz,zmg,ind,tx,ty,tz,nn)
cbsetas2d-mpi-modified.f
c     call para(nn)

      call parainit
c
cf. funcsb
      x(1)=log(w1)
      x(2)=log(w2)
      x(3)=log(a3)
c     x(4)=log(aa4/2)
      x(4)=log(aa4)
      x(5)=log(a5)
c     x(5)=log(a5-1)
      x(6)=log(a6)
      x(7)=log(a7)
c     x(7)=log(a7-1)


    1 format(1h ,5f10.3/(1h ,5f10.3))
c  Murata's routine follows.
c
      call simplx(x, 7, funcsb, dist, eps, xm)
c     call simplx(x, 2, funcsb, dist, eps, xm)
c
      w1=exp(x(1))
      w2=exp(x(2))
      a3=exp(x(3))
      aa4=exp(x(4))
      a5=exp(x(5))
c     a5=exp(x(5))+1
      a6=exp(x(6))
      a7=exp(x(7))
c     a7=exp(x(7))+1


      if(myrank.eq.0) then
      open(11,file='hist-etas-mk.upda')
c     write(11,5) w1,w2,abic,m
      write(11,5) w1,w2,fmin,m
    5 format(2e15.5,f11.2,i7)
c      write(11,*) (x(i),i=1,2),abic,npext
      write(11,2) a1,a2,a3,aa4,a5,a6,a7
       write(11,2) (xd(i),i=1,m)
c   2 format(5e15.5)
    2 format(4e20.12)
      close(11)
      do 210 i=1,m/2
        xdd(i)=log(a1)+xd(i)
        xdd(i+m/2)=a2*exp(xd(i+m/2))
  210 continue
      open(21,file='hist-etas-mk.omap')
c     open(21,file='delopara7pa.fig')
      write(21,5) (x(i),i=1,2),abic,m
      write(21,2) a1,a2,a3,aa4,a5,a6,a7
      write(21,2) (xdd(i),i=1,m)
      close(21)
      end if
c     call mpi_finalize(ierr)
c     time2=dclock()
c     write(6,*) 'time',time,'sec'
c     write(6,*) 'time',time2-time1,'sec'
      stop
      end
c***********************************************************************
      subroutine redata(xx,yy,zz,zmg,ind,tx,ty,tz,nn)
      implicit real * 8 (a-h,o-z)
      parameter(npa=22222,nda=22222,mdl=43333,nr=020,nmax=6000)
      parameter(npall=npa*5)
      common /coef/ aa(nda),bb(nda),cc(nda)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /delink/ nzu(npa),iu(nr,npa)
      common /xad/xd(npall),m
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      common /jfun/rmax(nda),jxy(nda,nmax),jnn(nda)
      common /range/tstart,ntstar
      dimension xx(1),yy(1),zz(1),zmg(1), ind(1)
      character*8 fmt(10)
      character*80 maindata
      open(1,file='hist-etas-mk.conf')
      read(1,112) maindata
  112 format(a)
!     maindata='delone2.out'
c     read(1,1) nfunct
      nfunct=17
    1 format(i10)
c     read(1,2) fmt
    2 format(10a8)
      read(1,*) tx,ty,tz,nn
    3 format(3f10.1,i10)
      read(1,*) xmin,ymin,xmg0,zmin,tstart,bi2
c     if(myrank.eq.0) then
      write(6,112) maindata
c     write(6,2) fmt
      write(6,3) tx,ty,tz,nn,tstart
      write(6,8) xmin,ymin,xmg0,zmin,tstart,bi2
!     write(9,112) maindata
c     write(9,2) fmt
!     write(9,3) tx,ty,tz,nn
!     write(9,8) xmin,ymin,xmg0,zmin,tstart,bi2
c     end if
    8 format(8f10.2)
      in=10
      open(in,file=maindata)
      read(in,*) np,npext,ndelo,tx,ty
    7 format(1h ,'data set  ',10a8)
      m=npext*2
      nnc=0
c
*     open(13,file='../delo2d/alljp95m50.config')
*     read(13,*) 
c
      do 20 i=1,npext
c     read(in,fmt,end=25) xx(i),yy(i),zmg(i),zz(i),dd,ee,ff
      read(in,*,end=25) idum,xx(i),yy(i),zmg(i),zz(i),dd,ee,ff
c     if(myrank.eq.0) write(6,*) i,xx(i),yy(i),zmg(i),zz(i),dd,ee,ff
*     read(13,*) idummy,xx(i),yy(i)
      if(i.gt.np) go to 20
      aa(i)=ee/dd/sqrt(1-ff**2)
      bb(i)=ff   /sqrt(1-ff**2)
      cc(i)=dd/ee/sqrt(1-ff**2)
c     write(6,*) i,myrank,xx(i)
   11 format(i5,4f8.4)
   20 continue
      i=nn+1
   25 nn=i-1
c     write(6,*) 'nn',nn,xx(nn),yy(nn),zz(nn)
c     read(in,*) (idummy,(ldt(i,j),j=1,3),adt(i), i=1,ndelo)
      do 509 i=1,ndelo
      read(in,*) idummy,(ldt(i,j),j=1,3),adt(i)
*     read(13,*) idummy,(ldt(i,j),j=1,3),adt(i)
  509 continue
      do 520 i = 1, npext
        read(in,*) idummy, nzu(i), (iu(j,i), j=1,nzu(i))
*       read(13,*) idummy, nzu(i), (iu(j,i), j=1,nzu(i))
  520 continue
      ii=0
      do 10 i=1,nn
*     if(xx(i).gt.tx+xmin) go to 10
*     if(xx(i).lt.   xmin) go to 10
*     if(yy(i).gt.ty+ymin) go to 10
*     if(yy(i).lt.   ymin) go to 10
*     if(zz(i).gt.tz+zmin) go to 10
*     if(zz(i).lt.   zmin) go to 10
      ind(i)=1
      if(zz(i).gt.tz+zmin) ind(i)=0
      if(zz(i).lt.   zmin) ind(i)=0
      ii=ii+1
!     write(6,*) 'chk',ii,zz(i),ntstar
      if(zz(i).lt.tstart) ntstar=ii
c  2007/06/22
c     xx(ii)=xx(i)-xmin
c     yy(ii)=yy(i)-ymin
*     zz(ii)=zz(i)-zmin
      xx(ii)=xx(i)
      yy(ii)=yy(i)
      zz(ii)=zz(i)
      zmg(ii)=zmg(i)-xmg0
       aa(ii)=aa(i)
       bb(ii)=bb(i)
       cc(ii)=cc(i)
c      if(aa(ii).ne.1.d0) write(6,*) ii,aa(ii),bb(ii),cc(ii)
       if(zz(ii).gt.tstart) nnc=nnc+1
       xx(ii)=xx(ii)*cos((ymin+ty/2)*3.14159/180.)
c      write(6,*) ii,xx(ii),ymin,ty/2,cos((ymin+ty/2)*3.14159/180.)
   10 continue
      tx=tx*cos((ymin+ty/2)*3.14159/180.)
      do 40 i=1,ndelo
        adt(i)=adt(i)*cos((ymin+ty/2)*3.14159/180.)
   40 continue
c
      nn=ii
      if(myrank.eq.0) then
      write(6,4) in
!     write(9,4) in
    4 format(1h ,'input device',i10)
      write(6,*) 'nfunct=',nfunct
      write(6,*) 'nn,ntstar,nnc',ii,ntstar,nnc
!     write(9,*) 'nn,ntstar,nnc',ii,ntstar,nnc
      write(6,5) tx,ty,tz,xmin,ymin,xmg0,zmin,tstart,nn,nnc
!     write(9,*) 'nfunct=',nfunct
!     write(9,5) tx,ty,tz,xmin,ymin,xmg0,zmin,tstart,nn,nnc
    5 format(1h ,'tx,ty,tz,xmin,ymin,xmg0,zmin,tsta'/8f10.3,
     &/' nn =',i5,' nnc =',i5)
c     write(6,6) ((i,xx(i),yy(i),zmg(i),zz(i)),i=1,nn)
      end if
    6 format(///1h ,i5,4f13.5/(1h ,i5,4f13.5))
c
      do i=1,nn
      do j=1,nmax
      jxy(i,j)=0
      end do
      jnn(i)=0
      end do
      jmax=0
c j-select-redata
      do 500 i=2,nn
      do 510 j=1,i-1
      xxx=max(10**(0.5*(zmg(j)+xmg0)-1),100.d0)*bi2
c     if(zz(i)-zz(j).gt.xxx) go to 510
      rmax(j)=(0.015*10.**(0.5*(zmg(j)+xmg0)-2.)+0.3)*bi2
      xxr=abs(xx(i)-xx(j))
      yyr=abs(yy(i)-yy(j))
      if(xxr**2+yyr**2 .gt. rmax(j)**2) go to 510
      jnn(i)=jnn(i)+1
      ji=jnn(i)
      if(ji.gt.nmax) then
      write(6,*) 'nmax over',ji
      stop 119
      end if
      jxy(i,ji)=j
c     go to 510
c 512 write(6,*) i,j,rmax(j)**2,xxr**2+yyr**2,xxr,yyr
  510 continue
      ji=jnn(i)
      if(jmax.lt.ji) jmax=ji
c     do j=1,ji
c     ij=jxy(i,j)
c     write(6,*) i,ji,ij,zmg(ij)+xmg0
c     end do
  511 format(2i6,10f5.1/(12x,10f5.1))
  500 continue
      write(6,*) 'jmax',jmax
!     write(9,*) 'jmax',jmax
c     stop 'to show cluster'
      return
      end
c***********************************************************************
      subroutine parainit
      implicit real * 8 (a-h,o-z)
      parameter(npa=22222,nda=22222,mdl=43333,nr=020)
c 00/08/29
      parameter(npall=5*npa)
      common /swgrad/ w1,w2,w4,w5,w7,ptdet,w00,w01
      common /fixpar/ a1,a2,a3,aa4,a5,a6,a7
c 00/08/29
c     common /xad/xd(npa),m
      common /xad/xd(npall),m
      common /data/ xx(nda),yy(nda),zz(nda),zmg(nda),
     &              nn,mm,tx,ty,tz,ind(nda)
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      common /abicc/abic,fmin
      common /initsi/inits,initf
      common /simpdiep/dist,eps
      common /hesapprox/nhesapp
      common     / ccc /  isw , ipr
      common /davtau/tau1,tau2
      common /daveps/eps1,eps2
      external func17
      dimension x(50)
      character*60 delocoeff

c     include 'mpif.h'

c Choice of initial coefficients for the delone system
      fmin=1.d60
continuation
c      init=0
c     init=1
      read(1,*) init
      read(1,*) inits
      read(1,*) initf
      if(init.eq.0) then
c       w1=1.d0
c       w2=1.d0
c

c        Set up initial values for the parameter surface
        do 10 i=1,m
c         xd(i)=i*1.0d0/m
          xd(i)=i*0.0d0/m
   10   continue


        read(1,112) delocoeff
        read(1,*) w00,w1,w2
  112   format(a)
      else
        read(1,112) delocoeff
        read(1,*) w00,w1,w2
        open(2, file=delocoeff)
        read(2,*) w1,w2,fmin,m
        read(2,*) a1,a2,a3,a4,a5,a6,a7
        read(2,*) (xd(i),i=1,m)
*       read(2,9) (xd(i),i=1,npext)
        close(2)
      fmin=1.d60
    9   format(4f20.0)
c       if(m.ne.npext) then
c         write(6,*) 'two data sets are inconsitent.'
c         stop 44444
c       endif
      endif
c
      read(1,1000) n
      read(1,*) (x(i),i=1,n)
      write(6,*) 'n=',n
      write(9,*) 'n=',n
      write(6,1030)  (x(i),i=1,n)
      read(1,*) ipr
      write(6,*) 'linear ipr',ipr
c     read(1,*) tau1,tau2
c     write(6,*) 'davidn tau',tau1,tau2
c     read(1,*) eps1,eps2
c     write(6,*) 'davidn eps',eps1,eps2
      read(1,*) nhesapp,dist,eps
      write(6,*) 'nhesapp,simplex(dist,eps)',nhesapp,dist,eps
      if(myrank.eq.0) then
      write(6,*) 'n=',n
!     write(9,*) 'n=',n
c     write(6,1030)  (x(i),i=1,n)
      end if
      if(init.eq.0) then
        a1=x(1)
        a2=x(2)
        a3=x(3)
        a4=x(4)
        a5=x(5)
        a6=x(6)
        a7=x(7)
      endif
c     aa4=a4*2
      aa4=a4

      return
 1000 format(3i10,2f15.6)
 1010 format(8f10.4)
 1030 format(1h ,'para-init=',7e14.5)
      end
c***********************************************************************

      subroutine func17(n,b,f,h,ifg)
c-----------------------------------------------------------------------
c     likelihood function of modified oomori type
c     etas point process (with trend: not realized)
c     the space-time <<anisotropic>> version.
c \lambda(t,x,y)=\mu+\sum k/(t-t_j+c)^p
c                /( (x-x_j)^2+(y-y_j)^2+d )^q/\exp(\alpham_j)^q
c---------------------------------------------------------------------
      implicit real * 8 (a-h,o-z)
      parameter(npa=22222,nda=22222,mdl=43333,nr=020,nmax=6000)
      parameter(npall=5*npa)
      common /data/ xx(nda),yy(nda),zz(nda),zmg(nda),
     &             nn,mm,tx,ty,tz,ind(nda)
      common /coef/ aa(nda),bb(nda),cc(nda)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /delink/ nzu(npa),iu(nr,npa)
      common /fixpar/ a1,a2,a3,aa4,a5,a6,a7
c     common /fixpar/ a1,a2,a3,a4,a5,a6,a7
      common /swgrad/ w1,w2,w4,w5,w7,ptdet,w00,w01
      common /parts / rmd(nda),bx(4096),by(4096),ft
      common /sublin/ilin
      common /ellips/qrrd1(nda,128,4),qrrd2(nda,128,4),qrr1(nda,128,4),
     &               qrr2(nda,128,4),qrr3(4),mbx,mby
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      common /jfun/rmax(nda),jxy(nda,nmax),jnn(nda)
      common /range/tstart,ntstar
      data pi /3.14159265358979d0/
      dimension b(npall),g(npall),h(npall)
      dimension b1(npa),d1ff(npa),d1ft(npa),g1(npa),gp1(npa)
      dimension b2(npa),d2ff(npa),d2ft(npa),g2(npa),gp2(npa)
      dimension rmd2i(nda),rmd4i(nda),rmd5i(nda),rmd7i(nda)
      dimension bp1(npa),bp2(npa),bp4(npa),bp5(npa),bp7(npa)
      data index /0/
      save index
      q(ai,bi,ci,x,y)=ai*x**2-2*bi*x*y+ci*y**2
c 
      a4=aa4/2      
c     write(6,*) 'a1-7',a1,a2,a3,a4,a5,a6,a7
      if(index.eq.0) then
        mbx=2**3
        mby=2**3
        do 30 i=1,mbx+1
          bx(i)=(i-1)*tx/mbx
   30   continue
        do 40 i=1,mby+1
          by(i)=(i-1)*ty/mby
   40   continue
      end if
c
      sum=0
      do 70 i=1,n
        if(b(i).gt.50.) b(i)=50d0
      sum=sum+b(i)
   70 continue
      b1min=1.0d99
      b2min=1.0d99
      b1max=0.0
      b2max=0.0
      do 80 i=1,n/2
        b1(i)=a1*exp(b(i))
        b2(i)=a2*exp(b(i+n/2))
        d1ff(i)=0.0
        d2ff(i)=0.0
        if(b1(i).lt.b1min) b1min=b1(i)
        if(b2(i).lt.b2min) b2min=b2(i)
        if(b1(i).gt.b1max) b1max=b1(i)
        if(b2(i).gt.b2max) b2max=b2(i)
   80 continue
      ifg=0
      ff=0.0
      ramdai=b1(1)
      if(ramdai.gt.0.0) then
         ff=log(ramdai)*ind(1)
         d1ff(1)=1.0d0/ramdai*ind(1)
      endif
      rmd(1)=0.0
      if(ilin.eq.1) then
        ftemp=0
chistory
        do 90 i=1,nn
          if(ind(i).eq.1) then
            rmdi=0.0
            ji=jnn(i)
            do 100 j1=1,ji
              j=jxy(i,j1)
              rmd2i(j)=0.0
              xxij=(xx(i)-xx(j))*exp(-a4*zmg(j))
              yyij=(yy(i)-yy(j))*exp(-a4*zmg(j))
              space=(aa(j)*xxij**2-2*bb(j)*xxij*yyij+
     &              cc(j)*yyij**2+a6)**(-a7)
              rmdi=rmdi+b2(j)/(zz(i)-zz(j)+a3)**a5 *space
  100       continue
chistory
            if(zz(i).lt.tstart) go to 90
c
            ramdai=b1(i)+rmdi
            if(ramdai.le.0.0) go to 50
            if(i.gt.1)ff = ff +log(ramdai) 
          endif
   90   continue
      else
chistory
        do 10 i=2,nn
          if(ind(i).eq.1)then
            rmdi=0.0
            ji=jnn(i)
            do 20 j1=1,ji
              j=jxy(i,j1)
              rmd2i(j)=0.0
              xxij=(xx(i)-xx(j))*exp(-a4*zmg(j))
              yyij=(yy(i)-yy(j))*exp(-a4*zmg(j))
              space=(aa(j)*xxij**2-2*bb(j)*xxij*yyij+
     &          cc(j)*yyij**2+a6)**(-a7)
c       if(i.ge.3350) write(6,*) 'ee',i,j,rmdi
c       if(i.ge.3350) write(6,*) 'ef',i,j,rmdi
c    &                               ,b2(j)/(zz(i)-zz(j)+a3)**a5*space
c       if(i.ge.3350) write(6,*) 'eg',i,j
c    &                               ,b2(j),(zz(i)-zz(j)+a3)**a5,space
c       if(i.ge.3350) write(6,*) 'ehhhhh',i,j
c    &       ,aa(j),xxij**2,2*bb(j),xxij*yyij,cc(j),yyij**2,a6
              rmdi=rmdi+b2(j)/(zz(i)-zz(j)+a3)**a5 *space
ccccccccccccccccccccccccc
c       if(i.ge.3350.and.j.ge.2112) write(6,*) 'ff',i,j,rmdi
c       if(i.ge.3350.and.j.ge.2112) write(6,*)'g',b2(j)/(zz(i)-zz(j)+a3)**a5*space
c       if(i.ge.3350.and.j.ge.2112) stop 3351
              rmd2i(j)=rmd2i(j)+1/(zz(i)-zz(j)+a3)**a5 *space
   20       continue
chistory
            if(zz(i).lt.tstart) go to 10
c
            ramdai=b1(i)+rmdi
            if(ramdai.le.0.0) write(6,*) b1(i), rmdi
            ff=ff+log(ramdai)     
c       write(6,*) 'ff',i,ff
            d1ff(i)=d1ff(i) + 1.0d0/ramdai
            ji=jnn(i)
            do 60 j1=1,ji
              j=jxy(i,j1)
              d2ff(j)=d2ff(j)+rmd2i(j)/ramdai
   60       continue
          endif
   10   continue
      end if
c
      do 110 i=1,n/2
        d1ft(i)=0.
        d2ft(i)=0.0
  110 continue
      ft=0.0
      if(index.eq.0) then
c
        do 120 i=1,nn
          do 130 j=1,mby
            if(xx(i).le.0.0) go to 130
            jy=j
            k=1
            qrrd1(i,j,k)=q(aa(i),bb(i),cc(i), xx(i),by(jy)-yy(i))
            qrrd2(i,j,k)=q(aa(i),bb(i),cc(i), xx(i),by(jy+1)-yy(i))
            qrr1(i,j,k)=q(1.0d0,0.0d0,1.0d0, xx(i),by(jy)-yy(i))
            qrr2(i,j,k)=q(1.0d0,0.0d0,1.0d0, xx(i),by(jy+1)-yy(i))
            qrr3(k)=(ty/mby)**2
  130     continue
          do 140 j=mby+1,mbx+mby
            if(yy(i).ge.ty) go to 140
            jx=j-mby
            k=2
            qrrd1(i,j,k)=q(aa(i),bb(i),cc(i), bx(jx)-xx(i),ty-yy(i))
            qrrd2(i,j,k)=q(aa(i),bb(i),cc(i), bx(jx+1)-xx(i),ty-yy(i))
            qrr1(i,j,k)=q(1.0d0,0.0d0,1.0d0, bx(jx)-xx(i),ty-yy(i))
            qrr2(i,j,k)=q(1.0d0,0.0d0,1.0d0, bx(jx+1)-xx(i),ty-yy(i))
            qrr3(k)=(tx/mbx)**2
  140     continue
          do 150 j=mbx+mby+1,mbx+2*mby
            if(xx(i).ge.tx) go to 150
            jy=j-mbx-mby
            k=3
            qrrd1(i,j,k)=q(aa(i),bb(i),cc(i),tx-xx(i),by(jy)-yy(i))
            qrrd2(i,j,k)=q(aa(i),bb(i),cc(i), tx-xx(i),by(jy+1)-yy(i))
            qrr1(i,j,k)=q(1.0d0,0.0d0,1.0d0, tx-xx(i),by(jy)-yy(i))
            qrr2(i,j,k)=q(1.0d0,0.0d0,1.0d0, tx-xx(i),by(jy+1)-yy(i))
            qrr3(k)=(ty/mby)**2
  150     continue
          do 160 j=mbx+2*mby+1,2*(mbx+mby)
            if(yy(i).le.0.0) go to 160
            jx=j-mbx-2*mby
            k=4
            qrrd1(i,j,k)=q(aa(i),bb(i),cc(i), bx(jx)-xx(i),yy(i))
            qrrd2(i,j,k)=q(aa(i),bb(i),cc(i), bx(jx+1)-xx(i),yy(i))
            qrr1(i,j,k)=q(1.0d0,0.0d0,1.0d0, bx(jx)-xx(i),yy(i))
            qrr2(i,j,k)=q(1.0d0,0.0d0,1.0d0, bx(jx+1)-xx(i),yy(i))
            qrr3(k)=(tx/mbx)**2
  160     continue
  120   continue
      end if
      if(ilin.eq.1) then
chistory
        do 170 i=1,nn
          if(ind(i).eq.1)then
            do 200 j=1,mby
              rrd1=qrrd1(i,j,1)
              rrd2=qrrd2(i,j,1)
              rr1=qrr1(i,j,1)
              rr2=qrr2(i,j,1)
              rr3=qrr3(1)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft+b2(i)/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  200       continue
            do 210 j=mby+1,mbx+mby
              rrd1=qrrd1(i,j,2)
              rrd2=qrrd2(i,j,2)
              rr1=qrr1(i,j,2)
              rr2=qrr2(i,j,2)
              rr3=qrr3(2)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft+b2(i)/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  210       continue
            do 220 j=mbx+mby+1,mbx+2*mby
              rrd1=qrrd1(i,j,3)
              rrd2=qrrd2(i,j,3)
              rr1=qrr1(i,j,3)
              rr2=qrr2(i,j,3)
              rr3=qrr3(3)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft+b2(i)/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  220       continue
            do 240 j=mbx+2*mby+1,2*(mbx+mby)
              rrd1=qrrd1(i,j,4)
              rrd2=qrrd2(i,j,4)
              rr1=qrr1(i,j,4)
              rr2=qrr2(i,j,4)
              rr3=qrr3(4)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft+b2(i)/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  240       continue
          endif
  170   continue
chistory
        do 175 i=1,ntstar
          if(ind(i).eq.1)then
            do 205 j=1,mby
              rrd1=qrrd1(i,j,1)
              rrd2=qrrd2(i,j,1)
              rr1=qrr1(i,j,1)
              rr2=qrr2(i,j,1)
              rr3=qrr3(1)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft-b2(i)/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  205       continue
            do 215 j=mby+1,mbx+mby
              rrd1=qrrd1(i,j,2)
              rrd2=qrrd2(i,j,2)
              rr1=qrr1(i,j,2)
              rr2=qrr2(i,j,2)
              rr3=qrr3(2)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft-b2(i)/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  215       continue
            do 225 j=mbx+mby+1,mbx+2*mby
              rrd1=qrrd1(i,j,3)
              rrd2=qrrd2(i,j,3)
              rr1=qrr1(i,j,3)
              rr2=qrr2(i,j,3)
              rr3=qrr3(3)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft-b2(i)/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  225       continue
            do 245 j=mbx+2*mby+1,2*(mbx+mby)
              rrd1=qrrd1(i,j,4)
              rrd2=qrrd2(i,j,4)
              rr1=qrr1(i,j,4)
              rr2=qrr2(i,j,4)
              rr3=qrr3(4)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft-b2(i)/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  245       continue
          endif
  175   continue
c
      else
chistory
        do 250 i=1,nn
          if(ind(i).eq.1)then
            do 260 j=1,mby
              rrd1=qrrd1(i,j,1)
              rrd2=qrrd2(i,j,1)
              rr1=qrr1(i,j,1)
              rr2=qrr2(i,j,1)
              rr3=qrr3(1)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft+b2(i)/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
              d2ft(i)=d2ft(i)+1/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  260       continue
            do 270 j=mby+1,mbx+mby
              rrd1=qrrd1(i,j,2)
              rrd2=qrrd2(i,j,2)
              rr1=qrr1(i,j,2)
              rr2=qrr2(i,j,2)
              rr3=qrr3(2)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft+b2(i)/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
              d2ft(i)=d2ft(i)+1/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  270       continue
            do 280 j=mbx+mby+1,mbx+2*mby
              rrd1=qrrd1(i,j,3)
              rrd2=qrrd2(i,j,3)
              rr1=qrr1(i,j,3)
              rr2=qrr2(i,j,3)
              rr3=qrr3(3)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft+b2(i)/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
              d2ft(i)=d2ft(i)+1/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  280       continue
            do 290 j=mbx+2*mby+1,2*(mbx+mby)
              rrd1=qrrd1(i,j,4)
              rrd2=qrrd2(i,j,4)
              rr1=qrr1(i,j,4)
              rr2=qrr2(i,j,4)
              rr3=qrr3(4)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft+b2(i)/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
              d2ft(i)=d2ft(i)+1/(1.d0-a5)*de2pi
     &          *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  290       continue
          endif
  250   continue
chistory
        do 255 i=1,ntstar
          if(ind(i).eq.1)then
            do 265 j=1,mby
              rrd1=qrrd1(i,j,1)
              rrd2=qrrd2(i,j,1)
              rr1=qrr1(i,j,1)
              rr2=qrr2(i,j,1)
              rr3=qrr3(1)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft-b2(i)/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
              d2ft(i)=d2ft(i)-1/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  265       continue
            do 275 j=mby+1,mbx+mby
              rrd1=qrrd1(i,j,2)
              rrd2=qrrd2(i,j,2)
              rr1=qrr1(i,j,2)
              rr2=qrr2(i,j,2)
              rr3=qrr3(2)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft-b2(i)/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
              d2ft(i)=d2ft(i)-1/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  275       continue
            do 285 j=mbx+mby+1,mbx+2*mby
              rrd1=qrrd1(i,j,3)
              rrd2=qrrd2(i,j,3)
              rr1=qrr1(i,j,3)
              rr2=qrr2(i,j,3)
              rr3=qrr3(3)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft-b2(i)/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
              d2ft(i)=d2ft(i)-1/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  285       continue
            do 295 j=mbx+2*mby+1,2*(mbx+mby)
              rrd1=qrrd1(i,j,4)
              rrd2=qrrd2(i,j,4)
              rr1=qrr1(i,j,4)
              rr2=qrr2(i,j,4)
              rr3=qrr3(4)
              ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
              ddd=sign(min(abs(ddd),1.0d0),ddd)
              delta=dacos(ddd)
              rr=(rrd1+rrd2)/2
              scl=exp(2*a4*zmg(i))
              spac2=pi/(1-a7)*scl*( (rr/scl+a6)**(1-a7)-a6**(1-a7) )
              de2pi=delta/(2*pi)
              ft=ft-b2(i)/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
              d2ft(i)=d2ft(i)-1/(1.d0-a5)*de2pi
     &          *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))*spac2
  295       continue
          endif
  255   continue
c
      endif
      vf=0.0
      do 180 id=1,ndelo
        do 190 k=1,3
          idk=ldt(id,k)
          vf=vf+b1(idk)/3*adt(id)*(tz-tstart)
          d1ft(idk)=d1ft(idk)+1.d0/3*adt(id)*(tz-tstart)
  190   continue
  180 continue
      ft=vf+ft
      f=ff-ft
      do 440 i=1,n/2
        g1(i)=d1ff(i)-d1ft(i)
        g2(i)=d2ff(i)-d2ft(i)
  440 continue
      f=-f
      do 330 k=1,n/2
        g(k      )=-g1(k)*b1(k)
        g(k+n/2  )=-g2(k)*b2(k)
  330 continue
    3 format(1h ,110x,d18.10)
    1 format(1h ,8d13.6)
c
c     return
c
      fp=0.0
c 'n' is defined in subroutine funcsb.
c the following bp1-bp7 are consistent with the originals
      do 230 i=1,n/2
        bp1(i)=b(i)
        bp2(i)=b(i+n/2)
  230 continue
      call penalt(w1,w00,bp1,fp1,gp1)
      call penalt(w2,w00,bp2,fp2,gp2)
c     write(6,*) f,fp1,fp2,f+fp1+fp2
      f=f+fp1+fp2
      if(fp.lt.0.0) write(6,*) 'fp < 0 !',  fp
      if(fp.lt.0.0) write(6,1) (b(i),i=1,n)
!     if(fp.lt.0.0) write(9,*) 'fp < 0 !',  fp
!     if(fp.lt.0.0) write(9,1) (b(i),i=1,n)
c
      sum=0.
      do 350 ia=1,n/2
        h(ia      )=g(ia      )+gp1(ia)
        h(ia  +n/2)=g(ia  +n/2)+gp2(ia)
        sum=sum+h(ia)**2
        sum=sum+h(ia+n/2)**2
  350 continue
      index=1
      return
   50 continue
      ifg=1
      if(myrank.eq.0) write(6,*) 'jumped down',i,1+i*5/n,b(i)
!     if(myrank.eq.0) write(9,*) 'jumped down',i,1+i*5/n,b(i)
      f=1.0d30
      return
      end
c***********************************************************************
      subroutine hesian(n,b,au,ad)
      implicit real * 8 (a-h,o-z)
      parameter(npa=22222,nda=22222,mdl=43333,nr=020,nmax=6000)
      parameter(npall=5*npa)
      common /data/ xx(nda),yy(nda),zz(nda),zmg(nda),
     &              nn,mm,tx,ty,tz, ind(nda)
      common /coef/ aa(nda),bb(nda),cc(nda)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /delink/ nzu(npa),iu(nr,npa)
      common /delin5/ nzu5(npall),iu5(nr,npall)
      common /fixpar/ a1,a2,a3,aa4,a5,a6,a7
      common /sublin/ilin
      common /ellips/qrrd1(nda,128,4),qrrd2(nda,128,4),qrr1(nda,128,4),
     &               qrr2(nda,128,4),qrr3(4),mbx,mby
      common /hesapprox/nhesapp
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      common /jfun/rmax(nda),jxy(nda,nmax),jnn(nda)
      common /range/tstart,ntstar
      dimension b(npall)
      dimension b1(npa),d1ff(npa),d1ft(npa),g1(npa)
      dimension b2(npa),d2ff(npa),d2ft(npa),g2(npa)
      dimension rmd2i(nda)
      dimension h22i(nda)
      dimension h11ff(nda),h22ff(nda)
      dimension h11ft(nda),h22ft(nda)
      dimension h11(nda),h22(nda)
      dimension au(nr,npall),ad(npall),aup(nr,npall),adp(npall)
      common /prihe5/ adp1(npa),aup1(nr,npa),adp2(npa),aup2(nr,npa),
     & adp4(npa),aup4(nr,npa),adp5(npa),aup5(nr,npa),
     & adp7(npa),aup7(nr,npa)
      data pi /3.14159265358979d0/
c
      do 410 i=1,n/2
        nzu5(i      )=nzu(i)
        nzu5(i+  n/2)=nzu(i)
        do 420 j=1,nzu(i)
          iu5(j,i      )=iu(j,i)
          iu5(j,i+  n/2)=iu(j,i)+  n/2
  420   continue
  410 continue
      a4=aa4/2      
      do 210 i=1,n/2
        b1(i)=a1*exp(b(i))
        b2(i)=a2*exp(b(i+n/2))
        d1ff(i)=0.0
        d2ff(i)=0.0
        h11ff(i)=0.0
        h22ff(i)=0.0
  210 continue
      ff=0.0
      ramdai=b1(1)
      if(ramdai.gt.0.0) then
        ff=log(ramdai)
        d1ff(1)=1.0d0/ramdai
        h11ff(1)=-1.0d0/ramdai**2
      endif
chistory
      do 10 i=2,nn
        if(ind(i).eq.1) then
          rmdi=0.0
          ji=jnn(i)
          do 20 j1=1,ji
            j=jxy(i,j1)
            rmd2i(j)=0.0
            h22i(j)=0.0
            xxij=(xx(i)-xx(j))*exp(-a4*zmg(j))
            yyij=(yy(i)-yy(j))*exp(-a4*zmg(j))
            spac0=aa(j)*xxij**2-2*bb(j)*xxij*yyij+cc(j)*yyij**2+a6
            space=spac0**(-a7)
            rmdi=rmdi+b2(j)/(zz(i)-zz(j)+a3)**a5*space
            rmd2i(j)=rmd2i(j)+(zz(i)-zz(j)+a3)**(-a5)*space
            h22i(j)=h22i(j)+0.d0/(zz(i)-zz(j)+a3)**a5*space
   20     continue
chistory
          if(zz(i).lt.tstart) go to 10
c
          ramdai=b1(i)+rmdi
          ft=ft+log(ramdai)
          d1ff(i)=d1ff(i)+1.0d0/ramdai
          h11ff(i)=h11ff(i)-1.0d0/ramdai**2
          ji=jnn(i)
          do 60 j1=1,ji
            j=jxy(i,j1)
            d2ff(j)=d2ff(j)+rmd2i(j)/ramdai
            h22ff(j)=h22ff(j)+h22i(j)/ramdai-rmd2i(j)*rmd2i(j)/ramdai**2
   60     continue
        endif
   10 continue
c
      ft=0.0
      do 30 i=1,n/2
        d1ft(i)=0.0
        d2ft(i)=0.0
        h22ft(i)=0.0
   30 continue
chistory
      do 70 i=1,nn
        if(ind(i).eq.1)then
          do 80 j=1,mby
            rrd1=qrrd1(i,j,1)
            rrd2=qrrd2(i,j,1)
            rr1=qrr1(i,j,1)
            rr2=qrr2(i,j,1)
            rr3=qrr3(1)
            ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
            ddd=sign(min(abs(ddd),1.0d0),ddd)
            delta=dacos(ddd)
            rr=(rrd1+rrd2)/2
            scl=exp(2*a4*zmg(i))
            rrs6=rr/scl+a6
            spainte=pi/(1-a7)*scl*(rrs6**(1-a7)-a6**(1-a7))
            ft=ft+b2(i)/(1-a5)*spainte*delta/(2*pi)
     &        *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
            d2ft(i)=d2ft(i)+1/(1-a5)*spainte*delta/(2*pi)
     &        *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
   80     continue
          do 90 j=mby+1,mbx+mby
            rrd1=qrrd1(i,j,2)
            rrd2=qrrd2(i,j,2)
            rr1=qrr1(i,j,2)
            rr2=qrr2(i,j,2)
            rr3=qrr3(2)
            ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
            ddd=sign(min(abs(ddd),1.0d0),ddd)
            delta=dacos(ddd)
            rr=(rrd1+rrd2)/2
            scl=exp(2*a4*zmg(i))
            rrs6=rr/scl+a6
            ft=ft+b2(i)/(1-a5)*spainte*delta/(2*pi)
     &        *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
            d2ft(i)=d2ft(i)+1/(1-a5)*spainte*delta/(2*pi)
     &        *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
   90     continue
          do 100 j=mbx+mby+1,mbx+2*mby
            rrd1=qrrd1(i,j,3)
            rrd2=qrrd2(i,j,3)
            rr1=qrr1(i,j,3)
            rr2=qrr2(i,j,3)
            rr3=qrr3(3)
            ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
            ddd=sign(min(abs(ddd),1.0d0),ddd)
            delta=dacos(ddd)
            rr=(rrd1+rrd2)/2
            scl=exp(2*a4*zmg(i))
            rrs6=rr/scl+a6
            spainte=pi/(1-a7)*scl*(rrs6**(1-a7)-a6**(1-a7))
            ft=ft+b2(i)/(1-a5)*spainte*delta/(2*pi)
     &        *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
            d2ft(i)=d2ft(i)+1/(1-a5)*spainte*delta/(2*pi)
     &        *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
  100     continue
          do 110 j=mbx+2*mby+1,2*(mbx+mby)
            rrd1=qrrd1(i,j,4)
            rrd2=qrrd2(i,j,4)
            rr1=qrr1(i,j,4)
            rr2=qrr2(i,j,4)
            rr3=qrr3(4)
            ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
            ddd=sign(min(abs(ddd),1.0d0),ddd)
            delta=dacos(ddd)
            rr=(rrd1+rrd2)/2
            scl=exp(2*a4*zmg(i))
            rrs6=rr/scl+a6
            spainte=pi/(1-a7)*scl*(rrs6**(1-a7)-a6**(1-a7))
            ft=ft+b2(i)/(1-a5)*spainte*delta/(2*pi)
     &        *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
            d2ft(i)=d2ft(i)+1/(1-a5)*spainte*delta/(2*pi)
     &        *((tz-zz(i)+a3)**(1-a5)-a3**(1-a5))
  110     continue
        endif
   70 continue
chistory
      do 75 i=1,ntstar
        if(ind(i).eq.1)then
          do 85 j=1,mby
            rrd1=qrrd1(i,j,1)
            rrd2=qrrd2(i,j,1)
            rr1=qrr1(i,j,1)
            rr2=qrr2(i,j,1)
            rr3=qrr3(1)
            ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
            ddd=sign(min(abs(ddd),1.0d0),ddd)
            delta=dacos(ddd)
            rr=(rrd1+rrd2)/2
            scl=exp(2*a4*zmg(i))
            rrs6=rr/scl+a6
            spainte=pi/(1-a7)*scl*(rrs6**(1-a7)-a6**(1-a7))
            ft=ft-b2(i)/(1-a5)*spainte*delta/(2*pi)
     &        *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))
            d2ft(i)=d2ft(i)-1/(1-a5)*spainte*delta/(2*pi)
     &        *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))
   85     continue
          do 95 j=mby+1,mbx+mby
            rrd1=qrrd1(i,j,2)
            rrd2=qrrd2(i,j,2)
            rr1=qrr1(i,j,2)
            rr2=qrr2(i,j,2)
            rr3=qrr3(2)
            ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
            ddd=sign(min(abs(ddd),1.0d0),ddd)
            delta=dacos(ddd)
            rr=(rrd1+rrd2)/2
            scl=exp(2*a4*zmg(i))
            rrs6=rr/scl+a6
            ft=ft-b2(i)/(1-a5)*spainte*delta/(2*pi)
     &        *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))
            d2ft(i)=d2ft(i)-1/(1-a5)*spainte*delta/(2*pi)
     &        *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))
   95     continue
          do 105 j=mbx+mby+1,mbx+2*mby
            rrd1=qrrd1(i,j,3)
            rrd2=qrrd2(i,j,3)
            rr1=qrr1(i,j,3)
            rr2=qrr2(i,j,3)
            rr3=qrr3(3)
            ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
            ddd=sign(min(abs(ddd),1.0d0),ddd)
            delta=dacos(ddd)
            rr=(rrd1+rrd2)/2
            scl=exp(2*a4*zmg(i))
            rrs6=rr/scl+a6
            spainte=pi/(1-a7)*scl*(rrs6**(1-a7)-a6**(1-a7))
            ft=ft-b2(i)/(1-a5)*spainte*delta/(2*pi)
     &        *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))
            d2ft(i)=d2ft(i)-1/(1-a5)*spainte*delta/(2*pi)
     &        *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))
  105     continue
          do 115 j=mbx+2*mby+1,2*(mbx+mby)
            rrd1=qrrd1(i,j,4)
            rrd2=qrrd2(i,j,4)
            rr1=qrr1(i,j,4)
            rr2=qrr2(i,j,4)
            rr3=qrr3(4)
            ddd=(rr1+rr2-rr3)/( 2*sqrt(rr1)*sqrt(rr2) )
            ddd=sign(min(abs(ddd),1.0d0),ddd)
            delta=dacos(ddd)
            rr=(rrd1+rrd2)/2
            scl=exp(2*a4*zmg(i))
            rrs6=rr/scl+a6
            spainte=pi/(1-a7)*scl*(rrs6**(1-a7)-a6**(1-a7))
            ft=ft-b2(i)/(1-a5)*spainte*delta/(2*pi)
     &        *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))
            d2ft(i)=d2ft(i)-1/(1-a5)*spainte*delta/(2*pi)
     &        *((tstart-zz(i)+a3)**(1-a5)-a3**(1-a5))
  115     continue
        endif
   75 continue
c
      vf=0.0
      do 180 id=1,ndelo
        do 190 k=1,3
          idk=ldt(id,k)
          vf=vf+b1(idk)/3*adt(id)*(tz-tsatrt)
          d1ft(idk)=d1ft(idk)+1.d0/3 *adt(id)*(tz-tsatrt)
          h11ft(idk)=h11ft(idk)+0.0
  190   continue
  180 continue
      do 40 i=1,n/2
        g1(i)=d1ff(i)-d1ft(i)
        g2(i)=d2ff(i)-d2ft(i)
        if(nhesapp.eq.0) then
        g1(i)=0.0
        g2(i)=0.0
        endif
        h11(i)=h11ff(i)-h11ft(i)
        h22(i)=h22ff(i)-h22ft(i)
   40 continue
      do 50 k=1,n/2
        ad(k      )=-h11(k)*b1(k)**2-g1(k)*b1(k)
        ad(k+  n/2)=-h22(k)*b2(k)**2-g2(k)*b2(k)      
   50 continue
      do 330 k=1,npext
        adp(k        )=adp1(k)
        adp(k+npext  )=adp2(k)
        do 340 j=1,nzu(k)
          aup(j,k        )=aup1(j,k)
          aup(j,k+npext  )=aup2(j,k)
  340   continue
  330 continue
c
      sum1=0
      sum2=0
      do 220 i=1,n
        ad(i)= ad(i)+adp(i)
        sum1=sum1+ad(i)
        do 230 j=1,nzu5(i)
          au(j,i)=aup(j,i)
          sum2=sum2+au(j,i)
  230   continue
  220 continue
c     write(6,*) 'hes',sum1,sum2
      return
      end
c***********************************************************************
      subroutine davidn( x,n,funct,ff )
c       ----------------------------------------------------------------
c       the following subroutines are directly called by this subroutine
c             funct
c             hesian
c             linear
c       ----------------------------------------------------------------
c          inputs:
c             x:       vector of initial values
c             k:       dimension of the vector x
c          output:
c             x:       vector of minimizing solution
      implicit  real * 8  ( a-h , o-z )
      parameter(nda=22222,npa=22222,mdl=43333,nr=020)
      parameter(npall=npa*5)
      dimension x(npall),dx(npall),g(npall),g0(npall),y(npall)
      dimension  hd(npall),wrk(npall),s(npall)
      dimension  sneg(npall)
      dimension  x1(npall),x2(npall),g1(npall)
      common /abicc/abic,fmin
      common /swgrad/ w1,w2,w4,w5,w7,ptdet,w00,w01
      common /print/ itp
      common /delin5/ nzu5(npall),iu5(nr,npall)
      common /delink/ nzu(npa),iu(nr,npa)
      integer nzl(npall),il(nr,npall)
      real*8 au(nr,npall),ad(npall),al(nr,npall)
      real*8 uu(nr,npall),ul(nr,npall),ud(npall)
      real*8 r(npall),p(npall),ru(npall),ap(npall)
      real*8 aut(nr,npa),adt(npa),alt(nr,npa)
      real*8 uut(nr,npa),ult(nr,npa),udt(npa)
      real*8 rt(npa),pt(npa),rut(npa),apt(npa)
      real*8 yt(npa),wrkt(npa),gt(npa),snegt(npa)
      integer iv(npall)
      integer ivt(npa)
      integer nzlt(npa),ilt(nr,npa)
      common /diagnl/ idet
      common /sublin/ ilin
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      common /gradient/ igsum
      common /davtau/tau1,tau2
      common /daveps/eps1,eps2
c
      external funct
c     data  tau1 , tau2  /  0.1d-4 , 0.1d-4  /
c     data  eps1 , eps2  / 0.1d-5 , 0.1d-5  /
      data pi / 3.14159265358979d0 /
c     include'mpif.h'
      tau1=0.1d-4
      tau2=0.1d-4
      eps1=0.1d-4
      eps2=0.1d-4
      igsum=0
      iter=0
         lpr=0
c        lpr=2
c        lpr=3
      if(itp.eq.0) lpr=2
      ramda = 0.5d0
      const1 = 1.0d-70
c
c          initial estimate of inverse of hessian
c
11112 continue
c
      call  funct(n,x,xm,g,ig )
      call hesian(n,x,au,ad)
c
  337 continue
      do 335 i=1,n
        x1(i)=x(i)
  335 continue
      eps0=eps1
*
      go to 339
*
      do 338 i=1,n/2

        x1(i)=x(i)+eps0
        call  funct( n,x1,xm1,g,ig )
        write(6,*) i,'1',xm1,xm1-xm,g(i)*eps0
        x1(i)=x(i)

        x1(i+1*n/2)=x(i+1*n/2)+eps0
        call  funct( n,x1,xm1,g,ig )
        write(6,*) i,'2',xm1,xm1-xm,g(i+1*n/2)*eps0
        x1(i+1*n/2)=x(i+1*n/2)

  338 continue
  339 continue
*
      go to 334
*
      write(6,*) 'input?'
      eps0=eps1
      eps0=0.1d-2
      eps0=0.1d-5
      eps0=0.1d-9
      read(5,*) k
      xmsum=0.0
      xm1sum=0.0
      dxmsum=0.0
      gsum=0.0
      do 336 i=(k-1)*n/2+1,k*n/2
        x1(i)=x(i)+eps0
        x1(i-1)=x(i-1)
        call  funct( n,x1,xm1,g,ig )
        write(6,5) k,i,xm1-xm-g(i)*eps0,ad(i)/2*eps0**2,
     &                 xm1-xm-g(i)*eps0+ad(i)/2*eps0**2
    5 format(2i6,4d13.4)
        xmsum=xmsum+xm
        xm1sum=xm1sum+xm1
        dxmsum=dxmsum+(xm1-xm)/eps0
        gsum=gsum+g(i)
  336 continue
      write(6,*) (xm1sum-xmsum)/eps0,dxmsum,gsum
      go to 337
c      call mpi_finalize(ierr)
      stop 12345
  334 continue
c
      do  20   i=1,n
        s(i) = 0.0d00
        dx(i) = 0.0d00
        hd(i)=1.d0/ad(i)
c       hd(i)=1.d0
c     if(i.le.n/5) write(6,*) i,hd(i)
   20 continue
c     stop 1234
c
      call  funct( n,x,xm,g,ig )
c
      if(myrank.eq.0) then
      if(lpr.ge.2) write( 6,345 ) w1,w2
      if(lpr.ge.2) write(6,*)'Initial Penalized log likelihood =',xm
!     if(lpr.ge.2) write( 9,345 ) w1,w2
!     if(lpr.ge.2) write(9,*)'Initial Penalized log likelihood =',xm
c     if(lpr.ge.2) write( 6,340 ) w1,w2
      end if
c
      icc = 0
c      iteration
11110 continue
      icc = icc + 1
      do  11111   ic=1,n
      if( ic .eq. 1 .and. icc .eq. 1 )     go to 120
c
      do  40   i=1,n
   40 y(i) = g(i) - g0(i)
c
      iter=iter+1
c
      if(mod(iter,0001).eq.0) then
        call hesian(n,x,au,ad)
c       call iccg (npall, n, nr, nzu5, nzl, iu5, il,
c    $           au, al, ad, y, eps1i, wrk,
c    $           uu, ul, ud, r, p, ru, ap, iv)
c--------------------------------------------------
     
      do myrk =0, 1
          
          do 46 k=1,n/2
              adt(k)=ad(k+n/2*myrk)
              yt(k)=y(k+n/2*myrk)
              wrkt(k)=wrk(k+n/2*myrk)
              do 46 j=1,nzu(k)
                  aut(j,k)=au(j,k+n/2*myrk)
   46     continue
          call iccg (npa, n/2, nr, nzu, nzlt, iu, ilt,
     $           aut, alt, adt, yt, eps1i, wrkt,
     $           uut, ult, udt, rt, pt, rut, apt, ivt)
     
          do k=1, n/2
            y(k+n/2*myrk)= yt(k)
            wrk(k+n/2*myrk)= wrkt(k)
          enddo

       enddo

      else

        do  60   i=1,n
          wrk(i) = y(i) * hd(i)
   60   continue
      end if
c
   55 continue
      s1 = 0.0d00
      s2 = 0.0d00
      do  70   i=1,n
      s1 = s1 + wrk(i) * y(i)
   70 s2 = s2 + dx(i) * y(i)
      if( s1.le.const1 .or. s2.le.const1 ) go to 900
*     if( s1.le.const1 .or. s2.le.const1 ) then
*       write(6,*) 's1 or s2 =<const1', s1,s2,const1
*       go to 11110
*       go to 11112
*       go to 900
*     endif
      if( s1 .le. s2 ) go to 100
c
c          update the inverse of hessian matrix
c
c     ---  davidon-fletcher-powell type correction  ---
c                  omitted
c     ---  fletcher type correction  ---
c                  omitted
c
  100 continue
c
  120 continue
c
      if(iter.gt.0 .and. mod(iter,0001).eq.0) then
c
c       call hesian(n,x,au,ad)
c
c       call iccg (npall, n, nr, nzu5, nzl, iu5, il,
c    $           au, al, ad, g, eps, sneg,
c    $           uu, ul, ud, r, p, ru, ap, iv)
       do myrk=0, 1
          do 146 k=1,n/2
             adt(k)=ad(k+n/2*myrk)
             gt(k)=g(k+n/2*myrk)
             snegt(k)=sneg(k+n/2*myrk)
             do 146 j=1,nzu(k)
                aut(j,k)=au(j,k+n/2*myrk)
  146     continue
          call iccg (npa, n/2, nr, nzu, nzlt, iu, ilt,
     $           aut, alt, adt, gt, eps, snegt,
     $           uut, ult, udt, rt, pt, rut, apt, ivt)
          do k=1, n/2
             g(k+n/2*myrk)=gt(k)
             sneg(k+n/2*myrk)=snegt(k)
          enddo

       enddo 
   
       myrk = 0
 


        ss=0.0
        do 155 i=1,n
          s(i) = -sneg(i)
          ss=ss+s(i)*s(i)
  155   continue
c     end if
      else
c
      ss = 0.0d00
      do  150   i=1,n
      sum = hd(i)*g(i)
      ss = ss + sum * sum
  150 s(i) = -sum
      end if
c
      s1 = 0.0d00
      s2 = 0.0d00
      do  170   i=1,n
      s1 = s1 + s(i)*g(i)
  170 s2 = s2 + g(i)*g(i)
      ds2 = sqrt(s2)
      gtem = abs(s1) / ds2
      if( gtem .le. tau1  .and.  ds2 .le. tau2 ) go to 900
      if( s1 .lt. 0.0d00 ) go to  200
c
      call hesian(n,x,au,ad)
c
      do  190   i=1,n
        hd(i) = 1.d0/ad(i)
        s(i) = -s(i)
  190 continue
  200 continue
c
      ed = xm
c
c          linear  search
c
      ilin=1
      call linear( x,s,ramda,ed,n,ig,funct )
      ilin=0
c
      if(myrank.eq.0) then
        if(lpr.ge.2) write( 6,330 )     ramda , ed , s1 , s2
!       if(lpr.ge.2) write( 9,330 )     ramda , ed , s1 , s2
      end if
      prints2=s2
 1010 continue
c
      s1 = 0.0d00
      do  210   i=1,n
      dx(i) = s(i) * ramda
      s1 = s1 + dx(i) * dx(i)
      g0(i) = g(i)
  210 x(i) = x(i) + dx(i)
      xmb = xm
c
      call  funct( n,x,xm,g,ig )
c
      s2 = 0.d0
      do  220     i=1,n
  220 s2 = s2 + g(i)*g(i)
      if( sqrt(s2) .gt. tau2 )   go to  11111
*     if( sqrt(s2) .gt. tau2 ) then
*       write(6,*) 'sqrt(s2)>tau2 ',sqrt(s2),tau2
*       write(6,*) 'go to 11111, icc =', icc
*       go to 11111
*     endif 
      if( xmb/xm-1.d0 .lt. eps1  .and.  sqrt(s1) .lt. eps2 )  go to 900
*     if( xmb/xm-1.d0 .lt. eps1  .or.  sqrt(s1) .lt. eps2 ) then
*       write(6,*) 'xmb/xm-1.d0 =',xmb/xm-1.d0
*       write(6,*) 'sqrt(s1) =',sqrt(s1)
*       go to 900
*     endif 
11111 continue
      if( icc .ge. 5 )     go to 900
      go to 11110
  900 continue
      if(prints2.le.0.1d-3) igsum=1
      ff=xm
      if(myrank.eq.0) then
      if(lpr.ge.2) write(6,350) ff
!     if(lpr.ge.2) write(9,350) ff
  350 format(1h ,'penalized log likelihood =',d23.15)
      if(lpr.ge.3) write( 6,600 )
      if(lpr.ge.3) write( 6,610 )     (x(i),i=1,n)
      if(lpr.ge.3) write( 6,601 )
      if(lpr.ge.3) write( 6,610 )     (g(i),i=1,n)
      end if
c
      ff=xm-ptdet/2+(n-1)*log(sqrt(2*pi))
c
      call hesian(n,x,au,ad)
c
c     call icdata(n,nr,h,au,ad,y,wrk)
c
      idet=1
      call iccg (npall, n, nr, nzu5, nzl, iu5, il,
     $           au, al, ad, y, tdet, wrk,
     $           uu, ul, ud, r, p, ru, ap, iv)
      idet=0
c
      pn=tdet
      abic=2*ff+pn-(n-1)*log(2*pi)
      if(myrank.eq.0) then
      write(6,599) w1,w2
      write(6,598) abic,ff,pn
!     write(9,599) w1,w2
!     write(9,598) abic,ff,pn
      end if
  599 format(1h ,'#e: w1,w2 =',2d11.3)
  598 format(1h ,' abic =',d17.10 ,' -l =',d17.10,' pn =',d17.10)
c 598 format(1h ,' abic =',d17.10 ,' -l =',d17.10,' pn =',d17.10/)
c
c     if(mod(iter,5).eq.0) go to 1010
      return
  330 format( 1h ,'lambda =',d15.7,4x,'pell =',d23.15,2x,d9.2,2x,d9.2)
  340 format( 1h ,'w1,w2 =',2d11.3,3x,'pell =',d23.15)
  345 format( 1h ,'#s: w1,w2 =',2d11.3)
  600 format( 1h ,'-----  x  -----' )
  601 format( 1h ,'***  gradient  ***' )
  610 format( 1h ,10d13.5 )
      end
************************************************************************
      subroutine  linear( x,h,ram,ee,k,ig,funct )
c
c     this subroutine performs the linear search along the direction spe
c     by the vector h
c       ----------------------------------------------------------------
c       the following subroutine is directly called by this subroutine:
c             funct
c       ----------------------------------------------------------------
c
c     inputs:
c        x:       vector of position
c        h:       search direction
c        k:       dimension of vector x
c
c     outputs:
c        ram:     optimal step width
c        e2:      minimum function value
c        ig:      error code
c
      implicit  real  *8 ( a-h,o-z )
      parameter(nda=22222,npa=22222,mdl=43333,nr=020)
      parameter(npall=npa*5)
      integer  return,sub
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      common /ccc/ isw,ipr
      dimension  x(npall),h(npall),x1(npall)
      dimension  g(npall)
      external funct
c
c     ipr = 0
c     ipr = 7
      if( ram .le. 1.0d-30 )  ram = 0.01d0
      const2 = 1.0d-60
      hnorm = 0.d0
      do 10  i=1,k
   10 hnorm = hnorm + h(i)**2
      hnorm = sqrt( hnorm )
c
      ram2 = ram
      e1 =ee
      ram1 = 0.d0
c
      do 20  i=1,k
   20 x1(i) = x(i) + ram2*h(i)
      call  funct( k,x1,e2,g,ig )
      if(ipr.ge.7)  write(6,2)  ram2,e2
c
      if( ig .eq. 1 )  go to  50
      if( e2 .gt. e1 )  go to 50
   30 ram3 = ram2*2.d0
      do 40  i=1,k
   40 x1(i) = x(i) + ram3*h(i)
      call  funct( k,x1,e3,g,ig )
      if( ig.eq.1 )  go to  500
      if( ipr.ge.7 )  write(6,3)  ram3,e3
      if( e3 .gt. e2 )  go to 70
      ram1 = ram2
      ram2 = ram3
      e1 = e2
      e2 = e3
      go to 30
c
   50 ram3 = ram2
      e3 = e2
      ram2 = ram3*0.1d0
      if( ram2*hnorm .lt. const2 )  go to  400
      do 60  i=1,k
   60 x1(i) = x(i) + ram2*h(i)
      call  funct( k,x1,e2,g,ig )
      if(ipr.ge.7)  write(6,4)  ram2,e2
c     write(6,8) e1,e2,e3,ee
    8 format(1h ,5d20.10)
      if( e2.gt.e1 )  go to 50
c
   70 assign 80 to return
      go to 200
c
   80 do 90  i=1,k
   90 x1(i) = x(i) + ram*h(i)
      call  funct( k,x1,ee,g,ig )
      if(ipr.ge.7)  write(6,5)  ram,ee
c
      ifg = 0
      assign  300 to  sub
      assign 200 to sub
   95 assign 130 to return
      if( ram .gt. ram2 )  go to 110
      if( ee .ge. e2 )  go to 100
      ram3 = ram2
      ram2 = ram
      e3 =e2
      e2 =ee
      go to  sub,( 200,300 )
c
  100 ram1 = ram
      e1 = ee
      go to  sub,( 200,300 )
c
  110 if( ee .le. e2 )  go to 120
      ram3 = ram
      e3 = ee
      go to  sub,( 200,300 )
c
  120 ram1 = ram2
      ram2 = ram
      e1 = e2
      e2 = ee
      go to  sub,( 200,300 )
c
  130 do 140  i=1,k
  140 x1(i) = x(i) + ram*h(i)
      call  funct( k,x1,ee,g,ig )
      if( ipr.ge.7 )  write(6,6)  ram,ee
      assign 200 to sub
      ifg = ifg+1
      ifg = 0
      if( ifg .eq. 1 )  go to 95
c
      if( e2 .lt. ee )  ram = ram2
      if( e2 .lt. ee )  ee = e2
      return
c
c      -------  internal subroutine sub1  -------
  200 a1 = (ram3-ram2)*e1
      a2 = (ram1-ram3)*e2
      a3 = (ram2-ram1)*e3
      b2 = (a1+a2+a3)*2.d0
      b1 = a1*(ram3+ram2) + a2*(ram1+ram3) + a3*(ram2+ram1)
      if( b2 .eq. 0.d0 )  go to 210
      ram = b1 /b2
      go to return ,( 80,130 )
c
  210 ig = 1
      ram = ram2
      return
c
c      -------  internal subroutine sub2  -------
c
  300 if( ram3-ram2 .gt. ram2-ram1 )  go to 310
      ram = (ram1+ram2)*0.5d0
      go to return ,( 80,130 )
c
  310 ram = (ram2+ram3)*0.5d0
      go to return ,( 80,130 )
c ------------------------------------------------------------
c
  400 ram = 0.d0
      return
c ------------------------------------------------------------
c
  500 ram = (ram2+ram3)*0.5d0
  510 do 520  i=1,k
  520 x1(i) = x(i) + ram*h(i)
      call  funct( k,x1,e3,g,ig )
      if( ipr.ge.7 )  write(6,7)  ram,e3
      if( ig.eq.1 )  go to 540
      if( e3.gt.e2 )  go to 530
      ram1 = ram2
      ram2 = ram
      e1 = e2
      e2 = e3
      go to 500
c
  530 ram3 = ram
      go to 70
c
  540 ram = (ram2+ram)*0.5d0
      go to 510
c
c ------------------------------------------------------------
    1 format( 1h ,'lambd1 =',d18.10, 10x,'e1 =',d25.17 )
    2 format( 1h ,'lambd2 =',d18.10, 10x,'e2 =',d25.17 )
    3 format( 1h ,'lambd3 =',d18.10, 10x,'e3 =',d25.17 )
    4 format( 1h ,'lambd4 =',d18.10, 10x,'e4 =',d25.17 )
    5 format( 1h ,'lambd5 =',d18.10, 10x,'e5 =',d25.17 )
    6 format( 1h ,'lambd6 =',d18.10, 10x,'e6 =',d25.17 )
    7 format( 1h ,'lambd7 =',d18.10, 10x,'e7 =',d25.17 )
      end
************************************************************************
      subroutine funcsb(n,xx,f)
      implicit real*8(a-h,o-z)
      parameter(nda=22222,npa=22222,mdl=43333,nr=020)
      parameter(npall=npa*5)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /print/ itp
      common /abicc/abic,fmin
      common /swgrad/ w1,w2,w4,w5,w7,ptdet,w00,w01
      common /prihe5/ adp1(npa),aup1(nr,npa),adp2(npa),aup2(nr,npa),
     & adp4(npa),aup4(nr,npa),adp5(npa),aup5(nr,npa),
     & adp7(npa),aup7(nr,npa)
      common /xad/xd(npall),m
      common /sublin/ ilin
      common /prupd/ iupd
      common /delink/ nzu(npa),iu(nr,npa)
      common /fixpar/ a1,a2,a3,aa4,a5,a6,a7
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      common /hesapprox/nhesapp
      integer nzl(npa),il(nr,npa)
      real*8 uu(nr,npa),ul(nr,npa),ud(npa),al(nr,npa)
      real*8 r(npa),p(npa),ru(npa),ap(npa),y(npa),wrk(npa)
      integer iv(npa)
      dimension xdd(npall)
      common /diagnl/ idet
      common /gradient/ igsum
c
      external func17
      dimension xx(7),x(npall)
      data nzl/npa*0/
c     data fmin/1.d60/
c     save fmin
c     ifg=0
c     eps=0.1d-2
c
c     if(myrank.eq.0) write(6,*) 'f',(xx(i),i=1,7)
      w1=exp(xx(1))
      w2=exp(xx(2))
c     if(myrank.eq.0) write(6,*) w1,w2
cf.  see main.
      a3=exp(xx(3))
      aa4=exp(xx(4))
      a5=exp(xx(5))
c     a5=exp(xx(5))+1
      a6=exp(xx(6))
      a7=exp(xx(7))
c     a7=exp(xx(7))+1

c
      idet=1
      call prihes(npext,w1,w00,adp1,aup1)
      call iccg (npa, m/2-1, nr, nzu, nzl, iu, il,
     $           aup1, al, adp1, y, ptdet1, wrk,
     $           uu, ul, ud, r, p, ru, ap, iv)
      call prihes(npext,w2,w00,adp2,aup2)
      call iccg (npa, m/2-1, nr, nzu, nzl, iu, il,
     $           aup2, al, adp2, y, ptdet2, wrk,
     $           uu, ul, ud, r, p, ru, ap, iv)
      ptdet=ptdet1+ptdet2

      if(myrank.eq.0) then
      write(6,*) ptdet1,ptdet2
      write(6,3) ptdet
!     write(9,*) ptdet1,ptdet2
!     write(9,3) ptdet
      end if
    3 format(1h ,'ptdet =',d20.13)
c
      idet=0
      call prihes(npext,w1,w00,adp1,aup1)
      call prihes(npext,w2,w00,adp2,aup2)

      itp=0
      do 60 i=1,m
   60 x(i)=xd(i)
c     do 20 i=1,1
      do 20 i=1,33
        call davidn( x,m,func17,xm )
        if(myrank.eq.0) write(6,*) 'repeated davidn =', i
!       if(myrank.eq.0) write(9,*) 'repeated davidn =', i
        if(igsum.eq.1) go to 25
   20 continue
   25 continue
ckkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk
c     nhesapp=1
      if(fmin.lt.abic) go to 50
      iupd=iupd+1
      fmino=fmin
      if(iupd.eq.1) fmino=abic
      fmin=abic


      sss1=0
      sss2=0
*
      do i=1, m/2
          sss1 = sss1 + x(i)
          sss2 = sss2 + x(i + m / 2)
      enddo
      if(myrank.eq.0) then
        write(*,*)'Surface Sliding: Old a1, a2=', a1, a2
!       write(9,*)'Surface Sliding: Old a1, a2=', a1, a2
      endif
      a1=a1*exp(sss1/(m/2))
      a2=a2*exp(sss2/(m/2))
      if(myrank.eq.0) then
          write(*,*)'Surface Sliding: sss1, sss2=', sss1, sss2
          write(*,*)'Surface Sliding: New a1, a2=', a1, a2
!         write(9,*)'Surface Sliding: sss1, sss2=', sss1, sss2
!         write(9,*)'Surface Sliding: New a1, a2=', a1, a2
      endif
*
      do 10 i=1,m/2
          xd(i)=x(i)-sss1/m*2.0d0
   10   xd(i + m /2)=x( i + m/2) - sss2/m*2d0
      if(myrank.eq.0) then
      write(6,617) '-------- xd -------',iupd,abic,abic-fmino
!     write(9,617) '-------- xd -------',iupd,abic,abic-fmino
  617 format(a,i6,e20.12,f11.3)
      write(6,615) a1,a2,a3,aa4,a5,a6,a7
!     write(9,615) a1,a2,a3,aa4,a5,a6,a7
  615 format(1h ,'a1-7',7e10.3)
      write(6,616) w1,w2
!     write(9,616) w1,w2
  616 format(1h ,'w1-2',5e10.3)
      write(6,*)
!     write(9,*)
c     write(6,610) (xd(i),i=1,m)
  610 format(1h ,10f7.2)
      open(11,file='hist-etas-mk.upda')
      write(11,5) w1,w2,abic,m
      write(11,6) a1,a2,a3,aa4,a5,a6,a7
      write(11,6) (xd(i),i=1,m)
      call flush(11)
      close(11)
      do 210 i=1,m/2
        xdd(i)=log(a1)+xd(i)
        xdd(i+m/2)=a2*exp(xd(i+m/2))
  210 continue
      open(21,file='hist-etas-mk.omap')
      write(21,5) w1,w2,abic,m
      write(21,6) a1,a2,a3,aa4,a5,a6,a7
      write(21,6) (xdd(i),i=1,m)
      close(21)
      end if
c     stop 119
    5 format(2e15.5,f11.2,i7)
    6 format(4e20.12)
   50 continue 
c     if(myrank.eq.0)then
c        write(*,9931) w1, w2, abic
c        write(*,9932) a1, a2, a3, aa4, a5, a6, a7
c        write(9,9931) w1, w2, abic
c        write(9,9932) a1, a2, a3, aa4, a5, a6, a7
c     endif
c9931 format('w1, w2, abic=', 3d15.7)
c9932 format(' a1 ... a7=', 7d15.7)
      f=abic
c     if(ilin.eq.1) return
      return
c 100 continue
c     f=1.d30
c     ifg=1
      end
************************************************************************
      subroutine prihes(n,w1,w0,adp,aup)
c
c     hesian matrix for gaussian smoothness prior: 3:3(2)
c
c  inputs:
c        np: number of points within a rectangular
c      npex: number of points including the minimum # of mirror image.
c      ndelo: number of the delaunay triangulars.
c        tx: horizontal length of the rectangular.
c        ty: vertical lengh of the rectangular.
c       ipt: identification of a triangular: i=ipt in the following.
c  ldt(i,.): identification of the verteces of the i-th triangular.
cx ndt(i,.): id of the neighboring triangular of the i-th one.
c    adt(i): area of the i-th triangular
cx   bdt(i): = adt(i)**na
c
c   input
c        n : number of parameters
c        a : parameters
c   output
c        h : hessian matrix
c
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222,mdl=43333,nr=020)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /data/ xx(nda),yy(nda),zz(nda),zmg(nda),
     &              nn,mm,tx,ty,tz, ind(nda)
c     common /swgrad/ w1,wx,wy,ptdet
      common /delink/ nzu(npa),iu(nr,npa)
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
      dimension adp(npa),aup(nr,npa)
      dimension xi(3),yi(3)
c     dimension ldi(3),xi(3),yi(3)
c
c     return
c
c     if(myrank.eq.0) write(6,*) 'npex',npex,ndelo,w1
      wx = w1
      wy = w1
c     write(6,1) w1
    1 format(1h ,'#3: w1 =',d15.8)
c
*     do 10 ia=1,npex
*     do 10 ib=1,npex
*  10   h(ia,ib)=0.0
*     do 10 i=1,n
      do 10 i=1,npa
        adp(i)=0.0
        do 10 j=1,nr
          aup(j,i)=0.0
   10 continue
c
      do 20 ipt=1,ndelo
        do 30 i=1,3
c         ldi(i)=ldt(ipt,i)
          xi(i)=xx(ldt(ipt,i))
          yi(i)=yy(ldt(ipt,i))
c         write(6,*) ipt,i,ldt(ipt,i),xi(i),yi(i)
   30     continue
        dn0 = xi(1)*(yi(2)-yi(3)) + xi(2)*(yi(3)-yi(1)) +
     &        xi(3)*(yi(1)-yi(2))
c     if(myrank.eq.0) write(6,*) 'dn0',ipt,dn0,(xi(i),i=1,3)
        if(dn0.eq.0.0) go to 20
        do 40 i=1,3
        do 40 j=1,3
          inx=i+1
          if(inx.eq.4) inx=1
          jnx=j+1
          if(jnx.eq.4) jnx=1
          inx2=inx+1
          if(inx2.eq.4) inx2=1
          jnx2=jnx+1
          if(jnx2.eq.4) jnx2=1
          il=ldt(ipt,i)
          jl=ldt(ipt,j)
*         h(il,jl) = h(il,jl) +
*    &      adt(ipt)/dn0**2 * (
*    &      wx*(yi(inx)-yi(inx2))*(yi(jnx)-yi(jnx2)) +
*    &      wy*(xi(inx)-xi(inx2))*(xi(jnx)-xi(jnx2))   )
          if(jl.lt.il) go to 50
          if(jl.eq.il) then
            adp(il)=adp(il) + adt(ipt)/dn0**2 * (
     &      wx*(yi(inx)-yi(inx2))*(yi(jnx)-yi(jnx2)) +
     &      wy*(xi(inx)-xi(inx2))*(xi(jnx)-xi(jnx2))  )
c           if(il.eq.1) write(6,*) adp(il),dn0,ipt,adt(ipt)
          else
            do 60 k=1,nzu(il)
              if(iu(k,il).eq.jl) then
                aup(k,il)=aup(k,il) + adt(ipt)/dn0**2 * (
     &          wx*(yi(inx)-yi(inx2))*(yi(jnx)-yi(jnx2)) +
     &          wy*(xi(inx)-xi(inx2))*(xi(jnx)-xi(jnx2))  )
              endif
   60       continue
          endif
   50   continue
   40   continue
c
   20 continue
cyo      do 70 i=1,npa
      do 70 i=nn+1,npa
c       fp=fp+(a(i)*w0)**2
c       g(i)=g(i)+2*a(i)*w0**2
        adp(i)=adp(i)+2*w0**2
   70 continue
c     if(myrank.eq.0) then
c     do 1111 i=1,n
c       write(6,*) i,adp(i),(aup(k,i),k=1,nr)
c1111 continue
c 111 format(17d14.6)
c     end if
c     stop 1111
      return
      end
c*********************************************************************
      subroutine penalt(w1,w0,a,fp,g)
c     subroutine penalt(n,w1,a,fp,g)
c
c     penalty function for gaussian smoothness prior
c
c  inputs:
c        np: number of points within a rectangular
c      npex: number of points including the minimum # of mirror image.
c     ndelo: number of the delaunay triangulars.
c        tx: horizontal length of the rectangular.
c        ty: vertical lengh of the rectangular.
c       ipt: identification of a triangular: i=ipt in the following.
c  ldt(i,.): identification of the verteces of the i-th triangular.
c    adt(i): area of the i-th triangular
cx   bdt(i): = adt(i)**na
c
c   input
c        n : number of parameters
c        a : parameters
c   output
c       fp : penalty value
c
      implicit real * 8 (a-h,o-z)
      parameter(nda=22222,npa=22222,mdl=43333,nr=020)
      common /triang/ adt(mdl),ldt(mdl,3),ndt(mdl,3),np,npext,ndelo
      common /data/ xx(nda),yy(nda),zz(nda),zmg(nda),
     &               nn,mm,tx,ty,tz, ind(nda)
c     common /swgrad/ w1,wx,wy,ptdet,adp(npa),aup(nr,npa)
      dimension a(npa),g(npa),gx(npa),gy(npa)
      dimension xi(3),yi(3)
c     dimension ldi(3),xi(3),yi(3)
c     data jj /0/
c
      wx = w1
      wy = w1
c     write(6,1) w1
    1 format(1h ,'#3: w1 =',d15.8)
c
      fp=0.0
      do 10 i=1,npa
      g(i)=0.0
   10 continue
      do 20 ipt=1,ndelo
        fpx=0.0
        fpy=0.0
!       do 25 i=1,npa
!         gx(i)=0.0
!         gy(i)=0.0
!  25   continue
        do 30 i=1,3
c         ldi(i)=ldt(ipt,i)
          xi(i)=xx(ldt(ipt,i))
          yi(i)=yy(ldt(ipt,i))
   30   continue
        dn0 = xi(1)*(yi(2)-yi(3)) + xi(2)*(yi(3)-yi(1)) +
     &        xi(3)*(yi(1)-yi(2))
        if(dn0.eq.0.0) go to 20
c
        i11=ldt(ipt,1)
        i22=ldt(ipt,2)
        i33=ldt(ipt,3)
        gx(i11)=0
        gy(i11)=0
        gx(i22)=0
        gy(i22)=0
        gx(i33)=0
        gy(i33)=0
c
        do 40 i=1,3
          inx=i+1
          if(inx.eq.4) inx=1
          inx2=inx+1
          if(inx2.eq.4) inx2=1
          il=ldt(ipt,i)
          fpx = fpx + a(il)*(yi(inx)-yi(inx2))/dn0
          fpy = fpy + a(il)*(xi(inx)-xi(inx2))/dn0
          gx(il) = gx(il) + (yi(inx)-yi(inx2))/dn0
          gy(il) = gy(il) + (xi(inx)-xi(inx2))/dn0
   40   continue
        fp = fp + adt(ipt) * (wx*fpx**2+wy*fpy**2) /2
        do 45 i=1,3
          il=ldt(ipt,i)
          g(il) = g(il) + adt(ipt) * ( wx*fpx*gx(il)+wy*fpy*gy(il) )
   45   continue
c       write(6,*) 'penalt',ipt,adt(ipt),fpx,fpy,wx,wy
   20 continue
cyo      do 50 i=1,npa
      do 50 i=nn+1,npa
        fp=fp+(a(i)*w0)**2
        g(i)=g(i)+2*a(i)*w0**2
   50 continue
      return
      end
************************************************************************
c I'm sending a new simplex program.  I have changed the subroutines
c 'simplx' and 'first'.  the interface between main and subroutine
c simplexwas not changed. in order to restart, after calling the funct,
c the program saves all the values to a file named 'simplex.root'.
c at first of calculation, this program trys to read the file
c 'simplex.root'.  if you want to start the iteration initially,
c please remove the file.
c y.murata thu jun  3 15:48:35 jst 1993
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine simplx(xinit, n, funct, dist, eps, fr)
c -------------------------------------------------------------------- c
c -- parameter: ------------------------------------------------------ c
c ---- i/o xinit:      vector of initial and final x ----------------- c
c ---- i   n:          dimension of the vector x --------------------- c
c ---- i   funct:      minimize function ----------------------------- c
c ---- i   dist:       size of first simplex ------------------------- c
c ---- i   eps:        epsilon --------------------------------------- c
c ---- o   fr:         final value ----------------------------------- c
c -------------------------------------------------------------------- c
c -- this subroutine call next subroutine and function. -------------  c
c --------- first, minmax, newsim, center, reduce, update, epslon ---- c
c -------------------------------------------------------------------- c
      parameter   (maxh=9, maxh5=maxh+10)
c
      external   funct
c
      integer    n
      real*8     xinit(maxh), dist, eps, fr
c
      integer    ipoint, iter, xh, xs, xl, x0, xr, xe, xc, i
      real*8     f(maxh5), x(maxh5, maxh)
      real*8     epslon
      real*8     eps1, alpha, beta, gamma
      common /cont/ init
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
c
      data       alpha/-1.0d0/, beta/0.5d0/, gamma/2.0d0/
c -------------------------------------------------------------------- c
      ipoint=0
cyo
      init=0
c      init=1
      if(init.ne.0) then
c     open(1, file='simplex.root', form='unformatted')
      open(1, file='simplex.root')
*     open(1)
      if(myrank.eq.0) write(6, *) 'read start'
!     if(myrank.eq.0) write(9, *) 'read start'
c     read(1, end=5) ipoint, iter, xh, xs, xl, x0, xr, xe, xc, i,
      read(1,100, end=5) ipoint, iter, xh, xs, xl, x0, xr, xe, xc, i,
     &               xinit, dist, f, x, eps1
  100 format(10i8/(3d25.15))
      if(myrank.eq.0) write(6, '(10(a, i10))')
     &  'start point = ', ipoint, ' iter = ', iter, ' first = ', i
!     if(myrank.eq.0) write(9, '(10(a, i10))')
!    &  'start point = ', ipoint, ' iter = ', iter, ' first = ', i
      close(1)
      end if
c
      if ( ipoint .eq. 11 ) go to 11
      if ( ipoint .eq. 12 ) go to 12
      if ( ipoint .eq. 13 ) go to 13
      if ( ipoint .eq. 14 ) go to 14
      if ( ipoint .eq. 15 ) go to 15
c
5     continue
      x0 = n + 2
      xr = n + 3
      xe = n + 4
      xc = n + 5
      xh = 0
      xs = 0
      xl = 0
c                                            ****   first simplex   ****
      i = 0
4     i = i + 1
      if ( i .gt. n+1 ) go to 3
      call first(i, n, f, x, funct, xinit, dist)
c     open(1, file='simplex.root', form='unformatted')
      open(1, file='simplex.root')
*     open(1)
c     write(1) 11, 0, xh, xs, xl, x0, xr, xe, xc, i,
      write(1,100) 11, 0, xh, xs, xl, x0, xr, xe, xc, i,
     &         xinit, dist, f, x, eps1
      close(1)
11    continue
      if ( dist .eq. 0.0d0 ) return
      go to 4
3     continue
c -------------------------------------------------------------------- c
c                                                ****   main loop   ****
      iter = 0
1     continue
c                                           ****   min max search   ****
      call minmax(n, f, xh, xs, xl)
c                                            ****   epsilon check   ****
      eps1 = epslon(n, f)
      if(myrank.eq.0) then
      write(6, '(31h ####  iteration, f, epsilon = , i5, 2d16.8)')
     &          iter, f(xl), eps1
      write(6, '(5h x = , 4d16.8)') (x(xl, i), i=1, n)
!     write(9, '(31h ####  iteration, f, epsilon = , i5, 2d16.8)')
!    &          iter, f(xl), eps1
!     write(9, '(5h x = , 4d16.8)') (x(xl, i), i=1, n)
      end if
c
      if ( iter .gt. 2 .and. eps1 .lt. eps ) go to 2
c                                               ****   center set   ****
      call center(n, x, xh, x0)
c                                               ****   reflection   ****
      call newsim(n, f, x, funct, x0, xh, alpha, xr)
*     open(1, file='simplex.root', form='unformatted')
      open(1, file='simplex.root')
*     open(1)
c     write(1) 12, iter, xh, xs, xl, x0, xr, xe, xc, i,
      write(1,100) 12, iter, xh, xs, xl, x0, xr, xe, xc, i,
     &         xinit, dist, f, x, eps1
      close(1)
12    continue
c
      if ( f(xr) .gt. f(xs) ) go to 6
c
      if ( f(xr) .ge. f(xl) ) then
        call update(n, f, x, xr, xh)
        go to 7
        end if
c                                                ****   expansion   ****
      call newsim(n, f, x, funct, x0, xr, gamma, xe)
c     open(1, file='simplex.root', form='unformatted')
      open(1, file='simplex.root')
c     open(1)
c     write(1) 13, iter, xh, xs, xl, x0, xr, xe, xc, i,
      write(1,100) 13, iter, xh, xs, xl, x0, xr, xe, xc, i,
     &         xinit, dist, f, x, eps1
      close(1)
13    continue
c
      if ( f(xe) .lt. f(xr) ) then
        call update(n, f, x, xe, xh)
      else
        call update(n, f, x, xr, xh)
        end if
      go to 7
c
6     continue
c
      if ( f(xr) .lt. f(xh) ) call update(n, f, x, xr, xh)
c                                              ****   contraction   ****
      call newsim(n, f, x, funct, x0, xh, beta, xc)
c     open(1, file='simplex.root', form='unformatted')
      open(1, file='simplex.root')
c     open(1)
c     write(1) 14, iter, xh, xs, xl, x0, xr, xe, xc, i,
      write(1,100) 14, iter, xh, xs, xl, x0, xr, xe, xc, i,
     &         xinit, dist, f, x, eps1
      close(1)
14    continue
c
      if ( f(xc) .lt. f(xh) ) then
        call update(n, f, x, xc, xh)
        go to 15
        end if
c                                                   ****   reduce   ****
      call reduce(n, f, x, funct, xl)
c     open(1, file='simplex.root', form='unformatted')
      open(1, file='simplex.root')
c     open(1)
c     write(1) 15, iter, xh, xs, xl, x0, xr, xe, xc, i,
      write(1,100) 15, iter, xh, xs, xl, x0, xr, xe, xc, i,
     &         xinit, dist, f, x, eps1
      close(1)
15    continue
c
7     iter = iter + 1
c
c     if ( iter .gt. 20 ) go to 2
      if ( iter .gt. 1000 ) go to 2
c
      go to 1
c -------------------------------------------------------------------- c
2     continue
c
      do 10 i = 1, n
        xinit(i) = x(xl, i)
10      continue
      fr = f(xl)
c -------------------------------------------------------------------- c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine minmax(n, f, xh, xs, xl)
c
      parameter   (maxh=9, maxh5=maxh+10)
c
      integer    n, xh, xs, xl
      real*8     f(maxh5)
c
      integer    i
      real*8     fmin1, fmax1, fmax2
c -------------------------------------------------------------------- c
      if ( f(1) .gt. f(2) ) then
        fmax1 = f(1)
        xh = 1
        fmax2 = f(2)
        xs = 2
        fmin1 = f(2)
        xl = 2
      else
        fmax1 = f(2)
        xh = 2
        fmax2 = f(1)
        xs = 1
        fmin1 = f(1)
        xl = 1
        end if
c
      do 10 i = 3, n+1
c
        if ( f(i) .gt. fmax1 ) then
          fmax2 = fmax1
          xs = xh
          fmax1 = f(i)
          xh = i
c
        else if ( f(i) .gt. fmax2 ) then
          fmax2 = f(i)
          xs = i
          end if
c
        if ( f(i) .lt. fmin1 ) then
          fmin1 = f(i)
          xl = i
          end if
c
10      continue
c -------------------------------------------------------------------- c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine first(i, n, f, x, funct, xinit, dist)
c
      parameter   (maxh=9, maxh5=maxh+10)
c
      external   funct
c
      integer    i, n
      real*8     f(maxh5), x(maxh5, maxh), xinit(maxh), dist
c
      integer    j
      real*8     xx(maxh)
c -------------------------------------------------------------------- c
      do 10 j = 1, n
        x(i, j) = xinit(j)
        xx(j) = x(i, j)
10      continue
c
      if ( i .ne. 1 ) then
        x(i, i-1) = x(i, i-1) + dist
        xx(i-1) = x(i, i-1)
        end if
c
      call funct(n, xx, f(i))
c -------------------------------------------------------------------- c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine center(n, x, xh, x0)
c
      parameter   (maxh=9, maxh5=maxh+10)
c
      integer    n, xh, x0
      real*8     x(maxh5, maxh)
c
      integer    i, j
c -------------------------------------------------------------------- c
      do 10 j = 1, n
        x(x0, j) = 0.0d0
        do 20 i = 1, n+1
          if ( i .ne. xh ) x(x0, j) = x(x0, j) + x(i, j)
20        continue
        x(x0, j) = x(x0, j) / dfloat(n)
10      continue
c -------------------------------------------------------------------- c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      function   epslon(n, f)
c
      parameter   (maxh=9, maxh5=maxh+10)
c
      integer    n
      real*8     f(maxh5)
c
      integer    i
      real*8     epslon, emean
c -------------------------------------------------------------------- c
      emean = 0.0d0
      do 10 i = 1, n+1
        emean = emean + f(i)
10      continue
      emean = emean / dfloat(n+1)
c
      epslon = 0.0d0
      do 20 i = 1, n+1
        epslon = epslon + (f(i) - emean) ** 2
20      continue
      epslon = dsqrt(epslon) / dfloat(n+1)
c -------------------------------------------------------------------- c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine newsim(n, f, x, funct, x0, in, prm, out)
c
      parameter   (maxh=9, maxh5=maxh+10)
c
      external   funct
c
      integer    n, x0, in, out
      real*8     f(maxh5), x(maxh5,maxh), xx(maxh), prm
c
      integer    j
c -------------------------------------------------------------------- c
      do 10 j = 1, n
        x(out, j) = prm * x(in, j) + (1.0d0 - prm) * x(x0, j)
        xx(j) = x(out,j)
10      continue
      call funct(n, xx, f(out))
c -------------------------------------------------------------------- c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine update(n, f, x, in, out)
c
      parameter   (maxh=9, maxh5=maxh+10)
c
      integer    n, in, out
      real*8     f(maxh5), x(maxh5, maxh)
c
      integer    j
c -------------------------------------------------------------------- c
      do 10 j = 1 , n
        x(out, j) = x(in, j)
10      continue
c
      f(out) = f(in)
c -------------------------------------------------------------------- c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine reduce(n, f, x, funct, xl)
c
      parameter   (maxh=9, maxh5=maxh+10)
c
      external   funct
c
      integer    n, xl
      real*8     f(maxh5), x(maxh5, maxh), xx(maxh)
c
      integer    i, j
c -------------------------------------------------------------------- c
      do 10 i = 1, n+1
c
        if ( i .ne. xl ) then
c
          do 20 j = 1, n
            x(i, j) = (x(i, j) + x(xl, j)) / 2.0d0
            xx(j) = x(i, j)
20          continue
          call funct(n, xx, f(i))
c
          end if
c
10      continue
c -------------------------------------------------------------------- c
      return
      end
      subroutine r6(jx,nox,jy,noy,ram,rrr)
      implicit real * 8 (a-h,o-z)
      if(jx.eq.1 .or. jx.eq.nox) ram=rrr/2
      if(jy.eq.1 .or. jy.eq.noy) ram=rrr/2
      if(jx.eq.1 .and. jy.eq.1) ram=rrr/3
      if(jx.eq.nox .and. jy.eq.noy) ram=rrr/3
      if(jx.eq.1 .and. jy.eq.noy) ram=rrr/6
      if(jx.eq.nox .and. jy.eq.1) ram=rrr/6
      return
      end
************************************************************************
      subroutine iccg (nd, n, nr, nzu, nzl, iu, il,
     $           au, al, ad, b, eps, x,
     $           uu, ul, ud, r, p, ru, ap, iv)
******************************************************
*                                                    *
*     solution of system of linear equations         *
*             a * x = b  (a : symmetric)             *
*     by incomplete cholesky decomposition and       *
*          conjugate gradient (iccg) method          *
*                                                    *
*     ==== input ====                                *
*                                                    *
*     n...order of original matrix a                 *
*     nr...maximum number of nonzero elements in one *
*          row and in one column in upper-half       *
*          off-diagonal part                         *
*                                                    *
*     the elements of a must be given as follows :   *
*                                                    *
*     nzu(i)...number of nonzero off-diagonal        *
*              elements in the i-th row in upper-    *
*              half part of a                        *
*     iu(j,i)...column number of the j-th nonzero    *
*               element in the i-th row              *
*               (count rightwards from diagonal+1)   *
*     au(j,i)...j-th nonzero element in the i-th row *
*     ad(i)...i-th diagonal element                  *
*                                                    *
*     b(i)...right hand side vector                  *
*     x(i)...initial guess for solution              *
*     eps...absolute tolerance for residual vector   *
*                                                    *
*     ==== output ====                               *
*                                                    *
*     x(i)...solution x                              *
*                                                    *
*     ==== work arrays ====                          *
*                                                    *
*     nzl,il,al are arranged from nzu,iu,au          *
*               by calling crconv                    *
*                                                    *
*     nzl(i)...number of nonzero off-diagonal        *
*              elements in the i-th column in upper- *
*              half part of a                        *
*     il(j,i)...row number of the j-th nonzero       *
*               element in the i-th column           *
*               (count downwards from the 1-st row)  *
*     al(j,i)...j-th nonzero element of i-th column  *
*                                                    *
*     uu,ul and ud are lt,l and d of incompletely    *
*        cholesky decomposed matrix ldlt of a        *
*        obtained by calling decomp and solve        *
*              (lt = transpose of l)                 *
*        icdcmp computes ul and ud                   *
*        crconv arranges uu from ul                  *
*                                                    *
*     uu(j,i)...j-th nonzero element in the i-th row *
*               of decomposed matrix lt of ldlt      *
*     ul(j,i)...j-th nonzero element in the i-th     *
*               column of decomposed matrix lt of    *
*               ldlt  or                             *
*               j-th nonzero element in the i-th row *
*               of decomposed matrix l of ldlt       *
*     ud(i)...i-th diagonal element of decomposed    *
*             matrix d of ldlt                       *
*                                                    *
*     iv...for easy vectorization for compiler       *
*                                                    *
*     r(i)...residual vector (in cg)                 *
*     p(i)...direction vector for search (in cg)     *
*     ru(i)...(ldlt)**(-1)*r (in cg)                 *
*     ap(i)...product a*p (in cg)                    *
*                                                    *
*     ==== subroutines called ====                   *
*                                                    *
*     cg, prodmv, crconv, icdcmp, icsolv             *
*                                                    *
******************************************************
      implicit real*8(a-h,o-z)
      integer nzu(nd),nzl(nd),iu(nr,nd),il(nr,nd)
      real*8 au(nr,nd),al(nr,nd),ad(nd)
      real*8 uu(nr,nd),ul(nr,nd),ud(nd)
      real*8 b(nd),x(nd)
      integer iv(nd)
      real*8 r(nd),p(nd),ru(nd),ap(nd)
      common /diagnl/ idet
      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised
*
c     do i=1,nd
c     nzl(i)=0
c     end do
*
      call crconv (nd, n, nr, nzu, nzl, iu, il, au, al, iv)
      call icdcmp (nd, n, nr, nzl, il, al, ad, ul, ud)
c     write(6,1) ud
c   3 format(9f8.1)
c
      if(idet.eq.1) then
        eps=0.0
        do 10 i=1,n
c         write(6,*) i,ud(i),eps
          if(ud(i).gt.0.0) then
            eps=eps-log(ud(i))
          else
           if(myrank.eq.0)
     &      write(6,*) 'non-pos diag.',i,ud(i),eps
!          if(myrank.eq.0)
!    &      write(9,*) 'non-pos diag.',i,ud(i),eps
          endif
   10   continue
        return
      else
c
      call crconv (nd, n, nr, nzl, nzu, il, iu, ul, uu, iv)
c     write(6,1) ud
c     write(6,3) uu
c     stop 555
*
      call cg (nd, n, nr, nzu, nzl, iu, il,
     $         au, al, ad, b, eps, x,
     $         uu, ul, ud, r, p, ru, ap)
c
      endif
*
      return
      end
************************************************************************
      subroutine crconv (nd, n, nr, nz1, nz2, i1, i2,
     $                   a1, a2, iv)
************************************************
*                                              *
*     convert row information into column      *
*         information or vice versa            *
*                                              *
*     nz1 to nz2,   i1 to i2,   a1 to a2       *
*                                              *
*     when row to column :                     *
*     nzu to nzl,   iu to il,   au to al       *
*                               uu to ul       *
*     when column to row :                     *
*     nzl to nzu,   il to iu,   al to au       *
*                               ul to uu       *
*                                              *
*     ==== work array ====                     *
*                                              *
*     iv...for easy vectorization for compiler *
*                                              *
************************************************
      implicit real*8(a-h,o-z)
      integer nz1(nd),nz2(nd),i1(nr,nd),i2(nr,nd)
      real*8 a1(nr,nd),a2(nr,nd)
      integer iv(nd)
*
c     do 540 i = 1, n
      do 540 i=1,n+1
        nz2(i) = 0
        do 550 j = 1, nr
          i2(j,i) = 0
          a2(j,i) = 0.0
  550   continue
  540 continue
*
      do 10 i = 1, n
        do 520 nu = 1, nz1(i)
          iv(i1(nu,i)) = nz2(i1(nu,i)) + 1
  520   continue
        do 530 nu = 1, nz1(i)
          nz2(i1(nu,i)) = iv(i1(nu,i))
          i2(iv(i1(nu,i)),i1(nu,i)) = i
          a2(iv(i1(nu,i)),i1(nu,i)) = a1(nu,i)
  530   continue
   10 continue
*
      return
      end
************************************************************************
      subroutine icdcmp (nd, n, nr, nzl, il, al, ad,
     $                   ul, ud)
******************************************************
*                                                    *
*       incomplete cholesky decomposition of a       *
*             into ldlt in iccg method               *
*              (lt = transpose of l)                 *
*                                                    *
*     relation between l(k,k) and d(k) is            *
*           l(k,k) * d(k) = 1   (iccg version)       *
*                                                    *
******************************************************
      implicit real*8(a-h,o-z)
      integer nzl(nd),il(nr,nd)
      real*8 al(nr,nd),ad(nd),ul(nr,nd),ud(nd)
*
c     real*8 t,eps
      real*8 t
*
c     call maceps (eps)
      eps=1.0d-6
*
      do 10 k = 1, n
*
        do 20 mu = 1, nzl(k)
*
          i = il(mu,k)
          mui = 1
          muk = 1
          ul(mu,k) = al(mu,k)
*
    1     continue
*
          if (il(mui,i) .gt. il(muk,k)) then
            muk = muk + 1
          else if (il(mui,i) .lt. il(muk,k)) then
            mui = mui + 1
          else
            ul(mu,k) = ul(mu,k)
     $        - ud(il(muk,k)) * ul(mui,i) * ul(muk,k)
            mui = mui + 1
            muk = muk + 1
          end if
*
c         if ((il(mui,i) .le. i-1) .and.
c    $        (il(muk,k) .le. i-1)) go to 1
          if ((mui .le. nzl(i)) .and.
     $        (muk .le. nzl(k))) go to 1
*
   20   continue
*
        t = ad(k)
        do 530 mu = 1, nzl(k)
          t = t - ud(il(mu,k)) * dble(ul(mu,k))**2
  530   continue
        if (abs(t) .le. eps) t = eps
*
        ud(k) = 1 / t
*
   10 continue
*
      return
      end
************************************************************************
      subroutine maceps (epsmac)
***************************************************
*                                                 *
*     get machine epsilon in single precision     *
*                                                 *
***************************************************
      implicit real*8(a-h,o-z)
      epsmac = 1.0
  100 continue
      if (1.0 + epsmac .le. 1.0) then
        epsmac = 2 * epsmac
        return
      end if
      epsmac = epsmac * 0.5
      go to 100
*
      end
************************************************************************
      subroutine cg (nd, n, nr, nzu, nzl, iu, il,
     $               au, al, ad, b, eps, x,
     $               uu, ul, ud, r, p, ru, ap)
***************************************************
*                                                 *
*       conjugate gradient method for iccg        *
*                                                 *
***************************************************
      implicit real*8(a-h,o-z)
      integer nzu(nd),nzl(nd),iu(nr,nd),il(nr,nd)
      real*8 au(nr,nd),al(nr,nd),ad(nd)
      real*8 uu(nr,nd),ul(nr,nd),ud(nd)
      real*8 b(nd),x(nd),xxmin(5*22222)
c                   in the above 22222=nd
      real*8 r(nd),p(nd),ru(nd),ap(nd)
      real*8 rur0,rur1,pap,res2
*
      kend = n + 10
c     res22b2=1.0d77
      resdbmin=1.0d77
      eps=1.0d-6
      eps2 = eps**2
*
*     ---- initialization ----
*
      call prodmv (nd, n, nr, nzu, nzl, iu, il, au, al,
     $         ad, x, ap)
*
      do 510 i = 1, n
        r(i) = b(i) - ap(i)
  510 continue
*
      call icsolv (nd, n, nr, nzu, nzl, iu, il,
     $             uu, ul, ud, r, ru)
*
      rur0 = 0.0d0
      do 520 i = 1, n
        p(i) = ru(i)
        rur0 = rur0 + r(i) * dble(ru(i))
  520 continue
*
      b2 = 0.0
      do 530 i = 1, n
        b2 = b2 + b(i)**2
  530 continue
*
*     ---- iterative loop ----
*
      do 10 k = 1, kend
*
        call prodmv (nd, n, nr, nzu, nzl, iu, il, au, al,
     $           ad, p, ap)
*
        pap = 0.0d0
        do 540 i = 1, n
          pap = pap + p(i) * dble(ap(i))
  540   continue
*
        alpha = rur0 / pap
*
        do 550 i = 1, n
          x(i) = x(i) + alpha * p(i)
          r(i) = r(i) - alpha * ap(i)
  550   continue
c     write(6,*) kend,k,x(1),x(n)
*
        call icsolv (nd, n, nr, nzu, nzl, iu, il,
     $               uu, ul, ud, r, ru)
*
        rur1 = 0.0d0
        res2 = 0.0d0
*
        do 560 i = 1, n
          rur1 = rur1 + r(i) * dble(ru(i))
          res2 = res2 + r(i)**2
  560   continue
c
c       if(k.ge.10 .and. res2/b2.gt.res22b2) then
c         write(6,*) 'cgres_1',k,res2,res2/b2
c         return
c       else
c         res22b2=res2/b2
c       endif
c
        if(res2/b2.lt.resdbmin) then
          kmin=k
          resdbmin=res2/b2
          do 561 i=1,n
            xxmin(i)=x(i)
  561     continue
        endif
        if(k.ge.kend) then
          do 562 i=1,n
            x(i)=xxmin(i)
  562     continue
          write(6,*) 'cgres_1',kmin,res2,res2/b2
!         write(9,*) 'cgres_1',kmin,res2,res2/b2
          return
        endif
*
*       write(6,*) k,res2/b2,res2,b2,eps2
c       if (res2/b2 .le. eps2) then
        if (res2.le.eps2**2 .or. res2/b2.le.eps2) then
c         write(6,*) 'cgres_0',k,res2,res2/b2,b2
c         stop 66666
*
*         ---- converged ----
*
          return
        end if
*
*       ---- not converged ----
*
        beta = rur1 / rur0
        rur0 = rur1
*
        do 570 i = 1, n
          p(i) = ru(i) + beta * p(i)
  570   continue
*
   10 continue
*
      return
      end
************************************************************************
      subroutine prodmv (nd, n, nr, nzu, nzl, iu, il,
     $               au, al, ad, x, ap)
***************************************************
*                                                 *
*               compute product                   *
*          ap = (matrix a) * (vector x)           *
*                                                 *
***************************************************
      implicit real*8(a-h,o-z)
      integer nzu(nd),nzl(nd),iu(nr,nd),il(nr,nd)
      real*8 au(nr,nd),al(nr,nd),ad(nd)
      real*8 x(nd),ap(nd)
*
      real*8 apv
*
      do 10 i = 1, n
*
        apv = 0.0d0
        do 520 mu = 1, nzl(i)
          apv = apv + al(mu,i) * dble(x(il(mu,i)))
  520   continue
        apv = apv + ad(i) * dble(x(i))
        do 530 nu = 1, nzu(i)
          apv = apv + au(nu,i) * dble(x(iu(nu,i)))
  530   continue
        ap(i) = apv
*
   10 continue
*
      return
      end
************************************************************************
      subroutine icsolv (nd, n, nr, nzu, nzl, iu, il,
     $                   uu, ul, ud, b, x)
***************************************************
*                                                 *
*     solve system of linear equations            *
*              (l * d * lt) * x = b               *
*     in iccg method (lt = transpose of l)        *
*     ldlt is incomplete cholesky decomposition   *
*             arranged by icdcmp                  *
*                                                 *
***************************************************
      implicit real*8(a-h,o-z)
      integer nzu(nd),nzl(nd),iu(nr,nd),il(nr,nd)
      real*8 uu(nr,nd),ul(nr,nd),ud(nd),b(nd),x(nd)
*
      real*8 t
*
*     ---- forward substitution ----
*
      do 10 k = 1, n
        t = b(k)
        do 520 mu = 1, nzl(k)
          t = t - ul(mu,k) * dble(x(il(mu,k)))
  520   continue
        x(k) = ud(k) * t
   10 continue
*
*     ---- backward substitution ----
*
      do 30 k = n, 1, -1
        t = 0.0
        do 540 nu = 1, nzu(k)
          t = t + uu(nu,k) * dble(x(iu(nu,k)))
  540   continue
        x(k) = x(k) - ud(k) * t
   30 continue
*
      return
      end


c********************************************************************
      subroutine para(nk)

      common/mpi/nprocs, myrank, ista, iend, ista2, iend2,ised            


c     write(*,*)'para nn', nk,myrank
      if(myrank.eq.0)ista=1
      if(myrank.ne.0)ista=sqrt(float(nk)*float(nk)*
     &         float(myrank)/float(nprocs))+1 
      if(myrank.ne.nprocs-1)iend=sqrt(float(nk)*float(nk)*
     &    float(myrank+1)/float(nprocs))
      if(myrank.eq.nprocs-1)iend=nk

      ista2=nk*myrank/nprocs+1
      iend2=nk*(myrank+1)/nprocs      
  
c     print*,'in process',myrank,'data:', ista, iend, 'of ',nk
c     print*,'in process',myrank,'data:', ista2, iend2, 'of ',nk


      return
      end


c**********************************************************************
