cprocess subchk,argchk
c  Modified from delo2d.iccgap.f: we "read(5," iccgap.in or use deloP.sh.
c  copied from #ogataxx.delone.fort at 97/05/07.
c  simplified penalty only consider to suppress the gradient of
c  function on the 2-dimensional delonay trianglar.
c  data set is made by @tane20p.cntl and @plptpat.fort
      program main
      implicit real * 8 (a-h,o-z)
      parameter(nda=18888,npa=18888,mdl=22222)
      common /data/ xx(nda),yy(nda),zz(nda,5),sd(nda),nn,mp(nda)
      common /abicc/abic,dabic
      common /xad/xd(nda),m
      common /files/ indata,outfile,an
      common /triang/ tx,ty,adt(mdl),bdt(mdl),cdt(mdl),idt(nda),
     &                ldt(mdl,3),ndt(mdl,3),np,npex,nd
      common     / ccc /  isw , ipr
      character*72 indata,outfile
      dimension x(30)
      external funcs,funcsb
c
c     open(1,file='deloPoi.in')
c     read(5,111) indata
c     read(5,*) an
c     read(5,111) outfile
c     read(1,111) indata
c     read(1,*) an
c     read(1,111) outfile
c     close(1)
c     write(6,*) indata
c     write(6,*) an
c     write(6,*) outfile
      indata='delone2.out'
      outfile='delo2d-poisson.omap'
      an=1
  111 format(a)
      open(4,file='delo2d-poisson.conf')
      read(4,*) xmin,ymin,threshmag
      write(6,*) 'xmin,ymin,threshmag=',xmin,ymin,threshmag
      read(4,*) w1
      write(6,*) 'weight=', w1
      read(4,*) ipr
      write(6,*) 'linear ipr',ipr
      close(4)
      call redata(xx,yy,mp,zz,sd,nn,dens,xmin,ymin,threshmag)
c
      do 10 i=1,nda
        xd(i)=log(dens)
   10 continue
      xd(1)=0.01
      x(1)=sqrt(w1)
c     x(2)=sqrt(0.10000000d+2)
c\m murata's routine follows: two more lines c\m should be changed.
      dist=1.0d0
      eps=1.d-3
      call simplx(x, 1, funcsb, dist, eps, xm)
c     call simplx(x, 6, funcsb, dist, eps, xm)
c
c     open(12)
c     read(12,6) w1,w2,abic,upd
c     read(12,6) (xd(i),i=1,95)
c     close(12)
c   6 format(4e20.12)
c     x(1)=sqrt(w1)
c     x(2)=sqrt(w2)
c
c     write(6,1) (xd(i),i=1,m)
    1 format(1h ,5f10.3/(1h ,5f10.3))
      x(1)=x(1)**2
c     x(2)=x(2)**2
c10   x(1)=exp(x(1))
c10   x(2)=exp(x(2))
      write(6,4) (x(i),i=1,1), abic,m
*     open(11,file='work')
      open(11,file=outfile)
c     open(11)
      write(11,4) (x(i),i=1,1), abic,m
      write(11,2) (xd(i),i=1,m)
    2 format(5e15.5)
    4 format(2e15.5,i7)
      close(11)
c
      stop
      end
c*********************************************************************
      subroutine redata(xx,yy,mp,zz,sd,nn,dens,xmin,ymin,threshmag)
c  inputs:
c         np: number of points within a rectangular
c       npex: number of points including the minimum # of mirror image.
c         nd: number of the delaunay triangulars.
c         tx: horizontal length of the rectangular.
c         ty: vertical lengh of the rectangular.
c       dens: density=np/(tx*ty)
c        ipt: identification of a triangular: i=ipt in the following.
c   ldt(i,.): identification of the verteces of the i-th triangular.
c   ndt(i,.): id of the neighboring triangular of the i-th one.
c     adt(i): area of the i-th triangular
c
      implicit real * 8 (a-h,o-z)
      parameter(nda=18888,npa=18888,mdl=22222,nr=016)
      common /triang/ tx,ty,adt(mdl),bdt(mdl),cdt(mdl),idt(nda),
     &                ldt(mdl,3),ndt(mdl,3),np,npex,nd
      common /delink/ nzu(npa),iu(nr,npa)
      dimension xx(nda),yy(nda),zz(nda,5),sd(nda),mp(nda)
      common /files/ indata,outfile,an
      character*72 indata,outfile
c
c reading configulation of the delaunay tassilation
c     read(5,*) np,npex,nd,tx,ty,dens
*     read(5,*) np,npex,nd,tx,ty
*     read(5,*) (idum,xx(ipt),yy(ipt),ipt=1,npex)
      open(2,file=indata)
      read(2,*) np,npex,nd,tx,ty
      dens=npex/tx/ty
      yd=ymin+ty/2
