	subroutine manxct(calcfg)
c     fletchers	implementation of marquardts good idea
c
	include '/u1/gcf/e350/crystal/cb/hotsun'
c
c ****	  meaning of variables in /hotsun/ *************************************
c ....	  npar	    * number of	parameters
c ....	  ndata	    * number of	data points
c ....	  zeromn    * minimum chisq
c ....	  zerocr    * current chisq
c ....	  idata	    * used by accum to count function values
c ....	  numit	    * counts iterations
c
c ....	  param	    * parameter	values
c ....	  perr	    *  estimated errors	in parameter values
c ....	  pbest	    * saved copy of param for best solution
c ....	  first	    * first derivative vector
c ....	  firbes    * saved copy of first for best solution
c ....	  func	    * residuals	of each	data point
c ....	  fbest	    * saved copy of func  for best solution
c
c ....	  secpck    * packed second derivative matrix
c ....	  secbes    * saved copy of secpck for best solution
c ....	  lllpck    * packed matrix l where secpck= l *	l-transpose
c		    * packed matrix is stored a(1,1) / a(2,1) a(2,2) /
c					      a(3,1) a(3,2) a(3,3) etc.
c
c ....	  bnderr    * number of	consequitive boundary violations
c ....	  tcalc	    * time spent in calcfg
c ....	  tfit	    * time spent outside calcfg
c
c ....	  nbadgo    * stop if nbadgo tries produces less change	than dellim
c ....	  maxit	    * allowed number of	iterations
c ....	  dellim    * stop when	predicted change of chisq le dellim
c ....	  errcal    * =0 dont calculate	errors in fitted parameters * if nonzero
c		    * on input do and return errcal=1 if possible *
c					     errcal=-1 if impossible
c ....	  iprint    * 0	no print * =1 full print * =-1 final print only
c ....	  timmax    * stop if total time taken gt timmax
c
c ....	  endres    * reason for stopping fit
c		    * =1 iteration limit
c		    * =2 converged correctly
c		    * =3 time cut
c		    * =4 no progress made in nbadgo iterations
c		    * =5 boundary vioation
c		    * =6 matrix	singular in spite of adding lambda
c
c ****	  meaning of local variables *******************************************
c ....	  materr    * if nonzero , fit aborted due to singular matrix
c ....	  rho	    * fletchers	parameter rho
c ....	  sigma	    * fletchers	parameter sigma
c ....	  isaved    * =0 current best solt is stored in	current	(param,func)
c			      type arrays
c		    * =1 current best solt stored in save(secbes etc) arrays
c		    * =2 current best solt stored in current storage except
c		    *	for residuals which are	in fbest
c ....	  succ	    * =true, current fit reduced chisq , =false	it increased
c ....	  tbegin    * initial time
c ....	  q	    * marquardt	lambda parameter
c ....	  qlow	    * lower limit on q (i.e. estimate of minimum eigenvalue
c		    * of second	derivative matrix )
c ....	  qhigh	    * upper limit on q (i.e. estimate of maximum eigenvalue
c		    * of second	derivative matrix )
c ....	  qgcf	    * gcf's estimate of	best q
c ....	  qgood	    * value of q on last succesful iteration
c ....	  qscale    * force all	q changes to be	at least factor	of sqrt(qscale)
c ....	  xnorm	    * sum of xshift(i)**2 $ diag(i)
c ....	  isdtry    $ =0 steepest descent not tried * =1 try next go * =2 tried
c
c ....	  expchg    * expected chisq change of current fit
c ....	  pred1	    * predicted	chisq reduction	from linear term in expansion
c ....	  pred2	    * predicted	chisq reduction	from quadratic term in expansion
c
c ....	  dgsave    * save diagonal elements of	current	secpck
c ....	  diag	    * initial diagonal elements	of secpck
c ....	  sqdiag    * sqrt(diag)
c ....	  dginv	    * inverse of diag
c ....	  sqdinv    * inverse of sqdiag
c ....	  xshift    * minus predicted change in	parameter values
c ....	  chsave    * circular array to	save previous chisqs
c ....	  ichsav    * counts circularly	from 1..nbadgo positions in chsave
c
	dimension	diag(24),sqdiag(24),dgsave(24),xshift(24),chsave(6),
