--- /dev/null
+*CMZ : 26/11/93 11.31.36 by Jonathan Butterworth
+*-- Author :
+ 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, NEAR
+ 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/
+ SAVE
+* write(*,*) 'SOL',sol
+* write(*,*) 'FNLIN',fnlin
+ NFCNT=0
+ FDIF=FMAJOR-FMINOR
+ PROD=1.0D+0
+ DO 10 I=1,N
+C ** JMB
+ IF (ABS(DELPLS(I)).LT.1.D-16) DELPLS(I)=0.D0
+ IF (ABS(DELNEG(I)).LT.1.D-16) DELNEG(I)=0.D0
+C -- JMB
+ PROD=PROD*(DELPLS(I)+DELNEG(I))
+ 10 CONTINUE
+C ** JMB
+ IF (VOL.NE.0.D0) THEN
+ GAMMA=PROD/VOL
+ FGAM=GAMMA*FMAJOR+(1.0D+0-GAMMA)*FMINOR
+ ELSE
+ GAMMA=1.0D0
+ FGAM = 0.0D0
+ ENDIF
+C -- JMB
+ 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
+ 60 YDMIN= 1.0D+30
+ NEAR=0
+ DO 70 I=1,NCUT
+ YDI=ABS(FMAJOR-FNEW(I))
+ IF(YDI.GT.YDMIN) GOTO 70
+ YDMIN=YDI
+ NEAR=I
+ 70 CONTINUE
+ IF(NEAR.EQ.0) RETURN
+ IF((LMAX.AND.FNEW(NEAR).LT.FGAM).OR.
+ 1(.NOT.LMAX.AND.FNEW(NEAR).GT.FGAM)) GOTO 130
+ ISAVE=ICUT(NEAR)
+ IF(NEAR.EQ.NCUT.OR.NCUT.EQ.1) GOTO 90
+ NCUTM1=NCUT-1
+ DO 80 I=NEAR,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)
+* write(*,*) 'A:delpls(isave),isave,xup(isave),x(isave)'
+* & ,delpls(isave),isave,xup(isave),x(isave)
+ GOTO 110
+ 100 IAB=ABS(ISAVE)
+ DELNEG(IAB)=X(IAB)-XLOW(IAB)
+* write(*,*) 'delneg(iab),iab,x(iab),xlow(iab)',delneg(iab),iab
+* & ,x(iab),xlow(iab)
+ 110 PROD=1.0D+0
+ DO 120 I=1,N
+C ** JMB
+ IF (ABS(DELPLS(I)).LT.1.D-16) DELPLS(I)=0.D0
+ IF (ABS(DELNEG(I)).LT.1.D-16) DELNEG(I)=0.D0
+C -- JMB
+ PROD=PROD*(DELPLS(I)+DELNEG(I))
+ 120 CONTINUE
+C ** JMB
+ IF (VOL.NE.0.D0) THEN
+ GAMMA=PROD/VOL
+ FGAM=GAMMA*FMAJOR+(1.0D+0-GAMMA)*FMINOR
+ ELSE
+ GAMMA=1.0D0
+ FGAM = 0.0D0
+ ENDIF
+C -- JMB
+ 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)
+ IF (DEL.NE.0.D0) THEN
+ DFORIG(I)=(FNEW(I)-FMAJOR)/DEL
+ ELSE
+ DFORIG(I)=0.D0
+ ENDIF
+ 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)
+ IF (GAMMA.NE.0) THEN
+ IF (ABS(DELPLS(II)+DELNEG(II)).NE.0.D0) THEN
+ FNLROW(I)=-GAMMA*FDIF/(DELPLS(II)+DELNEG(II))
+ ELSE
+* write(*,440) 'DELPLS(II)+DELNEG(II)='
+* write(*,*) DELPLS(II)+DELNEG(II)
+ FNLROW(I)=0
+ ENDIF
+ ELSE
+ FNLROW(I)=0.0
+ ENDIF
+ 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))) THEN
+ IF (DIAGJ(I).NE.0.D0) THEN
+ XMULT=FNLROW(I)/DIAGJ(I)
+ ELSE
+* write(*,440) 'diag(j)=0!'
+ XMULT=0.D0
+ ENDIF
+ ENDIF
+ 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))) THEN
+ IF (FNLROW(NCUT).NE.0.D0) THEN
+ SOL(NCUT)=FNLIN(NCUT)/FNLROW(NCUT)
+ ELSE
+ SOL(NCUT)=0.D0
+ ENDIF
+ ENDIF
+* write(*,*) 'FNLIN(NCUT),ncut,fnlrow(ncut)',FNLIN(NCUT),ncut
+* & ,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
+* write(*,*) 'sol(iback),iback,val',sol(iback),iback,val
+ IF(ABS(DIAGJ(IBACK)).LT.BIG*ABS(VAL)) THEN
+ IF (DIAGJ(IBACK).NE.0.D0) THEN
+ SOL(IBACK)=VAL/DIAGJ(IBACK)
+ ELSE
+ SOL(IBACK)=0.D0
+ ENDIF
+ ENDIF
+ 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)
+* write(*,*) 'A:delpls(iab),iab,sol(i),i,delmax,delmin,',delpls(iab)
+* & ,iab,sol(i),i,delmax,delmin
+ GOTO 310
+ 300 DELNEG(IAB)=DELNEG(IAB)+SOL(I)
+* write(*,*) 'A:delneg(iab),iab,sol(i),i',delneg(iab),iab
+* & ,sol(i),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 CONTINUE
+* write(*,*) 'iab,x(iab),z(iab),delpls,delneg',iab,x(iab),z(iab)
+* & ,delpls,delneg
+ 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
+C ** JMB
+ IF (ABS(DELPLS(I)).LT.1.D-16) DELPLS(I)=0.D0
+ IF (ABS(DELNEG(I)).LT.1.D-16) DELNEG(I)=0.D0
+C -- JMB
+ PROD=PROD*(DELPLS(I)+DELNEG(I))
+ 350 CONTINUE
+ IF (VOL.NE.0.D0) THEN
+ GAMNEW=PROD/VOL
+ ELSE
+ GAMNEW=0.D0
+ ENDIF
+ 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)) THEN
+ IF (SOL(I).NE.0.D0) THEN
+ TSTVAL=DFNEW/SOL(I)
+ ELSE
+ TSTVAL=0.D0
+ ENDIF
+ ENDIF
+ 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)
+* write(*,*) 'B:delpls(iab),iab,sol(i),i',delpls(iab),iab
+* & ,sol(i),i
+ ELSE
+ DELNEG(IAB)=DELNEG(IAB)-SOL(I)
+* write(*,*) 'B:delneg(iab),iab,sol(i),i',delneg(iab),iab
+* & ,sol(i),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))
+ 440 FORMAT(A)
+ END
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff