c----------------------------------------------------------------------------
c    File: edro9903.for
c
c    Sample main program to test subroutine box9903
c
c    This is the ``Hard spheres problem''
c
c    It consists of finding npun points in the unitary sphere of 
c    \R^{ndim} maximizing the minimimum distance between any pair
c    of them.
c
c    In other words, finding npun points in the unitary sphere of 
c    \R^{ndim} minimizing the maximum scalar product between any
c    pair of them
c
c
c
c                            Minimize z
c
c                            subject to
c 
c      z \geq <x_i, x_j> for all different pair of indices i, j
c
c                            ||x_i||^2 = 1    for all i=1, npun
c
c     The goal is to find an objective value less than 0.5 (This means
c     that the npun points stored belong to the sphere and every distance
c     between two of them is greater than 1.
c
c     Introducing slack variables w_{ij}, the inequality constraints can
c     be written as
c
c              <x_i, x_j> + w_{ij} - z = 0     for all i, j =1, npun
c
c     and                w{ij} \geq 0      for all i, j =1, npun
c
c
c     So, we have n = ndim * npun + npun (npun-1)/2 + 1 variables
c
c     and         nonlc = npun * (npun-1)/2 + npun    nonlinear constraints
c
c     The variables of the problem are organized in the array x in
c     the following way: the first ndim*npun variables are the coordinates
c     of the npun points;  the following npun (npun-1)/2 variables are
c     the slack variables. The last coordinate of x corresponds to z

	parameter(nnpun=40)
	parameter(nndim=5)
	parameter(nn=nndim*nnpun+nnpun*(nnpun-1)/2+1)
	parameter(nnonlc=nnpun*(nnpun-1)/2+nnpun)
	parameter(nnbd=nn+1)
	implicit double precision (a-h,o-z)
	double precision l(nn), lint(nn), lambda(1), mu(nnonlc)
	dimension x(nn), u(nn), g(nn), hess(1),  uint(nn)
        dimension s(nn), gq(nn), gc(nn), bd(nnbd), xn(nn), gn(nn)
	dimension d(nn), ingc(nn), infree(nn), inda(1)
	dimension a(1), b(1)
	dimension res(1), hres(nnonlc)
	dimension inhess(1)
	dimension dismat(nnpun, nnpun)
	dimension rainic(nnonlc), ra(nnonlc), hresan(nnonlc)
	character arq1*20
	character arq2*20
	character arq3*20
	character arq4*20
	character arq5*20
	character arq6*20

	external phi, mapro, guess
	dimension v1(3), v2(3), pv(3), difer(100), indif(100)
	Real*4 Dtime, tarray(2)
	Real*8 Tini, Tfim
	common/idoble/ indo(nnpun,nnpun), jump
	common/bolas/npun,ndim,n2,nlc
	common/gues/diago

	write(*,*)' Dear user: you must define six output files.'
	write(*,*)

	write(*,*)' Name of first output file:'
	read(*,'(a20)') arq1
	open (unit=10,file=arq1)

	write(*,*)' Name of second output file:'
	read(*,'(a20)') arq2
	open (unit=20,file=arq2)


	write(*,*)' Name of third output file:'
	read(*,'(a20)') arq3
	open (unit=30,file=arq1)

	write(*,*)' Name of fourth output file:'
	read(*,'(a20)') arq4
	open (unit=10,file=arq1)

	write(*,*)' Name of fifth output file:'
	read(*,'(a20)') arq5
	open (unit=10,file=arq5)

	write(*,*)' Name of sixth output file:'
	read(*,'(a20)') arq6
	open (unit=10,file=arq6)





	disbet = 0.d0
	ineq = 0
	write(10,*)
	write(20,*)
	write(30,*)
	write(40,*)
	write(50,*)
	write(60,*)
	write(10,*)' Output of Box9903 with Kissing problem:'
	write(10,*)
	write(*,*)' ilag = (0 if penalty, 1 if Augmented Lagrangian)'
	read(*,*) ilag
	write(*,*)' iguess = (Example: 0) '
	read(*,*) iguess
	if(iguess.eq.1)then
	write(*,*)' Element added on the diagonal of L: (Example: 10)'
	read(*,*) diago
	write(10,*)' Guess option used. Element added to diagonal:',
     *               diago
	endif
	if(ilag.eq.0)then
	write(10,*)' Pure external penalty method'
	else
	write(10,*)' Augmented Lagrangian method'
	endif

	write(*,*)' If you want to use a simplified Hessian type 1:'
	read(*,*) jump

	if(jump.eq.1) then
	write(10,*) ' Gauss-Newton simplified Hessian used'
	else
	write(10,*) ' Full Hessian used'
	endif

	write(10,*)
	write(*,*)' Dimension of the space (Example: 3):'
	read(*,*) ndim
	write(*,*)' Number of points (Example: 12):'
	read(*,*) npun
	k=0
	do 22 i=1,npun-1
	do 22 j=i+1,npun
	if(i.eq.j)then
	indo(i,j)=0
	else
	k=k+1
	indo(i,j)=k
	indo(j,i)=k
	endif
22	continue
	
	write(10,*)' Dimension of the space:', ndim,
     *     '  Number of points:', npun
	write(20,*)' Dimension of the space:', ndim,
     *     '  Number of points:', npun
	write(30,*)' Dimension of the space:', ndim,
     *     '  Number of points:', npun
	write(40,*)' Dimension of the space:', ndim,
     *     '  Number of points:', npun
	write(50,*) ndim,npun
	write(*,*)' Maximum number of initial points tried: 
     *             (Example:1)'
	read(*,*) maxini
	n2 = npun*(npun-1)/2
	nonlc = n2 + npun
	n=ndim*npun+n2+1
	nlc=nonlc
	write(10,*)' Number of pairs of different points:',n2
	write(10,*)' Number of nonlinear constraints:',nonlc
	write(10,*)' Number of variables (including slack):',n
	write(*,*)' Number of pairs of different points:',n2
	write(*,*)' Number of nonlinear constraints:',nonlc
	write(10,*)' Number of variables of the nonlinear program:',n
	write(*,*)' Number of variables of the nonlinear program:',n
	isem=237
	write(10,*)' Seed for the initial point:',isem
	do 1 i=1,npun*ndim
	l(i)=-1000.d0
1	u(i)=1000.d0
	do 4 i=npun*ndim+1, npun*ndim+n2
	l(i)=0.d0
4	u(i)=1000.d0
	l(n)=-1000.d0
	u(n)=1000.d0

	accuracy=0.1d0
	acumenor = 0.5d0
	write(10,*)' Accuracy in the resolution of the subproblem:',
     *       accuracy
	write(10,*)' Acumenor:', acumenor
	relarm = 1.d-4
	absarm = 1.d-6
	write(10,*)' relarm =', relarm,'  absarm = ', absarm

	write(*,*)' eps inicial, eps final, eps-div:'
	write(*,*)' Example: 0.1, 1.d-5, 10'
	read(*,*) epsinic, epsfin, epsdiv

	write(10,*)' Epsilon for projected gradient (initial/final):'
	write(10,*)' Epsinic:', epsinic,'  epsfin:', epsfin
	write(10,*) ' If Initial epsilon is not Final epsilon, '
	write(10,*)' Epsilon is divided by ', epsdiv,' after each outer'
	write(10,*) ' iteration.'
	write(10,*)' When it reaches final epsilon it is not divided '
	write(10,*) ' anymore.'
	write(*,*)' If new  feasibility < ratfea * Old feasibility'
	write(*,*)' we do not increase penalty parameter.'
	write(*,*)' ratfea = (Example: 0.01)'
	read(*,*) ratfea
	write(*,*)' Delta - min = (example: 0.001)'
	read(*,*) deltamin
	write(10,*)' Delta - min = ', deltamin
	
	write(10,*)' If new  feasibility < ratfea * Old feasibility'
	write(10,*)' we do not increase penalty parameter.'
	write(10,*)' ratfea =', ratfea

	epsd=1.d-8
	write(10,*)' epsd=', epsd
	itmax = 100
	write(10,*)' Maximum number of iterations at each call to Box:'
     *     ,    itmax
	write(*,*)' iprint = (Example: 1)'
	read(*,*) iprint

	nafmax=10*itmax
	par=4.d0
	delta0=10.d0
	dont=0.d0
	kdont=100
	eta=0.95d0
	factor=10.d0
	mitqu= n
	write(10,*)' Maximum number of iterations for quacan calls:',
     *                mitqu
	maxmv=10*mitqu
	write(*,*)' iprbox, iprquacan: (Example: -1  -1)'
	read(*,*) ipr, iprq

	bont = 0.d0
	kbont = 1000
	facra = 10.d0
	rai = 10.d0
	kauts = 1000
	indivi = 0
	epsnon = 1.d-8
	maxout = 20
	m = 0

	write(10,*) ' bont:', bont,'   kbont:', kbont
	write(10,*)' Initial penalty parameter:', rai
	write(10,*)' Augmenting factor facra:', facra
	write(10,*)' Parameter kauts:', kauts
	write(10,*)' Parameter indivi:', indivi
c
c   Set the initial point
c
	do 46 initial =1,maxini
	do 2 i=1,n
2	x(i) = 2.d0 * rnd(isem) - 1.d0

	do 3 i = 1, nonlc
3	rainic(i) = rai

12      continue

c	call gettim(ihour,imin,isec,ihund)

	Tini = DTime(tarray)

c       o tempo de CPU esta armazenado em tarray(1)


	call box9903 (n, m, nonlc, roinic, ra, rainic, indivi, 
     *  facro, facra, ilag, iguess, epslin, epsnon, ratfea,
     *  maxout, konout, konint, nefint, itqint, mvpint, ierla, iprint, 
     *  x, l , u, fx, hres, hresan, res, g,
     *  accuracy, acumenor, relarm, absarm, epsinic, epsdiv, epsfin,
     *  epsd, nafmax, itmax, par, delta0, deltamin, istop, 
     *  bont, kbont, kauts, dont, kdont, eta, factor, 
     *  mitqu, maxmv, phi, mapro, guess, hess, inhess, ipr, 
     *  iprq, lint, uint, s, gq, gc, bd, xn, gn, d, ingc, infree,
     *  mu, lambda, ispar, nelem, inda, a, b)


	Tfim = Dtime(tarray)

c       o tempo de CPU esta armazenado em tarray(1)


c	call gettim(ihour1,imin1,isec1,ihund1)
c	texec=3.6e+5*(ihour1-ihour)+6.e+3*(imin1-imin)+
c    *        1.e+2 *(isec1-isec)+(ihund1-ihund)
c	texec = texec/100.d0
	texec = tarray(1)
	write(50,51)konout,konint,nefint,itqint,mvpint,texec
51	format(1x,i3,i7,i9,i10,i10,f8.1)	
	write(10,*)' CPU time (seconds) = ',texec
	write(*,*)' CPU time CPU (seconds) = ',texec

c    Normalization
	do 42 i=1,npun
	z = 0.d0
	do 43 j=1,ndim
43	z = z + x((i-1)*ndim+j)**2
	z = dsqrt(z)
	do 45 j = 1, ndim
45	x((i-1)*ndim+j) = x((i-1)*ndim+j)/z
42	continue


c    Compute the minimum distance between two points 
	dist= 1.d30
	dismaxi=0.d0
	do 14 j=1,npun-1
	do 15 i=j+1,npun
	dismat(i,j)=0.d0
	do 16 k=1,ndim
16	dismat(i,j)=dismat(i,j)+ (x((i-1)*ndim+k)-x((j-1)*ndim+k))**2
	dismat(i,j)=dsqrt(dismat(i,j))
	dismaxi = dmax1(dismaxi, dismat(i,j))
	dist = dmin1(dist, dismat(i,j))
15	dismat(j,i)=dismat(i,j)
14	dismat(j,j)=0.d0
	
	if(dist.gt.disbet) then
	disbet = dist
	write(60,*)
	write(60,*)' Initial approx. number ', initial
	write(60,*)' gave best distance up to now.'
	write(60,*)' The points are:'
	do 44 i=1,npun
	write(60,*) (x((i-1)*ndim+j), j=1,ndim)
44	continue
	write(60,*)
	write(60,*)' The minimum distance is :', disbet
	endif

	write(10,*)' Initial  point number ', initial
	write(*,*)' Initial  point number ', initial
	write(10,*)' Minimum distance between two points:',dist
	write(*,*)' Minimum distance between two points:',dist
c	if(ndim.eq.3)then
	write(20,41)
41	format(/////)
	write(40,*)
	write(30,*) ' New Polytope:'
	write(40,*) ' New Polytope:'
	write(30,*)
	write(20,*)
	write(20,*)' Initial point: ', initial,'  seed:', isem
	write(40,*)' Initial point: ', initial,'  seed:', isem
	write(20,*)' Minimum distance between two points:', dist
	write(20,*)' Maximum distance between two points:', dismaxi
	write(40,*)' Minimum distance between two points:', dist
	write(40,*)' Maximum distance between two points:', dismaxi
	write(20,*)' Matriz de distancias:'
	
	do 17 i=1,npun
	write(20,18)(dismat(i,j),j=1,npun)
17	continue
18	format(1x, 13f5.2)
	write(20,*)
	write(20,*)'    The points are:'
	do 19 i=1,npun
	write(20,20) i, (x((i-1)*ndim+k), k=1,ndim)
20	format(i3, 3x, 5f7.2)
19	continue
	ndifer=0
	do 36 i = 1, npun-1
	do 36 j = i+1, npun
	if(ndifer.eq.0) then
	difer(1)=dismat(1,2)
	ndifer = 1
	indif(1) = 1 
	else
	do 37 k = 1, ndifer
	if(dabs(dismat(i,j)-difer(k)).le.0.01d0)then
	indif(k) = indif(k)+1
	goto36
	endif
37	continue
	ndifer = ndifer+1
	difer(ndifer) = dismat(i,j)
	indif(ndifer) = 1
	endif
36	continue

	do 38 k=1,ndifer
	write(20,39) difer(k), indif(k)
	write(40,39) difer(k), indif(k)
	
39	format(' Distance ', f6.2,'  appears ', i3,'  times')
38	continue

	if(ndim.eq.3)then
	
	
	ncara = 0
	do 21 i=1,npun-2
	do 21 j=i+1, npun-1
	do 21 k=j+1, npun
	do 24 ll=1,3
	v1(ll) = x((j-1)*ndim+ll) - x((i-1)*ndim+ll)
24	v2(ll) = x((k-1)*ndim+ll) - x((i-1)*ndim+ll)
	pv(1) = v1(2)*v2(3) - v1(3)*v2(2)
	pv(2) = - v1(1)*v2(3) + v1(3) * v2(1)
	pv(3) = v1(1)*v2(2) - v1(2)*v2(1)
	pesa=0.d0
	do 25 np=1,npun
	if(np.eq.i.or.np.eq.j.or.np.eq.k)goto25
	pes=0.d0
	do 26 ll=1,3
26	pes=pes + pv(ll) * (x((np-1)*ndim+ll)-x((i-1)*ndim+ll))
	if(pesa*pes.lt.0.d0)goto21
	pesa=pes
25	continue
	ncara = ncara+1
	write(20,27)i, j, k, ncara
27	format(' The vertices ', i3, i3, i3, ' determin face', i3)
	write(20,40) i,j,dismat(i,j)
40	format(' distance between vertices',i4,' and ',i4,' is',f6.2)
	write(20,40) i,k,dismat(i,k)
	write(20,40) j,k,dismat(j,k)

	
	write(30,*)' Polygon[{'
	write(30, 33) x((i-1)*ndim+1),x((i-1)*ndim+2),x((i-1)*ndim+3)
	write(30, 33) x((j-1)*ndim+1),x((j-1)*ndim+2),x((j-1)*ndim+3)
33	format(' {', f7.2,',',f7.2,',',f7.2,'},')

	do 28 np=1,npun
	if(np.eq.i.or.np.eq.j.or.np.eq.k)goto28
	pes=0.d0
	z1=0.d0
	z2=0.d0
	do 29 ll=1,3
	pes=pes + pv(ll) * (x((np-1)*ndim+ll)-x((i-1)*ndim+ll))
	z1=z1+ pv(ll)**2
29	z2=z2+ (x((np-1)*ndim+ll)- x((i-1)*ndim+ll))**2
	cosen= pes/dsqrt(z1*z2)
c	write(20,31)np, pes, cosen
31	format(' Ponto ', i3,' pes:', f12.4,'  cosen:', f12.4)
	if(dabs(cosen).le.0.1d0)then
	write(20,30)np
30	format(' The point ',i3, ' probably belongs to the face above.')
	write(30, 33) x((np-1)*ndim+1),x((np-1)*ndim+2),x((np-1)*ndim+3)
	endif
28	continue	

	write(30, 34) x((k-1)*ndim+1),x((k-1)*ndim+2),x((k-1)*ndim+3)
34	format(' {', f7.2,',',f7.2,',',f7.2,'}}],')


21	continue
	write(40,*)' Number of triangular faces detected for this'
	write(40,*)' polytope:', ncara
	endif

46	continue
	close (10)
	close (20)
	close (30)
	close (40)
	close (50)
	close (60)

	stop
	end


	subroutine phi (n, nonlc, x, mu, ra, f, fx, g, hres,  
     *   hess, inhess, modo )
	implicit double precision (a-h,o-z)
	parameter(nnpun=40)
	parameter(nndim=5)
	parameter(nn=nndim*nnpun+nnpun*(nnpun-1)/2+1)
	parameter(nnonlc=nnpun*(nnpun-1)/2+nnpun)
	parameter(nnbd=nndim*nnpun+nnpun*(nnpun-1)/2+2)
	double precision mu(*)
	dimension x(n), g(n), hres(*), hess(*), inhess(*), ra(*)
	common/bolas/npun,ndim,n2,nlc
	common/ra/rara(nnonlc)
	common/info/xx(nn),hhres(nnonlc),amu(nnonlc)
	if(modo.eq.1)then
	fx=x(n)
	n1=ndim*npun
	nrnl=0
	do 1 i=1,npun-1
	do 1 j=i+1, npun
	pes=0.d0
	do 2 k=1,ndim
2	pes=pes+ x((i-1)*ndim+k)*x((j-1)*ndim+k)
	nrnl=nrnl+1
1	hres(nrnl) = pes + x(n1+nrnl) - x(n)
c
	do 3 i = 1, npun
	nrnl=nrnl+1
	xnor=0.d0
	do 4 j=1,ndim
4	xnor=xnor + x(ndim*(i-1)+j)**2
	hres(nrnl) = xnor - 1.d0
3	continue
c
	hnor=0.d0
	do 5 i=1,nonlc
5	hnor=hnor+hres(i)**2
	tlan=0.d0
	do 6 i=1,nonlc
6	tlan=tlan+mu(i)*hres(i)
c
	f = fx + tlan
	do 14 i = 1, nonlc
14	f = f + ra(i)*hres(i)**2 / 2.d0


	return
	else
c
c    Attention, we are using the fact that callings to fun with
c    modo=2 always is preceeded by a call with the same point
c    with modo=1. 
c
	do 15 i=1,nonlc
	hhres(i)=hres(i)
15	rara(i)=ra(i)
	g(n)=1.d0
	do 7 i=1,n-1
7	g(i)=0.d0
	nrnl=0
	do 8 i=1,npun-1
	do 8 j=i+1,npun
	nrnl=nrnl+1
	amult= mu(nrnl) + ra(nrnl)*hres(nrnl)
	amu(nrnl) = amult
	g(n1+nrnl) = g(n1+nrnl) + amult
	g(n) = g(n) - amult
	do 9 k=1,ndim
	g((i-1)*ndim+k) = g((i-1)*ndim+k) + amult *x((j-1)*ndim+k)
	g((j-1)*ndim+k) = g((j-1)*ndim+k) + amult *x((i-1)*ndim+k)
9	continue
8	continue
	do 12 i=1,npun
	nrnl=n2+i
	amult=mu(nrnl) + ra(nrnl)*hres(nrnl)
	amu(nrnl)=amult
	do 13 j=1,ndim
13	g((i-1)*ndim+j) = g((i-1)*ndim+j)+
     *     amult * 2.d0 * x((i-1)*ndim+j)
12	continue
c
c
	do 11 i=1,n
11	xx(i)=x(i)

	endif
	return
	end





	subroutine mapro (n, nind, ind, u, v, hess, inhess )
	implicit double precision (a-h,o-z)
	parameter(nnpun=40)
	parameter(nndim=5)
	parameter(nn=nndim*nnpun+nnpun*(nnpun-1)/2+1)
	parameter(nnonlc=nnpun*(nnpun-1)/2+nnpun)
	parameter(nnbd=nndim*nnpun+nnpun*(nnpun-1)/2+2)

	dimension ind(nind), u(n), v(n), hess(*), inhess(*)
	common/bolas/npun,ndim,n2,nonlc
	common/ra/ra(nnonlc)
	common/info/x(nn),hres(nnonlc),amu(nnonlc)
	common/idoble/ indo(nnpun,nnpun), jump

	dimension aux(nn)
c
c     Product h'(x) u

c     First p(p-1)/2 rows of h'(x) u
c
c     Product aux = P u + something
c
	ndipun = ndim*npun	

	k=0
	do 1 i = 1, npun-1
	do 1 j = i+1, npun
	k=k+1
	aux(k) = u(ndipun+k) - u(n) 
	do 2 l = 1, ndim
	aux(k)=aux(k)+x((j-1)*ndim+l)*u((i-1)*ndim+l)
     *         + x((i-1)*ndim+l)*u((j-1)*ndim+l)
2	continue
1	continue

c    Last p rows of h'(x) u
	
	do 3 i = 1, npun
	k = k+1
	aux(k)=0.d0
	do 4 l = 1, ndim
4	aux(k) = aux(k) + 2.d0 * x((i-1)*ndim+l) * u((i-1)*ndim+l)
3	continue
c
	do 5 i=1, nonlc
5	aux(i) = aux(i) * ra(i)
c
c    Product h'(x)^T aux
c
c    First block of ndim*npun rows
c
c    Product v = P^T aux
c
	do 7 i=1,ndipun
7	v(i) = 0.d0
	k = 0
	do 8 i=1,npun-1
	do 8 j= i+1, npun
	k = k+1
	do 9 l = 1, ndim
	v((i-1)*ndim+l) = v((i-1)*ndim+l) + x((j-1)*ndim+l)*aux(k)
	v((j-1)*ndim+l) = v((j-1)*ndim+l) + x((i-1)*ndim+l)*aux(k)
9	continue
8	continue

	do 10 j=1,npun
	do 11 l = 1, ndim
11	v((j-1)*ndim+l) = v((j-1)*ndim+l) + 2.d0 * x((j-1)*ndim+l)
     *       * aux(n2 + j)   
10	continue




c    Second block of p(p-1)/2 rows
c
	do 6 i=1,n2
6	v(ndipun+i) = aux(i)
c

c    Last component

	v(n) = 0.d0
	do 12 j=1,n2
12	v(n) = v(n) -  aux(j)

	if(jump.eq.1) return


c    Product (Matriz S) x u
	do 13 i = 1, npun
	do 14 l = 1, ndim
14	aux((i-1)*ndim+l) =0.d0
	do 15 j = 1, npun
	if(i.eq.j)then
	do 16 l = 1, ndim
16	aux((i-1)*ndim+l) = aux((i-1)*ndim+l) 
     *    + 2.d0 * amu(n2+i) * u((j-1)*ndim+l)
	else
	k = indo(i,j)
	do 17 l = 1, ndim
17	aux((i-1)*ndim+l) = aux((i-1)*ndim+l) 
     *    + amu(k) * u((j-1)*ndim+l)
	endif
15	continue
13	continue

	do 18 i = 1, ndipun
18	v(i) = v(i) + aux(i)

	return
	end


	subroutine guess(n, x, grad, ro, ra, lambda, mu, lint, uint, s)
	parameter(nnpun=40)
	parameter(nndim=5)
	parameter(nn=nndim*nnpun+nnpun*(nnpun-1)/2+1)
	parameter(nnonlc=nnpun*(nnpun-1)/2+nnpun)
	parameter(nndipun = nnpun * nndim)
	implicit double precision (a-h,o-z)
	common/bolas/npun,ndim,n2,nlc
	common/gues/diago
	double precision x(n), grad(n), lint(n), s(n)
	dimension uint(n)
	dimension ind(1), hess(1), inhess(1)
	dimension ra(*)
	dimension smat (nndim,nndim), aux(nndim), auxbig(nn)
	dimension al1u(nndim)
	double precision lambda(*),mu(*)

	ndipun = ndim*npun

c    Compute the matrix Smat = \Sum x_i x_i^T  +  diago  Identity

	do 1 i=1,ndim
	do 2 j=1,i
	smat(i,j)=0.d0
	do 3 k=1,npun
3	smat(i,j) = smat(i,j) + x((k-1)*ndim+i) * x((k-1)*ndim+j)
2	continue
	smat(i,i) = smat(i,i) + diago
1	continue

c   Compute the Cholesky decomposition of matrix Smat
	call chole(nndim, ndim, smat)

	do 4 i=1,npun
	call sicho(nndim, ndim, smat, s((i-1)*ndim+1),
     *      grad((i-1)*ndim+1), aux)
c
c    Apply small Sherman-Morrison

	call sicho(nndim, ndim, smat, al1u, x((i-1)*ndim+1), aux)
	pescanum=0.d0
	pescaden=0.d0
	do 5 j=1,ndim
	pescanum = pescanum + al1u(j)*grad((i-1)*ndim+j)
5	pescaden = pescaden + al1u(j)*x((i-1)*ndim+j)
	
	denom = 1.d0 + 2.d0 * pescaden
	factor = 2.d0 * pescanum / denom

	do 6 j=1,ndim
	s((i-1)*ndim+j) = 
     *  s((i-1)*ndim+j) -  al1u(j) * factor
6	continue
4	continue


c    Apply big Sherman-Morrison

c    Compute ``big A^{-1} u'' in auxbig
	
	do 40 i=1,npun

	call sicho(nndim, ndim, smat, auxbig((i-1)*ndim+1),
     *      x((i-1)*ndim+1), aux)
c
c    Apply small Sherman-Morrison

	call sicho(nndim, ndim, smat, al1u, x((i-1)*ndim+1), aux)
	pescanum=0.d0
	pescaden=0.d0
	do 50 j=1,ndim
	pescanum = pescanum + al1u(j)*x((i-1)*ndim+j)
50	pescaden = pescaden + al1u(j)*x((i-1)*ndim+j)
	
	denom = 1.d0 + 2.d0 * pescaden
	factor = 2.d0 * pescanum / denom

	do 60 j=1,ndim
	auxbig((i-1)*ndim+j) = 
     *  auxbig((i-1)*ndim+j) - al1u(j) * factor
60	continue
40	continue

	pescanum = 0.d0
	do 7 i=1,ndipun
	pescanum = pescanum + auxbig(i) * grad(i)
7	pescaden = pescaden + auxbig(i) * x(i)
	denom = 1.d0 + pescaden
	factor = pescanum/denom
	do 8 i=1,ndipun
8	s(i) = - s(i) + factor * auxbig(i)

	pesca=0.d0
	do 9 i=1,ndipun
9	pesca = pesca +  s(i) * grad(i)
	write(10,*)' Scalar product <s (guess), gradient> = ', pesca


	break=1.d30
	do 10 i=1,ndipun
	if(s(i).eq.0.d0)goto10
	step = uint(i)/s(i)
	if(step.gt.0.d0)then
	break = dmin1(break,step)
	goto10
	endif
	step = lint(i)/s(i)
	if(step.gt.0.d0)then
	break = dmin1(break,step)
	endif
10	continue


	do 12 i=ndipun+1, n
12	s(i) = 0.d0


	call mapro(n, nind, ind, s, auxbig, hess, inhess )
	denom = 0.d0
	do 13 i=1,n
13	denom = denom + s(i) * auxbig(i)

	step = -pesca/denom

	if(step.gt.break) step = break

	do 14 i=1,n
14	s(i) = step * s(i)

	return
	end




	subroutine chole (nlin, n, a)
	implicit double precision (a-h,o-z)
	dimension a(nlin,n)
	a(1,1) = dsqrt(a(1,1))
	if(n.eq.1)return
	do 1 i=2,n
	do 2 j=1,i-1
	z = 0.d0
	if(j.gt.1)then
	do 3 k=1,j-1
3	z = z + a(i,k) * a(j,k)
	endif
	a(i,j) = (a(i,j) - z)/a(j,j)
2	continue
	z = 0.d0
	do 4 j=1,i-1
4	z = z + a(i,j)**2
	a(i,i) = dsqrt( a(i,i) - z )
1	continue
	return
	end

	subroutine sicho (nlin, n, a, x, b, aux)
	implicit double precision (a-h,o-z)
	dimension a(nlin,n), x(n), b(n), aux(n)
	aux(1) = b(1)/a(1,1)
	if(n.gt.1)then
	do 1 i=2,n
	z = 0.d0
	do 2 j=1,i-1
2	z = z + a(i,j)*aux(j)
	aux(i) = (b(i) - z) / a(i,i)
1	continue
	endif
	x(n) = aux(n)/a(n,n)
	if(n.eq.1)return
	do 3 i=n-1,1,-1
	z = 0.d0
	do 4 j=i+1,n
4	z = z + a(j,i)*x(j)
	x(i) = (aux(i) - z)/a(i,i)
3	continue
	return
	end