& elabel(2,6),dginv(24),sqdinv(24),iout1(2),method(2,2)
	character*10 elabel,ibl,method,iout1
	character*12 slabel
	logical violat,succ
	
	external calcfg
c
c ****	  data statements for elements of /hotsun/
	data timmax/100000./,dellim/.05/,maxit/20/,errcal/0/,
& iprint/-1/,nbadgo/6/,qscale/2.5/
c
c ****	  data statements for local variables
	data rho/.25/,sigma/.75/
	data elabel/10h iteration,6h limit,
&		  10h converged,3h ok,
&		  10h time cut ,1h ,
&		  10h no progre,2hss,
&		  10h boundary ,10h violation,
&		  10h matrix si,6hngular /
	data slabel/7hmanxcat/
	data ibl/1h /
	data method/10hmarquardt ,10halgorithm ,10h steepest ,7hdescent/
c
 1    format(29h *** start manxcat ****	maxit,i3,7h dellim,f8.3,
     1 7h nbadgo,i2,7h timmax,f8.2,4h rho,f8.3,6h sigma,f8.3,
     2 7h qscale,f8.2,11h **********)
 2    format(5h iter,i3,6h chisq,g12.2,10h marquardt,g11.3,
     1 6h q-low,g11.3,7h q-high,g11.3,11h q-gcf pred,g11.3,
     2 10h last ok q,g11.3/
     3 18h pred	chisq change,g12.3,7h actual,g12.3,
     4 9h fit-time,f10.3,5h calc,f10.3,4h ** ,2a10,4h ** ,2a10)
 3    format(11h parameters,1x,10g12.4/(1h ,11x,10g12.4))
 4    format(' resids'20f5.1/(1h ,20f5.1))
 5    format(9h	  errors,3x,10g12.4/(1h	,11x,10g12.4))
 6    format(54h ******************************	end manxcat fit	reason,
     1 i5,2x,2a10,49h ************************************************ )
 7    format(33h consequitive boundary violations,i5)
c
c     error if parameters out of allowed range
      if (npar.gt.24) call error(-npar,slabel)
      if (ndata.gt.500)	call error(ndata,slabel)
c
c     set initial time
      call second(tbegin)
      if(iprint.ne.0) print 1,maxit,dellim,nbadgo,timmax,rho,sigma,
     1 qscale
c
c     first call to user routine ** ignore violat
      call zrword(first,npar)
      nenter=npar*(npar+1)/2
      call zrword(secpck,nenter)
      zerocr=0.
      idata=0
      call calcfg(violat)
c
c     set intial parameters
      call second(tnow)
      tcalc=tnow-tbegin
      tfit=0.
      ichsav=1
      chsave(1)=zerocr
      expchg=0.
      delchi=0.
      isaved=0
      q=0.
      qlow=1.
      qhigh=1.
      qgcf=1.
      qgood=-1.
      isdtry=0
      numit=1
      bnderr=0
c
c     print intial trial if requested
      if(iprint.ne.1) go to 22
      print 2,numit,zerocr
      print 3,(param(i),i=1,npar)
      print 4,(func(i),i=1,ndata)
c
c     set up scaling from diagonal elements of initial second deriv matrix
 22   ipos=1
      do 10 i=1,npar
      diag(i)=abs(secpck(ipos))
      ipos=ipos+i+1
      sqdiag(i)=sqrt(diag(i))
      dginv(i)=1./diag(i)
 10   sqdinv(i)=1./sqrt(diag(i))
c
c     start of section used each iteration when	previous iteration successful
c     save diagonal elements of	secpck as marquardt addition will overwrite
 70   ipos=1
      do 31 i=1,npar
      dgsave(i)=secpck(ipos)
 31   ipos=ipos+i+1
c
c     start section used each iteration	whether	or not last iteration
c     succeeded
 35   materr=0
c
c     start section used if matrix failure in ludecp
c     set up modified chisq matrix with	marquardt addition
 75   if(isdtry.eq.1) go to 142
      if(isdtry.eq.2) isdtry=3
      ipos=1
      do 36 i=1,npar
      secpck(ipos)=dgsave(i)+q*diag(i)
 36   ipos=ipos+i+1
