share/hompack/fortran/fixpqs.f

   1         SUBROUTINE FIXPQS(N,Y,IFLAG,ARCRE,ARCAE,ANSRE,ANSAE,TRACE,
   2      &     A,NFE,ARCLEN,YP,YOLD,YPOLD,QR,LENQR,PIVOT,PP,RHOVEC,Z0,DZ,
   3      &     T,WORK,SSPAR,PAR,IPAR)
   4 C
   5 C SUBROUTINE  FIXPQS  FINDS A FIXED POINT OR ZERO OF THE
   6 C N-DIMENSIONAL VECTOR FUNCTION  F(X), OR TRACKS A ZERO CURVE OF A
   7 C GENERAL HOMOTOPY MAP  RHO(A,X,LAMBDA).  FOR THE FIXED POINT PROBLEM
   8 C F(X) IS ASSUMED TO BE A C2 MAP OF SOME BALL INTO ITSELF.  THE
   9 C EQUATION  X=F(X)  IS SOLVED BY FOLLOWING THE ZERO CURVE OF THE
  10 C HOMOTOPY MAP
  11 C
  12 C  LAMBDA*(X - F(X)) + (1 - LAMBDA)*(X - A),
  13 C
  14 C STARTING FROM  LAMBDA = 0, X = A.   THE CURVE IS PARAMETERIZED
  15 C BY ARC LENGTH  S, AND IS FOLLOWED BY SOLVING THE ORDINARY
  16 C DIFFERENTIAL EQUATION  D(HOMOTOPY MAP)/DS = 0  FOR
  17 C Y(S) = (X(S),LAMBDA(S)).  THIS IS DONE BY USING A HERMITE CUBIC
  18 C PREDICTOR AND A CORRECTOR WHICH RETURNS TO THE ZERO CURVE IN A
  19 C HYPERPLANE PERPENDICULAR TO THE TANGENT TO THE ZERO CURVE AT THE
  20 C MOST RECENT POINT.
  21 C
  22 C FOR THE ZERO FINDING PROBLEM  F(X)  IS ASSUMED TO BE A C2 MAP SUCH
  23 C THAT FOR SOME  R > 0,  X*F(X) >= 0  WHENEVER  NORM(X) = R.
  24 C THE EQUATION  F(X) = 0  IS SOLVED BY FOLLOWING THE ZERO CURVE OF
  25 C THE HOMOTOPY MAP
  26 C
  27 C  LAMBDA*F(X) + (1 - LAMBDA)*(X - A)
  28 C
  29 C EMANATING FROM  LAMBDA = 0, X = A.
  30 C
  31 C A  MUST BE AN INTERIOR POINT OF THE ABOVE MENTIONED BALLS.
  32 C
  33 C FOR THE CURVE TRACKING PROBLEM RHO(A,X,LAMBDA) IS ASSUMED TO
  34 C BE A C2 MAP FROM  E**M X [0,1) X E**N  INTO  E**N, WHICH FOR
  35 C ALMOST ALL PARAMETER VECTORS  A  IN SOME NONEMPTY OPEN SUBSET
  36 C OF E**M SATISFIES
  37 C
  38 C  RANK [D RHO(A,X,LAMBDA)/D LAMBDA, D RHO(A,X,LAMBDA)/DX] = N
  39 C
  40 C FOR ALL POINTS  (X,LAMBDA)  SUCH THAT  RHO(A,X,LAMBDA) = 0.  IT IS
  41 C FURTHER ASSUMED THAT
  42 C
  43 C         RANK [ D RHO(A,X0,0)/DX ] = N.
  44 C
  45 C WITH  A  FIXED, THE ZERO CURVE OF  RHO(A,X,LAMBDA)  EMANATING FROM
  46 C LAMBDA = 0, X = X0  IS TRACKED UNTIL  LAMBDA = 1  BY SOLVING THE
  47 C ORDINARY DIFFERENTIAL EQUATION  D RHO(A,X(S),LAMBDA(S))/DS = 0
  48 C FOR  Y(S) = (X(S),LAMBDA(S)), WHERE  S  IS ARC LENGTH ALONG THE
  49 C ZERO CURVE.  ALSO THE HOMOTOPY MAP  RHO(A,X,LAMBDA)  IS ASSUMED TO
  50 C BE CONSTRUCTED SUCH THAT
  51 C
  52 C         D LAMBDA(0)/DS > 0.
  53 C
  54 C FOR THE FIXED POINT AND ZERO FINDING PROBLEMS, THE USER MUST SUPPLY
  55 C A SUBROUTINE  F(X,V)  WHICH EVALUATES F(X) AT X AND RETURNS THE
  56 C VECTOR F(X) IN  V,  AND A SUBROUTINE  FJACS(X,QR,LENQR,PIVOT)  WHICH
  57 C EVALUATES THE (SYMMETRIC) JACOBIAN MATRIX OF F(X) AT X, AND RETURNS
  58 C THE SYMMETRIC JACOBIAN MATRIX IN PACKED SKYLINE STORAGE FORMAT IN  QR.
  59 C LENQR  AND  PIVOT  DESCRIBE THE DATA STRUCTURE IN  QR.  FOR THE CURVE
  60 C TRACKING PROBLEM, THE USER MUST SUPPLY A SUBROUTINE
  61 C  RHO(A,LAMBDA,X,V,PAR,IPAR)  WHICH EVALUATES THE HOMOTOPY MAP RHO
  62 C AT (A,X,LAMBDA) AND RETURNS THE VECTOR RHO(A,X,LAMBDA) IN  V,
  63 C AND A SUBROUTINE  RHOJS(A,LAMBDA,X,QR,LENQR,PIVOT,PP,PAR,IPAR)  WHICH
  64 C RETURNS IN  QR  THE SYMMETRIC N X N JACOBIAN MATRIX [D RHO/DX]
  65 C EVALUATED AT (A,X,LAMBDA) AND STORED IN PACKED SKYLINE FORMAT,
  66 C AND RETURNS IN  PP  THE VECTOR  -(D RHO/D LAMBDA)  EVALUATED AT
  67 C (A,X,LAMBDA).  LENQR  AND  PIVOT  DESCRIBE THE DATA STRUCTURE IN
  68 C QR.
  69 C *** NOTE THE MINUS SIGN IN THE DEFINITION OF  PP. ***
  70 C
  71 C
  72 C FIXPQS DIRECTLY OR INDIRECTLY USES THE SUBROUTINES  D1MACH, F
  73 C (OR RHO), FJACS (OR RHOJS), GMFADS, MULTDS, PCGQS, ROOTQS, STEPQS,
  74 C SOLVDS,  AND THE BLAS ROUTINES  DAXPY, DCOPY, DDOT, DNRM2, AND  DSCAL.
  75 C ONLY  D1MACH  CONTAINS MACHINE DEPENDENT CONSTANTS.  NO OTHER
  76 C MODIFICATIONS BY THE USER ARE REQUIRED.
  77 C
  78 C
  79 C ON INPUT:
  80 C
  81 C N  IS THE DIMENSION OF X, F(X), AND RHO(A,X,LAMBDA).
  82 C
  83 C Y(1:N+1)  CONTAINS THE STARTING POINT FOR TRACKING THE HOMOTOPY MAP.
  84 C    (Y(1),...,Y(N)) = A  FOR THE FIXED POINT AND ZERO FINDING
  85 C    PROBLEMS.  (Y(1),...,Y(N)) = X0  FOR THE CURVE TRACKING PROBLEM.
  86 C    Y(N+1)  NEED NOT BE DEFINED BY THE USER.
  87 C
  88 C IFLAG  CAN BE -2, -1, 0, 2, OR 3.  IFLAG  SHOULD BE 0 ON THE FIRST
  89 C    CALL TO  FIXPQS  FOR THE PROBLEM  X=F(X), -1 FOR THE PROBLEM
  90 C    F(X)=0, AND -2 FOR THE PROBLEM  RHO(A,X,LAMBDA)=0.   IN CERTAIN
  91 C    SITUATIONS  IFLAG  IS SET TO 2 OR 3 BY  FIXPQS, AND  FIXPQS  CAN
  92 C    BE CALLED AGAIN WITHOUT CHANGING  IFLAG.
  93 C
  94 C ARCRE, ARCAE  ARE THE RELATIVE AND ABSOLUTE ERRORS, RESPECTIVELY,
  95 C    ALLOWED THE ITERATION ALONG THE ZERO CURVE.  IF
  96 C    ARC?E .LE. 0.0  ON INPUT, IT IS RESET TO  .5*SQRT(ANS?E).
  97 C    NORMALLY  ARC?E  SHOULD BE CONSIDERABLY LARGER THAN  ANS?E.
  98 C
  99 C ANSRE, ANSAE  ARE THE RELATIVE AND ABSOLUTE ERROR VALUES USED FOR
 100 C    THE ANSWER AT  LAMBDA = 1.  THE ACCEPTED ANSWER  Y = (X,LAMBDA)
 101 C    SATISFIES
 102 C
 103 C      |Y(1) - 1| .LE. ANSRE + ANSAE      .AND.
 104 C
 105 C      ||DZ|| .LE. ANSRE*||Y|| + ANSAE      WHERE
 106 C
 107 C      DZ IS THE NEWTON STEP TO Y.
 108 C
 109 C TRACE  IS AN INTEGER SPECIFYING THE LOGICAL I/O UNIT FOR
 110 C    INTERMEDIATE OUTPUT.  IF  TRACE .GT. 0  THE POINTS COMPUTED ON
 111 C    THE ZERO CURVE ARE WRITTEN TO I/O UNIT  TRACE .
 112 C
 113 C A(1:*)  CONTAINS THE PARAMETER VECTOR  A.  FOR THE FIXED POINT
 114 C    AND ZERO FINDING PROBLEMS,  A  NEED NOT BE INITIALIZED BY THE
 115 C    USER, AND IS ASSUMED TO HAVE LENGTH  N.  FOR THE CURVE
 116 C    TRACKING PROBLEM,  A  MUST BE INITIALIZED BY THE USER.
 117 C
 118 C YP(1:N+1)  IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO THE
 119 C    ZERO CURVE AT THE CURRENT POINT  Y.
 120 C
 121 C YOLD(1:N+1) IS A WORK ARRAY CONTAINING THE PREVIOUS POINT FOUND
 122 C    ON THE ZERO CURVE.
 123 C
 124 C YPOLD(1:N+1) IS A WORK ARRAY CONTAINING THE TANGENT VECTOR TO
 125 C    THE ZERO CURVE AT  YOLD.
 126 C
 127 C QR(1:LENQR)  IS A WORK ARRAY CONTAINING THE N X N SYMMETRIC
 128 C    JACOBIAN MATRIX WITH RESPECT TO X STORED IN PACKED SKYLINE
 129 C    STORAGE FORMAT.  LENQR  AND  PIVOT  DESCRIBE THE DATA
 130 C    STRUCTURE IN  QR.  (SEE SUBROUTINE  PCGQS  FOR A DESCRIPTION
 131 C    OF THIS DATA STRUCTURE).
 132 C
 133 C LENQR  IS THE LENGTH OF THE N-DIMENSIONAL ARRAY  QR.  I.E.
 134 C    IT IS THE NUMBER OF NON-ZERO ENTRIES IN THE JACOBIAN
 135 C    MATRIX [DF/DX] (OR [D RHO/DX]).
 136 C
 137 C PIVOT(1:N+2)  IS A WORK ARRAY WHOSE FIRST N+1 COMPONENTS CONTAIN
 138 C    THE INDICES OF THE DIAGONAL ELEMENTS OF THE N X N SYMMETRIC
 139 C    JACOBIAN MATRIX (WITH RESPECT TO X) WITHIN  QR.
 140 C
 141 C PP(1:N)  IS A WORK ARRAY CONTAINING THE NEGATIVE OF THE LAST COLUMN
 142 C    OF THE JACOBIAN MATRIX  -[D RHO/D LAMBDA].
 143 C
 144 C RHOVEC(1:N+1), Z0(1:N+1), DZ(1:N+1), T(1:N+1)  ARE ALL WORK ARRAYS
 145 C    USED BY  STEPQS, TANGQS, AND ROOTQS  TO CALCULATE THE TANGENT
 146 C    VECTORS AND NEWTON STEPS.
 147 C
 148 C WORK(1:8*(N+1)+LENQR)  IS A WORK ARRAY USED BY THE CONJUGATE GRADIENT
 149 C    ALGORITHM TO SOLVE LINEAR SYSTEMS.
 150 C
 151 C SSPAR(1:4) =  (HMIN, HMAX, BMIN, BMAX)  IS A VECTOR OF PARAMETERS
 152 C    USED FOR THE OPTIMAL STEP SIZE ESTIMATION.  A DEFAULT VALUE
 153 C    CAN BE SPECIFIED FOR ANY OF THESE FOUR PARAMETERS BY SETTING IT
 154 C    .LE. 0.0  ON INPUT.  SEE THE COMMENTS IN  STEPQS  FOR MORE
 155 C    INFORMATION ABOUT THESE PARAMETERS.
 156 C
 157 C PAR(1:*)  AND  IPAR(1:*)  ARE ARRAYS FOR (OPTIONAL) USER PARAMETERS,
 158 C    WHICH ARE SIMPLY PASSED THROUGH TO THE USER WRITTEN SUBROUTINES
 159 C    RHO, RHOJS.
 160 C
 161 C
 162 C ON OUTPUT:
 163 C
 164 C N , TRACE , A , LENQR  ARE UNCHANGED.
 165 C
 166 C Y(N+1) = LAMBDA, (Y(1),...,Y(N)) = X, AND  Y  IS AN APPROXIMATE
 167 C    ZERO OF THE HOMOTOPY MAP.  NORMALLY  LAMBDA = 1  AND  X  IS A
 168 C    FIXED POINT OR ZERO OF  F(X).   IN ABNORMAL SITUATIONS,  LAMBDA
 169 C    MAY ONLY BE NEAR 1 AND  X  NEAR A FIXED POINT OR ZERO.
 170 C
 171 C IFLAG =
 172 C
 173 C   1   NORMAL RETURN.
 174 C
 175 C   2   SPECIFIED ERROR TOLERANCE CANNOT BE MET.  SOME OR ALL OF
 176 C       ARCRE, ARCAE, ANSRE, ANSAE  HAVE BEEN INCREASED TO
 177 C       SUITABLE VALUES.  TO CONTINUE, JUST CALL  FIXPQS  AGAIN
 178 C       WITHOUT CHANGING ANY PARAMETERS.
 179 C
 180 C   3   STEPQS  HAS BEEN CALLED 1000 TIMES.  TO CONTINUE, CALL
 181 C       FIXPQS  AGAIN WITHOUT CHANGING ANY PARAMETERS.
 182 C
 183 C   4   JACOBIAN MATRIX DOES NOT HAVE FULL RANK.  THE ALGORITHM
 184 C       HAS FAILED (THE ZERO CURVE OF THE HOMOTOPY MAP CANNOT BE
 185 C       FOLLOWED ANY FURTHER).
 186 C
 187 C   5   THE TRACKING ALGORITHM HAS LOST THE ZERO CURVE OF THE
 188 C       HOMOTOPY MAP AND IS NOT MAKING PROGRESS.  THE ERROR
 189 C       TOLERANCES  ARC?E  AND  ANS?E  WERE TOO LENIENT.  THE PROBLEM
 190 C       SHOULD BE RESTRARTED BY CALLING  FIXPQS  WITH SMALLER ERROR
 191 C       TOLERANCES AND  IFLAG = 0 (-1, -2).
 192 C
 193 C   6   THE NEWTON ITERATION IN  STEPQS  OR  ROOTQS  FAILED TO
 194 C       CONVERGE.  THE ERROR TOLERANCES  ANS?E  MAY BE TOO STRINGENT.
 195 C
 196 C   7   ILLEGAL INPUT PARAMETERS, A FATAL ERROR.
 197 C
 198 C ARCRE, ARCAE, ANSRE, ANSAE  ARE UNCHANGED AFTER A NORMAL RETURN
 199 C    (IFLAG = 1).  THEY ARE INCREASED TO APPROPRIATE VALUES ON THE
 200 C    RETURN  IFLAG = 2.
 201 C
 202 C NFE  IS THE NUMBER OF JACOBIAN EVALUATIONS.
 203 C
 204 C ARCLEN  IS THE APPROXIMATE LENGTH OF THE ZERO CURVE.
 205 C
 206 C ***** DECLARATIONS *****
 207 C
 208 C     FUNCTION DECLARATIONS
 209 C
 210         DOUBLE PRECISION D1MACH, DNRM2
 211 C
 212 C     LOCAL VARIABLES
 213 C
 214         DOUBLE PRECISION ABSERR, H, HOLD, RELERR, S, WK
 215         INTEGER IFLAGC, ITER, JW, LIMITD, LIMIT, NP1, PCGWK
 216         LOGICAL CRASH, START
 217 C
 218 C     SCALAR ARGUMENTS
 219 C
 220         DOUBLE PRECISION ARCRE, ARCAE, ANSRE, ANSAE, ARCLEN
 221         INTEGER N, IFLAG, TRACE, NFE, LENQR
 222 C
 223 C     ARRAY DECLARATIONS
 224 C
 225         DOUBLE PRECISION A(N), Y(N+1), YP(N+1), YOLD(N+1), YPOLD(N+1),
 226      &    QR(LENQR), PP(N), RHOVEC(N+1), Z0(N+1), DZ(N+1), T(N+1),
 227      &    WORK(8*(N+1)+LENQR), SSPAR(4), PAR(1)
 228         INTEGER PIVOT(N+2), IPAR(1)
 229 C
 230         SAVE
 231 C
 232 C ***** END OF DECLARATIONS *****
 233 C
 234 C LIMITD  IS AN UPPER BOUND ON THE NUMBER OF STEPS.  IT MAY BE
 235 C CHANGED BY CHANGING THE FOLLOWING PARAMETER STATEMENT:
 236         PARAMETER (LIMITD =1000)
 237 C
 238 C ***** FIRST EXECUTABLE STATEMENT *****
 239 C
 240 C CHECK IFLAG
 241 C
 242         IF (N .LE. 0 .OR. ANSRE .LE. 0.0 .OR. ANSAE .LT. 0.0)
 243      &    IFLAG = 7
 244         IF (IFLAG .GE. -2 .AND. IFLAG .LE. 0) GO TO 10
 245         IF (IFLAG .EQ. 2) GO TO 50
 246         IF (IFLAG .EQ. 3) GO TO 40
 247 C
 248 C ONLY VALID INPUT FOR IFLAG IS -2, -1, 0, 2, 3.
 249 C
 250         IFLAG = 7
 251         RETURN
 252 C
 253 C ***** INITIALIZATION BLOCK  *****
 254 C
 255   10    ARCLEN = 0.0
 256         IF (ARCRE .LE. 0.0) ARCRE = .5*SQRT(ANSRE)
 257         IF (ARCAE .LE. 0.0) ARCAE = .5*SQRT(ANSAE)
 258         NFE=0
 259         IFLAGC = IFLAG
 260         NP1=N+1
 261         PCGWK = 2*N+3
 262 C
 263 C SET INITIAL CONDITIONS FOR FIRST CALL TO STEPQS.
 264 C
 265         START=.TRUE.
 266         CRASH=.FALSE.
 267         RELERR = ARCRE
 268         ABSERR = ARCAE
 269         HOLD=1.0
 270         H=0.1
 271         S=0.0
 272         YPOLD(NP1) = 1.0
 273         Y(NP1) = 0.0
 274         DO 20 JW=1,N
 275           YPOLD(JW)=0.0
 276   20    CONTINUE
 277 C
 278 C SET OPTIMAL STEP SIZE ESTIMATION PARAMETERS.
 279 C
 280 C     MINIMUM STEP SIZE HMIN
 281         IF (SSPAR(1) .LE. 0.0) SSPAR(1)= (SQRT(N+1.0)+4.0)*D1MACH(4)
 282 C     MAXIMUM STEP SIZE HMAX
 283         IF (SSPAR(2) .LE. 0.0) SSPAR(2)= 1.0
 284 C     MINIMUM STEP REDUCTION FACTOR BMIN
 285         IF (SSPAR(3) .LE. 0.0) SSPAR(3)= 0.1
 286 C     MAXIMUM STEP EXPANSION FACTOR BMAX
 287         IF (SSPAR(4) .LE. 0.0) SSPAR(4)= 7.0
 288 C
 289 C LOAD  A  FOR THE FIXED POINT AND ZERO FINDING PROBLEMS.
 290 C
 291         IF (IFLAGC .GE. -1) THEN
 292           CALL DCOPY(N,Y,1,A,1)
 293         ENDIF
 294 C
 295   40    LIMIT=LIMITD
 296 C
 297 C ***** END OF INITIALIZATION BLOCK. *****
 298 C
 299 C ***** MAIN LOOP. *****
 300 C
 301   50    DO 400 ITER=1,LIMIT
 302         IF (Y(NP1) .LT. 0.0) THEN
 303           ARCLEN = S
 304           IFLAG = 5
 305           RETURN
 306         END IF
 307 C
 308 C TAKE A STEP ALONG THE CURVE.
 309 C
 310         CALL STEPQS(N,NFE,IFLAGC,LENQR,START,CRASH,HOLD,H,WK,RELERR,
 311      &    ABSERR,S,Y,YP,YOLD,YPOLD,A,QR,PIVOT,PP,RHOVEC,Z0,DZ,T,
 312      &    WORK,SSPAR,PAR,IPAR)
 313 C
 314 C PRINT LATEST POINT ON CURVE IF REQUESTED.
 315 C
 316       IF (TRACE .GT. 0) THEN
 317         WRITE (TRACE,217) ITER,NFE,S,Y(NP1),(Y(JW),JW=1,N)
 318 217     FORMAT(/' STEP',I5,3X,'NFE =',I5,3X,'ARC LENGTH =',F9.4,3X,
 319      &  'LAMBDA =',F7.4,5X,'X vector:'/1P,(1X,6E12.4))
 320       ENDIF
 321 C
 322 C CHECK IF THE STEP WAS SUCCESSFUL.
 323 C
 324         IF (IFLAGC .GT. 0) THEN
 325           ARCLEN=S
 326           IFLAG=IFLAGC
 327           RETURN
 328         END IF
 329 C
 330         IF (CRASH) THEN
 331 C
 332 C         RETURN CODE FOR ERROR TOLERANCE TOO SMALL.
 333 C
 334           IFLAG=2
 335 C
 336 C         CHANGE ERROR TOLERANCES.
 337 C
 338           IF (ARCRE .LT. RELERR) THEN
 339             ARCRE=RELERR
 340             ANSRE=RELERR
 341           ENDIF
 342           IF (ARCAE .LT. ABSERR) ARCAE=ABSERR
 343 C
 344 C         CHANGE LIMIT ON NUMBER OF ITERATIONS.
 345 C
 346           LIMIT = LIMIT - ITER
 347           RETURN
 348         END IF
 349 C
 350 C IF  LAMBDA >= 1.0,  USE  ROOTQS  TO FIND SOLUTION.
 351 C
 352         IF (Y(NP1) .GE. 1.0) GOTO 500
 353 C
 354   400   CONTINUE
 355 C
 356 C ***** END OF MAIN LOOP *****
 357 C
 358 C DID NOT CONVERGE IN  LIMIT  ITERATIONS, SET  IFLAG  AND RETURN.
 359 C
 360         ARCLEN = S
 361         IFLAG = 3
 362         RETURN
 363 C
 364 C ***** FINAL STEP -- FIND SOLUTION AT LAMBDA=1 *****
 365 C
 366 C SAVE  YOLD  FOR ARC LENGTH CALCULATION LATER.
 367 C
 368   500   CALL DCOPY(NP1,YOLD,1,T,1)
 369 C
 370 C FIND SOLUTION.
 371 C
 372         CALL ROOTQS(N,NFE,IFLAGC,LENQR,ANSRE,ANSAE,Y,YP,YOLD,YPOLD,
 373      &    A,QR,PIVOT,PP,RHOVEC,Z0,DZ,WORK(PCGWK),PAR,IPAR)
 374 C
 375 C CHECK IF SOLUTION WAS FOUND AND SET  IFLAG  ACCORDINGLY.
 376 C
 377         IFLAG=1
 378 C
 379 C     SET ERROR FLAG IF ROOTQS COULD NOT GET THE POINT ON THE ZERO
 380 C     CURVE AT  LAMBDA = 1.0 .
 381 C
 382         IF (IFLAGC .GT. 0) IFLAG=IFLAGC
 383 C
 384 C CALCULATE FINAL ARC LENGTH.
 385 C
 386         CALL DCOPY(NP1,Y,1,DZ,1)
 387         WK=-1.0
 388         CALL DAXPY(NP1,WK,T,1,DZ,1)
 389         ARCLEN=S - HOLD + DNRM2(NP1,DZ,1)
 390 C
 391 C ***** END OF FINAL STEP *****
 392 C
 393         RETURN
 394 C
 395 C ***** END OF SUBROUTINE FIXPQS *****
 396         END