c     read(2,*) (idum,xx(ipt),yy(ipt),ipt=1,npex)
c2002/02/07
      do 100 ipt=1,npex
      read(2,*) idum,xx(ipt),yy(ipt)
      xx(ipt)=xx(ipt)*cos(yd*3.14159/180.)
  100 continue
      tx=tx*cos(yd*3.14159/180.)
c2002/02/07
      write(6,*) np,npex,nd,tx,ty,dens
c     read(5,*) nn,np,npex,nd,tx,ty,dens
c     read(5,*) (xx(i),yy(i),mp(i),(zz(i,k),k=1,mp(i)), i=1,npex)
c     read(5,*) (i,(ldt(i,j),j=1,3),(ndt(i,j),j=1,3),adt(i), i=1,nd)
*     read(5,*) (idum,(ldt(i,j),j=1,3),adt(i), i=1,nd)
      read(2,*) (idum,(ldt(i,j),j=1,3),adt(i), i=1,nd)
      do 510 i = 1, npex
*       read(5,*) idummy, nzu(i), (iu(j,i), j=1,nzu(i))
        read(2,*) idummy, nzu(i), (iu(j,i), j=1,nzu(i))
  510 continue
c
*     an=0
*     an=1
      write(6,*) 'an =',an
      do 1111 i=1,nd
c     adt(i)=adt(i)*cos(38.5*3.14159/180.)
      adt(i)=adt(i)*cos(yd*3.14159/180.)
      bdt(i)=adt(i)**an
      cdt(i)=log(adt(i))
 1111 continue
c
      nn=np
c     open(10, file='inoue.dat')
c     read(10,7)
c     do 10 i=1,nn
c  10   read(10,*) xx(i),yy(i),zz(i),sd(i)
c   7 format(4f10.4)
c     close(10)
c
c     write(6,5) (xx(i),yy(i),i=1,nn)
    5 format(///1h ,9f13.5/(1h ,9f13.5))
c     do 20 i=1,nn
c     sd(i)=1.d0
c  20 continue
      write(6,4) tx,ty,nn,np,npex,nd
    4 format(1h ,'tx,ty ',2f10.3,'     nn,np,npex,nd =',4i7)
c     do 30 i=1,npex
c  30    write(6,3) i,xx(i),yy(i),mp(i),(zz(i,k),k=1,mp(i))
c  30    write(6,6) xx(i),yy(i),zz(i),sd(i)
    6 format(1h ,4f13.5)
    3 format(1h ,i5,2f10.5,i3,5f5.1)
cremember!
c     cutoff=2.85
      do 40 i=1,npex
      if(i.gt.np) go to 60
      do 50 k=1,mp(i)
c     zz(i,k)=zz(i,k)-threshmag
      zz(i,k)=1.d0
c     if(i.le.17) zz(i,k)=-1.d0
   50 continue
   40 continue
   60 continue
      return
      end
c*********************************************************************
      subroutine  davidn( x,n,ihes,funct,ff )
c
c          minimization by davidon-fletcher-powell procedure
c
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             ihes:    =0   inverse of hessian matrix is not available
c                      =1   inverse of hessian matrix is available
c
c          output:
c             x:       vector of minimizing solution
c
      implicit  real * 8  ( a-h , o-z )
      parameter(nda=18888,npa=18888,mdl=22222,nr=016)
      dimension  x(nda) , dx(nda) , g(nda) , g0(nda) , y(nda)
*     dimension  h(nda,nda) , wrk(nda) , s(nda)
      dimension  hd(   nda) , wrk(nda) , s(nda)
*     dimension  h0(nda,nda),w(nda),sneg(nda),ip(nda)
      dimension              w(nda),sneg(nda),ip(nda)
      common /abicc/abic,dabic
*     common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,ph(npa,npa)
      common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,adp(npa),aup(nr,npa)
      common /print/ itp
      common /dvfnhs/ rss,rss1,gh(nda),fpdum
c
      common /delink/ nzu(npa),iu(nr,npa)
      integer nzl(npa),il(nr,npa)
      real*8 au(nr,npa),ad(npa),al(nr,npa)
      real*8 uu(nr,npa),ul(nr,npa),ud(npa)
      real*8 r(npa),p(npa),ru(npa),ap(npa)
      integer iv(npa)
      common /diagnl/ idet
c
      external funct
      data  tau1 , tau2  /  0.1d-2 , 0.1d-2  /
      data  eps1 , eps2  / 0.1d-3 , 0.1d-3  /
c?    nr=50
      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
*     do  20   i=1,n
*     do  10   j=1,n
*  10 h(i,j) = 0.0d00
*     s(i) = 0.0d00
*     dx(i) = 0.0d00
*  20 h(i,i) = 1.0d00
      call hes(n,x,au,ad)
      do  20   i=1,n
        s(i) = 0.0d00
        dx(i) = 0.0d00
        hd(i)=1.d0/ad(i)
   20 continue
c
      call  funct( n,x,xm,g,ig )
c
      if(lpr.ge.2) write( 6,345 ) w1
      if(lpr.ge.2) write( 6,340 ) wx,wxx, xm
c
c          inverse of hessian computation (if available)
c
      if( ihes .eq. 1 )   call  hesian( x,n,h )
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_down
      iter=iter+1
c
      if(mod(iter,001).eq.0) then
*       call hes(n,x,h)
*       call invdet(h,tdet,n,nda)
*       call hes(n,x,h0)
*       call icdata(n,nr,h0,au,ad,y,wrk)
        call hes(n,x,au,ad)
        call iccg (n, n, nr, nzu, nzl, iu, il,
     $           au, al, ad, y, eps1i, wrk,
     $           uu, ul, ud, r, p, ru, ap, iv)
        write(6,*) '#iteration=',iter
*       write(6,*) (y(j),j=1,9)
*       write(6,*) (wrk(j),j=1,9)
c       sumy=0.
c       sumw=0.
c       do 111 j=1,n
c       sumy=sumy+y(j)
c       sumw=sumw+wrk(j)
c 111 continue
c       write(6,*) n,sumy,sumw
c       stop 222
*     end if
      else
c_up
*     do  60   i=1,n
*     sum = 0.0d00
*     do  50   j=1,n
*  50 sum = sum + y(j) * h(i,j)
*  60 wrk(i) = sum
c_down
      do  60   i=1,n
   60 wrk(i) = y(i) * hd(i)
      end if
c_up
   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. s2 )     go to 100
c
c          update the inverse of hessian matrix
c
c     ---  davidon-fletcher-powell type correction  ---
c
*     do  90   i=1,n
*     do  90   j=i,n
*     h(i,j) = h(i,j) + dx(i)*dx(j)/s2 - wrk(i)*wrk(j)/s1
*  90 h(j,i) = h(i,j)
*     go to  120
c
c     ---  fletcher type correction  ---
c
  100 continue
*     stem = s1 / s2 + 1.0d00
*     do  110   i=1,n
*     do  110   j=i,n
*     h(i,j) = h(i,j)- (dx(i)*wrk(j)+wrk(i)*dx(j)-dx(i)*dx(j)*stem)/s2
* 110 h(j,i) = h(i,j)
c
  120 continue
c_down
      if(iter.gt.0 .and. mod(iter,001).eq.0) then
*       call hes(n,x,h)
*       call invdet(h,tdet,n,nda)
c       call hes(n,x,h0)
c       call decomp(nda,n,h0,w,ip,tdet)
*       call solve(nda,n,h0,g,sneg,ip)
*       call icdata(n,nr,h0,au,ad,g,sneg)
c       call hes(n,x,au,ad)
        call iccg (n, n, nr, nzu, nzl, iu, il,
     $           au, al, ad, g, eps, sneg,
     $           uu, ul, ud, r, p, ru, ap, iv)
*       write(6,*) '#iteratio2=',iter
*       write(6,*) (g(j),j=1,5)
*       write(6,*) (sneg(j),j=1,5)
c       sumg=0.
c       sums=0.
c       do 222 j=1,n
c       sumg=sumg+g(j)
c       sums=sums+sneg(j)
c 222 continue
c       write(6,*) n,sumg,sums
c       stop 222
c
        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_up
*     ss = 0.0d00
*     do  150   i=1,n
*     sum = 0.0d00
*     do  140   j=1,n
* 140 sum = sum + h(i,j)*g(j)
*     ss = ss + sum * sum
* 150 s(i) = -sum
c_down
      ss = 0.0d00
      do  150   i=1,n
      sum = hd(i)*g(i)
      ss = ss + sum * sum
  150 s(i) = -sum
      end if
c_up
      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
*     do  190   i=1,n
*     do  180   j=1,n
* 180 h(i,j) = 0.0d00
*     h(i,i) = 1.0d00
* 190 s(i) = -s(i)
      call hes(n,x,au,ad)
      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
      call linear( x,s,ramda,ed,n,ig,funct )
c
      if(lpr.ge.2) write( 6,330 )     ramda , ed , s1 , s2
c     iter=iter+1
c     if(mod(iter,10).eq.0) then
c       call hes(n,x,h)
c       call invdet(h,tdet,n,nda)
c       if(mod(iter,10).eq.0) write(6,*) '#iteration=',iter
c     end if
c   1 format(1h ,16f8.0)
 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( xmb/xm-1.d0 .lt. eps1  .and.  sqrt(s1) .lt. eps2 )  go to 900
