c	the subroutine below computes the values of the shrink
c	coefficents for use in the non-negative garrotte. 
c	It assumes that all variables have been centered at zero by
c	subtracting off their sample means.  M is the number of 
c	variables, SS is the X'X matric, and XY is the X'Y vector.
c	Both of these are destroyed in the subroutine.  YY is the
c	sum of squares of the y-values, and BOLS is the vector
c	containing the ols regression coefficients.  YY is altered
c	in the subroutine, but not BOLS.

c	SUMB is the desired sum of the non-negative garotte coefficents.
c	On exit from the subroutine, SHRCF is the vector carrying
c	the non-negative garotte coeficients that sum to SUMB.  The
c	coeficient of the kth variable in the nn-garotte regression
c	is SHRCF(K)*BOLS(K).

c	note the parameter setting of MMAX=60 in the BPLUS subroutine.
c	if M>60, then MMAX must be increased to be >=M.  
c	The algorithm is a modification of the Lawson-Hansen non-
c	negative least squares algorithm.  The sum constraint is
c	gotten by a barrier method where the size of the barrier
c	is set by the constant CON. The sum is checked and if its
c	not in the bounds (1+-TOL)*SUMB then CON is increased and 
c	the process repeated.

	SUBROUTINE GAROTTE(SS,XY,YY,BOLS,M,SUMB,SHRCF)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	PARAMETER(TOL=.001)
	DIMENSION SS(M,M),XY(M),SHRCF(M),BOLS(M)

	CON=1000
	
	DO 10 I=1,M
	XY(I)=BOLS(I)*XY(I)/SUMB
10	CONTINUE
	
	DO 20 I=1,M
	DO 30 J=1,M
	SS(I,J)=BOLS(I)*BOLS(J)*SS(I,J)
30	CONTINUE
20 	CONTINUE
	
	YY=YY/(SUMB*SUMB)
	
100	DO 40 I=1,M
	XY(I)=XY(I)+(CON/SUMB)
	DO 50 J=1,M
	SS(I,J)=SS(I,J)+(CON/SUMB)
50	CONTINUE
40	CONTINUE
	YY=YY+(CON/SUMB)
	
	CALL BPLUS(SS,XY,YY,M,SHRCF)
	
	SSUM=0
	DO 60 K=1,M
	SSUM=SSUM+SHRCF(K)
	SHRCF(K)=SUMB*SHRCF(K)
60	CONTINUE
	IF (DABS(SSUM-1.0).GT.TOL) THEN
	CON=10*CON
	GOTO 100
	END IF
	END
	
		
	SUBROUTINE BPLUS(SS,XY,YY,M,BT)
	IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	PARAMETER(MMAX=200)
	DIMENSION SS(M,M),XY(M),BT(M),U(MMAX+1,MMAX+1),PT(MMAX),
     *	W(MMAX),BS(MMAX),IZ(MMAX)
	DO 10 K=1,M
	DO 20 J=1,M
	U(J,K)=SS(J,K)
20	CONTINUE
10	CONTINUE
	DO 30 J=1,M
	U(J,M+1)=XY(J)
	U(M+1,J)=XY(J)
30	CONTINUE
	U(M+1,M+1)=YY
	
	DO 40 K=1,M
	IZ(K)=1
	BT(K)=0
40	CONTINUE
	ITER=0
	ITER1=0
c	start on Lawson and Hansen algorithm

1	DO 50 I=1,M
	PT(I)=0
	DO 60 J=1,M
	PT(I)=PT(I)+SS(I,J)*BT(J)
60	CONTINUE
50	CONTINUE
	
	DO 70 K=1,M
	W(K)=XY(K)-PT(K)
70	CONTINUE
	NZ=0
	NX=0
	WMAX=0
	DO 80 K=1,M
	IF (IZ(K).EQ.1) THEN
		NZ=NZ+1
		IF (W(K).GT.WMAX) THEN
			NX=K
			WMAX=W(K)
		END IF
	END IF
80	CONTINUE
	IF (NZ.EQ.0.OR.NX.EQ.0) RETURN
	IZ(NX)=0
	CALL SWEEP(U,M+1,NX)
2	BMIN=1E20
	DO 90 K=1,M
	IF (IZ(K).EQ.1) THEN
	BS(K)=.0
	ELSE
	BS(K)=U(K,M+1)
	BMIN=DMIN1(BS(K),BMIN)
	END IF
90	CONTINUE
	
	IF (BMIN.GT.0.) THEN
	DO 100 K=1,M
	BT(K)=BS(K)
100	CONTINUE
	ITER1=ITER1+1
	GOTO 1
	END IF
	
	AL=1E20
	IQ=0
	DO 110 K=1,M
	IF (IZ(K).EQ.0.AND.BS(K).LE.0.) THEN
		IF (BT(K)-BS(K).EQ.0) THEN
		WRITE(*,*) "DENOMINATOR PROBLEM"
		STOP
		END IF
		RT=BT(K)/(BT(K)-BS(K))
		IF (RT.LT.AL) THEN
			IQ=K
			AL=RT
		END IF
	END IF
110	CONTINUE
	IF (IQ.EQ.0) THEN
	WRITE(*,*) "IQ FAILURE"
	STOP
	END IF
	IF (AL.EQ.0) THEN
	WRITE(*,*) "AL ZERO"
	STOP
	ENDIF
	IZ(IQ)=1
	CALL SWEEP(U,M+1,IQ)
	
	DO 120 K=1,M
	BT(K)=BT(K)+AL*(BS(K)-BT(K))
120	CONTINUE
	ITER=ITER+1
	IF (ITER.GT.100000) THEN
	Write(*,*) "TO0 MANY ITERATIONS"
	STOP
	END IF
	GOTO 2
	
	END
	
	
c	this is the SWEEP operation
	SUBROUTINE SWEEP(V,M,LL)
	PARAMETER(MMAX=200)
	DOUBLE PRECISION V(MMAX+1,MMAX+1),TD,CT
	TD=V(LL,LL)
	IF (TD.LT.1.0E-10) THEN
		WRITE(*,*) "SMALL"
		END IF
	DO 10 J=1,M
		V(LL,J)=V(LL,J)/TD
10	CONTINUE 
	
	DO 20 I=1,M
		IF (I.EQ.LL) THEN
			GOTO 20
		END IF
		CT=V(I,LL)
		DO 30 J= 1,M
			V(I,J)=V(I,J)-CT*V(LL,J)
30		CONTINUE 
		V(I,LL)=-CT/TD
20	CONTINUE 
	V(LL,LL)=1/TD
	RETURN
	END
	  
	  
	
	
