[u/mrichter/AliRoot.git] / MINICERN / mathlib / gen / divon / delslv.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1996/04/01 15:03:25  mclareni
* Mathlib gen
*
*
#include "gen/pilot.h"
      SUBROUTINE DELSLV(N,FMAJOR,FMINOR,LMAX,FRACT,X,XLOW,XUP,VOL,NCUT,N
     1CDIM,ICUT,DELPLS,DELNEG,REGTOL,FTOL,FORIG,DFORIG,FNEW,FNLIN,FNLROW
     2,DIAGJ,SPDIAG,SOL,Z,NFCNT)
      INTEGER N, NCUT, NCDIM, NFCNT
      INTEGER ICUT(NCDIM)
      DOUBLE PRECISION FMAJOR, FMINOR, FRACT, REGTOL, FTOL
      DOUBLE PRECISION X(N), XUP(N), XLOW(N), DELPLS(N), DELNEG(N)
      DOUBLE PRECISION FORIG(NCDIM), FNLIN(NCDIM)
      DOUBLE PRECISION DFORIG(NCDIM), FNEW(NCDIM), FNLROW(NCDIM)
      DOUBLE PRECISION DIAGJ(NCDIM), SPDIAG(NCDIM), SOL(NCDIM), Z(N)
      LOGICAL LMAX
      INTEGER I, IAB, IBACK, II, ISAVE, ITRY,
     1        NCUTM1, NNEAR
      DOUBLE PRECISION BIG, DEL, DELMAX, DELMIN, DFNEW, DFUN,
     1        FDIF, FGAM, FNORM, FNRMNW, FOMX, FRAT, FZ, GAMMA,
     2        GAMNEW, PROD, RATGAM, REGINV, SINGTL, TSTVAL,
     3        VAL, VOL, XMULT, YDI, YDMIN
CMM   INTEGER IABS
      DATA SINGTL/ 1.0D-4/
      DATA BIG/ 1.0D+10/
      NFCNT=0
      FDIF=FMAJOR-FMINOR
      PROD=1.0D+0
      DO 10 I=1,N
      PROD=PROD*(DELPLS(I)+DELNEG(I))
 10   CONTINUE
      GAMMA=PROD/VOL
      FGAM=GAMMA*FMAJOR+(1.0D+0-GAMMA)*FMINOR
      DO 20 I=1,N
      Z(I)=X(I)
 20   CONTINUE
      DO 50 I=1,NCUT
      II=ICUT(I)
      IAB=ABS(II)
      IF(II.LT.0) GOTO 30
      Z(IAB)=X(IAB)+DELPLS(IAB)
      GOTO 40
 30   Z(IAB)=X(IAB)-DELNEG(IAB)
 40   FNEW(I)=DFUN(N,Z)
      Z(IAB)=X(IAB)
      NFCNT=NFCNT+1
 50   CONTINUE
#if defined(CERNLIB_IBM)||defined(CERNLIB_SINGLE)
 60   YDMIN=1.0D+40
#endif
#if (!defined(CERNLIB_IBM))&&(defined(CERNLIB_DOUBLE))
   60 YDMIN= 1.0D+30