c
c     invert modified chisq matrix
c     determinant d1 d2	information ignored
      call ludecp(secpck,lllpck,npar,d1,d2,ier)
      if(ier.eq.0) go to 40
c
c     singular matrix *	increase lambda
      if(materr.gt.1) go to 82
      materr=materr+1
      q=2.*q
      if(q.lt.1.e-5) q=1.
      go to 75
c
c     steepest descent algorithm requested * note that lllpck is *not* set
c     under these circumstances
 142   do 143 i=1,npar
 143   xshift(i)=dginv(i)*first(i)
      z=0.
      ipos=1
      do 144 i=1,npar
      a=vdot1(i,secpck(ipos),1,xshift,1)
      if(i.eq.npar) go to 145
      a=a+vdot2(npar-i,xshift(i+1),1,secpck(ipos+i+i-1),i+1)
 145  z=z+a*xshift(i)
 144  ipos=ipos+i
c
      z=vdot1(npar,xshift,1,first,1)/z
      do 147 i=1,npar
 147  xshift(i)=z*xshift(i)
      isdtry=2
      qsd=q
      go to 141
c
c     matrix was inverted ok * solve chisq equation to find shift
 40   materr=0
      call luelmp(lllpck,first,npar,xshift)
c
c     find predicted chisq changes
 141  pred1=0.
      pred2=0.
      xnorm=0.
      ipos=1
      do 37 i=1,npar
      xnorm=xnorm+xshift(i)*xshift(i)*diag(i)
      pred1=pred1+2.*xshift(i)*first(i)
      a=dgsave(i)*xshift(i)
      i1=i-1
      if(i1.le.0) go to	38
      a=vdot1(i1,secpck(ipos),1,xshift,1)+a
 38   if(i.eq.npar) go to 39
      a=a+vdot2(npar-i,xshift(i+1),1,secpck(ipos+i+i-1),i+1)
 39   pred2=pred2-a*xshift(i)
 37   ipos=ipos+i
      expchg=pred1+pred2
c
c     if necessary * save current best parameters * note we store true
c     second derivative	matrix * not one overwritten by	marquardt
      if(isaved.ne.0) go to 41
      zeromn=zerocr
      call mvword(pbest,param,npar)
      call mvword(fbest,func,ndata)
      call mvword(secbes,secpck,nenter)
      ipos=1
      do 21 i=1,npar
      secbes(ipos)=dgsave(i)
 21   ipos=ipos+i+1
      call mvword(firbes,first,npar)
c
c     set up arrays for	calcfg
 41   call zrword(first,npar)
      call zrword(secpck,nenter)
      do 42 i=1,npar
 42   param(i)=param(i)-xshift(i)
      zerocr=0.
c
c     save current best	chisq in circular list
      ichsav=ichsav+1
      if(ichsav.gt.nbadgo) ichsav=1
      chsave(ichsav)=zeromn
c
c     call user	routine
      call second(temp)
      idata=0
      call calcfg(violat)
c
c     store boundary violation information
      if(violat) go to 13
      bnderr=0
      go to 12
 13   bnderr=bnderr+1
      xnorm=0.
      do 130 i=1,npar
      xshift(i)=pbest(i)-param(i)
 130  xnorm=xnorm+xshift(i)*xshift(i)*diag(i)
c
c     increment	iteration pointer and set change in chisq
 12   numit=numit+1
      delchi=zeromn-zerocr
c
c     set up timing information
      call second(tnow)
      taken=tnow-temp
      tcalc=tcalc+taken
      tfit=tnow-tbegin-tcalc
c
c     print iteration if requested
      if(iprint.ne.1) go to 23
      if(violat) go to 118
      iout1(1)=ibl
      iout1(2)=ibl
      go to 119
 118  iout1(1)=elabel(1,5)
      iout1(2)=elabel(2,5)
 119  ipos=1
      if(isdtry.eq.2) ipos=2
