1************************-*-Mode: Fortran -*- ************************** 2* 3* thin.f -- These are the Algae sparse matrix routines. 4* 5* Copyright (C) 1994-97 K. Scott Hunziker. 6* Copyright (C) 1990-94 The Boeing Company. 7* 8* See the file COPYING for license, warranty, and permission details. 9* 10*********************************************************************** 11 12* $Id: thin.f,v 1.4 2001/10/12 00:37:52 ksh Exp $ 13 14 SUBROUTINE DGSADD( AI, AJ, AN, BI, BJ, BN, N, M, 15 $ CI, CJ, CN, W ) 16 INTEGER AI(*), AJ(*), BI(*), BJ(*), N, M 17 INTEGER CI(*), CJ(*) 18 DOUBLE PRECISION AN(*), BN(*), CN(*), W(M) 19 20* DGSADD adds two double precision matrices in the full, unordered, 21* sparse form. The CI and CJ vectors should be previously set by a 22* call to XGSADD. This routine does not change the order of CI and 23* CJ, so one can call XGSORD to order them before calling this 24* routine. 25 26* Input: AI, AJ, AN The first matrix. 27* BI, BJ, BN The second matrix. 28* CI, CJ The structure of the result. 29* N The number of rows. 30* M The number of columns. 31 32* Output: CN The resulting matrix. 33* W Workspace of length M. 34 35 INTEGER I, J, CI1, CI2 36 37* Loop over all the rows. 38 39 DO 100 I = 1, N 40 41 CI1 = CI(I) 42 CI2 = CI(I+1) - 1 43 44* Any nonzeros in this row? 45 46 IF ( CI1 .GT. CI2 ) GO TO 100 47 48* Set result columns to zero. 49 50 DO 10 J = CI1, CI2 51 W(CJ(J)) = 0.0D0 52 10 CONTINUE 53 54* Assign nonzeros from left matrix. 55 56 DO 20 J = AI(I), AI(I+1) - 1 57 W(AJ(J)) = AN(J) 58 20 CONTINUE 59 60* Add in nonzeros from right matrix. 61 62 DO 30 J = BI(I), BI(I+1) - 1 63 W(BJ(J)) = W(BJ(J)) + BN(J) 64 30 CONTINUE 65 66* Collect the results. 67 68 DO 40 J = CI1, CI2 69 CN(J) = W(CJ(J)) 70 40 CONTINUE 71 72 100 CONTINUE 73 74 RETURN 75 END 76 77*********************************************************************** 78 79 SUBROUTINE DGSMUL( AI, AJ, AN, BI, BJ, BN, NA, MB, 80 $ CI, CJ, CN, W ) 81 INTEGER AI(*), AJ(*), BI(*), BJ(*), NA, MB 82 INTEGER CI(*), CJ(*) 83 DOUBLE PRECISION AN(*), BN(*), CN(*), W(MB) 84 85* DGSMUL multiplies two matrices in the full, unordered, sparse 86* form. The structure of the result must have already been formed, 87* likely by calling XGSMUL. This routine does not change the order 88* of CI and CJ, so one can call XGSORD to order them before calling 89* this routine. 90 91* Input: AI, AJ, AN The first matrix. 92* BI, BJ, BN The second matrix. 93* CI, CJ The results matrix structure. 94* NA The number of rows of A. 95* MB The number of columns of B. 96 97* Output: CN The resulting matrix. 98* W Workspace of length MB. 99 100 INTEGER I, J, K, CI1, CI2 101 102* Loop over the rows of A. 103 104 DO 100 I = 1, NA 105 106 CI1 = CI(I) 107 CI2 = CI(I+1) - 1 108 109* Any nonzeros in this row? 110 111 IF ( CI1 .GT. CI2 ) GO TO 100 112 113* Set result columns to zero. 114 115 DO 10 J = CI1, CI2 116 W(CJ(J)) = 0.0D0 117 10 CONTINUE 118 119* Step through the nonzeros in this row of A. 120 121 DO 30 J = AI(I), AI(I+1) - 1 122 123* Work through the nonzeros in the corresponding row of B. 124 125 DO 20 K = BI(AJ(J)), BI(AJ(J)+1) - 1 126 W(BJ(K)) = W(BJ(K)) + AN(J)*BN(K) 127 20 CONTINUE 128 129 30 CONTINUE 130 131* Collect the results. 132 133 DO 40 J = CI1, CI2 134 CN(J) = W(CJ(J)) 135 40 CONTINUE 136 137 100 CONTINUE 138 139 RETURN 140 END 141 142*********************************************************************** 143 144 SUBROUTINE DGSSUB( AI, AJ, AN, BI, BJ, BN, N, M, 145 $ CI, CJ, CN, W ) 146 INTEGER AI(*), AJ(*), BI(*), BJ(*), N, M 147 INTEGER CI(*), CJ(*) 148 DOUBLE PRECISION AN(*), BN(*), CN(*), W(M) 149 150* DGSSUB subtracts two double precision matrices in the full, 151* unordered, sparse form. The CI and CJ vectors should be 152* previously set by a call to XGSADD. This routine does not change 153* the order of CI and CJ, so one can call XGSORD to order them 154* before calling this routine. 155 156* Input: AI, AJ, AN The first matrix. 157* BI, BJ, BN The second matrix. 158* CI, CJ The structure of the result. 159* N The number of rows. 160* M The number of columns. 161 162* Output: CN The resulting matrix. 163* W Workspace of length M. 164 165 INTEGER I, J, CI1, CI2 166 167* Loop over all the rows. 168 169 DO 100 I = 1, N 170 171 CI1 = CI(I) 172 CI2 = CI(I+1) - 1 173 174* Any nonzeros in this row? 175 176 IF ( CI1 .GT. CI2 ) GO TO 100 177 178* Set result columns to zero. 179 180 DO 10 J = CI1, CI2 181 W(CJ(J)) = 0.0D0 182 10 CONTINUE 183 184* Assign nonzeros from left matrix. 185 186 DO 20 J = AI(I), AI(I+1) - 1 187 W(AJ(J)) = AN(J) 188 20 CONTINUE 189 190* Subtract nonzeros from right matrix. 191 192 DO 30 J = BI(I), BI(I+1) - 1 193 W(BJ(J)) = W(BJ(J)) - BN(J) 194 30 CONTINUE 195 196* Collect the results. 197 198 DO 40 J = CI1, CI2 199 CN(J) = W(CJ(J)) 200 40 CONTINUE 201 202 100 CONTINUE 203 204 RETURN 205 END 206 207*********************************************************************** 208 209 SUBROUTINE DGSTRN( AI, AJ, AN, N, M, ATI, ATJ, ATN ) 210 INTEGER AI(*), AJ(*), ATI(*), ATJ(*), N, M 211 DOUBLE PRECISION AN(*), ATN(*) 212 213* DGSTRN transposes a double precision matrix in the full, 214* unordered, sparse form. 215 216* Input: AI, AJ, AN The original matrix. 217* N The number of rows. 218* M The number of columns. 219 220* Output: ATI, ATJ, ATN The transposed matrix. 221 222 INTEGER I, J, K, JJ 223 224 DO 10 I = 2, M+1 225 ATI(I) = 0 226 10 CONTINUE 227 228 DO 20 I = 1, AI(N+1) - 1 229 J = AJ(I) + 2 230 IF ( J .LE. M+1 ) ATI(J) = ATI(J) + 1 231 20 CONTINUE 232 233 ATI(1) = 1 234 ATI(2) = 1 235 236 DO 30 I = 3, M+1 237 ATI(I) = ATI(I) + ATI(I-1) 238 30 CONTINUE 239 240 DO 50 I = 1, N 241 242 DO 40 J = AI(I), AI(I+1) - 1 243 244 JJ = AJ(J) + 1 245 K = ATI(JJ) 246 247 ATJ(K) = I 248 ATN(K) = AN(J) 249 ATI(JJ) = K + 1 250 251 40 CONTINUE 252 253 50 CONTINUE 254 255 RETURN 256 END 257 258*********************************************************************** 259 260 SUBROUTINE IGSADD( AI, AJ, AN, BI, BJ, BN, N, M, 261 $ CI, CJ, CN, W ) 262 INTEGER AI(*), AJ(*), BI(*), BJ(*), N, M 263 INTEGER CI(*), CJ(*) 264 INTEGER AN(*), BN(*), CN(*), W(M) 265 266* IGSADD adds two integer matrices in the full, unordered, sparse 267* form. The CI and CJ vectors should be previously set by a call 268* to XGSADD. This routine does not change the order of CI and CJ, 269* so one can call XGSORD to order them before calling this routine. 270 271* Input: AI, AJ, AN The first matrix. 272* BI, BJ, BN The second matrix. 273* CI, CJ The structure of the result. 274* N The number of rows. 275* M The number of columns. 276 277* Output: CN The resulting matrix. 278* W Workspace of length M. 279 280 INTEGER I, J, CI1, CI2 281 282* Loop over all the rows. 283 284 DO 100 I = 1, N 285 286 CI1 = CI(I) 287 CI2 = CI(I+1) - 1 288 289* Any nonzeros in this row? 290 291 IF ( CI1 .GT. CI2 ) GO TO 100 292 293* Set result columns to zero. 294 295 DO 10 J = CI1, CI2 296 W(CJ(J)) = 0 297 10 CONTINUE 298 299* Assign nonzeros from left matrix. 300 301 DO 20 J = AI(I), AI(I+1) - 1 302 W(AJ(J)) = AN(J) 303 20 CONTINUE 304 305* Add in nonzeros from right matrix. 306 307 DO 30 J = BI(I), BI(I+1) - 1 308 W(BJ(J)) = W(BJ(J)) + BN(J) 309 30 CONTINUE 310 311* Collect the results. 312 313 DO 40 J = CI1, CI2 314 CN(J) = W(CJ(J)) 315 40 CONTINUE 316 317 100 CONTINUE 318 319 RETURN 320 END 321 322*********************************************************************** 323 324 SUBROUTINE IGSMUL( AI, AJ, AN, BI, BJ, BN, NA, MB, 325 $ CI, CJ, CN, W ) 326 INTEGER AI(*), AJ(*), BI(*), BJ(*), NA, MB 327 INTEGER CI(*), CJ(*) 328 INTEGER AN(*), BN(*), CN(*), W(MB) 329 330* IGSMUL multiplies two matrices in the full, unordered, sparse 331* form. The structure of the result must have already been formed, 332* likely by calling XGSMUL. 333 334* Input: AI, AJ, AN The first matrix. 335* BI, BJ, BN The second matrix. 336* CI, CJ The results matrix structure. 337* NA The number of rows of A. 338* MB The number of columns of B. 339 340* Output: CN The resulting matrix. 341* W Workspace of length MB. 342 343 INTEGER I, J, K, CI1, CI2 344 345* Loop over the rows of A. 346 347 DO 100 I = 1, NA 348 349 CI1 = CI(I) 350 CI2 = CI(I+1) - 1 351 352 IF ( CI1 .GT. CI2 ) GO TO 100 353 354* Set result columns to zero. 355 356 DO 10 J = CI1, CI2 357 W(CJ(J)) = 0 358 10 CONTINUE 359 360* Step through the nonzeros in this row of A. 361 362 DO 30 J = AI(I), AI(I+1) - 1 363 364* Work through the nonzeros in the corresponding row of B. 365 366 DO 20 K = BI(AJ(J)), BI(AJ(J)+1) - 1 367 W(BJ(K)) = W(BJ(K)) + AN(J)*BN(K) 368 20 CONTINUE 369 370 30 CONTINUE 371 372* Collect the results. 373 374 DO 40 J = CI1, CI2 375 CN(J) = W(CJ(J)) 376 40 CONTINUE 377 378 100 CONTINUE 379 380 RETURN 381 END 382 383*********************************************************************** 384 385 SUBROUTINE IGSSUB( AI, AJ, AN, BI, BJ, BN, N, M, 386 $ CI, CJ, CN, W ) 387 INTEGER AI(*), AJ(*), BI(*), BJ(*), N, M 388 INTEGER CI(*), CJ(*) 389 INTEGER AN(*), BN(*), CN(*), W(M) 390 391* IGSSUB subtracts two integer matrices in the full, unordered 392* sparse form. The CI and CJ vectors should be previously set by a 393* call to XGSADD. This routine does not change the order of CI and 394* CJ, so one can call XGSORD to order them before calling this 395* routine. 396 397* Input: AI, AJ, AN The first matrix. 398* BI, BJ, BN The second matrix. 399* CI, CJ The structure of the result. 400* N The number of rows. 401* M The number of columns. 402 403* Output: CN The resulting matrix. 404* W Workspace of length M. 405 406 INTEGER I, J, CI1, CI2 407 408 DO 100 I = 1, N 409 410 CI1 = CI(I) 411 CI2 = CI(I+1) - 1 412 413* Any nonzeros in this row? 414 415 IF ( CI1 .GT. CI2 ) GO TO 100 416 417* Set result columns to zero. 418 419 DO 10 J = CI1, CI2 420 W(CJ(J)) = 0 421 10 CONTINUE 422 423* Assign nonzeros from left matrix. 424 425 DO 20 J = AI(I), AI(I+1) - 1 426 W(AJ(J)) = AN(J) 427 20 CONTINUE 428 429* Subtract nonzeros from right matrix. 430 431 DO 30 J = BI(I), BI(I+1) - 1 432 W(BJ(J)) = W(BJ(J)) - BN(J) 433 30 CONTINUE 434 435* Collect the results. 436 437 DO 40 J = CI1, CI2 438 CN(J) = W(CJ(J)) 439 40 CONTINUE 440 441 100 CONTINUE 442 443 RETURN 444 END 445 446*********************************************************************** 447 448 SUBROUTINE IGSTRN( AI, AJ, AN, N, M, ATI, ATJ, ATN ) 449 INTEGER AI(*), AJ(*), ATI(*), ATJ(*), N, M 450 INTEGER AN(*), ATN(*) 451 452* IGSTRN transposes an integer matrix in the full, unordered, sparse 453* form. 454 455* Input: AI, AJ, AN The original matrix. 456* N The number of rows. 457* M The number of columns. 458 459* Output: ATI, ATJ, ATN The transposed matrix. 460 461 INTEGER I, J, K, JJ 462 463 DO 10 I = 2, M+1 464 ATI(I) = 0 465 10 CONTINUE 466 467 DO 20 I = 1, AI(N+1) - 1 468 J = AJ(I) + 2 469 IF ( J .LE. M+1 ) ATI(J) = ATI(J) + 1 470 20 CONTINUE 471 472 ATI(1) = 1 473 ATI(2) = 1 474 475 DO 30 I = 3, M+1 476 ATI(I) = ATI(I) + ATI(I-1) 477 30 CONTINUE 478 479 DO 50 I = 1, N 480 481 DO 40 J = AI(I), AI(I+1) - 1 482 483 JJ = AJ(J) + 1 484 K = ATI(JJ) 485 486 ATJ(K) = I 487 ATN(K) = AN(J) 488 ATI(JJ) = K + 1 489 490 40 CONTINUE 491 492 50 CONTINUE 493 494 RETURN 495 END 496 497*********************************************************************** 498 499 SUBROUTINE XGSADD( AI, AJ, BI, BJ, N, M, CI, CJ, IW ) 500 INTEGER AI(*), AJ(*), BI(*), BJ(*), N, M 501 INTEGER CI(*), CJ(*), IW(M) 502 503* XGSADD symbolicly adds two matrices in the full, unordered sparse 504* form. 505 506* Input: AI, AJ The first matrix structure. 507* BI, BJ The second matrix structure. 508* N The number of rows. 509* M The number of columns. 510 511* Output: CI, CJ The resulting matrix structure. 512* IW Integer workspace of length M. 513 514 INTEGER I, J, K 515 516 K = 1 517 518 DO 10 I = 1, M 519 IW(I) = 0 520 10 CONTINUE 521 522 DO 40 I = 1, N 523 524 CI(I) = K 525 526 DO 20 J = AI(I), AI(I+1) - 1 527 CJ(K) = AJ(J) 528 IW(AJ(J)) = I 529 K = K + 1 530 20 CONTINUE 531 532 DO 30 J = BI(I), BI(I+1) - 1 533 IF ( IW(BJ(J)) .EQ. I ) GO TO 30 534 CJ(K) = BJ(J) 535 K = K + 1 536 30 CONTINUE 537 538 40 CONTINUE 539 540 CI(N+1) = K 541 542 RETURN 543 END 544 545*********************************************************************** 546 547 SUBROUTINE XGSMUL( AI, AJ, BI, BJ, NA, MB, CI, CJ, MAXCJ, IW ) 548 INTEGER AI(*), AJ(*), BI(*), BJ(*), NA, MB 549 INTEGER CI(*), CJ(MAXCJ), IW(MB) 550 551* XGSMUL symbolicly multiplies two matrices in full, unordered 552* sparse form. If CJ is not dimensioned large enough, we return 553* with CI(1)=0. 554 555* Input: AI, AJ The first matrix structure. 556* BI, BJ The second matrix structure. 557* NA The number of rows of A. 558* MB The number of columns of B. 559 560* Output: CI, CJ The resulting matrix structure. 561* MAXCJ The dimensioned length of CJ. 562* IW Integer workspace of length MB. 563 564 INTEGER I, J, K, JJ 565* 566 JJ = 1 567 568* Initialize 569 570 DO 10 I = 1, MB 571 IW(I) = 0 572 10 CONTINUE 573 574* Loop over rows of left matrix. 575 576 DO 40 I = 1, NA 577 578 CI(I) = JJ 579 580 DO 30 J = AI(I), AI(I+1) - 1 581 582 DO 20 K = BI(AJ(J)), BI(AJ(J)+1) - 1 583 584 IF ( IW(BJ(K)) .EQ. I ) GO TO 20 585 586 IF ( JJ .GT. MAXCJ ) THEN 587 CI(1) = 0 588 GO TO 50 589 ENDIF 590 591 CJ(JJ) = BJ(K) 592 JJ = JJ + 1 593 IW(BJ(K)) = I 594 595 20 CONTINUE 596 30 CONTINUE 597 40 CONTINUE 598 599 CI(NA+1) = JJ 600 601 50 CONTINUE 602 RETURN 603 END 604 605*********************************************************************** 606 607 SUBROUTINE XGSTRN( AI, AJ, N, M, ATI, ATJ ) 608 INTEGER AI(*), AJ(*), ATI(*), ATJ(*), N, M 609 610* XGSTRN symbolically transposes a matrix in the full, unordered, 611* sparse form. In the process, the rows of the matrix are ordered, 612* so that two calls to this routine will order the matrix. 613 614* Input: AI, AJ The original matrix structure. 615* N The number of rows. 616* M The number of columns. 617 618* Output: ATI, ATJ The transposed matrix structure. 619 620 INTEGER I, J 621 622 DO 10 I = 2, M+1 623 ATI(I) = 0 624 10 CONTINUE 625 626 DO 20 I = 1, AI(N+1) - 1 627 J = AJ(I) + 2 628 IF ( J .LE. M+1 ) ATI(J) = ATI(J) + 1 629 20 CONTINUE 630 631 ATI(1) = 1 632 ATI(2) = 1 633 634 DO 30 I = 3, M+1 635 ATI(I) = ATI(I) + ATI(I-1) 636 30 CONTINUE 637 638 DO 50 I = 1, N 639 640 DO 40 J = AI(I), AI(I+1) - 1 641 ATJ( ATI( AJ(J)+1 ) ) = I 642 ATI( AJ(J)+1 ) = ATI( AJ(J)+1 ) + 1 643 40 CONTINUE 644 645 50 CONTINUE 646 647 RETURN 648 END 649 650*********************************************************************** 651 652 SUBROUTINE ZGSADD( AI, AJ, AN, BI, BJ, BN, N, M, 653 $ CI, CJ, CN, W ) 654 INTEGER AI(*), AJ(*), BI(*), BJ(*), N, M 655 INTEGER CI(*), CJ(*) 656 DOUBLE PRECISION AN(*), BN(*), CN(*), W(*) 657 658* ZGSADD adds matrices in the full, unordered, sparse form. The CI 659* and CJ vectors should be previously set by a call to XGSADD. This 660* routine does not change the order of CI and CJ, so one could call 661* XGSORD to order them before calling this routine. 662 663* Double precision complex data is stored as two consecutive double 664* precision numbers. (This is equivalent to the COMPLEX*16 data 665* type that most FORTRAN compilers provide but which is not part 666* of the standard language.) 667 668* Input: AI, AJ, AN The first matrix. 669* BI, BJ, BN The second matrix. 670* CI, CJ The structure of the result. 671* N The number of rows. 672* M The number of columns. 673 674* Output: CN The resulting matrix. 675* W Workspace of length M. 676 677 INTEGER I, J, K, CI1, CI2 678 679 DO 50 I = 1, N 680 681 CI1 = CI(I) 682 CI2 = CI(I+1) - 1 683 684 IF ( CI1 .GT. CI2 ) GO TO 50 685 686 DO 10 K = CI1, CI2 687 W(2*CJ(K)-1) = 0.0D0 688 W(2*CJ(K)) = 0.0D0 689 10 CONTINUE 690 691 DO 20 K = AI(I), AI(I+1) - 1 692 W(2*AJ(K)-1) = AN(2*K-1) 693 W(2*AJ(K)) = AN(2*K) 694 20 CONTINUE 695 696 DO 30 K = BI(I), BI(I+1) - 1 697 J = BJ(K) 698 W(2*J-1) = W(2*J-1) + BN(2*K-1) 699 W(2*J) = W(2*J) + BN(2*K) 700 30 CONTINUE 701 702 DO 40 K = CI1, CI2 703 CN(2*K-1) = W(2*CJ(K)-1) 704 CN(2*K) = W(2*CJ(K)) 705 40 CONTINUE 706 707 50 CONTINUE 708 709 RETURN 710 END 711 712*********************************************************************** 713 714 SUBROUTINE ZGSMUL( AI, AJ, AN, BI, BJ, BN, NA, MB, 715 $ CI, CJ, CN, W ) 716 INTEGER AI(*), AJ(*), BI(*), BJ(*), NA, MB 717 INTEGER CI(*), CJ(*) 718 DOUBLE PRECISION AN(*), BN(*), CN(*), W(MB) 719 720* ZGSMUL multiplies two matrices in the full, unordered, sparse 721* form. The structure of the result must have already been formed, 722* likely by calling XGSMUL. 723 724* Input: AI, AJ, AN The first matrix. 725* BI, BJ, BN The second matrix. 726* CI, CJ The results matrix structure. 727* NA The number of rows of A. 728* MB The number of columns of B. 729 730* Output: CN The resulting matrix. 731* W Workspace of length MB. 732 733 INTEGER I, J, K, KK, CI1, CI2 734 DOUBLE PRECISION ARE, AIM 735 736 DO 50 I = 1, NA 737 738 CI1 = CI(I) 739 CI2 = CI(I+1) - 1 740 741 IF ( CI1 .GT. CI2 ) GO TO 50 742 743 DO 10 J = CI1, CI2 744 W(2*CJ(J)-1) = 0.0D0 745 W(2*CJ(J)) = 0.0D0 746 10 CONTINUE 747 748 DO 30 J = AI(I), AI(I+1) - 1 749 750 ARE = AN(2*J-1) 751 AIM = AN(2*J) 752 753 DO 20 KK = BI(AJ(J)), BI(AJ(J)+1) - 1 754 K = BJ(KK) 755 W(2*K-1) = W(2*K-1) + ARE*BN(2*KK-1)-AIM*BN(2*KK) 756 W(2*K) = W(2*K) + ARE*BN(2*KK)+AIM*BN(2*KK-1) 757 20 CONTINUE 758 759 30 CONTINUE 760 761 DO 40 J = CI1, CI2 762 CN(2*J-1) = W(2*CJ(J)-1) 763 CN(2*J) = W(2*CJ(J)) 764 40 CONTINUE 765 766 50 CONTINUE 767 768 RETURN 769 END 770 771*********************************************************************** 772 773 SUBROUTINE ZGSSUB( AI, AJ, AN, BI, BJ, BN, N, M, 774 $ CI, CJ, CN, W ) 775 INTEGER AI(*), AJ(*), BI(*), BJ(*), N, M 776 INTEGER CI(*), CJ(*) 777 DOUBLE PRECISION AN(*), BN(*), CN(*), W(*) 778 779* ZGSSUB subtracts double precision complex matrices in the full, 780* unordered, sparse form. The CI and CJ vectors should be 781* previously set by a call to XGSADD. This routine does not change 782* the order of CI and CJ, so one could call XGSORD to order them 783* before calling this routine. 784 785* Input: AI, AJ, AN The first matrix. 786* BI, BJ, BN The second matrix. 787* CI, CJ The structure of the result. 788* N The number of rows. 789* M The number of columns. 790 791* Output: CN The resulting matrix. 792* W Workspace of length M. 793 794 INTEGER I, J, K, CI1, CI2 795 796 DO 50 I = 1, N 797 798 CI1 = CI(I) 799 CI2 = CI(I+1) - 1 800 801 IF ( CI1 .GT. CI2 ) GO TO 50 802 803 DO 10 J = CI1, CI2 804 W(2*CJ(J)-1) = 0.0D0 805 W(2*CJ(J)) = 0.0D0 806 10 CONTINUE 807 808 DO 20 J = AI(I), AI(I+1) - 1 809 W(2*AJ(J)-1) = AN(2*J-1) 810 W(2*AJ(J)) = AN(2*J) 811 20 CONTINUE 812 813 DO 30 J = BI(I), BI(I+1) - 1 814 K = BJ(J) 815 W(2*K-1) = W(2*K-1) - BN(2*J-1) 816 W(2*K) = W(2*K) - BN(2*J) 817 30 CONTINUE 818 819 DO 40 J = CI1, CI2 820 CN(2*J-1) = W(2*CJ(J)-1) 821 CN(2*J) = W(2*CJ(J)) 822 40 CONTINUE 823 824 50 CONTINUE 825 826 RETURN 827 END 828 829*********************************************************************** 830 831 SUBROUTINE ZGSTRN( AI, AJ, AN, N, M, ATI, ATJ, ATN ) 832 INTEGER AI(*), AJ(*), ATI(*), ATJ(*), N, M 833 DOUBLE PRECISION AN(*), ATN(*) 834 835* DGSTRN transposes a double precision complex matrix in the 836* full, unordered sparse form. 837 838* Input: AI, AJ, AN The original matrix. 839* N The number of rows. 840* M The number of columns. 841 842* Output: ATI, ATJ, ATN The transposed matrix. 843 844 INTEGER I, J, K, JP 845 846 DO 10 I = 2, M+1 847 ATI(I) = 0 848 10 CONTINUE 849 850 DO 20 I = 1, AI(N+1) - 1 851 J = AJ(I) + 2 852 IF ( J .LE. M+1 ) ATI(J) = ATI(J) + 1 853 20 CONTINUE 854 855 ATI(1) = 1 856 ATI(2) = 1 857 858 DO 30 I = 3, M+1 859 ATI(I) = ATI(I) + ATI(I-1) 860 30 CONTINUE 861 862 DO 50 I = 1, N 863 DO 40 JP = AI(I), AI(I+1) - 1 864 J = AJ(JP) + 1 865 K = ATI(J) 866 ATJ(K) = I 867 ATN(2*K-1) = AN(2*JP-1) 868 ATN(2*K) = AN(2*JP) 869 ATI(J) = K + 1 870 40 CONTINUE 871 50 CONTINUE 872 873 RETURN 874 END 875 876*********************************************************************** 877 878 SUBROUTINE DPSSOL( IU, JU, UN, DI, N, B, X ) 879 INTEGER IU(*), JU(*), N 880 DOUBLE PRECISION UN(*), DI(N), B(N), X(N) 881 882* DPSSOL solves a system of linear equations "U'*D*U*x=b" for "x", 883* where "U" is an upper triangular, double precision sparse matrix 884* in upper, unordered form. The "U" matrix likely comes from calls 885* to XPSFAC, XGSORD, and DPSFAC. "b" is a full input vector, and 886* "x" is the full output vector. 887 888* Input: IU, JU, UN The upper triangular matrix "U". Its 889* diagonal is assumed to be identity, 890* and is not used. 891* DI The inverses of the diagonal elements 892* of D. 893* N The order of the system. 894* B The right-hand-side vector. 895 896* Output: X The solution vector. 897 898 INTEGER I, K, NM 899 DOUBLE PRECISION XX 900 901 NM = N - 1 902 DO 10 I = 1, N 903 X(I) = B(I) 904 10 CONTINUE 905 DO 40 K = 1, NM 906 XX = X(K) 907 DO 20 I = IU(K), IU(K+1) - 1 908 X(JU(I)) = X(JU(I)) - UN(I)*XX 909 20 CONTINUE 910 X(K) = XX*DI(K) 911 40 CONTINUE 912 X(N) = X(N)*DI(N) 913 K = NM 914 50 CONTINUE 915 XX = X(K) 916 DO 60 I = IU(K), IU(K+1) - 1 917 XX = XX - UN(I)*X(JU(I)) 918 60 CONTINUE 919 X(K) = XX 920 K = K - 1 921 IF ( K .GT. 0 ) GO TO 50 922 923 RETURN 924 END 925 926*********************************************************************** 927 928 SUBROUTINE DPSFAC( AI, AJ, AN, AD, N, IU, JU, UN, DI, IP, IUP ) 929 INTEGER AI(*), AJ(*), N 930 INTEGER IU(*), JU(*), IP(N), IUP(N) 931 DOUBLE PRECISION AN(*), AD(N), UN(*), DI(N) 932 933* DPSFAC performs a triangular factorization of a symmetric, 934* positive definite, sparse matrix in upper, unordered form. The 935* structure of the result must already exist, and it must be in 936* RR(U)O form. Usually, one calls XPSFAC, then XGSORD, and then 937* DPSFAC. The factorization satisfies A=U'*D*U. The matrix 938* returned is in RR(U)O form, has 1/D on the diagonal and the 939* off-diag elements of U in the upper triangle. (The diagonal 940* elements of U are all equal to one.) 941 942* Input: AI, AJ, AN, AD The matrix to be factored. 943* N The order of the matrix. 944* IU, JU The matrix U structure. 945 946* Output: UN The matrix U off-diagonal elements. 947* DI The inverse of matrix D. 948* IP Integer workspace of length N. 949* IUP Integer workspace of length N. The 950* value of IUP(1) is zero on a successful 951* return or nonzero on an error. 952 953 INTEGER I, J, L, JJ, IH, IUA, IUB, AI1, AI2, LAST, LN 954 INTEGER IUC, IUD 955 DOUBLE PRECISION UM 956 957 DO 10 J = 1, N 958 IP(J) = 0 959 10 CONTINUE 960 DO 130 I = 1, N 961 IH = I + 1 962 IUA = IU(I) 963 IUB = IU(IH) - 1 964 IF ( IUB .LT. IUA ) GO TO 40 965 DO 20 J = IUA, IUB 966 DI(JU(J)) = 0.0D0 967 20 CONTINUE 968 AI1 = AI(I) 969 AI2 = AI(IH) - 1 970 IF ( AI2 .LT. AI1 ) GO TO 40 971 DO 30 J = AI1, AI2 972 DI(AJ(J)) = AN(J) 973 30 CONTINUE 974 40 CONTINUE 975 DI(I) = AD(I) 976 LAST = IP(I) 977 IF ( LAST .EQ. 0 ) GO TO 90 978 LN = IP(LAST) 979 50 CONTINUE 980 L = LN 981 LN = IP(L) 982 IUC = IUP(L) 983 IUD = IU(L+1) - 1 984 UM = UN(IUC)*DI(L) 985 DO 60 J = IUC, IUD 986 JJ = JU(J) 987 DI(JJ) = DI(JJ) - UN(J)*UM 988 60 CONTINUE 989 UN(IUC) = UM 990 IUP(L) = IUC + 1 991 IF ( IUC .EQ. IUD ) GO TO 80 992 J = JU(IUC+1) 993 JJ = IP(J) 994 IF ( JJ .EQ. 0 ) GO TO 70 995 IP(L) = IP(JJ) 996 IP(JJ) = L 997 GO TO 80 998 70 CONTINUE 999 IP(J) = L 1000 IP(L) = L 1001 80 CONTINUE 1002 IF ( L .NE. LAST ) GO TO 50 1003 90 CONTINUE 1004 IF ( 1.0D0 + DI(I) .EQ. 1.0D0 ) THEN 1005 IUP(1) = -1 1006 RETURN 1007 ENDIF 1008 DI(I) = 1.0D0/DI(I) 1009 IF ( IUB .LT. IUA ) GO TO 120 1010 DO 100 J = IUA, IUB 1011 UN(J) = DI(JU(J)) 1012 100 CONTINUE 1013 J = JU(IUA) 1014 JJ = IP(J) 1015 IF ( JJ .EQ. 0 ) GO TO 110 1016 IP(I) = IP(JJ) 1017 IP(JJ) = I 1018 GO TO 120 1019 110 CONTINUE 1020 IP(J) = I 1021 IP(I) = I 1022 120 CONTINUE 1023 IUP(I) = IUA 1024 130 CONTINUE 1025 1026 IUP(1) = 0 1027 RETURN 1028 END 1029 1030*********************************************************************** 1031 1032 SUBROUTINE XPSFAC( AI, AJ, N, NN, IU, JU, IP ) 1033 INTEGER AI(*), AJ(*), N, NN 1034 INTEGER IU(*), JU(*), IP(N) 1035 1036* XPSFAC performs a symbolic triangular factorization of a 1037* symmetric, positive definite, sparse matrix in upper, unordered 1038* form. 1039 1040* Input: AI, AJ The matrix to be factored. 1041* N The order of the matrix. 1042* NN The length of JU. 1043 1044* Output: IU, JU The resulting matrix structure. 1045* IP Integer workspace of length N. 1046 1047* If the required length of JU is greater than NN, then the routine 1048* returns with IU(1)=0 and should be called again with a longer JU. 1049 1050 INTEGER I, JJ, L, NM, NH, JP, JPI, JPP, MN, AI1, AI2, LAST 1051 INTEGER LH, IUA, IUB 1052 1053 NM = N - 1 1054 NH = N + 1 1055 DO 10 I = 1, N 1056 IU(I) = 0 1057 IP(I) = 0 1058 10 CONTINUE 1059 JP = 1 1060 DO 90 I = 1, NM 1061 JPI = JP 1062 JPP = N + JP - 1 1063 MN = NH 1064 AI1 = AI(I) 1065 AI2 = AI(I+1) - 1 1066 IF ( AI2 .LT. AI1 ) GO TO 30 1067 DO 20 J = AI1, AI2 1068 JJ = AJ(J) 1069 IF ( JP .GT. NN ) GO TO 100 1070 JU(JP) = JJ 1071 JP = JP + 1 1072 IF ( JJ .LT. MN ) MN = JJ 1073 IU(JJ) = I 1074 20 CONTINUE 1075 30 CONTINUE 1076 LAST = IP(I) 1077 IF ( LAST .EQ. 0 ) GO TO 60 1078 L = LAST 1079 40 CONTINUE 1080 L = IP(L) 1081 LH = L + 1 1082 IUA = IU(L) 1083 IUB = IU(LH) - 1 1084 IF ( LH .EQ. I ) IUB = JPI - 1 1085 IU(I) = I 1086 DO 50 J = IUA, IUB 1087 JJ = JU(J) 1088 IF ( IU(JJ) .EQ. I ) GO TO 50 1089 IF ( JP .GT. NN ) GO TO 100 1090 JU(JP) = JJ 1091 JP = JP + 1 1092 IU(JJ) = I 1093 IF ( JJ .LT. MN ) MN = JJ 1094 50 CONTINUE 1095 IF ( JP .EQ. JPP ) GO TO 70 1096 IF ( L .NE. LAST ) GO TO 40 1097 60 CONTINUE 1098 IF ( MN .EQ. NH ) GO TO 90 1099 70 CONTINUE 1100 L = IP(MN) 1101 IF ( L .EQ. 0 ) GO TO 80 1102 IP(I) = IP(L) 1103 IP(L) = I 1104 GO TO 90 1105 80 CONTINUE 1106 IP(MN) = I 1107 IP(I) = I 1108 90 IU(I) = JPI 1109 IU(N) = JP 1110 IU(NH) = JP 1111 RETURN 1112 1113 100 CONTINUE 1114 IU(1) = 0 1115 RETURN 1116 1117 END 1118 1119*********************************************************************** 1120 1121 SUBROUTINE ISSTRF( N, M, AI, AJ, AN, AD, P, C, W ) 1122 INTEGER N, AI(*), AJ(*) 1123 INTEGER AN(*), AD(*), P(N,M), C(M,M), W(N,M) 1124 1125* This is the "transform" routine. It multiplies P'*A*P to form the 1126* result C, where P and C are dense and A is an INTEGER, sparse, 1127* symmetric matrix in upper, unordered form. 1128 1129* Input: N The order of A. 1130* M The number of columns of P. 1131* AI,AJ,AN,AD The N by N matrix A. 1132* P The N by M matrix P. 1133* W A work matrix, N by M. 1134 1135* Output: C The M by M result. 1136 1137 DO 20 I = 1, N 1138 DO 10 J = 1, M 1139 W(I,J) = AD(I)*P(I,J) 1140 10 CONTINUE 1141 20 CONTINUE 1142 1143 DO 100 I = 1, N 1144 DO 50 K = AI(I), AI(I+1)-1 1145 DO 40 J = 1, M 1146 W(I,J) = W(I,J) + AN(K)*P(AJ(K),J) 1147 W(AJ(K),J) = W(AJ(K),J) + AN(K)*P(I,J) 1148 40 CONTINUE 1149 50 CONTINUE 1150 100 CONTINUE 1151 1152 DO 200 I = 1, M 1153 DO 180 J = I, M 1154 C(I,J) = 0 1155 DO 160 K = 1, N 1156 C(I,J) = C(I,J) + P(K,I)*W(K,J) 1157 160 CONTINUE 1158 C(J,I) = C(I,J) 1159 180 CONTINUE 1160 200 CONTINUE 1161 1162 RETURN 1163 END 1164 1165*********************************************************************** 1166 1167 SUBROUTINE DSSTRF( N, M, AI, AJ, AN, AD, P, C, W ) 1168 INTEGER N, AI(*), AJ(*) 1169 DOUBLE PRECISION AN(*), AD(*), P(N,M), C(M,M), W(N,M) 1170 1171* This is the "transform" routine. It multiplies P'*A*P to form the 1172* result C, where P and C are dense and A is a REAL*8, sparse, 1173* symmetric matrix in upper, unordered form. 1174 1175* Input: N The order of A. 1176* M The number of columns of P. 1177* AI,AJ,AN,AD The N by N matrix A. 1178* P The N by M matrix P. 1179* W A work matrix, N by M. 1180 1181* Output: C The M by M result. 1182 1183 DO 20 I = 1, N 1184 DO 10 J = 1, M 1185 W(I,J) = AD(I)*P(I,J) 1186 10 CONTINUE 1187 20 CONTINUE 1188 1189 DO 100 I = 1, N 1190 DO 50 K = AI(I), AI(I+1)-1 1191 DO 40 J = 1, M 1192 W(I,J) = W(I,J) + AN(K)*P(AJ(K),J) 1193 W(AJ(K),J) = W(AJ(K),J) + AN(K)*P(I,J) 1194 40 CONTINUE 1195 50 CONTINUE 1196 100 CONTINUE 1197 1198 DO 200 I = 1, M 1199 DO 180 J = I, M 1200 C(I,J) = 0.0 1201 DO 160 K = 1, N 1202 C(I,J) = C(I,J) + P(K,I)*W(K,J) 1203 160 CONTINUE 1204 C(J,I) = C(I,J) 1205 180 CONTINUE 1206 200 CONTINUE 1207 1208 RETURN 1209 END 1210 1211*********************************************************************** 1212 1213 SUBROUTINE ZHSTRF( N, M, AI, AJ, AN, AD, P, C, W ) 1214 INTEGER N, AI(*), AJ(*) 1215 COMPLEX*16 AN(*), AD(*), P(N,M), C(M,M), W(N,M) 1216 1217* This is the "transform" routine. It multiplies P'*A*P to form the 1218* result C, where P and C are dense and A is a COMPLEX*16, sparse, 1219* hermitian matrix in upper, unordered form. 1220 1221* Input: N The order of A. 1222* M The number of columns of P. 1223* AI,AJ,AN,AD The N by N matrix A. 1224* P The N by M matrix P. 1225* W A work matrix, N by M. 1226 1227* Output: C The M by M result. 1228 1229 DO 20 I = 1, N 1230 DO 10 J = 1, M 1231 W(I,J) = AD(I)*P(I,J) 1232 10 CONTINUE 1233 20 CONTINUE 1234 1235 DO 100 I = 1, N 1236 DO 50 K = AI(I), AI(I+1)-1 1237 DO 40 J = 1, M 1238 W(I,J) = W(I,J) + AN(K)*P(AJ(K),J) 1239 W(AJ(K),J) = W(AJ(K),J) + CONJG(AN(K))*P(I,J) 1240 40 CONTINUE 1241 50 CONTINUE 1242 100 CONTINUE 1243 1244 DO 200 I = 1, M 1245 DO 180 J = I, M 1246 C(I,J) = ( 0.0, 0.0 ) 1247 DO 160 K = 1, N 1248 C(I,J) = C(I,J) + CONJG(P(K,I))*W(K,J) 1249 160 CONTINUE 1250 C(J,I) = CONJG(C(I,J)) 1251 180 CONTINUE 1252 C(I,I) = DBLE(C(I,I)) 1253 200 CONTINUE 1254 1255 RETURN 1256 END 1257