#endif
      NNEAR=0
      DO 70 I=1,NCUT
      YDI=ABS(FMAJOR-FNEW(I))
      IF(YDI.GT.YDMIN) GOTO 70
      YDMIN=YDI
      NNEAR=I
 70   CONTINUE
      IF(NNEAR.EQ.0) RETURN
      IF((LMAX.AND.FNEW(NNEAR).LT.FGAM).OR.
     1(.NOT.LMAX.AND.FNEW(NNEAR).GT.FGAM)) GOTO 130
      ISAVE=ICUT(NNEAR)
      IF(NNEAR.EQ.NCUT.OR.NCUT.EQ.1) GOTO 90
      NCUTM1=NCUT-1
      DO 80 I=NNEAR,NCUTM1
      ICUT(I)=ICUT(I+1)
      FNEW(I)=FNEW(I+1)
 80   CONTINUE
 90   NCUT=NCUT-1
      IF(ISAVE.LT.0) GOTO 100
      DELPLS(ISAVE)=XUP(ISAVE)-X(ISAVE)
      GOTO 110
 100  IAB=ABS(ISAVE)
      DELNEG(IAB)=X(IAB)-XLOW(IAB)
 110  PROD=1.0D+0
      DO 120 I=1,N
      PROD=PROD*(DELPLS(I)+DELNEG(I))
 120  CONTINUE
      GAMMA=PROD/VOL
      FGAM=GAMMA*FMAJOR+(1.0D+0-GAMMA)*FMINOR
      IF(NCUT.EQ.0) RETURN
      GOTO 60
 130  DO 160 I=1,NCUT
      II=ICUT(I)
      IAB=ABS(II)
      IF(II.LT.0) GOTO 140
      DEL=DELPLS(IAB)
      GOTO 150
 140  DEL=DELNEG(IAB)
 150  FORIG(I)=FNEW(I)
      DFORIG(I)=(FNEW(I)-FMAJOR)/DEL
 160  CONTINUE
      CALL FEQN(NCUT,FORIG,FGAM,FNLIN)
      CALL RLEN(NCUT,FNLIN,FNORM)
 170  DO 180 I=1,NCUT
      FNLIN(I)=-FNLIN(I)
 180  CONTINUE
      IF(NCUT.EQ.1) GOTO 200
      DIAGJ(1)=DFORIG(1)
      SPDIAG(1)=-DFORIG(2)
      NCUTM1=NCUT-1
      DO 190 I=1,NCUTM1
      DIAGJ(I)=DFORIG(I)
      SPDIAG(I)=-DFORIG(I+1)
 190  CONTINUE
 200  DO 210 I=1,NCUT
      II=ICUT(I)
      II=ABS(II)
      FNLROW(I)=-GAMMA*FDIF/(DELPLS(II)+DELNEG(II))
 210  CONTINUE
      FNLROW(1)=DFORIG(1)+FNLROW(1)
      IF(NCUT.EQ.1) GOTO 230
      DO 220 I=1,NCUTM1
      XMULT=0.0D+0
      IF(ABS(FNLROW(I)).LT.BIG*ABS(DIAGJ(I))) XMULT=FNLROW(I)/DIAGJ(I
     1)
      FNLROW(I+1)=FNLROW(I+1)-XMULT*SPDIAG(I)
      FNLIN(NCUT)=FNLIN(NCUT)-XMULT*FNLIN(I)
 220  CONTINUE
 230  SOL(NCUT)=FNLIN(NCUT)
      IF(ABS(FNLROW(NCUT)).LT.BIG*ABS(FNLIN(NCUT))) SOL(NCUT)=FNLIN(N
     1CUT)/FNLROW(NCUT)
      IF(NCUT.EQ.1) GOTO 250
      DO 240 I=2,NCUT
      IBACK=NCUT-I+1
      VAL=FNLIN(IBACK)-SOL(IBACK+1)*SPDIAG(IBACK)
      SOL(IBACK)=VAL
      IF(ABS(DIAGJ(IBACK)).LT.BIG*ABS(VAL)) SOL(IBACK)=VAL/DIAGJ(IBAC
     1K)
 240  CONTINUE
 250  ITRY=0
      DO 280 I=1,NCUT
      II=ICUT(I)
      IAB=ABS(II)
      IF(II.LT.0) GOTO 260
      DELMAX=FRACT*(XUP(IAB)-X(IAB)-DELPLS(IAB))
      DELMIN=-DELPLS(IAB)
      GOTO 270
 260  DELMAX=FRACT*(X(IAB)-XLOW(IAB)-DELNEG(IAB))
      DELMIN=-DELNEG(IAB)
 270  IF(SOL(I).GT.DELMAX) SOL(I)=0.75D+0*DELMAX
      IF(SOL(I).LT.DELMIN) SOL(I)=0.75D+0*DELMIN
 280  CONTINUE
 290  DO 310 I=1,NCUT
      II=ICUT(I)
      IAB=ABS(II)
      IF(II.LT.0) GOTO 300
      DELPLS(IAB)=DELPLS(IAB)+SOL(I)
      GOTO 310
 300  DELNEG(IAB)=DELNEG(IAB)+SOL(I)
 310  CONTINUE
      DO 340 I=1,NCUT
      II=ICUT(I)
      IAB=ABS(II)
      IF(II.LT.0) GOTO 320
      Z(IAB)=X(IAB)+DELPLS(IAB)
      GOTO 330
 320  Z(IAB)=X(IAB)-DELNEG(IAB)
 330  FZ=DFUN(N,Z)
      NFCNT=NFCNT+1
      FNEW(I)=FZ
      Z(IAB)=X(IAB)
 340  CONTINUE
      PROD=1.0D+0
      DO 350 I=1,N
      PROD=PROD*(DELPLS(I)+DELNEG(I))
 350  CONTINUE
      GAMNEW=PROD/VOL
      FGAM=GAMNEW*FMAJOR+(1.0D+0-GAMNEW)*FMINOR
      CALL FEQN(NCUT,FNEW,FGAM,FNLIN)
      CALL RLEN(NCUT,FNLIN,FNRMNW)
      IF(FNRMNW.GT.FNORM) GOTO 380
      FOMX=0.0D+0
      DO 360 I=1,NCUT
      IF(ABS(FNEW(I)).GT.FOMX) FOMX=ABS(FNEW(I))
      DFNEW=FNEW(I)-FORIG(I)
      TSTVAL=1.0D+0
      IF(ABS(SOL(I)).LT.BIG*ABS(DFNEW)) TSTVAL=DFNEW/SOL(I)
      IF(ABS(TSTVAL).LT.SINGTL*ABS(DFORIG(I))) TSTVAL=SINGTL*DFORIG(I
     1)
      DFORIG(I)=TSTVAL
      FORIG(I)=FNEW(I)
 360  CONTINUE
      FNORM=FNRMNW
      FOMX=MAX(FOMX,ABS(FGAM))
      FRAT=FNORM/(1.0D+0+FOMX)
      REGINV=1.0D+0/REGTOL
      RATGAM=GAMNEW/GAMMA
      GAMMA=GAMNEW
C---  Activate to do debugging
C     WRITE(6,420) FNORM,FOMX,FRAT
C     WRITE(6,430) FGAM,RATGAM
      FGAM=GAMMA*FMAJOR+(1.0D+0-GAMMA)*FMINOR
      IF(RATGAM.GT.REGTOL.AND.RATGAM.LT.REGINV) GOTO 370
      IF(FRAT.LT.FTOL) GOTO 370
      GOTO 170
 370  RETURN
 380  ITRY=ITRY+1
      IF(ITRY.GT.2) RETURN
      DO 410 I=1,NCUT
      II=ICUT(I)
      IAB=ABS(II)
      IF(II.GE.0) THEN
         DELPLS(IAB)=DELPLS(IAB)-SOL(I)
      ELSE
         DELNEG(IAB)=DELNEG(IAB)-SOL(I)
      ENDIF
      SOL(I)=SOL(I)*0.25D+0
 410  CONTINUE
      GOTO 290
 420  FORMAT(' FNORM, FOMX, FRAT', 3(1PD15.5))
 430  FORMAT(' FGAM, RATGAM', 2(1PD15.5))
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1996/04/01 15:03:25 mclareni
	6	* Mathlib gen
	7	*
	8	*
	9	#include "gen/pilot.h"
	10	SUBROUTINE DELSLV(N,FMAJOR,FMINOR,LMAX,FRACT,X,XLOW,XUP,VOL,NCUT,N
	11	1CDIM,ICUT,DELPLS,DELNEG,REGTOL,FTOL,FORIG,DFORIG,FNEW,FNLIN,FNLROW
	12	2,DIAGJ,SPDIAG,SOL,Z,NFCNT)
	13	INTEGER N, NCUT, NCDIM, NFCNT
	14	INTEGER ICUT(NCDIM)
	15	DOUBLE PRECISION FMAJOR, FMINOR, FRACT, REGTOL, FTOL
	16	DOUBLE PRECISION X(N), XUP(N), XLOW(N), DELPLS(N), DELNEG(N)
	17	DOUBLE PRECISION FORIG(NCDIM), FNLIN(NCDIM)
	18	DOUBLE PRECISION DFORIG(NCDIM), FNEW(NCDIM), FNLROW(NCDIM)
	19	DOUBLE PRECISION DIAGJ(NCDIM), SPDIAG(NCDIM), SOL(NCDIM), Z(N)
	20	LOGICAL LMAX
	21	INTEGER I, IAB, IBACK, II, ISAVE, ITRY,
	22	1 NCUTM1, NNEAR
	23	DOUBLE PRECISION BIG, DEL, DELMAX, DELMIN, DFNEW, DFUN,
	24	1 FDIF, FGAM, FNORM, FNRMNW, FOMX, FRAT, FZ, GAMMA,
	25	2 GAMNEW, PROD, RATGAM, REGINV, SINGTL, TSTVAL,
	26	3 VAL, VOL, XMULT, YDI, YDMIN
	27	CMM INTEGER IABS
	28	DATA SINGTL/ 1.0D-4/
	29	DATA BIG/ 1.0D+10/
	30	NFCNT=0
	31	FDIF=FMAJOR-FMINOR
	32	PROD=1.0D+0
	33	DO 10 I=1,N
	34	PROD=PROD*(DELPLS(I)+DELNEG(I))
	35	10 CONTINUE
	36	GAMMA=PROD/VOL
	37	FGAM=GAMMAFMAJOR+(1.0D+0-GAMMA)FMINOR
	38	DO 20 I=1,N
	39	Z(I)=X(I)
	40	20 CONTINUE
	41	DO 50 I=1,NCUT
	42	II=ICUT(I)
	43	IAB=ABS(II)
	44	IF(II.LT.0) GOTO 30
	45	Z(IAB)=X(IAB)+DELPLS(IAB)
	46	GOTO 40
	47	30 Z(IAB)=X(IAB)-DELNEG(IAB)
	48	40 FNEW(I)=DFUN(N,Z)
	49	Z(IAB)=X(IAB)
	50	NFCNT=NFCNT+1
	51	50 CONTINUE
	52	#if defined(CERNLIB_IBM)\|\|defined(CERNLIB_SINGLE)
	53	60 YDMIN=1.0D+40
	54	#endif
	55	#if (!defined(CERNLIB_IBM))&&(defined(CERNLIB_DOUBLE))
	56	60 YDMIN= 1.0D+30
	57	#endif
	58	NNEAR=0
	59	DO 70 I=1,NCUT
	60	YDI=ABS(FMAJOR-FNEW(I))
	61	IF(YDI.GT.YDMIN) GOTO 70
	62	YDMIN=YDI
	63	NNEAR=I
	64	70 CONTINUE