c
      print 2,numit,zerocr,q,qlow,qhigh,qgcf,qgood,
     1 expchg,delchi,tfit,tcalc,iout1,(method(i,ipos),i=1,2)
      print 3,(param(i),i=1,npar)
      print 4,(func(i),i=1,ndata)
c
 23   qgcf=2.*(expchg-delchi)/xnorm
      if(delchi.ge.0.) go to 43
      succ=.false.
      isaved=1
      go to 44
 43   succ=.true.
      isaved=0
      if(isdtry.eq.2) isdtry=0
      if(isdtry.ne.3) go to 44
      if(q.lt.qsd-1.e-10) isdtry=0
c
c     decide whether to	give up	the ghost or not ***********************
 44   if(tfit+tcalc.gt.timmax) go to 17
c
c     test on iteration	count
      if(numit.ge.maxit) go to 15
c
c     test on expected chisq change
      if(expchg.le.dellim) go to 16
c
c     test on boundary violation
      if(bnderr.gt.3) go to 19
c
c     stop if too little progress made in nbadgo steps * do this by
c     comparing	current	best chisq with	that found nbadgo iterations before
      if(numit.le.nbadgo)  go to 14
      itest=ichsav+1
      if(itest.gt.nbadgo) itest=1
      if(chsave(itest)-zeromn.lt.dellim) go to 18
c
c     another iteration	needed * decide	on strategy
 14   if(delchi.ge.rho*expchg) go to 60
      igood=0
c
c     unexpectedly small change	* increase marquardt parameter
c     reset limits qlow	and qhigh on range of q	values
c     note that	to save	time , qlow,qhigh are only reset if last q=0
c     this is important	for qlow calculation which is a) harder	and b)
c     assumes lllpck contains inverse info for true secpck not secpck plus
c     q	times diagonal matrix
      if(q.gt.1.e-10) go to 58
      qlow=0.
      qhigh=0.
      trh=0.
      trl=0.
      ipos=1
      jpos=1
      if(succ) go to 68
c
c     first case * previous go was a failure
c     find qlow	 * previous go failed
      call cmpinv(secpck)
      do 53 i=1,nenter
 53   secpck(i)=abs(secpck(i))
      do 54 i=1,npar
      trl=trl+secpck(ipos+i-1)*diag(i)
      z=vdot1(i,secpck(ipos),1,sqdiag,1)
      if(i.eq.npar) go to 56
      z=z+ vdot2(npar-i,sqdiag(i+1),1,secpck(ipos+i+i-1),i+1)
 56   z=z*sqdiag(i)
      if(z.gt.qlow) qlow=z
 54   ipos=ipos+i
c
c     find qhigh * previous go failed
      do 103 i=1,nenter
 103  secpck(i)=abs(secbes(i))
      do 104 i=1,npar
      trh=trh+secpck(jpos+i-1)*dginv(i)
      z=vdot1(i,secpck(jpos),1,sqdinv,1)
      if(i.eq.npar) go to 106
      z=z+vdot2(npar-i,sqdinv(i+1),1,secpck(jpos+i+i-1),i+1)
 106  z=z*sqdinv(i)
      if(z.gt.qhigh) qhigh=z
 104  jpos=jpos+i
      go to 69
c
c     second case * previous go	was a success
c     find qlow	 * previous go worked
 68   call cmpinv(secbes)
      do 63 i=1,nenter
 63   secbes(i)=abs(secbes(i))
      do 64 i=1,npar
      trl=trl+secbes(ipos+i-1)*diag(i)
      z=vdot1(i,secbes(ipos),1,sqdiag,1)
      if(i.eq.npar) go to 66
      z=z+ vdot2(npar-i,sqdiag(i+1),1,secbes(ipos+i+i-1),i+1)
 66   z=z*sqdiag(i)
      if(z.gt.qlow) qlow=z
 64   ipos=ipos+i
c
c     find qhigh * previous go worked
      do 113 i=1,nenter
 113  secbes(i)=abs(secpck(i))
      do 114 i=1,npar
      trh=trh+secbes(jpos+i-1)*dginv(i)
      z=vdot1(i,secbes(jpos),1,sqdinv,1)
      if(i.eq.npar) go to 116
      z=z+vdot2(npar-i,sqdinv(i+1),1,secbes(jpos+i+i-1),i+1)
 116  z=z*sqdinv(i)
      if(z.gt.qhigh) qhigh=z
 114  jpos=jpos+i
