MINICERN/mathlib/gen/c/snleq64.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1996/04/01 15:01:52  mclareni
   6 * Mathlib gen
   7 *
   8 *
   9 #include "gen/pilot.h"
  10 #if defined(CERNLIB_DOUBLE)
  11       SUBROUTINE DSNLEQ(N,X,F,FTOL,XTOL,MAXF,IPRT,INFO,SUB,W)
  12 #include "gen/imp64.inc"
  13 #endif
  14 #if !defined(CERNLIB_DOUBLE)
  15       SUBROUTINE  RSNLEQ(N,X,F,FTOL,XTOL,MAXF,IPRT,INFO,SUB,W)
  16 #endif
  17
  18 C     Based on   J.J. More  and  M.Y. Cosnard
  19 C
  20 C       ALGORITHM 554 BRENTM, A Fortran Subroutine for the
  21 C       Numerical Solution of Systems of Nonlinear Equations [C5]
  22 C
  23 C     ACM Trans. Math. Software 6 (1980) 240-251.
  24
  25       DIMENSION X(N),F(N),W(N,*),MPT(288)
  26       LOGICAL LCV
  27
  28       PARAMETER (Z1 = 1, SCALE = 10, P05 = 5*Z1/100)
  29 C**** EPS = SQRT(SMALLEST FP.NUMBER)
  30 C     EPS = 1 / SQRT( 16D0**13 )
  31 #if !defined(CERNLIB_DOUBLE)
  32       PARAMETER (EPS =  0.84293 69702 17878 97282 52636 392E-07)
  33 #endif
  34 #if defined(CERNLIB_IBM)
  35       PARAMETER (EPS =  0.14901 16119 38476 562D-07)
  36 #endif
  37 #if defined(CERNLIB_VAX)
  38       PARAMETER (EPS =  0.37252 90298 46191 40625D-08)
  39 #endif
  40 #if (defined(CERNLIB_UNIX))&&(defined(CERNLIB_DOUBLE))
  41       PARAMETER (EPS =  0.14901 16119 38476 600D-07)
  42 #endif
  43       DATA (MPT(I),I=1,288)
  44      1/1* 1,1* 2,3* 3,3* 4,4* 5,4* 6,4* 7,4* 8,5* 9,5*10,5*11,5*12,
  45      2 5*13,5*14,6*15,6*16,5*17,6*18,6*19,6*20,7*21,6*22,6*23,7*24,
  46      3 6*25,7*26,6*27,7*28,7*29,7*30,7*31,7*32,7*33,7*34,7*35,7*36,
  47      4 8*37,7*38,7*39,8*40,7*41,8*42,7*43,8*44,8*45,7*46,8*47,8*48/
  48
  49       INFO=0
  50       IF(N .LE. 0 .OR. FTOL .LE. 0 .OR. XTOL .LE. 0) RETURN
  51 C
  52 C     Find optimal MOPT for iterative refinement
  53 C
  54       IF(N .LE. 288) THEN
  55        MOPT=MPT(N)
  56       ELSE
  57        H=0
  58        DO 1 I = 49,N
  59        TEMP=LOG(I+Z1)/(N+2*I+1)
  60        IF(TEMP .LT. H) THEN
  61         MOPT=I-1
  62         GO TO 2
  63        ENDIF
  64     1  H=TEMP
  65       ENDIF
  66
  67     2 IFLAG=0
  68       NUMF=0
  69       NFCALL=0
  70
  71       NIER6=-1
  72       NIER7=-1
  73       NIER8=0
  74       FNORM=0
  75       DIFIT=0
  76       XNORM=0
  77       DO 10 I = 1,N
  78    10 XNORM=MAX(XNORM,ABS(X(I)))
  79       DELTA=SCALE*XNORM
  80       IF(XNORM .EQ. 0) DELTA=SCALE
  81
  82    20 IF(IPRT .NE. 0) WRITE(6,'(1X,I5,D25.14)') (I,X(I),I=1,N)
  83
  84       NSING=N
  85       FNORM1=FNORM
  86       DIFIT1=DIFIT
  87       FNORM=0
  88 C
  89 C     Compute step H for the divided difference which approximates
  90 C     the K-th row of the Jacobian matrix
  91 C
  92       H=EPS*XNORM
  93       IF(H .EQ. 0) H=EPS
  94       DO 40 J = 1,N
  95       DO 30 I = 1,N
  96    30 W(I,J+3)=0
  97       W(J,J+3)=H
  98    40 W(J,2)=X(J)
  99 C
 100 C     Enter a subiteration
 101 C
 102       DO 150 K = 1,N
 103       IFLAG=K
 104       CALL SUB(N,W(1,2),F,IFLAG)
 105       FKY=F(K)
 106       NFCALL=NFCALL+1
 107       NUMF=NFCALL/N
 108       IF(IFLAG .LT. 0) GO TO 230
 109       FNORM=MAX(FNORM,ABS(FKY))
 110 C
 111 C     Compute the K-th row of the Jacobian matrix
 112 C
 113       DO 60 J = K,N
 114       DO 50 I = 1,N
 115    50 W(I,3)=W(I,2)+W(I,J+3)
 116       CALL SUB(N,W(1,3),F,IFLAG)
 117       FKZ=F(K)
 118       NFCALL=NFCALL+1
 119       NUMF=NFCALL/N
 120       IF(IFLAG .LT. 0) GO TO 230
 121    60 W(J,1)=FKZ-FKY
 122       F(K)=FKY
 123 C
 124 C     Compute the Householder transformation to reduce the K-th row
 125 C     of the Jacobian matrix to a multiple of the K-th unit vector
 126 C
 127       ETA=0
 128       DO 70 I = K,N
 129    70 ETA=MAX(ETA,ABS(W(I,1)))
 130       IF(ETA .EQ. 0) GO TO 150
 131       NSING=NSING-1
 132       SKNORM=0
 133       DO 80 I = K,N
 134       W(I,1)=W(I,1)/ETA
 135    80 SKNORM=SKNORM+W(I,1)**2
 136       SKNORM=SQRT(SKNORM)
 137       IF(W(K,1) .LT. 0) SKNORM=-SKNORM
 138       W(K,1)=W(K,1)+SKNORM
 139 C
 140 C     Apply the transformation
 141 C
 142       DO 90 I = 1,N
 143    90 W(I,3)=0
 144       DO 100 J = K,N
 145       DO 100 I = 1,N
 146   100 W(I,3)=W(I,3)+W(J,1)*W(I,J+3)
 147       DO 120 J = K,N
 148       TEMP=W(J,1)/(SKNORM*W(K,1))
 149       DO 120 I = 1,N
 150   120 W(I,J+3)=W(I,J+3)-TEMP*W(I,3)
 151 C
 152 C     Compute the subiterate
 153 C
 154       W(K,1)=SKNORM*ETA
 155       TEMP=FKY/W(K,1)
 156       IF(H*ABS(TEMP) .GT. DELTA) TEMP=SIGN(DELTA/H,TEMP)
 157       DO 140 I = 1,N
 158   140 W(I,2)=W(I,2)+TEMP*W(I,K+3)
 159   150 CONTINUE
 160 C
 161 C     Compute the norms of the iterate and correction vector
 162 C
 163       XNORM=0
 164       DIFIT=0
 165       DO 160 I = 1,N
 166       XNORM=MAX(XNORM,ABS(W(I,2)))
 167       DIFIT=MAX(DIFIT,ABS(X(I)-W(I,2)))
 168   160 X(I)=W(I,2)
 169 C
 170 C     Update the bound on the correction vector
 171 C
 172       DELTA=MAX(DELTA,SCALE*XNORM)
 173 C
 174 C     Determine the progress of the iteration
 175 C
 176       LCV=FNORM .LT. FNORM1 .AND. DIFIT .LT. DIFIT1 .AND. NSING .EQ. 0
 177       NIER6=NIER6+1
 178       NIER7=NIER7+1
 179       NIER8=NIER8+1
 180       IF(LCV) NIER6=0
 181       IF(FNORM .LT. FNORM1 .OR. DIFIT .LT. DIFIT1) NIER7=0
 182       IF(DIFIT .GT. EPS*XNORM) NIER8=0
 183 C
 184 C     Tests for convergence
 185 C
 186       IF(FNORM .LE. FTOL) INFO=1
 187       IF(DIFIT .LE. XTOL*XNORM .AND. LCV) INFO=2
 188       IF(FNORM .LE. FTOL .AND. INFO .EQ. 2) INFO=3
 189       IF(INFO .NE. 0) GO TO 230
 190 C
 191 C     Tests for termination
 192 C
 193       IF(NUMF .GE. MAXF) INFO=4
 194       IF(NSING .EQ. N) INFO=5
 195       IF(NIER6 .EQ. 5) INFO=6
 196       IF(NIER7 .EQ. 3) INFO=7
 197       IF(NIER8 .EQ. 4) INFO=8
 198       IF(INFO .NE. 0) GO TO 230
 199       IF(.NOT.LCV .OR. DIFIT .GT. P05*XNORM) GO TO 20
 200 C
 201 C     Iterative refinement  (if the iteration is converging)
 202 C
 203       DO 210 M = 2,MOPT
 204       FNORM1=FNORM
 205       FNORM=0
 206       DO 190 K = 1,N
 207       IFLAG=K
 208       CALL SUB(N,W(1,2),F,IFLAG)
 209       FKY=F(K)
 210       NFCALL=NFCALL+1
 211       NUMF=NFCALL/N
 212       IF(IFLAG .LT. 0) GO TO 230
 213       FNORM=MAX(FNORM,ABS(FKY))
 214 C
 215 C     Iterative refinement is terminated if it does not give a
 216 C     reduction on residuals
 217 C
 218       IF(FNORM .GE. FNORM1) THEN
 219        FNORM=FNORM1
 220        GO TO 20
 221       ENDIF
 222       TEMP=FKY/W(K,1)
 223       DO 180 I = 1,N
 224   180 W(I,2)=W(I,2)+TEMP*W(I,K+3)
 225   190 CONTINUE
 226 C
 227 C     Compute the norms of the iterate and correction vector
 228 C
 229       XNORM=0
 230       DIFIT=0
 231       DO 200 I = 1,N
 232       XNORM=MAX(XNORM,ABS(W(I,2)))
 233       DIFIT=MAX(DIFIT,ABS(X(I)-W(I,2)))
 234   200 X(I)=W(I,2)
 235 C
 236 C     Stopping criteria for iterative refinement
 237 C
 238       IF(FNORM .LE. FTOL) INFO=1
 239       IF(DIFIT .LE. XTOL*XNORM) INFO=2
 240       IF(FNORM .LE. FTOL .AND. INFO .EQ. 2) INFO=3
 241       IF(NUMF .GE. MAXF .AND. INFO .EQ. 0) INFO=4
 242       IF(INFO .NE. 0) GO TO 230
 243   210 CONTINUE
 244       GO TO 20
 245
 246   230 IF(IFLAG .LT. 0) INFO=IFLAG
 247       RETURN
 248       END