c     if( sqrt(s1) .lt. eps2 )  go to 900
11111 continue
      if( icc .ge. 5 )     go to 900
      go to 11110
  900 continue
      ff=xm
      if(lpr.ge.2) write(6,350) ff,rss1
  350 format(1h ,'penalized log likelihood =',d23.15,5x,'rss1 =',d12.5)
      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)
c
c
cremember!
      pi=3.14159265358979d0
cex   ff=xm-ptdet/2+n*log(sqrt(2*pi))
      ff=xm-ptdet/2+(n-1)*log(sqrt(2*pi))
cex   ff=xm-ptdet/2+(n-3)*log(sqrt(2*pi))
*     call hes(n,x,h)
cex   call ltinv(h,tdet,n,nda)
c     call clock
c     call ltinv(h,tdet,n,nda)
c     call decomp(nda,n,h,w,ip,tdet)
c
      idet=1
      call hes(n,x,au,ad)
c     call icdata(n,nr,h,au,ad,y,wrk)
      call iccg (n, n, nr, nzu, nzl, iu, il,
     $           au, al, ad, y, tdet, wrk,
     $           uu, ul, ud, r, p, ru, ap, iv)
      idet=0
c
c     call clock(tin)
c     write(6,*) 'ltinv-time',tin
      pn=tdet
cex   abic=2*ff+pn-n*log(2*pi)
      abic=2*ff+pn-(n-1)*log(2*pi)
      write(6,599) w1
      write(6,598) abic,ff,pn
  599 format(1h ,'#2: w1 =',d15.8)
  598 format(1h ,' abic =',d17.10,x,' -l =',d17.10,x,' 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 ,'wx,wy=',d12.5,'wxx,wyy=',d12.5,3x,'pell =',d23.15)
  345 format( 1h ,'#1: w1 = ',d15.8)
  600 format( 1h ,'-----  x  -----' )
  601 format( 1h ,'***  gradient  ***' )
  610 format( 1h ,10d13.5 )
      end
c*********************************************************************
      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=18888,npa=18888,mdl=22222)
      integer  return,sub
      dimension  x(npa) , h(npa) , x1(npa)
      dimension  g(npa)
      common     / ccc /  isw , ipr
      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
c*********************************************************************
      subroutine hesian( x,n,h )
      return
      end
c*********************************************************************
      subroutine invdet(x,tdet,mm,mj)
c-----------------------------------------------------------------------
c       the inverse and determinant of x computation
c----------------------------------------------------------------------
c       inputs:
c          x:     mm*mm square matrix
c          mm:    dimension of x
c          mj:    absolute dimension of x in the main program
c
c       outputs:
c          x:     inverse of x
c          xdet:  determinant of x
c
      implicit  real * 8  ( a-h , o-z )
      parameter(nda=18888,npa=18888,mdl=22222)
      dimension x(mj,mj)
      dimension  ids(nda)
c
c     return
c
      xdet = 1.0d00
      tdet=0.0d0
      do 10 l=1,mm
c     pivoting at l-th stage
      xmaxp=0.10000d-10
      maxi=0
      do 110 i=l,mm
    1 if( abs(xmaxp) .ge. abs(x(i,l)) )     go to 110
      xmaxp=x(i,l)
      maxi=i
  110 continue
      ids(l)=maxi
      if(maxi.eq.l) go to 120
      if(maxi.gt.0) go to 121
      xdet = 0.0d00
      go to 140
c     row interchange
  121 do 14 j=1,mm
      xc=x(maxi,j)
      x(maxi,j)=x(l,j)
   14 x(l,j)=xc
      xdet=-xdet
  120 tdet=tdet+log(abs(xmaxp))
c 120 xdet=xdet*xmaxp
c 120 continue
      xc = 1.0d00 / xmaxp
      x(l,l)=1.0d00
      do 11 j=1,mm
   11 x(l,j)=x(l,j)*xc
      do 12 i=1,mm
      if(i.eq.l) go to 12
      xc=x(i,l)
      x(i,l) = 0.0d00
      do 13 j=1,mm
   13 x(i,j)=x(i,j)-xc*x(l,j)
   12 continue
   10 continue
      if(mm.gt.1) go to 123
      go to 140
c     column interchange
  123 mm1=mm-1
      do 130 j=1,mm1
      mmj=mm-j
      jj=ids(mmj)
      if(jj.eq.mmj) go to 130
      do 131 i=1,mm
      xc=x(i,jj)
      x(i,jj)=x(i,mmj)
  131 x(i,mmj)=xc
  130 continue
  140 return
      end
c*********************************************************************
      subroutine funcsb(n,xx,f)
      implicit real*8(a-h,o-z)
      parameter(nda=18888,npa=18888,mdl=22222,nr=016)
      common /triang/ tx,ty,adt(mdl),bdt(mdl),cdt(mdl),idt(nda),
     &                ldt(mdl,3),ndt(mdl,3),np,npex,nd
      common /print/ itp
      common /abicc/abic,dabic
