Trailing-Edge
-
PDP-10 Archives
-
decuslib10-01
-
43,50127/matrix.f4
There is 1 other file named matrix.f4 in the archive. Click here to see a list.
00010 C EXECUTIVE ROUTINE FOR MATRIX INTERPRETIVE PROGRAM
00020 INTEGER FREZ,FRET
00030 INTEGER OPLIST(60),WORD(15),ERROR,WORD1,BLANK,SAVNM(2)
00040 COMMON /FF/ NFMX,ITAPE,NWMX
00050 COMMON /XX/ ERROR
00060 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00070 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1),NME(4),NSIZE
00080 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00090 INTEGER CHAR(0/9)
00100 DIMENSION FNUMB(10)
00110 C * * * * * * * * * * * * * * LIST OF OPERATIONS * * *
00120 DATA (OPLIST(I),I=1,5)/'ZERO','REMAR','PRINT','MAP','ENVEL'/
00130 DATA (OPLIST(I),I=6,10)/'LOAD','DUPL','DELET','SCALE','MSCAL'/
00140 DATA OPLIST(11) /'SQREL'/
00150 DATA (OPLIST(I),I=12,16)/'INVEL','LOG','ADD','SUB','TRANS'/
00160 DATA (OPLIST(I),I=17,20)/'MULT','TRMPY','SYMCK','SYMNV'/
00170 DATA (OPLIST(I),I=21,24)/'INVRT','STODG','RMVDG','IDADD'/
00180 DATA (OPLIST(I),I=25,28)/'DGADD','DGSUB','DGPRE','DGPST'/
00190 DATA (OPLIST(I),I=29,32)/'DGMPY','FORMK','FORMD','MERGE'/
00200 DATA (OPLIST(I),I=33,36)/'STOSM','ADDSM','RMVSM','DELRC'/
00210 DATA (OPLIST(I),I=37,40)/'SELCT','CHOL1','CHOL2','CHOL3'/
00220 DATA (OPLIST(I),I=41,44)/'EIGEN','FUNGN','RESPN','LIST'/
00230 DATA (OPLIST(I),I=45,48)/'EDIT','GET','SAVE','FEM'/
00240 DATA OPLIST(49) /'RANK'/
00250 DATA (CHAR(I),I=0,9) /'0','1','2','3','4','5','6','7',
00260 1 '8','9'/
00270 DATA BLANK/' '/
00280 EXTERNAL SUBERR
00290 C * * * * * * * * * * * * * * INITIALIZE BY ZEROING ARRAY * * *
00300 C SET PSIZE TO 44 AND MAXIMUM IN MLIST TO 30K FOR MATRXB
00310 CALL PSIZE(23)
00320 C CALL PSIZE(44)
00330 CALL EXSIZE(SUBERR)
00340 NNUM=10
00350 NSIZE=0
00360 NUM=0
00370 DO 2 I=1,4
00380 IADDR=FREZ(10)
00390 2 NME(I)=IADDR
00400 INC1=NME(1)-LOC(NAME)
00410 INC2=NME(2)-LOC(NCOL)
00420 INC3=NME(3)-LOC(NROW)
00430 INC4=NME(4)-LOC(NN)
00440 NOPER=49
00450 TYPE 4
00460 4 FORMAT(' MATRIX REV. 7 JULY 1969'///)
00470 5 NF=0
00480 NW=1
00490 ITAPE=0
00500 NFMX=10
00510 NWMX=15
00520 ERROR=0
00530 DO 50 I=1,5
00540 WORD(I)=BLANK
00550 FNUMB(I)=0
00560 50 CONTINUE
00570 DO 55 I=6,10
00580 FNUMB(I)=0
00590 55 CONTINUE
00600 TYPE 6
00610 6 FORMAT(' *'$)
00620 CALL FREAD(NW,WORD,NF,FNUMB)
00630 WORD1=WORD(1)
00640 I=1
00650 8 IF(NOPER-I)9,10,10
00660 9 TYPE 9001
00670 9001 FORMAT(' ILLEGAL COMMAND'/)
00680 GO TO 5
00690 10 IF(OPLIST(I).EQ.WORD1) GO TO 15
00700 I=I+1
00710 GO TO 8
00720 15 IDITO=I
00730 GO TO (100,5,120,140,160,180,200,220,240,260,280,300,
00740 1 320,340,350,380,400,420,440,460,480,500,520,540,560,
00750 2 580,600,620,640,660,680,700,720,740,760,780,800,
00760 3 820,840,860,880,1000,1020,2000,190,2020,2040,2060,2080), IDITO
00770 C * * * * * * ZEROS A MATRIX * * * * * *
00780 100 NAMA=WORD(2)
00790 IF ((NW.NE.2).OR.(NF.NE.2)) GO TO 995
00800 NR=FNUMB(1)
00810 NC=FNUMB(2)
00820 MPYTAG=0
00830 CALL MLIST(NAMA)
00840 IF(ERROR.GT.0) GO TO 990
00850 DO 110 I=1,NR*NC
00860 A(I+INCA)=0.0
00870 110 CONTINUE
00880 GO TO 5
00890 C * * * * * * * * * * * * * * MATRIX PRINT - PRINT* * *
00900 120 NAMA=WORD(2)
00910 CALL MFIND(NAMA)
00920 IF(ERROR.GT.0) GO TO 990
00930 CALL PRINT(NAMA)
00940 GO TO 5
00950 C * * * * * * * * * * * * * * MATRIX MAP - MAP* * *
00960 140 NAMA=WORD(2)
00970 SCAL1=FNUMB(1)
00980 CALL MFIND(NAMA)
00990 IF(ERROR.GT.0) GO TO 990
01000 CALL MTXMAP(SCAL1)
01010 GO TO 5
01020 C * * * * * * * * * * * * * * ENVELOPE VALUES OF MATRIC - ENVEL* * *
01030 160 SCAL1=FNUMB(2)
01040 IF(NF.EQ.1) TYPE 1601
01050 1601 FORMAT(' INTERVAL SCALAR OMITTED'/)
01060 IF ((NW.LT.2).OR.(NF.LT.2)) GO TO 995
01070 NAMA=WORD(2)
01080 KCOL=FNUMB(1)
01090 CALL MFIND(NAMA)
01100 IF(ERROR.GT.0) GO TO 990
01110 NA=INCA
01120 IF(KCOL) 998,163,164
01130 163 NRA=NR
01140 NCA=NC
01150 GO TO 165
01160 164 NCA=NR
01170 NRA=NC
01180 165 NAMB=WORD(3)
01190 IF(NAMB.EQ.BLANK) GO TO 175
01200 NR=1
01210 NC=NRA
01220 CALL MLIST(NAMB)
01230 IF(ERROR.GT.0) GO TO 990
01240 NBB=1
01250 GO TO 176
01260 175 NBB=0
01270 176 NR=NRA
01280 NC=NCA
01290 CALL ENVEL(NA,NBB,SCAL1,KCOL)
01300 GO TO 5
01310 C * * * * * * * * * * * * * * LOAD MATRIX - LOAD * * *
01320 180 IF (NW.LT.2.OR.NF.NE.2) GO TO 995
01330 IF (NW.EQ.3) GO TO 184
01340 NAMA=WORD(2)
01350 GO TO 185
01360 184 IF (WORD(2).NE.'TAPE') GO TO 9
01370 NAMA=WORD(3)
01380 ITAPE=1
01390 185 NR=FNUMB(1)
01400 NC=FNUMB(2)
01410 CALL MLIST(NAMA)
01420 IF(ERROR.GT.0) GO TO 990
01430 CALL LOAD
01440 GO TO 5
01450 C * * * * * * * EDIT AN ELEMENT OF A MATRIX * * * * * * * *
01460 190 IF (NW.NE.2.OR.NF.NE.1) GO TO 995
01470 NAMA=WORD(2)
01480 CALL MFIND(NAMA)
01490 IF (ERROR.GT.0) GO TO 990
01500 NTIM=FNUMB(1)
01510 IF (NTIM.LT.1) GO TO 996
01520 TYPE 1900
01530 1900 FORMAT(5X,'ROW-COL COORDINATES - VALUE'/)
01540 DO 195 I=1,NTIM
01550 NW=-3
01560 NF=3
01570 NWMX=1
01580 NFMX=10
01590 CALL FREAD(NW,WORD,NF,FNUMB)
01600 NREL=FNUMB(1)
01610 NCEL=FNUMB(2)
01620 IF (NREL.GT.NR.OR.NCEL.GT.NC) GO TO 901
01630 NA=(NREL-1)*NC+NCEL+INCA
01640 A(NA)=FNUMB(3)
01650 195 CONTINUE
01660 GO TO 5
01670 C * * * * * * * * * * * * * * DUPLICATE MATRIX - DUPL * * *
01680 200 NAMA=WORD(2)
01690 IF (NW.NE.3) GO TO 995
01700 NAMB=WORD(3)
01710 IGO=1
01720 CALL MFIND(NAMA)
01730 IF(ERROR.GT.0) GO TO 990
01740 NB=INCA
01750 CALL MLIST(NAMB)
01760 IF(ERROR.GT.0) GO TO 990
01770 NS=NR*NC
01780 DO 210 I=1,NS
01790 A(I+INCA)=A(I+NB)
01800 210 CONTINUE
01810 GO TO 5
01820 C * * * * * * * * * * * * * * DELETE MATRIX - DELETE * * *
01830 220 NTIME=NW-3+1
01840 DO 225 I=1,NTIME
01850 NAMA=WORD(2+I)
01860 CALL DELETE(NAMA)
01870 IF(ERROR.GT.0) GO TO 990
01880 225 CONTINUE
01890 GO TO 5
01900 C * * * * * * * * * * * * * * SCALAR MULTIPLICATION - SCALE * * *
01910 240 NAMA=WORD(2)
01920 IF ((NW.NE.2).OR.(NF.NE.1)) GO TO 995
01930 SCAL1=FNUMB(1)
01940 241 CALL MFIND(NAMA)
01950 IF(ERROR.GT.0) GO TO 990
01960 DO 245 I=1,NR*NC
01970 A(I+INCA)=SCAL1*A(I+INCA)
01980 245 CONTINUE
01990 GO TO 5
02000 C * * * * * * * * * * * * * * MULTIPLY A*B(NR,NC) - MSCALE * * *
02010 260 IFLAG=1
02020 IF ((NW.NE.3).OR.(NF.NE.2)) GO TO 995
02030 NAMA=WORD(2)
02040 NAMB=WORD(3)
02050 NR=FNUMB(1)
02060 NC=FNUMB(2)
02070 IF(NR*NC) 909,272,273
02080 272 NR=1
02090 NC=1
02100 273 NRB=NR
02110 NCB=NC
02120 CALL MFIND(NAMB)
02130 IF(ERROR.GT.0) GO TO 990
02140 IF((NRB.GT.NR).OR.(NCB.GT.NC)) GO TO 909
02150 IF(NR*NC.EQ.0) GO TO 909
02160 SCAL1=A(INCA+(NRB-1)*NC+NCB)
02170 GO TO 241
02180 C * * * * * * * * * * * * * * SQUARE ROOT OF EACH ELEMENT - SQREL * * *
02190 280 IFLAG=1
02200 281 NAMA=WORD(2)
02210 CALL MFIND(NAMA)
02220 IF(ERROR.GT.0) GO TO 990
02230 CALL ELEMNT(IFLAG)
02240 IF (ERROR.NE.0) GO TO 990
02250 GO TO 5
02260 C * * * * * * * * * * * * * * INVERSION OF EACH ELEMENT - INVEL * * *
02270 300 IFLAG=2
02280 GO TO 281
02290 C * * * * * * * * * * * * * * LOG OF MATRIX - LOG * * *
02300 320 IFLAG=3
02310 GO TO 281
02320 C * * * * * * * * * * * * * * ADD OR SUBTRACT MATRICES - ADD , SUB * * *
02330 340 TAG=1.0
02340 GO TO 351
02350 350 TAG=-1.
02360 351 NAMA=WORD(2)
02370 NAMB=WORD(3)
02380 IF (NW.NE.3) GO TO 995
02390 CALL MFIND(NAMA)
02400 IF(ERROR.GT.0) GO TO 990
02410 NA=INCA
02420 NRA=NR
02430 NCA=NC
02440 CALL MFIND(NAMB)
02450 IF(ERROR.GT.0) GO TO 990
02460 IF ((NR.NE.NRA).OR.(NC.NE.NCA)) GO TO 901
02470 NS=NR*NC
02480 DO 360 I=1,NS
02490 A(NA+I)=A(NA+I)+TAG*A(INCA+I)
02500 360 CONTINUE
02510 GO TO 5
02520 C * * * * * * * * * * * * * * MATRIX TRANSPOSE - TRANS * * *
02530 380 NAMA=WORD(2)
02540 IF (NW.NE.3) GO TO 995
02550 NAMB=WORD(3)
02560 CALL MFIND(NAMA)
02570 IF(ERROR.GT.0) GO TO 990
02580 NRA=NR
02590 NR=NC
02600 NB=INCA
02610 NC=NRA
02620 CALL MLIST(NAMB)
02630 IF(ERROR.GT.0) GO TO 990
02640 CALL TRANS(NB)
02650 GO TO 5
02660 C * * * * * * * * * * * * * * MATRIX RIGHT MULTIPLY - MULT * * *
02670 400 MPYTAG=FNUMB(1)
02680 ITRMPY=0
02690 410 NAMA=WORD(2)
02700 CALL MFIND(NAMA)
02710 IF(ERROR.GT.0) GO TO 990
02720 NA=INCA
02730 NRA=NR
02740 NCA=NC
02750 IF (NW.NE.4) GO TO 995
02760 NAMB=WORD(3)
02770 NAMC=WORD(4)
02780 CALL MFIND(NAMB)
02790 IF(ERROR.GT.0) GO TO 990
02800 NRB=NR
02810 NCB=NC
02820 NB=INCA
02830 NR=NRA
02840 IF(ITRMPY.EQ.0) GO TO 411
02850 NR=NCA
02860 NCA=NRA
02870 411 IF (NCA.NE.NRB) GO TO 901
02880 413 CALL MLIST(NAMC)
02890 IF(ERROR.GT.0) GO TO 990
02900 CALL MULT(NA,NB,NCA,NCB,ITRMPY)
02910 MPYTAG=0
02920 ITRMPY=0
02930 GO TO 5
02940 C * * * * * * * * * * PREMULTIPLY BY TRANSPOSE (AND ADD) -TRMPY
02950 420 MPYTAG=FNUMB(1)
02960 ITRMPY=1
02970 GO TO 410
02980 C * * * * * * * * * * * * * * *SYMMETRY CHECK * * * *
02990 440 L=1
03000 441 NAMA=WORD(2)
03010 SCAL1=FNUMB(1)
03020 442 CALL MFIND(NAMA)
03030 IF(ERROR.GT.0) GO TO 990
03040 IF (NR.NE.NC) GO TO 945
03050 IF(SCAL1) 443,443,445
03060 443 SCAL1=ABS(A(INCA+1)*1.0E-06)
03070 445 IFLAG=0
03080 CALL SYMCHK(SCAL1,IFLAG)
03090 IF(IFLAG) 908,455,450
03100 450 IF(L-1) 908,5,901
03110 455 GO TO (5,465,822,885), L
03120 C * * * * * * * * * * * *INVERSION OF SYMETRIC MATRIX * *
03130 460 L=2
03140 GO TO 441
03150 465 NA=INCA
03160 NRA=NR
03170 467 NAMB=WORD(3)
03180 IF(NAMB.EQ.BLANK) GO TO 470
03190 CALL MFIND(NAMB)
03200 IF(ERROR.GT.0) GO TO 990
03210 IF(NR.EQ.NRA) GO TO 470
03220 ERROR=1
03230 GO TO 990
03240 470 NC=0
03250 471 CALL SYMINV(NA,SCAL1)
03260 IF(ERROR.GT.0) GO TO 990
03270 GO TO 5
03280 C * * * * * * * * * * * * * *MATRIX INVERSION * * *
03290 480 NAMA=WORD(2)
03300 SCAL1=FNUMB(1)
03310 482 CALL MFIND(NAMA)
03320 IF(ERROR.GT.0) GO TO 990
03330 NA=INCA
03340 NRA=NR
03350 IF (NR.NE.NC) GO TO 945
03360 NAMB=WORD(3)
03370 IF(NAMB.EQ.BLANK) GO TO 490
03380 CALL MFIND(NAMB)
03390 IF(ERROR.GT.0) GO TO 990
03400 IF(NR.EQ.NRA) GO TO 492
03410 ERROR=1
03420 GO TO 990
03430 490 NC=0
03440 492 CALL INVERT(NA,SCAL1)
03450 IF(ERROR.GT.0) GO TO 990
03460 GO TO 5
03470 C * * * * * * * * * * * * * STORE ROW (B) ON DIAGONAL OF (A) -STODG
03480 500 NAMA=WORD(2)
03490 NAMB=WORD(3)
03500 IF (NW.NE.3) GO TO 995
03510 CALL MFIND(NAMA)
03520 IF(ERROR.GT.0) GO TO 990
03530 K=INCA+1
03540 IF (NR.NE.NC) GO TO 945
03550 NCA=NC
03560 CALL MFIND(NAMB)
03570 IF(ERROR.GT.0) GO TO 990
03580 IF (NR.NE.1) GO TO 901
03590 IF (NC.NE.NCA) GO TO 901
03600 DO 510 I=1,NC
03610 A(K)=A(I+INCA)
03620 K=K+NCA+1
03630 510 CONTINUE
03640 GO TO 5
03650 C * * * * * * *REMOVE ROW(B) FROM DIAGONAL OF (A) - RMVDG
03660 520 NAMA=WORD(2)
03670 NAMB=WORD(3)
03680 IF (NW.NE.3) GO TO 995
03690 CALL MFIND(NAMA)
03700 IF(ERROR.GT.0) GO TO 990
03710 IF (NR.NE.NC) GO TO 945
03720 NR=1
03730 CALL MLIST(NAMB)
03740 IF(ERROR.GT.0) GO TO 990
03750 NB=INCA
03760 CALL MFIND(NAMA)
03770 K=INCA+1
03780 DO 525 I=1,NC
03790 A(I+NB)=A(K)
03800 K=K+NC+1
03810 525 CONTINUE
03820 GO TO 5
03830 C * * * * * *ADDITION OF SCALAR TIMES IDENTITY MATRIX - IDADD
03840 540 NAMA=WORD(2)
03850 IF ((NW.NE.2).OR.(NF.NE.1)) GO TO 995
03860 SCAL1=FNUMB(1)
03870 CALL MFIND(NAMA)
03880 IF(ERROR.GT.0) GO TO 990
03890 IF (NR.NE.NC) GO TO 945
03900 K=1+INCA
03910 DO 545 I=1,NR
03920 A(K)=A(K)+SCAL1
03930 K=K+NC+1
03940 545 CONTINUE
03950 GO TO 5
03960 C * * * * ADDITION OR SUBTRACTION OF DIAGONAL MATRIX - DGADD,DGSUB
03970 560 TAG=1.0
03980 561 NAMA=WORD(2)
03990 NAMB=WORD(3)
04000 IF (NW.NE.3) GO TO 995
04010 CALL MFIND(NAMA)
04020 IF(ERROR.GT.0) GO TO 990
04030 NA=INCA
04040 IF (NR.NE.NC) GO TO 945
04050 NCA=NC
04060 CALL MFIND(NAMB)
04070 IF(ERROR.GT.0) GO TO 990
04080 IF (NR.NE.1) GO TO 901
04090 IF (NC.NE.NCA) GO TO 901
04100 K=1+NA
04110 DO 565 I=1,NC
04120 A(K)=A(K)+TAG*A(I+INCA)
04130 K=K+NC+1
04140 565 CONTINUE
04150 GO TO 5
04160 580 TAG=-1.0
04170 GO TO 561
04180 C * * * *PREMULTIPLICATION BY DIAGONAL MATRIX - DGPRE * * *
04190 600 IFLAG=1
04200 601 NAMA=WORD(2)
04210 NAMB=WORD(3)
04220 IF (NW.NE.3) GO TO 995
04230 CALL MFIND(NAMA)
04240 IF(ERROR.GT.0) GO TO 990
04250 NA=INCA
04260 NRA=NR
04270 NCA=NC
04280 CALL MFIND(NAMB)
04290 IF(ERROR.GT.0) GO TO 990
04300 IF (NR.NE.1) GO TO 901
04310 GO TO (610,611,612),IFLAG
04320 610 IF(NC-NRA) 901,615,901
04330 611 IF(NC-NCA) 901,615,901
04340 612 IF(NRA-1) 901,611,901
04350 615 CALL DGMPY(NA,NRA,NCA,IFLAG)
04360 GO TO 5
04370 C * * * * * * * * *POSTMULTIPLICATION BY DIAGONAL MATRIX DGPST
04380 620 IFLAG=2
04390 GO TO 601
04400 C * * * * * * * * * *MULTIPLY TWO DIAGONAL MATRICES * * * *
04410 640 IFLAG=3
04420 GO TO 601
04430 C * * * *FORM TWO DIMENSIONAL BEAM STIFFNESS MATRIX(2X2,3X3,4X4)
04440 660 NAMA=WORD(2)
04450 IF ((NW.NE.2).OR.(NF.LT.1)) GO TO 995
04460 NR=FNUMB(1)
04470 NC=FNUMB(2)
04480 IF((NR.LT.0).OR.(NR.GT.4)) GO TO 903
04490 IF(NR-1) 671,671,672
04500 671 NR=2
04510 672 NTYPE=NR
04520 IF(NC) 908,673,674
04530 673 NC=1
04540 674 N=NC
04550 NR=N*NTYPE
04560 NC=NR
04570 CALL MLIST(NAMA)
04580 IF(ERROR.GT.0) GO TO 990
04590 CALL FORMK(NTYPE,N)
04600 IF(ERROR) 990,5,990
04610 C * * * *FORM BEAM MATRIX FOR DIRECT STIFFNESS MERGING(4X4,6X6)
04620 680 NAMA=WORD(2)
04630 IF ((NW.NE.2).OR.(NF.NE.2)) GO TO 995
04640 NTYPE=FNUMB(1)
04650 N=FNUMB(2)
04660 IF(.NOT.(NTYPE.EQ.4.OR.NTYPE.EQ.6.AND.N.GE.0)) GO TO 995
04670 IF(N.LE.0) N=1
04680 NR=NTYPE
04690 NC=NR
04700 DO 685 I=1,5
04710 NPLAC=I
04720 CALL GET(NAMA,I,IC)
04730 IF(IC.EQ.BLANK) GO TO 686
04740 685 CONTINUE
04750 686 NDIG=1
04760 IF(N.GE.10) NDIG=2
04770 IF((NPLAC+NDIG).GT.5) NPLAC=6-NDIG
04780 J=0
04790 NEND=N
04800 NTIME=N/10
04810 IF(NTIME.EQ.0) GO TO 688
04820 NTIME=NTIME+1
04830 NPLAC=NPLAC+1
04840 687 J=J+1
04850 CALL PUT(NAMA,NPLAC-1,CHAR(J-1))
04860 NEND=N-(J-1)*10
04870 IF(NEND.GT.9) NEND=9
04880 IF(J.EQ.1) GO TO 688
04890 CALL PUT(NAMA,NPLAC,CHAR(0))
04900 CALL MLIST(NAMA)
04910 IF(ERROR.GT.0) GO TO 990
04920 IF(NEND.EQ.0) GO TO 692
04930 688 DO 690 I=1,NEND
04940 ICHAR=CHAR(I)
04950 CALL PUT(NAMA,NPLAC,ICHAR)
04960 CALL MLIST(NAMA)
04970 IF(ERROR.GT.0) GO TO 990
04980 690 CONTINUE
04990 IF(NTIME.EQ.0) GO TO 692
05000 IF(J-NTIME) 687,692,692
05010 692 CALL FORMKD(NTYPE,N)
05020 IF(ERROR)990,5,990
05030 C * * * * * * * * * * * * * * * * *MATRIX MERGE * * *
05040 700 IF(.NOT.(NW.EQ.3.AND.NF.EQ.1)) GO TO 995
05050 K=FNUMB(1)
05060 NAMA=WORD(2)
05070 NAMB=WORD(3)
05080 IF(K.EQ.0) GO TO 903
05090 CALL MFIND(NAMA)
05100 IF(ERROR.GT.0) GO TO 990
05110 NA=INCA
05120 NRA=NR
05130 NCA=NC
05140 CALL MFIND(NAMB)
05150 IF(ERROR.GT.0) GO TO 990
05160 IF((NR.GT.100).OR.(NC.GT.100)) GO TO 901
05170 CALL MERGE(NA,NRA,NCA,K)
05180 IF(ERROR) 990,5,990
05190 C * * * * * * * * * * * * * STORE SUBMATRIX * * * * * *
05200 720 NTAG=1
05210 721 IF ((NW.NE.3).OR.(NF.NE.2)) GO TO 995
05220 J=FNUMB(1)
05230 K=FNUMB(2)
05240 NAMB=WORD(3)
05250 CALL MFIND(NAMB)
05260 IF(ERROR.GT.0) GO TO 990
05270 GO TO 770
05280 C * * * * * * * * * * * * * * * * ADD SUBMATRIX * * * *
05290 740 NTAG=-1
05300 GO TO 721
05310 C * * * * * * * * * * * * * * * * * * REMOVE SUBMATRIX * * *
05320 760 IF ((NW.NE.3).OR.(NF.NE.4)) GO TO 995
05330 NAMB=WORD(3)
05340 J=FNUMB(1)
05350 K=FNUMB(2)
05360 NR=FNUMB(3)
05370 NC=FNUMB(4)
05380 NTAG=0
05390 CALL MLIST(NAMB)
05400 IF(ERROR.GT.0) GO TO 990
05410 770 NB=INCA
05420 NRB=NR
05430 NCB=NC
05440 NAMA=WORD(2)
05450 CALL MFIND(NAMA)
05460 IF(ERROR.GT.0) GO TO 990
05470 CALL STOSM(NB,NRB,NCB,J,K,NTAG)
05480 IF(ERROR.GT.0) GO TO 990
05490 GO TO 5
05500 C * * * * * DELETION OF NR ROWS AND NC COLUMNS FROM MATRIX *
05510 780 IF(.NOT.(NW.EQ.3.AND.(NF.EQ.1.OR.NF.EQ.2))) GO TO 995
05520 NAMA=WORD(2)
05530 NAMB=WORD(3)
05540 NRD=FNUMB(1)
05550 IF(NF.EQ.2) GO TO 790
05560 K=1
05570 NCD=NRD
05580 GO TO 795
05590 790 NCD=FNUMB(2)
05600 K=0
05610 795 IF((NRD.GT.100).OR.(NCD.GT.100)) GO TO 901
05620 CALL DELRC(NAMA,NRD,NCD,NAMB,K)
05630 IF(ERROR) 990,5,990
05640 C * * * * * * * * SELECT EVERY NTH VALUE IN A ROW MATRIX
05650 800 IF ((NW.NE.3).OR.(NF.LE.0)) GO TO 995
05660 N=FNUMB(1)
05670 M=FNUMB(2)
05680 NAMA=WORD(2)
05690 NAMB=WORD(3)
05700 IF(N.LE.0.OR.M.LE.0) GO TO 908
05710 CALL MFIND(NAMA)
05720 IF(ERROR.GT.0) GO TO 990
05730 NA=INCA+M
05740 IF(NR-1) 901,813,811
05750 811 IF (NC.NE.1) GO TO 901
05760 NTMP=NR
05770 NR=NC
05780 NC=NTMP
05790 813 NC=1+(NC-M)/N
05800 IF(NC.LT.M) GO TO 901
05810 CALL MLIST(NAMB)
05820 IF(ERROR.GT.0) GO TO 990
05830 DO 815 I=1,NC
05840 A(I+INCA)=A(NA)
05850 NA=NA+N
05860 815 CONTINUE
05870 GO TO 5
05880 C = = = = = = EIGENVALUE,RESPONSE,AND SUPPORTING OPERATIONS = = =
05890 C * * * * CHOLESKI RIGHT DECOMPOSITION (UPPER TRIANGULAR) * * *
05900 820 L=3
05910 IF (NW.NE.3) GO TO 995
05920 NAMB=WORD(3)
05930 GO TO 441
05940 822 CALL MLIST(NAMB)
05950 IF(ERROR.GT.0) GO TO 990
05960 NB=INCA
05970 CALL MFIND(NAMA)
05980 CALL DECOM(NB)
05990 IF(ERROR) 990,5,990
06000 C * * * * * * * * * * * * * * * * * * * CHOLESKI REDUCTION * * :
06010 840 IFLAG=0
06020 841 IF(.NOT.(NW.EQ.2.OR.NW.EQ.3)) GO TO 995
06030 SCAL1=FNUMB(1)
06040 NAMA=WORD(2)
06050 CALL MFIND(NAMA)
06060 IF(ERROR.GT.0) GO TO 990
06070 NA=INCA
06080 IF (NR.NE.NC) GO TO 945
06090 NRA=NR
06100 IF(NW.EQ.2) GO TO 845
06110 NAMB=WORD(3)
06120 CALL MFIND(NAMB)
06130 IF(ERROR.GT.0) GO TO 990
06140 IF (NRA.NE.NR) GO TO 901
06150 GO TO 846
06160 845 NC=0
06170 846 CALL CHSOLV(NA,IFLAG,SCAL1)
06180 IF(ERROR)5,5,990
06190 C * * * * * * * * CHOLESKI BACK SUBSITUTION * * *
06200 860 IFLAG=1
06210 GO TO 841
06220 C * * * * * * * EIGENVALUES AND EIGENVECTORS * * * * * * * * * * *
06230 880 L=4
06240 IF ((NW.NE.5).OR.(NF.LE.0)) GO TO 995
06250 M=FNUMB(1)
06260 SCAL1=FNUMB(2)
06270 NAMA=WORD(2)
06280 GO TO 442
06290 885 N=NR
06300 NR=1
06310 NC=N
06320 NAMD=WORD(5)
06330 CALL MLIST(NAMD)
06340 ND=INCA
06350 IF(ERROR.GT.0) GO TO 990
06360 NR=IABS(M)
06370 NAMC=WORD(4)
06380 CALL MLIST(NAMC)
06390 IF(ERROR.GT.0) GO TO 990
06400 NE=INCA
06410 CALL MFIND(NAMA)
06420 NA=INCA
06430 NAMB=WORD(3)
06440 CALL MFIND(NAMB)
06450 IF(ERROR.GT.0) GO TO 990
06460 IF (NR.NE.1) GO TO 901
06470 IF (N.NE.NC) GO TO 901
06480 CALL EIGEN(A(NA+1),A(NE+1),A(ND+1),N,M,N)
06490 IF(ERROR) 5,5,990
06500 C * * * * * * * * FUNCTION GENERATION * * * * * * * * *
06510 1000 IF ((NW.NE.3).OR.(NF.NE.2)) GO TO 995
06520 NAMA=WORD(2)
06530 NAMB=WORD(3)
06540 N=FNUMB(1)
06550 SCAL1=FNUMB(2)
06560 NR=1
06570 NC=N
06580 CALL MLIST(NAMB)
06590 IF(ERROR.GT.0) GO TO 990
06600 NB=INCA
06610 CALL MFIND(NAMA)
06620 IF(ERROR.GT.0) GO TO 990
06630 IF (NC.NE.2) GO TO 901
06640 NS=NR*NC
06650 CALL FUNGN(NB,N,NS,SCAL1)
06660 IF(ERROR)5,5,990
06670 C * * * * * * * * * RESPONSE* * * * * * * * * * * * * * * * *
06680 1020 IF ((NW.NE.5).OR.(NF.NE.4)) GO TO 995
06690 NAMA=WORD(2)
06700 NAMB=WORD(3)
06710 NAMC=WORD(4)
06720 NAMD=WORD(5)
06730 N=FNUMB(1)
06740 K=FNUMB(2)
06750 SCAL1=FNUMB(3)
06760 SCAL2=FNUMB(4)
06770 IF((K.LT.0).OR.(K.GT.2)) GO TO 995
06780 CALL MFIND(NAMA)
06790 IF(ERROR.GT.0) GO TO 990
06800 NA=INCA
06810 M=NC
06820 IF (NR.NE.1) GO TO 901
06830 CALL MFIND(NAMC)
06840 IF(ERROR.GT.0) GO TO 990
06850 L=NC
06860 NTAG=1
06870 IF (NR.EQ.1) NTAG=0
06880 NC=(NC-1)/N
06890 NR=M
06900 NRA=INCA
06910 CALL MLIST(NAMD)
06920 IF(ERROR.GT.0) GO TO 990
06930 ND=INCA
06940 CALL MFIND(NAMB)
06950 IF(ERROR.GT.0) GO TO 990
06960 IF ((NC.NE.1).OR.(NR.NE.M)) GO TO 901
06970 CALL RESPON(A(NA+1),A(NRA+1),A(ND+1),M,N,L,SCAL1,SCAL2,NTAG,K)
06980 IF(ERROR) 5,5,990
06990 C * * * * * * * * * * * * * LIST NAMES OF MATRICES
07000 2000 NS=NSIZE
07010 TYPE 2001, NUM,NS
07020 2001 FORMAT(1X,I3,' MATRICES DEFINED, USING A TOTAL OF',
07030 1 I6,' WORDS OF STORAGE'/)
07040 IF (NUM.EQ.0) GO TO 5
07050 DO 2010 I=1,NUM
07060 NAMA=NAME(I+INC1)
07070 CALL MFIND(NAMA)
07080 TYPE 2002,NAMA,NR,NC
07090 2002 FORMAT(2X,A5,'(',I3,',',I3,')')
07100 2010 CONTINUE
07110 TYPE 2012
07120 2012 FORMAT(/)
07130 GO TO 5
07140 C * * * * * * * * * * * GET MATRIX FILE * * * * * * * * * * *
07150 2020 IF ((NW.NE.2).OR.(NF.NE.0)) GO TO 995
07160 CALL MGET(WORD(2))
07170 IF (ERROR.NE.0) GO TO 990
07180 GO TO 5
07190 C * * * * * * * * * SAVE ALL MATRICES * * * * * * * * * * * * *
07200 2040 IF ((NW.EQ.3).AND.(NF.EQ.0)) GO TO 2050
07210 IF ((NW.EQ.2).AND.(NF.EQ.1)) GO TO 2042
07220 GO TO 995
07230 C * * * * * * * * SAVE "N" MATRICES * * * * * * * * * * *
07240 2042 NV=FNUMB(1)
07250 CALL MSAVE(WORD(2),NV)
07260 IF (ERROR.NE.0) GO TO 990
07270 GO TO 5
07280 2050 IF (WORD(2).NE.'ALL') GO TO 995
07290 NV=NUM
07300 CALL MSAVE(WORD(3),NV)
07310 GO TO 5
07320 C * * * * * * * * GENERATION OF FIXED END MOMENTS MATRIX * * *
07330 2060 IF(NW.NE.2.OR.NF.GT.1) GO TO 995
07340 IF(NF.EQ.0) FNUMB(1)=1.
07350 NAMA=WORD(2)
07360 NC=1
07370 NTIME=FNUMB(1)
07380 NR=4*NTIME
07390 CALL MLIST(NAMA)
07400 IF(ERROR.GT.0) GO TO 990
07410 CALL FIXEM(NTIME)
07420 GO TO 5
07430 C * * * * * * * * * * * RANK OF MATRIX * * * * * * * * * *
07440 2080 IF(NW.NE.2) GO TO 995
07450 NAMA=WORD(2)
07460 SCAL1=-99999.
07470 GO TO 482
07480 C * * * * * * * * * * END OF COMMANDS FOR NOW * * : : * * *
07490 C * * * * * * * * * * * *ERROR MESSAGES * * * * *
07500 901 ERROR=1
07510 GO TO 990
07520 903 ERROR=2
07530 GO TO 990
07540 908 ERROR=8
07550 GO TO 990
07560 909 ERROR=9
07570 GO TO 990
07580 945 TYPE 946, NAMA
07590 946 FORMAT(2X,' MATRIX ',A5,' IS NOT SQUARE'/)
07600 GO TO 5
07610 990 CALL SUBERR
07620 GO TO 5
07630 995 TYPE 996
07640 996 FORMAT(2X,' ERROR IN COMMAND STRING'/)
07650 GO TO 5
07660 998 TYPE 999
07670 999 FORMAT(2X,' ERROR IN MATRIX LOGIC, COMMAND STOPPED'/)
07680 GO TO 5
07690 END
00010 C***********************************************************************
00020 C FREE FORMAT READ SUBROUTINE
00030 C***********************************************************************
00040 SUBROUTINE FREAD(NA,ALPHA,NFLT,FNUMB)
00050 COMMON /FF/ NFMX,IPT,NWMX
00060 DIMENSION LINE(73),ALPHA(1),FNUMB(1)
00070 INTEGER A,BLANK,E,POINT,Z,COMMA,PLUS
00080 INTEGER ALPHA, ZERO
00090 LOGICAL DIGIT
00100 DATA A/'A'/
00110 DATA Z/'Z'/
00120 DATA E/'E'/
00130 DATA BLANK/' '/
00140 DATA POINT/'.'/
00150 DATA MINUS/'-'/
00160 DATA COMMA/','/
00170 DATA ZERO/'0'/
00180 DATA NINE/'9'/
00190 DATA PLUS/'+'/
00200 DIGIT(N) = N .GE. ZERO .AND. N .LE. NINE
00210 IF (NFMX.EQ.0) NFMX=100
00220 IF (NWMX.EQ.0) NWMX=5
00230 8 NI=0
00240 NF=0
00250 10 I=0
00260 IF(NA.GE.0.OR.IPT.GT.0) GO TO 120
00270 TYPE 111
00280 111 FORMAT('+#'$)
00290 120 CALL BCOUT(LINE,73)
00300 ACCEPT 20, LINE
00310 20 FORMAT(73A1)
00320 1 SIGN=1.
00330 PSIGN=1.0
00340 X=0
00350 15 I=I+1
00360 IC=LINE(I)
00370 IF(IC.EQ.BLANK.OR.IC.EQ.COMMA) GO TO 1000
00380 IF(IC.GE.A.AND.IC.LE.Z) GO TO 100
00390 IF(DIGIT(IC)) GO TO 300
00400 IF(IC.EQ.MINUS) GO TO 200
00410 IF(IC.EQ.POINT) GO TO 350
00420 11 IF(IPT.GT.0) GO TO 15
00430 GO TO 1003
00440 C ALPHA CHARACTER ENVOUNTERED - BUILD WORD (S)
00450 100 IF(NA) 11,11,101
00460 101 J=I
00470 NI=NI+1
00480 102 I=I+1
00490 IC=LINE(I)
00500 IF(IC.EQ.BLANK.OR.IC.EQ.COMMA) GO TO 110
00510 GO TO 102
00520 110 K=I-1
00530 JJ = 1
00540 DO 115 II = J,K
00550 IF ((NI+JJ/5).GT.(NWMX+1)) GO TO 116
00560 CALL PUT(ALPHA(NI),JJ,LINE(II))
00570 115 JJ = JJ+1
00580 116 NI=NI+(K-J)/5
00590 IF(ALPHA(1).EQ.'REMAR') GO TO 999
00600 GO TO 1000
00610 C MINUS ENCOUNTERED - SET SIGN NEGATIVE
00620 200 SIGN=-1.
00630 I=I+1
00640 IC=LINE(I)
00650 IF(DIGIT(IC)) GO TO 300
00660 IF(IC.EQ.POINT) GO TO 350
00670 GO TO 11
00680 300 IF(NFLT) 11,301,301
00690 C DIGIT ENCOUNTERED - BUILD NUMBER
00700 301 J=I
00710 302 I=I+1
00720 IC=LINE(I)
00730 IF(DIGIT(IC)) GO TO 302
00740 IF(IC.EQ.POINT) GO TO 324
00750 IF(.NOT.(IC.EQ.BLANK.OR.IC.EQ. COMMA)) GO TO 11
00760 IGOTO=1
00770 GO TO 325
00780 324 IGOTO=2
00790 325 K=I-1
00800 CALL ATSI(LINE,J,K,X)
00810 X=X*SIGN
00820 GO TO (900,400),IGOTO
00830 350 X=0
00840 I=I+1
00850 IC=LINE(I)
00860 GO TO 500
00870 400 I=I+1
00880 IC=LINE(I)
00890 IF(IC.EQ.BLANK) GO TO 900
00900 IF(IC.EQ.E) GO TO 700
00910 500 IF(.NOT.DIGIT(IC)) GO TO 11
00920 J=I
00930 502 I=I+1
00940 IC=LINE(I)
00950 IF(DIGIT(IC)) GO TO 502
00960 IF(IC.EQ.E) GO TO 524
00970 IF(.NOT.(IC.EQ.BLANK.OR.IC.EQ.COMMA)) GO TO 11
00980 IGOTO=1
00990 GO TO 525
01000 524 IGOTO=2
01010 525 K=I-1
01020 CALL ATSI(LINE,J,K,DEC)
01030 X=X+SIGN*(DEC/10.**(K-J+1))
01040 GO TO (900,700),IGOTO
01050 700 I=I+1
01060 IC=LINE(I)
01070 IF(IC.EQ.MINUS) GO TO 701
01080 IF(IC.NE.PLUS) GO TO 720
01090 I=I+1
01100 IC=LINE(I)
01110 GO TO 720
01120 701 PSIGN=-1.
01130 I=I+1
01140 IC=LINE(I)
01150 720 J=I
01160 IF(IC.EQ.'+') GO TO 722
01170 IF(.NOT.DIGIT(IC)) GO TO 11
01180 722 I=I+1
01190 IC=LINE(I)
01200 IF(DIGIT(IC)) GO TO 722
01210 IF(.NOT.(IC.EQ.BLANK.OR.IC.EQ.COMMA)) GO TO 11
01220 K=I-1
01230 IF(K-J-1) 725,725,11
01240 725 CALL ATSI(LINE,J,K,PW)
01250 NP=PW*PSIGN
01260 X=X*10.**NP
01270 900 NF=NF+1
01280 IF (NF.GT.NFMX) GO TO 1000
01290 FNUMB(NF)=X
01300 1000 IF(I-70) 1,1,1001
01310 1001 IF(NF-NFLT) 1002,999,1005
01320 1002 IF(IPT.GT.0) GO TO 10
01330 1003 NLT=NFLT-NF
01340 TYPE 1004, NF,NLT
01350 1004 FORMAT(1X,I4,' VALUES STORED ** ENTER',I4,/)
01360 GO TO 10
01370 1005 IF(NA.GE.0.OR.IPT.GT.0) GO TO 999
01380 NLT=NF-NFLT
01390 TYPE 1006, NLT
01400 1006 FORMAT(1X,I3,' TOO MANY VALUES ENTER & OMITTED'/)
01410 999 NFLT=NF
01420 NA=NI
01430 NFMX=0
01440 NWMX=0
01450 RETURN
01460 END
00010 SUBROUTINE ATSI(IN,IBEG,IEND,VALUE)
00020 DIMENSION IN(73),NUM(10)
00030 DATA(NUM(I),I=1,5)/'1','2','3','4','5'/
00040 DATA(NUM(I),I=6,10)/'6','7','8','9','0'/
00050 VALUE=0
00060 IF(IEND.LT.IBEG) GO TO 200
00070 DO 100 I=IBEG,IEND
00080 IC=IN(I)
00090 DO 10 J=1,10
00100 IF(.NOT.IC.EQ.NUM(J)) GO TO 10
00110 IF(J-10)9,8,8
00120 8 J=0
00130 9 VALUE=VALUE*10+J
00140 GO TO 100
00150 10 CONTINUE
00160 100 CONTINUE
00170 GO TO 999
00180 200 VALUE=0
00190 999 RETURN
00200 END
00010 C**********************************************************
00020 C SUBROUTINE TO BLANK OUT INPUT ARRAY
00030 C*************************************************************
00040 SUBROUTINE BCOUT(LINE,IC)
00050 INTEGER BLANK
00060 DIMENSION LINE(1)
00070 DATA BLANK/' '/
00080 DO 100 I=1,IC
00090 LINE(I)=BLANK
00100 100 CONTINUE
00110 RETURN
00120 END
00010 SUBROUTINE MFIND(MM)
00020 C***********************************************************************
00030 C SUBROUTINE TO FIND LOCATION AND SIZE OF MATRIX MM
00040 C***********************************************************************
00050 INTEGER ERROR
00060 COMMON /XX/ ERROR
00070 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00080 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1),NME(4),NSIZE
00090 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00100 I=1
00110 5 IF (NUM-I) 40,10,10
00120 10 IF (NAME(I+INC1)-MM) 20,30,20
00130 20 I=I+1
00140 GO TO 5
00150 30 NR=NROW(I+INC3)
00160 NC=NCOL(I+INC2)
00170 NX=NN(I+INC4)
00180 INCA=NX-LOC(A)
00190 GO TO 101
00200 40 ERROR=2
00210 101 RETURN
00220 END
00010 SUBROUTINE MLIST(MM)
00020 C***********************************************************************
00030 C SUBROUTINE TO LIST MATRIX
00040 C***********************************************************************
00050 INTEGER ERROR,FREZ
00060 COMMON /XX/ ERROR
00070 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00080 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1),NME(4),NSIZE
00090 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00100 I=1
00110 5 IF (NUM-I) 30,10,10
00120 10 IF (NAME(I+INC1)-MM) 20,25,20
00130 20 I=I+1
00140 GO TO 5
00150 C***********************************************************************
00160 C CHECK FOR POSSIBLE MULTIPLY AND ADD
00170 C***********************************************************************
00180 25 IF(MPYTAG) 52,29,26
00190 26 TYPE 100
00200 IF(NR-NROW(I+INC3)) 54,27,54
00210 27 IF(NC-NCOL(I+INC2)) 54,28,54
00220 28 NX=NN(I+INC4)
00230 GO TO 60
00240 29 ERROR=4
00250 GO TO 60
00260 30 NUM=NUM+1
00270 IF (NUM-101) 45,40,45
00280 40 ERROR=6
00290 GO TO 60
00300 45 NSIZE=NSIZE+NR*NC
00310 C SET LIMIT TO 30K FOR MATRXB
00320 C IF(30000-NSIZE) 50,65,65
00330 IF (9000-NSIZE) 50,65,65
00340 50 ERROR=5
00350 NSIZE=NSIZE-NR*NC
00360 NUM=NUM-1
00370 GO TO 60
00380 65 IF (NUM.LE.NNUM) GO TO 31
00390 CALL MSTOR
00400 31 NX=FREZ(NR*NC)
00410 NROW(NUM+INC3)=NR
00420 NCOL(NUM+INC2)=NC
00430 NAME(NUM+INC1)=MM
00440 NN(NUM+INC4)=NX
00450 INCA=NX-LOC(A)
00460 GO TO 60
00470 MPYTAG=0
00480 52 ERROR=8
00490 GO TO 60
00500 54 ERROR=1
00510 GO TO 60
00520 60 RETURN
00530 100 FORMAT(' THE MULTIPLICATION RESULT WILL BE ADDED TO THE MATRIX',
00540 1' IN FIELD C'///)
00550 END
00010 C***********************************************************************
00020 C SUBROUTINE TO FIND MAXIMUM VALUES IN ROWS ( OR COLUMNS) OF MATRIX
00030 C***********************************************************************
00040 SUBROUTINE ENVEL(NA,NBB,SCAL1,KCOL)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ A(1)
00070 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00080 DO 80 I=1,NR
00090 AMAX=0.0
00100 TIM =SCAL1
00110 IF(KCOL) 10,10,15
00120 10 K=(I-1)*NC+1
00130 GO TO 20
00140 15 K=I
00150 20 DO 50 J=1,NC
00160 XA=ABS(A(NA+K))
00170 IF(AMAX-XA) 25,30,30
00180 25 AMAX=XA
00190 TMAX=TIM
00200 30 IF(KCOL) 35,35,40
00210 35 K=K+1
00220 GO TO 50
00230 40 K=K+NR
00240 50 TIM = TIM + SCAL1
00250 IF(NBB) 60,60,55
00260 55 A(INCA+I)=AMAX
00270 60 IF(KCOL) 65,65,70
00280 65 TYPE 100, I, TMAX, AMAX
00290 GO TO 80
00300 70 TYPE 101, I, TMAX, AMAX
00310 80 CONTINUE
00320 TYPE 102
00330 102 FORMAT(/)
00340 RETURN
00350 100 FORMAT(' ROW NUMBER',I4,3X,'INTERVAL = ',1PE12.5,3X,
00360 1 'MAXIMUM VALUE = ',1PE12.5)
00370 101 FORMAT(14H COLUMN NUMBER,I4,2X,10HINTERVAL = ,1PE11.4,2X,16HMAXIMU
00380 1M VALUE = ,1PE12.5)
00390 END
00010 SUBROUTINE DELETE(MM)
00020 C***********************************************************************
00030 C SUBROUTINE TO DELETE MATRIX MM
00040 C***********************************************************************
00050 INTEGER ERROR,FRET
00060 COMMON /XX/ ERROR
00070 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00080 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1),NME(4),NSIZE
00090 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00100 I=1
00110 5 IF (NUM-I) 60,10,10
00120 10 IF ( NAME(I+INC1)-MM) 20,30,20
00130 20 I=I+1
00140 GO TO 5
00150 30 NX=NN(I+INC4)
00160 NUM=NUM-1
00170 NS=NROW(I+INC3)*NCOL(I+INC2)
00180 NSIZE=NSIZE-NS
00190 IF (NUM.EQ.0) GO TO 50
00200 IF (I.GT.NUM) GO TO 50
00210 DO 40 J=I,NUM
00220 NAME(J+INC1)=NAME(J+INC1+1)
00230 NROW(J+INC3)=NROW(J+INC3+1)
00240 NCOL(J+INC2)=NCOL(J+INC2+1)
00250 40 NN(J+INC4)=NN(J+INC4+1)
00260 50 CALL FRET(NS,NX)
00270 90 RETURN
00280 60 ERROR=2
00290 GO TO 90
00300 END
00010 C***********************************************************************
00020 C MATRIX TRANSPOSE
00030 C***********************************************************************
00040 SUBROUTINE TRANS(NB)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ A(1)
00070 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00080 K=1
00090 DO 100 I=1,NR
00100 L=I
00110 DO 100 J=1,NC
00120 A(K+INCA)=A(L+NB)
00130 L=L+NR
00140 100 K=K+1
00150 RETURN
00160 END
00010 C***********************************************************************
00020 C MATRIX MULTIPLICATION
00030 C***********************************************************************
00040 SUBROUTINE MULT(NA,NB,NCA,NCB,ITRMPY)
00050 INTEGER ERROR
00060 COMMON /X/ ERROR
00070 COMMON /YY/ INCA
00080 COMMON /Z/ A(1)
00090 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00100 K=1
00110 DO 200 I=1,NR
00120 DO 200 J=1,NCB
00130 MM=J
00140 C***********************************************************************
00150 C CHECK FOR POSSIBLE PREMULTIPLICATION BY TRANSPOSE
00160 C***********************************************************************
00170 IF(ITRMPY) 998,30,40
00180 30 NN=(I-1)*NCA+1
00190 GO TO 50
00200 40 NN=I
00210 C***********************************************************************
00220 C CHECK FOR POSSIBLE MULTIPLY AND ADD
00230 C***********************************************************************
00240 50 IF(MPYTAG) 998,60,70
00250 60 A(INCA+K)=0.0
00260 70 DO 100 L=1,NCA
00270 A(INCA+K)=A(INCA+K)+A(NA+NN)*A(NB+MM)
00280 IF(ITRMPY) 998,80,90
00290 80 NN=NN+1
00300 GO TO 100
00310 90 NN=NN+NR
00320 100 MM=MM+NCB
00330 200 K=K+1
00340 GO TO 999
00350 998 ERROR=8
00360 999 RETURN
00370 END
00010 SUBROUTINE SUBERR
00020 C***********************************************************************
00030 C ERROR PRINT OUT SUBROUTINE
00040 C***********************************************************************
00050 INTEGER ERROR
00060 COMMON /XX/ ERROR
00070 COMMON /ZZ/ NUM
00080 DATA NCT/0/
00090 IF (ERROR.EQ.0) GO TO 5
00100 GO TO (1,2,3,4,5,6,7,8,9,10),ERROR
00110 1 TYPE 101
00120 GO TO 100
00130 2 TYPE 102
00140 GO TO 100
00150 3 TYPE 103
00160 GO TO 100
00170 4 TYPE 104
00180 GO TO 100
00190 5 TYPE 105
00200 IF (NCT.NE.0) GO TO 60
00210 NCT=NCT+1
00220 TYPE 113
00230 TYPE 111
00240 CALL PSIZE(24)
00250 GO TO 100
00260 60 TYPE 112
00270 STOP
00280 6 TYPE 106
00290 GO TO 100
00300 7 TYPE 107
00310 GO TO 100
00320 8 TYPE 108
00330 GO TO 100
00340 9 TYPE 109
00350 GO TO 100
00360 10 TYPE 110
00370 100 RETURN
00380 101 FORMAT (2X' MATRICES INCOMPATIBLE'/)
00390 102 FORMAT (2X' MATRIX UNDEFINED'/)
00400 103 FORMAT (2X' OPERATION UNDEFINED'/)
00410 104 FORMAT (2X' MATRIX PREVIOUSLY DEFINED'/)
00420 105 FORMAT (2X'STORAGE EXCEEDED'/)
00430 106 FORMAT (2X' OVER 100 MATRICES IN CORE'/ )
00440 107 FORMAT (27H ILLEGAL TAPE USED )
00450 108 FORMAT (2X' STRICTLY POSITIVE NUMBER HAS TAKEN ON NEGATIVE',
00460 1 ' VALUE'/)
00470 109 FORMAT (2X' SCALAR UNDEFINED'/)
00480 110 FORMAT(2X,' INVALID ARITHMETIC OPERATION ATTEMPTED'/)
00490 111 FORMAT(2X,'PROTECT MATRICES WITH A "SAVE" '/)
00500 112 FORMAT(2X,'ERROR EXIT'/)
00510 113 FORMAT(2X,'DELETE ALL MATRICES NOT IN USE IMMEDIATELY!')
00520 END
00010 C***********************************************************************
00020 C PRINT OPERATION
00030 C***********************************************************************
00040 SUBROUTINE PRINT(NAME)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ A(1)
00070 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00080 DIMENSION NCOLL(5)
00090 TYPE 102,NAME,NR,NC
00100 DO 50 M=1,NC,5
00110 NN=NC-M
00120 IF (NN-4) 40,40,35
00130 35 NN=4
00140 40 MM=NN+1
00150 DO 45 N=1,MM
00160 45 NCOLL(N)=M-1+N
00170 TYPE 101, (NCOLL(N),N=1,MM)
00180 DO 50 N=1,NR
00190 NL=M+(N-1)*NC
00200 NH=NL+NN
00210 50 TYPE 100, N,(A(I+INCA),I=NL,NH)
00220 TYPE 103
00230 RETURN
00240 100 FORMAT(I4,5(1X,1PE12.5))
00250 101 FORMAT(/1X,I11,5I13/)
00260 102 FORMAT(2X,A5,'(',I3,',',I3,')')
00270 103 FORMAT(//)
00280 END
00010 C*************************************************************
00020 C SUBROUTINE TO LOAD MATRIX
00030 C*****************************************************************
00040 SUBROUTINE LOAD
00050 COMMON /FF/ NFMX,ITAPE,NWMX
00060 COMMON /YY/ INCA
00070 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1),NME(4),NSIZE
00080 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00090 TYPE 10
00100 10 FORMAT(2X,' ENTER MATRIX BY ROWS'/)
00110 K=1
00120 DO 100 J=1,NR
00130 NA=-NC
00140 NF=NC
00150 NWMX=1
00160 NFMX=NC
00170 CALL FREAD(NA,NA,NF,A(K+INCA))
00180 K=K+NC
00190 100 CONTINUE
00200 RETURN
00210 END
00010 C***********************************************************************
00020 C SUBROUTINE TO PRINT MAP OF NON-ZERO ELEMENTS IN A MATRIX
00030 C***********************************************************************
00040 SUBROUTINE MTXMAP(SCAL1)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ A(1)
00070 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00080 DIMENSION NCOLL(20)
00090 DATA NULL/'.'/
00100 DATA NZERO/'*'/
00110 DO 50 M=1,NC,20
00120 NN=NC-M
00130 IF(NN-19) 40,40,35
00140 35 NN=19
00150 40 MM=NN+1
00160 DO 41 N=1,MM
00170 41 NCOLL(N)=M-1+N
00180 TYPE 101, (NCOLL(N),N=1,MM)
00190 DO 50 I=1,NR
00200 NL=M+(I-1)*NC
00210 DO 45 J=1,MM
00220 ND=NL+J-1
00230 IF(ABS(A(INCA+ND))-SCAL1) 42,42,43
00240 42 NCOLL(J)=NULL
00250 GO TO 45
00260 43 NCOLL(J)=NZERO
00270 45 CONTINUE
00280 50 TYPE 100, I,(NCOLL(J),J=1,MM)
00290 TYPE 102
00300 RETURN
00310 100 FORMAT(I5,20(2X,A1))
00320 101 FORMAT(1H0,4X,20I3)
00330 102 FORMAT(//)
00340 END
00010 C***********************************************************************
00020 C SUBROUTINE TO PERFORM OPERATION ON EACH ELEMENT OF A MATRIX
00030 C***********************************************************************
00040 SUBROUTINE ELEMNT(IFLAG)
00050 INTEGER ERROR
00060 COMMON /XX/ ERROR
00070 COMMON /YY/ INCA
00080 COMMON /Z/ A(1)
00090 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00100 NS=NR*NC
00110 IF (IFLAG-2) 1,2,3
00120 1 DO 11 I=1,NS
00130 11 IF (A(I+INCA).LT.0) GO TO 13
00140 GO TO 3
00150 2 DO 21 I=1,NS
00160 21 IF (A(INCA+I).EQ.0) GO TO 13
00170 3 DO 10 I=1,NS
00180 IF(IFLAG-2) 4,6,8
00190 4 A(I+INCA)=SQRT(A(I+INCA))
00200 GO TO 10
00210 6 A(I+INCA)=1.0/A(I+INCA)
00220 GO TO 10
00230 8 A(I+INCA)=ALOG(A(I+INCA))
00240 10 CONTINUE
00250 GO TO 15
00260 13 ERROR=10
00270 15 RETURN
00280 END
00010 C***********************************************************************
00020 C SYMMETRICAL MATRIX INVERSION
00030 C***********************************************************************
00040 SUBROUTINE SYMINV(NA,SINGLR)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ A(1)
00070 COMMON /ZZ/ NUM,NNUM,NX,NRA,NCB
00080 INTEGER FREZ,FRET
00090 DIMENSION C(1),D(1)
00100 C************************************************************
00110 C GET FREE STORAGE FOR C AND D ARRAYS
00120 C****************************************************************
00130 LCPT=FREZ(2*NRA)
00140 LC=LCPT-LOC(C)
00150 LD=LCPT+NRA-LOC(D)
00160 DO 140 I=1,NRA
00170 NR=(I-1)*NRA
00180 DO 100 J=1,NRA
00190 K=NR+J
00200 100 C(J+LC)=A(K+NA)
00210 DMAX=C(I+LC)
00220 IF(SINGLR-ABS(DMAX)) 105,990,990
00230 C***********************************************************************
00240 C DIVIDE ROW BY ELEMENT ON DIAGONAL
00250 C***********************************************************************
00260 105 DO 110 J=1,NRA
00270 110 D(J+LD)=-C(J+LC)/DMAX
00280 IF(NCB) 998,116,112
00290 112 NLB=(I-1)*NCB+1
00300 NHB=NLB+NCB-1
00310 DO 115 K=NLB,NHB
00320 115 A(K+INCA)=A(K+INCA)/DMAX
00330 116 L=1
00340 C***********************************************************************
00350 C REDUCE REMAINING ROWS AND COLUMNS
00360 C***********************************************************************
00370 DO 125 J=1,NRA
00380 M=L
00390 DO 120 K=J,NRA
00400 A(L+NA)=A(L+NA)+C(J+LC)*D(K+LD)
00410 A(M+NA)=A(L+NA)
00420 M=M+NRA
00430 120 L=L+1
00440 125 L=L+J
00450 IF(NCB) 998,139,130
00460 130 LB=1
00470 TAG=-1.0
00480 DO 135 J=1,NRA
00490 IF(I-J) 132,131,132
00500 131 LB=LB+NCB
00510 TAG=1.0
00520 GO TO 135
00530 132 DO 133 K=NLB,NHB
00540 A(LB+INCA)=A(LB+INCA)-TAG*C(J+LC)*A(K+INCA)
00550 133 LB=LB+1
00560 135 CONTINUE
00570 139 D(I+LD)=-1.0/DMAX
00580 M=I
00590 C***********************************************************************
00600 C REPLACE ELEMENTS OF A WITH CALCULATED INVERSE ELEMENTS
00610 C***********************************************************************
00620 DO 140 J=1,NRA
00630 K=NR+J
00640 A(K+NA)=D(J+LD)
00650 A(M+NA)=D(J+LD)
00660 140 M=M+NRA
00670 NS=NRA*NRA
00680 DO 150 J=1,NS
00690 150 A(J+NA)=-A(J+NA)
00700 GO TO 999
00710 990 TYPE 991
00720 991 FORMAT(' AN ACCURATE SINGLE PRECISION INVERSE CANNOT',
00730 1 ' BE COMPUTED'/)
00740 GO TO 999
00750 998 ERROR=8
00760 999 CALL FRET(2*NRA,LCPT)
00770 RETURN
00780 END
00010 C***********************************************************************
00020 C SUBROUTINE TO PERFORM VARIOUS DIAGONAL MATRIX MULTIPLICATIONS
00030 C***********************************************************************
00040 SUBROUTINE DGMPY(NA,N1,N2,IFLAG)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ A(1)
00070 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00080 DO 10 I=1,N1
00090 DO 10 J=1,N2
00100 K=(I-1)*N2+J+NA
00110 IF(IFLAG-2) 4,6,8
00120 4 A(K)=A(K)*A(I+INCA)
00130 GO TO 10
00140 6 A(K)=A(K)*A(J+INCA)
00150 GO TO 10
00160 8 A(J+NA)=A(J+NA)*A(J+INCA)
00170 10 CONTINUE
00180 RETURN
00190 END
00010 C***********************************************************************
00020 C SUBROUTINE TO CHECK ON SYMMETRY OF A MATRIX
00030 C***********************************************************************
00040 SUBROUTINE SYMCHK(SCAL1,IFLAG)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ A(1)
00070 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00080 DO 10 I=1,NR
00090 DO 10 J=I,NR
00100 IJ=(I-1)*NR+J
00110 JI=(J-1)*NR+I
00120 TEST=ABS(A(INCA+IJ)-A(INCA+JI))
00130 IF(TEST-SCAL1) 10,10,5
00140 5 TYPE 100, I,J,A(INCA+IJ),J,I,A(INCA+JI)
00150 IFLAG=1
00160 10 CONTINUE
00170 RETURN
00180 100 FORMAT(3X,'(',I3,',',I3,') =',
00190 1 1PE14.7,5X,'(',I3,',',I3,') =',1PE14.7/)
00200 END
00010 C***********************************************************************
00020 C GENERAL MATRIX INVERSION SUBROUTINE
00030 C***********************************************************************
00040 SUBROUTINE INVERT(NA,SINGLR)
00050 INTEGER ERROR,FREZ,FRET
00060 COMMON /XX/ ERROR
00070 COMMON /YY/ INCA
00080 COMMON /Z/ A(1)
00090 COMMON /ZZ/ NUM,NNUM,NX,N,NCB
00100 DIMENSION M(1),C(1)
00110 INTEGER RANK
00120 LOGICAL IRANK
00130 IRANK=.FALSE.
00140 IF(SINGLR.GT.-99998.) GO TO 50
00150 IRANK=.TRUE.
00160 SINGLR=0.
00170 50 DET=1.
00180 MPT=FREZ(N)
00190 MC=MPT-LOC(M)
00200 LCPT=FREZ(N)
00210 LC=LCPT-LOC(C)
00220 DO 90 I=1,N
00230 90 M(I+MC)=-I
00240 DO 140 I=1,N
00250 RANK=I
00260 C***********************************************************************
00270 C LOCATE LARGEST ELEMENT
00280 C***********************************************************************
00290 D=0.0
00300 DO 105 L=1,N
00310 IF(M(L+MC)) 100,100,105
00320 100 J=L
00330 DO 104 K=1,N
00340 IF(M(K+MC)) 101,101,103
00350 101 IF (ABS(D)-ABS(A(J+NA))) 102,103,103
00360 102 LD=L
00370 KD=K
00380 D=A(J+NA)
00390 103 J=J+N
00400 104 CONTINUE
00410 105 CONTINUE
00420 IF(I-1) 998,107,110
00430 107 IF(SINGLR) 998,108,109
00440 108 SINGLR=ABS(D*1.0E-06)
00450 109 TYPE 992, SINGLR
00460 110 IF(SINGLR-ABS(D)) 111,990,990
00470 C***********************************************************************
00480 C INTERCHANGE COLUMNS
00490 C***********************************************************************
00500 111 TEMP=-M(LD+MC)
00510 M(LD+MC)=M(KD+MC)
00520 M(KD+MC)=TEMP
00530 L=LD
00540 K=KD
00550 DO 114 J=1,N
00560 C(J+LC)=A(L+NA)
00570 A(L+NA)=A(K+NA)
00580 A(K+NA)=C(J+LC)
00590 L=L+N
00600 114 K=K+N
00610 C***********************************************************************
00620 C DIVIDE ROW BY LARGEST ELEMENT
00630 C***********************************************************************
00640 DET=DET*D
00650 NLA=(KD-1)*N+1
00660 NHA=NLA+N-1
00670 DO 115 K=NLA,NHA
00680 115 A(K+NA)=-A(K+NA)/D
00690 IF(NCB) 998,120,116
00700 116 NLB=(KD-1)*NCB+1
00710 NHB=NLB+NCB-1
00720 DO 117 K=NLB,NHB
00730 117 A(K+INCA)=A(K+INCA)/D
00740 C***********************************************************************
00750 C REDUCE REMAINING ROWS AND COLUMNS
00760 C***********************************************************************
00770 120 L=1
00780 LB=1
00790 DO 137 J=1,N
00800 IF (J-KD) 130,125,130
00810 125 L=L+N
00820 LB=LB+NCB
00830 GO TO 137
00840 130 DO 131 K=NLA,NHA
00850 A(L+NA)=A(L+NA)+C(J+LC)*A(K+NA)
00860 131 L=L+1
00870 IF(NCB) 998,137,132
00880 132 IF(M(J+MC)) 134,134,133
00890 133 TAG=-1.0
00900 GO TO 135
00910 134 TAG= 1.0
00920 135 DO 136 K=NLB,NHB
00930 A(LB+INCA)=A(LB+INCA)-TAG*C(J+LC)*A(K+INCA)
00940 136 LB=LB+1
00950 137 CONTINUE
00960 C***********************************************************************
00970 C REDUCE COLUMN
00980 C***********************************************************************
00990 C(KD+LC)=1.0
01000 J=KD
01010 DO 139 K=1,N
01020 A(J+NA)=C(K+LC)/D
01030 139 J=J+N
01040 140 CONTINUE
01050 C***********************************************************************
01060 C INTERCHANGE ROWS
01070 C***********************************************************************
01080 DO 200 I=1,N
01090 L=0
01100 150 L=L+1
01110 IF(M(L+MC)-I) 150,160,150
01120 160 K=(L-1)*N+1
01130 J=(I-1)*N+1
01140 KB=(L-1)*NCB+1
01150 JB=(I-1)*NCB+1
01160 M(L+MC)=M(I+MC)
01170 M(I+MC)=I
01180 DO 170 LL=1,N
01190 TEMP=A(K+NA)
01200 A(K+NA)=A(J+NA)
01210 A(J+NA)=TEMP
01220 J=J+1
01230 170 K=K+1
01240 IF(NCB) 998,200,180
01250 180 DO 190 LL=1,NCB
01260 TEMP=A(KB+INCA)
01270 A(KB+INCA)=A(JB+INCA)
01280 A(JB+INCA)=TEMP
01290 JB=JB+1
01300 190 KB=KB+1
01310 200 CONTINUE
01320 IF(IRANK) GO TO 995
01330 TYPE 201, DET
01340 201 FORMAT(' DETERMINANT =',1PE15.7/)
01350 C***********************************************************************
01360 GO TO 999
01370 C***********************************************************************
01380 990 IF(IRANK) GO TO 994
01390 TYPE 991, KD,LD,D
01400 991 FORMAT(2X,' (',I3,',',I3,') =',1PE14.7,3X,'WILL NOT',
01410 1 ' PRODUCE'/' AN ACCURATE SINGLE PRECISION INVERSION'/)
01420 992 FORMAT (9H SCALAR= 1PE15.7/)
01430 GO TO 999
01440 994 RANK=RANK-1
01450 995 TYPE 996, RANK
01460 996 FORMAT(' RANK =',I5/)
01470 GO TO 999
01480 998 ERROR=8
01490 999 CALL FRET(N,MPT)
01500 CALL FRET(N,LCPT)
01510 RETURN
01520 END
00010 C***********************************************************************
00020 C SUBROUTINE TO FORM MEMBER STIFFNESS MATRICES
00030 C***********************************************************************
00040 SUBROUTINE FORMK(NTYPE,N)
00050 COMMON /FF/ NFMX,ITAPE,NWMX
00060 COMMON /YY/ INCA
00070 COMMON /Z/ A(1)
00080 DIMENSION DATA(5)
00090 NRC=NTYPE*N
00100 NS=NRC*NRC
00110 C***********************************************************************
00120 C INITIALIZE MATRIX, READ IN DATA, AND CALCULATE STIFFNESS CONSTANTS
00130 C***********************************************************************
00140 DO 1 I=1,NS
00150 1 A(INCA+I)=0.0
00160 DO 2 I=1,5
00170 DATA(I)=0.0
00180 2 CONTINUE
00190 TYPE 100
00200 100 FORMAT(' ENTER THE FOLLOWING DATA FOR EACH MEMBER'/)
00210 IF(NTYPE.EQ.3) GO TO 10
00220 TYPE 101
00230 101 FORMAT(1X,' I ,SHEAR AREA, L , E'/)
00240 NF=4
00250 GO TO 11
00260 10 TYPE 102
00270 102 FORMAT(1X,' I,SHEAR AREA, L , E ,AREA'/)
00280 NF=5
00290 11 DO 50 I=1,N
00300 NA=-NF
00310 NO=NF
00320 NWMX=1
00330 NFMX=5
00340 CALL FREAD(NA,NA,NO,DATA)
00350 XI=DATA(1)
00360 ABAR=DATA(2)
00370 XL=DATA(3)
00380 E=DATA(4)
00390 AREA=DATA(5)
00400 G=E/2.4
00410 IF(ABAR) 98,15,20
00420 15 BETA=0.0
00430 GO TO 25
00440 20 BETA=(6.0*E*XI)/((XL**2)*ABAR*G)
00450 25 XKMM=(2.0*E*XI*(2.0+BETA))/(XL*(1.0+2.0*BETA))
00460 XKMV=(2.0*E*XI*(1.0-BETA))/(XL*(1.0+2.0*BETA))
00470 K=(I-1)*NTYPE*(NRC+1)+1+INCA
00480 IF(NTYPE-3) 35,30,40
00490 C***********************************************************************
00500 C 3 X 3 ELEMENT ( BENDING PLUS ONE AXIAL)
00510 C***********************************************************************
00520 30 A(K)=AREA*E/XL
00530 K=K+NRC+1
00540 C***********************************************************************
00550 C 2 X 2 ELEMENT ( BENDING ONLY)
00560 C***********************************************************************
00570 35 A(K)=XKMM
00580 K=K+1
00590 A(K)=XKMV
00600 K=K+NRC-1
00610 A(K)=XKMV
00620 K=K+1
00630 A(K)=XKMM
00640 GO TO 50
00650 C***********************************************************************
00660 C 4 X 4 ELEMENT ( BENDING AND TRANSVERSE SHEAR)
00670 C***********************************************************************
00680 40 XKVV=(12.0*E*XI)/((XL**3)*(1.0+2.0*BETA))
00690 XKVM=0.5*XL*XKVV
00700 A(K)=XKVV
00710 L=K
00720 K=K+1
00730 A(K)=-XKVM
00740 K=K+1
00750 A(K)=-XKVV
00760 K=K+1
00770 A(K)=-XKVM
00780 K=K+NRC-NTYPE+1
00790 A(K)=-XKVM
00800 K=K+1
00810 A(K)=XKMM
00820 K=K+1
00830 A(K)=XKVM
00840 K=K+1
00850 A(K)=XKMV
00860 K=K+NRC-NTYPE+1
00870 A(K)=-XKVV
00880 K=K+1
00890 A(K)=XKVM
00900 K=K+1
00910 A(K)=XKVV
00920 K=K+1
00930 A(K)=XKVM
00940 K=K+NRC-NTYPE+1
00950 A(K)=-XKVM
00960 K=K+1
00970 A(K)=XKMV
00980 K=K+1
00990 A(K)=XKVM
01000 K=K+1
01010 A(K)=XKMM
01020 50 CONTINUE
01030 GO TO 99
01040 98 ERROR=8
01050 99 RETURN
01060 END
00010 C***********************************************************************
00020 C SUBROUTINE TO FORM COMPLETE MEMBER STIFFNESS MATRICES
00030 C***********************************************************************
00040 SUBROUTINE FORMKD(NTYPE,N)
00050 COMMON /FF/ NFMX,ITAPE,NWMX
00060 COMMON /YY/ INCA, INC1
00070 COMMON /Z/ A(1),NAME(1)
00080 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00090 DIMENSION DATA(9)
00100 NRC=NTYPE
00110 NS=NRC*NRC
00120 C***********************************************************************
00130 C INITIALIZE MATRIX, READ IN DATA, AND CALCULATE STIFFNESS CONSTANTS
00140 C***********************************************************************
00150 DO 2 I=1,9
00160 DATA(I)=0.0
00170 2 CONTINUE
00180 TYPE 100
00190 100 FORMAT(' ENTER THE FOLLOWING DATA FOR EACH MEMBER'/)
00200 IF(NTYPE.EQ.6) GO TO 6
00210 TYPE 103
00220 NF=6
00230 GO TO 15
00240 6 TYPE 101
00250 101 FORMAT(3X,'I , SHEAR AREA,'$)
00260 TYPE 103
00270 103 FORMAT('+ E , CROSS-AREA, X1 , Y1 , X2 , Y2'//)
00280 NF=8
00290 15 DO 90 I=1,N
00300 NAMA=NAME(NUM-N+I+INC1)
00310 CALL MFIND(NAMA)
00320 DO 18 J=1,NS
00330 18 A(INCA+J)=0.0
00340 NO=NF
00350 NA=-NF
00360 NWMX=1
00370 NFMX=9
00380 CALL FREAD(NA,NA,NO,DATA)
00390 IF(NTYPE.EQ.6) GO TO 20
00400 XI=0
00410 ABAR=0
00420 E=DATA(1)
00430 AREA=DATA(2)
00440 X1=DATA(3)
00450 Y1=DATA(4)
00460 X2=DATA(5)
00470 Y2=DATA(6)
00480 GO TO 25
00490 20 XI=DATA(1)
00500 ABAR=DATA(2)
00510 E=DATA(3)
00520 AREA=DATA(4)
00530 X1=DATA(5)
00540 Y1=DATA(6)
00550 X2=DATA(7)
00560 Y2=DATA(8)
00570 25 G=E/2.4
00580 DX=X2-X1
00590 DY=Y2-Y1
00600 XL=SQRT(DX**2+DY**2)
00610 S=DY/XL
00620 C=DX/XL
00630 K=INCA+1
00640 XKAA=AREA*E/XL
00650 TYPE 105, XL
00660 105 FORMAT(5X,'LENGTH =',1PE14.7//)
00670 IF(NTYPE-4) 98,10,50
00680 C***********************************************************************
00690 C 4 X 4 ELEMENT (AXIAL ONLY)
00700 C***********************************************************************
00710 10 A(K)=XKAA*C**2
00720 L=K
00730 K=K+1
00740 A(K)=XKAA*C*S
00750 LP1=L+1
00760 K=K+1
00770 A(K)=-A(L)
00780 K=K+1
00790 A(K)=-A(LP1)
00800 K=K+1
00810 A(K)=A(LP1)
00820 K=K+1
00830 A(K)=XKAA*S**2
00840 LPN1=LP1+NRC
00850 K=K+1
00860 A(K)=-A(LP1)
00870 K=K+1
00880 A(K)=-A(LPN1)
00890 GO TO 75
00900 C***********************************************************************
00910 C 6 X 6 ELEMENT (BENDING, TRANSVERSE SHEAR, AND AXIAL)
00920 C***********************************************************************
00930 50 IF(ABAR) 98,55,60
00940 55 BETA=0.0
00950 GO TO 65
00960 60 BETA=(6.0*E*XI)/((XL**2)*ABAR*G)
00970 65 XKVV=(12.0*E*XI)/((XL**3)*(1.0+2.0*BETA))
00980 XKVM=0.5*XL*XKVV
00990 XKMM=(2.0*E*XI*(2.0+BETA))/(XL*(1.0+2.0*BETA))
01000 XKMV=(2.0*E*XI*(1.0-BETA))/(XL*(1.0+2.0*BETA))
01010 A(K)=XKVV*S**2+XKAA*C**2
01020 L=K
01030 K=K+1
01040 A(K)=(-XKVV+XKAA)*C*S
01050 K=K+1
01060 LP1=L+1
01070 A(K)=XKVM*S
01080 K=K+1
01090 A(K)=-A(L)
01100 K=K+1
01110 A(K)=-A(LP1)
01120 K=K+1
01130 A(K)=XKVM*S
01140 K=K+1
01150 A(K)=A(LP1)
01160 K=K+1
01170 A(K)=XKVV*C**2+XKAA*S**2
01180 LPN1=LP1+NRC
01190 K=K+1
01200 A(K)=-XKVM*C
01210 K=K+1
01220 A(K)=-A(LP1)
01230 K=K+1
01240 A(K)=-A(LPN1)
01250 K=K+1
01260 A(K)=-XKVM*C
01270 K=K+1
01280 A(K)=XKVM*S
01290 K=K+1
01300 A(K)=-XKVM*C
01310 K=K+1
01320 A(K)=XKMM
01330 K=K+1
01340 A(K)=-XKVM*S
01350 K=K+1
01360 A(K)=XKVM*C
01370 K=K+1
01380 A(K)=XKMV
01390 C***********************************************************************
01400 C DUPLICATE THE ROWS CORRESPONDING TO THE R1 AND R2 FREEDOMS
01410 C***********************************************************************
01420 75 K=K+1
01430 DO 80 J=1,2
01440 K=K+(J-1)*NRC
01450 KPN=K+NTYPE-1
01460 LL=L+(J-1)*NRC
01470 DO 80 JJ=K,KPN
01480 A(JJ)=-A(LL)
01490 80 LL=LL+1
01500 IF(NTYPE-4) 98,90,85
01510 85 K=K+(J-1)*NRC
01520 A(K)=XKVM*S
01530 K=K+1
01540 A(K)=-XKVM*C
01550 K=K+1
01560 A(K)=XKMV
01570 K=K+1
01580 A(K)=-XKVM*S
01590 K=K+1
01600 A(K)=XKVM*C
01610 K=K+1
01620 A(K)=XKMM
01630 90 CONTINUE
01640 GO TO 99
01650 98 ERROR=8
01660 99 RETURN
01670 END
00010 C***********************************************************************
00020 C SUBROUTINE TO PERFORM DIRECT STIFFNESS MERGE
00030 C***********************************************************************
00040 SUBROUTINE MERGE(NBIG,NRA,NCA,K)
00050 INTEGER ERROR,FREZ,FRET
00060 COMMON /FF/ NFMX,ITAPE,NWMX
00070 COMMON /XX/ ERROR
00080 COMMON /YY/ INCA
00090 COMMON /Z/ A(1)
00100 COMMON /ZZ/ NUM,NNUM,NX,NRB,NCB
00110 DIMENSION NROW(1),NCOL(1),W1(1)
00120 C************************************************************
00130 C GET FREE STORAGE
00140 C*******************************************************
00150 NS=3*NRB
00160 NP=FREZ(NS)
00170 N1=NP-LOC(NROW)
00180 N2=NP+NRB-LOC(NCOL)
00190 N3=NP+2*NRB-LOC(W1)
00200 C***********************************************************************
00210 C CHECK FOR SYMMETRICAL MERGE
00220 C***********************************************************************
00230 IF(IABS(K)-1) 98,1,5
00240 1 IF(NRB-NCB)91,2,91
00250 2 TYPE 100
00260 100 FORMAT(2X,' SYMMETRIC MERGE - ENTER ROW-COL COORDINATES'/)
00270 NF=NRB
00280 NA=-NF
00290 NFMX=NRB
00300 NWMX=1
00310 CALL FREAD(NA,NA,NF,W1(N3+1))
00320 DO 3 I=1,NRB
00330 NROW(I+N1)=W1(I+N3)
00340 NCOL(I+N2)=W1(I+N3)
00350 3 CONTINUE
00360 GO TO 6
00370 5 TYPE 101
00380 101 FORMAT(2X,' ENTER ROW COORDINATES'/)
00390 NF=NRB
00400 NA=-NF
00410 NFMX=NRB
00420 NWMX=1
00430 CALL FREAD(NA,NA,NF,W1(N3+1))
00440 DO 4 I=1,NRB
00450 NROW(I+N1)=W1(I+N3)
00460 4 CONTINUE
00470 TYPE 102
00480 102 FORMAT(2X,' ENTER COLUMN COORDINATES'/)
00490 NF=NCB
00500 NA=-NF
00510 NFMX=NRB
00520 NWMX=1
00530 CALL FREAD(NA,NA,NF,W1(N3+1))
00540 DO 50 I=1,NCB
00550 NCOL(I+N2)=W1(I+N3)
00560 50 CONTINUE
00570 6 DO 8 I=1,NRB
00580 IF(NRA-IABS(NROW(I+N1))) 91,8,8
00590 8 CONTINUE
00600 DO 10 J=1,NCB
00610 IF(NCA-IABS(NCOL(J+N2))) 91,10,10
00620 10 CONTINUE
00630 DO 21 I=1,NRB
00640 IF(NROW(I+N1)) 11,21,12
00650 11 TAGI=-1.0
00660 NROW(I+N1)=IABS(NROW(I+N1))
00670 GO TO 13
00680 12 TAGI=1.0
00690 13 DO 20 J=1,NCB
00700 IF(NCOL(J+N2)) 14,20,15
00710 14 IF(TAGI) 15,98,16
00720 15 TAGJ=1.0
00730 NCOL(J+N2) =IABS(NCOL(J+N2))
00740 GO TO 17
00750 16 TAGJ=-1.0
00760 NCOL(J+N2) =IABS(NCOL(J+N2))
00770 17 NNA=NCA*(NROW(I+N1)-1)+NCOL(J+N2)
00780 NNB=NCB*(I-1)+J
00790 IF(K) 18,98,19
00800 18 A(NBIG+NNA)=TAGI*TAGJ*A(NNB+INCA)
00810 GO TO 20
00820 19 A(NNA+NBIG)=A(NNA+NBIG)+TAGI*TAGJ*A(INCA+NNB)
00830 20 CONTINUE
00840 21 CONTINUE
00850 GO TO 99
00860 91 ERROR=1
00870 GO TO 99
00880 98 ERROR=8
00890 99 CALL FRET(NS,NP)
00900 RETURN
00910 END
00010 C***********************************************************************
00020 C SUBROUTINE TO STORE OR REMOVE A SUBMATRIX
00030 C***********************************************************************
00040 SUBROUTINE STOSM(NB,NRB,NCB,JJ,KK,NTAG)
00050 INTEGER ERROR
00060 COMMON /XX/ ERROR
00070 COMMON /YY/ INCA
00080 COMMON /Z/ A(1)
00090 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00100 IF (NR+1-JJ-NRB) 60,50,50
00110 50 IF (NC+1-KK-NCB) 60,70,70
00120 60 ERROR=1
00130 GO TO 900
00140 70 DO 100 I=1,NRB
00150 K=(KK-1)+(JJ+I-2)*NC
00160 L=(I-1)*NCB
00170 DO 100 J=1,NCB
00180 N=K+J
00190 M=L+J
00200 IF (NTAG) 75,90,80
00210 75 A(INCA+N)=A(INCA+N)+A(NB+M)
00220 GO TO 100
00230 80 A(INCA+N)=A(NB+M)
00240 GO TO 100
00250 90 A(NB+M)=A(INCA+N)
00260 100 CONTINUE
00270 900 RETURN
00280 END
00010 C***********************************************************************
00020 C SUBROUTINE TO DELETE ROWS AND/OR COLUMNS FROM A MATRIX
00030 C***********************************************************************
00040 SUBROUTINE DELRC(NAMA,NRD,NCD,NAMB,K)
00050 INTEGER ALPHA(1),ERROR,FREZ,FRET
00060 COMMON /FF/ NFMX,ITAPE,NWMX
00070 DIMENSION MROW(1),MCOL(1),W1(1)
00080 COMMON /XX/ ERROR
00090 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00100 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1),NME(4),NSIZE
00110 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00120 C***********************************************************************
00130 C CHECK FOR SYMMETRICAL DELETION
00140 C***********************************************************************
00150 NMX=MAX0(NRD,NCD)
00160 NW1=FREZ(NMX)
00170 IW1=NW1-LOC(W1)
00180 NMR=FREZ(NRD+NCD)
00190 IMR=NMR-LOC(MROW)
00200 IMC=NMR+NRD-LOC(MCOL)
00210 IF(K) 98,5,1
00220 1 IF(NRD-NCD) 91,2,91
00230 2 IF(NRD) 98,98,3
00240 3 NF=NRD
00250 NA=-NF
00260 TYPE 300
00270 300 FORMAT(' ENTER ROWS-COLUMNS TO BE DELETED'/)
00280 NFMX=NMX
00290 NWMX=1
00300 CALL FREAD(NA,ALPHA,NF,W1(IW1+1))
00310 DO 4 I=1,NRD
00320 MROW(I+IMR)=W1(I+IW1)
00330 MCOL(I+IMC)=W1(I+IW1)
00340 4 CONTINUE
00350 GO TO 10
00360 5 IF(NRD)98,7,6
00370 6 NF=NRD
00380 NA=-NF
00390 TYPE 600
00400 600 FORMAT(' ENTER ROWS TO BE DELETED'/)
00410 NFMX=NMX
00420 NWMX=1
00430 CALL FREAD(NA,ALPHA,NF,W1(IW1+1))
00440 DO 70 I=1,NRD
00450 MROW(I+IMR)=W1(I+IW1)
00460 70 CONTINUE
00470 7 IF(NCD) 98,10,8
00480 8 NF=NCD
00490 NA=-NF
00500 TYPE 800
00510 800 FORMAT(' ENTER COLUMNS TO BE DELETED'/)
00520 NFMX=NMX
00530 NWMX=1
00540 CALL FREAD(NA,ALPHA,NF,W1(IW1+1))
00550 DO 90 I=1,NCD
00560 MCOL(I+IMC)=W1(I+IW1)
00570 90 CONTINUE
00580 10 CALL MFIND(NAMA)
00590 IF(ERROR.GT.0) GO TO 99
00600 NA=INCA+1
00610 NRA=NR
00620 NCA=NC
00630 IF(NRD) 98,17,15
00640 15 DO 16 I=1,NRD
00650 IF(NRA-MROW(I+IMR)) 91,16,16
00660 16 CONTINUE
00670 17 IF(NCD) 98,20,18
00680 18 DO 19 J=1,NCD
00690 IF(NCA-MCOL(J+IMC)) 91,19,19
00700 19 CONTINUE
00710 20 NR=NRA-NRD
00720 NC=NCA-NCD
00730 CALL MLIST(NAMB)
00740 IF(ERROR.GT.0) GO TO 99
00750 NB=INCA+1
00760 DO 30 II=1,NRA
00770 DO 28 JJ=1,NCA
00780 IF(NRD) 98,23,21
00790 21 DO 22 I=1,NRD
00800 IF(II-MROW(I+IMR)) 22,29,22
00810 22 CONTINUE
00820 23 IF(NCD) 98,26,24
00830 24 DO 25 J=1,NCD
00840 IF(JJ-MCOL(J+IMC)) 25,28,25
00850 25 CONTINUE
00860 26 A(NB)=A(NA)
00870 NB=NB+1
00880 28 NA=NA+1
00890 GO TO 30
00900 29 NA=NA+NCA
00910 30 CONTINUE
00920 GO TO 99
00930 91 ERROR=1
00940 GO TO 99
00950 98 ERROR=8
00960 99 CALL FRET(NRD+NCD,NMR)
00970 CALL FRET(NMX,NW1)
00980 RETURN
00990 END
00010 C***********************************************************************
00020 C SUBROUTINE TO PERFORM CHOLESKI RIGHT DECOMPOSTITION
00030 C***********************************************************************
00040 SUBROUTINE DECOM(NB)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ A(1)
00070 COMMON /ZZ/ NUM,NNUM,NX,NR,NC
00080 NS=NR*NR
00090 DO 1 I=1,NS
00100 1 A(I+NB)=0.0
00110 NU=1
00120 SRA=A(INCA+NU)
00130 IF(SRA.LT.0) GO TO 90
00140 A(NU+NB)=SQRT(SRA)
00150 DO 2 K=2,NR
00160 2 A(K+NB)=A(K+INCA)/A(NU+NB)
00170 DO 10 I=2,NR
00180 IP1=I+1
00190 IM1=I-1
00200 5 NU=(I-1)*NR+I
00210 A(NU+NB)=A(NU+INCA)
00220 DO 6 II=1,IM1
00230 NP=NU-II*NR
00240 6 A(NU+NB)=A(NU+NB)-(A(NP+NB)**2)
00250 SRA=A(NU+NB)
00260 IF(SRA.LT.0) GO TO 90
00270 A(NU+NB)=SQRT(SRA)
00280 IF(I-NR) 7,10,10
00290 7 NUSTR=NU
00300 DO 9 J=IP1,NR
00310 NU=NU+1
00320 A(NU+NB)=A(NU+INCA)
00330 DO 8 JJ=1,IM1
00340 NP=NU-JJ*NR
00350 NNP=NUSTR-JJ*NR
00360 8 A(NU+NB)=A(NU+NB)-A(NP+NB)*A(NNP+NB)
00370 9 A(NU+NB)=A(NU+NB)/A(NUSTR+NB)
00380 10 CONTINUE
00390 99 RETURN
00400 90 TYPE 100
00410 100 FORMAT(2X,' MATRIX IS NOT POSITIVE DEFINITE'/)
00420 END
00010 C***********************************************************************
00020 C SUBROUTINE TO PERFORM FORWARD OR BACK SUBSTITUTION AND INVERSION
00030 C***********************************************************************
00040 SUBROUTINE CHSOLV(NA,IFLAG,SINGLR)
00050 INTEGER ERROR,FREZ,FRET
00060 COMMON /XX/ ERROR
00070 COMMON /YY/ INCA
00080 COMMON /Z/ A(1)
00090 COMMON /ZZ/ NUM,NNUM,NX,N,M
00100 DIMENSION C(1)
00110 C*************************************************
00120 C GET FREE STORAGE FOR ARRAY C
00130 C***************************************************
00140 NP=FREZ(N)
00150 LC=NP-LOC(C)
00160 NM1=N-1
00170 DO 10 I=1,NM1
00180 IP1=I+1
00190 IF(IFLAG) 98,1,2
00200 1 II=I
00210 GO TO 3
00220 2 II=N-I+1
00230 3 DO 10 J=IP1,N
00240 IF(IFLAG) 98,6,7
00250 6 JJ=(II-1)*N+J
00260 GO TO 8
00270 7 JJ=II*N-J+1
00280 8 IF(ABS(A(JJ+NA))-SINGLR) 10,10,97
00290 10 CONTINUE
00300 DO 70 I=1,N
00310 IM1=I-1
00320 IF(IFLAG) 98,12,13
00330 12 II=I
00340 GO TO 14
00350 13 II=N-I+1
00360 14 IF(IM1) 98,15,30
00370 15 K=(II-1)*N+II
00380 IF(SINGLR) 98,16,17
00390 16 SINGLR=ABS(A(K+NA)*1.0E-06)
00400 17 TYPE 101, SINGLR
00410 101 FORMAT(' SCALAR =',1PE14.7/)
00420 IF(M) 98,25,18
00430 18 JJ=(II-1)*M
00440 DO 20 J=1,M
00450 JJ=JJ+1
00460 20 A(JJ+INCA)=A(JJ+INCA)/A(K+NA)
00470 25 A(K+NA)=1.0/A(K+NA)
00480 GO TO 70
00490 30 DO 35 J=1,I
00500 IF(IFLAG) 98,32,33
00510 32 JJ=(II-1)*N+J
00520 GO TO 35
00530 33 JJ=II*N-J+1
00540 35 C(J+LC)=A(JJ+NA)
00550 IF(SINGLR-ABS(C(I+LC))) 38,96,96
00560 38 IF(M) 98,55,40
00570 40 DO 50 J=1,M
00580 JJ=(II-1)*M+J
00590 IF(IFLAG) 98,42,43
00600 42 JL=J
00610 GO TO 45
00620 43 JL=(N-1)*M+J
00630 45 DO 49 L=1,IM1
00640 A(JJ+INCA)=A(JJ+INCA)-C(L+LC)*A(JL+INCA)
00650 IF(IFLAG) 98,46,47
00660 46 JL=JL+M
00670 GO TO 49
00680 47 JL=JL-M
00690 49 CONTINUE
00700 50 A(JJ+INCA)=A(JJ+INCA)/C(I+LC)
00710 55 K=(II-1)*N+II
00720 A(K+NA)=1.0/C(I+LC)
00730 DO 65 J=1,IM1
00740 IF(IFLAG) 98,57,58
00750 57 KK=(II-1)*N+J
00760 GO TO 59
00770 58 KK=II*N-J+1
00780 59 A(KK+NA)=0.0
00790 DO 64 L=J,IM1
00800 IF(IFLAG) 98,61,62
00810 61 LL=(L-1)*N+J
00820 GO TO 64
00830 62 LL=N*N-N*(L-1)-J+1
00840 64 A(KK+NA)=A(KK+NA)-C(L+LC)*A(LL+NA)
00850 65 A(KK+NA)=A(KK+NA)/C(I+LC)
00860 70 CONTINUE
00870 GO TO 99
00880 96 TYPE 102, II,II,C(I)
00890 GO TO 99
00900 97 ERROR=1
00910 GO TO 99
00920 98 ERROR=8
00930 99 CALL FRET(N,NP)
00940 RETURN
00950 102 FORMAT(2X,' (',I3,',',I3,') =',1PE14.7,3X,'WILL NOT PRODUCE '
00960 1 ,/' AN ACCURATE SINGLE PRECISION INVERSION'/)
00970 END
00010 C***********************************************************************
00020 SUBROUTINE EIGEN(R,V,E,LP,M,NM)
00030 COMMON /YY/ INCA
00040 COMMON /Z/ A(1)
00050 INTEGER FREZ,FRET
00060 DIMENSION R(1),E(1),V(1),AA(1),B(1),W1(1),W2(1)
00070 NAA=FREZ(LP)
00080 NB=FREZ(LP)
00090 NW1=FREZ(LP)
00100 NW2=FREZ(LP)
00110 INA=NAA-LOC(AA)+1
00120 INB=NB-LOC(B)+1
00130 IW1=NW1-LOC(W1)+1
00140 IW2=NW2-LOC(W2)+1
00150 IF(LP-1) 20,14,1
00160 1 DO 2 I=1,LP
00170 2 A(I+INCA)=1./(SQRT(A(I+INCA)))
00180 DO 3 I=1,LP
00190 DO 3 J=1,LP
00200 JJ=I+(J-1)*NM
00210 3 R(JJ)=A(I+INCA)*R(JJ)*A(J+INCA)
00220 C***********************************************************************
00230 C TRI-DIAGONALIZE MATRIX
00240 C***********************************************************************
00250 4 CALL TRIDI(LP,NM,R,AA(INA),B(INB),W1(IW1),W2(IW2))
00260 IF(ERROR) 20,40,20
00270 C***********************************************************************
00280 C FIND EIGEN VALUES
00290 C***********************************************************************
00300 40 CALL EIGVAL(LP,NM,M,E,AA(INA),B(INB),W1(IW1),W2(IW2))
00310 IF (M) 5,14,5
00320 C***********************************************************************
00330 C FIND EIGEN VECTORS
00340 C***********************************************************************
00350 5 K=IABS(M)
00360 J=1
00370 DO 7 I=1,K
00380 CALL EIGVEC(LP,NM,R,AA(INA),B(INB),E(I),V(J),W1(IW1),W2(IW2))
00390 7 J=J+NM
00400 C***********************************************************************
00410 C MODIFY EIGEN VECTORS
00420 C***********************************************************************
00430 9 DO 15 I=1,LP
00440 DO 15 J=1,K
00450 JJ=I+(J-1)*NM
00460 15 V(JJ)=A(I+INCA)*V(JJ)
00470 GO TO 20
00480 C***********************************************************************
00490 C SOLUTION FOR ORDER 1
00500 C***********************************************************************
00510 14 E(1)=R(1)/E(1)
00520 V(1)=1.0
00530 20 CALL FRET(LP,NAA)
00540 CALL FRET(LP,NB)
00550 CALL FRET(LP,NW1)
00560 CALL FRET(LP,NW2)
00570 RETURN
00580 END
00010 C***********************************************************************
00020 C TRI-DIAGONALIZATION SUBROUTINE DWM 1517-UB
00030 C***********************************************************************
00040 C RETURN ORIGNAL R IN UPPER TRIANGULAR HALF INCLUDING DIAGONAL
00050 C RETURN MODIFIED W MATRICES IN LOWER HALF OF R MATRIX
00060 C RETURN NEW DIAGONAL IN A
00070 C RETURN NEW FIRST OFF DIAGONAL IN B
00080 C***********************************************************************
00090 SUBROUTINE TRIDI(LP,NM,R,A,B,W,Q)
00100 INTEGER ERROR
00110 COMMON /XX/ ERROR
00120 DIMENSION R(1),A(1),B(1),Q(1),W(1)
00130 LP1=LP-1
00140 LP2=LP1*NM+LP
00150 LPP=LP2-NM
00160 NM1=NM+1
00170 C***********************************************************************
00180 C STORE ORIGINAL DIAGONAL
00190 C***********************************************************************
00200 L=0
00210 DO 10 I=1,LP2,NM1
00220 L=L+1
00230 10 A(L)=R(I)
00240 C***********************************************************************
00250 B(1)=0.
00260 C***********************************************************************
00270 IF(LP-2)99,65,11
00280 11 SUM=0.
00290 DO 12 I=1,LP-2
00300 NR=(I-1)*LP
00310 DO 12 J=1,LP
00320 IF((J-I).LT.2) GO TO 12
00330 SUM=SUM+ABS(R(NR+J))
00340 12 CONTINUE
00350 IF(SUM.GT.0) GO TO 15
00360 SMALL=R(1)*1.E-8
00370 IF(SMALL.EQ.0) SMALL=1.E-8
00380 R(LP)=SMALL
00390 R(LP*LP-LP+1)=SMALL
00400 15 KK=0
00410 DO 50 K=2,LP1
00420 KL=KK+K
00430 KU=KK+LP
00440 KJ=K+1
00450 C***********************************************************************
00460 C MAKE SURE THERE IS AT LEAST ONE NON-ZERO ELEMENT IN THE FIRST ROW
00470 C***********************************************************************
00480 IF(K-2) 99,16,19
00490 16 DO 17 J=2,KU
00500 IF(R(J)-0.0) 19,17,19
00510 17 CONTINUE
00520 18 TYPE 101
00530 ERROR=-1
00540 GO TO 99
00550 C***********************************************************************
00560 C CALCULATE AND STORE MODIFIED COLUMN MATRIX W
00570 C***********************************************************************
00580 19 SUM=0.0
00590 DO 20 J=KL,KU
00600 20 SUM=SUM+R(J)**2
00610 S=SQRT(SUM)
00620 B(K)=SIGN(S,-R(KL))
00630 S=1./S
00640 W(K)=SQRT(ABS(R(KL))*S+1.)
00650 X=SIGN(S/W(K),R(KL))
00660 R(KL)=W(K)
00670 DO 30 I=KJ,LP
00680 JJ=I+KK
00690 W(I)=X*R(JJ)
00700 30 R(JJ)=W(I)
00710 C***********************************************************************
00720 C CALCULATE NEW R MATRIX WITH ROW K-1 NOW HAVING ZEROS OFF 2ND DIAGONAL
00730 C***********************************************************************
00740 DO 35 J=K,LP
00750 JJ=J+1
00760 Q(J)=0.0
00770 L=KK+J
00780 DO 33 I=K,J
00790 L=L+NM
00800 33 Q(J)=Q(J)+R(L)*W(I)
00810 IF(JJ-LP)34,34,36
00820 34 DO 35 I=JJ,LP
00830 L=L+1
00840 35 Q(J)=Q(J)+R(L)*W(I)
00850 36 X=0.0
00860 DO 40 J=K,LP
00870 40 X=X+W(J)*Q(J)
00880 X=.5*X
00890 DO 45 I=K,LP
00900 45 Q(I)=X*W(I)-Q(I)
00910 LL=KK
00920 KK=KK+NM
00930 DO 50 I=K,LP
00940 LL=LL+NM
00950 DO 50 J=I,LP
00960 L=LL+J
00970 50 R(L)=R(L)+Q(I)*W(J)+Q(J)*W(I)
00980 C***********************************************************************
00990 C SORT OUTPUT
01000 C***********************************************************************
01010 L=1
01020 DO 60 I=1,LP
01030 X=A(I)
01040 A(I)=R(L)
01050 R(L)=X
01060 60 L=L+NM1
01070 65 B(LP)=R(LPP)
01080 99 RETURN
01090 101 FORMAT(48H4NO NON-ZERO OFF-DIAGONAL ELEMENTS IN FIRST ROW./
01100 1' DIAGONAL ELEMENT IS AN EIGENVALUE--'/' NEED ONLY PERFORM EIGEN '
01110 2' OPERATION ON REMAINING SUBMATRIX'/)
01120 END
00010 C***********************************************************************
00020 C SUBROUTINE TO GENERATE EIGENVECTORS
00030 C***********************************************************************
00040 SUBROUTINE EIGVEC (LP,NM,R,A,B,E,V,P,Q)
00050 DIMENSION R(1),A(1),B(1),V(1),P(1),Q(1)
00060 C***********************************************************************
00070 C SET UP SIMULTANEOUS EQUATIONS FOR EIGEN VECTOR WITH EIGEN VALUE E
00080 C***********************************************************************
00090 X=A(1)-E
00100 Y=B(2)
00110 LP1=LP-1
00120 DO 10 I=1,LP1
00130 IF(ABS(X)-ABS(B(I+1))) 4,6,8
00140 4 P(I)=B(I+1)
00150 Q(I)=A(I+1)-E
00160 V(I)=B(I+2)
00170 Z=-X/P(I)
00180 X=Z*Q(I)+Y
00190 IF(LP1-I)5,10,5
00200 5 Y=Z*V(I)
00210 GO TO 10
00220 6 IF(X)8,7,8
00230 7 X=1.0E-10
00240 8 P(I)=X
00250 Q(I)=Y
00260 V(I)=0.
00270 X=A(I+1)-(B(I+1)/X*Y+E)
00280 Y=B(I+2)
00290 10 CONTINUE
00300 C***********************************************************************
00310 C SOLVE SIMULTANEOUS EQUATIONS FOR EIGEN VECTOR OF TRI-DIAGONAL MATRIX
00320 C***********************************************************************
00330 20 IF(X) 21,28,21
00340 21 V(LP)=1./X
00350 22 I=LP1
00360 V(I)=(1.-Q(I)*V(LP))/P(I)
00370 X=V(LP)**2+V(I)**2
00380 25 I=I-1
00390 IF(I) 26,30,26
00400 26 V(I)=(1.-(Q(I)*V(I+1)+V(I)*V(I+2)))/P(I)
00410 X=X+V(I)**2
00420 GO TO 25
00430 28 V(LP)=1.0E10
00440 GO TO 22
00450 30 X=SQRT(X)
00460 DO 31 I=1,LP
00470 31 V(I)=V(I)/X
00480 C***********************************************************************
00490 C TRANSFORM EIGEN VECTOR TO SOLUTION OF ORIGINAL MATRIX
00500 C***********************************************************************
00510 J=LP1*NM-NM
00520 K=LP
00530 GO TO 42
00540 32 K=K-1
00550 J=J-NM
00560 Y=0.0
00570 DO 35 I=K,LP
00580 L=J+I
00590 35 Y=Y+V(I)*R(L)
00600 DO 40 I=K,LP
00610 L=J+I
00620 40 V(I)=V(I)-Y*R(L)
00630 42 IF(J)32,44,32
00640 44 RETURN
00650 END
00010 C***********************************************************************
00020 C EIGEN VALUE SUBROUTINE FOR TRI-DIAGONAL MATRICES DWM 1517-UB
00030 C***********************************************************************
00040 SUBROUTINE EIGVAL(LP,NM,M,E,A,B,F,W)
00050 DIMENSION E(1),A(1),B(1),F(1),W(1)
00060 EQUIVALENCE (S1,IS1),(S2,IS2)
00070 C***********************************************************************
00080 C FIND UPPER AND LOWER BOUNDS AND NORMALIZE INPUT
00090 C***********************************************************************
00100 BD=ABS(A(1))
00110 DO 5 I=2,LP
00120 5 BD=AMAX1(BD,ABS(A(I))+B(I)**2)
00130 BD=BD+1.
00140 DO 6 I=1,LP
00150 A(I)=A(I)/BD
00160 B(I)=B(I)/BD
00170 W(I)=1.
00180 6 E(I)=-1.
00190 DO 50 K=1,LP
00200 8 IF((W(K)-E(K))/AMAX1(ABS(W(K)),ABS(E(K)),1.0E-9)-1.0E-7) 50,50,10
00210 10 X=(W(K)+E(K))*.5
00220 C***********************************************************************
00230 C FIND NUMBER OF EIGEN VALUES ,N, GREATER THAN OR EQUAL TO X
00240 C***********************************************************************
00250 IS2=1
00260 F(1)=A(1)-X
00270 IF(F(1))102,104,104
00280 102 IS1=-1
00290 N=0
00300 GO TO 105
00310 104 IS1=1
00320 N=1
00330 105 DO 120 I=2,LP
00340 IF(B(I))106,113,106
00350 106 IF(B(I-1))107,114,107
00360 107 IF(ABS(F(I-1))+ABS(F(I-2))-1.0E-15) 110,112,112
00370 110 F(I-1)=F(I-1)*1.0E15
00380 F(I-2)=F(I-2)*1.0E15
00390 112 F(I)=(A(I)-X)*F(I-1)-B(I)**2*F(I-2)
00400 GO TO 115
00410 113 F(I)=(A(I)-X)*SIGN(1.,S1)
00420 GO TO 115
00430 114 F(I)=(A(I)-X)*F(I-1)-SIGN(B(I)**2,S2)
00440 115 S2=S1
00450 IF(F(I))116,117,116
00460 116 S1=SIGN(S1,F(I))
00470 IF(IS2+IS1)117,120,117
00480 117 N=N+1
00490 120 CONTINUE
00500 C***********************************************************************
00510 C TRAP EIGEN VALUES IN SMALLER AND SMALLER BOUNDS
00520 C***********************************************************************
00530 N=LP-N
00540 IF(N-K)20,12,12
00550 12 DO 15 J=K,N
00560 15 W(J)=X
00570 20 N=N+1
00580 IF(LP-N)8,24,24
00590 24 DO 26 J=N,LP
00600 IF(X-E(J))8,8,26
00610 26 E(J)=X
00620 GO TO 8
00630 50 CONTINUE
00640 C***********************************************************************
00650 C RESTORE INPUT AND ORDER EIGEN VALUES
00660 C***********************************************************************
00670 DO 60 I=1,LP
00680 A(I)=A(I)*BD
00690 B(I)=B(I)*BD
00700 60 F(I)=(W(I)+E(I))*BD*.5
00710 J=LP
00720 K=1
00730 DO 70 I=1,LP
00740 IF(ABS(F(K))-ABS(F(J))) 63,63,65
00750 63 E(I)=F(J)
00760 J=J-1
00770 GO TO 70
00780 65 E(I)=F(K)
00790 K=K+1
00800 70 CONTINUE
00810 IF(ISIGN(1,M)) 75,80,80
00820 75 DO 77 I=1,LP
00830 77 F(I)=E(I)
00840 J=LP
00850 DO 79 I=1,LP
00860 E(I)=F(J)
00870 79 J=J-1
00880 80 CONTINUE
00890 RETURN
00900 END
00010 C***********************************************************************
00020 C SUBROUTINE TO GENERATE FUNCTION FROM DISCRETE DATA POINTS
00030 C***********************************************************************
00040 SUBROUTINE FUNGN(NB,N,M,SCAL1)
00050 INTEGER ERROR
00060 COMMON /XX/ ERROR
00070 COMMON /YY/ INCA
00080 COMMON /Z/ A(1)
00090 K=1+NB
00100 L=1+INCA
00110 TIM =A(L)
00120 10 IF(M+INCA-L-3) 30,20,20
00130 20 DT=A(L+2)-A(L)
00140 DP=A(L+3)-A(L+1)
00150 L=L+2
00160 IF (DT) 30,10,40
00170 30 ERROR=1
00180 GO TO 70
00190 40 SLOPE=DP/DT
00200 50 IF (A(L)-TIM) 10,10,60
00210 60 A(K)=A(L-1)+(TIM -A(L-2))*SLOPE
00220 TIM =TIM +SCAL1
00230 K=K+1
00240 IF(N+NB-K) 70,65,50
00250 65 IF (A(L)-TIM) 10,66,60
00260 66 A(K)=A(L+1)
00270 70 RETURN
00280 END
00010 C***********************************************************************
00020 C RESPONSE PROGRAM
00030 C***********************************************************************
00040 SUBROUTINE RESPON(W,P,X,NM,NT,L,DAMP,DT,NTAG,NTYPE)
00050 COMMON /YY/ INCA
00060 COMMON /Z/ XM(1)
00070 DIMENSION W(1),P(1),X(1)
00080 NTYPE=NTYPE+1
00090 K=1
00100 C1=DT/2.
00110 C2=C1*DT/3.
00120 C3=C2*2.
00130 DO 160 N=1,NM
00140 NOUT=NT+1
00150 C4=W(N)**2
00160 C5=2.*DAMP*W(N)
00170 F=1.+C1*C5+C2*C4
00180 DISP=0.0
00190 VEL=0.0
00200 IL=L*(N-1)*NTAG+1
00210 ACEL=P(IL)*XM(N+INCA)
00220 IL=IL+1
00230 IH=IL+L-2
00240 DO 160 I=IL,IH
00250 A=VEL+C1*ACEL
00260 B=DISP+DT*VEL+C3*ACEL
00270 ACEL=(P(I)*XM(N+INCA)-C5*A-C4*B)/F
00280 VEL=A+C1*ACEL
00290 DISP=B+C2*ACEL
00300 IF(NOUT-I) 160,150,160
00310 150 GO TO (151,152,153),NTYPE
00320 151 X(K)=DISP
00330 GO TO 154
00340 152 X(K)=VEL
00350 GO TO 154
00360 153 X(K)=ACEL
00370 154 K=K+1
00380 NOUT=NOUT+NT
00390 160 CONTINUE
00400 RETURN
00410 END
00010 SUBROUTINE MSTOR
00020 INTEGER FREC,FRET,ERROR
00030 COMMON /XX/ ERROR
00040 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00050 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1),NME(4),NSIZE
00060 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00070 NNUM=NNUM+10
00080 DO 10 I=1,4
00090 IADDR=NME(I)
00100 ITMP=FREC(NNUM,IADDR)
00110 CALL FRET(NNUM-10,IADDR)
00120 10 NME(I)=ITMP
00130 INC1=NME(1)-LOC(NAME)
00140 INC2=NME(2)-LOC(NCOL)
00150 INC3=NME(3)-LOC(NROW)
00160 INC4=NME(4)-LOC(NN)
00170 RETURN
00180 END
00010 C * * * * *GETS MATRICES FROM DISK FILE * * * * *
00020 SUBROUTINE MGET(NAM)
00030 INTEGER FREZ,FRET,DISKIO,ERROR
00040 COMMON /XX/ ERROR
00050 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00060 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1),NME(4),NSIZE
00070 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00080 DIMENSION FBUF(1),NFILN(2)
00090 DATA NFILN(2) /5H .MAT/
00100 C * * * * * REQUEST FREE STORAGE * * * *
00110 SLASH='//'
00120 IPT=FREZ(128)
00130 INC=IPT-LOC(FBUF)
00140 NFILN(1)=NAM
00150 NMB=0
00160 NREC=1
00170 ISET=0
00180 NP=1
00190 IERR=DISKIO(0,NFILN,NREC,FBUF(1+INC))
00200 IF (IERR.NE.0) GO TO 50
00210 10 IF ((NP+2).LT.129) GO TO 20
00220 11 NREC=NREC+1
00230 IERR=DISKIO(0,NFILN,NREC,FBUF(INC+1))
00240 IF (IERR.NE.0) GO TO 50
00250 NP=1
00260 IF (ISET.EQ.1) GO TO 45
00270 C * * * * * LIST MATRIX * * * * *
00280 20 DO 21 J=1,5
00290 CALL GET(FBUF(NP+INC),J,IC)
00300 CALL PUT(NAMA,J,IC)
00310 21 CONTINUE
00320 NR=FBUF(NP+1+INC)
00330 NC=FBUF(NP+2+INC)
00340 NP=NP+2
00350 NS=NR*NC
00360 CALL MLIST(NAMA)
00370 IF (ERROR.EQ.0) GO TO 25
00380 IF (ERROR.NE.4) GO TO 60
00390 ERROR=0
00400 NP=NP-2
00410 CALL DELETE(NAMA)
00420 TYPE 22,NAMA
00430 22 FORMAT(1H ,'MATRIX ',A5,' REPLACED'/)
00440 GO TO 20
00450 25 K=0
00460 NMB=NMB+1
00470 DO 30 I=1,NS
00480 IF ((NP+I-K).LT.129) GO TO 40
00490 NREC=NREC+1
00500 IERR=DISKIO(0,NFILN,NREC,FBUF(1+INC))
00510 NP=1
00520 K=I
00530 40 A(I+INCA)=FBUF(INC+NP+I-K)
00540 30 CONTINUE
00550 NP=NP+I-K+1
00560 IF (NP.LT.129) GO TO 45
00570 ISET=1
00580 GO TO 11
00590 C * * * IF THERE ARE MORE MATRICES, READ & LIST THEM * * *
00600 45 ISET=0
00610 IF (FBUF(INC+NP).EQ.SLASH) GO TO 60
00620 GO TO 10
00630 50 TYPE 100,IERR
00640 100 FORMAT(1H ,'ERROR ',I3/)
00650 IF (IERR.EQ.2) TYPE 300,NFILN
00660 300 FORMAT(1H ,'FILE ',2A5,' NOT FOUND'/)
00670 IERR=0
00680 GO TO 70
00690 60 IERR=DISKIO(2,NFILN,NREC,FBUF(INC+1))
00700 IF (IERR.NE.0) GO TO 50
00710 70 CALL FRET(128,IPT)
00720 TYPE 200,NMB,NFILN
00730 200 FORMAT(1H ,I3,' MATRICES RETRIEVED FROM DISK FILE ',2A5/)
00740 RETURN
00750 END
00010 C * * * * * * SUBROUTINE TO SAVE CURRENT MATRICES * * * * * *
00020 SUBROUTINE MSAVE(NAM,NV)
00030 INTEGER FREZ,FRET,DISKIO,ERROR,FRES
00040 COMMON /FF/ NFMX,ITAPE,NWMX
00050 COMMON /XX/ ERROR
00060 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00070 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1)
00080 COMMON /ZZ/ NUM,NNUM,NX,NR,NC,MPYTAG
00090 DIMENSION NFILN(2),FBUF(1),NW(1),FNUMB(1)
00100 DATA NFILN(2) /5H .MAT/,IBL/' '/
00110 IPT=FREZ(128)
00120 INC=IPT-LOC(FBUF)
00130 NREC=1
00140 NP=1
00150 NFILN(1)=NAM
00160 IF (NV.EQ.NUM) GO TO 9
00170 C * * * * * * SAVES "N" MATRICES * * * * * *
00180 TYPE 1
00190 1 FORMAT(1H ,'ENTER MATRICES TO SAVE:'/)
00200 NWP=FRES(NV,IBL)
00210 INW=NWP-LOC(NW)
00220 3 NA=NV
00230 NF=0
00240 NWMX=NV
00250 NFMX=1
00260 CALL FREAD(NA,NW(1+INW),NF,FNUMB)
00270 IF ((NA.EQ.NV).AND.(NF.EQ.0)) GO TO 5
00280 TYPE 4,NV
00290 4 FORMAT(1H ,'ERROR. ENTER ',I3,' MATRICES TO BE SAVED.'/)
00300 GO TO 3
00310 5 DO 7 I=1,NV
00320 DO 6 J=1,NUM
00330 6 IF (NAME(INC1+J).EQ.NW(I+INW)) GO TO 8
00340 ERROR=2
00350 GO TO 20
00360 8 CALL SAVE(J,FBUF(INC+1),NFILN(1),NREC,NP)
00370 7 CONTINUE
00380 GO TO 12
00390 C * * * * * * SAVES ALL CURRENT MATRICES * * * * * *
00400 9 DO 10 I=1,NUM
00410 CALL SAVE(I,FBUF(INC+1),NFILN(1),NREC,NP)
00420 10 CONTINUE
00430 12 FBUF(INC+NP)='//'
00440 NP=NP+1
00450 IERR=DISKIO(-1,NFILN,NREC,FBUF(1+INC))
00460 IF (IERR.NE.0) GO TO 30
00470 IERR=DISKIO(-2,NFILN,NREC,FBUF(1+INC))
00480 IF (IERR.NE.0) GO TO 30
00490 TYPE 200,NV,NFILN
00500 200 FORMAT(1H ,I3,' MATRICES STORED IN DISK FILE ',2A5/)
00510 GO TO 20
00520 30 TYPE 100,IERR
00530 100 FORMAT(1H ,'ERROR ',I3)
00540 20 IF (NV.NE.NUM) CALL FRET(NV,NWP)
00550 22 CALL FRET(128,IPT)
00560 RETURN
00570 END
00010 C * * * * SAVES ONE MATRIX IN DISK FILE NFILN * * * *
00020 SUBROUTINE SAVE(I,FBUF,NFILN,NREC,NP)
00030 COMMON /YY/ INCA,INC1,INC2,INC3,INC4
00040 COMMON /Z/ A(1),NAME(1),NCOL(1),NROW(1),NN(1)
00050 DIMENSION FBUF(1),NFILN(2)
00060 INTEGER DISKIO
00070 ISET=0
00080 IF ((NP+2).LT.129) GO TO 30
00090 5 IERR=DISKIO(-1,NFILN,NREC,FBUF)
00100 IF (IERR.NE.0) GO TO 20
00110 DO 10 J=1,128
00120 10 FBUF(J)=0.0
00130 NREC=NREC+1
00140 NP=1
00150 IF (ISET.EQ.1) GO TO 60
00160 GO TO 30
00170 20 TYPE 100,IERR
00180 100 FORMAT(1H ,'I/O ERROR ',I3)
00190 GO TO 60
00200 C * * * * SAVE NAME, # OF ROWS, # OF COLS * * * *
00210 30 DO 31 J=1,5
00220 CALL GET(NAME(I+INC1),J,IC)
00230 CALL PUT(FBUF(NP),J,IC)
00240 31 CONTINUE
00250 FBUF(NP+1)=NROW(I+INC3)
00260 FBUF(NP+2)=NCOL(I+INC2)
00270 INCA=NN(I+INC4)-LOC(A)
00280 NP=NP+2
00290 NS=NROW(I+INC3)*NCOL(I+INC2)
00300 K=0
00310 C * * * * SAVE MATRIX ELEMENTS * * * *
00320 DO 50 J=1,NS
00330 IF ((NP+J-K).LT.129) GO TO 40
00340 IERR=DISKIO(-1,NFILN,NREC,FBUF)
00350 IF (IERR.NE.0) GO TO 20
00360 DO 35 L=1,128
00370 35 FBUF(L)=0.0
00380 NP=1
00390 NREC=NREC+1
00400 K=J
00410 40 FBUF(NP+J-K)=A(INCA+J)
00420 50 CONTINUE
00430 NP=NP+J-K+1
00440 IF (NP.LT.129) GO TO 60
00450 ISET=1
00460 GO TO 5
00470 60 ISET=0
00480 RETURN
00490 END
00010 C SUBROUTINE FIXEM -FIXED END MOMENTS FOR MATRIX 20,11
00020 SUBROUTINE FIXEM(NTIME)
00030 COMMON /YY/ INCA
00040 COMMON /Z/ A(1)
00050 DIMENSION DATA(3)
00060 TYPE 10
00070 10 FORMAT(2X,'LOAD,',1X,'LENGTH'/)
00080 DO 100 I=1,NTIME
00090 20 NA=-2
00100 NF=2
00110 CALL FREAD(NA,NA,NF,DATA)
00120 W=DATA(1)
00130 XL=DATA(2)
00140 MEMBR=(I-1)*4+INCA
00150 WXL=W*XL/2.
00160 WXL2=W*XL**2./12.
00170 A(MEMBR+1)=-WXL
00180 A(MEMBR+2)=WXL2
00190 A(MEMBR+3)=-WXL
00200 A(MEMBR+4)=-WXL2
00210 GO TO 100
00220 GO TO 20
00230 100 CONTINUE
00240 RETURN
00250 END