c
c     at this stage trl	is trace of a-1	, trh trace of a
c     qlow is l-infinity norm of a-1 and qhigh l-infinity norm of a
 69   if(trh.lt.qhigh) qhigh=trh
      if(trl.lt.qlow) qlow=trl
      qlow=1./qlow
      q=sqrt(qlow*qhigh)
      if(qhigh.gt.qscale*q) go to 120
      go to 117
c
c     qlow,qhigh already set * decide on new marqardt parameter
 58   if(qhigh.le.qscale*q) go to 121
      if(qgood.le.qscale*q) go to 127
      q=sqrt(q*qgood)
      go to 120
 127  q=sqrt(q*qhigh)
      go to 120
c
c     q	is getting comparable to qhigh * things	look really bad
c     find z where lambda-new( i.e. q new) = lambda-old	* z
 121  if(isdtry.ne.0) go to 128
 117  isdtry=1
      go to 120
 128  z=pred1-delchi
      if(z.gt.0.) go to	126
      z=2.
      go to 129
 126  z=z*2./pred1
      if(z.lt.2.) z=2.
      if(z.gt.10.) z=10.
 129  q=q*z
c
c     q	reset *	reset current parms if last go failed
 120  if(succ) go to 70
      isaved=2
      call mvword(first,firbes,npar)
      call mvword(secpck,secbes,nenter)
      call mvword(param,pbest,npar)
      go to 35
c
c     reasonable success * decrease marquardt parameter	if going very well
 60   qgood=q
      igood=igood+1
      if(q.lt.1.e-10) go to 70
      if(delchi.le.sigma*expchg) go to 70
      if(q.le.qscale*qlow) go to 122
      q=sqrt(q*qlow)
      go to 70
 122  if(igood.lt.3) go	to 70
      q=0.
      go to 70
c     end determination	of strategy for	next iteration *****************
c
c     end section * error as singular matrix
 82   endres=6
      ipos=1
      do 83 i=1,npar
      secpck(ipos)=dgsave(i)
 83   ipos=ipos+i+1
      call second(tnow)
      tfit=tnow-tbegin-tcalc
      if(isaved.eq.2) call mvword(func,fbest,ndata)
      go to 81
c
c     set termination parameter
 15   endres=1
      go to 80
 16   endres=2
      go to 80
 17   endres=3
      go to 80
 18   endres=4
      go to 80
 19   endres=5
c
c     end of fit * first restore current parameters if necessary
 80   if(isaved.eq.0) go to 81
      call mvword(func,fbest,ndata)
      zerocr=zeromn
      if(isaved.eq.2) go to 81
      call mvword(first,firbes,npar)
      call mvword(secpck,secbes,nenter)
      call mvword(param,pbest,npar)
c
c     calculate	errors if requested
c     when errors found	* secpck holds in packed form full error matrix
c     while perr has sqrts of diagonal elements
 81   if(errcal.eq.0) go to 84
      call ludecp(secpck,lllpck,npar,d1,d2,ier)
      if(ier.eq.0) go to 85
c
c     errors cannot be calculated
      errcal=-1
      do 86 i=1,npar
 86   perr(i)=-1.
      go to 84
c
c     errors can be calculated
 85   errcal=1
      call cmpinv(secpck)
      ipos=1
      do 87 i=1,npar
      perr(i)=sqrt(abs(secpck(ipos)))
 87   ipos=ipos+i+1
c
c     final print out if requested
 84   if(iprint.eq.0) return
      if(materr.ne.0) return
c
      print 2,numit,zerocr,q,qlow,qhigh,qgcf,qgood,
     1 expchg,delchi,tfit,tcalc
      print 3,(param(i),i=1,npar)
      if(errcal.eq.1) print 5,(perr(i),i=1,npar)
      write(6,4) (func(i),i=1,ndata)
      if(bnderr.ne.0) print 7,bnderr
245   print 6,endres,(elabel(i,endres),i=1,2)
      return
c
      end