*     common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,ph(npa,npa)
      common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,adp(npa),aup(nr,npa)
      common /xad/xd(nda),m
      common /sublin/ ilin
      common /prupd/ iupd
      common /files/ indata,outfile,an
      character*72 indata,outfile
c
      common /delink/ nzu(npa),iu(nr,npa)
      integer nzl(npa),il(nr,npa)
*     real*8 au(nr,npa),ad(npa),al(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)
      integer iv(npa)
      common /diagnl/ idet
c
      external funcs
      dimension xx(nda),x(nda)
      dimension  w(nda),ip(nda)
c     data xd/.01,299*0.0/
      data fmin/1.d60/
      data upd/0.0/
      m=npex
      ifg=0
      eps=0.1d-2
c
      w1=xx(1)**2
      w2=w1
c     w2=xx(2)**2
c10   if(xx(1).gt. 30.0 .or. xx(2).gt. 30.0) go to 100
c10   if(xx(1).lt.-30.0 .or. xx(2).lt.-30.0) go to 100
c10   w1=exp(xx(1))
c10   w2=exp(xx(2))
      w4=1.0d0
*     call prihes(m,ph)
      call prihes(m)
cremember! 3 cases following:
cex   call ltinv(ph,p4265./178./212.tde,m-1,nda)
c.....<<for the case where w1=0.0>>.........
c     write(6,*) 'ph'
c     call ltinv(ph,ptdet,m-1,nda)
c     call decomp(nda,m-1,ph,w,ip,ptdet)
c
      idet=1
*     call icdata(m-1,nr,ph,au,ad,y,wrk)
      call iccg (m, m-1, nr, nzu, nzl, iu, il,
     $           aup, al, adp, y, ptdet, wrk,
     $           uu, ul, ud, r, p, ru, ap, iv)
      idet=0
c
      write(6,3) ptdet
c     stop 222
    3 format(1h ,'ptdet =',d20.13)
c...........................................
    1 format(1h ,10f8.4)
*     call prihes(m,ph)
      call prihes(m)
      itp=0
      do 60 i=1,m
   60 x(i)=xd(i)
      do 20 i=1,1
   20 call davidn( x,m,0,funcs,xm )
      if(fmin.lt.abic) go to 50
      upd=upd+1
      iupd=upd+0.5
      fmin=abic
      do 10 i=1,m
   10 xd(i)=x(i)
      write(6,*) '-------- xd -------',upd,w1,abic
c     write(6,610) (xd(i),i=1,m)
  610 format(1h ,10f7.2)
*     open(11,file='work')
      open(11,file=outfile)
c     open(11)
      write(11,4) w1,abic,m
      write(11,6) (xd(i),i=1,m)
      close(11)
    6 format(5e15.5)
    4 format(2e15.5,i7)
   50 f=abic
c     if(ilin.eq.1) return
c     stop 119
      return
  100 continue
      f=1.d30
      ifg=1
      end
c*********************************************************************
*     subroutine prihes(n,h)
      subroutine prihes(n)
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        nd: 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
c    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=18888,npa=18888,mdl=22222,nr=016)
      common /triang/ tx,ty,adt(mdl),bdt(mdl),cdt(mdl),idt(nda),
     &                ldt(mdl,3),ndt(mdl,3),np,npex,nd
      common /data/ xx(nda),yy(nda),zz(nda,5),sd(nda),nn,mp(nda)
*     common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,ph(npa,npa)
      common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,adp(npa),aup(nr,npa)
      common /delink/ nzu(npa),iu(nr,npa)
*     dimension h(npa,npa)
      dimension ldi(3),xi(3),yi(3)
c
c     return
c
c     write(6,*) 'npex',npex,nd
      wx = w1
      wy = w2
c     write(6,1) w1,w2
    1 format(1h ,'#3: w1,w2 =',2d15.8)
c
*     do 10 ia=1,npex
*     do 10 ib=1,npex
*  10   h(ia,ib)=0.0
      do 10 i=1,n
        adp(i)=0.0
        do 10 j=1,nr
          aup(j,i)=0.0
   10 continue
c
      do 20 ipt=1,nd
        do 30 i=1,3
          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
        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) +
*    &      bdt(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) + bdt(ipt)/dn0**2 * (
     &      wx*(yi(inx)-yi(inx2))*(yi(jnx)-yi(jnx2)) +
     &      wy*(xi(inx)-xi(inx2))*(xi(jnx)-xi(jnx2))  )
          else
            do 60 k=1,nzu(il)
              if(iu(k,il).eq.jl) then
                aup(k,il)=aup(k,il) + bdt(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
c         if(il.eq.4243.and.jl.eq.4243) then
c         write(6,*) 'ph4167'
c         write(6,*) ipt,i,j,inx,inx2,jnx,jnx2
c         write(6,*) h(il,jl),bdt(ipt),dn0,wx,wy
c         write(6,*) xi(inx)-xi(inx2),xi(jnx)-xi(jnx2)
c         write(6,*) yi(inx)-yi(inx2),yi(jnx)-yi(jnx2)
c         end if
   40   continue
c
   20 continue
      return
      end
c*********************************************************************
      subroutine funcs(n,a,f,g,ifg)
c-----------------------------------------------------------------------
c      log likelihood function for intensity estimation.
c-----------------------------------------------------------------------
c   input
c        n : number of parameters
c        a : parameters
c   output
c        f : value of the function
c        g : gradient vecters
c        ifg : index for restrictions
c
      implicit real * 8 (a-h,o-z)
      parameter(nda=18888,npa=18888,mdl=22222,nr=016)
      common /triang/ tx,ty,adt(mdl),bdt(mdl),cdt(mdl),idt(nda),
     &                ldt(mdl,3),ndt(mdl,3),np,npex,nd
      dimension a(nda),g(nda),gp(nda)
      common /data/ xx(nda),yy(nda),zz(nda,5),sd(nda),nn,mp(nda)
*     common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,ph(npa,npa)
      common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,adp(npa),aup(nr,npa)
      common /dvfnhs/ rss,rss1,gh(nda),fpdum
      common /prupd/ iupd
      ifg=0
c
      do 10 ia=1,n
c       if(a(ia).gt.170.0) go to 70
        if(a(ia).gt.170.0) write(6,*) ia,a(ia)
   10   g(ia)=0.0
      f=0.0
c
      do 20 k=1,np
*     do 20 k=1,npex
c       do 25 ki=1,mp(k)  multiplication of data at a point.
          f=f+a(k)
          g(k)=g(k)+1.d0
c  25     continue
   20   continue
        do 80 id=1,nd
          ff=0.0
          do 90 k=1,3
            idk=ldt(id,k)
c           ff=ff+exp(a(idk))/3*adt(id)
c           g(idk)=g(idk) - exp(a(idk))/3*adt(id)
*           if(a(idk)+cdt(id).gt.170.0) go to 70
            ff=ff+exp(a(idk)+cdt(id))/3
            g(idk)=g(idk) - exp(a(idk)+cdt(id))/3
   90       continue
   80     f=f-ff
c
      f=-f
c     do 30 k=1,np
      do 30 k=1,npex
   30   g(k)=-g(k)
c
c     return
c
      fp=0.0
c 'n' is defined in subroutine funcsb.
      call penalt(n,a,fp,gp)
c     do 40 ia=1,n
c     do 40 ib=1,n
c       fp=fp+a(ia)*ph(ia,ib)*a(ib)/2
c  40   continue
      f=f+fp
      if(fp.lt.0.0) write(6,*) 'fp < 0 !',  fp
      if(fp.lt.0.0) write(6,1) (a(i),i=1,n)
c
      do 50 ia=1,n
c       pha=0.0
c       do 60 ib=1,n
c  60     pha=pha+ph(ia,ib)*a(ib)
c       g(ia)=g(ia)+pha
        g(ia)=g(ia)+gp(ia)
   50   continue
c     write(6,1) (a(i),i=1,n),f
c     write(6,1) (g(i),i=1,n)
    1 format(1h ,7f12.6)
      return
   70 continue
      f=1.d30
      ifg=1
      return
      end
c*********************************************************************
*     subroutine hes(n,a,h)
      subroutine hes(n,a,au,ad)
c----------------------------------------------------------------------
c     hesian matrix of the log-likelihood for intensity estimation.
c----------------------------------------------------------------------
c   input
c         nn: data length
c          n: number of parameters = npex
c       a(i): parameters
c   output
c     h(.,.): hessian matrix
c
      implicit real * 8 (a-h,o-z)
      parameter(nda=18888,npa=18888,mdl=22222,nr=016)
      common /triang/ tx,ty,adt(mdl),bdt(mdl),cdt(mdl),idt(nda),
     &                ldt(mdl,3),ndt(mdl,3),np,npex,nd
      common /data/ xx(nda),yy(nda),zz(nda,5),sd(nda),nn,mp(nda)
*     common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,ph(nr,npa)
      common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,adp(npa),aup(nr,npa)
*     dimension a(nda),h(npa,npa)
      dimension a(npa)
      common /delink/ nzu(npa),iu(nr,npa)
      dimension au(nr,npa),ad(npa)
c
c     return
c
*     do 10 i=1,n
*     do 10 j=1,n
*  10   h(i,j)=0.0
      do 10 i=1,n
        ad(i)=0.0
   10 continue
c
      do 80 id=1,nd
        do 90 k=1,3
          idk=ldt(id,k)
*         h(idk,idk)=h(idk,idk) - exp(a(idk))/3*adt(id)
c         ad(idk)=ad(idk)       - exp(a(idk))/3*adt(id)
          ad(idk)=ad(idk)       - exp(a(idk)+cdt(id))/3
   90   continue
   80  continue
c
*     do 30 ib=1,n
*     do 30 ia=1,n
*  30   h(ia,ib)=-h(ia,ib)+ph(ia,ib)
      do 20 i=1,n
        ad(i)=-ad(i)+adp(i)
        do 30 j=1,nzu(i)
          au(j,i)=aup(j,i)
   30   continue
   20 continue
      return
      end
c*********************************************************************
      subroutine penalt(n,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        nd: 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
c    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=18888,npa=18888,mdl=22222,nr=016)
      common /triang/ tx,ty,adt(mdl),bdt(mdl),cdt(mdl),idt(nda),
     &                ldt(mdl,3),ndt(mdl,3),np,npex,nd
      common /data/ xx(nda),yy(nda),zz(nda,5),sd(nda),nn,mp(nda)
*     common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,ph(npa,npa)
      common /swgrad/ w1,w2,w3,w4,wx,wxx,ptdet,adp(npa),aup(nr,npa)
      dimension a(nda),g(nda),gx(nda),gy(nda),gx2(nda),gy2(nda)
      dimension ldi(3),xi(3),yi(3)
c     data jj /0/
c
      wx = w1
      wy = w2
c     write(6,1) w1,w2
    1 format(1h ,'#3: w1,w2 =',2d15.8)
c
      fp=0.0
      do 10 i=1,npa
      g(i)=0.0
   10 continue
      do 20 ipt=1,nd
        fpx=0.0
        fpy=0.0
      do 25 i=1,npa
      gx(i)=0.0
      gy(i)=0.0
      gx2(i)=0.0
      gy2(i)=0.0
   25 continue
        do 30 i=1,3
          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
        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 + bdt(ipt) * (wx*fpx**2+wy*fpy**2) /2
        do 45 i=1,3
          il=ldt(ipt,i)
          g(il) = g(il) + bdt(ipt) * ( wx*fpx*gx(il)+wy*fpy*gy(il) )
   45     continue
   20 continue
      return
      end
c main:    test program for sample1 and sample2 function.
c funcs1:  sample1 function.
c funcs2:  sample2 function.
c simplx:  simplex minimization subroutine.
c minmax:  called by subroutine simplx.
c first:             "
c centor:            "
c newsim:            "
c update:            "
c reduce:            "
c epslon:            "
c
c define parameter.
c maxh:    max dimension of vector x.
c
c procedure to use this program.
c 1. save this mail to a file. (save file-name.f)
c 2. delete the comments. (vi file name.f)
c 3. change the define parameter and function.
c
c program starts. ------------------------------------------------------
c
c     parameter   (maxh=6, maxh5=maxh+5)
c
c     integer    n
c     real*8     xinit(maxh), dist, eps, f
c     external   funcs1, funcs2
c                                                     **** sample 1 ****
c     n = 2
c
c     xinit(1) = 10.0d0
c     xinit(2) = 10.0d0
c
c     dist = 1.0d0
c     eps = 10.0d-10
c
c     call simplx(xinit, n, funcs1, dist, eps, f)
c                                                     **** sample 2 ****
c     n = 6
c
c     xinit(1) = 0.0d0
c     xinit(2) = 0.0d0
c     xinit(3) = 0.0d0
c     xinit(4) = 0.0d0
c     xinit(5) = 0.0d0
c     xinit(6) = 0.0d0
c
c     dist = 1.0d0
c     eps = 10.0d-10
c
c     call simplx(xinit, n, funcs2, dist, eps, f)
c
c     stop
c     end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
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
c
      parameter   (maxh=6, maxh5=maxh+5)
c
      external   funct
c
      integer    n
      real*8     xinit(maxh), dist, eps, fr
c
      integer    iter, xh, xs, xl, x0, xr, xe, xc, j
      real*8     f(maxh5), x(maxh5,maxh)
      real*8     epslon
      real*8     eps1, alpha, beta, gamma
c
      data       alpha/-1.0d0/, beta/0.5d0/, gamma/2.0d0/
c
      x0 = n + 2
      xr = n + 3
      xe = n + 4
      xc = n + 5
c
      call first(n, f, x, funct, xinit, dist)
c -------------------------------------------------------------------- c
ctemp open(2, file='simplex.line')
      iter = 0
c
c                                                ****   main loop   ****
1     continue
c
c                                           ****   min max search   ****
      call minmax(n, f, xh, xs, xl)
c
c                                            ****   epsilon check   ****
      eps1 = epslon(n, f)
      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)