65	IF(NNEAR.EQ.0) RETURN
66	IF((LMAX.AND.FNEW(NNEAR).LT.FGAM).OR.
67	1(.NOT.LMAX.AND.FNEW(NNEAR).GT.FGAM)) GOTO 130
68	ISAVE=ICUT(NNEAR)
69	IF(NNEAR.EQ.NCUT.OR.NCUT.EQ.1) GOTO 90
70	NCUTM1=NCUT-1
71	DO 80 I=NNEAR,NCUTM1
72	ICUT(I)=ICUT(I+1)
73	FNEW(I)=FNEW(I+1)
74	80 CONTINUE
75	90 NCUT=NCUT-1
76	IF(ISAVE.LT.0) GOTO 100
77	DELPLS(ISAVE)=XUP(ISAVE)-X(ISAVE)
78	GOTO 110
79	100 IAB=ABS(ISAVE)
80	DELNEG(IAB)=X(IAB)-XLOW(IAB)
81	110 PROD=1.0D+0
82	DO 120 I=1,N
83	PROD=PROD*(DELPLS(I)+DELNEG(I))
84	120 CONTINUE
85	GAMMA=PROD/VOL
86	FGAM=GAMMAFMAJOR+(1.0D+0-GAMMA)FMINOR
87	IF(NCUT.EQ.0) RETURN
88	GOTO 60
89	130 DO 160 I=1,NCUT
90	II=ICUT(I)
91	IAB=ABS(II)
92	IF(II.LT.0) GOTO 140
93	DEL=DELPLS(IAB)
94	GOTO 150
95	140 DEL=DELNEG(IAB)
96	150 FORIG(I)=FNEW(I)
97	DFORIG(I)=(FNEW(I)-FMAJOR)/DEL
98	160 CONTINUE
99	CALL FEQN(NCUT,FORIG,FGAM,FNLIN)
100	CALL RLEN(NCUT,FNLIN,FNORM)
101	170 DO 180 I=1,NCUT
102	FNLIN(I)=-FNLIN(I)
103	180 CONTINUE
104	IF(NCUT.EQ.1) GOTO 200
105	DIAGJ(1)=DFORIG(1)
106	SPDIAG(1)=-DFORIG(2)
107	NCUTM1=NCUT-1
108	DO 190 I=1,NCUTM1
109	DIAGJ(I)=DFORIG(I)
110	SPDIAG(I)=-DFORIG(I+1)
111	190 CONTINUE
112	200 DO 210 I=1,NCUT
113	II=ICUT(I)
114	II=ABS(II)
115	FNLROW(I)=-GAMMA*FDIF/(DELPLS(II)+DELNEG(II))
116	210 CONTINUE
117	FNLROW(1)=DFORIG(1)+FNLROW(1)
118	IF(NCUT.EQ.1) GOTO 230
119	DO 220 I=1,NCUTM1
120	XMULT=0.0D+0
121	IF(ABS(FNLROW(I)).LT.BIG*ABS(DIAGJ(I))) XMULT=FNLROW(I)/DIAGJ(I
122	1)
123	FNLROW(I+1)=FNLROW(I+1)-XMULT*SPDIAG(I)
124	FNLIN(NCUT)=FNLIN(NCUT)-XMULT*FNLIN(I)
125	220 CONTINUE
126	230 SOL(NCUT)=FNLIN(NCUT)
127	IF(ABS(FNLROW(NCUT)).LT.BIG*ABS(FNLIN(NCUT))) SOL(NCUT)=FNLIN(N
128	1CUT)/FNLROW(NCUT)
129	IF(NCUT.EQ.1) GOTO 250
130	DO 240 I=2,NCUT
131	IBACK=NCUT-I+1
132	VAL=FNLIN(IBACK)-SOL(IBACK+1)*SPDIAG(IBACK)
133	SOL(IBACK)=VAL
134	IF(ABS(DIAGJ(IBACK)).LT.BIG*ABS(VAL)) SOL(IBACK)=VAL/DIAGJ(IBAC
135	1K)
136	240 CONTINUE
137	250 ITRY=0
138	DO 280 I=1,NCUT
139	II=ICUT(I)
140	IAB=ABS(II)
141	IF(II.LT.0) GOTO 260
142	DELMAX=FRACT*(XUP(IAB)-X(IAB)-DELPLS(IAB))
143	DELMIN=-DELPLS(IAB)
144	GOTO 270
145	260 DELMAX=FRACT*(X(IAB)-XLOW(IAB)-DELNEG(IAB))
146	DELMIN=-DELNEG(IAB)
147	270 IF(SOL(I).GT.DELMAX) SOL(I)=0.75D+0*DELMAX
148	IF(SOL(I).LT.DELMIN) SOL(I)=0.75D+0*DELMIN
149	280 CONTINUE
150	290 DO 310 I=1,NCUT
151	II=ICUT(I)
152	IAB=ABS(II)
153	IF(II.LT.0) GOTO 300
154	DELPLS(IAB)=DELPLS(IAB)+SOL(I)
155	GOTO 310
156	300 DELNEG(IAB)=DELNEG(IAB)+SOL(I)
157	310 CONTINUE
158	DO 340 I=1,NCUT
159	II=ICUT(I)
160	IAB=ABS(II)
161	IF(II.LT.0) GOTO 320
162	Z(IAB)=X(IAB)+DELPLS(IAB)
163	GOTO 330
164	320 Z(IAB)=X(IAB)-DELNEG(IAB)
165	330 FZ=DFUN(N,Z)
166	NFCNT=NFCNT+1
167	FNEW(I)=FZ
168	Z(IAB)=X(IAB)
169	340 CONTINUE
170	PROD=1.0D+0
171	DO 350 I=1,N
172	PROD=PROD*(DELPLS(I)+DELNEG(I))
173	350 CONTINUE
174	GAMNEW=PROD/VOL
175	FGAM=GAMNEWFMAJOR+(1.0D+0-GAMNEW)FMINOR
176	CALL FEQN(NCUT,FNEW,FGAM,FNLIN)
177	CALL RLEN(NCUT,FNLIN,FNRMNW)
178	IF(FNRMNW.GT.FNORM) GOTO 380
179	FOMX=0.0D+0
180	DO 360 I=1,NCUT
181	IF(ABS(FNEW(I)).GT.FOMX) FOMX=ABS(FNEW(I))
182	DFNEW=FNEW(I)-FORIG(I)
183	TSTVAL=1.0D+0
184	IF(ABS(SOL(I)).LT.BIG*ABS(DFNEW)) TSTVAL=DFNEW/SOL(I)
185	IF(ABS(TSTVAL).LT.SINGTLABS(DFORIG(I))) TSTVAL=SINGTLDFORIG(I
186	1)
187	DFORIG(I)=TSTVAL
188	FORIG(I)=FNEW(I)
189	360 CONTINUE
190	FNORM=FNRMNW
191	FOMX=MAX(FOMX,ABS(FGAM))
192	FRAT=FNORM/(1.0D+0+FOMX)
193	REGINV=1.0D+0/REGTOL
194	RATGAM=GAMNEW/GAMMA
195	GAMMA=GAMNEW
196	C--- Activate to do debugging
197	C WRITE(6,420) FNORM,FOMX,FRAT
198	C WRITE(6,430) FGAM,RATGAM
199	FGAM=GAMMAFMAJOR+(1.0D+0-GAMMA)FMINOR
200	IF(RATGAM.GT.REGTOL.AND.RATGAM.LT.REGINV) GOTO 370
201	IF(FRAT.LT.FTOL) GOTO 370
202	GOTO 170
203	370 RETURN
204	380 ITRY=ITRY+1
205	IF(ITRY.GT.2) RETURN
206	DO 410 I=1,NCUT
207	II=ICUT(I)
208	IAB=ABS(II)
209	IF(II.GE.0) THEN
210	DELPLS(IAB)=DELPLS(IAB)-SOL(I)
211	ELSE
212	DELNEG(IAB)=DELNEG(IAB)-SOL(I)
213	ENDIF
214	SOL(I)=SOL(I)*0.25D+0
215	410 CONTINUE
216	GOTO 290
217	420 FORMAT(' FNORM, FOMX, FRAT', 3(1PD15.5))
218	430 FORMAT(' FGAM, RATGAM', 2(1PD15.5))
219	END