share/hompack/fortran/pcgqs.f

   1       SUBROUTINE PCGQS(NN,AA,LENAA,MAXA,PP,YP,RHO,START,WORK,IFLAG)
   2 C
   3 C THIS SUBROUTINE SOLVES A SYSTEM OF EQUATION USING THE METHOD OF
   4 C CONJUGATE GRADIENTS.  THE SYSTEM TO BE SOLVED IS IN THE FORM
   5 C
   6 C   (AUG)*X = B,
   7 C
   8 C WHERE
   9 C
  10 C            +--          --+        +-    -+
  11 C            |        |     |        |      |
  12 C            |   AA   | -PP |        |      |
  13 C      AUG = |        |     | ,  B = | -RHO |
  14 C            +--------------+        |      |
  15 C            |      YP      |        |      |
  16 C            +--          --+        +-    -+
  17 C
  18 C
  19 C
  20 C THE SYSTEM IS SOLVED BY SPLITTING  AUG  INTO TWO MATRICES
  21 C AUG = M + L,  WHERE
  22 C
  23 C       +-        -+                                +-     -+
  24 C       |      |   |                                |       |
  25 C       |  AA  | C |                                | -PP-C |
  26 C   M = |      |   |  ,  L = U * [E(NN+1)**T],  U = |       | ,
  27 C       +------+---+                                +-------+
  28 C       |  C   | D |                                |   0   |
  29 C       +-        -+                                +-     -+
  30 C
  31 C E(NN+1) IS THE (NN+1) X 1 VECTOR CONSISTING OF ALL ZEROS EXCEPT FOR
  32 C A '1' IN ITS LAST POSITION.
  33 C
  34 C THE FINAL SOLUTION VECTOR,  X,  IS GIVEN BY
  35 C
  36 C            +-                                    -+
  37 C            |           [SOL(U)]*[E(NN+1)**T]      |
  38 C       X =  | I  -  -----------------------------  | * SOL(B)
  39 C            |        {[(SOL(U))**T]*E(NN+1)}+1.0   |
  40 C            +-                                    -+
  41 C
  42 C WHERE SOL(A)=[M**(-1)]*A.  THE TWO SYSTEMS  (MZ=U, MZ=B) ARE SOLVED
  43 C BY A PRECONDITIONED CONJUGATE GRADIENT ALGORITHM.
  44 C
  45 C
  46 C ON INPUT:
  47 C
  48 C NN  = THE DIMENSION OF THE MATRIX PACKED IN AA.
  49 C
  50 C AA(1:LENAA)  CONTAINS THE MATRIX  AA,  STORED IN PACKED SKYLINE
  51 C    FORMAT.  LENAA  AND  MAXA  DESCRIBE THE DATA STRUCTURE.
  52 C
  53 C LENAA  = THE LENGTH OF THE ONE-DIMENSIONAL ARRAY  AA.
  54 C
  55 C MAXA(1:NN+2)  IN ITS FIRST  N+1  COMPONENTS CONTAINS THE INDICES OF
  56 C    THE DIAGONAL ELEMENTS OF THE MATRIX STORED IN  AA.
  57 C    MAXA(NN+2)  IS ASSIGNED THE VALUE  LENNAA + NN + 2.
  58 C
  59 C    AS AN EXAMPLE OF THE PACKED SKYLINE STORAGE FORMAT, CONSIDER THE
  60 C    SYMMETRIC MATRIX
  61 C
  62 C
  63 C                +--             --+
  64 C                |  1  3  0  0  0  |
  65 C                |  3  2  0  7  0  |
  66 C                |  0  0  4  6  0  |
  67 C                |  0  7  6  5  9  |
  68 C                |  0  0  0  9  8  |
  69 C                +--             --+
  70 C
  71 C    THIS WOULD RESULT IN  NN=5, LENAA=9, MAXA=(1,2,4,5,8,10,16),
  72 C    AND  AA=(1,2,3,4,5,6,7,8,9).
  73 C
  74 C PP(1:NN)  = THE NEGATIVE OF THE LAST COLUMN OF  AUG.
  75 C
  76 C YP(1:NN+1)  = THE LAST ROW OF  AUG.
  77 C
  78 C RHO(1:NN+1)  = THE NEGATIVE OF THE RIGHT HAND SIDE OF THE EQUATION TO
  79 C    BE SOLVED.
  80 C
  81 C WORK(1:6*(NN+1)+LENAA)  IS A WORK ARRAY DIVIDED UP AS FOLLOWS:
  82 C
  83 C    WORK(1:NN+1) = TEMPORARY WORKING STORAGE.
  84 C
  85 C    WORK(NN+2:2*NN+2) = INTERMEDIATE SOLUTION VECTOR Z FOR  MZ=U.
  86 C       THE INPUT VALUE IS USED AS THE INITIAL ESTIMATE FOR Z.
  87 C
  88 C    WORK(2*NN+3:3*NN+3) = INTERMEDIATE SOLUTION VECTOR Z FOR  MZ=B.
  89 C
  90 C    WORK(3*NN+4:4*NN+4) = STORAGE FOR THE RESIDUAL VECTORS.
  91 C
  92 C    WORK(4*NN+5:5*NN+5) = STORAGE FOR THE DIRECTION VECTORS.
  93 C
  94 C    WORK(5*NN+6:6*NN+6+LENAA) = STORAGE FOR THE PRECONDITIONING
  95 C       MATRIX Q.
  96 C
  97 C
  98 C ON OUTPUT:
  99 C
 100 C NN, AA, LENAA, MAXA, PP, YP, AND RHO  ARE UNCHANGED.
 101 C
 102 C START(1:N+1)  CONTAINS THE SOLUTION VECTOR  X.
 103 C
 104 C IFLAG  IS UNCHANGED UNLESS THE CONJUGATE GRADIENT ITERATION
 105 C    FAILS TO CONVERGE IN 10*(NN+1) ITERATIONS (MOST LIKELY DUE
 106 C    TO A SINGULAR JACOBIAN MATRIX).  IN THIS CASE,  PCGQS  RETURNS
 107 C    IFLAG = 4, AND DOES NOT COMPUTE  X.
 108 C
 109 C
 110 C CALLS  D1MACH, DAXPY, DCOPY, DDOT, DNRM2, DSCAL, GMFADS, MULTDS,
 111 C    SOLVDS.
 112 C
 113 C ***** DECLARATIONS *****
 114 C
 115 C     FUNCTION DECLARATIONS
 116 C
 117         DOUBLE PRECISION D1MACH, DDOT, DNRM2
 118 C
 119 C     LOCAL VARIABLES
 120 C
 121         DOUBLE PRECISION AB, AU, BB, BU, DZNRM, PBNPRD, PUNPRD,
 122      $     RBNPRD, RBTOL, RNPRD, RUNPRD, RUTOL, TEMP, ZLEN, ZTOL
 123         INTEGER DIR, IMAX, J, LENQ, NP1, Q, RES, ZB, ZU
 124         LOGICAL STILLU, STILLB
 125 C
 126 C     SCALAR ARGUMENTS
 127 C
 128         INTEGER NN, LENAA, IFLAG
 129 C
 130 C     ARRAY DECLARATIONS
 131 C
 132         DOUBLE PRECISION AA(LENAA), PP(NN), YP(NN+1), RHO(NN+1),
 133      $     START(NN+1), WORK(6*(NN+1)+LENAA)
 134         INTEGER MAXA(NN+2)
 135 C
 136 C ***** END OF DECLARATIONS *****
 137 C
 138 C ***** FIRST EXECUTABLE STATEMENT *****
 139 C
 140 C SET UP BASES FOR VECTORS STORED IN WORK ARRAY.
 141 C
 142       NP1=NN+1
 143       ZU=NN+2
 144       ZB=(2*NN)+3
 145       RES=(3*NN)+4
 146       DIR=(4*NN)+5
 147       Q=(5*NN)+6
 148 C
 149 C INITIALIZE  PRECONDITIONING MATRIX  Q,  SET VALUES OF  MAXA(NN+1)
 150 C AND  MAXA(NN+2),  COMPUTE PRECONDITIONER.
 151 C
 152       CALL DCOPY(LENAA,AA,1,WORK(Q),1)
 153       CALL DCOPY(NP1,YP,1,WORK(Q+LENAA),1)
 154       MAXA(NN+1)=LENAA+1
 155       MAXA(NN+2)=LENAA+NN+2
 156       LENQ = MAXA(NN+2)-1
 157       CALL GMFADS(NP1,WORK(Q),LENQ,MAXA)
 158 C
 159 C COMPUTE ALL TOLERANCES NEEDED FOR EXIT CRITERIA.
 160 C
 161       CALL DCOPY(NN,PP,1,WORK,1)
 162       WORK(NP1)=0.0
 163       CALL DAXPY(NP1,1.0D0,YP,1,WORK,1)
 164       IMAX=10*NP1
 165       STILLU=.TRUE.
 166       STILLB=.TRUE.
 167       ZTOL=100.0*D1MACH(4)
 168       RUTOL=ZTOL*DNRM2(NP1,WORK,1)
 169       RBTOL=ZTOL*DNRM2(NP1,RHO,1)
 170 C
 171 C ***** END OF INITIALIZATION *****
 172 C
 173 C ***** SOLVE SYSTEM  M Z = U *****
 174 C
 175 C COMPUTE INITIAL RESIDUAL VECTOR FOR THE SYSTEM M Z = U .
 176 C     RES = (Q**(-1))*(U - M*Z.)
 177 *
 178       CALL MULTDS(WORK(RES),AA,WORK(ZU),MAXA,NN,LENAA)
 179       WORK(RES+NN)= DDOT(NN,YP,1,WORK(ZU),1)
 180       CALL DAXPY(NP1,WORK(ZU+NN),YP,1,WORK(RES),1)
 181       CALL DSCAL(NP1,-1.0D0,WORK(RES),1)
 182       CALL DAXPY(NN,-1.0D0,PP,1,WORK(RES),1)
 183       CALL DAXPY(NN,-1.0D0,YP,1,WORK(RES),1)
 184       CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK(RES))
 185 C
 186 C COMPUTE INITIAL DIRECTION VECTOR.
 187 C     DIR = (A**T)*(Q**(-T))*RES.
 188 C
 189       CALL DCOPY(NP1,WORK(RES),1,WORK,1)
 190       CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
 191       CALL MULTDS(WORK(DIR),AA,WORK,MAXA,NN,LENAA)
 192       WORK(DIR+NN)=DDOT(NN,YP,1,WORK,1)
 193       CALL DAXPY(NP1,WORK(NP1),YP,1,WORK(DIR),1)
 194 C
 195 C COMPUTE INITIAL INNER PRODUCTS.
 196 C
 197       RUNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
 198       PUNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
 199 C
 200 C REPEAT UNTIL CONVERGENCE OR TOO MANY ITERATIONS.
 201 C
 202       J=1
 203 C
 204 C     DO WHILE ((STILLU) .AND. (J .LE. IMAX))
 205 100   IF (.NOT. ((STILLU) .AND. (J .LE. IMAX)) ) GO TO 200
 206 C
 207 C        IF ||RESIDUAL|| IS STILL NOT SMALL ENOUGH, CONTINUE.
 208 C
 209          IF (SQRT(RUNPRD) .GT. RUTOL) THEN
 210 C
 211 C           IF DIRECTION VECTOR IS ZERO, THEN RE-COMPUTE RESIDUAL,
 212 C           DIRECTION VECTOR, AND INNER PRODUCTS FROM SCRATCH
 213 C           (RATHER THAN FROM UPDATES OF PREVIOUS VALUES).
 214 C
 215             IF (PUNPRD .EQ. 0.0) THEN
 216 C
 217 C              COMPUTE RESIDUAL.
 218 C
 219                CALL MULTDS(WORK(RES),AA,WORK(ZU),MAXA,NN,LENAA)
 220                WORK(RES+NN)= DDOT(NN,YP,1,WORK(ZU),1)
 221                CALL DAXPY(NP1,WORK(ZU+NN),YP,1,WORK(RES),1)
 222                CALL DSCAL(NP1,-1.0D0,WORK(RES),1)
 223                CALL DAXPY(NN,-1.0D0,PP,1,WORK(RES),1)
 224                CALL DAXPY(NN,-1.0D0,YP,1,WORK(RES),1)
 225                CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK(RES))
 226 C
 227 C              COMPUTE DIRECTION VECTOR.
 228 C
 229                CALL DCOPY(NP1,WORK(RES),1,WORK,1)
 230                CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
 231                CALL MULTDS(WORK(DIR),AA,WORK,MAXA,NN,LENAA)
 232                WORK(DIR+NN)=DDOT(NN,YP,1,WORK,1)
 233                CALL DAXPY(NP1,WORK(NP1),YP,1,WORK(DIR),1)
 234 C
 235 C              COMPUTE INNER PRODUCTS
 236 C
 237                RUNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
 238                PUNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
 239 C
 240 C              CHECK FOR CONVERGENCE.
 241 C
 242                IF (SQRT(RUNPRD) .LE. RUTOL) THEN
 243                   STILLU=.FALSE.
 244                ENDIF
 245             ENDIF
 246             IF (STILLU) THEN
 247 C
 248 C              UPDATE SOLUTION VECTOR.
 249 C                Z = Z + AU*DIR, WHERE AU= RUNPRD/PUNPRD.
 250 C
 251                AU=RUNPRD/PUNPRD
 252                CALL DAXPY(NP1,AU,WORK(DIR),1,WORK(ZU),1)
 253 C
 254 C              COMPUTE RELATIVE CHANGE IN THE SOLUTION.
 255 C
 256                DZNRM=AU*SQRT(PUNPRD)
 257                ZLEN=DNRM2(NP1,WORK(ZU),1)
 258 C
 259 C              IF RELATIVE CHANGE IN SOLUTIONS IS SMALL ENOUGH, EXIT.
 260 C
 261                IF ( (DZNRM/ZLEN) .LT. ZTOL) STILLU=.FALSE.
 262             ENDIF
 263          ELSE
 264             STILLU=.FALSE.
 265          ENDIF
 266 C
 267 C        IF NO EXIT CRITERIA FOR  MZ=U  HAVE BEEN MET, UPDATE RESIDUAL,
 268 C          DIRECTION VECTORS, AND INNER PRODUCTS FOR NEXT ITERATION.
 269 C
 270          IF (STILLU) THEN
 271 C
 272 C           UPDATE RESIDUAL VECTOR; COMPUTE <RES,RES>.
 273 C             RES = RES - AU*(Q**(-1))*M*DIR.
 274 C
 275             CALL MULTDS(WORK,AA,WORK(DIR),MAXA,NN,LENAA)
 276             WORK(NP1)=DDOT(NN,YP,1,WORK(DIR),1)
 277             CALL DAXPY(NP1,WORK(DIR+NN),YP,1,WORK,1)
 278             CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
 279             CALL DAXPY(NP1,-AU,WORK,1,WORK(RES),1)
 280             RNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
 281 C
 282 C           UPDATE DIRECTION VECTOR; COMPUTE <DIR,DIR>.
 283 C             DIR = (M**T)*(Q**(-T))*RES + BU*DIR,
 284 C             WHERE  BU = RNPRD/RUNPRD. (NOTE: START IS USED AS
 285 C             A WORK ARRAY HERE).
 286 C
 287             BU=RNPRD/RUNPRD
 288             RUNPRD=RNPRD
 289             CALL DCOPY(NP1,WORK(RES),1,WORK,1)
 290             CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
 291             CALL MULTDS(START,AA,WORK,MAXA,NN,LENAA)
 292             START(NP1)=DDOT(NN,YP,1,WORK,1)
 293             CALL DAXPY(NP1,WORK(NP1),YP,1,START,1)
 294             CALL DAXPY(NP1,BU,WORK(DIR),1,START,1)
 295             CALL DCOPY(NP1,START,1,WORK(DIR),1)
 296             PUNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
 297          ENDIF
 298 C
 299          J=J+1
 300       GO TO 100
 301 200   CONTINUE
 302 C     END DO
 303 C
 304 C SET ERROR FLAG IF THE CONJUGATE GRADIENT ITERATION DID NOT CONVERGE.
 305 C
 306       IF (J .GT. IMAX) THEN
 307          IFLAG=4
 308          RETURN
 309       ENDIF
 310 C
 311 C ***** END OF M Z = U SYSTEM *****
 312 C
 313 C ***** SOLVE SYSTEM M Z = B *****
 314 C
 315 C
 316 C COMPUTE INITIAL RESIDUAL VECTOR FOR THE SYSTEM M Z = B .
 317 C
 318       CALL MULTDS(WORK(RES),AA,WORK(ZB),MAXA,NN,LENAA)
 319       WORK(RES+NN)=DDOT(NN,YP,1,WORK(ZB),1)
 320       CALL DAXPY(NP1,WORK(ZB+NN),YP,1,WORK(RES),1)
 321       CALL DSCAL(NP1,-1.0D0,WORK(RES),1)
 322       CALL DAXPY(NP1,-1.0D0,RHO,1,WORK(RES),1)
 323       CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK(RES))
 324 C
 325 C COMPUTE INITIAL DIRECTION VECTOR.
 326 C
 327       CALL DCOPY(NP1,WORK(RES),1,WORK,1)
 328       CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
 329       CALL MULTDS(WORK(DIR),AA,WORK,MAXA,NN,LENAA)
 330       WORK(DIR+NN)=DDOT(NN,YP,1,WORK,1)
 331       CALL DAXPY(NP1,WORK(NP1),YP,1,WORK(DIR),1)
 332 C
 333 C COMPUTE INITIAL INNER PRODUCTS.
 334 C
 335       RBNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
 336       PBNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
 337 C
 338 C REPEAT UNTIL CONVERGENCE, OR TOO MANY ITERATIONS.
 339 C
 340       J=1
 341 C
 342 C     DO WHILE ( STILLB  .AND. (J .LE. IMAX) )
 343 300   IF (.NOT. ( STILLB  .AND. (J .LE. IMAX) ) ) GO TO 400
 344 C
 345 C        IF ||RESIDUAL|| IS STILL NOT SMALL ENOUGH, CONTINUE.
 346 C
 347          IF (SQRT(RBNPRD) .GT. RBTOL) THEN
 348 C
 349 C           IF DIRECTION VECTOR IS ZERO, RE-COMPUTE RESIDUAL,
 350 C             DIRECTION VECTOR, AND INNER PRODUCTS FROM SCRATCH.
 351 C
 352             IF (PBNPRD .EQ. 0.0) THEN
 353 C
 354 C              COMPUTE RESIDUAL.
 355 C
 356                CALL MULTDS(WORK(RES),AA,WORK(ZB),MAXA,NN,LENAA)
 357                WORK(RES+NN)=DDOT(NN,YP,1,WORK(ZB),1)
 358                CALL DAXPY(NP1,WORK(ZB+NN),YP,1,WORK(RES),1)
 359                CALL DSCAL(NP1,-1.0D0,WORK(RES),1)
 360                CALL DAXPY(NP1,-1.0D0,RHO,1,WORK(RES),1)
 361                CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK(RES))
 362 C
 363 C              COMPUTE DIRECTION VECTOR.
 364 C
 365                CALL DCOPY(NP1,WORK(RES),1,WORK,1)
 366                CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
 367                CALL MULTDS(WORK(DIR),AA,WORK,MAXA,NN,LENAA)
 368                WORK(DIR+NN)=DDOT(NN,YP,1,WORK,1)
 369                CALL DAXPY(NP1,WORK(NP1),YP,1,WORK(DIR),1)
 370 C
 371 C              COMPUTE INNER PRODUCTS.
 372 C
 373                RBNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
 374                PBNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
 375 C
 376 C              CHECK FOR CONVERGENCE.
 377 C
 378                IF (SQRT(RBNPRD) .LE. RBTOL) THEN
 379                   STILLB=.FALSE.
 380                ENDIF
 381             ENDIF
 382             IF (STILLB) THEN
 383 C
 384 C              UPDATE SOLUTION VECTOR.
 385 C                Z = Z + AB*DIR, WHERE AB=RBNPRD/PBNPRD.
 386 C
 387                AB=RBNPRD/PBNPRD
 388                CALL DAXPY(NP1,AB,WORK(DIR),1,WORK(ZB),1)
 389 C
 390 C              COMPUTE RELATIVE CHANGE IN SOLUTIONS.
 391 C
 392                DZNRM=AB*SQRT(PBNPRD)
 393                ZLEN=DNRM2(NP1,WORK(ZB),1)
 394 C
 395 C              IF RELATIVE CHANGE IN SOLUTIONS IS SMALL ENOUGH, EXIT.
 396 C
 397                IF ( (DZNRM/ZLEN) .LT. ZTOL) STILLB=.FALSE.
 398             ENDIF
 399          ELSE
 400             STILLB=.FALSE.
 401          ENDIF
 402 C
 403 C        IF NO EXIT CRITERIA FOR  MZ=B  HAVE BEEN MET, UPDATE RESIDUAL,
 404 C          DIRECTION VECTORS, AND INNER PRODUCTS.
 405 C
 406          IF (STILLB) THEN
 407 C
 408 C           UPDATE RESIDUAL VECTOR; COMPUTE <RES,RES>.
 409 C             RES = RES - AB*(Q**(-1))*M*DIR.
 410 C
 411             CALL MULTDS(WORK,AA,WORK(DIR),MAXA,NN,LENAA)
 412             WORK(NP1)=DDOT(NN,YP,1,WORK(DIR),1)
 413             CALL DAXPY(NP1,WORK(DIR+NN),YP,1,WORK,1)
 414             CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
 415             CALL DAXPY(NP1,-AB,WORK,1,WORK(RES),1)
 416             RNPRD=DDOT(NP1,WORK(RES),1,WORK(RES),1)
 417 C
 418 C           UPDATE DIRECTION VECTOR; COMPUTE <DIR,DIR>.
 419 C             DIR = (M**T)*(Q**(-T))*RES + BB*DIR,
 420 C             WHERE BB=RNPRD/RBNPRD.
 421 C           (NOTE:  START IS USED AS A WORK ARRAY HERE).
 422 C
 423             BB=RNPRD/RBNPRD
 424             RBNPRD=RNPRD
 425             CALL DCOPY(NP1,WORK(RES),1,WORK,1)
 426             CALL SOLVDS(NP1,WORK(Q),LENQ,MAXA,WORK)
 427             CALL MULTDS(START,AA,WORK,MAXA,NN,LENAA)
 428             START(NP1)=DDOT(NN,YP,1,WORK,1)
 429             CALL DAXPY(NP1,WORK(NP1),YP,1,START,1)
 430             CALL DAXPY(NP1,BB,WORK(DIR),1,START,1)
 431             CALL DCOPY(NP1,START,1,WORK(DIR),1)
 432             PBNPRD=DDOT(NP1,WORK(DIR),1,WORK(DIR),1)
 433          ENDIF
 434 C
 435          J=J+1
 436       GO TO 300
 437 400   CONTINUE
 438 C     END DO
 439 C
 440 C SET ERROR FLAG IF THE CONJUGATE GRADIENT ITERATION DID NOT CONVERGE.
 441 C
 442       IF (J .GT. IMAX) THEN
 443          IFLAG=4
 444          RETURN
 445       ENDIF
 446 C
 447 C ***** END OF  M Z = B  SYSTEM *****
 448 C
 449 C COMPUTE FINAL SOLUTION VECTOR  X,  AND RETURN IT IN START.
 450 C
 451       TEMP=-WORK(ZB+NN)/(1.0D0+WORK(ZU+NN))
 452       CALL DCOPY(NP1,WORK(ZB),1,START,1)
 453       CALL DAXPY(NP1,TEMP,WORK(ZU),1,START,1)
 454 C
 455       RETURN
 456       END