C ALGORITHM 601 COLLECTED ALGORITHMS FROM ACM. C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 3, C SEP., 1983, P. 344-345. C DRIVER PROGRAM TO TEST SPARSE MATRIX PACKAGE FOR MAN 10 C SPECIAL CASES...J.M.MCNAMEE...NOV 1981. MAN 20 C NAMELY :- MAN 30 C MUSFMX... MULTIPLY SPARSE MATRIX BY FULL MATRIX MAN 40 C MUFSMX... MULTIPLY FULL MATRIX BY SPARSE MATRIX MAN 50 C ADSFMX... ADD SPARSE MATRIX TO FULL MATRIX MAN 60 C MUSMFV... MULTIPLY SPARSE MATRIX BY FULL VECTOR MAN 70 C MUFVSM... MULTIPLY FULL VECTOR BY SPARSE MATRIX MAN 80 C MUFMSV... MULTIPLY FULL MATRIX BY SPARSE VECTOR MAN 90 C MUSVFM... MULTIPLY SPARSE VECTOR BY FULL MATRIX MAN 100 C SPEMSV... MULTIPLY SPARSE ELEMENTARY MATRIX BY SPARSE VECTOR MAN 110 C MSMGUS... MULTIPLY TWO SPARSE MATRICES BY GUSTAFSON'S FAST METHOD MAN 120 C TRSMGU... TRANSPOSE SPARSE MATRIX. MAN 130 C FOR DESCRIPTION OF INPUT AND OUTPUT SEE COMMENTS AT START MAN 140 C OF EACH SUBROUTINE. MAN 150 C USES ANCILLARY ROUTINES AS FOLLOWS:- MAN 160 C READSM...READ SPARSE MATRIX MAN 170 C READFM...READ FULL MATRIX MAN 180 C PRNTFM...PRINT FULL MATRIX MAN 190 C READFV...READ FULL VECTOR MAN 200 C PRNTFV...PRINT FULL VECTOR MAN 210 C READSV...READ SPARSE VECTOR MAN 220 C READEM... READ ELEMENTARY MATRIX MAN 230 C PRNTSV...PRINT SPARSE VECTOR MAN 240 C PRNTSM... PRINT SPARSE MATRIX MAN 250 C CNSFMX... CONVERT SPARSE TO FULL MATRIX MAN 260 C MUFUMX... MULTIPLY TWO FULL MATRICES MAN 270 C CMFUMX... COMPARE TWO FULL MATRICES MAN 280 C ADFUMX... ADD TWO FULL MATRICES MAN 290 C MUFMXV... MULTIPLY FULL MATRIX BY FULL VECTOR MAN 300 C CMFUVC... COMPARE TWO FULL VECTORS MAN 310 C MUFVMX... MULTIPLY FULL VECTOR BY FULL MATRIX MAN 320 C CNSFVC... CONVERT SPARSE VECTOR TO FULL VECTOR MAN 330 C CNEFMX... CONVERT SPARSE ELEMENTARY MATRIX TO FULL MATRIX MAN 340 C MAN 350 REAL A(20), B(10,10), C(10,10), X(10), Y(10), D(20), E(20), ET(20)MAN 360 REAL F(10,10), G(10,10), H(10,10) MAN 370 C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- MAN 380 C INTEGER*2 MA(20),MX(10),MY(10) ,MD(20),ME(30),IA(10),ID(10), MAN 390 C * IE(10),MET(20),IET(10) MAN 400 INTEGER MA(20), MX(10), MY(10), MD(20), ME(30), IA(10), ID(10), MAN 410 * IE(10), MET(20), IET(10) MAN 420 C SET DIMENSIONAL PARAMETERS. MAN 430 DATA NA, NM, MRB, MCB, MRC, MCC, NX, NY /20,20,10,10,10,10,10,10/ MAN 440 DATA NDE, NDME /20,30/ MAN 450 DATA MRF, MCF, MRG, MCG, MRH, MCH, ND /10,10,10,10,10,10,20/ MAN 460 DATA IR, LP /5,6/ MAN 470 C TEST ON MUSFMX...MULTIPLY SPARSE MATRIX IN A,MA ON RIGHT MAN 480 C BY FULL MATRIX IN B. MAN 490 DO 30 I=1,3 MAN 500 WRITE (LP,99999) MAN 510 99999 FORMAT (1H0, 70X, 40H***MULTIPLY SPARSE MATRIX BY FULL MATRIX) MAN 520 CALL READSM(A, NA, MA, NM, 1) MAN 530 CALL READFM(B, NRB, NCB, MRB, MCB, 2) MAN 540 CALL MUSFMX(A, NA, MA, NM, B, NRB, NCB, MRB, MCB, C, NRC, NCC, MAN 550 * MRC, MCC, IERROR) MAN 560 IF (IERROR.EQ.0) GO TO 10 MAN 570 WRITE (LP,99998) MAN 580 99998 FORMAT (1H , 70X, 30HA AND B INCOMPATIBLE IN MUSFMX) MAN 590 GO TO 30 MAN 600 10 CALL PRNTFM(C, NRC, NCC, MRC, MCC) MAN 610 CALL CNSFMX(A, NA, MA, NM, F, NRF, NCF, MRF, MCF) MAN 620 CALL MUFUMX(F, NRF, NCF, MRF, MCF, B, NRB, NCB, MRB, MCB, G, MAN 630 * NRG, NCG, MRG, MCG) MAN 640 CALL CMFUMX(G, NRG, NCG, MRG, MCG, C, NRC, NCC, MRC, MCC, MAN 650 * IYESNO) MAN 660 IF (IYESNO.EQ.0) GO TO 20 MAN 670 WRITE (LP,99997) I MAN 680 99997 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON MUSFMX UNSUCCE,MAN 690 * 8HSSFUL***) MAN 700 GO TO 30 MAN 710 20 WRITE (LP,99996) I MAN 720 99996 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON MUSFMX SUCCESSFUL)MAN 730 30 CONTINUE MAN 740 C TEST ON MUFSMX...MULTIPLY FULL MATRIX IN B ON RIGHT BY MAN 750 C SPARSE MATRIX IN A,MA. MAN 760 DO 60 I=1,4 MAN 770 WRITE (LP,99995) MAN 780 99995 FORMAT (1H0, 70X, 40H***MULTIPLY FULL MATRIX BY SPARSE MATRIX) MAN 790 CALL READFM(B, NRB, NCB, MRB, MCB, 1) MAN 800 CALL READSM(A, NA, MA, NM, 2) MAN 810 CALL MUFSMX(B, NRB, NCB, MRB, MCB, A, NA, MA, NM, C, NRC, NCC, MAN 820 * MRC, MCC, IERROR) MAN 830 IF (IERROR.EQ.0) GO TO 40 MAN 840 WRITE (LP,99994) MAN 850 99994 FORMAT (1H , 70X, 30HA AND B INCOMPATIBLE IN MUFSMX) MAN 860 GO TO 60 MAN 870 40 CALL PRNTFM(C, NRC, NCC, MRC, MCC) MAN 880 CALL CNSFMX(A, NA, MA, NM, F, NRF, NCF, MRF, MCF) MAN 890 CALL MUFUMX(B, NRB, NCB, MRB, MCB, F, NRF, NCF, MRF, MCF, G, MAN 900 * NRG, NCG, MRG, MCG) MAN 910 CALL CMFUMX(G, NRG, NCG, MRG, MCG, C, NRC, NCC, MRC, MCC, MAN 920 * IYESNO) MAN 930 IF (IYESNO.EQ.0) GO TO 50 MAN 940 WRITE (LP,99993) I MAN 950 99993 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON MUFSMX UNSUCCE,MAN 960 * 8HSSFUL***) MAN 970 GO TO 60 MAN 980 50 WRITE (LP,99992) I MAN 990 99992 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON MUFSMX SUCCESSFUL)MAN 1000 60 CONTINUE MAN 1010 C TEST ON ADSFMX...ADD SPARSE MATRIX IN A,MA TO FULL MATRIX MAN 1020 C IN B,LEAVING RESULT IN B. MAN 1030 DO 110 I=1,2 MAN 1040 WRITE (LP,99991) MAN 1050 99991 FORMAT (1H0, 70X, 28H***ADD SPARSE TO FULL MATRIX) MAN 1060 CALL READSM(A, NA, MA, NM, 1) MAN 1070 CALL READFM(B, NRB, NCB, MRB, MCB, 2) MAN 1080 DO 80 K=1,NRB MAN 1090 DO 70 J=1,NCB MAN 1100 H(K,J) = B(K,J) MAN 1110 70 CONTINUE MAN 1120 80 CONTINUE MAN 1130 CALL ADSFMX(A, NA, MA, NM, B, NRB, NCB, MRB, MCB, IERROR) MAN 1140 IF (IERROR.EQ.0) GO TO 90 MAN 1150 WRITE (LP,99990) MAN 1160 99990 FORMAT (1H , 70X, 30HA AND B INCOMPATIBLE IN ADSFMX) MAN 1170 GO TO 110 MAN 1180 90 CALL PRNTFM(B, NRB, NCB, MRB, MCB) MAN 1190 CALL CNSFMX(A, NA, MA, NM, F, NRF, NCF, MRF, MCF) MAN 1200 CALL ADFUMX(H, NRB, NCB, MRH, MCH, F, NRF, NCF, MRF, MCF, G, MAN 1210 * NRG, NCG, MRG, MCG) MAN 1220 CALL CMFUMX(G, NRG, NCG, MRG, MCG, B, NRB, NCB, MRB, MCB, MAN 1230 * IYESNO) MAN 1240 IF (IYESNO.EQ.0) GO TO 100 MAN 1250 WRITE (LP,99989) I MAN 1260 99989 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON ADSFMX UNSUCCE,MAN 1270 * 8HSSFUL***) MAN 1280 GO TO 110 MAN 1290 100 WRITE (LP,99988) I MAN 1300 99988 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON ADSFMX SUCCESSFUL)MAN 1310 110 CONTINUE MAN 1320 C TEST ON MUSMFV...MULTIPLY SPARSE MATRIX IN A,MA ON RIGHT MAN 1330 C BY FULL VECTOR IN X. MAN 1340 DO 140 I=1,2 MAN 1350 WRITE (LP,99987) MAN 1360 99987 FORMAT (1H0, 70X, 40H***MULTIPLY SPARSE MATRIX BY FULL VECTOR) MAN 1370 CALL READSM(A, NA, MA, NM, 1) MAN 1380 CALL READFV(X, NRX, NX, 2) MAN 1390 CALL MUSMFV(A, NA, MA, NM, X, NRX, NX, Y, NRY, NY, IERROR) MAN 1400 IF (IERROR.EQ.0) GO TO 120 MAN 1410 WRITE (LP,99986) MAN 1420 99986 FORMAT (1H , 70X, 30HA AND X INCOMPATIBLE IN MUSMFV) MAN 1430 GO TO 140 MAN 1440 120 CALL PRNTFV(Y, NRY, NY) MAN 1450 CALL CNSFMX(A, NA, MA, NM, F, NRF, NCF, MRF, MCF) MAN 1460 CALL MUFMXV(F, NRF, NCF, MRF, MCF, X, NRX, NX, D, NRD, ND) MAN 1470 CALL CMFUVC(D, NRD, ND, Y, NRY, NY, IYESNO) MAN 1480 IF (IYESNO.EQ.0) GO TO 130 MAN 1490 WRITE (LP,99985) I MAN 1500 99985 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON MUSMFV UNSUCCE,MAN 1510 * 8HSSFUL***) MAN 1520 GO TO 140 MAN 1530 130 WRITE (LP,99984) I MAN 1540 99984 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON MUSMFV SUCCESSFUL)MAN 1550 140 CONTINUE MAN 1560 C TEST ON MUFVSM...MULTIPLY FULL VECTOR IN X ON RIGHT BY MAN 1570 C SPARSE MATRIX IN A,MA. MAN 1580 DO 170 I=1,2 MAN 1590 WRITE (LP,99983) MAN 1600 99983 FORMAT (1H0, 70X, 40H***MULTIPLY FULL VECTOR BY SPARSE MATRIX) MAN 1610 CALL READFV(X, NCX, NX, 1) MAN 1620 CALL READSM(A, NA, MA, NM, 2) MAN 1630 CALL MUFVSM(X, NCX, NX, A, NA, MA, NM, Y, NCY, NY, IERROR) MAN 1640 IF (IERROR.EQ.0) GO TO 150 MAN 1650 WRITE (LP,99982) MAN 1660 99982 FORMAT (1H , 70X, 30HX AND A INCOMPATIBLE IN MUFVSM) MAN 1670 GO TO 170 MAN 1680 150 CALL PRNTFV(Y, NCY, NY) MAN 1690 CALL CNSFMX(A, NA, MA, NM, F, NRF, NCF, MRF, MCF) MAN 1700 CALL MUFVMX(X, NCX, NX, F, NRF, NCF, MRF, MCF, D, NCD, ND) MAN 1710 CALL CMFUVC(D, NCD, ND, Y, NCY, NY, IYESNO) MAN 1720 IF (IYESNO.EQ.0) GO TO 160 MAN 1730 WRITE (LP,99981) I MAN 1740 99981 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON MUFVSM UNSUCCE,MAN 1750 * 8HSSFUL***) MAN 1760 GO TO 170 MAN 1770 160 WRITE (LP,99980) I MAN 1780 99980 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON MUFVSM SUCCESSFUL)MAN 1790 170 CONTINUE MAN 1800 C TEST ON MUFMSV..MULTIPLY FULL MATRIX IN B ON RIGHT BY MAN 1810 C SPARSE VECTOR IN X,MX. MAN 1820 DO 200 I=1,3 MAN 1830 WRITE (LP,99979) MAN 1840 99979 FORMAT (1H0, 70X, 40H***MULTIPLY FULL MATRIX BY SPARSE VECTOR) MAN 1850 CALL READFM(B, NRB, NCB, MRB, MCB, 1) MAN 1860 CALL READSV(X, NRX, NEX, NX, MX, 2) MAN 1870 CALL MUFMSV(B, NRB, NCB, MRB, MCB, X, NRX, NEX, NX, MX, Y, NRY, MAN 1880 * NY, IERROR) MAN 1890 IF (IERROR.EQ.0) GO TO 180 MAN 1900 WRITE (LP,99978) MAN 1910 99978 FORMAT (1H , 70X, 30HB AND X INCOMPATIBLE IN MUFMSV) MAN 1920 GO TO 200 MAN 1930 180 CALL PRNTFV(Y, NRY, NY) MAN 1940 CALL CNSFVC(X, NRX, NEX, NX, MX, D, NRD, ND) MAN 1950 CALL MUFMXV(B, NRB, NCB, MRB, MCB, D, NRD, ND, E, NRE, NDE) MAN 1960 CALL CMFUVC(Y, NRY, NY, E, NRE, NDE, IYESNO) MAN 1970 IF (IYESNO.EQ.0) GO TO 190 MAN 1980 WRITE (LP,99977) I MAN 1990 99977 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON MUFMSV UNSUCCE,MAN 2000 * 8HSSFUL***) MAN 2010 GO TO 200 MAN 2020 190 WRITE (LP,99976) I MAN 2030 99976 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON MUFMSV SUCCESSFUL)MAN 2040 200 CONTINUE MAN 2050 C TEST ON MUSVFM...MULTIPLY SPARSE VECTOR IN X,MX BY FULL MAN 2060 C MATRIX IN B. MAN 2070 DO 230 I=1,3 MAN 2080 WRITE (LP,99975) MAN 2090 99975 FORMAT (1H0, 70X, 40H***MULTIPLY SPARSE VECTOR BY FULL MATRIX) MAN 2100 CALL READSV(X, NCX, NEX, NX, MX, 1) MAN 2110 CALL READFM(B, NRB, NCB, MRB, MCB, 2) MAN 2120 CALL MUSVFM(X, NCX, NEX, NX, MX, B, NRB, NCB, MRB, MCB, Y, NCY, MAN 2130 * NY, IERROR) MAN 2140 IF (IERROR.EQ.0) GO TO 210 MAN 2150 WRITE (LP,99974) MAN 2160 99974 FORMAT (1H , 70X, 30HX AND B INCOMPATIBLE IN MUSVFM) MAN 2170 GO TO 230 MAN 2180 210 CALL PRNTFV(Y, NCY, NY) MAN 2190 CALL CNSFVC(X, NCX, NEX, NX, MX, D, NCD, ND) MAN 2200 CALL MUFVMX(D, NCD, ND, B, NRB, NCB, MRB, MCB, E, NCE, NDE) MAN 2210 CALL CMFUVC(Y, NCY, NY, E, NCE, NDE, IYESNO) MAN 2220 IF (IYESNO.EQ.0) GO TO 220 MAN 2230 WRITE (LP,99973) I MAN 2240 99973 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON MUSVFM UNSUCCE,MAN 2250 * 8HSSFUL***) MAN 2260 GO TO 230 MAN 2270 220 WRITE (LP,99972) I MAN 2280 99972 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON MUSVFM SUCCESSFUL)MAN 2290 230 CONTINUE MAN 2300 C TEST ON SPEMSV...MULTIPLY SPARSE ELEMENTARY MATRIX IN A,MA MAN 2310 C ON RIGHT BY SPARSE VECTOR IN X,MX. MAN 2320 DO 260 I=1,4 MAN 2330 WRITE (LP,99971) MAN 2340 99971 FORMAT (1H0, 60X, 39H***MULTIPLY SPARSE ELEMENTARY MATRIX BY, MAN 2350 * 14H SPARSE VECTOR) MAN 2360 C READ THE ELEMENTARY MATRIX. MAN 2370 CALL READEM(A, NRA, NEA, NA, MA, 1) MAN 2380 C READ THE INDEX OF THE NON-UNIT COLUMN. MAN 2390 READ (IR,99970) K MAN 2400 99970 FORMAT (I5) MAN 2410 WRITE (LP,99969) K MAN 2420 99969 FORMAT (28H0INDEX OF NON-UNIT COLUMN = , I5) MAN 2430 CALL READSV(X, NRX, NEX, NX, MX, 2) MAN 2440 CALL SPEMSV(A, NRA, NEA, NA, MA, K, X, NRX, NEX, NX, MX, Y, MAN 2450 * NRY, NEY, NY, MY, IERROR) MAN 2460 IF (IERROR.EQ.0) GO TO 240 MAN 2470 WRITE (LP,99968) MAN 2480 99968 FORMAT (1H , 70X, 30HA AND X INCOMPATIBLE IN SPEMSV) MAN 2490 GO TO 260 MAN 2500 240 CALL PRNTSV(Y, NRY, NEY, NY, MY) MAN 2510 CALL CNEFMX(A, NRA, NEA, NA, MA, K, F, NRF, NCF, MRF, MCF) MAN 2520 CALL CNSFVC(X, NRX, NEX, NX, MX, D, NRD, ND) MAN 2530 CALL MUFMXV(F, NRF, NCF, MRF, MCF, D, NRD, ND, E, NRE, NDE) MAN 2540 CALL CNSFVC(Y, NRY, NEY, NY, MY, ET, NRET, NDE) MAN 2550 CALL CMFUVC(ET, NRET, NDE, E, NRE, NDE, IYESNO) MAN 2560 IF (IYESNO.EQ.0) GO TO 250 MAN 2570 WRITE (LP,99967) I MAN 2580 99967 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON SPEMSV UNSUCCE,MAN 2590 * 8HSSFUL***) MAN 2600 GO TO 260 MAN 2610 250 WRITE (LP,99966) I MAN 2620 99966 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON SPEMSV SUCCESSFUL)MAN 2630 260 CONTINUE MAN 2640 C TEST ON MSMGUS... MULTIPLY SPARSE MATRIX IN A,MA BY SPARSE MAN 2650 C MATRIX IN D,MD. MAN 2660 DO 300 I=1,5 MAN 2670 WRITE (LP,99965) MAN 2680 99965 FORMAT (1H0, 70X, 42H***MULTIPLY SPARSE MATRIX BY SPARSE MATRIX)MAN 2690 CALL READSM(A, NA, MA, NM, 1) MAN 2700 CALL READSM(D, NA, MD, NM, 2) MAN 2710 CALL MSMGUS(A, MA, D, MD, ID, NDE, E, NDME, ME, IE, IERROR, X, MAN 2720 * MY) MAN 2730 IF (IERROR.EQ.0) GO TO 280 MAN 2740 IF (IERROR.EQ.1) GO TO 270 MAN 2750 WRITE (LP,99964) MAN 2760 99964 FORMAT (1H0, 70X, 30HA AND D INCOMPATIBLE IN MSMGUS) MAN 2770 GO TO 300 MAN 2780 270 WRITE (LP,99963) MAN 2790 99963 FORMAT (1H , 75X, 24HSPACE EXCEEDED IN MSMGUS) MAN 2800 GO TO 300 MAN 2810 280 CALL TRSMGU(E, ME, IE, ET, MET, IET) MAN 2820 CALL TRSMGU(ET, MET, IET, E, ME, IE) MAN 2830 CALL PRNTSM(E, ME) MAN 2840 CALL CNSFMX(A, NA, MA, NM, F, NRF, NCF, MRF, MCF) MAN 2850 CALL CNSFMX(D, ND, MD, NM, G, NRG, NCG, MRG, MCG) MAN 2860 CALL CNSFMX(E, NDE, ME, NDME, B, NRB, NCB, MRB, MCB) MAN 2870 CALL MUFUMX(F, NRF, NCF, MRF, MCF, G, NRG, NCG, MRG, MCG, H, MAN 2880 * NRH, NCH, MRH, MCH) MAN 2890 CALL CMFUMX(H, NRH, NCH, MRH, MCH, B, NRB, NCB, MRB, MCB, MAN 2900 * IYESNO) MAN 2910 IF (IYESNO.EQ.0) GO TO 290 MAN 2920 WRITE (LP,99962) I MAN 2930 99962 FORMAT (1H , 70X, 15H***TEST NUMBER , I5, 18H ON MSMGUS UNSUCCE,MAN 2940 * 8HSSFUL***) MAN 2950 GO TO 300 MAN 2960 290 WRITE (LP,99961) I MAN 2970 99961 FORMAT (1H , 70X, 12HTEST NUMBER , I5, 21H ON MSMGUS SUCCESSFUL)MAN 2980 300 CONTINUE MAN 2990 STOP MAN 3000 END MAN 3010 SUBROUTINE MUSFMX(A, NA, MA, NM, B, NRB, NCB, MRB, MCB, C, NRC, MUS 10 * NCC, MRC, MCC, IERROR) C********************************************************** C MULTIPLY SPARSE MATRIX STORED IN A,MA ON RIGHT BY FULL C MATRIX STORED IN B. PUT RESULT IN C. C*** INPUT C A A ONE-DIMENSIONAL ARRAY CONTAINING THE NON- C ZERO ELEMENTS OF FIRST(SPARSE)MATRIX,ARRANGED C ROW-BY-ROW,BUT NOT NECESSARILY IN ORDER OF C INCREASING COLUMN NUMBER. C NA THE DIMENSION OF A IN THE CALLING ROUTINE. C MA(1) NUMBER OF ROWS IN A. LATER CALLED NRA. C MA(2) NUMBER OF COLUMNS IN A. LATER CALLED NCA. C MA(3) BLANK C MA(4) TOTAL NUMBER OF NON-ZERO ELEMENTS IN A. C MA(5) NUMBER OF NON-ZERO ELEMENTS IN ROW 1 OF A. C ...MA(4+I)... NUMBER OF NON-ZERO ELEMENTS IN ROW I OF A. C ... C ..MA(4+NRA+1) COLUMN INDEX OF 1ST NON-ZERO ELEMENT OF A. C ... C MA(4+NRA+I) COLUMN INDEX OF I'TH NON-ZERO ELEMENT OF A. C C NM DIMENSION OF MA IN CALLING ROUTINE. C C C B A TWO-DIMENSIONAL ARRAY CONTAINING ALL THE C ELEMENTS OF THE SECOND MATRIX. C NRB ACTUAL NUMBER OF ROWS IN B. C NCB ACTUAL NUMBER OF COLUMNS IN B. C MRB ROW-DIMENSION OF B IN CALLING ROUTINE. C MCB COLUMN DIMENSION OF B IN CALLING ROUTINE. C C*** OUTPUT C C A TWO-DIMENSIONAL ARRAY CONTAINING ALL THE C ELEMENTS OF THE PRODUCT MATRIX. C NRC ACTUAL NUMBER OF ROWS IN C. C NCC ACTUAL NUMBER OF COLUMNS IN C. C MRC ROW DIMENSION OF C IN CALLING ROUTINE. C MCB COLUMN DIMENSION OF C IN CALLING ROUTINE. C IERROR =0 IF NO ERROR,=1 IF A AND B INCOMPATIBLE C REAL A(NA), B(MRB,MCB), C(MRC,MCC) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(NM) INTEGER MA(NM) C LP IS UNIT NUMBER OF LINE PRINTER. LP = 6 C PICK OUT CONTROL INFORMATION. C NRA,NCA ARE NUMBER OF ROWS,COLUMNS IN A. NRA = MA(1) NCA = MA(2) C NRC,NCC ARE NUMBER OF ROWS,COLUMNS IN C. NRC = NRA NCC = NCB C SET ERROR FLAG IERROR = 0 C CHECK FOR COMPATIBILITY UNDER MULTIPLICATION. IF (NCA.EQ.NRB) GO TO 10 C INCOMPATIBLE.RESET ERROR FLAG AND RETURN. IERROR = 1 RETURN C CLEAR C TO ZERO. 10 DO 30 I=1,NRC DO 20 J=1,NCC C(I,J) = 0.0 20 CONTINUE 30 CONTINUE C N2 WILL BE POINTER TO END OF ROW I OF A. N2 = 0 C I WILL BE ROW-INDEX FOR A. DO 60 I=1,NRA C PICK OUT NUMBER OF NON-ZERO ELEMENTS IN ROW I. NIR = MA(I+4) C IF NO NON-ZEROES SKIP PROCESSING OF ROW I OF A. IF (NIR.EQ.0) GO TO 60 C N1 POINTS TO START OF ROW I IN A,N2 TO END. N1 = N2 + 1 N2 = N2 + NIR C PROCESS ROW I OF A, I.E. FORM ALL PRODUCTS OF NON-ZERO A(I,L) WITH C B(L,J); PUT INTO C(I,J). C K POINTS TO NON-ZERO ELEMENTS IN ROW I OF A. DO 50 K=N1,N2 NRAK4 = NRA + K + 4 L = MA(NRAK4) T = A(K) DO 40 J=1,NCB C(I,J) = C(I,J) + T*B(L,J) 40 CONTINUE 50 CONTINUE 60 CONTINUE RETURN END SUBROUTINE MUFSMX(B, NRB, NCB, MRB, MCB, A, NA, MA, NM, C, NRC, MUF 10 * NCC, MRC, MCC, IERROR) C MULTIPLY A FULL MATRIX STORED IN B ON RIGHT BY A SPARSE C MATRIX STORED IN A,MA.PUT RESULT IN C. C*** INPUT C B A TWO-DIMENSIONAL ARRAY CONTAINING ALL THE C ELEMENTS OF THE FIRST MATRIX. C NRB ACTUAL NUMBER OF ROWS IN B. C NCB ACTUAL NUMBER OF COLUMNS IN B. C MRB ROW-DIMENSION OF B IN CALLING ROUTINE. C MCB COLUMN DIMENSION OF B IN CALLING ROUTINE. C C A A ONE-DIMENSIONAL ARRAY CONTAINING THE NON- C ZERO ELEMENTS OF SECOND (SPARSE) MATRIX,ARRANGED C ROW-BY-ROW,BUT NOT NECESSARILY IN ORDER OF C INCREASING COLUMN NUMBER. C NA THE DIMENSION OF A IN THE CALLING ROUTINE. C MA(1) NUMBER OF ROWS IN A. LATER CALLED NRA. C MA(2) NUMBER OF COLUMNS IN A. LATER CALLED NCA C MA(3) BLANK C MA(4) TOTAL NUMBER OF NON-ZERO ELEMENTS IN A. C MA(5) NUMBER OF NON-ZERO ELEMENTS IN ROW 1 OF A. C ...MA(4+I)... NUMBER OF NON-ZERO ELEMENTS IN ROW I OF A. C ... C ..MA(4+NRA+1) COLUMN INDEX OF 1ST NON-ZERO ELEMENT OF A. C ... C MA(4+NRA+I) COLUMN INDEX OF I'TH NON-ZERO ELEMENT OF A C C NM DIMENSION OF MA IN CALLING ROUTINE. C***OUTPUT C C A TWO-DIMENSIONAL ARRAY CONTAINING ALL THE C ELEMENTS OF THE PRODUCT MATRIX. C NRC ACTUAL NUMBER OF ROWS IN C. C NCC ACTUAL NUMBER OF COLUMNS IN C. C MRC ROW DIMENSION OF C IN CALLING ROUTINE. C MCB COLUMN DIMENSION OF C IN CALLING ROUTINE. C IERROR =1 IF A AND B INCOMPATIBLE,=0 IF NO ERROR. C REAL A(NA), B(MRB,MCB), C(MRC,MCC) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(NM) INTEGER MA(NM) C LP IS UNIT NUMBER OF LINE PRINTER. LP = 6 C PICK OUT CONTROL INFORMATION. C NRA,NCA ARE NUMBER OF ROWS,COLUMNS IN A. NRA = MA(1) NCA = MA(2) C NRC,NCC ARE NUMBER OF ROWS,COLUMNS IN C. NRC = NRB NCC = NCA C SET ERROR FLAG IERROR = 0 C CHECK FOR COMPATIBILITY UNDER MULTIPLICATION. IF (NCB.EQ.NRA) GO TO 10 C INCOMPATIBLE.RESET ERROR FLAG AND RETURN. IERROR = 1 RETURN C CLEAR C TO ZERO. 10 DO 30 I=1,NRC DO 20 J=1,NCC C(I,J) = 0. 20 CONTINUE 30 CONTINUE C N2 WILL BE POINTER TO END OF ROW K OF A. N2 = 0 C K WILL BE ROW-INDEX FOR A. DO 60 K=1,NRA C PICK OUT NUMBER OF NON-ZEROS IN ROW K OF A. K4 = K + 4 NIR = MA(K4) IF (NIR.EQ.0) GO TO 60 C IF NO NON-ZEROS SKIP PROCESSING OF ROW K. C N1 WILL POINT TO START OF ROW K IN A,N2 TO END OF ROW. N1 = N2 + 1 N2 = N2 + NIR C PROCESS ROW K OF A,I.E. C FORM ALL PRODUCTS OF B(I,K) WITH NON-ZERO A(K,J), C ADD INTO C(I,J) C L POINTS TO NON-ZERO ELEMENTS IN ROW K OF A. DO 50 L=N1,N2 NRAL4 = NRA + L + 4 J = MA(NRAL4) DO 40 I=1,NRB C(I,J) = C(I,J) + B(I,K)*A(L) 40 CONTINUE 50 CONTINUE 60 CONTINUE RETURN END SUBROUTINE ADSFMX(A, NA, MA, NM, B, NRB, NCB, MRB, MCB, IERROR) ADS 10 C ADD SPARSE MATRIX IN A,MA TO FULL MATRIX IN B. C SUM OVERWRITES B. C*** INPUT C A A ONE-DIMENSIONAL ARRAY CONTAINING THE NON- C ZERO ELEMENTS OF FIRST(SPARSE)MATRIX,ARRANGED C ROW-BY-ROW,BUT NOT NECESSARILY IN ORDER OF C INCREASING COLUMN NUMBER. C NA THE DIMENSION OF A IN THE CALLING ROUTINE. C MA(1) NUMBER OF ROWS IN A. LATER CALLED NRA. C MA(2) NUMBER OF COLUMNS IN A. LATER CALLED NCA. C MA(3) BLANK C MA(4) TOTAL NUMBER OF NON-ZERO ELEMENTS IN A. C MA(5) NUMBER OF NON-ZERO ELEMENTS IN ROW 1 OF A. C ...MA(4+I)... NUMBER OF NON-ZERO ELEMENTS IN ROW I OF A. C ... C ..MA(4+NRA+1) COLUMN INDEX OF 1ST NON-ZERO ELEMENT OF A. C ... C MA(4+NRA+I) COLUMN INDEX OF I'TH NON-ZERO ELEMENT OF A. C C NM DIMENSION OF MA IN CALLING ROUTINE. C C B A TWO-DIMENSIONAL ARRAY CONTAINING ALL THE C ELEMENTS OF THE SECOND MATRIX. C NRB ACTUAL NUMBER OF ROWS IN B. C NCB ACTUAL NUMBER OF COLUMNS IN B. C MRB ROW-DIMENSION OF B IN CALLING ROUTINE. C MCB COLUMN DIMENSION OF B IN CALLING ROUTINE. C C*** OUTPUT C B THE OUTPUT (SUM OF A AND B) OVERWRITES THE FULL MATRIX B. C IERROR =1 IF A AND B INCOMPATIBLE,=0 IF NO ERROR. C REAL A(NA), B(MRB,MCB) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(NM) INTEGER MA(NM) C LP IS UNIT NUMBER OF LINE PRINTER. LP = 6 C PICK OUT CONTROL INFORMATION. C NRA,NCA ARE NUMBER OF ROWS,COLUMNS IN A. NRA = MA(1) NCA = MA(2) C SET ERROR FLAG IERROR = 0 C CHECK COMPATIBILITY FOR ADDITION . IF (NRA.EQ.NRB .AND. NCA.EQ.NCB) GO TO 10 C INCOMPATIBLE.RESET ERROR FLAG AND RETURN. IERROR = 1 RETURN 10 CONTINUE C JR IS NUMBER OF CONTROL ITEMS IN A. JR = 4 + NRA C N1,N2 ARE POINTERS TO START AND END OF ROW I OF A. N2 = 0 C I COUNTS ROWS OF A. DO 30 I=1,NRA C NIR IS NUMBER OF NON-ZERO ELEMENTS IN ROW I OF A. I4 = I + 4 NIR = MA(I4) C IF NO NON-ZEROS IN ROW I SKIP PROCESSING OF ROW I. IF (NIR.EQ.0) GO TO 30 N1 = N2 + 1 N2 = N2 + NIR DO 20 K=N1,N2 JRK = JR + K J = MA(JRK) B(I,J) = B(I,J) + A(K) 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE MUSMFV(A, NA, M, NM, X, NRX, NX, Y, NRY, NY, IERROR) MUS 10 C TO MULTIPLY A SPARSE MATRIX IN A,MA ON THE RIGHT BY A FULL C VECTOR IN X. PUT THE RESULT IN Y. C*** INPUT C A A ONE-DIMENSIONAL ARRAY CONTAINING THE NON- C ZERO ELEMENTS OF FIRST(SPARSE)MATRIX,ARRANGED C ROW-BY-ROW,BUT NOT NECESSARILY IN ORDER OF C INCREASING COLUMN NUMBER. C NA THE DIMENSION OF A IN THE CALLING ROUTINE. C M(1) NUMBER OF ROWS IN A. LATER CALLED NR. C M(2) NUMBER OF COLUMNS IN A. LATER CALLED NC. C M(3) BLANK C M(4) TOTAL NUMBER OF NON-ZERO ELEMENTS IN A. C M(5) NUMBER OF NON-ZERO ELEMENTS IN ROW 1 OF A. C ...M(4+I) ... NUMBER OF NON-ZERO ELEMENTS IN ROW I OF A. C ... C .. M(4+NR +1) COLUMN INDEX OF 1ST NON-ZERO ELEMENT OF A. C ... C M(4+NR +I) COLUMN INDEX OF I'TH NON-ZERO ELEMENT OF A. C C NM DIMENSION OF M IN CALLING ROUTINE. C C X THE ELEMENTS OF THE FULL VECTOR IN A ONE- C DIMENSIONAL ARRAY. C NRX ACTUAL NUMBER OF ELEMENTS IN X. C NX,NY DIMENSION OF X,Y IN CALLING ROUTINE. C*** OUTPUT C Y ELEMENTS OF THE PRODUCT VECTOR. C NRY ACTUAL NUMBER OF ELEMENTS (ROWS) IN Y. C IERROR =1 IF A AND X INCOMPATIBLE,=0 IF NO ERROR. C REAL A(NA), X(NX), Y(NY) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 M(NM) INTEGER M(NM) C EXTRACT CONTROL INFORMATION. C NR,NC IS NUMBER OF ROWS,COLUMNS IN A. NR = M(1) NC = M(2) NRY = NR C JR IS NUMBER OF CONTROL ELEMENTS IN M. JR = 4 + NR C SET ERROR FLAG IERROR = 0 C CHECK FOR COMPATIBILITY BETWEEN A,X. IF (NC.EQ.NRX) GO TO 10 C INCOMPATIBLE.RESET ERROR FLAG AND RETURN. IERROR = 1 RETURN C N1,N2 ARE POINTERS TO START AND END OF ROW I OF A. 10 N2 = 0 C I COUNTS ROWS OF A. DO 40 I=1,NR C FORM INNER PRODUCT OF ROW I OF A WITH VECTOR X. S = 0. I4 = I + 4 NIR = M(I4) C IF NO NON-ZEROS IN ROW I OF A SKIP THE PRODUCT IF (NIR.EQ.0) GO TO 30 N1 = N2 + 1 N2 = N2 + NIR DO 20 K=N1,N2 JRK = JR + K J = M(JRK) S = S + A(K)*X(J) 20 CONTINUE 30 Y(I) = S 40 CONTINUE RETURN END SUBROUTINE MUFVSM(X, NCX, NX, A, NA, M, NM, Y, NCY, NY, IERROR) MUF 10 C MULTIPLY A FULL VECTOR IN X BY A SPARSE MATRIX STORED IN A,M. C PUT RESULT IN FULL VECTOR Y. C*** INPUT C X THE ELEMENTS OF THE FULL VECTOR IN A ONE- C DIMENSIONAL ARRAY. C NCX ACTUAL NUMBER OF ELEMENTS IN X. C NX,NY DIMENSION OF X,Y IN CALLING ROUTINE. C A A ONE-DIMENSIONAL ARRAY CONTAINING THE NON- C ZERO ELEMENTS OF SECOND (SPARSE) MATRIX A,ARRANGED C ROW-BY-ROW,BUT NOT NECESSARILY IN ORDER OF C INCREASING COLUMN NUMBER. C NA THE DIMENSION OF A IN THE CALLING ROUTINE. C M(1) NUMBER OF ROWS IN A. LATER CALLED NR. C M(2) NUMBER OF COLUMNS IN A. LATER CALLED NC. C M(3) BLANK C M(4) TOTAL NUMBER OF NON-ZERO ELEMENTS IN A. C M(5) NUMBER OF NON-ZERO ELEMENTS IN ROW 1 OF A. C ...M(4+I) ... NUMBER OF NON-ZERO ELEMENTS IN ROW I OF A. C ... C .. M(4+NR +1) COLUMN INDEX OF 1ST NON-ZERO ELEMENT OF A. C ... C M(4+NR +I) COLUMN INDEX OF I'TH NON-ZERO ELEMENT OF A. C C NM DIMENSION OF M IN CALLING ROUTINE. C C*** OUTPUT C Y ELEMENTS OF THE PRODUCT VECTOR. C NCY ACTUAL NUMBER OF ELEMENTS (COLS) IN Y. C IERROR =1 IF X AND A INCOMPATIBLE,=0 IF NO ERROR. C REAL A(NA), X(NX), Y(NY) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 M(NM) INTEGER M(NM) C LP IS UNIT NUMBER OF LINE PRINTER LP = 6 C EXTRACT CONTROL INFORMATION. C NR,NC IS NUMBER OF ROWS,COLUMNS IN A. NR = M(1) NC = M(2) NCY = NC C JR IS NUMBER OF CONTROL ELEMENTS IN M. JR = 4 + NR C SET ERROR FLAG IERROR = 0 C CHECK FOR COMPATIBILITY BETWEEN A,X. IF (NCX.EQ.NR) GO TO 10 C INCOMPATIBLE.RESET ERROR FLAG AND RETURN. IERROR = 1 RETURN C N1,N2 ARE POINTERS TO START AND END OF ROW I OF A. 10 N2 = 0 DO 20 I=1,NC Y(I) = 0. 20 CONTINUE C I COUNTS ROWS OF A. DO 40 I=1,NR I4 = I + 4 NIR = M(I4) IF (NIR.EQ.0) GO TO 40 T = X(I) N1 = N2 + 1 N2 = N2 + NIR DO 30 K=N1,N2 JRK = JR + K J = M(JRK) Y(J) = Y(J) + T*A(K) 30 CONTINUE 40 CONTINUE RETURN END SUBROUTINE MUFMSV(A, NRA, NCA, MRA, MCA, X, NRX, NEX, NX, MX, Y, MUF 10 * NRY, NY, IERROR) C MULTIPLY A FULL MATRIX A ON RIGHT BY A SPARSE COLUMN C VECTOR IN X,M. PUT RESULT IN FULL VECTOR Y. C*** INPUT C A A TWO-DIMENSIONAL ARRAY CONTAINING THE FULL MATRIX. C NRA ACTUAL NUMBER OF ROWS IN A. C NCA ACTUAL NUMBER OF COLS IN A. . C MRA ROW-DIMENSION OF A IN CALLING ROUTINE. C MCA COLUMN-DIMENSION OF A IN CALLING ROUTINE. C C X A ONE-DIMENSIONAL ARRAY CONTAINING THE NON-ZERO C ELEMENTS OF THE SPARSE VECTOR,NOT NECESSARILY C IN ORDER OF INCREASING ROW NUMBER. C NRX NUMBER OF ROWS IN X (INCLUDING ZEROS). C NEX NUMBER OF NON-ZERO ELEMENTS IN X. C NX,NY DIMENSIONS OF X,Y IN CALLING ROUTINE. C MX M(I) CONTAINS ROW-INDEX OF X(I) C*** OUTPUT C Y A ONE-DIMENSIONAL ARRAY CONTAINING ALL THE C ELEMENTS OF THE PRODUCT VECTOR. C NRY ACTUAL NUMBER OF ELEMENTS IN Y. C IERROR =0 IF NO ERRORS,=1 IF A AND X INCOMPATIBLE. C REAL A(MRA,MCA), X(NX), Y(NY) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MX(NX) INTEGER MX(NX) C SET ERROR FLAG. IERROR = 0 C CHECK COMPATIBILITY OF A,X FOR MULTIPLICATION. IF (NCA.EQ.NRX) GO TO 10 C INCOMPATIBLE.RESET ERROR FLAG. IERROR = 1 RETURN 10 NRY = NRA C CLEAR Y TO ZERO. DO 20 I=1,NRA Y(I) = 0.0 20 CONTINUE C IF X IDENTICALLY ZERO SKIP PRODUCT. IF (NEX.EQ.0) GO TO 50 C J COUNTS ELEMENTS IN X. DO 40 J=1,NEX C PICK OUT ROW INDEX OF J'TH ELEMENT OF X. K = MX(J) T = X(J) C MULTIPLY EACH ELEMENT OF COLUMN K OF A BY X(K) AND ADD TO Y(I) DO 30 I=1,NRA Y(I) = Y(I) + A(I,K)*T 30 CONTINUE 40 CONTINUE 50 RETURN END SUBROUTINE MUSVFM(X, NCX, NEX, NX, MX, A, NRA, NCA, MRA, MCA, Y, MUS 10 * NCY, NY, IERROR) C MULTIPLY A SPARSE ROW VECTOR IN X,M BY A FULL MATRIX IN A. C PUT RESULT IN FULL VECTOR Y. C*** INPUT C X A ONE-DIMENSIONAL ARRAY CONTAINING THE NON-ZERO C ELEMENTS OF THE SPARSE VECTOR,NOT NECESSARILY C IN ORDER OF INCREASING COLUMN NUMBER. C NCX NUMBER OF COLS IN X (INCLUDING ZEROS). C NEX NUMBER OF NON-ZERO ELEMENTS IN X. C NX,NY DIMENSIONS OF X,Y IN CALLING ROUTINE. C MX M(I) CONTAINS ROW-INDEX OF X(I) C C A A TWO-DIMENSIONAL ARRAY CONTAINING THE FULL MATRIX. C NRA ACTUAL NUMBER OF ROWS IN A. C NCA ACTUAL NUMBER OF COLS IN A. C MRA ROW-DIMENSION OF A IN CALLING ROUTINE. C MCA COLUMN-DIMENSION OF A IN CALLING ROUTINE. C*** OUTPUT C Y A ONE-DIMENSIONAL ARRAY CONTAINING ALL THE C ELEMENTS OF THE PRODUCT VECTOR. C NCY ACTUAL NUMBER OF ELEMENTS IN Y. C IERROR =0 IF NO ERRORS,=1 IF A AND X INCOMPATIBLE. C REAL A(MRA,MCA), X(NX), Y(NY) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MX(NX) INTEGER MX(NX) C SET ERROR FLAG. IERROR = 0 C CHECK COMPATIBILITY OF A,X FOR MULTIPLICATION. IF (NCX.EQ.NRA) GO TO 10 C INCOMPATIBLE.RESET ERROR FLAG. IERROR = 1 RETURN 10 NCY = NCA C CLEAR Y TO ZERO. DO 20 I=1,NCA Y(I) = 0.0 20 CONTINUE C IF X IDENTICALLY ZERO SKIP THE PRODUCT. IF (NEX.EQ.0) GO TO 50 C I COUNTS ELEMENTS OF X. DO 40 I=1,NEX C PICK OUT COLUMN INDEX OF I'TH ELEMENT. K = MX(I) T = X(I) C MULTIPLY EACH ELEMENT OF ROW K OF A BY T AND ADD TO Y(J). DO 30 J=1,NCA Y(J) = Y(J) + A(K,J)*T 30 CONTINUE 40 CONTINUE 50 RETURN END SUBROUTINE SPEMSV(A, NRA, NEA, NA, MA, K, X, NRX, NEX, NX, MX, Y, SPE 10 * NRY, NEY, NY, MY, IERROR) C MULTIPLY A SPARSE ELEMENTARY MATRIX (STORED SEQUENTIALLY C IN A,MA,K) ON THE RIGHT BY A SPARSE COLUMN VECTOR C STORED SEQUENTIALLY IN X,MX. C AN ELEMENTARY MATRIX IS ONE IDENTICAL TO A UNIT MATRIX C EXCEPT FOR ONE COLUMN,WHICH IS THE ONLY ONE HAVING C NON-ZERO OFF-DIAGONAL ELEMENTS AND HAVING A NON-UNITY C DIAGONAL ELEMENT. C*** INPUT C A THE NON-ZERO ELEMENTS OF THE NON-UNIT COLUMN OF A, C IN ORDER OF INCREASING ROW NUMBER. C NRA NUMBER OF ROWS IN A. C NEA THE NUMBER OF NON-ZERO ELEMENTS IN A IN THE C NON-UNIT COLUMN C NA THE DIMENSION OF A,MA IN THE CALLING ROUTINE. C MA MA(I) CONTAINS THE ROW-INDEX OF A(I) C K THE COLUMN INDEX OF THE NON-UNIT COLUMN OF A. C X THE NON-ZERO ELEMENTS OF THE COLUMN VECTOR, C IN ORDER OF INCREASING ROW NUMBER. C NRX NUMBER OF ROWS IN X. C NEX THE NUMBER OF NON-ZERO ELEMENTS IN X,MX C NX THE DIMENSION OF X,MX IN THE CALLING ROUTINE. C MX MX(I) CONTAINS THE ROW-INDEX OF X(I). C*** OUTPUT C Y THE NON-ZERO ELEMENTS OF THE PRODUCT VECTOR. C NRY NUMBER OF ROWS IN Y. C NEY THE NUMBER OF NON-ZERO ELEMENTS IN Y,MY. C NY THE DIMENSION OF Y,MY IN THE CALLING ROUTINE. C MY MY(I) CONTAINS THE ROW-INDEX OF Y(I). C IERROR =0 IF NO ERRORS,=1 IF A AND X INCOMPATIBLE. REAL A(NA), X(NX), Y(NY) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(NA),MX(NX),MY(NY) INTEGER MA(NA), MX(NX), MY(NY) C SET ERROR FLAG. IERROR = 0 C TEST FOR COMPATIBILITY. IF (NRA.EQ.NRX) GO TO 10 C INCOMPATIBLE. RESET ERROR FLAG. IERROR = 1 RETURN 10 NRY = NRA C PICK OUT X(K) DO 20 I=1,NEX IF (MX(I)-K) 20, 30, 100 20 CONTINUE GO TO 100 30 R = X(I) C SET COUNTERS FOR A,Y IA = 1 IY = 1 C MAIN LOOP COUNTS ON INDEX IX OF X-ARRAY. DO 80 IX=1,NEX C PICK OUT ROW-INDEX OF CURRENT X-ELEMENT JX = MX(IX) C SET UP TO AVOID ADDING X-ELEMENT IN K TH ROW T = X(IX) IF (JX.EQ.K) T = 0. C TEST FOR END OF A-ARRAY,IF SO USE X-ELEMENTS ONLY. 40 IF (IA.GT.NEA) GO TO 70 C PICK OUT ROW-INDEX OF CURRENT A-ELEMENT. JA = MA(IA) C COMPARE ROW-INDICES OF A- AND X-ELEMENTS. C PATH DEPENDS ON WHICH IS EARLIER. IF (JA-JX) 50, 60, 70 C IF A-ELEMENT IN EARLIER ROW,IGNORE X-ELEMENT. 50 Y(IY) = A(IA)*R MY(IY) = JA C UPDATE COUNTERS FOR A AND Y,USE SAME X-ELEMENT. IA = IA + 1 IY = IY + 1 GO TO 40 C C IF A- AND X-ELEMENT IN SAME ROW,USE BOTH. 60 Y(IY) = T + A(IA)*R MY(IY) = JA C UPDATE COUNTERS FOR A AND Y AND GO TO END OF LOOP C TO UPDATE COUNTER FOR X. IA = IA + 1 IY = IY + 1 GO TO 80 C C IF X-ELEMENT IN EARLIER ROW THAN A-ELEMENT, IGNORE C A-ELEMENT. IF ROW -INDEX IS K,IGNORE EVERYTHING C (X-ELEMENT WAS STORED IN T,EXCEPT FOR ROW K,WHEN T=0.) 70 IF (T.EQ.0.0) GO TO 80 Y(IY) = T MY(IY) = JX IY = IY + 1 C END OF MAIN LOOP 80 CONTINUE C ALL X-ELEMENTS PROCESSED. CHECK IF ANY A-ELEMENTS C REMAIN. IF SO,PROCESS THEM. 90 IF (IA.GT.NEA) GO TO 120 Y(IY) = A(IA)*R MY(IY) = MA(IA) IY = IY + 1 IA = IA + 1 GO TO 90 C CASE WHERE X(K) IS ZERO. A-ELEMENTS ALL IGNORED. 100 IY = 1 DO 110 I=1,NEX J = MX(I) Y(IY) = X(I) MY(IY) = J IY = IY + 1 110 CONTINUE C ALL ELEMENTS PROCESSED. RECORD NUMBER OF ELEMENTS IN Y. 120 NEY = IY - 1 RETURN END SUBROUTINE MSMGUS(A, MA, B, MB, IB, NDC, C, NDMC, MC, IC, IERROR, MSM 10 * X, XB) C MULTIPLY SPARSE MATRICES BY THE METHOD OF GUSTAFSON,ACM T.O.M.S. C VOL 4 (1978) P250. C*** INPUT C A A ONE-DIMENSIONAL ARRAY CONTAINING THE NON-ZERO ELEMENTS C OF THE FIRST MATRIX,ARRANGED ROW-WISE, BUT NOT C NECESSARILY IN ORDER WITHIN ROWS. C MA(1) NUMBER OF ROWS IN A,LATER CALLED P. C MA(2) NUMBER OF COLUMNS IN A,LATER CALLED Q. C MA(3) BLANK. C MA(4) NUMBER OF NON-ZERO ELEMENTS IN A. C MA(5) NUMBER OF NON-ZERO ELEMENTS IN ROW 1 OF A. C ...MA(4+I)... NUMBER OF NON-ZERO ELEMENTS IN ROW I OF A. C ... C ..MA(4+P+1) COLUMN INDEX OF 1ST NON-ZERO ELEMENT OF A. C ... C MA(4+P+I) COLUMN INDEX OF I'TH NON-ZERO ELEMENT OF A. C B AS A,BUT FOR SECOND MATRIX. C MB AS MA,BUT FOR SECOND MATRIX. C NDC DIMENSION OF C (PRODUCT MATRIX) IN CALLING ROUTINE. C NDMC DIMENSION OF MC (CONTROL AND INDEX MATRIX FOR PRODUCT) C IN CALLING ROUTINE. C*** OUTPUT C C A ONE-DIMENSIONAL ARRAY CONTAINING THE NON-ZERO ELEMENTS C OF THE PRODUCT MATRIX,ARRANGED ROW-BY-ROW,BUT NOT C USUALLY IN ORDER WITHIN ROWS. C ORDERING BY INCREASING COLUMN NUMBER MAY BE ATTAINED C BY SUBSEQUENTLY CALLING TRSMGU TWICE,I.E.: C CALL TRSMGU(C,MC,IC,CT,MCT,ICT) C CALL TRSMGU(CT,MCT,IC,C,MC,IC) C WHERE CT,MCT,ICT ARE WORKING STORAGE AREAS OF C SAME DIMENSION AS C,MC,IC RESPECTIVELY. C MC AS MA,BUT FOR PRODUCT MATRIX. C IERROR =1 IF SPACE EXCEEDED IN C,=2 IF A AND B INCOMPATIBLE, C =0 OTHERWISE. C*** WORKING STORAGE PARAMETERS. C IB IB(I) IS ADDRESS IN B OF FIRST ELEMENT OF ROW I OF B. C IB(NUMBER OF ROWS +1)=NUMBER OF ELEMENTS IN B,+1. C IC AS ABOVE,BUT FOR C. C X A ONE-DIMENSIONAL ARRAY OF SIZE GE NUMBER OF COLS OF C, C TO CONTAIN ELEMENTS OF CURRENT ROW OF C, C IN FULL,I.E. NON-SPARSE FORM. C XB AN ARRAY OF SAME SIZE AS X. XB(J)=I IF ELEMENT IN ROW I, C COLUMN J OF C IS NON-ZERO. REAL A(1), B(1), C(1), X(1) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(1),MB(1),IB(1),MC(1),IC(1),XB(1) INTEGER MA(1), MB(1), IB(1), MC(1), IC(1), XB(1) INTEGER P, Q, R, V, VP, PP4, QP4, VPPPP4 C CHECK COMPATIBILITY OF A,B IF (MA(2).EQ.MB(1)) GO TO 10 IERROR = 2 RETURN C UNPACK CONTROL INFORMATION FROM MA,MB C P=NUMBER OF ROWS IN A. C Q=NUMBER OF COLUMNS IN A(SHOULD EQUAL NUMBER OF ROWS IN B). C R=NUMBER OF COLUMNS IN B(AND C). 10 P = MA(1) Q = MA(2) R = MB(2) MC(1) = P MC(2) = R IB(1) = 1 DO 20 I=1,Q IB(I+1) = IB(I) + MB(I+4) 20 CONTINUE IERROR = 0 IP = 1 C INITIALIZE THE NON-ZERO -ELEMENT INDICATOR FOR ROW I OF C. DO 30 V=1,R XB(V) = 0 30 CONTINUE C PROCESS THE ROWS OF A. PP4 = P + 4 QP4 = Q + 4 C INEXT WILL POINT TO START OF NEXT ROW,I.E. ROW I+1 INEXT = 1 DO 80 I=1,P IC(I) = IP C ISTART POINTS TO START OF CURRENT ROW. ISTART = INEXT INEXT = INEXT + MA(I+4) I1 = ISTART I2 = INEXT - 1 IF (I1.GT.I2) GO TO 80 C PROCESS ROW I OF A. DO 60 JP=I1,I2 C J IS COLUMN-INDEX OF CURRENT ELEMENT OF A,I.E. ROW-INDEX OF ROW OF B C TO BE PROCESSED. JPPPP4 = JP + PP4 J = MA(JPPPP4) I3 = IB(J) I4 = IB(J+1) - 1 IF (I3.GT.I4) GO TO 60 C PROCESS ROW OF B. DO 50 KP=I3,I4 C K IS COLUMN INDEX OF CURRENT ELEMENT OF B. KPPQP4 = KP + QP4 K = MB(KPPQP4) C CHECK IF CONTRIBUTION ALREADY EXIXTS TO C(I,K) IF (XB(K).EQ.I) GO TO 40 C SET COLUMN-INDEX AND NON-ZERO INDICATOR FOR NEW ELEMENT OF C. IPPPP4 = IP + PP4 IF (IPPPP4.GT.NDMC) GO TO 100 MC(IPPPP4) = K IP = IP + 1 XB(K) = I X(K) = A(JP)*B(KP) GO TO 50 C ADD NEW CONTRIBUTION TO EXISTING ELEMENT OF C 40 X(K) = X(K) + A(JP)*B(KP) 50 CONTINUE 60 CONTINUE C CHECK FOR OVERFLOW IN C. IF ((IP-1).GT.NDC) GO TO 100 I5 = IC(I) I6 = IP - 1 C EXTRACT NON-ZEROS FROM CURRENT ROW OF C (STORED IN X). DO 70 VP=I5,I6 VPPPP4 = VP + PP4 V = MC(VPPPP4) C(VP) = X(V) 70 CONTINUE 80 CONTINUE C IC(P+1)= NUMBER OF NON-ZEROS IN C,+1. IC(P+1) = IP C EXTRACT CONTROL INFORMATION IN REQUIRED FORM FOR MC. DO 90 I=1,P MC(I+4) = IC(I+1) - IC(I) 90 CONTINUE MC(3) = 0 MC(4) = IP - 1 RETURN 100 IERROR = 1 RETURN END SUBROUTINE TRSMGU(A, MA, IA, AT, MAT, IAT) TRS 10 C TRANSPOSE A SPARSE MATRIX USING A DISTRIBUTION COUNT SORT. C*** INPUT C A A ONE-DIMENSIONAL ARRAY CONTAINING THE NON-ZERO C ELEMENTS OF MATRIX A,ARRANGED ROW-WISE,BUT IN C GENERAL NOT ORDERED WITHIN ROWS. C MA AS IN MSMGUS. C IA IA(I) IS ADDRESS IN A OF FIRST ELEMENT OF ROW I OF A. C IA(NUMBER OF ROWS + 1) IS NUMBER OF ELEMENTS IN A,+1. C*** OUTPUT C AT A ONE-DIMENSIONAL ARRAY CONTAINING NON-ZERO ELEMENTS C OF MATRIX A TRANSPOSED,ARRANGED ROW-WISE,AND IN C ORDER WITHIN ROWS. C MAT CONTAINS CONTROL INFORMATION AND COLUMN INDICES FOR AT, C IN SAME FORMAT AS MA...SEE MSMGUS. C*** WORKING STORAGE. C IAT IAT(I) IS ADDRESS IN AT OF FIRST ELEMENT IN ROW I OF AT. C IAT(S+1) CONTAINS NUMBER OF ELEMENTS IN AT,PLUS 1. REAL A(1), AT(1) C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(1),IA(1),MAT(1),IAT(1) INTEGER MA(1), IA(1), MAT(1), IAT(1) INTEGER R, S, RP4, SP4 C R,S = NUMBER OF ROWS,COLUMNS IN A. R = MA(1) S = MA(2) C C DETERMINE COLUMN COUNTS OF MATRIX A(I.E ROW COUNTS OF AT) C IN ARRAY IAT. C DO 10 I=1,S IAT(I) = 0 10 CONTINUE IAT(S+1) = 0 RP4 = R + 4 SP4 = S + 4 NE = IA(R+1) - 1 DO 20 I=1,NE IPRP4 = I + RP4 K = MA(IPRP4) IAT(K) = IAT(K) + 1 20 CONTINUE C C CALCULATE ROW POINTERS OF AT FROM COLUMN COUNTS OBTAINED ABOVE. C POINTER FOR ROW I STORED IN LOCATION I+1. C ITEMP1 = IAT(1) ITEMP2 = IAT(2) IAT(2) = 1 IF (S.LE.1) GO TO 40 DO 30 I=2,S ITEMP3 = IAT(I+1) IAT(I+1) = IAT(I) + ITEMP1 ITEMP1 = ITEMP2 ITEMP2 = ITEMP3 30 CONTINUE C C CALCULATE COLUMN INDICES(IN ARRAY MAT) AND NUMERICAL VALUES (IN ARRAY C AT) OF MATRIX A-TRANSPOSED USING THE LIST POINTERS IAT(I+1) WHICH C ALWAYS POINT TO THE NEXT ELEMENT TO BE ENTERED IN ROW I OF AT. 40 DO 60 I=1,R J1 = IA(I) J2 = IA(I+1) - 1 IF (J1.GT.J2) GO TO 60 DO 50 JP=J1,J2 JPPRP4 = JP + RP4 J = MA(JPPRP4) JPT = IAT(J+1) JPPSP4 = JPT + SP4 MAT(JPPSP4) = I AT(JPT) = A(JP) IAT(J+1) = JPT + 1 50 CONTINUE 60 CONTINUE C C NOW ALL ROW POINTERS IAT FOR AT HAVE CORRECT VALUES EXCEPT POSITION 1. C FIX IT. IAT(1) = 1 C EXTRACT NEEDED CONTROL INFORMATION IN MAT. DO 70 I=1,S MAT(I+4) = IAT(I+1) - IAT(I) 70 CONTINUE MAT(1) = S MAT(2) = R MAT(3) = 0 MAT(4) = NE RETURN END SUBROUTINE READSM(A, NA, MA, NM, J) REA 10 C READ AND PRINT SPARSE MATRIX A AND ITS ARRAY MA OF C CONTROL DATA AND COLUMN INDICES. C *** INPUT: C J ORDER OF MATRIX IN OVERALL INPUT. C NA DIMENSION OF A IN CALLING ROUTINE. C NM DIMENSION OF MA IN CALLING ROUTINE. C *** OUTPUT: C A ONE-DIMENSIONAL ARRAY CONTAINING NON- C ZERO ELEMENTS OF MATRIX,IN ORDER OF C INCREASING ROW NUMBER AND INCREASING C COLUMN NUMBER WITHIN ROWS. C MA(1) NUMBER OF ROWS IN MATRIX(LATER CALLED NR) C MA(2) NUMBER OF COLS IN MATRIX. C MA(3) NULL C MA(4) NUMBER OF NON-ZERO ELEMENTS IN MATRIX. C MA(5) NUMBER OF NON-ZERO ELEMENTS IN ROW 1. C ... C ... C MA(4+I) NUMBER OF NON-ZERO ELEMENTS IN ROW I. C ... C ... C MA(4+NR+1) COL-INDEX OF FIRST NON-ZERO IN ROW 1. C ... C MA(4+NR+I) COL-INDEX OF A(I). C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(NM) INTEGER MA(NM) REAL A(NA) DATA LP, IR /6,5/ READ (IR,99999) N, (A(I),I=1,N) 99999 FORMAT (I5/(7F10.5)) WRITE (LP,99998) J, (A(I),I=1,N) 99998 FORMAT (45H0NON-ZERO ELEMENTS OF SPARSE MATRIX NUMBERED , I5, * 12H IN INPUT = /(1X, 8F10.5/)) READ (IR,99997) NIM, (MA(I),I=1,NIM) 99997 FORMAT (I5/(14I5)) WRITE (LP,99996) (MA(I),I=1,NIM) 99996 FORMAT (34H0CONTROL AND COLUMN-INDEX ARRAY = /(1X, 16I5/)) RETURN END SUBROUTINE READFM(A, NRA, NCA, MRA, MCA, K) REA 10 C READ AND PRINT FULL MATRIX A. C *** INPUT: C MRA ROW-DIMENSION OF A IN CALLING ROUTINE. C MCA COL-DIMENSION OF A IN CALLING ROUTINE. C K ORDER OF MATRIX IN OVERALL INPUT. C *** OUTPUT: C A A TWO-DIMENSIONAL ARRAY CONTAINING C THE MATRIX IN STANDARD FORM. C NRA NUMBER OF ROWS IN A. C NCA NUMBER OF COLS INA. C DIMENSION A(MRA,MCA) DATA LP, IR /6,5/ READ (IR,99999) NRA, NCA 99999 FORMAT (2I5) WRITE (LP,99998) K 99998 FORMAT (22H0FULL MATRIX NUMBERED , I5, 3H = ) DO 10 I=1,NRA READ (IR,99997) (A(I,J),J=1,NCA) 99997 FORMAT ((7F10.5)) WRITE (LP,99996) (A(I,J),J=1,NCA) 99996 FORMAT ((1X, 8F10.5/)) 10 CONTINUE RETURN END SUBROUTINE PRNTFM(A, NRA, NCA, MRA, MCA) PRN 10 C PRINT A FULL MATRIX IN A. C *** INPUT: C A TWO-DIMENSIONAL ARRAY CONTAINING C THE MATRIX IN STANDARD FORM. C NRA NUMBER OF ROWS IN A. C NCA NUMBER OF COLS IN A. C MRA ROW-DIMENSION OF A IN CALLING ROUTINE. C MCA COL-DIMENSION OF A IN CALLING ROUTINE. C*** OUTPUT: C A ON THE LINE-PRINTER,ROW BY ROW. C DIMENSION A(MRA,MCA) DATA LP /6/ WRITE (LP,99999) 99999 FORMAT (41H0RESULTING FULL MATRIX FOLLOWS,ROW BY ROW) DO 10 I=1,NRA WRITE (LP,99998) (A(I,J),J=1,NCA) 99998 FORMAT ((1X, 8F10.5/)) 10 CONTINUE RETURN END SUBROUTINE READFV(X, NRX, NX, J) REA 10 C READ AND PRINT FULL VECTOR X. C*** INPUT: C J ORDER OF VECTOR IN OVERALL INPUT. C NX DIMENSION OF X IN CALLING ROUTINE. C***OUTPUT: C X A ONE DIMENSIONAL ARRAY CONTAINING ELEMENTS OF THE C FULL VECTOR C NRX NUMBER OF ROWS(OR COLS) IN X. DIMENSION X(NX) DATA LP, IR /6,5/ READ (IR,99999) NRX, (X(I),I=1,NRX) 99999 FORMAT (I5/(7F10.5)) WRITE (LP,99998) J, (X(I),I=1,NRX) 99998 FORMAT (21H0FULL VECTOR NUMBERED, I5, 3H = /(1X, 8F10.5/)) RETURN END SUBROUTINE PRNTFV(X, NRX, NX) PRN 10 C PRINT A FULL VECTOR IN X. C *** INPUT: C X ONE-DIMENSIONAL ARRAY CONTAINING C THE ELEMENTS OF THE FULL VECTOR. C NRX NUMBER OF ROWS(COLS) IN X. C NX DIMENSION OF X IN CALLING ROUTINE. C *** OUTPUT: C ELEMENTS OF X ON LINE-PRINTER. DIMENSION X(NX) DATA LP /6/ WRITE (LP,99998) NRX WRITE (LP,99999) (X(I),I=1,NRX) 99999 FORMAT (37H0ELEMENTS OF RESULTING FULL VECTOR = /(1X, 8F10.5/)) 99998 FORMAT (36H0NUMBER OF ROWS OR COLS IN RESULTING, 13H FULL VECTOR , C * 2H= , I5) RETURN END SUBROUTINE READSV(X, NRX, NEX, NX, MX, J) REA 10 C READ A SPARSE VECTOR INTO X,MX AND PRINT IT. C*** INPUT: C J ORDER OF VECTOR IN OVERALL INPUT. C NX DIMENSION OF X IN CALLING ROUTINE. C*** OUTPUT: C X AN ARRAY CONTAINING THE NON-ZERO C ELEMENTS OF THE SPARSE VECTOR, C IN ORDER OF INCREASING ROW NUMBER. C NRX NUMBER OF ROWS IN X(INCLUDING ZEROS). C NEX NUMBER OF NON-ZERO ELEMENTS IN X. C MX MX(I) CONTAINS ROW INDEX OF X(I). C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MX(NX) INTEGER MX(NX) REAL X(NX) DATA LP, IR /6,5/ READ (IR,99999) NRX, NEX 99999 FORMAT (2I5) IF (NEX.EQ.0) GO TO 10 READ (IR,99995) (X(I),I=1,NEX) WRITE (LP,99998) J, NRX, (X(I),I=1,NEX) 99998 FORMAT (15H0VECTOR NUMBER , I5, 10H CONTAINS , I5, 5H ROWS/ * 31H0NON-ZERO ELEMENTS OF VECTOR = , 5F10.5/(1X, 8F10.5/)) READ (IR,99997) (MX(I),I=1,NEX) 99997 FORMAT (14I5) WRITE (LP,99996) (MX(I),I=1,NEX) 99996 FORMAT (36H0ROW INDICES OF NON-ZERO ELEMENTS = , 10I5/(1X, 16I5/)) RETURN 99995 FORMAT ((7F10.5)) 10 WRITE (LP,99994) J 99994 FORMAT (15H0VECTOR NUMBER , I5, 8H IS NULL) RETURN END SUBROUTINE READEM(X, NRX, NEX, NX, MX, J) REA 10 C READ NON-UNIT COLUMN OF SPARSE ELEMENTARY MATRIX C INTO X,MX AND PRINT IT. C*** INPUT: C J ORDER OF VECTOR IN OVERALL INPUT. C NX DIMENSION OF X IN CALLING ROUTINE. C*** OUTPUT: C X AN ARRAY CONTAINING THE NON-ZERO C ELEMENTS OF THE NON-UNIT COLUMN, C IN ORDER OF INCREASING ROW NUMBER. C NRX NUMBER OF ROWS IN X(INCLUDING ZEROS). C NEX NUMBER OF NON-ZERO ELEMENTS IN X. C MX MX(I) CONTAINS ROW INDEX OF X(I). C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MX(NX) INTEGER MX(NX) REAL X(NX) DATA LP, IR /6,5/ READ (IR,99999) NRX, NEX, (X(I),I=1,NEX) 99999 FORMAT (2I5/(8F10.5)) WRITE (LP,99998) J, NRX, (X(I),I=1,NEX) 99998 FORMAT (26H0ELEMENTARY MATRIX NUMBER , I5, 10H CONTAINS , I5, * 5H ROWS/37H0NON-ZERO ELEMENTS OF NON-UNIT COL = , 5F10.5/(1X, * 8F10.5/)) READ (IR,99997) (MX(I),I=1,NEX) 99997 FORMAT (16I5) WRITE (LP,99996) (MX(I),I=1,NEX) 99996 FORMAT (36H0ROW INDICES OF NON-ZERO ELEMENTS = , 10I5/(1X, 16I5/)) RETURN END SUBROUTINE PRNTSV(X, NRX, NEX, NX, MX) PRN 10 C PRINT A VECTOR WHICH IS RESULT OF SOME SPARSE MATRIX C OPERATION. C *** INPUT: C X ARRAY CONTAINING NON-ZERO ELEMENTS OF C THE VECTOR C NEX NUMBER OF NON-ZERO ELEMENTS. C NX DIMENSION OF X IN THE CALLING ROUTINE C MX MX(I) CONTAINS ROW(COL) INDEX OF X(I). C *** OUTPUT: C X,MX ON THE LINE PRINTER. C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MX(NX) INTEGER MX(NX) REAL X(NX) DATA LP /6/ WRITE (LP,99997) NRX IF (NEX.EQ.0) GO TO 10 WRITE (LP,99999) (X(I),I=1,NEX) 99999 FORMAT (32H0ELEMENTS OF RESULTING VECTOR = , 5F10.5/(1X, 8F10.5/)) WRITE (LP,99998) (MX(I),I=1,NEX) 99998 FORMAT (27H0ELEMENTS OF INDEX ARRAY = , 10I5/(1X, 16I5/)) 99997 FORMAT (46H0NUMBER OF ROWS OR COLS IN RESULTING VECTOR = , I5) RETURN 10 WRITE (LP,99996) 99996 FORMAT (42H0ALL ELEMENTS OF RESULTING VECTOR ARE ZERO) RETURN END SUBROUTINE PRNTSM(A, MA) PRN 10 C PRINT SPARSE MATRIX IN A WITH CONTROL DATA AND COLUMN INDICES IN MA. C THE FOLLOWING MAY BE CHANGED ON IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(1) INTEGER MA(1) REAL A(1) DATA LP /6/ NE = MA(4) WRITE (LP,99999) (A(I),I=1,NE) 99999 FORMAT (20H0RESULTING MATRIX = /(1X, 8F10.5/)) NR = MA(1) N = NR + NE + 4 WRITE (LP,99998) (MA(I),I=1,N) 99998 FORMAT (33H0CONTROL DATA AND COL INDICES = /(1X, 16I5/)) RETURN END SUBROUTINE CNSFMX(A, NA, MA, NM, B, NRB, NCB, MRB, MCB) CNS 10 C CONVERT MATRIX STORED IN SPARSE MODE IN A,MA TO FULL C STORAGE MODE IN B. C *** INPUT: C A,MA AS IN USUAL SPARSE MODE. C MRB,MCB ROW AND COLUMN DIMENSION OF B IN CALLING ROUTINE. C *** OUTPUT C B TWO DIMENSIONAL ARRAY CONTAINING SAME MATRIX C AS IN A BUT STORED IN FULL MODE INCLUDING ZEROS. C NRB NUMBER OF ROWS IN B C NCB NUMBER OF COLUMNS IN B. REAL A(NA), B(MRB,MCB) C FOLLOWING MAY BE ALTERED FOR IBM/370 TYPE MACHINES BY:- C INTEGER*2 MA(NM) INTEGER MA(NM) C EXTRACT CONTROL INFORMATION. NRB = MA(1) NCB = MA(2) C CLEAR WHOLE OF B TO ZERO INITIALLY. DO 20 I=1,NRB DO 10 J=1,NCB B(I,J) = 0. 10 CONTINUE 20 CONTINUE C K2 WILL POINT TO END OF CURRENT ROW I. K2 = 0 C I COUNTS ROWS OF A,B. DO 40 I=1,NRB NIRI = MA(I+4) C IF NO NON-ZEROS IN ROW I SKIP PROCESSING OF THAT ROW. IF (NIRI.EQ.0) GO TO 40 C K1 POINTS TO START OF CURRENT ROW I. K1 = K2 + 1 K2 = K2 + NIRI DO 30 K=K1,K2 NRBK4 = NRB + K + 4 J = MA(NRBK4) B(I,J) = A(K) 30 CONTINUE 40 CONTINUE RETURN END SUBROUTINE MUFUMX(A, NRA, NCA, NRDA, NCDA, B, NRB, NCB, NRDB, MUF 10 * NCDB, C, NRC, NCC, NRDC, NCDC) C MULTIPLY FULL MATRICES A B ,RESULT IN C REAL A(NRDA,NCDA), B(NRDB,NCDB), C(NRDC,NCDC) NRC = NRA NCC = NCB IF (NCA.EQ.NRB) GO TO 10 LP = 6 WRITE (LP,99999) 99999 FORMAT (31H A AND B INCOMPATIBLE IN MUFUMX) RETURN 10 DO 40 I=1,NRA DO 30 J=1,NCB SUM = 0.0 DO 20 K=1,NCA SUM = SUM + A(I,K)*B(K,J) 20 CONTINUE C(I,J) = SUM 30 CONTINUE 40 CONTINUE RETURN END SUBROUTINE CMFUMX(A, NRA, NCA, NRDA, NCDA, B, NRB, NCB, NRDB, CMF 10 * NCDB, IYESNO) C COMPARE TWO FULL MATRICES IN A AND B. IF IDENTICAL,IYESNO=0, C OTHERWISE IYESNO=1. REAL A(NRDA,NCDA), B(NRDB,NCDB) IYESNO = 0 DO 20 I=1,NRA DO 10 J=1,NCA IF (A(I,J).EQ.B(I,J)) GO TO 10 IYESNO = 1 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE ADFUMX(A, NRA, NCA, NRDA, NCDA, B, NRB, NCB, NRDB, ADF 10 * NCDB, C, NRC, NCC, NRDC, NCDC) C ADD TWO FULL MATRICES AN A,B, RESULT IN C. REAL A(NRDA,NCDA), B(NRDB,NCDB), C(NRDC,NCDC) DO 20 I=1,NRA DO 10 J=1,NCA C(I,J) = A(I,J) + B(I,J) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE MUFMXV(A, NRA, NCA, NRDA, NCDA, X, NRX, NX, Y, NRY, NY)MUF 10 C MULTIPLY FULL MATRIX IN A BY FULL VECTOR IN X. C RESULT IN Y. REAL A(NRDA,NCDA), X(NX), Y(NY) NRY = NRA DO 20 I=1,NRA S = 0.0 DO 10 K=1,NCA S = S + A(I,K)*X(K) 10 CONTINUE Y(I) = S 20 CONTINUE RETURN END SUBROUTINE CMFUVC(X, NRX, NX, Y, NRY, NY, IYESNO) CMF 10 C COMPARE TWO FULL VECTORS IN X AND Y. IYESNO=1 IF THEY ARE DIFFERENT. REAL X(NX), Y(NY) IYESNO = 0 DO 10 I=1,NRX IF (X(I).EQ.Y(I)) GO TO 10 IYESNO = 1 10 CONTINUE RETURN END SUBROUTINE MUFVMX(X, NCX, NX, A, NRA, NCA, NRDA, NCDA, Y, NCY, NY)MUF 10 C MULTIPLY FULL VECTOR IN X BY FULL MATRIX IN A. RESULT IN Y. REAL X(NX), A(NRDA,NCDA), Y(NY) NCY = NCA DO 20 J=1,NCA S = 0.0 DO 10 K=1,NRA S = S + X(K)*A(K,J) 10 CONTINUE Y(J) = S 20 CONTINUE RETURN END SUBROUTINE CNSFVC(X, NRX, NEX, NX, MX, Y, NRY, NY) CNS 10 C CONVERT SPARSE VECTOR IN X TO FULL VECTOR IN Y. REAL X(NX), Y(NY) C FOLLOWING MAY BE ALTERED FOR IBM /370 TYPE MACHINES BY:- C INTEGER*2 MX(NX) INTEGER MX(NX) NRY = NRX DO 10 I=1,NRX Y(I) = 0 10 CONTINUE IF (NEX.EQ.0) GO TO 30 DO 20 I=1,NEX J = MX(I) Y(J) = X(I) 20 CONTINUE 30 RETURN END SUBROUTINE CNEFMX(A, NRA, NEA, NDA, MA, K, B, NRB, NCB, NRDB, CNE 10 * NCDB) C CONVERT ELEMENTARY MATRIX IN A (NON-ZERO COL K) TO FULL MATRIX IN B. REAL A(NDA), B(NRDB,NCDB) C THE FOLLOWING MAY BE ALTERED ON IBM/370 TYPE MACHINES TO:- C INTEGER*2 MA(NDA) INTEGER MA(NDA) NRB = NRA NCB = NRA DO 20 I=1,NRA DO 10 J=1,NRA B(I,J) = 0.0 IF (I.EQ.J .AND. I.NE.K) B(I,J) = 1.0 10 CONTINUE 20 CONTINUE DO 30 L=1,NEA I = MA(L) B(I,K) = A(L) 30 CONTINUE RETURN END 5 00000100 1. 2. 1. 1. 1. 00000200 13 00000300 4 4 0 5 2 1 2 0 1 4 2 1 3 00000400 4 4 00000500 1. 2. 3. 4. 00000600 4. 3. 2. 1. 00000700 2. 3. 4. 5. 00000800 5. 4. 3. 2. 00000900 5 00001000 1. 2. 3. 4. 5. 00001100 12 00001200 3 4 0 5 2 1 2 1 4 2 1 4 00001300 4 3 00001400 1. 2. 3. 00001500 4. 5. 6. 00001600 6. 5. 4. 00001700 3. 2. 0. 00001800 5 00001900 1. 2. 3. 4. 5. 00002000 12 00002100 3 3 0 5 2 1 2 1 3 2 1 3 00002200 4 3 00002300 1. 2. 3. 00002400 4. 5. 6. 00002500 6. 5. 4. 00002600 3. 2. 0. 00002700 4 4 00002800 1. 2. 3. 4. 00002900 4. 3. 2. 1. 00003000 2. 3. 4. 5. 00003100 5. 4. 3. 2. 00003200 5 00003300 1. 2. 1. 1. 1. 00003400 13 00003500 4 4 0 5 2 1 2 0 1 4 2 1 3 00003600 4 3 00003700 1. 2. 3. 00003800 4. 5. 6. 00003900 6. 5. 4. 00004000 3. 2. 0. 00004100 5 00004200 1. 2. 3. 4. 5. 00004300 12 00004400 3 4 0 5 2 1 2 1 4 2 1 4 00004500 3 3 00004600 1. 2. 3. 00004700 3. 2. 1. 00004800 1. 2. 3. 00004900 3 00005000 1. 2. 3. 00005100 10 00005200 3 4 0 3 2 0 1 1 3 4 00005300 3 3 00005400 1. 2. 3. 00005500 3. 2. 1. 00005600 1. 2. 3. 00005700 5 00005800 1. 2. 1. 1. 1. 00005900 13 00006000 4 4 0 5 2 1 2 0 1 4 2 1 3 00006100 3 00006200 1. 2. 3. 00006300 10 00006400 3 4 0 3 2 1 0 1 3 2 00006500 3 4 00006600 1. 2. 3. 4. 00006700 5. 6. 7. 8. 00006800 9. 8. 7. 6. 00006900 3 00007000 1. 2. 3. 00007100 10 00007200 3 4 0 3 2 1 0 1 3 2 00007300 4 4 00007400 1. 2. 3. 4. 00007500 5. 6. 7. 8. 00007600 9. 8. 7. 6. 00007700 9. 8. 7. 6. 00007800 3 00007900 1. 2. 3. 00008000 10 00008100 3 4 0 3 2 1 0 1 3 3 00008200 4 00008300 1. 2. 3. 4. 00008400 3 00008500 1. 2. 3. 00008600 10 00008700 3 4 0 3 2 1 0 1 3 3 00008800 3 00008900 1. 2. 3. 00009000 3 00009100 1. 2. 3. 00009200 3 00009300 1. 2. 3. 00009400 10 00009500 3 4 0 3 2 1 0 1 3 3 00009600 4 00009700 1. 2. 3. 4. 00009800 3 00009900 1. 2. 3. 00010000 10 00010100 3 4 0 3 2 1 0 1 3 3 00010200 2 3 00010300 1. 2. 3. 00010400 4. 5. 6. 00010500 3 2 00010600 1. 2. 00010700 1 3 00010800 2 3 00010900 1. 2. 3. 00011000 4. 5. 6. 00011100 2 2 00011200 1. 2. 00011300 1 2 00011400 2 3 00011500 1. 2. 3. 00011600 4. 5. 6. 00011700 3 0 00011800 3 2 00011900 1. 2. 00012000 1 3 00012100 3 2 00012200 1. 4. 00012300 2. 5. 00012400 3. 6. 00012500 3 2 00012600 1. 2. 00012700 1 3 00012800 2 2 00012900 1. 4. 00013000 2. 5. 00013100 3 0 00013200 3 2 00013300 1. 4. 00013400 2. 5. 00013500 3. 6. 00013600 5 2 00013700 2. 3. 00013800 1 3 00013900 3 00014000 5 3 00014100 5. 2. 4. 00014200 2 3 5 00014300 5 2 00014400 2. 3. 00014500 1 3 00014600 3 00014700 5 3 00014800 5. 2. 4. 00014900 2 4 5 00015000 5 2 00015100 2. 3. 00015200 1 4 00015300 3 00015400 5 2 00015500 5. 2. 00015600 2 3 00015700 5 2 00015800 2. 3. 00015900 1 4 00016000 3 00016100 4 2 00016200 5. 2. 00016300 2 3 00016400 4 00016500 1. 2. 3. 4. 00016600 11 00016700 3 4 0 4 2 2 0 1 3 2 4 00016800 2 00016900 2. 1. 00017000 10 00017100 4 2 0 2 0 1 1 0 2 1 00017200 4 00017300 1. 2. 3. 4. 00017400 11 00017500 3 4 0 4 2 2 0 1 3 2 4 00017600 6 00017700 1. 2. 3. 4. 5. 6. 00017800 14 00017900 4 4 0 6 1 2 3 0 3 2 4 1 3 4 00018000 4 00018100 1. 2. 3. 4. 00018200 11 00018300 3 4 0 4 2 2 0 1 3 2 4 00018400 1 00018500 1. 00018600 9 00018700 4 1 0 1 0 0 1 0 1 00018800 6 00018900 1. 2. 3. 4. 5. 6. 00019000 16 00019100 6 1 0 6 1 1 1 1 1 1 1 1 1 1 00019200 1 1 00019300 6 00019400 1. 2. 3. 4. 5. 6. 00019500 11 00019600 1 6 0 6 6 1 2 3 4 5 6 00019700 4 00019800 1.0 2.0 3.0 4.0 00019900 11 00020000 3 4 0 4 2 2 0 1 3 2 4 00020100 2 00020200 2.0 1.0 00020300 10 00020400 3 2 0 2 0 1 1 0 2 1 00020500