c
      if(iter .gt. 2 .and. eps1 .lt. eps) go to 2
c
c                                               ****   center set   ****
      call center(n, x, xh, x0)
c
c                                               ****   reflection   ****
      call newsim(n, f, x, funct, x0, xh, alpha, xr)
c
      if(f(xr) .le. f(xs)) then
c
        if(f(xr) .lt. f(xl)) then
c
c                                                ****   expansion   ****
          call newsim(n, f, x, funct, x0, xr, gamma, xe)
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
c
        else
          call update(n, f, x, xr, xh)
          end if
c
      else
c
        if(f(xr) .lt. f(xh)) call update(n, f, x, xr, xh)
c
c                                              ****   contraction   ****
        call newsim(n, f, x, funct, x0, xh, beta, xc)
c
        if(f(xc) .lt. f(xh)) then
          call update(n, f, x, xc, xh)
c
        else
c
c                                                   ****   reduce   ****
          call reduce(n, f, x, funct, xl)
c
          end if
c
        end if
c
      iter = iter + 1
c
      if(iter .gt. 1000 ) go to 2
c
c                                            ****   epsilon check   ****
      go to 1
c
2     continue
c
      do 10 j = 1, n
        xinit(j) = x(xl, j)
10      continue
      fr = f(xl)
c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine minmax(n, f, xh, xs, xl)
c
      parameter   (maxh=6, maxh5=maxh+5)
c
      integer    n, xh, xs, xl
      real*8     f(maxh5)
c
      integer    i
      real*8     fmin1, fmax1, fmax2
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
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine first(n, f, x, funct, xinit, dist)
c
      parameter   (maxh=6, maxh5=maxh+5)
c
      external   funct
c
      integer    n
      real*8     f(maxh5), x(maxh5,maxh), xinit(maxh), xx(maxh), dist
c
      integer    i, j
c
      do 10 i = 1, n + 1
c
        do 20 j = 1, n
          x(i,j) = xinit(j)
          xx(j) = x(i,j)
20        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
10      continue
c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine center(n, x, xh, x0)
c
      parameter   (maxh=6, maxh5=maxh+5)
c
      integer    n, xh, x0
      real*8     x(maxh5,maxh)
c
      integer    i, j
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
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      function   epslon(n, f)
c
      parameter   (maxh=6, maxh5=maxh+5)
c
      integer    n
      real*8     f(maxh5)
c
      integer    i
      real*8     epslon, emean
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
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine newsim(n, f, x, funct, x0, in, prm, out)
c
      parameter   (maxh=6, maxh5=maxh+5)
c
      external   funct
c
      integer    n, x0, in, out
      real*8     f(maxh5), x(maxh5,maxh), xx(maxh), prm
c
      integer    j
c
      do 10 j = 1, n
        x(out,j) = prm * x(in,j) + (1.0 - prm) * x(x0,j)
        xx(j) = x(out,j)
10      continue
      call funct(n, xx, f(out))
c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine update(n, f, x, in, out)
c
      parameter   (maxh=6, maxh5=maxh+5)
c
      integer    n, in, out
      real*8     f(maxh5), x(maxh5,maxh)
c
      integer    j
c
      do 10 j = 1 , n
        x(out,j) = x(in,j)
10      continue
c
      f(out) = f(in)
c
      return
      end
c -------------------------------------------------------------------- c
c -------------------------------------------------------------------- c
      subroutine reduce(n, f, x, funct, xl)
c
      parameter   (maxh=6, maxh5=maxh+5)
c
      external   funct
c
      integer    n, xl
      real*8     f(maxh5), x(maxh5,maxh), xx(maxh)
c
      integer    i, j
c
      do 10 i = 1, n+1
        if(i .ne. xl) then
          do 20 j = 1, n
            x(i,j) = (x(i,j) + x(xl,j)) / 2.0d0
            xx(j) = x(i,j)
20          continue
c
          call funct(n, xx, f(i))
          end if
10      continue
c
      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
*
      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
          if(ud(i).gt.0.0) then
            eps=eps-log(ud(i))
          else
            write(6,*) '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,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)
*
      do 540 i = 1, n
        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(22222)
      real*8 r(nd),p(nd),ru(nd),ap(nd)
*
      real*8 rur0,rur1,pap,res2
*
      kend = n + 10
      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.gt.999) then
          do 562 i=1,n
            x(i)=xxmin(i)
  562     continue
          write(6,*) 'cgres_1',kmin,res2,res2/b2
          return
        endif
*
c       if (res2/b2 .le. eps2) then
        if (res2.le.eps2**2 .or. res2/b2.le.eps2) then
          write(6,*) 'cgres_0',k,res2,res2/b2
*
*         ---- 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
