luksan/mssubs.for

   1 * SUBROUTINE MXBSBM                ALL SYSTEMS                92/12/01
   2 * PURPOSE :
   3 * MULTIPLICATION OF A BLOCKED SYMMETRIC MATRIX A BY A VECTOR X.
   4 *
   5 * PARAMETERS :
   6 * PARAMETERS :
   7 *  II  L  BLOCK DIMENSION.
   8 *  RI  ABL(L*(L+1)/2)  VALUES OF NONZERO ELEMENTS OF THE GIVEN BLOCK.
   9 *  II  JBL(L)  INDICES OF THE INDIVIDUAL BLOCKS
  10 *  RI  X(N)  UNPACKED INPUT VECTOR.
  11 *  RI  Y(N)  UNPACKED OR PACKED OUTPUT VECTOR EQUAL TO A*X.
  12 *  II  JOB  FORM OF THE VECTOR Y. JOB=1-UNPACKED FORM. JOB=2-PACKED
  13 *         FORM.
  14 *
  15       SUBROUTINE MXBSBM(L,ABL,JBL,X,Y,JOB)
  16       INTEGER L,JBL(*),JOB
  17       DOUBLE PRECISION ABL(*),X(*),Y(*)
  18       INTEGER I,J,IP,JP,K
  19       DOUBLE PRECISION TEMP
  20       DOUBLE PRECISION ZERO
  21       PARAMETER  (ZERO = 0.0D 0)
  22       DO 1 I=1,L
  23       IP=JBL(I)
  24       IF (IP.GT.0) THEN
  25       IF (JOB.EQ.1) THEN
  26       Y(IP)=ZERO
  27       ELSE
  28       Y(I)=ZERO
  29       END IF
  30       END IF
  31    1  CONTINUE
  32       K=0
  33       DO 4 I=1,L
  34       IP=JBL(I)
  35       IF (IP.GT.0) TEMP=X(IP)
  36       IF (JOB.EQ.1) THEN
  37       DO 2 J=1,I-1
  38       JP=JBL(J)
  39       K=K+1
  40       IF (IP.GT.0.AND.JP.GT.0) THEN
  41       Y(IP)=Y(IP)+ABL(K)*X(JP)
  42       Y(JP)=Y(JP)+ABL(K)*TEMP
  43       END IF
  44     2 CONTINUE
  45       K=K+1
  46       IF (IP.GT.0) Y(IP)=Y(IP)+ABL(K)*TEMP
  47       ELSE
  48       DO 3 J=1,I-1
  49       JP=JBL(J)
  50       K=K+1
  51       IF (IP.GT.0.AND.JP.GT.0) THEN
  52       Y(I)=Y(I)+ABL(K)*X(JP)
  53       Y(J)=Y(J)+ABL(K)*TEMP
  54       END IF
  55     3 CONTINUE
  56       K=K+1
  57       IF (IP.GT.0) Y(I)=Y(I)+ABL(K)*TEMP
  58       END IF
  59     4 CONTINUE
  60       RETURN
  61       END
  62 * SUBROUTINE MXBSBU                ALL SYSTEMS                92/12/01
  63 * PURPOSE :
  64 * CORRECTION OF A BLOCKED SYMMETRIC MATRIX A. THE CORRECTION IS DEFINED
  65 * AS A:=A+ALF*X*TRANS(X) WHERE ALF IS A GIVEN SCALING FACTOR AND X IS
  66 * A GIVEN VECTOR.
  67 *
  68 * PARAMETERS :
  69 *  II  L  BLOCK DIMENSION.
  70 *  RI  ABL(L*(L+1)/2)  VALUES OF NONZERO ELEMENTS OF THE GIVEN BLOCK.
  71 *  II  JBL(L)  INDICES OF THE INDIVIDUAL BLOCKS
  72 *  RI  ALF  SCALING FACTOR.
  73 *  RI  X(N)  UNPACKED OR PACKED INPUT VECTOR.
  74 *  II  JOB  FORM OF THE VECTOR X. JOB=1-UNPACKED FORM. JOB=2-PACKED
  75 *         FORM.
  76 *
  77       SUBROUTINE MXBSBU(L,ABL,JBL,ALF,X,JOB)
  78       INTEGER L,JBL(*),JOB
  79       DOUBLE PRECISION ABL(*),ALF,X(*)
  80       INTEGER I,J,IP,JP,K
  81       K=0
  82       IF (JOB.EQ.1) THEN
  83       DO 3 I=1,L
  84         IP=JBL(I)
  85         DO 2 J=1,I
  86           JP=JBL(J)
  87           K=K+1
  88           IF (IP.GT.0.AND.JP.GT.0) THEN
  89             ABL(K)=ABL(K)+ALF*X(IP)*X(JP)
  90           END IF
  91     2   CONTINUE
  92     3 CONTINUE
  93       ELSE
  94       DO 5 I=1,L
  95         IP=JBL(I)
  96         DO 4 J=1,I
  97           JP=JBL(J)
  98           K=K+1
  99           IF (IP.GT.0.AND.JP.GT.0) THEN
 100           ABL(K)=ABL(K)+ALF*X(I)*X(J)
 101           END IF
 102     4   CONTINUE
 103     5 CONTINUE
 104       END IF
 105       RETURN
 106       END
 107 * SUBROUTINE MXBSMI                ALL SYSTEMS                91/12/01
 108 * PURPOSE :
 109 * BLOCKS OF THE SYMMETRIC BLOCKED MATRIX ARE SET TO THE UNIT MATRICES.
 110 *
 111 * PARAMETERS :
 112 *  II  NBLKS  NUMBER OF BLOCKS OF THE MATRIX.
 113 *  RI  ABL(NBLA)  VALUES OF THE NONZERO ELEMENTS OF THE MATRIX.
 114 *  II  IBLBG(NBLKS+1)  BEGINNINGS OF THE BLOCKS IN THE MATRIX.
 115 *
 116 * SUBROUTINES USED :
 117 *  MXDSMI  DENSE SYMMETRIC MATRIX IS SET TO THE UNIT MATRIX.
 118 *
 119       SUBROUTINE MXBSMI(NBLKS,ABL,IBLBG)
 120       INTEGER NBLKS,IBLBG(*)
 121       DOUBLE PRECISION ABL(*)
 122       INTEGER I,K,KBEG,KLEN
 123       K=1
 124       DO 1 I=1,NBLKS
 125       KBEG=IBLBG(I)
 126       KLEN=IBLBG(I+1)-KBEG
 127       CALL MXDSMI(KLEN,ABL(K))
 128       K=K+KLEN*(KLEN+1)/2
 129    1  CONTINUE
 130       RETURN
 131       END
 132 * SUBROUTINE MXDCMD               ALL SYSTEMS                91/12/01
 133 * PURPOSE :
 134 * MULTIPLICATION OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A
 135 * BY A VECTOR X AND ADDITION OF THE SCALED VECTOR ALF*Y.
 136 *
 137 * PARAMETERS :
 138 *  II  N  NUMBER OF ROWS OF THE MATRIX A.
 139 *  II  M  NUMBER OF COLUMNS OF THE MATRIX A.
 140 *  RI  A(N*M)  RECTANGULAR MATRIX STORED COLUMNWISE IN THE
 141 *         ONE-DIMENSIONAL ARRAY.
 142 *  RI  X(M)  INPUT VECTOR.
 143 *  RI  ALF  SCALING FACTOR.
 144 *  RI  Y(N)  INPUT VECTOR.
 145 *  RO  Z(N)  OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
 146 *
 147 * SUBPROGRAMS USED :
 148 *  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 149 *  S   MXVSCL  SCALING OF A VECTOR.
 150 *
 151       SUBROUTINE MXDCMD(N,M,A,X,ALF,Y,Z)
 152       INTEGER N,M
 153       DOUBLE PRECISION A(*),X(*),ALF,Y(*),Z(*)
 154       INTEGER J,K
 155       CALL MXVSCL(N,ALF,Y,Z)
 156       K=0
 157       DO 1 J=1,M
 158       CALL MXVDIR(N,X(J),A(K+1),Z,Z)
 159       K=K+N
 160     1 CONTINUE
 161       RETURN
 162       END
 163 * SUBROUTINE MXDCMU               ALL SYSTEMS                91/12/01
 164 * PURPOSE :
 165 * UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A. THIS MATRIX
 166 * IS UPDATED BY THE RULE A:=A+ALF*X*TRANS(Y).
 167 *
 168 * PARAMETERS :
 169 *  II  N  NUMBER OF ROWS OF THE MATRIX A.
 170 *  II  M  NUMBER OF COLUMNS OF THE MATRIX A.
 171 *  RU  A(N*M)  RECTANGULAR MATRIX STORED COLUMNWISE IN THE
 172 *         ONE-DIMENSIONAL ARRAY.
 173 *  RI  ALF  SCALAR PARAMETER.
 174 *  RI  X(N)  INPUT VECTOR.
 175 *  RI  Y(M)  INPUT VECTOR.
 176 *
 177       SUBROUTINE MXDCMU(N,M,A,ALF,X,Y)
 178       INTEGER N,M
 179       DOUBLE PRECISION A(*),ALF,X(*),Y(*)
 180       DOUBLE PRECISION TEMP
 181       INTEGER I,J,K
 182       K=0
 183       DO 2 J=1,M
 184       TEMP=ALF*Y(J)
 185       DO 1 I=1,N
 186       A(K+I)=A(K+I)+TEMP*X(I)
 187     1 CONTINUE
 188       K=K+N
 189     2 CONTINUE
 190       RETURN
 191       END
 192 * SUBROUTINE MXDCMV               ALL SYSTEMS                91/12/01
 193 * PURPOSE :
 194 * RANK-TWO UPDATE OF A COLUMNWISE STORED DENSE RECTANGULAR MATRIX A.
 195 * THIS MATRIX IS UPDATED BY THE RULE A:=A+ALF*X*TRANS(U)+BET*Y*TRANS(V).
 196 *
 197 * PARAMETERS :
 198 *  II  N  NUMBER OF ROWS OF THE MATRIX A.
 199 *  II  M  NUMBER OF COLUMNS OF THE MATRIX A.
 200 *  RU  A(N*M)  RECTANGULAR MATRIX STORED COLUMNWISE IN THE
 201 *         ONE-DIMENSIONAL ARRAY.
 202 *  RI  ALF  SCALAR PARAMETER.
 203 *  RI  X(N)  INPUT VECTOR.
 204 *  RI  U(M)  INPUT VECTOR.
 205 *  RI  BET  SCALAR PARAMETER.
 206 *  RI  Y(N)  INPUT VECTOR.
 207 *  RI  V(M)  INPUT VECTOR.
 208 *
 209       SUBROUTINE MXDCMV(N,M,A,ALF,X,U,BET,Y,V)
 210       INTEGER N,M
 211       DOUBLE PRECISION A(*),ALF,X(*),U(*),BET,Y(*),V(*)
 212       DOUBLE PRECISION TEMPA,TEMPB
 213       INTEGER I,J,K
 214       K=0
 215       DO 2 J=1,M
 216       TEMPA=ALF*U(J)
 217       TEMPB=BET*V(J)
 218       DO 1 I=1,N
 219       A(K+I)=A(K+I)+TEMPA*X(I)+TEMPB*Y(I)
 220     1 CONTINUE
 221       K=K+N
 222     2 CONTINUE
 223       RETURN
 224       END
 225 * SUBROUTINE MXDPGB                ALL SYSTEMS                91/12/01
 226 * PURPOSE :
 227 * SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A DENSE SYMMETRIC
 228 * POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
 229 * OBTAINED BY THE SUBROUTINE MXDPGF.
 230 *
 231 * PARAMETERS :
 232 *  II  N ORDER OF THE MATRIX A.
 233 *  RI  A(N*(N+1)/2) FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
 234 *         SUBROUTINE MXDPGF.
 235 *  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
 236 *         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
 237 *         EQUATIONS.
 238 *  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN
 239 *         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
 240 *
 241 * METHOD :
 242 * BACK SUBSTITUTION
 243 *
 244       SUBROUTINE MXDPGB(N,A,X,JOB)
 245       INTEGER          JOB,N
 246       DOUBLE PRECISION A(*),X(*)
 247       INTEGER          I,II,IJ,J
 248       IF (JOB.GE.0) THEN
 249 *
 250 *     PHASE 1 : X:=L**(-1)*X
 251 *
 252           IJ = 0
 253           DO 20 I = 1,N
 254               DO 10 J = 1,I - 1
 255                   IJ = IJ + 1
 256                   X(I) = X(I) - A(IJ)*X(J)
 257    10         CONTINUE
 258               IJ = IJ + 1
 259    20     CONTINUE
 260       END IF
 261       IF (JOB.EQ.0) THEN
 262 *
 263 *     PHASE 2 : X:=D**(-1)*X
 264 *
 265           II = 0
 266           DO 30 I = 1,N
 267               II = II + I
 268               X(I) = X(I)/A(II)
 269    30     CONTINUE
 270       END IF
 271       IF (JOB.LE.0) THEN
 272 *
 273 *     PHASE 3 : X:=TRANS(L)**(-1)*X
 274 *
 275           II = N* (N-1)/2
 276           DO 50 I = N - 1,1,-1
 277               IJ = II
 278               DO 40 J = I + 1,N
 279                   IJ = IJ + J - 1
 280                   X(I) = X(I) - A(IJ)*X(J)
 281    40         CONTINUE
 282               II = II - I
 283    50     CONTINUE
 284       END IF
 285       RETURN
 286       END
 287 * SUBROUTINE MXDPGF                ALL SYSTEMS                89/12/01
 288 * PURPOSE :
 289 * FACTORIZATION A+E=L*D*TRANS(L) OF A DENSE SYMMETRIC POSITIVE DEFINITE
 290 * MATRIX A+E WHERE D AND E ARE DIAGONAL POSITIVE DEFINITE MATRICES AND
 291 * L IS A LOWER TRIANGULAR MATRIX. IF A IS SUFFICIENTLY POSITIVE
 292 * DEFINITE THEN E=0.
 293 *
 294 * PARAMETERS :
 295 *  II  N ORDER OF THE MATRIX A.
 296 *  RU  A(N*(N+1)/2)  ON INPUT A GIVEN DENSE SYMMETRIC (USUALLY POSITIVE
 297 *         DEFINITE) MATRIX A STORED IN THE PACKED FORM. ON OUTPUT THE
 298 *         COMPUTED FACTORIZATION A+E=L*D*TRANS(L).
 299 *  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
 300 *         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
 301 *         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
 302 *         INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
 303 *         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
 304 *         PROCESS.
 305 *  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON
 306 *         OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
 307 *         FACTORIZATION PROCESS (IF INF>0).
 308 *  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
 309 *
 310 * METHOD :
 311 * P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND
 312 * LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974)
 313 * PP. 311-350.
 314 *
 315       SUBROUTINE MXDPGF(N,A,INF,ALF,TAU)
 316       DOUBLE PRECISION ALF,TAU
 317       INTEGER          INF,N
 318       DOUBLE PRECISION A(*)
 319       DOUBLE PRECISION BET,DEL,GAM,RHO,SIG,TOL
 320       INTEGER          I,IJ,IK,J,K,KJ,KK,L
 321       L = 0
 322       INF = 0
 323       TOL = ALF
 324 *
 325 *     ESTIMATION OF THE MATRIX NORM
 326 *
 327       ALF = 0.0D0
 328       BET = 0.0D0
 329       GAM = 0.0D0
 330       TAU = 0.0D0
 331       KK = 0
 332       DO 20 K = 1,N
 333           KK = KK + K
 334           BET = MAX(BET,ABS(A(KK)))
 335           KJ = KK
 336           DO 10 J = K + 1,N
 337               KJ = KJ + J - 1
 338               GAM = MAX(GAM,ABS(A(KJ)))
 339    10     CONTINUE
 340    20 CONTINUE
 341       BET = MAX(TOL,BET,GAM/N)
 342 *      DEL = TOL*BET
 343       DEL = TOL*MAX(BET,1.0D0)
 344       KK = 0
 345       DO 60 K = 1,N
 346           KK = KK + K
 347 *
 348 *     DETERMINATION OF A DIAGONAL CORRECTION
 349 *
 350           SIG = A(KK)
 351           IF (ALF.GT.SIG) THEN
 352               ALF = SIG
 353               L = K
 354           END IF
 355           GAM = 0.0D0
 356           KJ = KK
 357           DO 30 J = K + 1,N
 358               KJ = KJ + J - 1
 359               GAM = MAX(GAM,ABS(A(KJ)))
 360    30     CONTINUE
 361           GAM = GAM*GAM
 362           RHO = MAX(ABS(SIG),GAM/BET,DEL)
 363           IF (TAU.LT.RHO-SIG) THEN
 364               TAU = RHO - SIG
 365               INF = -1
 366           END IF
 367 *
 368 *     GAUSSIAN ELIMINATION
 369 *
 370           A(KK) = RHO
 371           KJ = KK
 372           DO 50 J = K + 1,N
 373               KJ = KJ + J - 1
 374               GAM = A(KJ)
 375               A(KJ) = GAM/RHO
 376               IK = KK
 377               IJ = KJ
 378               DO 40 I = K + 1,J
 379                   IK = IK + I - 1
 380                   IJ = IJ + 1
 381                   A(IJ) = A(IJ) - A(IK)*GAM
 382    40         CONTINUE
 383    50     CONTINUE
 384    60 CONTINUE
 385       IF (L.GT.0 .AND. ABS(ALF).GT.DEL) INF = L
 386       RETURN
 387       END
 388 * SUBROUTINE MXDRMM               ALL SYSTEMS                91/12/01
 389 * PURPOSE :
 390 * MULTIPLICATION OF A ROWWISE STORED DENSE RECTANGULAR MATRIX A BY
 391 * A VECTOR X.
 392 *
 393 * PARAMETERS :
 394 *  II  N  NUMBER OF COLUMNS OF THE MATRIX A.
 395 *  II  M  NUMBER OF ROWS OF THE MATRIX A.
 396 *  RI  A(M*N)  RECTANGULAR MATRIX STORED ROWWISE IN THE
 397 *         ONE-DIMENSIONAL ARRAY.
 398 *  RI  X(N)  INPUT VECTOR.
 399 *  RO  Y(M)  OUTPUT VECTOR EQUAL TO A*X.
 400 *
 401       SUBROUTINE MXDRMM(N,M,A,X,Y)
 402       INTEGER N,M
 403       DOUBLE PRECISION A(*),X(*),Y(*)
 404       DOUBLE PRECISION TEMP
 405       INTEGER I,J,K
 406       DOUBLE PRECISION ZERO
 407       PARAMETER (ZERO=0.0D 0)
 408       K=0
 409       DO 2 J=1,M
 410       TEMP=ZERO
 411       DO 1 I=1,N
 412       TEMP=TEMP+A(K+I)*X(I)
 413     1 CONTINUE
 414       Y(J)=TEMP
 415       K=K+N
 416     2 CONTINUE
 417       RETURN
 418       END
 419 * SUBROUTINE MXDRCB               ALL SYSTEMS                91/12/01
 420 * PURPOSE :
 421 * BACKWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION OF
 422 * THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
 423 *
 424 * PARAMETERS :
 425 *  II  N  NUMBER OF ROWS OF THE MATRICES A AND B.
 426 *  II  M  NUMBER OF COLUMNS OF THE MATRICES A AND B.
 427 *  RI  A(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 428 *  RI  B(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 429 *  RI  U(M)  VECTOR OF SCALAR COEFFICIENTS.
 430 *  RO  V(M)  VECTOR OF SCALAR COEFFICIENTS.
 431 *  RU  X(N)  PREMULTIPLIED VECTOR.
 432 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 433 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 434 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 435 *         IX(I).EQ.-5.
 436 *
 437 * SUBPROGRAM USED :
 438 *  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 439 *  RF  MXUDOT  DOT PRODUCT OF VECTORS.
 440 *
 441 * METHOD :
 442 * H.MATTHIES, G.STRANG: THE SOLUTION OF NONLINEAR FINITE ELEMENT
 443 * EQUATIONS. INT.J.NUMER. METHODS ENGN. 14 (1979) 1613-1626.
 444 *
 445       SUBROUTINE MXDRCB(N,M,A,B,U,V,X,IX,JOB)
 446       INTEGER N,M,IX(*),JOB
 447       DOUBLE PRECISION A(*),B(*),U(*),V(*),X(*)
 448       DOUBLE PRECISION MXUDOT
 449       INTEGER I,K
 450       K=1
 451       DO 1 I=1,M
 452       V(I)=U(I)*MXUDOT(N,X,A(K),IX,JOB)
 453       CALL MXUDIR(N,-V(I),B(K),X,X,IX,JOB)
 454       K=K+N
 455     1 CONTINUE
 456       RETURN
 457       END
 458 * SUBROUTINE MXDRCF               ALL SYSTEMS                91/12/01
 459 * PURPOSE :
 460 * FORWARD PART OF THE STRANG FORMULA FOR PREMULTIPLICATION OF
 461 * THE VECTOR X BY AN IMPLICIT BFGS UPDATE.
 462 *
 463 * PARAMETERS :
 464 *  II  N  NUMBER OF ROWS OF THE MATRICES A AND B.
 465 *  II  M  NUMBER OF COLUMNS OF THE MATRICES A AND B.
 466 *  RI  A(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 467 *  RI  B(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 468 *  RI  U(M)  VECTOR OF SCALAR COEFFICIENTS.
 469 *  RI  V(M)  VECTOR OF SCALAR COEFFICIENTS.
 470 *  RU  X(N)  PREMULTIPLIED VECTOR.
 471 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
 472 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
 473 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
 474 *         IX(I).EQ.-5.
 475 *
 476 * SUBPROGRAM USED :
 477 *  S   MXUDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 478 *  RF  MXUDOT  DOT PRODUCT OF VECTORS.
 479 *
 480 * METHOD :
 481 * H.MATTHIES, G.STRANG: THE SOLUTION OF NONLINEAR FINITE ELEMENT
 482 * EQUATIONS. INT.J.NUMER. METHODS ENGN. 14 (1979) 1613-1626.
 483 *
 484       SUBROUTINE MXDRCF(N,M,A,B,U,V,X,IX,JOB)
 485       INTEGER N,M,IX(*),JOB
 486       DOUBLE PRECISION A(*),B(*),U(*),V(*),X(*)
 487       DOUBLE PRECISION TEMP,MXUDOT
 488       INTEGER I,K
 489       K=(M-1)*N+1
 490       DO 1 I=M,1,-1
 491       TEMP=U(I)*MXUDOT(N,X,B(K),IX,JOB)
 492       CALL MXUDIR(N,V(I)-TEMP,A(K),X,X,IX,JOB)
 493       K=K-N
 494     1 CONTINUE
 495       RETURN
 496       END
 497 * SUBROUTINE MXDRSU               ALL SYSTEMS                91/12/01
 498 * PURPOSE :
 499 * SHIFT OF COLUMNS OF THE RECTANGULAR MATRICES A AND B. SHIFT OF
 500 * ELEMENTS OF THE VECTOR U. THESE SHIFTS ARE USED IN THE LIMITED
 501 * MEMORY BFGS METHOD.
 502 *
 503 * PARAMETERS :
 504 *  II  N  NUMBER OF ROWS OF THE MATRIX A AND B.
 505 *  II  M  NUMBER OF COLUMNS OF THE MATRIX A AND B.
 506 *  RU  A(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 507 *  RU  B(N*M)  RECTANGULAR MATRIX STORED AS A ONE-DIMENSIONAL ARRAY.
 508 *  RU  U(M)  VECTOR.
 509 *
 510       SUBROUTINE MXDRSU(N,M,A,B,U)
 511       INTEGER N,M
 512       DOUBLE PRECISION A(*),B(*),U(*)
 513       INTEGER I,K,L
 514       K=(M-1)*N+1
 515       DO 1 I=M-1,1,-1
 516       L=K-N
 517       CALL MXVCOP(N,A(L),A(K))
 518       CALL MXVCOP(N,B(L),B(K))
 519       U(I+1)=U(I)
 520       K=L
 521     1 CONTINUE
 522       RETURN
 523       END
 524 * SUBROUTINE MXDSMI                ALL SYSTEMS                88/12/01
 525 * PURPOSE :
 526 * DENSE SYMMETRIC MATRIX A IS SET TO THE UNIT MATRIX WITH THE SAME
 527 * ORDER.
 528 *
 529 * PARAMETERS :
 530 *  II  N  ORDER OF THE MATRIX A.
 531 *  RO  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM
 532 *         WHICH IS SET TO THE UNIT MATRIX (I.E. A:=I).
 533 *
 534       SUBROUTINE MXDSMI(N,A)
 535       INTEGER N
 536       DOUBLE PRECISION A(*)
 537       INTEGER I, M
 538       DOUBLE PRECISION ZERO,ONE
 539       PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
 540       M = N * (N+1) / 2
 541       DO 1  I = 1, M
 542       A(I) = ZERO
 543     1 CONTINUE
 544       M = 0
 545       DO 2  I = 1, N
 546       M = M + I
 547       A(M) = ONE
 548     2 CONTINUE
 549       RETURN
 550       END
 551 * SUBROUTINE MXDSMS                ALL SYSTEMS                91/12/01
 552 * PURPOSE :
 553 * SCALING OF A DENSE SYMMETRIC MATRIX.
 554 *
 555 * PARAMETERS :
 556 *  II  N  ORDER OF THE MATRIX A.
 557 *  RU  A(N*(N+1)/2)  DENSE SYMMETRIC MATRIX STORED IN THE PACKED FORM
 558 *         WHICH IS SCALED BY THE VALUE ALF (I.E. A:=ALF*A).
 559 *  RI  ALF  SCALING FACTOR.
 560 *
 561       SUBROUTINE MXDSMS( N, A, ALF)
 562       INTEGER N
 563       DOUBLE PRECISION  A(*), ALF
 564       INTEGER I,M
 565       M = N * (N+1) / 2
 566       DO 1  I = 1, M
 567       A(I) = A(I) * ALF
 568     1 CONTINUE
 569       RETURN
 570       END
 571 * SUBROUTINE MXLIIM                ALL SYSTEMS                96/12/01
 572 * PURPOSE :
 573 * MATRIX MULTIPLICATION FOR LIMITED STORAGE INVERSE COLUMN UPDATE
 574 * METHOD.
 575 *
 576 * PARAMETERS :
 577 *  II  N  NUMBER OF VARIABLES.
 578 *  II  M  NUMBER OF QUASI-NEWTON STEPS.
 579 *  RI  D(N) DIAGONAL OF A DECOMPOSED TRIDIAGONAL MATRIX.
 580 *  RI  DL(N) SUBDIAGONAL OF A DECOMPOSED TRIDIAGONAL MATRIX.
 581 *  RI  DU(N) SUPERDIAGONAL OF A DECOMPOSED TRIDIAGONAL MATRIX.
 582 *  RI  DU2(N) SECOND SUPERDIAGONAL OF A DECOMPOSED TRIDIAGONAL MATRIX.
 583 *  II  ID(N)  PERMUTATION VECTOR.
 584 *  RI  XM(N*M)  SET OF VECTORS FOR INVERSE COLUMN UPDATE.
 585 *  RI  GM(M)  SET OF VALUES FOR INVERSE COLUMN UPDATE.
 586 *  II  IM(M)  SET OF INDICES FOR INVERSE COLUMN UPDATE.
 587 *  RI  U(N)  INPUT VECTOR.
 588 *  RO  V(N)  OUTPUT VECTOR.
 589 *
 590 * SUBPROGRAMS USED :
 591 *  S   MXSGIB  BACK SUBSTITUTION AFTER INCOMPLETE LU DECOMPOSITION.
 592 *  S   MXVCOP  COPYING OF A VECTOR.
 593 *  S   MXVDIR  VECTOR AUGMENTED BY THE SCALED VECTOR.
 594 *
 595       SUBROUTINE MXLIIM(N,M,A,IA,JA,IP,ID,XM,GM,IM,U,V,S)
 596       INTEGER          M,N
 597       DOUBLE PRECISION A(*),GM(*),S(*),U(*),V(*),XM(*)
 598       INTEGER          IA(*),ID(*),IM(*),IP(*),JA(*)
 599       INTEGER          I,L
 600       CALL MXVCOP(N,U,V)
 601       CALL MXSGIB(N,A,IA,JA,IP,ID,V,S,0)
 602       L = 1
 603       DO 10 I = 1,M
 604           CALL MXVDIR(N,U(IM(I))/GM(I),XM(L),V,V)
 605           L = L + N
 606    10 CONTINUE
 607       RETURN
 608       END
 609 * SUBROUTINE MXSCMD               ALL SYSTEMS                92/12/01
 610 * PURPOSE :
 611 * MULTIPLICATION OF A DENSE RECTANGULAR MATRIX A BY A VECTOR X AND
 612 * ADDITIOON OF THE SCALED VECTOR ALF*Y.
 613 *
 614 * PARAMETERS :
 615 *  II  N  NUMBER OF ROWS OF THE MATRIX A.
 616 *  II  NA NUMBER OF COLUMNS OF THE MATRIX A.
 617 *  II  MA  NUMBER OF ELEMENTS IN THE FIELD A.
 618 *  RI  A(MA)  RECTANGULAR MATRIX STORED AS A TWO-DIMENSIONAL ARRAY.
 619 *  II  IA(NA+1)  POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD A.
 620 *  II  JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD A.
 621 *  RI  X(NA)  INPUT VECTOR.
 622 *  RI  ALF  SCALING FACTOR.
 623 *  RI  Y(N)  INPUT VECTOR.
 624 *  RO  Z(N)  OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
 625 *
 626 * SUBPROGRAMS USED :
 627 *  S   MXVSCL  SCALING OF A VECTOR.
 628 *
 629       SUBROUTINE MXSCMD(N,NA,A,IA,JA,X,ALF,Y,Z)
 630       INTEGER N,NA,IA(*),JA(*)
 631       DOUBLE PRECISION A(*),X(*),ALF,Y(*),Z(*)
 632       INTEGER I,J,K,L,JP
 633       CALL MXVSCL(N,ALF,Y,Z)
 634       DO 2 I=1,NA
 635       K=IA(I)
 636       L=IA(I+1)-K
 637       DO 1 J=1,L
 638       JP=JA(K)
 639       IF (JP.GT.0) Z(JP)=Z(JP)+A(K)*X(I)
 640       K=K+1
 641     1 CONTINUE
 642     2 CONTINUE
 643       RETURN
 644       END
 645 * SUBROUTINE MXSCMM               ALL SYSTEMS                92/12/01
 646 * PURPOSE :
 647 * MULTIPLICATION OF A DENSE RECTANGULAR MATRIX A BY A VECTOR X.
 648 *
 649 * PARAMETERS :
 650 *  II  N  NUMBER OF ROWS OF THE MATRIX A.
 651 *  II  NA NUMBER OF COLUMNS OF THE MATRIX A.
 652 *  II  MA  NUMBER OF ELEMENTS IN THE FIELD A.
 653 *  RI  A(MA)  RECTANGULAR MATRIX STORED AS A TWO-DIMENSIONAL ARRAY.
 654 *  II  IA(NA+1)  POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD A.
 655 *  II  JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD A.
 656 *  RI  X(NA)  INPUT VECTOR.
 657 *  RO  Y(N)  OUTPUT VECTOR EQUAL TO A*X.
 658 *
 659 * SUBPROGRAMS USED :
 660 *  S   MXVSET  INITIATION OF A VECTOR.
 661 *
 662       SUBROUTINE MXSCMM(N,NA,A,IA,JA,X,Y)
 663       INTEGER N,NA,IA(*),JA(*)
 664       DOUBLE PRECISION A(*),X(*),Y(*)
 665       INTEGER I,J,K,L,JP
 666       CALL MXVSET(N,0.0D 0,Y)
 667       DO 2 I=1,NA
 668       K=IA(I)
 669       L=IA(I+1)-K
 670       DO 1 J=1,L
 671       JP=JA(K)
 672       IF (JP.GT.0) Y(JP)=Y(JP)+A(K)*X(I)
 673       K=K+1
 674     1 CONTINUE
 675     2 CONTINUE
 676       RETURN
 677       END
 678 * SUBROUTINE MXSGIB                ALL SYSTEMS                95/12/01
 679 * PURPOSE :
 680 * SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE UNSYMMETRIC
 681 * MATRIX A USING INCOMPLETE FACTORIZATION OBTAINED BY THE SUBROUTINE
 682 * MXSGIF.
 683 *
 684 * PARAMETERS :
 685 *  II  N  MATRIX DIMENSION.
 686 *  II  M  NUMBER OF MATRIX NONZERO ELEMENTS.
 687 *  RU  A(M)  NONZERO ELEMENTS OF THE MATRIX A.
 688 *  II  IA(N+1)  ROW POINTERS OF THE MATRIX A.
 689 *  II  JA(M)  COLUMN INDICES OF THE MATRIX A.
 690 *  IO  IP(N)  PERMUTATION VECTOR.
 691 *  II  ID(N) DIAGONAL POINTERS OF THE MATRIX A.
 692 *  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
 693 *         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
 694 *         EQUATIONS.
 695 *  RA  Y(N)  AUXILIARY VECTOR.
 696 *  JOB  OPTION. JOB=0 - SOLUTION WITH THE ORIGINAL MATRIX.
 697 *         JOB=1 - SOLUTION WITH THE MATRIX TRANSPOSE.
 698 *
 699       SUBROUTINE MXSGIB(N,A,IA,JA,IP,ID,X,Y,JOB)
 700       DOUBLE PRECISION CON
 701       PARAMETER        (CON=1.0D120)
 702       INTEGER          JOB,N
 703       DOUBLE PRECISION A(*),X(*),Y(*)
 704       INTEGER          IA(*),ID(*),IP(*),JA(*)
 705       DOUBLE PRECISION APOM
 706       INTEGER          J,JJ,JP,K,KSTOP,KSTRT
 707       IF (JOB.LE.0) THEN
 708           DO 20 K = 1,N
 709               KSTRT = IA(K)
 710               KSTOP = IA(K+1) - 1
 711               DO 10 JJ = KSTRT,KSTOP
 712                   J = JA(JJ)
 713                   JP = IP(J)
 714                   IF (JP.LT.K) THEN
 715                       X(K) = X(K) - A(JJ)*X(JP)
 716                       IF (ABS(X(K)).GE.CON) X(K) = SIGN(CON,X(K))
 717                   END IF
 718    10         CONTINUE
 719    20     CONTINUE
 720           DO 40 K = N,1,-1
 721               KSTRT = IA(K)
 722               KSTOP = IA(K+1) - 1
 723               DO 30 JJ = KSTRT,KSTOP
 724                   J = JA(JJ)
 725                   JP = IP(J)
 726                   IF (JP.GT.K) X(K) = X(K) - A(JJ)*X(JP)
 727                   IF (JP.EQ.K) APOM = A(JJ)
 728    30         CONTINUE
 729               X(K) = X(K)/APOM
 730    40     CONTINUE
 731           CALL MXVSFP(N,IP,X,Y)
 732       ELSE
 733           CALL MXVSBP(N,IP,X,Y)
 734           DO 60 K = 1,N
 735               X(K) = X(K)/A(ID(K))
 736               KSTRT = IA(K)
 737               KSTOP = IA(K+1) - 1
 738               DO 50 JJ = KSTRT,KSTOP
 739                   J = JA(JJ)
 740                   JP = IP(J)
 741                   IF (JP.GT.K) X(JP) = X(JP) - A(JJ)*X(K)
 742    50         CONTINUE
 743    60     CONTINUE
 744           DO 80 K = N,1,-1
 745               KSTRT = IA(K)
 746               KSTOP = IA(K+1) - 1
 747               DO 70 JJ = KSTRT,KSTOP
 748                   J = JA(JJ)
 749                   JP = IP(J)
 750                   IF (JP.LT.K) X(JP) = X(JP) - A(JJ)*X(K)
 751    70         CONTINUE
 752    80     CONTINUE
 753       END IF
 754       RETURN
 755       END
 756 * SUBROUTINE MXSGIF                ALL SYSTEMS                95/12/01
 757 * PURPOSE :
 758 * INCOMPLETE FACTORIZATION OF A GENERAL SPARSE MATRIX A.
 759 *
 760 * PARAMETERS :
 761 *  II  N  MATRIX DIMENSION.
 762 *  II  M  NUMBER OF MATRIX NONZERO ELEMENTS.
 763 *  RU  A(M)  NONZERO ELEMENTS OF THE MATRIX A.
 764 *  II  IA(N+1)  ROW POINTERS OF THE MATRIX A.
 765 *  II  JA(M)  COLUMN INDICES OF THE MATRIX A.
 766 *  IO  IP(N)  PERMUTATION VECTOR.
 767 *  IO  ID(N)  DIAGONAL POINTERS OF THE MATRIX A.
 768 *  RA  IW(N)  AUXILIARY VECTOR.
 769 *  RI  TOL  PIVOT TOLERANCE.
 770 *  IO  INF  INFORMATION.
 771 *
 772       SUBROUTINE MXSGIF(N,A,IA,JA,IP,ID,IW,TOL,INF)
 773       DOUBLE PRECISION ZERO,ONE,CON
 774       PARAMETER        (ZERO=0.0D0,ONE=1.0D0,CON=1.0D-30)
 775       DOUBLE PRECISION TOL
 776       INTEGER          INF,N
 777       DOUBLE PRECISION A(*)
 778       INTEGER          IA(*),ID(*),IP(*),IW(*),JA(*)
 779       DOUBLE PRECISION TEMP
 780       INTEGER          I,II,J,JJ,JSTOP,JSTRT,K,KK,KSTOP,KSTRT
 781       INF = 0
 782       DO 10 I = 1,N
 783           IF (IP(I).LE.0 .OR. IP(I).GT.N) THEN
 784               CALL MXVINP(N,IP)
 785               GO TO 20
 786           END IF
 787    10 CONTINUE
 788    20 CALL MXVINS(N,0,IW)
 789       DO 70 K = 1,N
 790           KSTRT = IA(K)
 791           KSTOP = IA(K+1) - 1
 792           ID(K) = 0
 793           DO 30 JJ = KSTRT,KSTOP
 794               J = JA(JJ)
 795               IW(J) = JJ
 796               IF (IP(J).EQ.K) ID(K) = JJ
 797    30     CONTINUE
 798           IF (ID(K).EQ.0) THEN
 799               INF = -45
 800               RETURN
 801           END IF
 802           IF (TOL.GT.ZERO) A(ID(K)) = (ONE+TOL)*A(ID(K))
 803           IF (ABS(A(ID(K))).LT.TOL) A(ID(K)) = A(ID(K)) +
 804      *                              SIGN(TOL,A(ID(K)))
 805           DO 50 JJ = KSTRT,KSTOP
 806               J = IP(JA(JJ))
 807               IF (J.LT.K) THEN
 808                   JSTRT = IA(J)
 809                   JSTOP = IA(J+1) - 1
 810                   TEMP = A(JJ)/A(ID(J))
 811                   A(JJ) = TEMP
 812                   DO 40 II = JSTRT,JSTOP
 813                       I = JA(II)
 814                       IF (IP(I).GT.J) THEN
 815                           KK = IW(I)
 816                           IF (KK.NE.0) A(KK) = A(KK) - TEMP*A(II)
 817                       END IF
 818    40             CONTINUE
 819               END IF
 820    50     CONTINUE
 821           KK = ID(K)
 822           IF (ABS(A(KK)).LT.CON) THEN
 823               INF = K
 824               IF (A(KK).EQ.ZERO) THEN
 825                   A(KK) = CON
 826               ELSE
 827                   A(KK) = SIGN(CON,A(KK))
 828               END IF
 829           END IF
 830           DO 60 JJ = KSTRT,KSTOP
 831               J = JA(JJ)
 832               IW(J) = 0
 833    60     CONTINUE
 834    70 CONTINUE
 835       RETURN
 836       END
 837 * SUBROUTINE MXSPCA                 ALL SYSTEMS                93/12/01
 838 * PURPOSE :
 839 * REWRITE SYMMETRIC MATRIX INTO THE PERMUTED FACTORIZED COMPACT SCHEME.
 840 * MOIDIFIED VERSION FOR THE USE WITH MXSPCJ.
 841 *
 842 * PARAMETERS:
 843 *  II  N  SIZE OF THE SYSTEM SOLVED.
 844 *  II  NB NUMBER OF NONZEROS IN THE UPPER TRIANGLE OF THE ORIGINAL
 845 *         MATRIX.
 846 *  II  ML SIZE OF THE COMPACT FACTOR.
 847 *  RU  A(MMAX) NUMERICAL VALUES OF THE SPARSE HESSIAN APPROXIMATION
 848 *              STORED AT THE POSITIONS 1, ...,NB.
 849 *  IU  JA(MMAX) INDICES OF THE NONZERO ELEMENTS OF THE HESSIAN MATRIX IN
 850 *             THE PACKED ROW FORM AT THE FIRST NB POSITIONS.
 851 *             NEW POSITIONS
 852 *             IN THE PERMUTED FACTOR STORED IN A(NB+1), ..., A(2*NB),
 853 *             INDICES OF COMPACT SCHEME IN A(2*NB+1), ..., A(2*NB+ML).
 854 *  II  PSL(N+1)  POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR
 855 *             FACTOR OF THE HESSIAN APPROXIMATION.
 856 *  RI  T  CORRECTION FACTOR THAT IS ADDED TO THE DIAGONAL.
 857 *
 858 *
 859       SUBROUTINE MXSPCA(N,NB,ML,A,IA,JA,T)
 860       INTEGER N,NB,ML,IA(*),JA(*)
 861       DOUBLE PRECISION A(*),T
 862       INTEGER I,J
 863       DO 100 I=1,N
 864         J=ABS(JA(IA(I)+NB+ML))
 865         A(NB+J)=A(NB+J)+T
 866  100  CONTINUE
 867       RETURN
 868       END
 869 * SUBROUTINE MXSPCB                ALL SYSTEMS                92/12/01
 870 * PURPOSE :
 871 * SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE SYMMETRIC
 872 * POSITIVE DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L)
 873 * STORED IN THE COMPACT SCHEME. THIS FACTORIZATION CAN BE OBTAINED
 874 * USING THE SUBROUTINE MXSPCF.
 875 *
 876 * PARAMETERS :
 877 *  II  N ORDER OF THE MATRIX A.
 878 *  RI  A(MMAX)  FACTORS L,D OF THE FACTORIZATION A+E=L*D*TRANS(L)
 879 *                STORED USING THE COMPACT SCHEME OF STORING.
 880 *  II  PSL(N+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX
 881 *  II  SL(MMAX)  ARRAY OF COLUMN INDICES OF THE FACTORS L AND D
 882 *         STORED USING THE COMPACT SCHEME.
 883 *  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
 884 *         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
 885 *         EQUATIONS.
 886 *  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN
 887 *         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
 888 *
 889 * METHOD :
 890 * BACK SUBSTITUTION
 891 *
 892       SUBROUTINE MXSPCB(N,A,PSL,SL,X,JOB)
 893       INTEGER N
 894       DOUBLE PRECISION  A(*),X(*)
 895       INTEGER PSL(*),SL(*),JOB
 896       INTEGER I,J,IS
 897 *
 898 *     FIRST PHASE
 899 *
 900       IF (JOB.GE.0) THEN
 901                     DO 300 I=1,N
 902                     IS=SL(I)+N+1
 903                     DO 200 J=PSL(I)+I,PSL(I+1)+I-1
 904                     X(SL(IS))=X(SL(IS)) - A(J)*X(I)
 905                     IS=IS+1
 906 200                 CONTINUE
 907 300                 CONTINUE
 908       END IF
 909 *
 910 *     SECOND PHASE
 911 *
 912       IF (JOB.EQ.0) THEN
 913                     DO 400 I=1,N
 914                     X(I) = X(I)/A(PSL(I)+I-1)
 915 400   CONTINUE
 916       END IF
 917 *
 918 *     THIRD PHASE
 919 *
 920       IF (JOB.LE.0) THEN
 921                     DO 600 I=N,1,-1
 922                     IS=SL(I)+N+1
 923                     DO 500 J=PSL(I)+I,PSL(I+1)+I-1
 924                     X(I)=X(I)-A(J)*X(SL(IS))
 925                     IS=IS+1
 926 500                 CONTINUE
 927 600                 CONTINUE
 928       END IF
 929       RETURN
 930       END
 931 * SUBROUTINE MXSPCC                ALL SYSTEMS               92/12/01
 932 * PURPOSE :
 933 * SPARSE MATRIX REORDER, SYMBOLIC FACTORIZATION, DATA STRUCTURES
 934 * TRANSFORMATION - INITIATION OF THE DIRECT SPARSE SOLVER.
 935 * MODIFIED VERSION WITH CHANGED DATA STRUCTURES.
 936 *
 937 *  PARAMETERS :
 938 *  II  N  ACTUAL NUMBER OF VARIABLES.
 939 *  II  NJA  NUMBER OF NONZERO ELEMENTS OF THE MATRIX.
 940 *  IO  ML  SIZE OF THE COMPACT STRUCTURE OF THE TRIANGULAR FACTOR
 941 *         OF THE HESSIAN APPROXIMATION.
 942 *  II  MMAX  SIZE OF THE ARRAYS JA,A.
 943 *  RO  A(MMAX)   NUMERICAL VALUES OF THE SPARSE HESSIAN APPROXIMATION
 944 *         STORED AT THE POSITIONS 1, ...,NJA. LOWER TRIANGULAR
 945 *         FACTOR OF THE HESSIAN APPROXIMATION STORED AT THE
 946 *         POSITIONS NJA+1, ..., NJA+MB.
 947 *  II  IA(N) POINTERS OF THE DIAGONAL ELEMENTS OF THE HESSIAN MATRIX.
 948 *  II  JA(MMAX)  INDICES OF THE NONZERO ELEMENTS OF THE HESSIAN MATRIX IN
 949 *         THE PACKED ROW FORM AT THE FIRST NJA POSITIONS. COMPACT
 950 *         STRUCTURE OF INDICES OF ITS TRIANGULAR FACTOR IS ROWWISE
 951 *         STORED.
 952 *  II  PSL(N+1)  POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR
 953 *         FACTOR OF THE HESSIAN APPROXIMATION.
 954 *  IO  PERM(N) PERMUTATION VECTOR.
 955 *  IO  INVP(N) INVERSE PERMUTATION VECTOR.
 956 *  IA  WN11(N) AUXILIARY VECTOR.
 957 *  IA  WN12(N) AUXILIARY VECTOR.
 958 *  IA  WN13(N) AUXILIARY VECTOR.
 959 *  IA  WN14(N) AUXILIARY VECTOR.
 960 *  IO  ITERM  TERMINATION INDICATOR. TERMINATION IF ITERM .NE. 0.
 961 *
 962 * COMMON DATA :
 963 *
 964 * SUBPROGRAMS USED :
 965 *  S   MXSTG1  WIDTHENING OF THE PACKED FORM OF THE SPARSE MATRIX.
 966 *  S   MXSSMN  SPARSE MATRIX REORDERING.
 967 *  S   MXVSIP  INVERSE PERMUTATION COMPUTING.
 968 *  S   MXSPCI  SYMBOLIC FACTORIZATION.
 969 *  S   MXSTL1  PACKING OF THE WIDTHENED FORM OF THE SPARSE MATRIX.
 970 *
 971       SUBROUTINE MXSPCC(N,NJA,ML,MMAX,A,IA,JA,PSL,PERM,INVP,WN11,WN12,
 972      *                  WN13,WN14,ITERM)
 973       INTEGER N,NJA,MMAX,ML,ITERM
 974       INTEGER PERM(*),INVP(*),WN11(*),WN12(*),WN13(*),WN14(*)
 975       INTEGER PSL(*),IA(*),JA(*)
 976       INTEGER JSTRT,JSTOP,I,J,K,L,NJASAVE
 977       INTEGER LL,LL1,NJABIG,KSTRT
 978       DOUBLE PRECISION A(*)
 979       IF (ML.GT.0) RETURN
 980       IF (2*NJA.GE.MMAX) THEN
 981         ITERM=-41
 982         GO TO 1000
 983       END IF
 984 *
 985 *     WIDTHENING OF THE PACKED FORM
 986 *
 987       NJASAVE=NJA
 988       CALL MXSTG1(N,NJA,IA,JA,WN12,WN11)
 989       NJABIG=NJA
 990 *
 991 *     REORDERING OF THE SPARSE MATRIX
 992 *
 993       CALL MXSSMN(N,IA,JA,PERM,WN11,WN12,WN13)
 994 *
 995 *     FIND THE INVERSE PERMUTATION VECTOR INVP
 996 *
 997       CALL MXVSIP(N,PERM,INVP)
 998 *
 999 *     SHRINK THE STRUCTURE
1000 *
1001       CALL MXSTL1(N,NJA,IA,JA,WN11)
1002       DO 40 I=1,N
1003         WN11(I)=0
1004         WN12(I)=0
1005  40   CONTINUE
1006 *
1007 *     WN11 CONTAINS BEGINNINGS OF THE FACTOR ROWS
1008 *
1009       DO 50 I=1,N
1010         K=PERM(I)
1011         JSTRT=IA(K)
1012         JSTOP=IA(K+1)-1
1013         DO 60 J=JSTRT,JSTOP
1014           L=JA(J)
1015           L=INVP(L)
1016           IF (L.GE.I) THEN
1017             WN12(I)=WN12(I)+1
1018           ELSE
1019             WN12(L)=WN12(L)+1
1020           END IF
1021  60     CONTINUE
1022  50   CONTINUE
1023       WN11(1)=1
1024       DO 69 I=1,N-1
1025        WN11(I+1)=WN11(I)+WN12(I)
1026  69   CONTINUE
1027 *
1028 *     CREATE UPPER TRIANGULAR STRUCTURE NECESSARY FOR THE TRANSFER
1029 *
1030       DO 300 I=1,N
1031         K=PERM(I)
1032         JSTRT=IA(K)
1033         JSTOP=IA(K+1)-1
1034         DO 200 J=JSTRT,JSTOP
1035           L=JA(J)
1036           L=INVP(L)
1037           IF (L.GE.I) THEN
1038             LL1=WN11(I)
1039             WN11(I)=LL1+1
1040             JA(NJABIG+LL1)=L
1041             A(J)=LL1
1042             A(NJA+LL1)=J
1043           ELSE
1044             LL1=WN11(L)
1045             WN11(L)=LL1+1
1046             JA(NJABIG+LL1)=I
1047             A(J)=LL1
1048             A(NJA+LL1)=J
1049           END IF
1050  200    CONTINUE
1051  300  CONTINUE
1052 *
1053 *     SORT INDICES IN THE PERMUTED UPPER TRIANGLE
1054 *
1055       DO 599 I=1,N
1056         WN11(I)=0
1057  599  CONTINUE
1058       WN11(1)=1
1059       WN14(1)=1
1060       DO 67 I=2,N+1
1061         WN11(I)=WN11(I-1)+WN12(I-1)
1062         WN14(I)=WN11(I)
1063  67   CONTINUE
1064       DO 602 I=1,N
1065         WN12(I)=0
1066  602  CONTINUE
1067        JSTOP=WN11(N+1)
1068        DO 499 I=N,1,-1
1069          JSTRT=WN11(I)
1070          CALL MXVSR5(JSTOP-JSTRT,JSTRT-1,JA(NJABIG+JSTRT),
1071      &               A,A(NJASAVE+JSTRT))
1072          JSTOP=JSTRT
1073  499   CONTINUE
1074 *
1075 *     WIDTHENING OF THE PERMUTED PACKED FORM.
1076 *
1077       NJASAVE=NJA
1078       CALL MXSTG1(N,NJA,IA,JA,WN12,WN11)
1079       NJABIG=NJA
1080 *
1081 *     SYMBOLIC FACTORIZATION.
1082 *
1083       CALL MXSPCI(N,ML,MMAX-2*NJA,IA,JA,PSL,A(2*NJASAVE+1),PERM,
1084      &            INVP,WN11,WN12,WN13,ITERM)
1085       IF (ITERM.NE.0) THEN
1086         ITERM=-42
1087         GO TO 1000
1088       END IF
1089 *
1090 *     RETRIEVE PARAMETERS
1091 *
1092       CALL MXSTL1(N,NJA,IA,JA,WN11)
1093 *
1094 *     SHIFT PERMUTED UPPER TRIANGLE.
1095 *
1096       DO 73 I=1,NJASAVE
1097         JA(NJA+I)=JA(NJABIG+I)
1098   73  CONTINUE
1099 *
1100 *     SHIFT STRUCTURE SL.
1101 *
1102       IF (2*NJASAVE+ML.GE.MMAX) THEN
1103         ITERM=-41
1104         GO TO 1000
1105       END IF
1106       DO 70 I=1,ML
1107         JA(2*NJASAVE+I)=A(2*NJASAVE+I)
1108   70  CONTINUE
1109 *
1110 *     SET POINTERS
1111 *
1112       DO 600 I=1,N
1113         WN12(I)=0
1114  600  CONTINUE
1115        LL1=PSL(N)+N-1
1116        JSTOP=WN14(N+1)
1117        DO 500 I=N,1,-1
1118          JSTRT=WN14(I)
1119          DO 700 J=JSTRT,JSTOP-1
1120            K=JA(NJASAVE+J)
1121            WN12(K)=J
1122            LL=A(NJASAVE+J)
1123            WN13(K)=LL
1124  700     CONTINUE
1125          JSTOP=JSTRT
1126          KSTRT=JA(2*NJASAVE+I)+N+1+2*NJASAVE
1127          DO 800 J=KSTRT+PSL(I+1)-PSL(I)-1,KSTRT,-1
1128            L=JA(J)
1129            IF (WN12(L).NE.0) THEN
1130              LL=WN13(L)
1131              A(LL)=LL1
1132              WN12(L)=0
1133            END IF
1134            LL1=LL1-1
1135  800     CONTINUE
1136          K=WN12(I)
1137          WN12(I)=0
1138          LL=WN13(I)
1139          A(LL)=LL1
1140          LL1=LL1-1
1141  500   CONTINUE
1142       DO 76 I=1,ML
1143         JA(NJASAVE+I)=JA(2*NJASAVE+I)
1144   76  CONTINUE
1145       DO 72 I=1,NJASAVE
1146         JA(ML+NJASAVE+I)=A(I)
1147   72  CONTINUE
1148  1000 CONTINUE
1149       RETURN
1150       END
1151 * SUBROUTINE MXSPCD                ALL SYSTEMS                92/12/01
1152 * PURPOSE :
1153 * COMPUTATION OF A DIRECTION OF NEGATIVE CURVATURE WITH RESPECT TO A
1154 * SPARSE SYMMETRIC MATRIX A USING THE FACTORIZATION A+E=L*D*TRANS(L)
1155 * STORED IN THE COMPACT SPARSE FORM.
1156 *
1157 * PARAMETERS :
1158 *  II  N ORDER OF THE MATRIX A.
1159 *  II  MMAX  LENGTH OF THE PRINCIPAL MATRIX VECTORS (SL,A).
1160 *  RI  A(MMAX)      FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
1161 *         SUBROUTINE MXSPGF.IT CONTAINS THE NUMERICAL VALUES OF THE
1162 *         FACTORS STORED IN THE COMPACT FORM ACCORDING TO THE
1163 *         INFORMATION IN THE VECTORS PSL,SL.
1164 *  II  PSL(N+1)  POINTER VECTOR OF THE FACTORIZED MATRIX A.
1165 *  II  SL(MMAX)  COMPACT SHEME OF THE FACTORIZED MATRIX A.
1166 *  RO  X(N)  COMPUTED DIRECTION OF NEGATIVE CURVATURE (I.E.
1167 *         TRANS(X)*A*X<0) IF IT EXISTS.
1168 *  II  INF  INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. THE
1169 *         DIRECTION OF NEGATIVE CURVATURE EXISTS ONLY IF INF>0.
1170 *
1171 * METHOD :
1172 * P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND
1173 * LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974)
1174 * PP. 311-350.
1175 *
1176       SUBROUTINE MXSPCD(N,A,PSL,SL,X,INF)
1177       INTEGER  N,INF,PSL(*),SL(*)
1178       DOUBLE PRECISION  A(*),X(*)
1179       INTEGER  I, J, IS
1180 *
1181 *     RIGHT HAND SIDE FORMATION
1182 *
1183       DO 100 I=1,N
1184       X(I) = 0.0D 0
1185 100   CONTINUE
1186       IF (INF .LE. 0) RETURN
1187       X(INF) = 1.0D 0
1188 *
1189 *     BACK SUBSTITUTION
1190 *
1191       DO 300 I=INF-1,1,-1
1192       IS=SL(I)+N+1
1193       DO 200 J=PSL(I)+I,PSL(I+1)+I-1
1194       X(I)=X(I)-A(J)*X(SL(IS))
1195       IS=IS+1
1196 200   CONTINUE
1197 300   CONTINUE
1198       RETURN
1199       END
1200 * SUBROUTINE MXSPCF                ALL SYSTEMS                90/12/01
1201 * PURPOSE :
1202 * NUMERICAL  FACTORIZATION A+E=L*D*TRANS(L) OF A SPARSE
1203 * SYMMETRIC POSITIVE DEFINITE MATRIX A+E WHERE D AND E ARE DIAGONAL
1204 * POSITIVE DEFINITE MATRICES AND L IS A LOWER TRIANGULAR MATRIX. IF
1205 * A IS SUFFICIENTLY POSITIVE DEFINITE THEN E=0. THE STRUCTURE ON
1206 * INPUT WAS OBTAINED BY THE SYMBOLIC FACTORIZATION AND IT MAKES
1207 * USE OF THE COMPACT SCHEME OF STORING THE SPARSE MATRIX IN THE
1208 * POINTER ARRAY PSL ,INDEX ARRAY SL. NUMERICAL VALUES OF THE FACTOR
1209 * CAN BE FOUND IN THE ARRAY A.
1210 *
1211 * PARAMETERS :
1212 *  II  N ORDER OF THE MATRIX A.
1213 *  RU  A(MMAX) ON INPUT NUMERICAL VALUES OF THE LOWER HALF OF THE
1214 *      MATRIX THAT IS BEEING FACTORIZED(INCLUDING THE DIAGONAL
1215 *      ELEMENTS. ON OUTPUT IT CONTAINS FACTORS L AND D AS IF THEY
1216 *      FORM THE LOWER HALF OF THE MATRIX.STRUCTURE INFORMATION
1217 *      IS SAVED IN THE COMPACT SCHEME IN THE PAIR OF VECTORS PSL
1218 *      AND SL.
1219 *  II  PSL(NF+1) POINTER VECTOR OF THE FACTOR
1220 *  II  SL(MMAX) STRUCTURE OF THE FACTOR IN THE COMPACT FORM
1221 *  IA  WN11(NF+1) AUXILIARY VECTOR.
1222 *  IA  WN12(NF+1) AUXILIARY VECTOR.
1223 *  RA  RN01(NF+1) AUXILIARY VECTOR.
1224 *  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
1225 *         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
1226 *         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
1227 *         INF>0 THEN THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
1228 *         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
1229 *         PROCESS.
1230 *  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON
1231 *         OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
1232 *         FACTORIZATION PROCESS (IF INF>0).
1233 *  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
1234 *
1235 * METHOD :
1236 * S.C.EISENSTAT,M.C.GURSKY,M.H.SCHULTZ,A.H.SHERMAN:YALE SPARSE MATRIX
1237 * PACKAGE I. THE SYMMETRIC CODES,YALE UNIV. RES. REPT.
1238 * NO.112,1977.
1239 *
1240       SUBROUTINE MXSPCF(N,A,PSL,SL,WN11,WN12,RN01,INF,ALF,TAU)
1241       INTEGER N,PSL(*),SL(*),WN11(*),WN12(*),INF
1242       DOUBLE PRECISION A(*),RN01(*),ALF
1243       DOUBLE PRECISION BET,GAM,DEL,RHO,SIG,TOL,TADD,TBDD,TAU
1244       INTEGER  I, J, K, L, II
1245       INTEGER ISTRT,ISTOP,NEWK,KPB,ISUB
1246       L = 0
1247       INF = 0
1248       TOL = ALF
1249       ALF = 0.0D 0
1250       BET = 0.0D 0
1251       GAM = 0.0D 0
1252       TAU = 0.0D 0
1253       DO 60 I=1,N
1254       BET=MAX(BET,ABS(A(PSL(I)+I-1)))
1255       DO 50 J=PSL(I)+I,PSL(I+1)+I-1
1256       GAM = MAX( GAM, ABS(A(J)) )
1257 50    CONTINUE
1258 60    CONTINUE
1259       BET = MAX(TOL,BET,GAM/N)
1260       DEL = TOL*BET
1261       DO 100 I=1,N
1262       WN11(I)=0
1263       RN01(I)=0.0D 0
1264 100   CONTINUE
1265       DO 600 J=1,N
1266 *
1267 *     DETERMINATION OF A DIAGONAL CORRECTION
1268 *
1269       SIG=A(PSL(J)+J-1)
1270       RHO=0.0D 0
1271       NEWK=WN11(J)
1272 200   K=NEWK
1273       IF (K.EQ.0) GO TO 400
1274       NEWK=WN11(K)
1275       KPB=WN12(K)
1276       TADD=A(KPB+K)
1277       TBDD=TADD*A(PSL(K)+K-1)
1278       RHO=RHO+TADD*TBDD
1279       ISTRT=KPB+1
1280       ISTOP=PSL(K+1)-1
1281       IF (ISTOP.LT.ISTRT) GO TO 200
1282       WN12(K)=ISTRT
1283       I=SL(K)+(KPB-PSL(K))+1
1284       ISUB=SL(N+1+I)
1285       WN11(K)=WN11(ISUB)
1286       WN11(ISUB)=K
1287       DO 300 II=ISTRT,ISTOP
1288       ISUB=SL(N+1+I)
1289       RN01(ISUB)=RN01(ISUB)+A(II+K)*TBDD
1290       I=I+1
1291 300   CONTINUE
1292       GO TO 200
1293 400   SIG=A(PSL(J)+J-1)-RHO
1294       IF (ALF.GT.SIG) THEN
1295         ALF=SIG
1296         L=J
1297       END IF
1298       GAM=0.0D 0
1299       ISTRT=PSL(J)
1300       ISTOP=PSL(J+1)-1
1301       IF (ISTOP.LT.ISTRT) GO TO 370
1302       WN12(J)=ISTRT
1303       I=SL(J)
1304       ISUB=SL(N+1+I)
1305       WN11(J)=WN11(ISUB)
1306       WN11(ISUB)=J
1307       DO 500 II=ISTRT,ISTOP
1308       ISUB=SL(N+1+I)
1309       A(II+J)=(A(II+J)-RN01(ISUB))
1310       RN01(ISUB)=0.0D 0
1311       I=I+1
1312 500   CONTINUE
1313       DO 350 K=PSL(J)+J,PSL(J+1)+J-1
1314       GAM=MAX(GAM,ABS(A(K)))
1315 350   CONTINUE
1316       GAM=GAM*GAM
1317 370   RHO=MAX(ABS(SIG),GAM/BET,DEL)
1318       IF (TAU.LT.RHO-SIG) THEN
1319                          TAU=RHO-SIG
1320                          INF=-1
1321       END IF
1322 *
1323 *     GAUSSIAN ELIMINATION
1324 *
1325       A(PSL(J)+J-1)=RHO
1326       DO 550 II=ISTRT,ISTOP
1327       A(II+J)=A(II+J)/RHO
1328 550   CONTINUE
1329 600   CONTINUE
1330       IF (L.NE.0.AND.ABS(ALF).GT.DEL) INF=L
1331       RETURN
1332       END
1333 * SUBROUTINE MXSPCI                ALL SYSTEMS                89/12/01
1334 * PURPOSE :
1335 * SYMBOLIC FACTORIZATION OF A SPARSE SYMMETRIC MATRIX GIVEN IN THE
1336 * NORMAL SCHEME PA,SA. ON OUTPUT WE HAVE POINTER VECTOR OF THE FACTOR
1337 * PSL AND VECTOR OF COLUMN INDICES SL. ML IS THE NUMBER OF THE INDICES
1338 * USED FOR THE VECTOR SL, WHERE WE HAVE FACTOR IN THE COMPACT FORM.
1339 *
1340 * PARAMETERS :
1341 *  II  N ORDER OF THE MATRIX A.
1342 *  IO  ML NUMBER OF THE NONZERO ELEMENTS IN THE FACTOR'S COMPACT SCHEME
1343 *  II  MMAX  LENGTH OF THE ARRAY SL. IN THE CASE OF THE
1344 *            INSUFFICIENT SPACE IT IS TO BE INCREASED.
1345 *  II  PA(N+1) POINTER VECTOR OF THE INPUT MATRIX
1346 *  II  SA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX
1347 *  IO  PSL(N+1) POINTER VECTOR OF THE FACTOR
1348 *  RO  SL(MMAX) COMPACT SCHEME OF THE INDICES OF THE FACTOR
1349 *  II  PERM(N) PERMUTATION VECTOR
1350 *  II  INVP(N) INVERSE PERMUTATION VECTOR
1351 *  IA  WN11(N+1) WORK VECTOR OF THE LENGTH N+1
1352 *  IA  WN12(N+1) WORK VECTOR OF THE LENGTH N+1
1353 *  IA  WN13(N+1) WORK VECTOR OF THE LENGTH N+1
1354 *  IO  ISPACE AN INFORMATION ON SPACE OBTAINED DURING THE PROCESS
1355 *       OF THE FACTORIZATION.
1356 *      ISPACE=0    THE FACTORIZATION HAS TERMINATED NORMALLY
1357 *      ISPACE=1    INSUFFICIENT SPACE AVAILABLE
1358 *
1359 * METHOD :
1360 * S.C.EISENSTAT,M.C.GURSKY,M.H.SCHULTZ,A.H.SHERMAN:YALE SPARSE MATRIX
1361 * PACKAGE I. THE SYMMETRIC CODES,ACM TRANS. ON MATH. SOFTWARE.
1362 *
1363 * NOTE: TYPE OF SL CHANGED FOR THE UFO APPLICATION.
1364 *
1365       SUBROUTINE MXSPCI(N,ML,MMAX,PA,SA,PSL,SL,PERM,INVP,
1366      &            WN11,WN12,WN13,ISPACE)
1367       INTEGER N,MMAX,PA(*),SA(*),PSL(*)
1368       INTEGER PERM(*),INVP(*),WN11(*),WN12(*),WN13(*)
1369       INTEGER ISPACE,I,INZ,J,JSTOP,JSTRT,K,KNZ,KXSUB,MRGK,LMAX,ML
1370       INTEGER NABOR,NODE,NP1,NZBEG,NZEND,RCHM,MRKFLG,M
1371       DOUBLE PRECISION SL(*)
1372       NZBEG=1
1373       NZEND=0
1374       PSL(1)=1
1375       DO 100 K=1,N
1376       WN11(K)=0
1377       WN13(K)=0
1378 100   CONTINUE
1379       NP1=N+1
1380       DO 1500 K=1,N
1381       KNZ=0
1382       MRGK=WN11(K)
1383       MRKFLG=0
1384       WN13(K)=K
1385       IF (MRGK.NE.0) WN13(K)=WN13(MRGK)
1386       SL(K)=NZEND
1387       NODE=PERM(K)
1388       JSTRT=PA(NODE)
1389       JSTOP=PA(NODE+1)-1
1390       IF (JSTRT.GT.JSTOP) GO TO 1500
1391       WN12(K)=NP1
1392       DO 300 J=JSTRT,JSTOP
1393       NABOR=SA(J)
1394       IF (NABOR.EQ.NODE) GO TO 300
1395       NABOR=INVP(NABOR)
1396       IF (NABOR.LE.K) GO TO 300
1397       RCHM=K
1398 200   M=RCHM
1399       RCHM=WN12(M)
1400       IF (RCHM.LE.NABOR) GO TO 200
1401       KNZ=KNZ+1
1402       WN12(M)=NABOR
1403       WN12(NABOR)=RCHM
1404       IF (WN13(NABOR).NE.WN13(K)) MRKFLG=1
1405 300   CONTINUE
1406       LMAX=0
1407       IF (MRKFLG.NE.0.OR.MRGK.EQ.0) GO TO 350
1408       IF (WN11(MRGK).NE.0) GO TO 350
1409       SL(K)=SL(MRGK)+1
1410       KNZ=PSL(MRGK+1)-(PSL(MRGK)+1)
1411       GO TO 1400
1412 350   I=K
1413 400   I=WN11(I)
1414       IF (I.EQ.0) GO TO 800
1415       INZ=PSL(I+1)-(PSL(I)+1)
1416       JSTRT=SL(I)+1
1417       JSTOP=SL(I)+INZ
1418       IF (INZ.LE.LMAX) GO TO 500
1419       LMAX=INZ
1420       SL(K)=JSTRT
1421 500   RCHM=K
1422       DO 700 J=JSTRT,JSTOP
1423       NABOR=SL(N+1+J)
1424 600   M=RCHM
1425       RCHM=WN12(M)
1426       IF (RCHM.LT.NABOR) GO TO 600
1427       IF (RCHM.EQ.NABOR) GO TO 700
1428       KNZ=KNZ+1
1429       WN12(M)=NABOR
1430       WN12(NABOR)=RCHM
1431       RCHM=NABOR
1432 700   CONTINUE
1433       GO TO 400
1434 800   IF (KNZ.EQ.LMAX) GO TO 1400
1435       IF (NZBEG.GT.NZEND) GO TO 1200
1436       I=WN12(K)
1437       DO 900 JSTRT=NZBEG,NZEND
1438       IF (SL(N+1+JSTRT)-I .GE.0) THEN
1439         IF (SL(N+1+JSTRT).EQ.I) THEN
1440           GO TO 1000
1441         ELSE
1442           GO TO 1200
1443         END IF
1444       END IF
1445 900   CONTINUE
1446       GO TO 1200
1447 1000  SL(K)=JSTRT
1448       DO 1100 J=JSTRT,NZEND
1449       IF (SL(N+1+J).NE.I) GO TO 1200
1450       I=WN12(I)
1451       IF (I.GT.N) GO TO 1400
1452 1100  CONTINUE
1453       NZEND=JSTRT-1
1454 1200  NZBEG=NZEND+1
1455       NZEND=NZEND+KNZ
1456 *
1457 *     A VARIANT IS USED WHEN CALLED SO THAT SL(X)=A(NB+X)
1458 *
1459       IF (NZEND.GE.MMAX-N-1) GO TO 1600
1460       I=K
1461       DO 1300 J=NZBEG,NZEND
1462       I=WN12(I)
1463       SL(N+1+J)=I
1464       WN13(I)=K
1465 1300  CONTINUE
1466       SL(K)=NZBEG
1467       WN13(K)=K
1468 1400     IF (KNZ.GT.1) THEN
1469             KXSUB=SL(K)
1470             I=SL(N+1+KXSUB)
1471             WN11(K)=WN11(I)
1472             WN11(I)=K
1473       END IF
1474       PSL(K+1)=PSL(K)+KNZ
1475 1500  CONTINUE
1476       SL(N)=SL(N)+1
1477       SL(N+1)=SL(N)
1478       ML=N+SL(N+1)
1479       ISPACE=0
1480       RETURN
1481 1600  ISPACE=1
1482       RETURN
1483       END
1484 * SUBROUTINE MXSPCM                ALL SYSTEMS                92/12/01
1485 * PURPOSE :
1486 * MULTIPLICATION OF A GIVEN VECTOR X BY A SPARSE SYMMETRIC POSITIVE
1487 * DEFINITE MATRIX A+E USING THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED
1488 * BY THE SUBROUTINE MXSPGN. FACTORS ARE STORED IN THE COMPACT FORM.
1489 *
1490 * PARAMETERS :
1491 *  II  N ORDER OF THE MATRIX A.
1492 *  RI  A(MMAX)      FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
1493 *         SUBROUTINE MXSPGN.IT CONTAINS THE NUMERICAL VALUES OF THE
1494 *         FACTORS STORED IN THE COMPACT FORM ACCORDING TO THE
1495 *         INFORMATION IN THE VECTORS PSL,SL.
1496 *  II  PSL(N+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX
1497 *  II  SL(MMAX) INDICES OF THE COMPACT SCHEME OF THE FACTORS.
1498 *  RU  X(N)  ON INPUT THE GIVEN VECTOR, ON OUTPUT THE RESULT
1499 *         OF MULTIPLICATION.
1500 *  RA  RN01(N) AUXILIARY VECTOR.
1501 *  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)*X. IF JOB>0 THEN
1502 *         X:=TRANS(L)*X. IF JOB<0 THEN X:=L*X.
1503 *
1504       SUBROUTINE MXSPCM(N,A,PSL,SL,X,RN01,JOB)
1505       INTEGER N
1506       INTEGER PSL(*),SL(*),JOB
1507       DOUBLE PRECISION  A(*),X(*),RN01(*),ZERO
1508       PARAMETER(ZERO=0.0D0)
1509       INTEGER I,J,IS
1510       DO 50 I=1,N
1511       RN01(I)=ZERO
1512  50   CONTINUE
1513 *
1514 *     FIRST PHASE:X=TRANS(L)*X
1515 *
1516       IF (JOB.GE.0) THEN
1517                     DO 300 I=1,N
1518                     IS=SL(I)+N+1
1519                     DO 200 J=PSL(I)+I,PSL(I+1)+I-1
1520                     X(I)=X(I)+A(J)*X(SL(IS))
1521                     IS=IS+1
1522 200                 CONTINUE
1523 300                 CONTINUE
1524       END IF
1525 *
1526 *     SECOND PHASE:X=D*X
1527 *
1528       IF (JOB.EQ.0) THEN
1529                     DO 400 I=1,N
1530                     X(I) = X(I)*A(PSL(I)+I-1)
1531 400   CONTINUE
1532       END IF
1533 *
1534 *     THIRD PHASE:X=L*X
1535 *
1536       IF (JOB.LE.0) THEN
1537                     DO 600 I=N,1,-1
1538                     IS=SL(I)+N+1
1539                     DO 500 J=PSL(I)+I,PSL(I+1)+I-1
1540                     RN01(SL(IS))=RN01(SL(IS))+A(J)*X(I)
1541                     IS=IS+1
1542 500                 CONTINUE
1543 600                 CONTINUE
1544                     DO 700 I=1,N
1545                     X(I)=RN01(I)+X(I)
1546 700                 CONTINUE
1547       END IF
1548       RETURN
1549       END
1550 * SUBROUTINE MXSPCN                ALL SYSTEMS               93/12/01
1551 * PURPOSE :
1552 *  ESTIMATION OF THE MINIMUM EIGENVALUE AND THE CORRESPONDING EIGENVECTOR
1553 *  OF A SPARSE SYMMETRIC POSITIVE DEFINITE MATRIX A+E USING THE
1554 *  FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXSPCF.
1555 *
1556 * PARAMETERS :
1557 *  II  N ORDER OF THE MATRIX A.
1558 *  RI  A(MMAX)  FACTORS L,D OF THE FACTORIZATION A+E=L*D*TRANS(L)
1559 *                STORED USING THE COMPACT SCHEME OF STORING.
1560 *  II  PSL(N+1) POINTER ARRAY OF THE FACTORIZED SPARSE MATRIX
1561 *  II  SL(MMAX)  ARRAY OF COLUMN INDICES OF THE FACTORS L AND D
1562 *         STORED USING THE COMPACT SCHEME.
1563 *         SUBROUTINE MXDPGF.
1564 *  RO  X(N)  ESTIMATED EIGENVECTOR.
1565 *  RO  ALF  ESTIMATED EIGENVALUE.
1566 *  II  JOB  OPTION. IF JOB=0 THEN ONLY ESTIMATED EIGENVALUE IS
1567 *         COMPUTED. IS JOB>0 THEN BOTH ESTIMATED EIGENVALUE AND
1568 *         ESTIMATED EIGENVECTOR ARE COMPUTED.
1569 *
1570       SUBROUTINE MXSPCN(N,A,PSL,SL,X,ALF,JOB)
1571       INTEGER N
1572       DOUBLE PRECISION A(*),X(*),ALF
1573       INTEGER PSL(*),SL(*),JOB
1574       DOUBLE PRECISION XP,XM,FP,FM,MXVDOT
1575       INTEGER I,K,IS
1576       DOUBLE PRECISION ZERO,ONE
1577       PARAMETER (ZERO=0.0D 0,ONE=1.0D 0)
1578 *
1579 *     COMPUTATION OF THE VECTOR V WITH POSSIBLE MAXIMUM NORM SUCH
1580 *     THAT  L*D**(1/2)*V=U  WHERE U HAS ELEMENTS +1 OR -1
1581 *
1582       DO 2 I=1,N
1583       X(I)=ZERO
1584     2 CONTINUE
1585       DO 6 K=1,N
1586       XP=-X(K)+ONE
1587       XM=-X(K)-ONE
1588       FP=ABS(XP)
1589       FM=ABS(XM)
1590       IS=SL(K)+N+1
1591       DO 3 I=PSL(K)+K,PSL(K+1)+K-1
1592       FP=FP+ABS(X(SL(IS))+A(I)*XP)
1593       FM=FM+ABS(X(SL(IS))+A(I)*XM)
1594       IS=IS+1
1595     3 CONTINUE
1596       IF (FP.GE.FM) THEN
1597       X(K)=XP
1598       IS=SL(K)+N+1
1599       DO 4 I=PSL(K)+K,PSL(K+1)+K-1
1600       X(SL(IS))=X(SL(IS))+A(I)*XP
1601       IS=IS+1
1602     4 CONTINUE
1603       ELSE
1604       X(K)=XM
1605       IS=SL(K)+N+1
1606       DO 5 I=PSL(K)+K,PSL(K+1)+K-1
1607       X(SL(IS))=X(SL(IS))+A(I)*XM
1608       IS=IS+1
1609     5 CONTINUE
1610       END IF
1611     6 CONTINUE
1612 *
1613 *     COMPUTATION OF THE VECTOR X SUCH THAT
1614 *     D**(1/2)*TRANS(L)*X=V
1615 *
1616       FM=ZERO
1617       DO 7 K=1,N
1618       IF (JOB.LE.0) THEN
1619       FP=SQRT(A(PSL(K)+K-1))
1620       X(K)=X(K)/FP
1621       FM=FM+X(K)*X(K)
1622       ELSE
1623       X(K)=X(K)/A(PSL(K)+K-1)
1624       END IF
1625     7 CONTINUE
1626       FP=DBLE(N)
1627       IF (JOB.LE.0) THEN
1628 *
1629 *     FIRST ESTIMATION OF THE MINIMUM EIGENVALUE BY THE FORMULA
1630 *     ALF=(TRANS(U)*U)/(TRANS(V)*V)
1631 *
1632       ALF=FP/FM
1633       RETURN
1634       END IF
1635       FM=ZERO
1636       DO 9 K=N,1,-1
1637       IS=SL(K)+N+1
1638       DO 8 I=PSL(K)+K,PSL(K+1)+K-1
1639       X(K)=X(K)-A(I)*X(SL(IS))
1640       IS=IS+1
1641     8 CONTINUE
1642       FM=FM+X(K)*X(K)
1643     9 CONTINUE
1644       FM=SQRT(FM)
1645       IF (JOB.LE.1) THEN
1646 *
1647 *     SECOND ESTIMATION OF THE MINIMUM EIGENVALUE BY THE FORMULA
1648 *     ALF=SQRT(TRANS(U)*U)/SQRT(TRANS(X)*X)
1649 *
1650       ALF=SQRT(FP)/FM
1651       ELSE
1652 *
1653 *     INVERSE ITERATIONS
1654 *
1655       DO 11 K=2,JOB
1656 *
1657 *     SCALING THE VECTOR X BY ITS NORM
1658 *
1659       DO 10 I=1,N
1660       X(I)=X(I)/FM
1661    10 CONTINUE
1662       CALL MXSPCB(N,A,PSL,SL,X,0)
1663       FM=SQRT(MXVDOT(N,X,X))
1664    11 CONTINUE
1665       ALF=ONE/FM
1666       END IF
1667 *
1668 *     SCALING THE VECTOR X BY ITS NORM
1669 *
1670       DO 12 I=1,N
1671       X(I)=X(I)/FM
1672    12 CONTINUE
1673       RETURN
1674       END
1675 * FUNCTION MXSPCP                  ALL SYSTEMS                92/12/01
1676 * PURPOSE :
1677 * COMPUTATION OF THE NUMBER MXSPCP=TRANS(X)*D**(-1)*X WHERE D IS A
1678 * DIAGONAL MATRIX IN THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
1679 * SUBROUTINE MXSPGN. THE FACTORS ARE STORED IN THE COMPACT FORM.
1680 *
1681 * PARAMETERS :
1682 *  II  N ORDER OF THE MATRIX A.
1683 *  RI  A(MMAX)      FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
1684 *         SUBROUTINE MXSPGN.IT CONTAINS THE NUMERICAL VALUES OF THE
1685 *         FACTORS STORED IN THE COMPACT FORM ACCORDING TO THE
1686 *         INFORMATION IN THE VECTORS PSL,SL.
1687 *  II  PSL(N+1)  POINTER VECTOR OF THE FACTORIZED MATRIX A.
1688 *  RI  X(N)  INPUT VECTOR
1689 *  RR  MXSPCP  COMPUTED NUMBER MXSPCP=TRANS(X)*D**(-1)*X
1690 *
1691       FUNCTION MXSPCP(N,A,PSL,X)
1692       INTEGER  N
1693       DOUBLE PRECISION  A(*), X(*), MXSPCP
1694       DOUBLE PRECISION  TEMP
1695       INTEGER  PSL(*),I
1696       TEMP = 0.0D 0
1697       DO  100 I=1,N
1698       TEMP = TEMP + X(I)*X(I)/A(PSL(I)+I-1)
1699 100   CONTINUE
1700       MXSPCP = TEMP
1701       RETURN
1702       END
1703 * FUNCTION MXSPCQ                  ALL SYSTEMS                92/12/01
1704 * PURPOSE :
1705 * COMPUTATION OF THE NUMBER MXSPCQ=TRANS(X)*D*X WHERE D IS A
1706 * DIAGONAL MATRIX IN THE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
1707 * SUBROUTINE MXSPGN. FACTORS ARE STORED IN THE COMPACT FORM.
1708 *
1709 * PARAMETERS :
1710 *  II  N ORDER OF THE MATRIX A.
1711 *  RI  A(MMAX)      FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
1712 *         SUBROUTINE MXSPGN.IT CONTAINS THE NUMERICAL VALUES OF THE
1713 *         FACTORS STORED IN THE COMPACT FORM ACCORDING TO THE
1714 *         INFORMATION IN THE VECTORS PSL,SL.
1715 *  II  PSL(N+1)  POINTER VECTOR OF THE FACTORIZED MATRIX A
1716 *  RI  X(N)  INPUT VECTOR
1717 *  RR  MXSPCQ  COMPUTED NUMBER MXSPCQ=TRANS(X)*D*X
1718 *
1719       FUNCTION MXSPCQ(N,A,PSL,X)
1720       INTEGER  N
1721       DOUBLE PRECISION  A(*), X(*), MXSPCQ
1722       DOUBLE PRECISION  TEMP
1723       INTEGER  PSL(N+1),I
1724       DOUBLE PRECISION ZERO
1725       PARAMETER (ZERO = 0.0D0)
1726       TEMP = ZERO
1727       DO  100 I=1,N
1728       TEMP = TEMP + X(I)*X(I)*A(PSL(I)+I-1)
1729 100   CONTINUE
1730       MXSPCQ = TEMP
1731       RETURN
1732       END
1733 * SUBROUTINE MXSPCT                 ALL SYSTEMS                92/12/01
1734 * PURPOSE :
1735 * REWRITE SYMMETRIC MATRIX INTO THE PERMUTED FACTORIZED COMPACT SCHEME.
1736 * MOIDIFIED VERSION FOR THE USE WITH MXSPCJ.
1737 *
1738 * PARAMETERS:
1739 *  II  N  SIZE OF THE SYSTEM SOLVED.
1740 *  II  NB NUMBER OF NONZEROS IN THE UPPER TRIANGLE OF THE ORIGINAL
1741 *         MATRIX.
1742 *  II  ML SIZE OF THE COMPACT FACTOR.
1743 *  II  MMAX DECLARED LENGTH OF THE ARRAYS JA,A.
1744 *  RU  A(MMAX) NUMERICAL VALUES OF THE SPARSE HESSIAN APPROXIMATION
1745 *              STORED AT THE POSITIONS 1, ...,NB.
1746 *  IU  JA(MMAX) INDICES OF THE NONZERO ELEMENTS OF THE HESSIAN MATRIX IN
1747 *             THE PACKED ROW FORM AT THE FIRST NB POSITIONS.
1748 *             NEW POSITIONS
1749 *             IN THE PERMUTED FACTOR STORED IN A(NB+1), ..., A(2*NB),
1750 *             INDICES OF COMPACT SCHEME IN A(2*NB+1), ..., A(2*NB+ML).
1751 *  II  PSL(N+1)  POINTER VECTOR OF THE COMPACT FORM OF THE TRIANGULAR
1752 *             FACTOR OF THE HESSIAN APPROXIMATION.
1753 *  IO  ITERM  ERROR FLAG. IF ITERM < 0 - LACK OF SPACE IN JA.
1754 *
1755 *
1756       SUBROUTINE MXSPCT(N,NB,ML,MMAX,A,JA,PSL,ITERM)
1757       INTEGER N,NB,ML,MMAX,JA(*)
1758       INTEGER PSL(*),ITERM
1759       DOUBLE PRECISION A(*)
1760       INTEGER I,J
1761 *
1762 *     WN11 CONTAINS BEGINNINGS OF THE FACTOR ROWS
1763 *
1764       ITERM=0
1765 *
1766 *     LACK OF SPACE
1767 *
1768       IF (MMAX.LE.NB+PSL(N+1)+N-1) THEN
1769         ITERM=-43
1770         RETURN
1771       END IF
1772       IF (MMAX.LE.2*NB+ML) THEN
1773         ITERM=-44
1774         RETURN
1775       END IF
1776       DO 50 I=NB+1,NB+PSL(N+1)+N-1
1777         A(I)=0.0D 0
1778  50   CONTINUE
1779       DO 100 I=NB+ML+1,2*NB+ML
1780         J=JA(I)
1781         A(NB+J)=A(I-NB-ML)
1782  100  CONTINUE
1783       RETURN
1784       END
1785 * SUBROUTINE MXSPTB                ALL SYSTEMS                94/12/01
1786 * PURPOSE :
1787 * SOLUTION OF A SYSTEM OF LINEAR EQUATIONS WITH A SPARSE SYMMETRIC
1788 * POSITIVE DEFINITE MATRIX A+E USING INCOMPLETE ILUT-TYPE FACTORIZATION
1789 * A+E=L*D*TRANS(L) OBTAINED BY THE SUBROUTINE MXSPTF.
1790 *
1791 * PARAMETERS :
1792 *  II  N ORDER OF THE MATRIX A.
1793 *  RI  A(MMAX) INCOMPLETE FACTORIZATION A+E=L*D*TRANS(L) OBTAINED BY THE
1794 *         SUBROUTINE MXSPTF.
1795 *  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
1796 *  II  JA(MMAX)  INDICES OF THE NONZERO ELEMENTS OF A.
1797 *  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
1798 *         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
1799 *         EQUATIONS.
1800 *  II  JOB  OPTION. IF JOB=0 THEN X:=(A+E)**(-1)*X. IF JOB>0 THEN
1801 *         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
1802 *
1803 * METHOD :
1804 * BACK SUBSTITUTION
1805 *
1806       SUBROUTINE MXSPTB(N,A,IA,JA,X,JOB)
1807       INTEGER N,IA(*),JA(*),JOB
1808       DOUBLE PRECISION A(*),X(*)
1809       INTEGER I,J,K
1810       DOUBLE PRECISION TEMP,SUM
1811 *
1812 *     FIRST PHASE
1813 *
1814       IF (JOB.GE.0) THEN
1815         DO 300 I=1,N
1816           K=IA(I)
1817           IF (K.LE.0) GO TO 300
1818           TEMP=X(I)*A(K)
1819           DO 200 J=IA(I)+1,IA(I+1)-1
1820             K=JA(J)
1821             IF (K.GT.0) X(K)=X(K)-A(J)*TEMP
1822 200       CONTINUE
1823           IF (JOB.EQ.0) X(I)=TEMP
1824 300     CONTINUE
1825       END IF
1826 *
1827 *     THIRD PHASE
1828 *
1829       IF (JOB.LE.0) THEN
1830         DO 600 I=N,1,-1
1831           K=IA(I)
1832           IF (K.LE.0) GO TO 600
1833           SUM=0.0D 0
1834           TEMP=A(K)
1835           DO 500 J=IA(I)+1,IA(I+1)-1
1836             K=JA(J)
1837             IF (K.GT.0) SUM=SUM+A(J)*X(K)
1838 500       CONTINUE
1839           SUM=SUM*TEMP
1840           X(I)=X(I)-SUM
1841 600     CONTINUE
1842       END IF
1843       RETURN
1844       END
1845 * SUBROUTINE MXSPTF                ALL SYSTEMS                03/12/01
1846 * PURPOSE :
1847 * INCOMPLETE CHOLESKY FACTORIZATION A+E=L*D*TRANS(L) OF A SPARSE
1848 * SYMMETRIC POSITIVE DEFINITE MATRIX A+E WHERE D AND E ARE DIAGONAL
1849 * POSITIVE DEFINITE MATRICES AND L IS A LOWER TRIANGULAR MATRIX.
1850 * METHOD IS BASED ON THE SIMPLE IC(0) ALGORITHM WITHOUT DIAGONAL
1851 * COMPENSATION. SPARSE RIGHT-LOOKING IMPLEMENTATION.
1852 *
1853 * PARAMETERS :
1854 *  II  N ORDER OF THE MATRIX A.
1855 *  RI  A(M)  SPARSE SYMMETRIC (USUALLY POSITIVE DEFINITE) MATRIX.
1856 *  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
1857 *  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
1858 *  IA  WN01(N+1)  AMXILIARY ARRAY.
1859 *  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
1860 *         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
1861 *         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
1862 *         INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
1863 *         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
1864 *         PROCESS.
1865 *  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS.
1866 *         ON OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
1867 *         FACTORIZATION PROCESS (IF INF>0).
1868 *  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
1869 *
1870 * METHOD :
1871 * P.E.GILL, W.MURRAY : NEWTON TYPE METHODS FOR UNCONSTRAINED AND
1872 * LINEARLY CONSTRAINED OPTIMIZATION, MATH. PROGRAMMING 28 (1974)
1873 * PP. 311-350.
1874 *
1875       SUBROUTINE MXSPTF(N,A,IA,JA,WN01,INF,ALF,TAU)
1876       INTEGER N,IA(*),JA(*),WN01(*),INF
1877       DOUBLE PRECISION A(*),ALF,TAU
1878       INTEGER I,J,K,L,II,LL,K1,L1,L2,JSTRT,JSTOP,IDIAG,KSTRT,KSTOP,NJA
1879       DOUBLE PRECISION PTOL,BET,GAM,TEMP,DEL,DIAG,NDIAG,INVDIAG
1880       PTOL=ALF
1881       NJA=IA(N+1)-1
1882 *
1883 *     INITIALIZE AMXILIARY VECTOR
1884 *
1885       INF=0
1886       CALL MXVINS(N,0,WN01)
1887 *
1888 *     GILL-MURRAY MODIFICATION
1889 *
1890       ALF=0.0D 0
1891       BET=0.0D 0
1892       GAM=0.0D 0
1893       TAU=0.0D 0
1894       DO 2 I=1,N
1895         IDIAG=IA(I)
1896         IF (JA(IDIAG).LE.0) GO TO 2
1897         TEMP=A(IDIAG)
1898         BET=MAX(BET,ABS(TEMP))
1899         DO 1 J=IA(I)+1,IA(I+1)-1
1900           IF (JA(J).LE.0) GO TO 1
1901           TEMP=A(J)
1902           GAM=MAX(GAM,ABS(TEMP))
1903     1   CONTINUE
1904     2 CONTINUE
1905       BET=MAX(PTOL,BET,GAM/DBLE(N))
1906       DEL=PTOL*BET
1907 *
1908 *     COMPUTE THE PRECONDITIONER
1909 *
1910       LL=0
1911       DO 8 K=1,N
1912         KSTRT=IA(K)
1913         KSTOP=IA(K+1)-1
1914         IF (JA(KSTRT).LE.0) GO TO 8
1915         DIAG=A(KSTRT)
1916         IF (ALF.GT.DIAG) THEN
1917           ALF=DIAG
1918           LL=K
1919         END IF
1920         GAM=0.0D 0
1921         DO 3 J=KSTRT+1,KSTOP
1922           IF (JA(J).LE.0) GO TO 3
1923           TEMP=A(J)
1924           GAM=MAX(GAM,ABS(TEMP))
1925     3   CONTINUE
1926         GAM=GAM*GAM
1927         INVDIAG=MAX(ABS(DIAG),GAM/BET,DEL)
1928         IF (TAU.LT.INVDIAG-DIAG) THEN
1929           TAU=INVDIAG-DIAG
1930           INF=-1
1931         END IF
1932         INVDIAG=1.0D 0/INVDIAG
1933         A(KSTRT)=INVDIAG
1934 *
1935 *     RIGHT-LOOKING UPDATE
1936 *
1937 *
1938 *     SET POINTERS
1939 *
1940         DO 4 II=KSTRT,KSTOP
1941           K1=JA(II)
1942           IF (K1.GT.0) WN01(K1)=II
1943     4   CONTINUE
1944 *
1945 *     INNER LOOP
1946 *
1947         DO 6 I=KSTRT+1,KSTOP
1948           J=JA(I)
1949           IF (J.LE.0) GO TO 6
1950           NDIAG=A(I)
1951           JSTRT=IA(J)
1952           JSTOP=IA(J+1)-1
1953           DO 5 L=JSTRT,JSTOP
1954             L1=JA(L)
1955             IF (L1.LE.0) GO TO 5
1956             L2=WN01(L1)
1957             IF (L2.NE.0) THEN
1958               A(L)=A(L)-(A(L2)*INVDIAG)*NDIAG
1959             END IF
1960     5     CONTINUE
1961     6   CONTINUE
1962 *
1963 *     CLEAR THE POINTERS
1964 *
1965         DO 7 II=KSTRT,KSTOP
1966           K1=JA(II)
1967           IF (K1.GT.0) WN01(K1)=0
1968     7   CONTINUE
1969     8 CONTINUE
1970       IF (LL.GT.0.AND.ABS(ALF).GT.DEL) INF=LL
1971       RETURN
1972       END
1973 * SUBROUTINE MXSRMD               ALL SYSTEMS                92/12/01
1974 * PURPOSE :
1975 * MULTIPLICATION OF TRANSPOSE OF A DENSE RECTANGULAR MATRIX A BY
1976 * A VECTOR X AND ADDITION OF A SCALED VECTOR ALF*Y.
1977 *
1978 * PARAMETERS :
1979 *  II  N  NUMBER OF ROWS OF THE MATRIX A.
1980 *  II  NA NUMBER OF COLUMNS OF THE MATRIX A.
1981 *  II  MA  NUMBER OF ELEMENTS IN THE FIELD A.
1982 *  RI  A(MA)  RECTANGULAR MATRIX STORED AS A TWO-DIMENSIONAL ARRAY.
1983 *  II  IA(NA+1)  POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD A.
1984 *  II  JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD A.
1985 *  RI  X(N)  INPUT VECTOR.
1986 *  RI  ALF  SCALING FACTOR.
1987 *  RI  Y(NA)  INPUT VECTOR.
1988 *  RO  Z(NA)  OUTPUT VECTOR EQUAL TO TRANS(A)*X+ALF*Y.
1989 *
1990       SUBROUTINE MXSRMD(NA,A,IA,JA,X,ALF,Y,Z)
1991       INTEGER NA,IA(*),JA(*)
1992       DOUBLE PRECISION A(*),X(*),ALF,Y(*),Z(*)
1993       DOUBLE PRECISION TEMP
1994       INTEGER I,J,K,L,JP
1995       DO 2 I=1,NA
1996       K=IA(I)
1997       L=IA(I+1)-K
1998       TEMP=ALF*Y(I)
1999       DO 1 J=1,L
2000       JP=JA(K)
2001       IF (JP.GT.0) TEMP=TEMP+A(K)*X(JP)
2002       K=K+1
2003     1 CONTINUE
2004       Z(I)=TEMP
2005     2 CONTINUE
2006       RETURN
2007       END
2008 * SUBROUTINE MXSRMM               ALL SYSTEMS                92/12/01
2009 * PURPOSE :
2010 * MULTIPLICATION OF TRANSPOSE OF A DENSE RECTANGULAR MATRIX A BY
2011 * A VECTOR X.
2012 *
2013 * PARAMETERS :
2014 *  II  N  NUMBER OF ROWS OF THE MATRIX A.
2015 *  II  NA NUMBER OF COLUMNS OF THE MATRIX A.
2016 *  II  MA  NUMBER OF ELEMENTS IN THE FIELD A.
2017 *  RI  A(MA)  RECTANGULAR MATRIX STORED AS A TWO-DIMENSIONAL ARRAY.
2018 *  II  IA(NA+1)  POSITION OF THE FIRST RORWS ELEMENTS IN THE FIELD A.
2019 *  II  JA(MA) COLUMN INDICES OF ELEMENTS IN THE FIELD A.
2020 *  RI  X(N)  INPUT VECTOR.
2021 *  RO  Y(M)  OUTPUT VECTOR EQUAL TO TRANS(A)*X.
2022 *
2023       SUBROUTINE MXSRMM(NA,A,IA,JA,X,Y)
2024       INTEGER NA,IA(*),JA(*)
2025       DOUBLE PRECISION A(*),X(*),Y(*)
2026       DOUBLE PRECISION TEMP
2027       INTEGER I,J,K,L,JP
2028       DO 2 I=1,NA
2029       K=IA(I)
2030       L=IA(I+1)-K
2031       TEMP=0.0D 0
2032       DO 1 J=1,L
2033       JP=JA(K)
2034       IF (JP.GT.0) TEMP=TEMP+A(K)*X(JP)
2035       K=K+1
2036     1 CONTINUE
2037       Y(I)=TEMP
2038     2 CONTINUE
2039       RETURN
2040       END
2041 * SUBROUTINE  MXSRSP               ALL SYSTEMS               95/12/01
2042 * PURPOSE : CREATE ROW PERMUTATIONS FOR OBTAINING DIAGONAL NONZEROS.
2043 *
2044 *  PARAMETERS :
2045 *  II  N  NUMBER OF COLUMNS OF THE MATRIX.
2046 *  II  M  NUMBER OF NONZEROS MEMBERS IN THE MATRIX.
2047 *  II  IA(M+1)  ROW POINTERS OF THE SPARSE MATRIX.
2048 *  II  JA(M)  COLUMN INDICES OF THE SPARSE MATRIX.
2049 *  IO  IP(N)  PERMUTATION VECTOR.
2050 *  II  NR  NUMBER OF STRUCTURALLY INDEPENDENT ROWS.
2051 *  IA  IW1(N) AMXILIARY VECTOR.
2052 *  IA  IW2(N) AMXILIARY VECTOR.
2053 *  IA  IW3(N) AMXILIARY VECTOR.
2054 *  IA  IW4(N) AMXILIARY VECTOR.
2055 *
2056       SUBROUTINE MXSRSP(N,IA,JA,IP,NR,IW1,IW2,IW3,IW4)
2057       INTEGER          N,NR
2058       INTEGER          IA(*),IP(*),IW1(*),IW2(*),IW3(*),IW4(*),JA(*)
2059       INTEGER          I,I1,I2,II,J,J1,K,KK,L
2060       DO 10 I = 1,N
2061           IW2(I) = IA(I+1) - IA(I) - 1
2062           IW3(I) = 0
2063           IP(I) = 0
2064    10 CONTINUE
2065       NR = 0
2066 *
2067 *     MAIN LOOP.
2068 *     EACH PASS ROUND THIS LOOP EITHER RESULTS IN A NEW ASSIGNMENT
2069 *     OR GIVES A ROW WITH NO ASSIGNMENT.
2070 *
2071       DO 100 L = 1,N
2072           J = L
2073           IW1(J) = -1
2074           DO 70 K = 1,L
2075 *
2076 *     LOOK FOR A CHEAP ASSIGNMENT
2077 *
2078               I1 = IW2(J)
2079               IF (I1.LT.0) GO TO 30
2080               I2 = IA(J+1) - 1
2081               I1 = I2 - I1
2082               DO 20 II = I1,I2
2083                   I = JA(II)
2084                   IF (IP(I).EQ.0) GO TO 80
2085    20         CONTINUE
2086 *
2087 *     NO CHEAP ASSIGNMENT IN ROW.
2088 *
2089               IW2(J) = -1
2090 *
2091 *     BEGIN LOOKING FOR ASSIGNMENT CHAIN STARTING WITH ROW J.
2092 *
2093    30         CONTINUE
2094               IW4(J) = IA(J+1) - IA(J) - 1
2095 *
2096 *     INNER LOOP.  EXTENDS CHAIN BY ONE OR BACKTRACKS.
2097 *
2098               DO 60 KK = 1,L
2099                   I1 = IW4(J)
2100                   IF (I1.LT.0) GO TO 50
2101                   I2 = IA(J+1) - 1
2102                   I1 = I2 - I1
2103 *
2104 *     FORWARD SCAN.
2105 *
2106                   DO 40 II = I1,I2
2107                       I = JA(II)
2108                       IF (IW3(I).EQ.L) GO TO 40
2109 *
2110 *     COLUMN I HAS NOT YET BEEN ACCESSED DURING THIS PASS.
2111 *
2112                       J1 = J
2113                       J = IP(I)
2114                       IW3(I) = L
2115                       IW1(J) = J1
2116                       IW4(J1) = I2 - II - 1
2117                       GO TO 70
2118    40             CONTINUE
2119 *
2120 *     BACKTRACKING STEP.
2121 *
2122    50             CONTINUE
2123                   J = IW1(J)
2124                   IF (J.EQ.-1) GO TO 100
2125    60         CONTINUE
2126    70     CONTINUE
2127 *
2128 *     NEW ASSIGNMENT IS MADE.
2129 *
2130    80     CONTINUE
2131           IP(I) = J
2132           IW2(J) = I2 - II - 1
2133           NR = NR + 1
2134           DO 90 K = 1,L
2135               J = IW1(J)
2136               IF (J.EQ.-1) GO TO 100
2137               II = IA(J+1) - IW4(J) - 2
2138               I = JA(II)
2139               IP(I) = J
2140    90     CONTINUE
2141   100 CONTINUE
2142 *
2143 *     IF MATRIX IS STRUCTURALLY SINGULAR, WE NOW COMPLETE THE
2144 *     PERMUTATION IP.
2145 *
2146       IF (NR.EQ.N) RETURN
2147       DO 110 I = 1,N
2148           IW2(I) = 0
2149   110 CONTINUE
2150       K = 0
2151       DO 130 I = 1,N
2152           IF (IP(I).NE.0) GO TO 120
2153           K = K + 1
2154           IW4(K) = I
2155           GO TO 130
2156   120     CONTINUE
2157           J = IP(I)
2158           IW2(J) = I
2159   130 CONTINUE
2160       K = 0
2161       DO 140 I = 1,N
2162           IF (IW2(I).NE.0) GO TO 140
2163           K = K + 1
2164           L = IW4(K)
2165           IP(L) = I
2166   140 CONTINUE
2167       RETURN
2168       END
2169 * SUBROUTINE MXSSDA                ALL SYSTEMS                91/12/01
2170 * PURPOSE :
2171 * A SPARSE SYMMETRIC MATRIX A IS AUGMENTED BY THE SCALED UNIT MATRIX
2172 * SUCH THAT A:=A+ALF*I (I IS THE UNIT MATRIX OF ORDER N).
2173 *
2174 * PARAMETERS :
2175 *  II  N  ORDER OF THE MATRIX A.
2176 *  RI  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
2177 *      PACKED FORM.
2178 *  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
2179 *  RI  ALF  SCALING FACTOR.
2180 *
2181       SUBROUTINE MXSSDA(N,A,IA,ALF)
2182       INTEGER N,IA(*)
2183       DOUBLE PRECISION  A(*), ALF
2184       INTEGER I
2185       DO 100 I=1,N
2186       A(IA(I))=A(IA(I))+ALF
2187 100   CONTINUE
2188       RETURN
2189       END
2190 * FUNCTION MXSSDL                  ALL SYSTEMS                88/12/01
2191 * PURPOSE :
2192 * DETERMINATION OF A MINIMUM DIAGONAL ELEMENT OF A SPARSE SYMMETRIC
2193 * MATRIX
2194 *
2195 * PARAMETERS :
2196 *  II  N  ORDER OF THE MATRIX A
2197 *  RI  A(MMAX) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
2198 *      USUAL FORM.
2199 *  II  IA(N+1)  POINTER VECTOR OF THE DIAGONAL OF THE SPARSE MATRIX.
2200 *  II  JA(MMAX)  INDICES OF NONZERO ELEMENTS OF THE SPARSE MATRIX.
2201 *  IO  INF  INDEX OF INIMUM DIAGONAL ELEMENT OF THE MATRIX A.
2202 *  RR  MXSSDL  MINIMUM DIAGONAL ELEMENT OF THE MATRIX A.
2203 *
2204       FUNCTION MXSSDL(N,A,IA,JA,INF)
2205       INTEGER  N,IA(*),JA(*),INF
2206       DOUBLE PRECISION  A(*),MXSSDL
2207       DOUBLE PRECISION  CON
2208       PARAMETER (CON=1.0D 60)
2209       INTEGER I,J
2210       INF=0
2211       MXSSDL = CON
2212       DO 100 I=1,N
2213       J=IA(I)
2214       IF (JA(J).GT.0.AND.MXSSDL.GT.A(J)) THEN
2215       INF=I
2216       MXSSDL=A(J)
2217       END IF
2218 100   CONTINUE
2219       RETURN
2220       END
2221 * SUBROUTINE MXSSMD                ALL SYSTEMS                93/12/01
2222 * PURPOSE :
2223 * MULTIPLICATION OF A SPARSE SYMMETRIC MATRIX A BY A VECTOR X
2224 * AND ADDITION OF A SCALED VECTOR ALF*Y.
2225 *
2226 * PARAMETERS :
2227 *  II  N  ORDER OF THE MATRIX A.
2228 *  RI  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
2229 *      PACKED FORM.
2230 *  II  IA(N)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
2231 *  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
2232 *  RI  X(N)  INPUT VECTOR.
2233 *  RI  ALF  SCALING FACTOR.
2234 *  RI  Y(NA)  INPUT VECTOR.
2235 *  RO  Z(NA)  OUTPUT VECTOR EQUAL TO A*X+ALF*Y.
2236 *
2237       SUBROUTINE MXSSMD(N,A,IA,JA,X,ALF,Y,Z)
2238       INTEGER N,IA(*),JA(*)
2239       DOUBLE PRECISION  A(*),X(*),ALF,Y(*),Z(*)
2240       INTEGER I,J,K,JSTRT,JSTOP
2241       DO 100 I=1,N
2242       Z(I)=ALF*Y(I)
2243 100   CONTINUE
2244       JSTOP=0
2245       DO 300 I=1,N
2246         JSTRT=JSTOP+1
2247         JSTOP=IA(I+1)-1
2248         IF (JA(JSTRT).GT.0) THEN
2249           DO 200 J=JSTRT,JSTOP
2250             K=JA(J)
2251             IF (J.EQ.JSTRT) THEN
2252               Z(I)=Z(I)+A(J)*X(I)
2253             ELSE IF (K.GT.0) THEN
2254               Z(K)=Z(K)+A(J)*X(I)
2255               Z(I)=Z(I)+A(J)*X(K)
2256             END IF
2257 200       CONTINUE
2258         END IF
2259 300   CONTINUE
2260       RETURN
2261       END
2262 * SUBROUTINE MXSSMG                ALL SYSTEMS                91/12/01
2263 * PURPOSE :
2264 *  GERSHGORIN BOUNDS OF THE EIGENVALUAE OF A DENSE SYMMETRIC MATRIX.
2265 *  AMIN .LE. ANY EIGENVALUE OF A .LE. AMAX.
2266 *
2267 * PARAMETERS :
2268 *  II  N  DIMENSION OF THE MATRIX A.
2269 *  RI  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
2270 *      PACKED FORM.
2271 *  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
2272 *  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
2273 *  RO  AMIN  LOWER BOUND OF THE EIGENVALUE OF A.
2274 *  RO  AMAX  UPPER BOUND OF THE EIGENVALUE OF A.
2275 *
2276       SUBROUTINE MXSSMG(N,A,IA,JA,AMIN,AMAX,X)
2277       INTEGER N,IA(*),JA(*)
2278       DOUBLE PRECISION  A(*),AMIN,AMAX,X(*)
2279       INTEGER I,J,K,JSTRT,JSTOP
2280       DOUBLE PRECISION CMAX
2281       PARAMETER (CMAX=1.0D 60)
2282       DO 1 I=1,N
2283       X(I)=0.0D 0
2284     1 CONTINUE
2285       JSTOP=0
2286       DO 3 I=1,N
2287       JSTRT=JSTOP+1
2288       JSTOP=IA(I+1)-1
2289       IF (JA(JSTRT).GT.0) THEN
2290       DO 2 K=JSTRT+1,JSTOP
2291       J=JA(K)
2292       IF (J.GT.0) THEN
2293       X(I)=X(I)+ABS(A(K))
2294       X(J)=X(J)+ABS(A(K))
2295       END IF
2296     2 CONTINUE
2297       END IF
2298     3 CONTINUE
2299       AMIN= CMAX
2300       AMAX=-CMAX
2301       DO 4 I=1,N
2302       K=IA(I)
2303       IF (K.GT.0) THEN
2304       AMAX=MAX(AMAX,A(K)+X(I))
2305       AMIN=MIN(AMIN,A(K)-X(I))
2306       END IF
2307     4 CONTINUE
2308       RETURN
2309       END
2310 * SUBROUTINE MXSSMI                ALL SYSTEMS                92/12/01
2311 * PURPOSE :
2312 * SPARSE SYMMETRIC MATRIX A IS SET TO THE UNIT MATRIX WITH THE SAME
2313 * ORDER.
2314 *
2315 * PARAMETERS :
2316 *  II  N  ORDER OF THE MATRIX A.
2317 *  RU  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
2318 *      PACKED FORM.
2319 *  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
2320 *  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
2321 *
2322       SUBROUTINE MXSSMI(N,A,IA)
2323       INTEGER N,IA(*)
2324       DOUBLE PRECISION  A(*)
2325       INTEGER I,K
2326       DO 100 I=1,IA(N+1)-1
2327       A(I)=0.0D 0
2328 100   CONTINUE
2329       DO 200 I=1,N
2330       K=ABS(IA(I))
2331       A(K)=1.0D 0
2332 200   CONTINUE
2333       RETURN
2334       END
2335 * SUBROUTINE MXSSMM                ALL SYSTEMS                92/12/01
2336 * PURPOSE :
2337 * MULTIPLICATION OF A SPARSE SYMMETRIC MATRIX BY A VECTOR X.
2338 *
2339 * PARAMETERS :
2340 *  II  N  ORDER OF THE MATRIX A.
2341 *  II  M  NUMBER OF NONZERO ELEMENTS OF THE MATRIX A.
2342 *  RI  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
2343 *      PACKED FORM.
2344 *  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
2345 *  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
2346 *  RI  X(N)  INPUT VECTOR.
2347 *  RO  Y(N)  OUTPUT VECTOR WHERE Y := A * X.
2348 *
2349       SUBROUTINE MXSSMM(N,A,IA,JA,X,Y)
2350       INTEGER N,IA(*),JA(*)
2351       DOUBLE PRECISION  A(*),X(*),Y(*)
2352       INTEGER I,J,K,JSTRT,JSTOP
2353       DO 100 I=1,N
2354       Y(I)=0.0D 0
2355 100   CONTINUE
2356       JSTOP=0
2357       DO 300 I=1,N
2358         JSTRT=JSTOP+1
2359         JSTOP=IA(I+1)-1
2360         IF (JA(JSTRT).GT.0) THEN
2361           DO 200 J=JSTRT,JSTOP
2362             K=JA(J)
2363             IF (J.EQ.JSTRT) THEN
2364               Y(I)=Y(I)+A(J)*X(I)
2365             ELSE IF (K.GT.0) THEN
2366               Y(K)=Y(K)+A(J)*X(I)
2367               Y(I)=Y(I)+A(J)*X(K)
2368             END IF
2369 200       CONTINUE
2370         END IF
2371 300   CONTINUE
2372       RETURN
2373       END
2374 * SUBROUTINE MXSSMN                ALL SYSTEMS                89/12/01
2375 * PURPOSE :
2376 * THIS SUBROUTINE FINDS THE PERMUTATION VECTOR PERM FOR THE
2377 * SPARSE SYMMETRIC MATRIX GIVEN IN THE VECTOR PAIR PA,SA.IT USES
2378 * THE SO-CALLED NESTED DISSECTION METHOD.
2379 *
2380 * PARAMETERS :
2381 *  II  N ORDER OF THE MATRIX A.
2382 *  II  MMAX  LENGTH OF THE PRINCIPAL MATRIX VECTORS.
2383 *  II  PA(N+1) POINTER VECTOR OF THE INPUT MATRIX.
2384 *  II  SA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX.
2385 *  IO  PERM(N) PERMUTATION VECTOR.
2386 *  IA  WN11(N+1) AMXILIARY VECTOR.
2387 *  IA  WN12(N+1) AMXILIARY VECTOR.
2388 *  IA  WN13(N+1) AMXILIARY VECTOR.
2389 *
2390 * METHOD :
2391 * NESTED DISSECTION METHOD
2392 *
2393       SUBROUTINE MXSSMN(N,PA,SA,PERM,WN11,WN12,WN13)
2394       INTEGER N
2395       INTEGER PA(*),SA(*),PERM(*)
2396       INTEGER WN11(*),WN12(*),WN13(*)
2397       INTEGER I,J,K,NUM,ROOT,NLVL,LVLEND,LBEGIN,ICS
2398       INTEGER NN,N1,MINDEG,N2,LVSIZE,NDEG,NUNLVL,MIDLVL
2399       INTEGER TEMP,NPUL,NSEP,I1,I2,I3,I4,J1,J2
2400       DO 100 I = 1, N
2401       WN11(I) = 1
2402   100 CONTINUE
2403       NUM= 0
2404       DO 2000 I = 1, N
2405   200 IF ( WN11(I) .EQ. 0 ) GO TO 2000
2406       ROOT = I
2407       WN11(ROOT) = 0
2408       WN13(1) = ROOT
2409       NLVL = 0
2410       LVLEND = 0
2411       ICS = 1
2412   300 LBEGIN = LVLEND + 1
2413       LVLEND = ICS
2414       NLVL = NLVL + 1
2415       WN12(NLVL) = LBEGIN
2416       DO 500 K = LBEGIN, LVLEND
2417       NN = WN13(K)
2418       DO 400 J=PA(NN),PA(NN+1)-1
2419       N2 = SA(J)
2420       IF (N2.EQ.NN) GO TO 400
2421       IF  (WN11(N2).EQ.0) GO TO 400
2422       ICS = ICS + 1
2423       WN13(ICS) = N2
2424       WN11(N2) = 0
2425   400 CONTINUE
2426   500 CONTINUE
2427       LVSIZE=ICS-LVLEND
2428       IF (LVSIZE.GT.0) GO TO 300
2429       WN12(NLVL+1) = LVLEND + 1
2430       DO 600 K = 1, ICS
2431       NN = WN13(K)
2432       WN11(NN) = 1
2433   600 CONTINUE
2434       ICS = WN12(NLVL+1) - 1
2435       IF  ( NLVL .EQ. 1 .OR. NLVL .EQ. ICS )GO TO 1470
2436   700 J1 = WN12(NLVL)
2437       MINDEG = ICS
2438       ROOT = WN13(J1)
2439       IF  ( ICS .EQ. J1 ) GO TO 1000
2440       DO 900 J = J1, ICS
2441       NN = WN13(J)
2442       NDEG = 0
2443       DO 800 K=PA(NN),PA(NN+1)-1
2444       N1 = SA(K)
2445       IF (N1.EQ.NN) GO TO 800
2446       IF  ( WN11(N1)  .GT. 0 ) NDEG = NDEG + 1
2447   800 CONTINUE
2448       IF  ( NDEG .GE. MINDEG  ) GO TO 900
2449       ROOT = NN
2450       MINDEG = NDEG
2451   900 CONTINUE
2452  1000 CONTINUE
2453       WN11(ROOT) = 0
2454       WN13(1) = ROOT
2455       NUNLVL = 0
2456       LVLEND = 0
2457       ICS = 1
2458  1100 LBEGIN = LVLEND + 1
2459       LVLEND = ICS
2460       NUNLVL = NUNLVL + 1
2461       WN12(NUNLVL) = LBEGIN
2462       DO 1300 K = LBEGIN, LVLEND
2463       NN = WN13(K)
2464       DO 1200 J=PA(NN),PA(NN+1)-1
2465       N2 = SA(J)
2466       IF (N2.EQ.NN) GO TO 1200
2467       IF  (WN11(N2).EQ.0) GO TO 1200
2468       ICS = ICS + 1
2469       WN13(ICS) = N2
2470       WN11(N2) = 0
2471  1200 CONTINUE
2472  1300 CONTINUE
2473       LVSIZE=ICS-LVLEND
2474       IF (LVSIZE.GT.0) GO TO 1100
2475       WN12(NUNLVL+1) = LVLEND + 1
2476       DO 1400 K = 1, ICS
2477       NN = WN13(K)
2478       WN11(NN) = 1
2479  1400 CONTINUE
2480       IF (NUNLVL .LE.NLVL) GO TO 1470
2481       NLVL=NUNLVL
2482       IF (NLVL.LT.ICS) GO TO 700
2483 1470  CONTINUE
2484       IF ( NLVL .GE. 3 ) GO TO 1600
2485       NSEP = WN12(NLVL+1) - 1
2486       DO 1500 K = 1, NSEP
2487       NN = WN13(K)
2488       PERM(NUM+K) = NN
2489       WN11(NN) = 0
2490  1500 CONTINUE
2491       GO TO 1950
2492  1600 MIDLVL = (NLVL+2)/2
2493       I3 = WN12(MIDLVL)
2494       I1 = WN12(MIDLVL + 1)
2495       I4 = I1 - 1
2496       I2 = WN12(MIDLVL+2) - 1
2497       DO 1700 K = I1, I2
2498       NN = WN13(K)
2499       PA(NN) = - PA(NN)
2500  1700 CONTINUE
2501       NSEP = 0
2502       DO 1800 K = I3, I4
2503       NN = WN13(K)
2504       J1 = PA(NN)
2505       J2 = IABS(PA(NN+1)) - 1
2506       DO 1750 J = J1, J2
2507       N2 = SA(J)
2508       IF (N2.EQ.NN) GO TO 1750
2509       IF ( PA(N2) .GT. 0 ) GO TO 1750
2510       NSEP = NSEP + 1
2511       PERM(NSEP+NUM) = NN
2512       WN11(NN) = 0
2513       GO TO 1800
2514  1750 CONTINUE
2515  1800 CONTINUE
2516       DO 1900 K = I1, I2
2517       NN = WN13(K)
2518       PA(NN) = - PA(NN)
2519  1900 CONTINUE
2520  1950 CONTINUE
2521       NUM = NUM + NSEP
2522       IF ( NUM .GE. N ) GO TO 2100
2523       GO TO 200
2524  2000 CONTINUE
2525  2100 CONTINUE
2526       IF (N.LT.2) GO TO 2300
2527       NPUL = N/2
2528       DO 2200 I=1,NPUL
2529       TEMP=PERM(I)
2530       PERM(I)=PERM(N-I+1)
2531       PERM(N-I+1)=TEMP
2532 2200  CONTINUE
2533 2300  CONTINUE
2534       RETURN
2535       END
2536 * FUNCTION MXSSMQ                  ALL SYSTEMS                92/12/01
2537 * PURPOSE :
2538 * VALUE OF A QUADRATIC FORM WITH A SPARSE SYMMETRIC MATRIX A.
2539 *
2540 * PARAMETERS :
2541 *  II  N  ORDER OF THE MATRIX A.
2542 *  RI  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
2543 *      PACKED FORM.
2544 *  II  IA(N+1)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
2545 *  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
2546 *  RI  X(N)  INPUT VECTOR.
2547 *  RI  Y(N)  INPUT VECTOR.
2548 *  RR  MXSSMQ  VALUE OF THE QUADRATIC FORM MXSSMQ=TRANS(Y)*A*X.
2549 *
2550       FUNCTION MXSSMQ(N,A,IA,JA,X,Y)
2551       INTEGER  N,IA(*),JA(*)
2552       DOUBLE PRECISION  A(*), X(*), Y(*), MXSSMQ
2553       DOUBLE PRECISION  TEMP1,TEMP2
2554       INTEGER  I,J,K,JSTRT,JSTOP
2555       JSTOP=0
2556       TEMP1=0.0D 0
2557       DO 300 I=1,N
2558         JSTRT=JSTOP+1
2559         JSTOP=IA(I+1)-1
2560         IF (JA(JSTRT).GT.0) THEN
2561           TEMP2=0.0D 0
2562           DO 200 J=JSTRT,JSTOP
2563             K=JA(J)
2564             IF (J.EQ.JSTRT) THEN
2565               TEMP2=TEMP2+A(J)*Y(I)
2566             ELSE IF (K.GT.0) THEN
2567               TEMP2=TEMP2+2*Y(K)*A(J)
2568           END IF
2569 200       CONTINUE
2570           TEMP1=TEMP1+X(I)*TEMP2
2571         END IF
2572 300   CONTINUE
2573       MXSSMQ=TEMP1
2574       RETURN
2575       END
2576 * SUBROUTINE MXSSMY                ALL SYSTEMS                93/12/01
2577 * PURPOSE :
2578 * CORRECTION OF A SPARSE SYMMETRIC MATRIX A. THE CORRECTION IS DEFINED
2579 * AS A:=A+SUM OF (HALF*(X*TRANS(Y)+Y*TRANS(X)))(I)/SIGMA(I) WHERE
2580 * SIGMA(I) IS A DOT PRODUCT TRANS(X)*X WHERE ONLY CONTRIBUTIONS
2581 * CORRESPONDING TO NONZEROS IN ROW I ARE SUMMED UP, X AND Y ARE GIVEN
2582 * VECTORS.
2583 *
2584 * PARAMETERS :
2585 *  II  N  ORDER OF THE MATRIX A.
2586 *  RI  A(M) ELEMENTS OF THE SPARSE SYMMETRIC MATRIX STORED IN THE
2587 *      PACKED FORM.
2588 *  II  IA(N)  POINTERS OF THE DIAGONAL ELEMENTS OF A.
2589 *  II  JA(M)  INDICES OF THE NONZERO ELEMENTS OF A.
2590 *  RA  XS(N) AMXILIARY VECTOR - USED FOR SIGMA(I).
2591 *  RI  X(N)  VECTOR IN THE CORRECTION TERM.
2592 *  RI  Y(N)  VECTOR IN THE CORRECTION TERM.
2593 *
2594       SUBROUTINE MXSSMY(N,A,IA,JA,XS,X,Y)
2595       INTEGER  N,IA(*),JA(*)
2596       DOUBLE PRECISION  A(*),X(*),Y(*),XS(*),SIGMA,TEMP
2597       INTEGER I,J,K,JSTRT,JSTOP
2598       CALL MXVSET(N,0.0D 0,XS)
2599 *
2600 *      COMPUTE SIGMA(I)
2601 *
2602       JSTOP=0
2603       DO 200 I=1,N
2604         JSTRT=JSTOP+1
2605         JSTOP=IA(I+1)-1
2606         IF (JA(JSTRT).GT.0) THEN
2607           SIGMA=0.0D 0
2608           DO 100 J=JSTRT,JSTOP
2609             K=JA(J)
2610             IF (K.GT.0) THEN
2611               SIGMA=SIGMA+Y(K)*Y(K)
2612               IF (K.NE.I) XS(K)=XS(K)+Y(I)*Y(I)
2613             END IF
2614 100       CONTINUE
2615           XS(I)=XS(I)+SIGMA
2616         END IF
2617 200   CONTINUE
2618 *
2619 *      UPDATE MATRIX
2620 *
2621       JSTOP=0
2622       DO 400 I=1,N
2623         JSTRT=JSTOP+1
2624         JSTOP=IA(I+1)-1
2625         IF (JA(JSTRT).GT.0) THEN
2626           IF (XS(I).EQ.0.0D 0) THEN
2627             TEMP=0.0D 0
2628           ELSE
2629             TEMP=X(I)/XS(I)
2630           END IF
2631           DO 300 J=JSTRT,JSTOP
2632             K=JA(J)
2633             IF (K.GT.0) THEN
2634               IF (XS(K).EQ.0.0D 0) THEN
2635               A(J)=A(J)+0.5D 0*TEMP*Y(K)
2636               ELSE
2637               A(J)=A(J)+0.5D 0*(TEMP*Y(K)+Y(I)*X(K)/XS(K))
2638               END IF
2639             END IF
2640 300       CONTINUE
2641         END IF
2642 400   CONTINUE
2643       RETURN
2644       END
2645 * SUBROUTINE MXSTG1                ALL SYSTEMS                89/12/01
2646 * PURPOSE :
2647 * WIDTHENING THE PACKED FORM OF THE VECTORS IA, JA OF THE SPARSE MATRIX
2648 *
2649 * PARAMETERS :
2650 *  II  N ORDER OF THE SPARSE MATRIX.
2651 *  IU  M NUMBER OF NONZERO ELEMENTS IN THE MATRIX.
2652 *  II  MMAX LENGTH OF THE ARRAY JA.
2653 *  II  IA(N+1) POINTER VECTOR OF THE INPUT MATRIX.
2654 *  II  JA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX.
2655 *  IA  PD(N+1) AMXILIARY VECTOR.
2656 *  IA  WN11(N+1) AMXILIARY VECTOR.
2657 *
2658       SUBROUTINE MXSTG1(N,M,IA,JA,PD,WN11)
2659       INTEGER N,M
2660       INTEGER IA(*),PD(*),JA(*),WN11(*)
2661       INTEGER I,J,L1,L,K
2662 *
2663 *     UPPER TRIANGULAR INFORMATION TO THE AMXILIARY ARRAY
2664 *
2665       L1=IA(1)
2666       DO 100 I=1,N
2667       L=L1
2668       L1=IA(I+1)
2669       WN11(I)=L1-L
2670 100   CONTINUE
2671 *
2672 *     LOWER TRIANGULAR INFORMATION TO THE AMXILIARY ARRAY
2673 *
2674       DO 300  I=1,N
2675         DO 200  J=IA(I)+1,IA(I+1)-1
2676         K=ABS(JA(J))
2677         WN11(K)=WN11(K)+1
2678 200     CONTINUE
2679 300   CONTINUE
2680 *
2681 *     BY PARTIAL SUMMING WE GET POINTERS OF THE WIDE STRUCTURE
2682 *     WN11(I) POINTS AT THE END OF THE ROW I
2683 *
2684       L=0
2685       DO 400  I=2,N
2686       WN11(I)=WN11(I)+WN11(I-1)
2687 400   CONTINUE
2688 *
2689 *     DEFINE LENGTH OF THE WITHENED STRUCTURE
2690 *
2691       M=WN11(N)
2692 *
2693 *     SHIFT OF UPPER TRIANGULAR ROWS
2694 *
2695       PD(1)=1
2696       DO 600  I=N,1,-1
2697       L=WN11(I)
2698       PD(I+1)=L+1
2699         DO 500  J=IA(I+1)-1,IA(I),-1
2700         JA(L)=JA(J)
2701         L=L-1
2702 500     CONTINUE
2703 600   CONTINUE
2704 *
2705 *     FORMING OF THE LOWER TRIANGULAR PART
2706 *
2707       DO 800  I=1,N
2708         DO 700  J=WN11(I)+IA(I)+2-IA(I+1),WN11(I)
2709         K=ABS(JA(J))
2710         JA(PD(K))=I
2711         PD(K)=PD(K)+1
2712 700     CONTINUE
2713 800   CONTINUE
2714       DO 900  I=1,N
2715       IA(I+1)=WN11(I)+1
2716 900   CONTINUE
2717       RETURN
2718       END
2719 * SUBROUTINE MXSTL1                ALL SYSTEMS                91/12/01
2720 * PURPOSE :
2721 * PACKING OF THE WIDTHENED FORM OF THE VECTORS IA, JA OF THE SPARSE
2722 * MATRIX
2723 *
2724 * PARAMETERS :
2725 *  II  N ORDER OF THE SPARSE MATRIX.
2726 *  IU  M NUMBER OF NONZERO ELEMENTS IN THE MATRIX.
2727 *  II  MMAX LENGTH OF THE ARRAY JA.
2728 *  IU  IA(N+1) POINTER VECTOR OF THE INPUT MATRIX.
2729 *  IU  JA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX.
2730 *  IA  PD(N+1) AMXILIARY VECTOR.
2731 *
2732       SUBROUTINE MXSTL1(N,M,IA,JA,PD)
2733       INTEGER N,M
2734       INTEGER IA(*),PD(*),JA(*)
2735       INTEGER I,J,L,JSTRT,JSTOP
2736       L=1
2737 *
2738 *     PD DEFINITION
2739 *
2740       JSTOP=0
2741       DO 60 I=1,N
2742         JSTRT=JSTOP+1
2743         JSTOP=IA(I+1)-1
2744         DO 50 J=JSTRT,JSTOP
2745           IF (ABS(JA(J)).EQ.I) THEN
2746             PD(I)=J
2747             GO TO 60
2748           END IF
2749    50    CONTINUE
2750    60 CONTINUE
2751 *
2752 *     REWRITE THE STRUCTURE
2753 *
2754       DO 200 I=1,N
2755         DO 100 J=PD(I),IA(I+1)-1
2756         JA(L)=JA(J)
2757         L=L+1
2758 100     CONTINUE
2759       IA(I+1)=L
2760 200   CONTINUE
2761       IA(1)=1
2762 *
2763 *     SET THE LENGTH OF THE PACKED STRUCTURE
2764 *
2765       M=L-1
2766       RETURN
2767       END
2768 * SUBROUTINE MXSTL2                ALL SYSTEMS                90/12/01
2769 * PURPOSE :
2770 * PACKING OF THE WIDTHENED FORM OF THE VECTORS A,IA,JA OF THE SPARSE
2771 * MATRIX
2772 *
2773 * PARAMETERS :
2774 *  II  N ORDER OF THE SPARSE MATRIX.
2775 *  IU  M NUMBER OF NONZERO ELEMENTS IN THE MATRIX.
2776 *  II  MMAX LENGTH OF THE ARRAY JA.
2777 *  RU  A(MMAX) VECTOR OF NUMERICAL VALUES OF THE MATRIX BEING SHRINKED.
2778 *  IU  IA(N+1) POINTER VECTOR OF THE INPUT MATRIX.
2779 *  IU  JA(MMAX) VECTOR OF THE COLUMN INDICES OF THE INPUT MATRIX.
2780 *  IA  PD(N+1) AMXILIARY VECTOR.
2781 *
2782       SUBROUTINE MXSTL2(N,M,A,IA,JA,PD)
2783       INTEGER N,M
2784       INTEGER IA(*),PD(*),JA(*)
2785       DOUBLE PRECISION A(*)
2786       INTEGER I,J,L,JSTRT,JSTOP
2787       L=1
2788 *
2789 *     PD DEFINITION
2790 *
2791       JSTOP=0
2792       DO 60 I=1,N
2793         JSTRT=JSTOP+1
2794         JSTOP=IA(I+1)-1
2795         DO 50 J=JSTRT,JSTOP
2796           IF (ABS(JA(J)).EQ.I) THEN
2797             PD(I)=J
2798             GO TO 60
2799           END IF
2800    50    CONTINUE
2801    60 CONTINUE
2802 *
2803 *     REWRITE THE STRUCTURE
2804 *
2805       DO 200 I=1,N
2806         DO 100 J=PD(I),IA(I+1)-1
2807         JA(L)=JA(J)
2808         A(L)=A(J)
2809         L=L+1
2810 100     CONTINUE
2811       IA(I+1)=L
2812 200   CONTINUE
2813       IA(1)=1
2814 *
2815 *     SET THE LENGTH OF THE PACKED STRUCTURE
2816 *
2817       M=L-1
2818       RETURN
2819       END
2820 * SUBROUTINE MXTPGB                ALL SYSTEMS                93/12/01
2821 * PURPOSE :
2822 * BACK SUBSTITUTION FOR A DECOMPOSED TRIDIAGONAL MATRIX.
2823 *
2824 * PARAMETERS :
2825 *  II  N  ORDER OF THE TRIDIAGONAL MATRIX T.
2826 *  RI  D(N)  ELEMENTS OF THE DIAGONAL MATRIX D IN THE DECOMPOSITION
2827 *         T=L*D*TRANS(L).
2828 *  RI  E(N)  SUBDIAGONAL ELEMENTS OF THE LOWER TRIANGULAR MATRIX L IN
2829 *         THE DECOMPOSITION T=L*D*TRANS(L).
2830 *  RU  X(N)  ON INPUT THE RIGHT HAND SIDE OF A SYSTEM OF LINEAR
2831 *         EQUATIONS. ON OUTPUT THE SOLUTION OF A SYSTEM OF LINEAR
2832 *         EQUATIONS.
2833 *  II  JOB  OPTION. IF JOB=0 THEN X:=T**(-1)*X. IF JOB>0 THEN
2834 *         X:=L**(-1)*X. IF JOB<0 THEN X:=TRANS(L)**(-1)*X.
2835 *
2836       SUBROUTINE MXTPGB(N,D,E,X,JOB)
2837       INTEGER N,JOB
2838       DOUBLE PRECISION  D(*),E(*),X(*)
2839       INTEGER I
2840       IF (JOB.GE.0) THEN
2841 *
2842 *     PHASE 1 : X:=L**(-1)*X
2843 *
2844       DO 1 I=2,N
2845       X(I)=X(I)-X(I-1)*E(I-1)
2846     1 CONTINUE
2847       END IF
2848       IF (JOB.EQ.0) THEN
2849 *
2850 *     PHASE 2 : X:=D**(-1)*X
2851 *
2852       DO 2 I=1,N
2853       X(I)=X(I)/D(I)
2854     2 CONTINUE
2855       END IF
2856       IF (JOB.LE.0) THEN
2857 *
2858 *     PHASE 3 : X:=TRANS(L)**(-1)*X
2859 *
2860       DO 3 I=N-1,1,-1
2861       X(I)=X(I)-X(I+1)*E(I)
2862     3 CONTINUE
2863       END IF
2864       RETURN
2865       END
2866 * SUBROUTINE MXTPGF                ALL SYSTEMS                03/12/01
2867 * PURPOSE :
2868 * CHOLESKI DECOMPOSITION OF A TRIDIAGONAL MATRIX.
2869 *
2870 * PARAMETERS :
2871 *  II  N  ORDER OF THE TRIDIAGONAL MATRIX T.
2872 *  RU  D(N)  ON INPUT DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX T.
2873 *         ON OUTPUT ELEMENTS OF THE DIAGONAL MATRIX D IN THE
2874 *         DECOMPOSITION T=L*D*TRANS(L).
2875 *  RU  E(N)  ON INPUT SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX T.
2876 *         ON OUTPUT SUBDIAGONAL ELEMENTS OF THE LOWER TRIANGULAR MATRIX L
2877 *         IN THE DECOMPOSITION T=L*D*TRANS(L).
2878 *  IO  INF  AN INFORMATION OBTAINED IN THE FACTORIZATION PROCESS. IF
2879 *         INF=0 THEN A IS SUFFICIENTLY POSITIVE DEFINITE AND E=0. IF
2880 *         INF<0 THEN A IS NOT SUFFICIENTLY POSITIVE DEFINITE AND E>0. IF
2881 *         INF>0 THEN A IS INDEFINITE AND INF IS AN INDEX OF THE
2882 *         MOST NEGATIVE DIAGONAL ELEMENT USED IN THE FACTORIZATION
2883 *         PROCESS.
2884 *  RU  ALF  ON INPUT A DESIRED TOLERANCE FOR POSITIVE DEFINITENESS. ON
2885 *         OUTPUT THE MOST NEGATIVE DIAGONAL ELEMENT USED IN THE
2886 *         FACTORIZATION PROCESS (IF INF>0).
2887 *  RO  TAU  MAXIMUM DIAGONAL ELEMENT OF THE MATRIX E.
2888 *
2889       SUBROUTINE MXTPGF(N,D,E,INF,ALF,TAU)
2890       INTEGER N,INF
2891       DOUBLE PRECISION  D(*),E(*),ALF,TAU
2892       DOUBLE PRECISION  DI,EI,BET,GAM,DEL,TOL
2893       INTEGER I,L
2894       DOUBLE PRECISION ZERO,ONE,TWO
2895       PARAMETER (ZERO=0.0D 0,ONE=1.0D 0,TWO=2.0D 0)
2896       L=0
2897       INF=0
2898       TOL=ALF
2899 *
2900 *     ESTIMATION OF THE MATRIX NORM
2901 *
2902       ALF=ZERO
2903       GAM=ZERO
2904       TAU=ZERO
2905       BET=ABS(D(1))
2906       DO 1 I=1,N-1
2907       BET=MAX(BET,ABS(D(I+1)))
2908       GAM=MAX(GAM,ABS(E(I)))
2909     1 CONTINUE
2910       BET=MAX(TOL,TWO*BET,GAM/MAX(ONE,DBLE(N-1)))
2911       DEL=TOL*MAX(BET,ONE)
2912       DO 2 I=1,N
2913       EI=D(I)
2914       IF (ALF.GT.EI) THEN
2915       ALF=EI
2916       L=I
2917       END IF
2918       GAM=ZERO
2919       IF (I.LT.N) GAM=E(I)**2
2920       DI=MAX(ABS(EI),GAM/BET,DEL)
2921       IF (TAU.LT.DI-EI) THEN
2922       TAU=DI-EI
2923       INF=-1
2924       END IF
2925 *
2926 *     GAUSSIAN ELIMINATION
2927 *
2928       D(I)=DI
2929       IF (I.LT.N) THEN
2930       EI=E(I)
2931       E(I)=EI/DI
2932       D(I+1)=D(I+1)-E(I)*EI
2933       END IF
2934     2 CONTINUE
2935       IF (L.GT.0.AND.ABS(ALF).GT.DEL) INF = L
2936       RETURN
2937       END
2938 * SUBROUTINE MXUCOP                ALL SYSTEMS                99/12/01
2939 * PURPOSE :
2940 * COPY OF THE VECTOR WITH INITIATION OF THE ACTIVE PART.
2941 *
2942 * PARAMETERS :
2943 *  II  N  VECTOR DIMENSION.
2944 *  RI  X(N)  INPUT VECTOR.
2945 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
2946 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
2947 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
2948 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
2949 *         IX(I).EQ.-5.
2950 *
2951       SUBROUTINE MXUCOP(N,X,Y,IX,JOB)
2952       INTEGER N,IX(*),JOB
2953       DOUBLE PRECISION X(*),Y(*)
2954       INTEGER I
2955       IF (JOB.EQ.0) THEN
2956       DO 1 I=1,N
2957       Y(I)=X(I)
2958     1 CONTINUE
2959       ELSE IF (JOB.GT.0) THEN
2960       DO 2 I=1,N
2961       IF (IX(I).GE. 0) THEN
2962       Y(I)=X(I)
2963       ELSE
2964       Y(I)=0.0D 0
2965       END IF
2966     2 CONTINUE
2967       ELSE
2968       DO 3 I=1,N
2969       IF (IX(I).NE.-5) THEN
2970       Y(I)=X(I)
2971       ELSE
2972       Y(I)=0.0D 0
2973       END IF
2974     3 CONTINUE
2975       END IF
2976       RETURN
2977       END
2978 * FUNCTION MXUDEL                  ALL SYSTEMS                99/12/01
2979 * PURPOSE :
2980 *  SQUARED NORM OF A SHIFTED VECTOR IN A BOUND CONSTRAINED CASE.
2981 *
2982 * PARAMETERS :
2983 *  II  N  VECTOR DIMENSION.
2984 *  RI  A  SCALING FACTOR.
2985 *  RI  X(N)  INPUT VECTOR.
2986 *  RI  Y(N)  INPUT VECTOR.
2987 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
2988 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
2989 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
2990 *         IX(I).EQ.-5.
2991 *  RR  MXUDEL SQUARED NORM OF Y+A*X.
2992 *
2993       FUNCTION MXUDEL(N,A,X,Y,IX,JOB)
2994       INTEGER N,IX(N),JOB
2995       DOUBLE PRECISION A,X(N),Y(N),MXUDEL
2996       INTEGER I
2997       DOUBLE PRECISION TEMP
2998       TEMP=0.0D 0
2999       IF (JOB.EQ.0) THEN
3000       DO 1 I=1,N
3001       TEMP=TEMP+(Y(I)+A*X(I))**2
3002     1 CONTINUE
3003       ELSE IF (JOB.GT.0) THEN
3004       DO 2 I=1,N
3005       IF (IX(I).GE. 0) TEMP=TEMP+(Y(I)+A*X(I))**2
3006     2 CONTINUE
3007       ELSE
3008       DO 3 I=1,N
3009       IF (IX(I).NE.-5) TEMP=TEMP+(Y(I)+A*X(I))**2
3010     3 CONTINUE
3011       END IF
3012       MXUDEL=TEMP
3013       RETURN
3014       END
3015 * SUBROUTINE MXUDIF                ALL SYSTEMS                99/12/01
3016 * PURPOSE :
3017 * VECTOR DIFFERENCE IN A BOUND CONSTRAINED CASE.
3018 *
3019 * PARAMETERS :
3020 *  II  N  VECTOR DIMENSION.
3021 *  RI  X(N)  INPUT VECTOR.
3022 *  RI  Y(N)  INPUT VECTOR.
3023 *  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X - Y.
3024 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
3025 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
3026 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
3027 *         IX(I).EQ.-5.
3028 *
3029       SUBROUTINE MXUDIF(N,X,Y,Z,IX,JOB)
3030       INTEGER N,IX(N),JOB
3031       DOUBLE PRECISION X(*),Y(*),Z(*)
3032       INTEGER I
3033       IF (JOB.EQ.0) THEN
3034       DO 1 I=1,N
3035       Z(I)=X(I)-Y(I)
3036     1 CONTINUE
3037       ELSE IF (JOB.GT.0) THEN
3038       DO 2 I=1,N
3039       IF (IX(I).GE. 0) Z(I)=X(I)-Y(I)
3040     2 CONTINUE
3041       ELSE
3042       DO 3 I=1,N
3043       IF (IX(I).NE.-5) Z(I)=X(I)-Y(I)
3044     3 CONTINUE
3045       END IF
3046       RETURN
3047       END
3048 * SUBROUTINE MXUDIR                ALL SYSTEMS                99/12/01
3049 * PURPOSE :
3050 * VECTOR AUGMENTED BY THE SCALED VECTOR IN A BOUND CONSTRAINED CASE.
3051 *
3052 * PARAMETERS :
3053 *  II  N  VECTOR DIMENSION.
3054 *  RI  A  SCALING FACTOR.
3055 *  RI  X(N)  INPUT VECTOR.
3056 *  RI  Y(N)  INPUT VECTOR.
3057 *  RO  Z(N)  OUTPUT VECTOR WHERE Z:= Y + A*X.
3058 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
3059 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
3060 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
3061 *         IX(I).EQ.-5.
3062 *
3063       SUBROUTINE MXUDIR(N,A,X,Y,Z,IX,JOB)
3064       INTEGER N,IX(*),JOB
3065       DOUBLE PRECISION  A, X(*), Y(*), Z(*)
3066       INTEGER I
3067       IF (JOB.EQ.0) THEN
3068       DO 1 I=1,N
3069       Z(I) = Y(I) + A*X(I)
3070     1 CONTINUE
3071       ELSE IF (JOB.GT.0) THEN
3072       DO 2 I=1,N
3073       IF (IX(I).GE. 0) Z(I) = Y(I) + A*X(I)
3074     2 CONTINUE
3075       ELSE
3076       DO 3 I=1,N
3077       IF (IX(I).NE.-5) Z(I) = Y(I) + A*X(I)
3078     3 CONTINUE
3079       END IF
3080       RETURN
3081       END
3082 * FUNCTION MXUDOT                  ALL SYSTEMS                99/12/01
3083 * PURPOSE :
3084 * DOT PRODUCT OF VECTORS IN A BOUND CONSTRAINED CASE.
3085 *
3086 * PARAMETERS :
3087 *  II  N  VECTOR DIMENSION.
3088 *  RI  X(N)  INPUT VECTOR.
3089 *  RI  Y(N)  INPUT VECTOR.
3090 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
3091 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
3092 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
3093 *         IX(I).EQ.-5.
3094 *  RR  MXUDOT  VALUE OF DOT PRODUCT MXUDOT=TRANS(X)*Y.
3095 *
3096       FUNCTION MXUDOT(N,X,Y,IX,JOB)
3097       INTEGER N,IX(*),JOB
3098       DOUBLE PRECISION X(*),Y(*),MXUDOT
3099       DOUBLE PRECISION TEMP
3100       INTEGER I
3101       DOUBLE PRECISION ZERO
3102       PARAMETER (ZERO = 0.0D 0)
3103       TEMP = ZERO
3104       IF (JOB.EQ.0) THEN
3105       DO 1 I=1,N
3106       TEMP=TEMP+X(I)*Y(I)
3107     1 CONTINUE
3108       ELSE IF (JOB.GT.0) THEN
3109       DO 2 I=1,N
3110       IF (IX(I).GE. 0) TEMP=TEMP+X(I)*Y(I)
3111     2 CONTINUE
3112       ELSE
3113       DO 3 I=1,N
3114       IF (IX(I).NE.-5) TEMP=TEMP+X(I)*Y(I)
3115     3 CONTINUE
3116       END IF
3117       MXUDOT=TEMP
3118       RETURN
3119       END
3120 * SUBROUTINE MXUNEG                ALL SYSTEMS                00/12/01
3121 * PURPOSE :
3122 * COPY OF THE VECTOR WITH INITIATION OF THE ACTIVE PART.
3123 *
3124 * PARAMETERS :
3125 *  II  N  VECTOR DIMENSION.
3126 *  RI  X(N)  INPUT VECTOR.
3127 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
3128 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
3129 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
3130 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
3131 *         IX(I).EQ.-5.
3132 *
3133       SUBROUTINE MXUNEG(N,X,Y,IX,JOB)
3134       INTEGER N,IX(*),JOB
3135       DOUBLE PRECISION X(*),Y(*)
3136       INTEGER I
3137       IF (JOB.EQ.0) THEN
3138       DO 1 I=1,N
3139       Y(I)=-X(I)
3140     1 CONTINUE
3141       ELSE IF (JOB.GT.0) THEN
3142       DO 2 I=1,N
3143       IF (IX(I).GE. 0) THEN
3144       Y(I)=-X(I)
3145       ELSE
3146       Y(I)=0.0D 0
3147       END IF
3148     2 CONTINUE
3149       ELSE
3150       DO 3 I=1,N
3151       IF (IX(I).NE.-5) THEN
3152       Y(I)=-X(I)
3153       ELSE
3154       Y(I)=0.0D 0
3155       END IF
3156     3 CONTINUE
3157       END IF
3158       RETURN
3159       END
3160 * FUNCTION  MXUNOR               ALL SYSTEMS                99/12/01
3161 * PURPOSE :
3162 * EUCLIDEAN NORM OF A VECTOR IN A BOUND CONSTRAINED CASE.
3163 *
3164 * PARAMETERS :
3165 *  II  N  VECTOR DIMENSION.
3166 *  RI  X(N)  INPUT VECTOR.
3167 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
3168 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
3169 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
3170 *         IX(I).EQ.-5.
3171 *  RR  MXUNOR  EUCLIDEAN NORM OF X.
3172 *
3173       FUNCTION MXUNOR(N,X,IX,JOB)
3174       INTEGER N,IX(*),JOB
3175       DOUBLE PRECISION X(*),MXUNOR
3176       DOUBLE PRECISION POM,DEN
3177       INTEGER I
3178       DOUBLE PRECISION ZERO
3179       PARAMETER (ZERO=0.0D 0)
3180       DEN=ZERO
3181       IF (JOB.EQ.0) THEN
3182       DO 1 I=1,N
3183       DEN=MAX(DEN,ABS(X(I)))
3184     1 CONTINUE
3185       ELSE IF (JOB.GT.0) THEN
3186       DO 2 I=1,N
3187       IF (IX(I).GE. 0) DEN=MAX(DEN,ABS(X(I)))
3188     2 CONTINUE
3189       ELSE
3190       DO 3 I=1,N
3191       IF (IX(I).NE.-5) DEN=MAX(DEN,ABS(X(I)))
3192     3 CONTINUE
3193       END IF
3194       POM=ZERO
3195       IF (DEN.GT.ZERO) THEN
3196       IF (JOB.EQ.0) THEN
3197       DO 4 I=1,N
3198       POM=POM+(X(I)/DEN)**2
3199     4 CONTINUE
3200       ELSE IF (JOB.GT.0) THEN
3201       DO 5 I=1,N
3202       IF (IX(I).GE. 0) POM=POM+(X(I)/DEN)**2
3203     5 CONTINUE
3204       ELSE
3205       DO 6 I=1,N
3206       IF (IX(I).NE.-5) POM=POM+(X(I)/DEN)**2
3207     6 CONTINUE
3208       END IF
3209       END IF
3210       MXUNOR=DEN*SQRT(POM)
3211       RETURN
3212       END
3213 * SUBROUTINE MXUZER                ALL SYSTEMS                99/12/01
3214 * PURPOSE :
3215 * VECTOR ELEMENTS CORRESPONDING TO ACTIVE BOUNDS ARE SET TO ZERO.
3216 *
3217 * PARAMETERS :
3218 *  II  N  VECTOR DIMENSION.
3219 *  RO  X(N)  OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
3220 *  II  IX(N)  VECTOR CONTAINING TYPES OF BOUNDS.
3221 *  II  JOB  OPTION. IF JOB.GT.0 THEN INDEX I IS NOT USED WHENEVER
3222 *         IX(I).LE.-1. IF JOB.LT.0 THEN INDEX I IS NOT USED WHENEVER
3223 *         IX(I).EQ.-5.
3224 *
3225       SUBROUTINE MXUZER(N,X,IX,JOB)
3226       INTEGER N,IX(*),JOB
3227       DOUBLE PRECISION X(*)
3228       INTEGER I
3229       IF (JOB.EQ.0) RETURN
3230       DO 1 I=1,N
3231       IF (IX(I).LT.0) X(I)=0.0D 0
3232     1 CONTINUE
3233       RETURN
3234       END
3235 * SUBROUTINE MXVCOP                ALL SYSTEMS                88/12/01
3236 * PURPOSE :
3237 * COPYING OF A VECTOR.
3238 *
3239 * PARAMETERS :
3240 *  II  N  VECTOR DIMENSION.
3241 *  RI  X(N)  INPUT VECTOR.
3242 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= X.
3243 *
3244       SUBROUTINE MXVCOP(N,X,Y)
3245       INTEGER          N
3246       DOUBLE PRECISION X(*),Y(*)
3247       INTEGER          I
3248       DO 10 I = 1,N
3249           Y(I) = X(I)
3250    10 CONTINUE
3251       RETURN
3252       END
3253 * SUBROUTINE MXVCOR                ALL SYSTEMS                93/12/01
3254 * PURPOSE :
3255 * CORRECTION OF A VECTOR.
3256 *
3257 * PARAMETERS :
3258 *  II  N  VECTOR DIMENSION.
3259 *  RI  A  CORRECTION FACTOR.
3260 *  RU  X(N)  CORRECTED VECTOR. ZERO ELEMENTS OF X ARE SET TO BE EQUAL A.
3261 *
3262       SUBROUTINE MXVCOR(N,A,X)
3263       INTEGER N
3264       DOUBLE PRECISION  A,X(*)
3265       DOUBLE PRECISION ZERO
3266       PARAMETER (ZERO=0.0D 0)
3267       INTEGER I
3268       DO 1 I=1,N
3269       IF (X(I).EQ.ZERO) X(I)=A
3270     1 CONTINUE
3271       RETURN
3272       END
3273 * SUBROUTINE MXVDIF                ALL SYSTEMS                88/12/01
3274 * PURPOSE :
3275 * VECTOR DIFFERENCE.
3276 *
3277 * PARAMETERS :
3278 *  RI  X(N)  INPUT VECTOR.
3279 *  RI  Y(N)  INPUT VECTOR.
3280 *  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X - Y.
3281 *
3282       SUBROUTINE MXVDIF(N,X,Y,Z)
3283       INTEGER          N
3284       DOUBLE PRECISION X(*),Y(*),Z(*)
3285       INTEGER          I
3286       DO 10 I = 1,N
3287           Z(I) = X(I) - Y(I)
3288    10 CONTINUE
3289       RETURN
3290       END
3291 * SUBROUTINE MXVDIR                ALL SYSTEMS                91/12/01
3292 * PURPOSE :
3293 * VECTOR AUGMENTED BY THE SCALED VECTOR.
3294 *
3295 * PARAMETERS :
3296 *  II  N  VECTOR DIMENSION.
3297 *  RI  A  SCALING FACTOR.
3298 *  RI  X(N)  INPUT VECTOR.
3299 *  RI  Y(N)  INPUT VECTOR.
3300 *  RO  Z(N)  OUTPUT VECTOR WHERE Z:= Y + A*X.
3301 *
3302       SUBROUTINE MXVDIR(N,A,X,Y,Z)
3303       DOUBLE PRECISION A
3304       INTEGER          N
3305       DOUBLE PRECISION X(*),Y(*),Z(*)
3306       INTEGER          I
3307       DO 10 I = 1,N
3308           Z(I) = Y(I) + A*X(I)
3309    10 CONTINUE
3310       RETURN
3311       END
3312 * FUNCTION MXVDOT                  ALL SYSTEMS                91/12/01
3313 * PURPOSE :
3314 * DOT PRODUCT OF TWO VECTORS.
3315 *
3316 * PARAMETERS :
3317 *  II  N  VECTOR DIMENSION.
3318 *  RI  X(N)  INPUT VECTOR.
3319 *  RI  Y(N)  INPUT VECTOR.
3320 *  RR  MXVDOT  VALUE OF DOT PRODUCT MXVDOT=TRANS(X)*Y.
3321 *
3322       FUNCTION MXVDOT(N,X,Y)
3323       INTEGER          N
3324       DOUBLE PRECISION X(*),Y(*),MXVDOT
3325       DOUBLE PRECISION TEMP
3326       INTEGER          I
3327       TEMP = 0.0D0
3328       DO 10 I = 1,N
3329           TEMP = TEMP + X(I)*Y(I)
3330    10 CONTINUE
3331       MXVDOT = TEMP
3332       RETURN
3333       END
3334 * SUBROUTINE MXVICP             ALL SYSTEMS                   93/12/01
3335 * PURPOSE :
3336 * COPYING OF AN INTEGER VECTOR.
3337 *
3338 * PARAMETERS :
3339 *  II  N DIMENSION OF THE INTEGER VECTOR.
3340 *  II  IX(N)  INPUT INTEGER VECTOR.
3341 *  IO  IY(N)  OUTPUT INTEGER VECTOR SUCH THAT IY(I):= IX(I) FOR ALL I.
3342 *
3343       SUBROUTINE MXVICP(N,IX,IY)
3344       INTEGER N,IX(*),IY(*)
3345       INTEGER I
3346       DO 1 I=1,N
3347       IY(I)=IX(I)
3348     1 CONTINUE
3349       RETURN
3350       END
3351 * SUBROUTINE MXVINB            ALL SYSTEMS                   91/12/01
3352 * PURPOSE :
3353 * UPDATE OF AN INTEGER VECTOR.
3354 *
3355 * PARAMETERS :
3356 *  II  N DIMENSION OF THE INTEGER VECTOR.
3357 *  II  M DIMENSION OF THE CHANGED INTEGER VECTOR.
3358 *  II  IX(N)  INTEGER VECTOR.
3359 *  IU  JA(M)  INTEGER VECTOR WHICH IS UPDATED SO THAT JA(I)=-JA(I)
3360 *         IF IX(JA(I)).LT.0.
3361 *
3362       SUBROUTINE MXVINB(M,IX,JA)
3363       INTEGER M,IX(*),JA(*)
3364       INTEGER I
3365       DO 1 I=1,M
3366       JA(I)=ABS(JA(I))
3367       IF (IX(JA(I)).LT.0) JA(I)=-JA(I)
3368     1 CONTINUE
3369       RETURN
3370       END
3371 * SUBROUTINE MXVINE             ALL SYSTEMS                   94/12/01
3372 * PURPOSE :
3373 * ELEMENTS OF THE INTEGER VECTOR ARE REPLACED BY THEIR ABSOLUTE VALUES.
3374 *
3375 * PARAMETERS :
3376 *  II  N DIMENSION OF THE INTEGER VECTOR.
3377 *  IU  IX(N)  INTEGER VECTOR WHICH IS UPDATED SO THAT IX(I):=ABS(IX(I))
3378 *         FOR ALL I.
3379 *
3380       SUBROUTINE MXVINE(N,IX)
3381       INTEGER N,IX(*)
3382       INTEGER I
3383       DO 1 I=1,N
3384       IX(I)=ABS(IX(I))
3385     1 CONTINUE
3386       RETURN
3387       END
3388 * SUBROUTINE MXVINI             ALL SYSTEMS                   99/12/01
3389 * PURPOSE :
3390 * ELEMENTS CORRESPONDING TO FIXED VARIABLES ARE SET TO -5.
3391 *
3392 * PARAMETERS :
3393 *  II  N DIMENSION OF THE INTEGER VECTOR.
3394 *  IU  IX(N)  INTEGER VECTOR WHICH IS UPDATED SO THAT IX(I):=ABS(IX(I))
3395 *         FOR ALL I.
3396 *
3397       SUBROUTINE MXVINI(N,IX)
3398       INTEGER N,IX(*)
3399       INTEGER I
3400       DO 1 I=1,N
3401       IF (ABS(IX(I)).EQ.5) IX(I)=-5
3402     1 CONTINUE
3403       RETURN
3404       END
3405 * SUBROUTINE MXVINP             ALL SYSTEMS                   91/12/01
3406 * PURPOSE :
3407 * INITIATION OF A INTEGER PERMUTATION VECTOR.
3408 *
3409 * PARAMETERS :
3410 *  II  N DIMENSION OF THE INTEGER VECTOR.
3411 *  IO  IP(N)  INTEGER VECTOR SUCH THAT IP(I)=I FOR ALL I.
3412 *
3413       SUBROUTINE MXVINP(N,IP)
3414       INTEGER          N
3415       INTEGER          IP(*)
3416       INTEGER          I
3417       DO 10 I = 1,N
3418           IP(I) = I
3419    10 CONTINUE
3420       RETURN
3421       END
3422 * SUBROUTINE MXVINS             ALL SYSTEMS                   90/12/01
3423 * PURPOSE :
3424 * INITIATION OF THE INTEGER VECTOR.
3425 *
3426 * PARAMETERS :
3427 *  II  N DIMENSION OF THE INTEGER VECTOR.
3428 *  II  IP  INTEGER PARAMETER.
3429 *  IO  IX(N)  INTEGER VECTOR SUCH THAT IX(I)=IP FOR ALL I.
3430 *
3431       SUBROUTINE MXVINS(N,IP,IX)
3432       INTEGER          IP,N
3433       INTEGER          IX(*)
3434       INTEGER          I
3435       DO 10 I = 1,N
3436           IX(I) = IP
3437    10 CONTINUE
3438       RETURN
3439       END
3440 * SUBROUTINE MXVLIN                ALL SYSTEMS                92/12/01
3441 * PURPOSE :
3442 * LINEAR COMBINATION OF TWO VECTORS.
3443 *
3444 * PARAMETERS :
3445 *  II  N  VECTOR DIMENSION.
3446 *  RI  A  SCALING FACTOR.
3447 *  RI  X(N)  INPUT VECTOR.
3448 *  RI  B  SCALING FACTOR.
3449 *  RI  Y(N)  INPUT VECTOR.
3450 *  RO  Z(N)  OUTPUT VECTOR WHERE Z:= A*X + B*Y.
3451 *
3452       SUBROUTINE MXVLIN(N,A,X,B,Y,Z)
3453       INTEGER N
3454       DOUBLE PRECISION  A, X(*), B, Y(*), Z(*)
3455       INTEGER I
3456       DO 1 I=1,N
3457       Z(I) = A*X(I) + B*Y(I)
3458     1 CONTINUE
3459       RETURN
3460       END
3461 * FUNCTION MXVMAX             ALL SYSTEMS                   91/12/01
3462 * PURPOSE :
3463 * L-INFINITY NORM OF A VECTOR.
3464 *
3465 * PARAMETERS :
3466 *  II  N  VECTOR DIMENSION.
3467 *  RI  X(N)  INPUT VECTOR.
3468 *  RR  MXVMAX  L-INFINITY NORM OF THE VECTOR X.
3469 *
3470       FUNCTION MXVMAX(N,X)
3471       INTEGER N
3472       DOUBLE PRECISION X(*),MXVMAX
3473       INTEGER I
3474       DOUBLE PRECISION ZERO
3475       PARAMETER (ZERO=0.0D 0)
3476       MXVMAX=ZERO
3477       DO 1 I=1,N
3478       MXVMAX=MAX(MXVMAX,ABS(X(I)))
3479     1 CONTINUE
3480       RETURN
3481       END
3482 * FUNCTION MXVMX1               ALL SYSTEMS                   91/12/01
3483 * PURPOSE :
3484 * L-INFINITY NORM OF A VECTOR WITH INDEX DETERMINATION.
3485 *
3486 * PARAMETERS :
3487 *  II  N  VECTOR DIMENSION.
3488 *  RI  X(N)  INPUT VECTOR.
3489 *  IO  K  INDEX OF ELEMENT WITH MAXIMUM ABSOLUTE VALUE.
3490 *  RR  MXVMX1  L-INFINITY NORM OF THE VECTOR X.
3491 *
3492       FUNCTION MXVMX1(N,X,K)
3493       INTEGER          K,N
3494       DOUBLE PRECISION X(*),MXVMX1
3495       INTEGER          I
3496       K = 1
3497       MXVMX1 = ABS(X(1))
3498       DO 10 I = 2,N
3499           IF (ABS(X(I)).GT.MXVMX1) THEN
3500               K = I
3501               MXVMX1 = ABS(X(I))
3502           END IF
3503    10 CONTINUE
3504       RETURN
3505       END
3506 * SUBROUTINE MXVMUL             ALL SYSTEMS                   89/12/01
3507 * PURPOSE :
3508 * VECTOR IS PREMULTIPLIED BY THE POWER OF A DIAGONAL MATRIX.
3509 *
3510 * PARAMETERS :
3511 *  II  N  VECTOR DIMENSION.
3512 *  RI  D(N)  DIAGONAL MATRIX STORED AS A VECTOR WITH N ELEMENTS.
3513 *  RI  X(N)  INPUT VECTOR.
3514 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:=(D**K)*X.
3515 *  II  K  INTEGER EXPONENT.
3516 *
3517       SUBROUTINE MXVMUL(N,D,X,Y,K)
3518       INTEGER          K,N
3519       DOUBLE PRECISION D(*),X(*),Y(*)
3520       INTEGER          I
3521       IF (K.EQ.0) THEN
3522           CALL MXVCOP(N,X,Y)
3523       ELSE IF (K.EQ.1) THEN
3524           DO 10 I = 1,N
3525               Y(I) = X(I)*D(I)
3526    10     CONTINUE
3527       ELSE IF (K.EQ.-1) THEN
3528           DO 20 I = 1,N
3529               Y(I) = X(I)/D(I)
3530    20     CONTINUE
3531       ELSE
3532           DO 30 I = 1,N
3533               Y(I) = X(I)*D(I)**K
3534    30     CONTINUE
3535       END IF
3536       RETURN
3537       END
3538 * SUBROUTINE MXVNEG                ALL SYSTEMS                88/12/01
3539 * PURPOSE :
3540 * CHANGE THE SIGNS OF VECTOR ELEMENTS.
3541 *
3542 * PARAMETERS :
3543 *  II  N  VECTOR DIMENSION.
3544 *  RI  X(N)  INPUT VECTOR.
3545 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= - X.
3546 *
3547       SUBROUTINE MXVNEG(N,X,Y)
3548       INTEGER          N
3549       DOUBLE PRECISION X(*),Y(*)
3550       INTEGER          I
3551       DO 10 I = 1,N
3552           Y(I) = -X(I)
3553    10 CONTINUE
3554       RETURN
3555       END
3556 * FUNCTION  MXVNOR               ALL SYSTEMS                91/12/01
3557 * PURPOSE :
3558 * EUCLIDEAN NORM OF A VECTOR.
3559 *
3560 * PARAMETERS :
3561 *  II  N  VECTOR DIMENSION.
3562 *  RI  X(N)  INPUT VECTOR.
3563 *  RR  MXVNOR  EUCLIDEAN NORM OF X.
3564 *
3565       FUNCTION MXVNOR(N,X)
3566       INTEGER N
3567       DOUBLE PRECISION X(*),MXVNOR
3568       DOUBLE PRECISION DEN,POM
3569       INTEGER I
3570       DOUBLE PRECISION ZERO
3571       PARAMETER (ZERO=0.0D0)
3572       DEN = ZERO
3573       DO 10 I = 1,N
3574           DEN = MAX(DEN,ABS(X(I)))
3575    10 CONTINUE
3576       POM = ZERO
3577       IF (DEN.GT.ZERO) THEN
3578           DO 20 I = 1,N
3579               POM = POM + (X(I)/DEN)**2
3580    20     CONTINUE
3581       END IF
3582       MXVNOR = DEN*SQRT(POM)
3583       RETURN
3584       END
3585 * SUBROUTINE MXVSAB             ALL SYSTEMS                   91/12/01
3586 * PURPOSE :
3587 * L-1 NORM OF A VECTOR.
3588 *
3589 * PARAMETERS :
3590 *  II  N  VECTOR DIMENSION.
3591 *  RI  X(N)  INPUT VECTOR.
3592 *  RR  MXVSAB  L-1 NORM OF THE VECTOR X.
3593 *
3594       FUNCTION MXVSAB(N,X)
3595       INTEGER N
3596       DOUBLE PRECISION X(N),MXVSAB
3597       INTEGER I
3598       DOUBLE PRECISION ZERO
3599       PARAMETER (ZERO=0.0D 0)
3600       MXVSAB=ZERO
3601       DO 1 I=1,N
3602       MXVSAB=MXVSAB+ABS(X(I))
3603     1 CONTINUE
3604       RETURN
3605       END
3606 * SUBROUTINE MXVSAV                ALL SYSTEMS                91/12/01
3607 * PORTABILITY : ALL SYSTEMS
3608 * 91/12/01 LU : ORIGINAL VERSION
3609 *
3610 * PURPOSE :
3611 * DIFFERENCE OF TWO VECTORS RETURNED IN THE SUBSTRACTED ONE.
3612 *
3613 * PARAMETERS :
3614 *  II  N  VECTOR DIMENSION.
3615 *  RI  X(N)  INPUT VECTOR.
3616 *  RU  Y(N)  UPDATE VECTOR WHERE Y:= X - Y.
3617 *
3618       SUBROUTINE MXVSAV(N,X,Y)
3619       INTEGER          N
3620       DOUBLE PRECISION X(*),Y(*)
3621       DOUBLE PRECISION TEMP
3622       INTEGER          I
3623       DO 10 I = 1,N
3624           TEMP = Y(I)
3625           Y(I) = X(I) - Y(I)
3626           X(I) = TEMP
3627    10 CONTINUE
3628       RETURN
3629       END
3630 * SUBROUTINE MXVSBP                ALL SYSTEMS                91/12/01
3631 * PURPOSE :
3632 * VECTOR X(N) IS PERMUTED ACCORDING TO THE FORMULA
3633 * X(PERM(I)):=X(I).
3634 *
3635 * PARAMETERS :
3636 *  II  N  LENGTH OF VECTORS.
3637 *  II  PERM(N)  INPUT PERMUTATION VECTOR.
3638 *  RU  X(N)  VECTOR THAT IS TO BE PERMUTED.
3639 *  RA  RN01(N)  AMXILIARY VECTOR.
3640 *
3641       SUBROUTINE MXVSBP(N,PERM,X,RN01)
3642       INTEGER N,PERM(*),I
3643       DOUBLE PRECISION RN01(*),X(*)
3644       DO 100 I=1,N
3645       RN01(PERM(I))=X(I)
3646 100   CONTINUE
3647       DO 200 I=1,N
3648       X(I)=RN01(I)
3649 200   CONTINUE
3650       RETURN
3651       END
3652 * SUBROUTINE MXVSCL                ALL SYSTEMS                88/12/01
3653 * PURPOSE :
3654 * SCALING OF A VECTOR.
3655 *
3656 * PARAMETERS :
3657 *  II  N  VECTOR DIMENSION.
3658 *  RI  X(N)  INPUT VECTOR.
3659 *  RI  A  SCALING FACTOR.
3660 *  RO  Y(N)  OUTPUT VECTOR WHERE Y:= A*X.
3661 *
3662       SUBROUTINE MXVSCL(N,A,X,Y)
3663       INTEGER N
3664       DOUBLE PRECISION  A, X(*), Y(*)
3665       INTEGER I
3666       DO 1 I=1,N
3667       Y(I) = A*X(I)
3668     1 CONTINUE
3669       RETURN
3670       END
3671 * SUBROUTINE MXVSET                ALL SYSTEMS                88/12/01
3672 * PURPOSE :
3673 * A SCALAR IS SET TO ALL THE ELEMENTS OF A VECTOR.
3674 *
3675 * PARAMETERS :
3676 *  II  N  VECTOR DIMENSION.
3677 *  RI  A  INITIAL VALUE.
3678 *  RO  X(N)  OUTPUT VECTOR SUCH THAT X(I)=A FOR ALL I.
3679 *
3680       SUBROUTINE MXVSET(N,A,X)
3681       DOUBLE PRECISION A
3682       INTEGER          N
3683       DOUBLE PRECISION X(*)
3684       INTEGER          I
3685       DO 10 I = 1,N
3686           X(I) = A
3687    10 CONTINUE
3688       RETURN
3689       END
3690 * SUBROUTINE MXVSFP                ALL SYSTEMS                91/12/01
3691 * PURPOSE :
3692 * VECTOR X(N) IS PERMUTED ACCORDING TO THE FORMULA
3693 * X(I)=X(PERM(I)).
3694 *
3695 * PARAMETERS :
3696 *  II  N  LENGTH OF VECTORS.
3697 *  II  PERM(N)  INPUT PERMUTATION VECTOR.
3698 *  RU  X(N)  VECTOR THAT IS TO BE PERMUTED.
3699 *  RA  RN01(N)  AMXILIARY VECTOR.
3700 *
3701       SUBROUTINE MXVSFP(N,PERM,X,RN01)
3702       INTEGER N,PERM(*),I
3703       DOUBLE PRECISION RN01(*),X(*)
3704 C
3705       DO 100 I=1,N
3706       RN01(I)=X(PERM(I))
3707 100   CONTINUE
3708       DO 200 I=1,N
3709       X(I)=RN01(I)
3710 200   CONTINUE
3711       RETURN
3712       END
3713 * SUBROUTINE MXVSIP                ALL SYSTEMS                91/12/01
3714 * PURPOSE :
3715 * THE VECTOR OF THE INVERSE PERMUTATION IS COMPUTED.
3716 *
3717 * PARAMETERS :
3718 *  II  N  LENGTH OF VECTORS.
3719 *  II  PERM(N)  INPUT PERMUTATION VECTOR.
3720 *  IO  INVP(N)  INVERSE PERMUTATION VECTOR.
3721 *
3722       SUBROUTINE MXVSIP(N,PERM,INVP)
3723       INTEGER N,PERM(*),INVP(*)
3724       INTEGER I,J
3725       DO 100 I=1,N
3726       J=PERM(I)
3727       INVP(J)=I
3728 100   CONTINUE
3729       RETURN
3730       END
3731 * SUBROUTINE MXVSR2                ALL SYSTEMS                92/12/01
3732 * PURPOSE :
3733 * RADIXSORT.
3734 *
3735 * PARAMETERS :
3736 *  II  MCOLS  NUMBER OF INTEGER VALUES OF THE SORTED ARRAY.
3737 *  RI  DEG(MCOLS)  VALUES OF THE SORTED ARRAY.
3738 *  RO  ORD(MCOLS)  SORTED OUTPUT.
3739 *  RA  RADIX(MCOLS+1)  AUXILIARY ARRAY.
3740 *  II  WN01(MCOLS)  INDICES OF THE SORTED ARRAY.
3741 *  II  LENGTH   NUMBER OF SORTED PIECES.
3742 *
3743       SUBROUTINE MXVSR2(MCOLS,DEG,ORD,RADIX,WN01,LENGTH)
3744       INTEGER MCOLS,WN01(*)
3745       DOUBLE PRECISION DEG(*),ORD(*),RADIX(*)
3746       INTEGER LENGTH,I,L,L1,L2
3747 *
3748 *     RADIX IS SHIFTED : 0-(MCOLS-1) --- 1-MCOLS
3749 *
3750       DO 50 I=1,MCOLS+1
3751         RADIX(I)=0
3752 50    CONTINUE
3753       DO 100 I=1,LENGTH
3754         L2=WN01(I)
3755         L=DEG(L2)
3756         RADIX(L+1)=RADIX(L+1)+1
3757 100   CONTINUE
3758 *
3759 *     RADIX COUNTS THE NUMBER OF VERTICES WITH DEG(I)>=L
3760 *
3761       L=0
3762       DO 200 I=MCOLS,0,-1
3763         L=RADIX(I+1)+L
3764         RADIX(I+1)=L
3765 200   CONTINUE
3766 *
3767 *     ARRAY ORD IS FILLED
3768 *
3769       DO 300 I=1,LENGTH
3770         L2=WN01(I)
3771         L=DEG(L2)
3772         L1=RADIX(L+1)
3773         ORD(L1)=L2
3774         RADIX(L+1)=L1-1
3775 300   CONTINUE
3776       RETURN
3777       END
3778 * SUBROUTINE MXVSR5                ALL SYSTEMS                92/12/01
3779 * PURPOSE :
3780 * SHELLSORT.
3781 *
3782 * PARAMETERS :
3783 *  II  K  NUMBER OF INTEGER VALUES OF THE SORTED ARRAY.
3784 *  II  L  CORRECTION FOR THE ABSOLUTE INDEX IN THE SORTED ARRAY
3785 *  IU  ARRAY(K)  INTEGER SORTED ARRAY.
3786 *  RO  ARRAYC(K)  REAL OUTPUT ARRAY.
3787 *  RU  ARRAYD(K)  REAL ARRAY WHICH IS PERMUTED IN THE SAME WAY
3788 *         AS THE INTEGER SORTED ARRAY.
3789 *
3790       SUBROUTINE MXVSR5(K,L,ARRAY,ARRAYC,ARRAYD)
3791       INTEGER K,L
3792       INTEGER ARRAY(*)
3793       DOUBLE PRECISION ARRAYC(*),ARRAYD(*)
3794       INTEGER IS,LA,LT,LS,LLS,I,J,JS,KHALF
3795       DOUBLE PRECISION LD
3796 *
3797 *     NOTHING TO BE SORTED
3798 *
3799       IF (K.LE.1) GO TO 400
3800 *
3801 *     SHELLSORT
3802 *
3803 *     L - CORRECTION FOR THE ABSOLUTE INDEX IN THE SORTED ARRAY
3804 *
3805       LS=131071
3806       KHALF=K/2
3807       DO 300 LT=1,17
3808         IF (LS.GT.KHALF) THEN
3809           LS=LS/2
3810           GO TO 300
3811         END IF
3812         LLS=K-LS
3813         DO 200 I=1,LLS
3814           IS=I+LS
3815           LA=ARRAY(IS)
3816           LD=ARRAYD(IS)
3817           J=I
3818           JS=IS
3819  100      IF (LA.GE.ARRAY(J)) THEN
3820             ARRAY(JS)=LA
3821             ARRAYD(JS)=LD
3822             ARRAYC(INT(LD))=JS+L
3823             GO TO 200
3824           ELSE
3825             ARRAY(JS)=ARRAY(J)
3826             ARRAYD(JS)=ARRAYD(J)
3827             ARRAYC(INT(ARRAYD(J)))=JS+L
3828             JS=J
3829             J=J-LS
3830           END IF
3831           IF (J.GE.1) GO TO 100
3832           ARRAY(JS)=LA
3833           ARRAYD(JS)=LD
3834           ARRAYC(INT(LD))=JS+L
3835  200    CONTINUE
3836         LS=LS/2
3837  300  CONTINUE
3838  400  CONTINUE
3839       RETURN
3840       END
3841 * SUBROUTINE MXVSR7               ALL SYSTEMS                94/12/01
3842 * PURPOSE :
3843 * SHELLSORT
3844 *
3845 * PARAMETERS :
3846 *  II  K LENGTH OF SORTED VECTOR.
3847 *  IU  ARRAY(K) SORTED ARRAY.
3848 *  IU  ARRAYB(K) SECOND SORTED ARRAY.
3849 *
3850       SUBROUTINE MXVSR7(K,ARRAY,ARRAYB)
3851       INTEGER K
3852       INTEGER ARRAY(*),ARRAYB(*)
3853       INTEGER IS,LA,LB,LT,LS,LLS,I,J,JS,KHALF
3854 *
3855 *     NOTHING TO BE SORTED
3856 *
3857       IF (K.LE.1) GO TO 400
3858 *
3859 *     SHELLSORT
3860 *
3861       LS=131071
3862       KHALF=K/2
3863       DO 300 LT=1,17
3864         IF (LS.GT.KHALF) THEN
3865           LS=LS/2
3866           GO TO 300
3867         END IF
3868         LLS=K-LS
3869         DO 200 I=1,LLS
3870           IS=I+LS
3871           LA=ARRAY(IS)
3872           LB=ARRAYB(IS)
3873           J=I
3874           JS=IS
3875  100      IF (LA.GE.ARRAY(J)) THEN
3876             ARRAY(JS)=LA
3877             ARRAYB(JS)=LB
3878             GO TO 200
3879           ELSE
3880             ARRAY(JS)=ARRAY(J)
3881             ARRAYB(JS)=ARRAYB(J)
3882             JS=J
3883             J=J-LS
3884           END IF
3885           IF (J.GE.1) GO TO 100
3886           ARRAY(JS)=LA
3887           ARRAYB(JS)=LB
3888  200    CONTINUE
3889         LS=LS/2
3890  300  CONTINUE
3891  400  CONTINUE
3892       RETURN
3893       END
3894 * SUBROUTINE MXVSRT                ALL SYSTEMS                91/12/01
3895 * PURPOSE :
3896 * SHELLSORT
3897 *
3898 * PARAMETERS :
3899 *  II  K LENGTH OF SORTED VECTOR.
3900 *  IU  ARRAY(K) SORTED ARRAY.
3901 *
3902       SUBROUTINE MXVSRT(K,ARRAY)
3903       INTEGER K
3904       INTEGER ARRAY(*)
3905       INTEGER IS,LA,LT,LS,LLS,I,J,JS,KHALF
3906 *
3907 *     NOTHING TO BE SORTED
3908 *
3909       IF (K.LE.1) GO TO 400
3910 *
3911 *     SHELLSORT
3912 *
3913       LS=131071
3914       KHALF=K/2
3915       DO 300 LT=1,17
3916         IF (LS.GT.KHALF) THEN
3917           LS=LS/2
3918           GO TO 300
3919         END IF
3920         LLS=K-LS
3921         DO 200 I=1,LLS
3922           IS=I+LS
3923           LA=ARRAY(IS)
3924           J=I
3925           JS=IS
3926  100      IF (LA.GE.ARRAY(J)) THEN
3927             ARRAY(JS)=LA
3928             GO TO 200
3929           ELSE
3930             ARRAY(JS)=ARRAY(J)
3931             JS=J
3932             J=J-LS
3933           END IF
3934           IF (J.GE.1) GO TO 100
3935           ARRAY(JS)=LA
3936  200    CONTINUE
3937         LS=LS/2
3938  300  CONTINUE
3939  400  CONTINUE
3940       RETURN
3941       END
3942 * SUBROUTINE MXVSUM                ALL SYSTEMS                88/12/01
3943 * PURPOSE :
3944 * SUM OF TWO VECTORS.
3945 *
3946 * PARAMETERS :
3947 *  II  N  VECTOR DIMENSION.
3948 *  RI  X(N)  INPUT VECTOR.
3949 *  RI  Y(N)  INPUT VECTOR.
3950 *  RO  Z(N)  OUTPUT VECTOR WHERE Z:= X + Y.
3951 *
3952       SUBROUTINE MXVSUM(N,X,Y,Z)
3953       INTEGER          N
3954       DOUBLE PRECISION X(*),Y(*),Z(*)
3955       INTEGER          I
3956       DO 10 I = 1,N
3957           Z(I) = X(I) + Y(I)
3958    10 CONTINUE
3959       RETURN
3960       END
3961 * FUNCTION MXVVDP                  ALL SYSTEMS                92/12/01
3962 * PURPOSE :
3963 * COMPUTATION OF THE NUMBER MXVVDP=TRANS(X)*D**(-1)*Y WHERE D IS A
3964 * DIAGONAL MATRIX STORED AS A VECTOR.
3965 *
3966 * PARAMETERS :
3967 *  II  N  VECTOR DIMENSION.
3968 *  RI  D(N)  DIAGONAL MATRIX STORED AS A VECTOR.
3969 *  RI  X(N)  INPUT VECTOR.
3970 *  RI  Y(N)  INPUT VECTOR.
3971 *  RR  MXVVDP  COMPUTED NUMBER MXVVDP=TRANS(X)*D**(-1)*Y.
3972 *
3973       FUNCTION MXVVDP(N,D,X,Y)
3974       INTEGER N
3975       DOUBLE PRECISION  D(*), X(*), Y(*), MXVVDP
3976       DOUBLE PRECISION  TEMP
3977       INTEGER  I
3978       DOUBLE PRECISION ZERO
3979       PARAMETER (ZERO = 0.0D 0)
3980       TEMP = ZERO
3981       DO  1  I = 1, N
3982       TEMP = TEMP + X(I)*Y(I)/D(I)
3983     1 CONTINUE
3984       MXVVDP = TEMP
3985       RETURN
3986       END
3987 * SUBROUTINE MXWDIR                ALL SYSTEMS                92/12/01
3988 * PURPOSE :
3989 * VECTOR AUGMENTED BY THE SCALED VECTOR IN THE PACKED CASE.
3990 *
3991 * PARAMETERS :
3992 *  II  L  PACKED VECTOR DIMENSION.
3993 *  II  N  VECTOR DIMENSION.
3994 *  II  JBL(L)  INDICES OF PACKED VECTOR.
3995 *  RI  A  SCALING FACTOR.
3996 *  RI  X(L)  PACKED INPUT VECTOR.
3997 *  RI  Y(N)  UNPACKED INPUT VECTOR.
3998 *  RO  Z(N)  UNPACKED OR PACKED OUTPUT VECTOR WHERE Z:= Y + A*X.
3999 *  II  JOB  FORM OF THE VECTOR Z. JOB=1-UNPACKED FORM. JOB=2-PACKED
4000 *         FORM.
4001 *
4002       SUBROUTINE MXWDIR(L,JBL,A,X,Y,Z,JOB)
4003       INTEGER L,JBL(*),JOB
4004       DOUBLE PRECISION  A, X(*), Y(*), Z(*)
4005       INTEGER I,IP
4006       IF (JOB.EQ.1) THEN
4007       DO 1 I=1,L
4008       IP=JBL(I)
4009       IF (IP.GT.0) Z(IP)=Y(IP)+A*X(I)
4010     1 CONTINUE
4011       ELSE
4012       DO 2 I=1,L
4013       IP=JBL(I)
4014       IF (IP.GT.0) Z(I)=Y(IP)+A*X(I)
4015     2 CONTINUE
4016       END IF
4017       RETURN
4018       END
4019 * FUNCTION MXWDOT                  ALL SYSTEMS                92/12/01
4020 * PURPOSE :
4021 * DOT PRODUCT OF TWO VECTORS IN THE PACKED CASE.
4022 *
4023 * PARAMETERS :
4024 *  II  L  PACKED OR UNPACKED VECTOR DIMENSION.
4025 *  II  N  UNPACKED VECTOR DIMENSION.
4026 *  II  JBL(L)  INDICES OF PACKED VECTOR.
4027 *  RI  X(L)  UNPACKED OR PACKED INPUT VECTOR.
4028 *  RI  Y(N)  UNPACKED INPUT VECTOR.
4029 *  II  JOB  FORM OF THE VECTOR X. JOB=1-UNPACKED FORM. JOB=2-PACKED
4030 *          FORM.
4031 *  RR  MXWDOT  VALUE OF DOT PRODUCT MXWDOT=TRANS(X)*Y.
4032 *
4033       FUNCTION MXWDOT(L,JBL,X,Y,JOB)
4034       INTEGER L,JBL(*),JOB
4035       DOUBLE PRECISION  X(*), Y(*), MXWDOT
4036       DOUBLE PRECISION TEMP
4037       INTEGER I,IP
4038       TEMP=0.0D0
4039       IF (JOB.EQ.1) THEN
4040       DO 1 I=1,L
4041       IP=JBL(I)
4042       IF (IP.GT.0) TEMP=TEMP+X(IP)*Y(IP)
4043     1 CONTINUE
4044       ELSE
4045       DO 2 I=1,L
4046       IP=JBL(I)
4047       IF (IP.GT.0) TEMP=TEMP+X(I)*Y(IP)
4048     2 CONTINUE
4049       END IF
4050       MXWDOT=TEMP
4051       RETURN
4052       END