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 906200 CONTINUE 907300 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) 915400 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 926500 CONTINUE 927600 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 1185100 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 1196200 CONTINUE 1197300 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)) ) 125750 CONTINUE 125860 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 1264100 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) 1272200 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 1291300 CONTINUE 1292 GO TO 200 1293400 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 1312500 CONTINUE 1313 DO 350 K=PSL(J)+J,PSL(J+1)+J-1 1314 GAM=MAX(GAM,ABS(A(K))) 1315350 CONTINUE 1316 GAM=GAM*GAM 1317370 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 1328550 CONTINUE 1329600 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 1378100 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 1398200 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 1405300 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 1412350 I=K 1413400 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 1421500 RCHM=K 1422 DO 700 J=JSTRT,JSTOP 1423 NABOR=SL(N+1+J) 1424600 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 1432700 CONTINUE 1433 GO TO 400 1434800 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 1445900 CONTINUE 1446 GO TO 1200 14471000 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 14521100 CONTINUE 1453 NZEND=JSTRT-1 14541200 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 14651300 CONTINUE 1466 SL(K)=NZBEG 1467 WN13(K)=K 14681400 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 14751500 CONTINUE 1476 SL(N)=SL(N)+1 1477 SL(N+1)=SL(N) 1478 ML=N+SL(N+1) 1479 ISPACE=0 1480 RETURN 14811600 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 1522200 CONTINUE 1523300 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) 1531400 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 1542500 CONTINUE 1543600 CONTINUE 1544 DO 700 I=1,N 1545 X(I)=RN01(I)+X(I) 1546700 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) 1699100 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) 1729100 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 1822200 CONTINUE 1823 IF (JOB.EQ.0) X(I)=TEMP 1824300 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) 1838500 CONTINUE 1839 SUM=SUM*TEMP 1840 X(I)=X(I)-SUM 1841600 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 2187100 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 2218100 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) 2243100 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 2257200 CONTINUE 2258 END IF 2259300 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 2328100 CONTINUE 2329 DO 200 I=1,N 2330 K=ABS(IA(I)) 2331 A(K)=1.0D 0 2332200 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 2355100 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 2369200 CONTINUE 2370 END IF 2371300 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 24831470 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 25322200 CONTINUE 25332300 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 2569200 CONTINUE 2570 TEMP1=TEMP1+X(I)*TEMP2 2571 END IF 2572300 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 2614100 CONTINUE 2615 XS(I)=XS(I)+SIGMA 2616 END IF 2617200 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 2640300 CONTINUE 2641 END IF 2642400 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 2670100 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 2678200 CONTINUE 2679300 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) 2687400 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 2702500 CONTINUE 2703600 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 2712700 CONTINUE 2713800 CONTINUE 2714 DO 900 I=1,N 2715 IA(I+1)=WN11(I)+1 2716900 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 2758100 CONTINUE 2759 IA(I+1)=L 2760200 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 2810100 CONTINUE 2811 IA(I+1)=L 2812200 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) 3646100 CONTINUE 3647 DO 200 I=1,N 3648 X(I)=RN01(I) 3649200 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(*) 3704C 3705 DO 100 I=1,N 3706 RN01(I)=X(PERM(I)) 3707100 CONTINUE 3708 DO 200 I=1,N 3709 X(I)=RN01(I) 3710200 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 3728100 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 375250 CONTINUE 3753 DO 100 I=1,LENGTH 3754 L2=WN01(I) 3755 L=DEG(L2) 3756 RADIX(L+1)=RADIX(L+1)+1 3757100 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 3765200 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 3775300 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 4053