Home

NASA Technical Memorandum X-2417

1. NP4 N M CALL TRANP G NG Z NP4 amp NNP 4 CALL MULT RI NRI Z NP4 NNP4 CALL 7 NNP4 IC II E 1 GN TO 204 CMLL MULT Q NO H4 NH 7 1 NNP1 CALL MULT Z NP3 NNP3 71 1 1 Z 2 NNP2 GO TO 205 4 19 72 205 NPATR MOD N 2 IF NPATIR FO 0 NPAIR TWO NN 1 N 2 1 GO TO 201 202 201 NO 300 Iz1 N 2 2 Nx N I 1 NTH3 TW ON 1 1 M 1 911 300 CALL EQUATE Z NP2 NN 7 NTH3 NN DO 302 Iz2 N 2 NP4zN 3 N 51 NTH2 T WU N N 1 302 CALL EQUATF Z NP4 NN 7 2 NN GO TO 202 203 NF A TR 202 D 301 I 2 N 2 NP2 1 1 MI LL 2 I VIC 2 M T1 1 N 1 201 CALL EOUATF Z NP2 NN Z NIH3 304 J 1 N 2 NPA42N 3xXN T 1 4 NTH2 TMOU EN CN 1 304 CALL EQOUATF Z NP4 NTH2 NN GO TO 203 201 NPAIR 203 DO 303 Iz1 N II I N NO 303 J 1 N JJ J N L1 NNZ J 1 L2 NNZ 1 1 JJ L3 N J 1 I 7 L1 F L3 303 L2 F L3 60 TO 1000 4 995 CALL LNCNT 2 MRITF 6 50 NF yNG yNRI NH NO 50 FORMAT DIMENSION ERROR IN AUG 135 NF1 9X NG 9X NRI 8X 1 NH 9 0 29 5 3 216 999 CALL ASPERR 1000 RETURN E ND LII SUBROUTINE RICAT PHI NPHI NCINCONT Ky PT KDB
2. STEP NUMBER 2 IN SAMPL P T MATRIX 3 ROWS 3 COLUMNS M 3 3 1702128D 00 1 00000000 00 5 3191 4850 01 1 0000000 00 1 00000000 00 0 0 5 23191489D 01 OO 5 4 8 2 3 49 D 949 F Q r P nF F F STEP NUMBER 2 IN SAMPL K T MATRIX 3 ROWS 2 COLUMNS 4 1054403D 01 5 27201350 02 ___ ______________________________________________ lt _ 3 05103760 01 1 52551880 00 T m m l I n STEP NUJMBERBZ IN SAMPI a a 22 22 P T MATRIX 3 ROWS 3 COLUMNS 3 378931995 086 3 000090068 06 6 40207528 03 o V 1 0000000D 00 1 0000000D 00 0 0 6 10207520 01 0 0 6 49186760 00 M M STEP NUMBER 3 IN SAMPL M A HA F T n EE e XII MATRIX 3 ROWS k UMS n nT 4 0964801D 01 4 8240037D 02 0 0 0 0 uices nnt STEP NUMBER 4 IN SAMPL P T MATRIX 3 ROMS 3 COLUMNS 3 18070400 00 1 0000000D 00 6 2633
3. a X B R gt _ _ _ END lt lt _ _ __ _ gt _ _ UNLOAD VASPMN ERASE SOURCE VASPMN SOURCE MNPG USERLIB VASPMN ERASE USEREtBOMNP8 3 RELEASE VASP1 DISPLAY ALL CONVERSATIONAL VASP PROGRAMS CLEARED ase YOU MAY RESTART WITH 2 PROCDEF CMPL AAA C s nn I cma 01980 5 1 1 5 02000 02100 02200 02500 02h00 02500 02600 02700 02800 02900 05000 4 as de A os ee M r DEFAULT SYSINX E EDIT SOURCE MNPG EXCERPT SOURCE VASPMN 1600 1700 800 1 s masse odis END DISPLAY MNPG NOW COMPILING DEFAULT LFMEN N UNLOAD MNPG ERASE USERLIB MNPG FTN MNPG Y s u Ink LOAD MNPG _ _ gt _3 _ 2 Figure 12 List of data set VASPPROC 03100 CALL MNPG 03200 DEFAULT LIME
4. 0028 CALL ETPHI FSTAR NFSTAR DELT2 A2 NA2 DUMMY KDUM a User s main program Figure 3 Example 3 UTE 0029 CALL LNCNT 2 onec WRITE 6 106 I 0031 106 FORMAT 10 TRANSITION MATRIX 0032 CALL A2 NA2 A2 1 003 CALL LNCNT 3 0035 WRITE 6 101 0036 CALkL MULTLCK NCK XT NXT DELTC NDELTC J 0037 DELTC 1 1 0 DELTC 1 MON 0038 CALL LNCNT 1 0039 WRITE 6 102 TT XT DELTC 0041 CALL MULT A2 NA2 XX NXX XT NXT 0042 CALL MULT CK NCK XT NXT DELTC NDELTC 0043 DELTCLI e 21 DELTCLA 2 2 0044 TT TT DFLT2 0045 WRITE 6 102 TT XT DELTC Faray EE L GET EHE 0047 GO TO 200 0048 END a User s main program Concluded Figure 3 Continued bb TEST PROGRAM 2A GENERATES TRANSIENT USING X I 1 EXP F T X I 0 0 1 0 0 01 3 5 L 3 3 0 2767 1 0 0 0372 17 0872 0 1785 12 1983 0 0 0 0 6 67 3 1 0 0 0 0 6 67 X0 3 1 1 0 0 0 0 0 3 3 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 Q 3 3 0 2 0 0 0 0 0 0 _ 0 2 0 0 0 0 0 0 0 0 R 1 1 3 3 0 0 0 0 0 0 0 0200 _ 0 0 0 0 0 0 0 0 b Data deck Figure 3 Continued N St TEST PROGRAM 2A GENERATES TRANSIENT USING X I 1 2EXP FxT Xxx I F MATRIX 3 ROWS 9 4610000D 01 1 0000000 5 00 1 7087200D 01 1 7850000D 01 0 0
5. Current Pivot AIT Figure 7 Information Systems Co flow chart subroutine PSEU A B C EE DEP JP D gt Rank 1 Yes Save Indow Qf Tris Pivot Set New Tnrethold DPRI 2x 105 Current Pivot No IR NR No ce Ye Quit 799 e Re invert B Here Compute Set IR Rank Matrix U From C E E Store into C Replace B Diagonal With Ones And Zeros If Zero In Diagonal Sat M 1 EE x Ct Re scan Diagonal Yes e E ts This DECOM Cali QDCM Rediagonalization No Sat JP3 No Iteration Call ANDRA B EE DRP JP 5 lt QR2 gt 2082 Error Error No C CALL TTRM EE EE BxC 691 692 693 606 Figure 7 Concluded TTRM NR C EE EE 61XC PSEUP Symmetric Use At as Needed to Form Non symmetric Pseudo inv in 8 Set Rank in Return Parameter Move Answer EE to B 29 BDNRM DESCRIPTION This subroutine computes the quantity norm QPQ Q norm Q where the values of P and Q are in the square arrays CT and EE or EE and CT depending on the sign of NR If P the return value is zero This routine can thus be used to test the quality of a pseudoinverse USAGE CALL BDNRM NR CT EE D KRV Input Arguments Matrices CT EE with dimensions NR X NR Constants NR size of
6. rr rs a mm ZL MATRIX 2 ROWS 2 COLUMNS 60 00 AAAA ASAA A 2 00000000 00 2 00000000 00 RN en es 0 0 0 0 0 0 0 0 EXR MATRIX 1 ROWS 1 COLUMNS 0 0 EX MATRIX 1 ROWS 1 COLUMNS 0 0 T MATRIX 2 ROWS 2 COLUMNS 1 00000000 00 0 0 A O00 0000D LL T T MATRIX 2 ROWS 2 COLUMNS 350050000050 UU 0 0 1 0000000D 00 _ R MATRIX 2 ROWS 2 COLUMNS 2 ee 1 0000000D 00 1 0000000 00 2 00000000 00 2 00000000 00 L suka u PS R1 MATRIX 2 ROWS 2 COLUMNS 1 00000000 00 1 0000000D 00 00 06 2 0090090990p 00 lt lt RERR MATRIX 2 ROWS 2 COLUMNS 463336396 4 4408921D 16 4 4408921D 16 d Case 1 Figure 5 Continued ee ma 1 MS SS 1 4142136D 00 2 8284271D 00 MATRIX 1 ROWS 2 COLUMNS 7 07106780 01 17 07106780 01 RANK MATRIX 1 ROWS COLUMNS re m aM M _ 1 0000000D 00 d Case 1 Concluded Figure 5 Continued 65 09 TEST PROGRAM FOR DECGEN AND DECOM CASE 4 2X kaNK2 VASP PROGRAM
7. GO CCAEE FOUATREITUSNZESRSNR 0028 IF DUM M7 NE 0 DO GO TO 50 0029 CALL DECSYM H 0020 M I program Figure 5 Example 5 VI N 0031 50 CALL DECGEN R NR G NG H NH DUM KDUM 0032 CALL MULT 1 1 0033 RyNRyRLyNRLyREyNRE CUCU gt 0034 CALL PRNT R NR R 1 0035 CALL PRNT RI NRI R1 1 0037 CALL PRNT H NH 1 1 0038 CALL G NG G6 1 1 0 0 3 9 NB H1 0040 1 0 0041 CALL PRNT DUM M6 ND RANK 1 0042 20 0232 E hp nM P nT J N O P a Main program Concluded Figure 2 Continued SS 95 n F n O OA nm n FV k 0001 SUBROUTINE ORTH T NT DUM KDUM 0002 DOUBLE PRECISION T 1 DUM 1 STH 000 3 DHAENS LON N E a 0004 LDM2 2 NT 1 XX2 1 0005 10 LDM1 NT 1 XX2 1 NEN S T 0007 CALL UNITY DUM LDM1 NT 0008 NM NT 1 1 0009 B0 20 LM 0010 JP J 1 001
8. 03300 DISPLAY w w COMPUTING DONE ww 05540 PROCDEF CONVASP 052580 C XCISE SS SS Se 03400 PROCDEF CONVASP 03500 PARAM N A B C W X Y Z 053600 DDEF VASP1 AP VASE 6 PT ONEUJOBE48 03700 JOBLIBS SYSULIB 02308 DISPLAY tee MATRTOES AVATLABEL AELCDIMENSTONED N ARE 0 05900 DISPLAY A 8 C W X Y 2 FOR INPUT OR COMPUTATIONS 04000 DISPLAY D1 D2 D3 D4 D5 D6 D7 FOR COMPUTATIONS ONLY 24H 099 DEFAULT SYSTNXCRE mn 04200 DEFAULT LIMEN N 04300 EDIT SOURCE VASPMN 01Hr06 EXCFSE EEAST nF 04500 INSERT 100 04600 IMPLICIT REAL 8 A H 0 Z ASP 04 709 OMMON AC 81 58CH07 C 8N7 04800 1 YC N 2 N D1C N D2 N D3 ND D amp N D5 N D6 ND D7 ND 04900 2 N A 2 N B 2 N C 2 N W 2 N X 2 05000 5 N Y 23 N 2 2Y ND1 2 ND2 2 ND5 2 NDt 2 ND5 2 ND6 2 ND7 2 05100 COMMON MAX MAXRC 05200 MAXRC N 5300 210 PRINT IS e ee eee 05400 15 FORMAT DATA 05500 CALL INPUT A A B B C C W W 05608 SX e 1 2 A Nt A S 7 05700 2 N 8 N B NSC N C NSW NSW NSX NSX 05800 5 N Y N Y N 2 NSZ 05809 EH2 FORMAT CHC IPEA
9. AASA gt AAA AERR gt A or Be ASAA gt A AA It be noted that the size of the numbers the AERR matrices is 1016 which is very good Ames Research Center National Aeronautics and Space Administration Moffett Field Calif 94035 July 12 1971 64 S9 000 1 DIMENSION A 50 B 50 VI 250 2 NB 2 yA1 50 2 50 NA1 2 1 NA2 2 LAB 2 NOO2 DOUBLE PRECTSTON A RoW Ad 2 0003 COMMON MAX MAXRC 0004 ios nn MEAXROSSEDS Lun 0006 5 CALL RDTITL 0007 CALL READ 1 D 0 0 8 NW 3 5 0 0009 CALL PSEUDO A NA NB W NW 0010 CALL PRNT 5 1 0011 CALL NA 0012 CALL 1 1 2 2 0013 CALL SCALE A NA A1 yNA1 1 D0 0015 CALL PRNT A2 NA2 AASA 1 0016 CALL PRNT Al NAl AERR 1 A ee CALL MUL T4 B4 NB 1 1 0018 CALL MULT A1 NA1 B NB A2 NA2 0019 CALL SCALE B NB A1 NA1 1 D0 0020 H GALL ADD AE NaH Ae He 0021 CALL 2 2 5 1 0022 CALL PRNT Al NA1 AERR 1 0023 unc a ee Sei Se es um ar usd 0024 END a Main program to check PSEUDO Figure 6 Example 6 rr oS 99 PSEUDO TFS
10. 0 R 1 1 0 K 2 5 1 Oe s 5 O 2 0 36 0 1 5 1 0 0 0 0 0 1 0 0 0 0 ls 0 0 Oe Oe 0 0 0 0 1 l 0 0 0 0 0 1 O 1 0 X 5 1 Le l 1 0 le T 4 1 5 2 3 5 0 b Data deck Figure 2 Continued Le TEST PROGRAM GENERATES TRANSIENT USING TRNSI VASP PROGRAM F MATRIX 5 ROWS 5 COLUMNS 0 0 0 0 0 0 5 0000000D 0 1 0 0 0 0 0 0 0 0 0 0 1 0000000D 00 100000000 00 0 0 0 0 0 0 0 0 1 0000000D 00 0 0 0 0 0 0 0 0 0 0 2 0000000D 00 METRIA 5 RUNS 7 2 8 PH VJ Lj 1 00000000 00 0 0 1 0000000D 00 1 00000000 00 0 0 0 0 1 00000000 00 0 0 0 0 1 0000000D 00 J MATRIX 2 ROWS 1 COLUMNS 0 0 0 0 R MATRIX 1 ROWS 1 COLUMNS 0 0 K MATRIX 2 ROWS 5 COLUMNS 1 0000000D 00 0 0 5 0000000D 01 0 0 2 0000000D 00 0 0 3 00000000 00 0 0 1 0000000D 00 0 0 MATRIX 7 ROWS 5 COLUMNS 1 00000000 00 0 0 0 0 0 0 0 0 0 0 1 00000000 00 0 0 0 0 0 0 0 0 0 0 1 0000000D 00 0 0 0 0 0 0 0 0 0 0 1 0000000D 00 0 0 0 0 0 0 0 0 0 0 1 0000000D 00 1 0000000D 00 0 0 0 0 0 0 0 0 0 0 1 00000000 00 0 0 1 0000000D 00 0 0 Output Figure 2 Continued C TEST PROGRAM GENERATFS TRANSIENT USING TRNSI VASP PROGRAM UN 1 0000000D 00 1 0000000D 00 1 0000000D 00 0 0 1 0000000D 00 MATRIX 4 ROWS 1 COLUMNS 5 00000000 01 2 00000000 00 3 5000000D 00 0 0
11. CALL EAT F NF TI DUMMY NNF DUMMY N4 NNF DUMMY 11 KDUM H 4 CALL PRNT DUMMY N3 NF 1 CALL EQUATE DUMMY N3 NF DUMMY 1 CALL PRNT DUMMY N4 NF INT 1 DEDE CALL MULT DUMMY N4 NF G NG DUMMY N2 CALL MULT J NJ RyNR DUMMY L 3 NUI CALL LNCNT 100 7 CALL LNCNT 3 WRITET6 501 50 FORMAT 1HO TRANSTENT RESPONSE INDICATES CONTROL CHANGES WRITE 6 51 11 1 51 FORMAT 1HO 4X TIME FIRST 13 EL NTS CONTAIN X NEXT I3 T 1 ELEMENTS CONTAIN Y HX LAST I3 ELEMENTS CONTAIN U JR 2 2 0 MU aQ IU 1 CALL MULT KyNKy X NX DUMMY L5 NU 9 DUMMYILS5S NU 203 CALL MULT H NH X NX DUMMY L4 NY IF IU NE 11 GO TO 205 WRITE 6 101 AR X TP 1 PJ 1 9 101 FORMAT 1HO A1 F8 2 1P7D15 7 10X 1P7015 7 GU 206 205 WRITE 6 102 IP 1 NXR DUMMY TP IP L4 LAST 102 FORMAT 1HO 1X F8 2 1P7D15 7 10X 1P 7D15 7 206 JU 11 10 IF DABS T 3 RETURN TT TT 202 CALL MULT DUMMY I1 NNF X NX DUMMY L6 NNX CALL MULT DUMMY N2 NG DUMMY L5 NU X CALL ADD X NX DUMMY L6 IF NO GO TO 201 1 60 203 C C DIMENSTON ERROR DIAGNOSTIC 999 WRITE 6 990 990 FORMAT 1HO DIMENSION ERROR TRNSI 25X COL SIZ
12. 2 7 2 o SUBROUTINE PRNT AR NAR NAM TP gt C SUBR PRNT PRINTS DOUBLE PRECISION MATRIX COMMON 6 2 6 COMMON L TNES NLP L IN TITLE 23 COMMON C NOTE NO LTNES PAGE VARIES WITH THE INSTALLATION DATA KZ KW KB 1H0 1H1 1H REAL 8 AR n DIMENSION AR 1 4 NAR 2 E NAM l HEADLINE SAME IF 2 HEADLINE NEW PAGE C IP 3 NO HFADLINE SAME PAGE IP 4 NO HEADLINE NEW 11 IP NR NAR 1 NC NAR 2 NLST NR MC IF NLST GT ORe NLST Lie 1 0R NR LT 1 GO TO 16 IF NAME NAME KB SKIP HEADLINE IF REQUESTED 7 10 CALL LNCNT 100 11 CALL LNCNT 2 3 WRITE 6 177 7 177 1 5 4 8 MATRIX4 5X 13 5H ROWS 5X 13 8H COLUMNS GO TO 13 12 CALL LNCNT 100 GO TO 13 132 CALL LNCNT 2 WRITE 6 891 891 FORMAT 1H0 C BELOW COMPUTE NR OF LINES ROW 13 J NC 1 NEPR 1 NLPM J JST 1 C COMPUTE LAST ROW POSITION 1 NLST NLST NR MN NC NC GT NFPR MN NFPR KLST NR MN 1 B V 56 91 CONTINUE 00 912 J JST NR CALL LNCNT NLPW KLST KLST 1 WRITE 6 1 AR N N J KLST IF _ NC LF NFPR 60 912 NLST NLST KNR KL ST NR 91 2 RETURN 16 CALL LNCNT 1 WRITE 6 916 916 FORMAT ERROR IN PRNT MATRIX A4 HAS NA 216 CALL ASPERR RETURN FND SUBRO
13. C6CM DEC lt A REAL GENERAL MATRIX SUBRRUTINF DECGEN R NR G NG H NH4DUM KDUM 0 DIMENSION NR 2 NH 2 NG 2 ND 2 NAC2 KP CA COMMON MAX MAXRC DOUBLE PRECISION R 1 G6 1 H C 1 DUR 1 DCR 3 NJ NR 1 NR 2 1 CALL TRANP R NR DUM NJ ND IF NR 1 GT NR 2 GU CALL MULT R NR DUM NJ ND DUM ICS 1 GO TO 30 10 CALL MULT DUM NJ ND R NR DUM ICS 2 GO TO 30 C ENTRY DECSYM DECOMPOSE A REAL SYMETRIC MATRIX M NG H N I ICS 0 IF NR 1 NE NR 2 TO 601 CALL EQUATE R NR DUM C DUMMY STORAGE MUST AT LEAST 2 30 M NA 1 MMS KDUM LT MM7 Gn TO 601 KP 1 20 KP 2 z0 KP 3 M KP 4 M NCM 1 0 DCM 23 0 MS M XM N2 MS N3 N2 MS N4 N3 MS INNING N6 N5 MS M 7 M6 M5 6 CALL TRANP DUM N3 NA DUM N6 ND IF ICS EQ 2 GOTO 100 ALL MULT DUM N2 NA DUM N6 ND NH 2 KP 2 IF ICS F0 1 GO TO 60 _ _ CALL TRANP HyNH G NG 60 150 60 CALL MULT 5 CALL TRANP 6 NG DUM N6 ND CALL MULT DUMING ND R NR G NG GO TO 150 100 CALL MULT DUM N2 NA DUM N6 ND H NH NG 2 KP 2 CALL TRANP H NH G NG CALL MULT G NG DUM N5 NA H NH CALL TRANP HzNH DUM N6 ND CALL MULT R NR DUM N6 ND H NH 150 DUM 6 KP 2
14. DI1 D2 D5 D5 D5 D6 D7 FOR COMPUTATIONS ONLY T CONVASP 0000200 DEFAULT SYSINX E CONVASP 0000250 DEFAULT LIMEN N Q00NVASP 0000300 EDTT SOURCE VASPMN S TT CONVASP 0000350 EXCISE 1 LAST CONVASP 0000380 INSERT 100 2 CONVASP 0000400 IMPLICIT REAL 8 A H 0 2 7 IN ISIKA CONVASP 0000430 COMMON ASP A N B N C N AWC N CONVASP 0000460 1 2 01 02 05 0 05 06 07 GNVASP 6606H96 N A 2 N 802 N CC20 N WC2 N XTZ CONVASP 0000520 5 N Y 2 N Z 2 ND1 2 ND2 2 ND3 2 ND4 2 ND5 2 ND6 2 ND7 2 CONVASP 0000550 COMMON MAX MAXRC CONVASP 0000580 MAXRC N gt CONVASP 0000700 10 PRINT 15 CONVASP 0000800 15 FORMAT DATA EC CONVASP 0000800 CALLE TNPUT TSA A TSBTSSB SCT ST 5W W CONVASP 0000880 1 X X 5 Y Y 2 2 N A N A CONVASP 0000920 2 N B N B N C N C NSW N X N X Q 00 CONVASP 0000950 3 NSY N Y N Z N Z CONVASP 0001000 13 FORMAT 1X 1 20 13 CONVASP 0001100 6 FORMAT 1X 1P6D13 6 m 2 0001200 RETURN acc CONVASP 0001300 END CONVASP 0001500 END ee CONVASP 0001450 ERASE USERLIB VASPMN i CONVASP 0001500 FTN VASPMN Y CONVASP 0001600 XLIST VASPMN
15. DOUBLE PRECISION A 25 DOUBLE PRECISION A 5 5 DOUBLE PRECISION A 3 10 DOUBLE PRECISION A 100 The important factor is the total number of cells reserved and the user may reserve more cells than the matrix requires or conversely he may put a smaller matrix than originally planned in a specific array The VASP program stores data in an array as a string of columns just as Fortran does The convention used here and throughout this report is that the name of a dimension array is obtained by prefixing the letter N to the matrix name 4 b However it stores the matrix according to the dimension given in NA whereas Fortran stores A according to the dimensions in the Fortran dimension statement 2 Consider the following example The Fortran statements are DOUBLE PRECISION 5 5 5 5 5 5 DIMENSION 2 NB 2 NC 2 CALL READ 3 A NA B NB C NC The first card in the data deck specifies NA 5 5 followed by cards with 25 data words for A then NB 4 4 followed by 16 data words of B finally NC 6 6 followed by 36 data words of C Since the storage of data in VASP is controlled by the VASP dimensioning the 25 words for A will exactly fill the reserved storage and the 16 words for B will fill the first 16 cells of the storage reserved for B The 36 words for C will completely fill the reserved storage for C and overflow into something else The user can prevent this overflow by set
16. RETURN 601 CALL LNCNT 1 WRITE 6 1601 MM7 1601 FOURMAT DIMENSION ERROR IN DECSYM DECGFN NRz 2I6 3X 1 KDUMz I4 3X KDUMI M IN 2 I4 CALL ASPERR RETURN END 0CI SUBROUTINE READ1 COMMON M X M DIMENSION 2 2 DOUBLE PRECISION IF 7 1 0 0 GO 410 NR NZ 1 MC M7 2 NLS T NR NC IF NLST MAXRC NLST LT I 0R NR LT 1 GO TO 16 DO 400 1 1 NR 400 READ 5 101 J J T NLST NR NA 1 zNR NA 2 NC 410 CALL PRNT 1 101 FORMAT 8F10 2 RETURN 16 CALL LNCNT 1 WRITE 6 916 NAM NR NC 916 FORMAT ERROR IN READI MATRIX A4 HAS 216 RETURN END OAL NI 100 SUBROUTINE ASPERR DATA 1 10 CALL TRACE ERRTRA IS THE 360 67 TRACE ROUTINE TRACE IS FOR TSS CALL ERRTRA THIS IS AN INSTALLATION DEPENDENT SUBROUTINE SUBROUTINE ERRTRA IS A SUBROUTINE SUPPLIED BY THE AMES OPERATING SYSTEM TO PROVIDE AN ERROR WALKBACK THE STATEMENT CALL ERRTRA SHOULD BE EITHER 1 CHANGED TO MATCH THE USERS OPERATING SYSTEM OR 2 DELETED ALTOGETHER 1 1 1 IF 1 61 0 RETURN 10 WRITE 6 100 FORMAT TOO MANY ERRORS EXIT CALLED CALL EXIT RETURN END cel COMMON 6 2 6 COMMON LINES NLP LIN TITLE 23 COMMON MAX MAXRC DATA MAXRC 6400 NOTE NLP LINES PAGE VARIES WITH THE INSTALLATION DATA LIN NLP 15 45 DATA
17. _ _ ___ _ _ e O0NVASP 0001700 LEOAD VASPMNSS vvv mU CONVASP 0001770 DEFAULT LIMEN N CONVASP 0001800 EDIT SOURCE MNPG 0 5 0001900_ 1 1 1457 CONVASP 0001950 INSERT 1 1 CONVASP 0002000 EXCERPT SOURCE VASPMN 100 1500 _ CONVASP 0002100 DISPLAY FORTRAN STATEMENTS CONVASP 0002150 DEFAULT LIMEN W CONVASP 0002200 INSERT 100 Figure 11 List of VASP procedures Continued 191 PROCDEF RECMPT DEFAULT LIMEN N CALL MNPG DEFAULT LIMEN W DISPLAY COMPUTING DONE PROCDEF REWRT PARAM LINE DEFAULT LIMEN W DEFAULT SYSINX E EDIT SOURCE MNPG EXCISE LINE LAST IF LINE 1100 DISPLAY FORTRAN STATEMENTS SLINE 11001 LIST 100 LAST DEFAULT SYSINX G INSERT LINE PROCDEF VASP PARAM INPUT OUTPUT DEFAULT N 9 A A B B C C WzW X X Y Y Z Z LINE 100 JBLB VASP VASP 0000300L0AD BLKDTA RECMPT 0000000 RECMPT 0000100 RECMPT 0000200 RECMPT 0000300 RECMPT 0000400 REWRT 0000000 REWRT 0000100 REWRT 0000200 REWRT 0000400 REWRT 0000500 REWRT 0000600 REWRT 0000700 REWRT 0000800 REWRT 0000900 REWRT 0001000 YASP 0000000 VASP 0000100 VASP 0000150 VASP 0000200 VASP 0000400 VASP 0000500 VASP 0000600 VASP 0000700 IF INPUT 7s DDEF S INPUTS DISPLAY t INPUT FROM DATA SET INPUT IF INPUT
18. smaller At the 6 At is used again form the pseudoinverse Square uses AtA The result 75 is always forced symmetric afterward even for symmetric entry ANDRA is called to diagonalize this result in B Most of the pivots are found and the steps made on the first call If not iterating this part is not repeated If singular rank of symmetric input not maximal a transforming matrix is computed A copy of the original symmetric B is transformed and reinverted by ANDRA The result is retransformed by premultiplication and postmultiplication If iterating the pivot tolerance is decreased and ANDRA is called to find one pivot at a time A side calculation is done to measure the quality of pseudoinverse formed at each step The routine backs up one step and stops with rank one less if it makes a bad step The result if singular is sent through the reinversion above The use of PSEUP by DECOM avoids reinverting in the singular case also it never uses a nonsquare input There is a find rank only option If PSEU is used without iteration four I O matrices are needed plus a dummy matrix Iteration uses four additional dummy matrices Iteration cannot be done during the reinversion Besides those mentioned entries BDNRM MULT and NORM are used for iteration TTRM is also used except in rank only case ANDRA diagonalization algorithm For a detailed description of the method see the documentation of ANDRA itself A mathem
19. DIMENSION MC 0II 31 2 yNC 2 MC 2 MXL21 M9L21 106IL21 DIMENSION PHI 1 C 1 K 1 PT 1 1 DOUBLE PRECISION PHI C K PT SUM SUMN AL Wa DET TWO 2 N NPHI 1 TWO 5 0 N N 504 4 50 NP1 1 2 M50 M01 NP3Z NSO4NP2 NP4 NSQ NP3 IF KDUM LT NSQ amp GN TO 600 IF NPHI 2 NE NPHI 1 0R NC 2 NE N UR NPT 1 NE N UR 2 1 IC 600 NPRINT NC ONT 1 NBL OCK NC 2 NPRINT NZ NC ONT 3 C REARRANGE PHI MATRIX NN 1 N NN 2 1 DO 300 I 1 N NTH1 TWU N 1 1 NTH3 NTH1 N NM1 1 1 1 MV 2 M M CALL EQUATE PHI NTH1 W NW1 NN 30 A QUATE PHT NTH3 NNM W NW2 NN NM 1 TWD N N NM 2 1 CALL EOUATE W C 1 NM PHT I NM DO 301 I 1 N NTH2 TWOXNX N I 1 1 NT H4 NTH2 N NW3 N IVIIIXM I 1 1 NW4z NM3 N N PHTI NTH NN MW NM3 NN 4 1 301 CALL EQUATE PHI NTH4 NN VI N WA 1 CALL EQUATF W NWW NM PHI1 NWW EVUATT VVVIS WW OEE TT MAIN COMPUT ATION CALL UNITY PT NPT DO 306 I 1 N dh rs 306 K I O DO 403 NBLQCK 202 201 DO 104 Js1 NPRINT CALL MULT PHI NP3 W 1 NPT CALL ADD PHI 1 1 Wii NPT CALL INV W 1 W NP2 CALL MULT PHI NP4 W NP2 CALL
20. DISPLAY INPUT FROM TERMINAL IF 0UTPUT 7 ait DDEF FT06F001 0UTPUT DISPLAY OUTPUT PLACED DATA SET OUTPUT IF OUTPUT eT DISPLAY OUTPUT TO TERMINAL Figure 11 List of VASP procedures Concluded 91 00050 00060 LOGON USERID 9 PROCDEF CHNGIN 4 MEN 00090 EXCISE 1 LAST 00100 PROCDEF CHNGIN 09200 PARAM 1NPUT t m mm 6 a A M 8 F 00500RELEASE 05 001 00500 00500 00540 14 INPUT pDEF FTO5F001 INPUT DISPLAY INPUT FROM DATA SET 1NPUT ti DISPLAY LINPUT FROM TERMINAL pm CHNGOUT M 00580_EXCISE 1 LAST 00600 PROCBEF CHNGOUT 00700 00800 ee oe PARAM 0UTPUT RELEASE 00960 1F OUTPUT s pDDEF FT0G6F001 0UTPUT DTISPEAY OUTPUT PLACED 1N DATA SET 0UTPUT 01000 01040_ 01680 FX TSE gt AS PROCDEF CLRVASP 01100 01200 01206 01400 01500 01600 01700 01800 01900 01940 DISPLAY IF OUTPUT PROCDEF CLRVASP DISPLAY OUTPUT TO TERMINAL
21. Ff MATRIX BROWS 5 COLUMNS 0 0 0 0 0 0 0 0 0 0 0 0 5 00000000 01 0 0 0 0 0 0 0 0 LI 0000000 00 590000000 00 0 0 0 0 0 0 1 00000000 00 0 0 0 0 0 0 0 0 0 0 2 00000000 00 MATRIX 5 ROWS 5 CULUMNS 1 0000000D 00 0 0 0 0 0 0 0 0 0 0 1 28402540 00 0 0 0 0 0 0 0 0 0 0 1 64872130 00 8 24360690 01 0 0 0 0 0 0 0 0 1 6487213D 00 0 0 0 0 0 0 0 0 0 0 2718281950 00 INT MATRIX 5 ROMS 5 COLUMNS 5 00000020 01 0 0 0 0 0 0 0 0 0 0 5 68050860 01 0 0 0 0 0 0 0 0 0 0 6 48 721310 01 1 15639380 01 0 0 0 0 0 0 0 0 6 4872131 0 01 0 0 0 0 0 0 0 0 0 0 8 59140970 01 Output Continued Figure 2 Continued q TEST PROGRAM GENERATES TRANSIENT USING TRNSI TRANSIENT RESPONSE INDICATES CONTROL CHANGES NEXT 7 1 0000000D 00 0000000D 00 1 0339835D 00 4085903D 0 7 8171846D 01 2 19452840 00 6 86126800 00 ELEMENTS CONTAIN Y 0 0 VASP PROGRAM HX LAST 2 1 00000000 00 22705246D 00 0000080 01 4 6788297D 1 00000000 00 50000000 ELEMENTS CONTAIN U zJR KX 1 00000000 00 00000000 00 1 40859030 01 7 50000080 01 2 40830520 00 DO 6 01398680 00 2 1945284D 00 2 50000020 00 6 7846557D 2 50000020 00 1 2798642D 01 3 5000000D 00 3 0000000D 00 1 21859130 01 8 5427698D 00 4 2500002D 00 1 2404001D TIME FIRST 5 ELEMENTS CONTAIN X 0 0 1 0000000D 00 1 00000000 00 1 00000000 00 0 0 0 50 7 5000008D 01 2 408
22. M _ 0 0 _ _ MATRIX ROW I COLUMNS eL LLL Ln 1 00000000 00 2 0000000D 00 34 0000000D 00 TIME RESPONSE MV DZ E RE XI 1 2 XT 3 YT 1 YT 2 YT 3 1 2 0 01 0 10201000 Ol 0 2040403D 01 0 3091364D 01 0 6151867D 01 0 2040403D 01 0 02 0 10404030 01 0 2081622D 01 0 3185510D 01 0 6307534D 01 0 20816220 01 0 03 x 0 1060909D O1 0 2123673D Ol 0 32825230 Ol 0 6467105 Ol 0 2123673 01 0404 Q 1081622D 01 0 2166574D Ol 0 3382491D 01 0 6630686D OL 0 2166574D 01 __ 0 05 0611025429 01 0 22103420 01 0 34855030 01 0 67983870 O1 0 22103420 01 000 0 06 0911236130 01 0 22549940 01 0 3591652D 01 69703190 OL 0 2254994D 01 0 07 0 11450160 01 0 23005480 01 0 3701034D 0 7146598D 01 0 2300548D Ol _ c Output Figure 1 Concluded The user s main program This program will be discussed statement by statement using the line numbers on figure 1 as a reference Lines 1 and 2 These two statements allocate the necessary storage for the variables to be used and define them as double precision Also the dimension arrays NF NG etc are allocated storage The dimensionality of F G etc could have been included in the double precision statement instead of the dimension statement and they could have been dimensioned as F 9 instead of F 3 3 The W array has been set up for du
23. NH NG Dummy Arguments Matrix DUM Constant KDUM Note KDUM contains the size of the DUM array which must be at least 7 min NR 1 NR 2 EXAMPLE USES OF VASP PROGRAM The examples given demonstrate directly the use of the principal subroutines EAT ETPHI AUG RICAT SAMPL DECGEN and PSEUDO In addition they exercise all of the subroutines except TRCE They can be used to indicate whether the programs are working properly They do not however provide an exhaustive test of the VASP program The first example discusses the user s main program in great detail to explain some of the system features The remainder of the examples simply state the problem and present the main program listing the data listings and the results Example 1 Transient Response A set of equations for a linear plant can be written as x t Fx t Gu t x 0 x y t Hx t where y are respectively the state control and observation vectors The system distribu tion and observation matrices F G and H respectively It is known that t x t ty eF t r G r u r dr O is the solution for x t If G and u are constant then 28 L x t eFtx_ LICL dr Gu By letting s t 7 the integral becomes O t f eFSd s I ds O Thus the solution to the system equation be written x t Bx CGu y t Hx t where eFt and t C f eF3 ds 0 It is des
24. NR DDI DABS B L IF FMX GE FDI GO TO 23 M L MX FDI QR1 L 0 12 692 SET TOLERANCE FOR ANDRA LIMIT OF SIZE OF DIAGONAL TOLERANCE OF ZERO IN ANDRA CALL DPR 1 DABS DP1 8 M C ASK FOR ALL ROWS DONE IN 1 CALL F 23 L I JP 1 I7 JP2 FIRST TIME INITIALIZATION FOR ANDRA JP 2 17 IF QIT 02 GO TO 561 JP 4 NR M 17 SOLD EFZ HAVE FINISHED PRELIM PART C INITIALIZATION FOR ANDRA DIAGONALIZATION NOW COMPLETED CALL ANDRA TO DIAGONALIZE SYMMETRIC MATRIX CALL ANDRA REDUCES ROWS BY MODIFIED GAUSS METHOD USING SORT PIVOT 30 C ONT TINUE TFIGOTIR U7 GU _ SAVE OLD VALUES CASE PIVOT IS REJECTED UNDER OPT DO 31 L 1 QNT J KRC L KRC2 L D J 8 L 31 D K C L 32 CALL ANDRA B C DPR JP JP 1 IR JP 3 C CHECK COMPLETION IS MATRIX ALL DANE IS MATRIX INVERTIBLE IFFOITR 17 ORe IR CEQ NR QIT 0 IZ GO TO 700 CHECK IF ITERATING WITH RHO TEST QR NOT C ALIT IF NO ITFRATION OR NO NEW PIVOT FOUND OMIT ITERATION CALCS IF NO NEW PIVOT DECRFASE TOLERANCE IF JP 5 F0 M GO TO 220 C COMPUTE RHO FOR ESTIMATING THRESHHOLD TO STOP SS IS RHO SS 6BDNRM NR C EE KRV BDNRM NR C EE
25. QS L FORCE PIVOT TO aACODEM FOR ONE B M DCC C SIGNAL NO LONGER FIRST TIME CALLED JP 2 1 s UPDATE EFFECTIVE RANK FOUND 1 DPR 2 DMX JP 5 NOW TEST IS THIS AN ITERATION TO DO ONLY 1 ROW AT A IF OKR OS GO TO 480 IF JP 1 IZ GO TO 490 C AT THIS POINT EITHER STOP WITH ONE ROW OR TRY NEXT C HERE GET READY TO RETURN RANK PARAMETER 4 480 JP 3 QKR RETURN IF ENOUGH TRIES TO D ALL ROWS PLUS 1 MORE QUIT 490 IF QITR ORI GO TO 480 QITR OITR GO TO 200 691 CALL LNCNT 1 WRITE 6 1691 QS QNT RM RETURN 692 CALL LNCNT 1 WRITE 65 1692 1692 FORMAT ERROR IN ANDRA FINDS NO P IVOTS CHECK FOR DIAGONAL ALLOWING NO PIVOTS 505 IF JP 2 EQ IZ ORI GO 692 GO TO 80 END TB 8v SUBROUTINE DECOM DCM KP SUBR DECDM FINDS 3 MATRICES FROM WHICH N GE COOP PL C INTO KRONECKER PRODUCT POSSIBLY USING A SEPARATE PSEU CALL Cee SURR DECOM INFORMATION SYSTEMS CU MAY 1969 J W ANDREWS IMPLICIT RFAL 8 D INTEGER 2 0 DOUBLE PRECISION A B C E Dy PIV COMMON MAX MA DIMENSION A1400 8 400 C 400 E 400 12000 2 JL 400 DCM 3 KP 4 NV 2 ND 2 EQUIVALENCE NV I PIV ND 1 DFZ0 I17 02 ND 2 J 1 ND 1 81 DATA 70 0 00 C SET PARAMETERS TO CALL PSEU FIND T TRANSFORMATION IN C IM OS KP 3 QNT 05
26. SUMMARY VASP is a Variable dimension Fortran version of the Automatic Synthesis Program a computer implementation of the Kalman filtering and control theory It consists of 31 subprograms for analy sis synthesis and optimization of complicated high order time variant problems associated with modern control theory These subprograms include operations of matrix algebra computation of the exponential of a matrix and its convolution integral solution of the matrix Ricatti equation and computation of dynamical response of a high order system Since VASP is programmed in Fortran the user has at his disposal not only VASP subprograms but all Fortran built in functions and any other programs written the Fortran language Allthe storage allocation is controlled by the user so the largest system that the program can handle is limited only by the size of the computer the complexity of the problem and the ingenuity of the user No accuracy was lost in converting the original machine language program to Fortran The principal part of this report contains a VASP dictionary and some example problems The dictionary contains a description of each subroutine and instructions on its use The example prob lems give the user a better perspective on the use of VASP for solving problems in modern control theory These example problems include dynamical response optimal control gain solution of the sampled data matrix Ricatti equation matrix deco
27. T 0015 10 CALL EAT F NF TT B NB C NC W 18 0016 CALL MULT B NB XO NXO V1 NV1 50 gpl 0018 CALL MULT AT NAL1 US NU V2 NV2 0019 CALL ADD VI NV1 V2 NV2 XT 0020 GM L MULTH HS NH XEN TT NY Fp 0021 CALL LNCNT 2 0022 WRITE 6 102 TT UXTUID 121 3 YTCI 4 121 2 BER a Ta PTaTRRBELTARL 0024 IF TT GT TFINAL STOP 0025 GO TO 10 Eno w v _ _ _ User s main program Figure 1 Example 1 GENERATES TRANSIENT USING EAT TEST PROGRAM 1 01 2 00 9 el oom e en e ojo Clo CI e e elm em ern m Ooo eio OIO 9 9 e e ec 9 Lj o ojo ojo Ijo gt HOO Xam b Data deck Figure 1 Continued 31 TEST PROGRAM 1 GENERATES TRANSIENT USING EAT VASP PROGRAM F MATRIX 3 ROWS 3 COLUMNS cc 0 0 2 00000002 00 0 0 0 0 0 0 3 0000000D 00 G MATRIX 3 ROWS 3 COLUMNS 1 0000000D 00 0 0 0 0 pup BODO OOD Oi ac M TR ai iu ME Lt er led IL es 0 0 0 0 1 0000000D 00 HL MATRIX 2 ROWS 3 COLUMNS 1 0000000D 00 1 0000000 00 1 00000000 0 0 1 0000000D 00 0 0 U MATRIX 73 ROWS 1 COLUMNS 1 00000000 00 _ 0 0 2
28. VASP PROGRAM The 92 characters are contained in the array TITLES which is in turn contained in the COMMON area LINES If N gt NLP the printer will automatically skip to the top of the next page and print the title line USAGE CALL LNCNT N Input Arguments Constant number of lines to be printed Output Arguments None 5 ADD DESCRIPTION This subroutine computes the matrix sum USAGE CALL ADD A NA B NB C NC Input Arguments Matrices A B Dimension arrays NA NB 1 Output Arguments Matrices C Dimension array Matrices share the same storage space or matrices share the same storage space 6 SUBT DESCRIPTION This subroutine computes the matrix difference USAGE CALL SUBT A NA B NB C NC Input Arguments Matrices Dimension Output Arguments Matrices C Dimension array NC REMARK Matrices A and C may share the same storage space or matrices B and C may share the same storage space 7 MULT DESCRIPTION This subroutine computes the matrix product C A B USAGE CALL MULT A NA B NB C NC Input Arguments Matrices A B Dimension arrays NA NB Output Arguments Matrix C Dimension II V 8 SCALE DESCRIPTION This subroutine multiplies every element of matrix by S and stores the resulting value in that is where S isa scalar USAGE CALL SCALE A NA B NB
29. lt F y Enn Oa END ll 601 SUBROUTINE JUXTC A NA B NB C NC DIMENSION A 1 B 1 5CC 1 NA 2 NB 2 NC 2 T DOUBLE PRECISION A B C COMMON IF NA 1 NE NB 1 GO TO 600 NC 1 NA 1 NC 2 NA 2 NB 2 L NA 1 2 NNC NC 1 NC 2 IF NA 1 LT 1 OR L LT 1 OR L GI MAXRC GO TO 600 IF NC 2 LT 1 OR NNC GT MAXRC GO TO 600 MS NA 1 XM 2 usu Orr a C 1 A T MBS NA 1 NB 2 J MS I 20 C J B T RETURN 600 CALL LNCNT 1 MRITE 6 1600 NA NB 1600 FORMAT DIMENSION ERROR IN JUXTC 2 6 5 216 CALL ASPERR RETURN END OTI SUBRDUTINE JUXTR A BP NB Cy 1 DIMENSION 1 1 1 2 2 2 DOURLE PRECISION COMMON X MAXRC IF NA 2 2 060 TO 600 2 2 NC 1 NA 1 NB 1 L NA 1 NA 2 NNC NC 1 2 IL 2 LT 1 OR NNC G1 MAXRC 2 MRA NA 1 MRB NB 1 MRC NC 1 DO 10 1 1 MCA DO 10 J 1 MRA KzJ MRA I 1 J MRC 1 1 10 C L A K DO 20 I 1 MCA DO 20 J 1 MRB X K J MRR 1 1 MR A J TMRC 1 20 C L R K RETURN 600 CALL LNCNT 1 WRITE 641600 600 FORMAT DIMENSION ERROR IN JUXTR CALL ASPERR RETURN 1 1 GO TO 600 GO TO 600 NAz 216 5X NB 21I6 SUBROUTIN
30. 0 0 1 0 0 3 3 3 0 1 0 0 0 190 1 0 0 0 0 0 1 0 2 2 1 0 1 0 1 0 2 0 b Data Figure 4 Example 4 L5 TEST PROGRAM FOR SAMPL CASE FROM ASP MANUAL P234 AND P244 VASP PRDGRAM PHI MATRIX 3 ROWS 3 COLUMNS c G0 0 0 0 0 0 0 0 0 0 0 2 0000000D 00 H MATRIX 2 ROWS 3 COLUMNS 0 0 2 0000000D 00 0 0 o 1 0000000 09 0 MATRIX 3 ROWS 3 COLUMNS FO 1 0000000D 00 1 0000000D 00 0 0 0 0 0 0 1 00000000 00 R MATRIX 2 ROWS 2 COLUMNS 1 0000000D 00 1 00000000 00 6 lt 21 _ 03 322 2 gt ___ ___ _____ STEP NUMBER O IN SAMPL K T MATRIX 3 ROWS 2 COLUMNS 4 28571430 01 1 428517140 01 0 0 y 1 42857140 01 7 14285710 01 SLEP NUMBERS TIIN P T MATRIX 3 ROWS 3 COLUMNS 3 1428571D 00 1 0000000D 00___ 24 8571429D201 h 2ML 1 00000000 00 1 00000000 00 0 0 2 85114290 01 0 0 3 5714286D 00 Output Figure 4 Continued STEP NUMBER 1 IN SAMPL _ MATRIX 23 ROMS 2 COLUMNS 4 1489362D 01 7 4468085D 02 0 0 0 0 2 6 95 745D 0 1 329 8 2D_ 00
31. 5 5 2 7755576D 17 1 00000000 00 4 1633363D 17 4 16333630 17 1 3877788D 17 4 16333630 17 1 00000000 00 0 0 5 55111510 17 4 163336359 17 0 0 90000009 006 Case 7 Figure 5 Continued 9 R MATRIX 4 ROWS 4 COLUMNS 1 5064996D 02 1 398 6860D 02 3 09032340 01 1 00423367 01 1 6179018D 00 1 5064996D 02 z Se 7 1972651D 01 6 1673338D 02 2 1264018 01 1 19726519 01 6 56133380 2 5 38436559 0 M1 MATRIX ROWS amp COLUMNS 3 09032340 01 1 6179018D 00 1 5064996D 02 1 0042336D 0 1 1 5064996D 02 1 39868010 02 5 Ty 7 19726510 01 6 1673628 02 20426401 80 01 1 19726510 01 6s 16T3628D 02 5 38415120501 4 ROWS 4 COLUMNS 5 55111510 17 RERR MATRIX 2 2204460D 16 1 73472350 18 4 1633363D 17 0 0 5 8972036D 08 8 3266727D 17 2 9047379D 07 1 38 88D 1 _ 1 5265567D 16 2 90271379p 07 3071633p ge MATRIX 4 ROWS 2 COLUMNS neg FE BAT Ch me Qa Dg hp 2 sh _ F nT 0 0 1 2719677D 00 1 17671260 01 1 18438510 02 4 671365399 03 5 658333 G MATRIX 2 ROWS 4 COLUMNS _ _ M e a FFE A t TeHp6b Gik TYt53 P rY _ v gt M 7 2 42956120 01 1 271967
32. 5 53529140 02 1 13136860 01 5 5352914D 02Z 7 951442 30 02 4 TTERATIONS K T MATRIX 1 ROWS 3 COLUMNS 54622 75D 01 2 6920418D 01 5 3036153 8 01 Output Continued Figure 3 Continued Lv P T MATRIX 3 ROWS 3 CULUMNS 22 5 16404040 01 2 11641710 02 1 1313684D 01 2 17641710 02 5 64703970 02 5 53529500 02 1 1313684D 01 5 53529500 02 7 95144720 02 A Re a A a s rd um A MuR 4 oov FSTR MATRIX 3 ROWS 3 COLUMNS 2 7670000D 01 1 00000000 00 3 72000000 02 7087200D 01 14 7850000D 01 1 2198300D 01 5 03333380 00 2 46259190 00 1 02075110 Gu RANSITION MNT 2 MATRIX 3 ROWS 3 COLUMNS 9 9640412D 01 9 96518910 03 9 4141918 04 1 6739033 01 9 9592474D 01 1 15736830 0 1 4 97750550 02 2 31282830 02 9 01578260 01 s x Output Continued Figure 3 Continued 8v TEST PROGRAM 2A TIRE TT 0 0 0 01 0 02 0 03 0 04 0 05 0 06 0 07 0 08 0 09 0 10 0 11 0 12 0 13 0 14 0 15 0 16 0 17 18 0 19 0 20 0 21 0 22 0 23 0 24 0 25 0 26 0 27 0 28 0 29 0 30 0 31 0 32 0 33 0 34 0 35 0 36 04 37 0 38 0 39 GENERATES TRANSIENT USING XT 1 0 1000000D 0 9964041D 0 9911999p 0 9844623 0 9762668 0 9666891D 0 9558051 0 9436909n 0 9304217 0 9160726D 0 90071
33. 59 DSUM 8 I XC M DSUM L K 1 NR LI EE L DSUM 60 CONTINUE C NOW RE ENTER MAIN SEQUENCE WITH PS INV IN EE 808 C GO FIX UP B PSUED0D INVERSE PRESUMABLY HAVE NIAGONALIZED HAVE DIAGONALIZED WITH ALL 1S IN UPPER LEFT C HERE WE HAVE FINISHED DIAGONALIZ WANT PSUEDO INV IN B C NEED TO SAVE DIAGONALIZED FOR USE BY DECOM CALL ODCM NEG FLAG DO 871 I 1 0 C A WAS SYMMETRIC JUST MOVE TO B RETURN FROM PSEUP ENTRY 871 EE T 60 877 80 CONTINUE C NOW FORM T TRP T APPRDX B SHRP PSUEDINV IN MATRX CALL T TR M NR C EE 808 IF QPS EQ 07 GO TO 870 IF QTP 819 819 818 C TRANS E TRP SHRP NRA LE NCA 818 DO 8181 I 1 0 DO 8181 J 1 NR DSUM DFZO QN J 1 NR DO 8182 K 1 L I 1 QR M ON K 8182 DSUM DSUM A L EE M L J 1 NCA I 8181 RB L DSUM 877 819 DO 8191 I 1 DO 8191 J 1 OR DSUM DFZO DO 8192 K 1 NR L K 1 QR J M K 1 NCA I 8192 DSUM DSUM A L EE M NOTE NCA IS USED BECAUSE A SHARP IS TRANSPOSED IN DIMENSTONS J 1 NCA I HERE EE A TRANS TRP A SHRP A TRANS NCA 819 8 L DSUM C HERE GET READY TO RETURN 877 C ONT TNUE C MOVE RANK TO RETURN PARAMETER IP 2 IR DEP 3 DPR 2 ABOVE RETURN PIVOT FROM ANDRA
34. 64 w COMPUTING DONE wwx 65 rewrt 150 LL or ee EE 66 0000100 CALL MULT A NA NX Y NY 67 0000150 call add x nx y ny w nw 68 0000250 call prnt rnt w nw w w 9 0000350 call prnt x nx x 1 70 0000450 71 MNPG NOW COMPILING wwx Figure 9 Example of conversational VASP _ Continued cen MM ee T a lt U M S M a MP L es ES 72 0000350 E www NOT FOUND WHERE REQUIRED 73 0000350 CALL PRNT X NX X 1 74 350 call prnt 1 _ s ___ wee 75 P 76 MODIFICATIONS 77 n LLL 78 DATA 79 80 VV MATRIX 3 ROWS E 81 1 5153413D 02 82 1 5705153D 02 83 1 6 768930 02 84 X MATRIX 3 ROWS 1 COLUMNS 85 1 0000000D 00 86 0 0 8 0 0 88 rrswCOMPUTING 89 Logoff 90 9 729 CPU SECONDS 05 11 71 AT 11 55 FST0 R3984 Figure 9 Example of conversational VASP Concluded being entered columnwise Do forget to input the matrix dimensions such as the example Data entry is ended with Execution of the program continues lines 18 through 27 display the requested output and line 28 indicates completion At this point line 29 the computer gives an underscore and the user may
35. ADD PHI NP2 W NP2 NPT W NP2 CALL 2 1 NPT SUMN 0 SUM 0 DO 202 IL 1 N ILP IL NP3 MIL IL 1 XM IL SUM SUM DARS PT NIL SUMN SUMN TDABS PT NIL W AL SUMN SUM NO 201 ILz1 N NIL IL 1 N IL ILP IL NP3 PT NIL 204 NO 104 M 2 N N1 M 104 L 1 N 13 2 MX M 1 PT L1 PT L1 PT L2 2 PT L2 PT L1 IF AL 00001 203 203 104 104 CONTINUE NTOT I NPRINT GO TO 305 203 NTOT NTOT J 305 CALL MULT 103 GO 403 400 401 402 NZ 400 CALL LNCNT 1 WRITE 6 50 NTOT 50 FORMAT 10X 14 14H ITERATIONS CALL PRNT 1 GN TN 403 401 CALL LNCNT 1 WRITE 6 50 NTOT CALL PRNT K NKy K T 1 GO TO 403 402 CALL LNCNT 1 WRITE 6 50 NTOT CALL PRNT K NK K T 1 CALL PT NPT P T 1 00001 210 210 403 403 CONTINUE C REARRANGE PHI MATRIX 210 CALL EQUATE PHI 1 NM W 1 NM 303 Iz1 N zTWOXN I 1 1 NTH3 NTHI N 1 1 1 1 NW2 NW1 N N CALL EQUATE W NW1 NN PHI NTH1 303 CALL EQUATE W NW2 NN PHT NTH3 NN CALL EQUATE PHINWW NM W NWW NM DO 304 Iz 9 NTH2 TWOXNX N I 1 1 NTH4 NTH2 N N TWOXN I 1 41 NW4 NMWM3 N N CALL EQUATE W NW3 NN PHT NTH2 NN
36. ANDRA CALLED BY PSEU J W ANDREWS INF SYSTEMS APRIL 1969 IMPLICIT REAL 8 D INTEGER 2 Q DOUBLE PRECISION DIMENSION B 400 T 400 DPR 2 JP 5 EQUIVALENCE 001 01 100 DMX FMX DRS TITS EQUIVALENCE DFZO FZRO IZ DATA ICC DFZO Z40000000 0 D0 C DPR1 IS MAGNITUDE THAT IS CONSIDERED ZRO PIVOT MUST BE NO SMALER DPR 2 IS TO RETURN FINAL PIVOT SO THAT USER MAY TEST SMALLNESS CC ANDRA CAN BE USED ALL BY ITSELF TO GET INV RANK OF POS SYMM TE THAT ORT HAS BE TAKEN PIVOTS ALONG THE DIAGONAL NOTE I AM DELIBERATELY ALLOWING SOME PARAMETERS TO CHANGE ON SUBSE C QUENT CALL DPR 1 CHANGES PIVOT SIZE A ROUGH TOLERANCE FOR ZRO PR C TEST IS THIS AN INTIALIZATION CALL TF JP 2 2 1 2 C INTIALIZE FORM IDENTITY MATRIX 1 QS JP 4 QNT QS QS IF QS 1 ONT GT 6400 GO TO 691 DO 18 I 1 ONT 18 ICI DFZO QR1 05 1810 I 1 QS T L 1 D0 1810 L 981 L DPR 2 DF7D C SET RANK ZRO TRIAL PIVOT VALUE TO ZERO QKR I7 C SET PIVOT CHOICE ITERATION O ALLOWANCE OF ROWS 1 ITER ucc UI a a 2 CONTINUE 200 CUNT INUE 2 f MAX DTA AND DI A MV 2 AR Srl FMX FZRO IF ALL DIAG ELEMENTS TESTED YET QR 30 IF I EQ 05 GO TO 40 I I L OR L DDI B L GET CURRENT DIAG ELEMENT FOR INTEGER SINGLE PREC TEST UPDATE L TQ GET DIAG ELEMENT TEST
37. C RE DO TO GET PS INV THAT MOVES ALL 15 OF DIAGONAL TO UPPER LEFT DIAG COMPUTE U MATRX 5 IN ASP FOR TRANSFORMING ORIG IN SINGULAR CASE L s DO 527 I 1 DO 525 J 1 NR J 1 NR I IF B L 521 522 521 522 C K C K C L DFZO 525 521 C K 0L70 III 14 il 861 C L 1 D0 525 CONTINUE 527 QR1 L C SAVE RANK SO FAR SHOULD BE SAME SIZE AFTER RE INVERSION QR2 IR DO 54 I 1 NR DO 54 1 NR DSUM DFZO K NR DO 53 J NR M I 1 NR J ON J 53 DSUM C M D L DSUM L QN I EE L DSUM 54 CONTINUE DO 56 I l NR DO 56 K 1 NR DSUM DFZO 1 DD 55 J J 1 NR I QN J DSUM EE L C M DSUM QN I B L DSUM 56 CONTINUE SET UP FOR SECONDARY ANDRA CALL NO ITERATION JP4 NR GO FIND LARGEST DIAG ELEMENT AGAIN gt GO TO 200 561 JP 3 17 CALL ANDRA B EE DPR JP IR JP 3 TEST FOR A CHANGE IN RANK eee ERROR IF QR2 IR 693 568 694 568 CALL T TRM NR D TRANSFORM C SHARP IN D BS U D U DO 58 I zl NR 6CI 00 58 DSUM DFZO K 1 NR DO 57 J 1 NR M J 1 NR J L ON J n 57 DSUM C M D L DSUM L I B L DSUM 58 CONT INUE DO 60 I 1 NR DO 60 K 1 NR DSUM DFZO DO 59 J 1 NR QN 3 1 M K L QN I
38. KDUM Programs 25 through 31 are called internally and need not be used by the programmer They are described in appendix A CALL READI A NA NZ NAM CALL ASPERR BLKDATA nonexecutable CALL PSEU A B C EE DEP IP D CALL PSEUP A B C EE DEP IP D FUNCTION BDNRM NR CT EE D KRV CALL CALL ANDRA B T DPR JP CALL DECOM A B C E JL DCM KP D The remainder of the items are procedures to facilitate the use of VASP on the Ames TSS VASP inputdsname outputdsname CHNGIN inputdsname CHNGOUT outputdsname CMPL CLRVASP CONVASP matrixsize A name B name C name W name X name Y7name Z name RECMPT REWRT n TABLE 2 STORAGE REQUIREMENTS AND EXTERNAL REFERENCES Storage Subroutines decimal bytes External references READ RDTITL PRNT LNCNT ADD SUBT MULT SCALE TRANP 10 INV 11 NORM 12 UNITY 13 TRACE 14 EQUATE 15 JUXTC 16 JUXTR 17 EAT 8 ETPHI 19 AUG 20 RICAT SAMPL 00 0 22 TRNSI 23 PSEUDO DECGEN 25 READI 26 ASPERR 27 BLKDATA 28 PSEU 29 BDNRM 30 ANDRA 31 DECOM TOTAL COMMON MAX COMMON LINES COMMON FORM READi PRNT LNCNT Z C C Yf amp W x ADD MULT SCALE NORM UNITY EQUATE ADD MULT SCALE NORM UNITY EQU
39. NEPR FMT 77 107016 7 gt 77 7 7 U DATA FMT2 3X 1P 7016 7 DATA TITLE 19 5 PROGRAM END 11 BR As DEP JP URE B V A Re ANDREWS 4 V Hs SUBR TO COMPUTE PSEUDO INVERSE CF GENERAL MATRIX RETURN FINAL PIVOT NOTE IMPLIT STATEMENTS MUST BE FIRST CAN BF REPLACED TYPE IMPLICIT REAL 8 D INTEGER 2 0 DOUBLE PRECISION IS ONLY THING ESSENTIAL INTEGER 2 M DOUBLE PRECISION D DIMENSION A 400 B 400 400 400 D 2000 1 KRV 4 2 DEP 3 DPR 2 4 JP 5 DATA ICC DFZO Z40000000 0 00 EQUIVALENCE 100 DMX FMX EQUIVALENCE DDT DSUM FZRO 12 07 OLL QR1 KRV 1 5KRC 5 KR KR KRV 3 4 KRC3 KRV 4 5 1 60 TO 1000 TRY PS D OPS IZ IP 4 IP 3 1000 CONTINUE DPl DEP 1 EF2 SNGL DEP 2 SET DEFAULT VALUES OF TOLERANCES IF DEP 1 DFZ0 2 D 6 IF EF2 FZRO EF2 1 0 NCA IP 4 C NUMBER OF ROWS OF ORIGNAL INPUT MATRIX QR IP 3 C SET SW FOR 0 DO ALL STEPS 0 THEN WANT RANK ONLY ONT QR NCA C TEST DIMENSIONS INPUT FOR REASONABLENESS LI 2 ORe ONT GT MAXRC UR QR LT 1 GO TO 691 IF DIMENSIONS ABSURD PSEU ERR EXIT 1 QDCM IP 1 QITR Q
40. NEXR 2 2 9 NGU2 2 9ND 2 NR1 2 __ _ ___ 0002 DOUBLE PRECISION ZL 100 T 100 EXR 100 R 100 G 100 H 100 1 DUM 700 R1 100 RE 100 EXC 100 f 0G002 A C OMMUN AM AX MARR E oa 0004 MAXRC 100 0005 e006 Me 0007 20 CALL RDTITL 0008 CALL READ 4 ZL NZL T NT EXR NEXR4 EXC NEXC4 T NT 490 9 AL E OR TH H 4 NT4 BUMKBUMJS IR DT MR DNE 0010 MzNT 1 0011 M2 M M 1 OO12 U 8 5C 7 gt _______ _ 0013 M3 M2 M5 0014 M4 M3 M5 OO DNS 0016 M6 M5 MS 0017 M7 M6 MS 0019 CALL MULT TNT ZL NZL DUM M3 ND 0020 CALL MULT DUM M3 ND DUM M2 ND ZL NZL 94923 FF NZtHOHHRNESNEXC OHHH88 T0 36 d ecco puede 0022 CALL JUXTC ZL NZL EXC NEXC 0023 60 TO 50 OP T Ff N7 2 J NF7NEXRT2YYG6E Tfr wy 0 0025 CALL JUXTR ZL 0026 60 50
41. S Input Arguments Matrix A Dimension array NA Scalar S Note If S isa constant it must be written as a double precision constant i e 2 00 0 D0 etc Output Arguments Matrix B Dimension array NB Note and can be the same matrix 9 TRANP DESCRIPTION This subroutine rearranges the elements of matrix A so that Bij Aii USAGE CALL TRANP A NA B NB Input Arguments Matrix A Dimension array NA Output Arguments Matrix B Dimension array NB 10 INV DESCRIPTION This subroutine computes the matrix inverse of A and stores this inverse in A that is Note that after the inversion is performed the values stored in the original matrix destroyed and replaced by the corresponding elements of its inverse USAGE CALL INV A NA DET DUM Input Arguments Matrix A Dimension array NA Output Arguments Matrix A the inverse of the original A Scalar DET the determinant of A Dummy Argument Matrix DUM work vector of length 2 NA 1 This subroutine is a slightly modified copy of the inverse routine given in the IBM scientific subroutine package 11 NORM DESCRIPTION This subroutine computes the norm of the matrix A as follows I Al min es 2 Ajj USAGE CALL NORM A NA ANORM Input Arguments Matrix A Dimension array NA Output Arguments Scalar ANORM 12 UNITY DESCRIPTION This subroutine computes the unit matrix USAGE
42. Such a technique works in a great many different Fortrans If a row is found to be marked by a pivot code it is skipped in steps 1 2 and 3 above The pivot position is forced to be exactly 1 before step 2 is applied The pivot code is actually tested for as an integer The pivot size is tested for in single precision These modifications are for speed A count is kept of the number of pivot searches If this count is one greater than the num ber of rows the routine always stops searching for pivots The result if B has maximum rank is matrix T such that TIT inverse of B The input B consists of 0 s everywhere except the diagonal which holds pivot codes 86 Error Messages From ANDRA Message Explanation Dimension error The total number of matrix elements was too large or too small The parameter JP 4 cannot be less than one nor more than MAXRC Finds no pivots ANDRA could find not a single pivot in its very first search of diagonal 31 DECOM SUMMARY Fortran IV subroutine DECOM generates four double precision output matrices from the symmetric non negative definite input matrix B One output is a matrix S such that if E isa unity matrix of rank r then B SEE s This matrix is obtained as the inverse of a matrix T by calling subroutine INV T comes from subroutine PSEU It is defined by E a diagonal matrix with elements 0 or 1 E is also returned along with a permutation matrix P such that PEP
43. X to obtain Y Line 21 Line 21 tells the counter we are going to print 2 lines If this will not fit on the present page LNCNT will advance to the next page print the title as on the first line of the first page of output and increment the line counter to allow for the paging and the two lines about to be printed Line 22 Line 22 prints the variables XT and YT skipping a line between each print line as required by the 1HO in FORMAT 102 Note that YT 3 is not printed Example 2 Transient Response Using TRNSI This example uses the same equations as Example 1 except that u is piecewise constant that is u t JR Kx t t itj lt lt It where i is a non negative integer and J R K are constant matrices The first term JR represents a forcing function and the second Kx is a feedback term See ASP manual p 121 for detailed explanation It is desired to generate the transient response of such a system in response to a given initial condition and a time varying control In particular given 0 I 1 0 0 1 0 1 M O O O O 34 1 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 1 1 01 O t 0 5 sec t 2 sec tg 3 5 t9 sec The system is monitored at intervals t while the control u t is changed only at sampling intervals t t must be a positive integral multiple of t Specifically the control
44. ZL MATRIX 2 ROWS 2 COLUMNS 1 00000005 00 9 9 0 0 2 0000000D 00 ui MATRIX RS COLUNNS _ ______ 1 99569850 01 6 00573110 01 6 00573110 01 7 9956985 0 01 MATRIX 1 ROWS 1 COLUMNS 0 0 C XC MATRIX 2 ROWS 1 COLUMNS 2 0000000D 00 3 9009009905 909 MATRIX 2 ROWS 2 COLUMNS n H nM n 5 6511539D 01 8 2501188D 01 MATRIX 2 ROWS 2 COLUMNS 1 00000000 00 0 0 0 0 1 0000000D 00 gt __ gt gt gt 3 _ R MATRIX 2 ROWS 3 COLUMNS 131935540 00 4 66226910 01 2 00000000 00 6622691 01 1 6 8064460 00 3 900000009 0900 R1 MATRIX 2 ROWS 3 COLUMNS 1 3493554D 00 4466226919 01 240000000B 99 n 4 6622691D 01 1 68064460 00 3 0000000D 0 Case 4 Figure 5 Continued L9 RERR MATE LX 2 3 EOL IM NS 2 22044600 16 2 49800180 16 4 44089210 16 9 7144515 17 222204460D 16 424408921D 16 H MATRIX 2 ROWS 2 COLUMNS 2 0493573D 00 1 3259717D 00 FB G MATRIX 2 ROWS 3 COLUMNS 7 20719060 01 5 140859658 01 4 165579 6 01 1 3435357D 01 4 84314830 01 8 64516200 01 RANK MATR
45. and space activities of the United States shall be conducted so as to contribute to the expansion of human knowl edge of phenomena in the atmosphere and space The Administration shall provide for the widest practicable and appropriate dissemination of information concerning its activities and the results thereof NATIONAL AERONAUTICS AND SPACE ACT 1928 A 2 NASA SCIENTIFIC AND _ _ PUBLICATION S RS TECHNICAL REPORTS Scientific and U TECHNICAL TRANSLATIONS Information tecHnical information considered impogtant published i in a foreign language considered c mplete and a lasting contribution oxisting sto merit NASA distribution in English knowledge 3 21 SPECIAL PUBLICATIONS Information TECHNICAL NOTES Informatica less broad ie derived from or of value to NASA activities in scope but nevertheless of importance asa 22 8 Publications include conference proceedings contribution to existing knowledge 107 52 27 monographs data compilations handbooks 2 sourcebooks and special bibliographies E TECHNICAL MEMORANDUMS Information receiving limited distribution TECHNOLOGY UTILIZATION because of preliminary data security PUBLICATIONS Information on technology tion or Other reasons R used by NASA that may be of particular See A interest in commercial and other non aerospace CONTRACTOR REPORTS Scientific and applications Publica
46. be symmetric c R must be invertible The Fortran symbol for is RI USAGE CALL AUG F NF G NG RI NRI H NH Q NQ C NC Z NZ ID Input Arguments Matrices F G RI H Q Dimension arrays NF NG NRI NH NQ Control constant II II 1 General case II 1 Special case H is not used in AUG Output Arguments Matrices CZ Dimension arrays NC NZ 20 RICAT DESCRIPTION This subroutine computes P t and K t by the following equations P t 7 021 0 0 0 P t K t CP t See ASP manual page 9 for reference MOTIVATION A standard control problem will be used to illustrate how this matrix Riccati equation arises Given the system equation x Fx Gu the output equation 21 and the performance index T J f x H QHx u Ru dt x T H S T Hx T 0 where Q R S symmetric matrices and is invertible We wish to find a control law which minimizes the performance index J Introducing the auxiliary variable A t into the system of equation we have the following Euler Lagrange equations 2 L GR G H QH F X X Z which have for a solution x t Xo Xo A t Xo optimal control law is u t R G X t Letting P t be linear transformation from the state variable x t to the auxiliary variable A t that 15 A t P t x t we obtain from the Euler Lagrange equation the following Riccati equation P F P PF PGR G P H OH where the
47. data set INPUT will also read the same data set Variables may also be set by Fortran arithmetic statements Output may be from the VASP subroutine PRNT or any Fortran WRITE or PRINT statement Two standard formats are available if desired for unlabeled output The program automatically dimensions 14 arrays to the desired size and the user may supply his own names to 7 of them Usage The use of conversational VASP is demonstrated by the accompanying figure fig 9 Lower case letters are input and upper case are the computer responses Detailed comments on the various statements follow To start the user calls VASP line 1 as for nonconversational usage II desired an output data set may be named Line 2 lists the DDNAME being used The next two lines lines 3 and 4 indicate where input and output are to reside The computer then gives an underscore after which the procedure is called The param eters of this procedure are first the total number of elements in a matrix followed by up to seven matrix names If the parameters are defaulted the system will select matrices with 9 ele ments and name the matrices A B C W X Y Z In addition 7 dummy matrices D1 through D7 are available for use In the figure all matrices are to be dimensioned 16 line 5 the second matrix is to be renamed F and the Z matrix is to be renamed FSTAR That is if you wish to rename a specific matrix put a dollar sign in front of
48. do anything he wishes In the example we are going to recompute with the same program using new data Accordingly the user asks for RECMPT line 29 The program is again executed and new data are asked for line 30 They are entered in lines 31 through 35 using a different style than in line 17 to show the flexibility available On completion of the data entry the results are printed in lines 36 through 44 At this point it is desired to rewrite the entire program so the user issues the command REWRT line 45 The system as at line 9 prompts the user with FORTRAN STATEMENTS and a line number lines 46 and 47 after which the user enters Fortran statements as desired In the example line 48 is entered incorrectly and then corrected lines 49 and 50 Following this a line 150 was inserted lines 51 and 52 Then a CMPL was issued line 53 to compile and execute the program New data were entered at line 56 and lines 57 through 59 are the output requested by the statement print 6 1 Note that all 16 elements of f are printed using one of the two FORMAT statements compiled into the program for convenience see lines 5900 and 6000 of VASPPROC fig 12 6 FORMAT 1X 1P6D13 6 13 FORMAT 1X 1P4D20 13 These statements request the output of a 6 decimal number or a 13 decimal number In the example we are printing a 6 decimal number The remainder of the output is then printed lines 60 through 63 Now it is desired to rewrit
49. for adequate accuracy that is all matrix and scalar variables except integers are assumed to be in double precision Universal Features The arguments in the subroutines are arranged in the following order Input arguments Output arguments Dummy arguments These are also arranged so that matrices occur before scalars An array of length two must be allocated by the user to store the dimensions of the matrix and this array must be included in the subroutine call statements For example the add subroutine is called by CALL ADD A NA B NB C NC and performs the matrix operation and NC arrays of length two which contain the dimensions of matrices B respectively In other words the numbers of rows columns of matrices and stored in NA NB and NC respectively Specifically the number of rows of is stored in NA 1 and the number of columns of A in NA 2 In general the dimension array associated with an input matrix contains input information to the subroutine while that associated with an output matrix contains output information The dic tionary shows the few cases where this rule is violated In the example above dimension arrays NA and NB are inputs since matrices A B are inputs and must be loaded before entering this sub routine On the other hand NC is an output since C is an output That is the values of NC 1 and NC 2 are computed in th
50. in this manual as example 6 shows errors determined by calculating AAT AAT LA on the order of 10715 or less The routine was also tested on the ill conditioned 7 X 7 matrix in the ASP manual NASA CR 475 p 150 The exact inverse is given on page 151 and the error obtained from the ASP program using the equivalent of the PSEUP entry p 152 is on the order of 1071 The error obtained using the VASP program and the PSEU entry was on the order of 10 or less Singular Case The routine forms a new inverse from the original symmetric matrix Since there are several steps more between the inverse and the original input A it is only natural that accuracy should fall off In many cases this inverse will give a diagonal matrix of 0 s and 175 when used as a left inverse of A or possibly as a right inverse The work of reinverting B requires no extra matrices it does destroy the usual values of C and D No iteration can be done in the stage after B is found to be singular It can be asked for in the starting stage Error exits are made if the rank changes during reinversion The smallest allowable pivot is redetermined Error Exits Messages The error exits are reasonably self explanatory Unless otherwise noted the errors cause a return from PSEU without completion of the calculations Subsequent calculations in other portions of the program are suspect Message Explanation Dimension error The total number of
51. matrix elements was too Jarge or too small Diagonal elements of matrix 0 Symmetric matrix B has no positive diagonal elements Check input A 74 Rank has decreased Rank has increased Rank greater than matrix size Timing Singular case Reinverting and the new rank is less than that of the earlier phase Singular case Reinverting and the new rank is greater than in the earlier phase Computation continues RANK returned from ANDRA 15 greater than maximalrank The ANDRA routine by itself is very fast The iteration mode is slower by large factor than the regular mode of subroutine PSEU The time estimates below in hundredths of a second are as used on the NASA Ames 360 50 High and low estimates are given because real time figures reflect an unknown percentage of time devoted to another CPU user Case PSEU 2 X 2 matrix PSEUP 4 X 4 matrix reinvert PSEU 7 X 7 no pivot rejection PSEU 7 X 7 rank 6 reinversion PSEU 7 X 7 iteration no tests PSEU 7 X 7 iteration one test PSEU 7 X 7 iteration some tests of pivots PSEU 7 X 7 iteration many pivot tests PSEU 7 X 7 iteration nearly all tests PSEU 4 X 2 reinversion PSEU 4 X 2 reinversion METHOD Summary of Method PSEU has two entry points The nonsymmetric entry forms AtA or whichever is High 2 14 42 103 53 182 253 501 607 Low I 10 30 62 30 118 170 286 324
52. subroutine PSEUP during diagonalization of input matrix A Integer control parameters just as for subroutine PSEUP Zero do not iterate in PSEUP One iterate in PSEUP Zero do all calculations Nonzero do rank only Changed to reflect the rank on output Should be set to zero or minus one before each call Note DECGEN uses and KP2 0 The row size of the matrix input The column size Note This parameter is forced negative as a signal if T cannot be inverted by INV D has five parts as does the dummy array in PSEUP Let these be denoted D D1 D2 D3 and D4 These five equal arrays must be included in the size of parameter D if iteration by PSEUP is selected If no iteration is used D2 D3 and D4 may be omitted D holds the inverse of output C DI holds the permutation matrix Note If rank only is computed D is computed but D is not A B C and DI are thus the only matrices with useful values returned rcm n METHOD The results from DECOM are an effective rank r matrices B and D which are used in further calculations to get a Kronecker decomposition or to see which variables are dependent and the permutation matrix P in This section describes the sequence to obtain the Kronecker decomposition in two different cases The goal is two matrices G and H DECOM does not produce these matrices they are produced either by DECGEN or by the user according to the following steps Le
53. the original name and then equate it to the desired name as in the example Fourteen arrays NA through ND7 used for dimension information are also defined and renamed to agree with the working matrices Lines 6 7 and 8 then define the matrices available Note that no l element variables are defined The user may define them in his program but they will not be available from one computation to the next The computer will then ask for FORTRAN STATEMENTS At this point a data set SOURCE MNPGS has been set up for editing and the necessary DIMENSION and other initial izing statements have been stored These statements are listed in figure 12 lines 4600 through 6000 The computer prompts the user with 100 and the user may enter any Fortran statements he wishes The full power of the text editor 15 available at this point In the example we have entered four statements lines 10 through 13 Note that we have defined a single variable t for use in the etphi statement The value of this variable will not be remembered by the system After completing the desired Fortran statements the user requests compilation by entering _CMPL line 14 The computer then indicates that compilation is proceeding line 15 and will give the usual error messages if the compile is unsuccessful After compilation the program is automatically executed and the first item in the execution is a request for data from the INPUT subroutine line 16 Data are entered fre
54. u t is updated by the equation u t JR Kx t 1 it S lt t t lt it It The x y vectors are computed at time intervals t and these vectors together with the time t and the control u for the subsequent time interval are printed out each time The problem terminates when the final time tr is reached The matrix T has elements LI 1 tg t in that order The user s main program to solve this problem is shown in figure 2 a the corresponding data deck is shown in figure 2 b and the results are presented in figure 2 c 190 DIMENSION F 54 5 NF 2 5G 552 6 2 291 2 1 1 NR 2 F K 25 200 15 NK 2 H 7 5 NH 2 X 5 NX 2 I 4 NT 2 DUMMY 100 300 DOUBLE PRECISION J ReKsH Xg T DUMMY 400 _ 1 6 FMT2 6 _ _ 500 CALL RDTITL 600 CALL READ 5 F NF G NG Jg NJ RS NRK NK 700 CALL READ 3 H NH IX NX T 5NT TRNS 7 F 6 6 T PUMMY 100 90N STOP 1000 END M KT Nn n n Y n m n nPT H s l t a User s main program Figure 2 Example 2 35 TEST PROGRAM GENERATES TRANSIENT USING TRNSI F 5 5 0 0 0 0 0 Oe 22 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 25 5 2 1 1 1 0 0 0 0 m J 2 1 0
55. x ALL CONVERSATIONAL VASP PROGRAMS CLEARED CLRVASP 0000700 DISPLAY x YOU MAY RESTART WITH CONVASPxssse xx ee ee CMPL 0000000 PROCDEF CMPL MPL 0000030 DEFAULT SYSINXSE CMPL 0000050 EDIT SOURCF MNPG CMPL 0000100 EXCERPT SOURCE VAS PMN 1600 1700 CMPL 0000150 END CMPL 0000200 DISPLAY MNPG NOW COMPILING w CMPL 0000220 DEFAULT LIMEN N CMPL 0000250 UNLOAD MNP6 Non CMPL 0000270 ERASE USERLIB MNPG CMPL 0000500 MNPG Y CMPL 0000400 LOAD MNPG CMPL 0000600 CALL MNPG CMPL 0000700 DEFAULT LIMEN M CMPL 0000800 DISPLAY COMPUTING M v M eg 12 m n X Figure 11 List of VASP procedures 091 ONVASP 0000000 PROCDEF CONVASP mmc SESS CONVASP 0000020 PARAM N A B C W X Y sz CONVASP 0000040 DDEF VASP1 VP VASP1 0PTION OBLIB 80 JOBLIBS SYSULIB gt CONVASP 0000110 DISPLAY MATRICES AVAILABLE ALL DIMENSIONED N ARE CONVASP 0000140 DISPLAY A B C W X Y Z FOR INPUT OR RTOS 2200002 0000170 DISPLAY
56. x DRS Reduce T Row Pivot Row B L DDM Reduce Pivot Column Symmetry B K DDM Update Next Element L Row Column Do Satisfied Yes 68 I Force Pivot to Unity B M 1 00 De 460 1 15 index QMP T Pivot Row Yes No Set Index Of Diagonal Element n Current Row DRS B K Test Part Of DRS For Integer Pivot Code 1CC K Index Element In Pivot Column To Eliminate DMM Take Negative Of Coeff Do 47 J 15 index Add Pivot Row Times OMM To Row 1 Ot B Add Pivot Row T Times DMM To Row Of T Do 460 Satisfied Set L Force Pivot Row To Zero Force Pivot To Code For One B M DCC No Longer First Cail JP2 1 Update Rank OKR 1 DPR2 JP5 QMR 7 Rank Maximal teration One Pivot Per Call JP1 0 480 Figure 8 Concluded Yes CT Of Tries gt Max Rank 1 Update Count of Tries 691 Error M055808 First Or Too Many Tries Call No Pivot 692 METHOD Mathematical The method is described in the ASP manual pages 137 to 139 The square matrix T is initialized to be the identity Step 1 The diagonal of B is scanned to find the largest pivot Pivots are only taken from the diagonal
57. 0 19 b Output Figure 6 Continued 9 PSEUDO TEST PROGRAM CASE 2 FROM ASP MANUAL PAGE 137 VASP PROGRAM A MATRIX 4 ROWS 4 COLUMNS 2 00000006 00 1 000000006 00 2 00000006 00 2 000000866 GO 1 0000000D 00 2 5000000D 01 8 00000000 6 0000000D 00 2 00000000 00 8 00000000 00 4 0000000D 00 0 0 2 0000000D 00 6 00000000 00 9 0 40000000P 00 APSE MATRIX 4 ROWS 4 COLUMNS _ 5 50970630 02 1 102609506 02 4 69110228 02 6 22831038 02 1 10260950 03 2 97370440 02 7 55120450 03 9 75642350 03 4 69110230 02 7 5512045D 03 4 2366935D 02 5 14551100 02 5 3283103D 02 9 7564235D 03 5 14551100 02 7 51110960 02 AASA MATRIX 4 ROWS 4 COLUMNS 1 0000000D 00 2 0000000D 00 2 00000009 00 100000000 00 2 5000000D 01 8 0000000D 00 6 00000000 00 2 00000000 00 8 00000000 00 4 0000000D 2 22044600 16 2 00000000 00 6 00000000 00 0 0 4 0000000D AERR MATRIX 4 ROWS 4 COLUMNS 4 4408921D 16 0 0 4 4408921 16 6 06133819 16 4 4408921D 16 3 55271370 15 2 22044600 16 4 4408921D 16 6 66133810 16 6 6613381D 16 6 66133810 16 2 22044600 16 2 2204460D 16 2 2204460D 16 0 0 4 4408921D 16 ASAA MATRIX 4 ROWS 4 COLUMNS 5450910630 02 1441 026095D 03 4 691102380 02 65 32831038 02 1 10260950 03 2 97370440 02 7 5512045D 03 9 7564235n 03 4 6911
58. 0 0 G 3 ROWS 0 0 2 6 6700000D 00 XO MATRIX 3 ROWS 1 00000000 00 0 0 JpDa0 2 L MATRIX 3 ROWS 1 00000000 00 0 0 0 0 1 00000000 00 0 0 0 0 0 MATRIX 3 RDNS 2 00000000 01 0 0 0 0 2 0000000D 01 0 0 0 0 R MATRIX 1 ROWS 1 0000000D 00 MATRIX 3 ROWS 0 0 0 0 0 0 0 0 0 0 0 0 ITERATIONS 3 COLUMNS 3 72000000 02 1 2198300D 01 6 6700000D 1 COLUMNS 1 CULUMNS 3 COLUMNS 0 0 0 0 1 0 0000000 00 3 COLUMNS oO O e d 00 1 COLUMNS 3 COLUMNS se S 4 9 Output Figure 3 Continued VASP PRAGRAM TEST PROGRAM 2A GENERATES TRANSIENT USING 1 1 CT VASP PROGRAM MATRIX 1 ROWS 3 COLUMNS 7 36037320 01 3 6359943 01 5 16454380 01 PLT MATRIX 8 ROWS 3 COLUMNS 6 67994270 01 2 1615347D 02 1 10350420 01 2 16153470 02 5 5929130D 02 5 45126590 02 1 1035042D 01 5 4512659D 02 7 14294420 02 2 ITERATIONS K T MATRIX 1 ROWS 3 COLUMNS 5463314D 01 3 691 988D 01 5 3032918D 01 P T MATRIX 3 ROWS 3 COLUMNS 6 1636913 01 2 1772154D 0 1 1313840D 01 2 1772154D 0 5 64663120 02 5 5349307D 02 1 1313840D 01 5 5349307D 02 7 9509622D 0 2 _ 5 K T MATRIX 1 ROWS 3 COLUMNS 7 5462283 01 3 6920394 01 5 3026120D 01 P T MATRIX 3 ROWS 3 COLUMNS 5 1640398D 01 2 12641915 202 1 13126860 01 2 17641910 02 5 64703680 02
59. 023D 02 7 55120450 03 4 23669350 02 5 14551100 02 5 32831030 02 9 7564235 03 5 14551100 02 17 5111096 02 AERR MATRIX 4 ROWS 4 COLUMNS 1 4 1451500 17 2 764271550 18 8 67361740 18 1 3977399 52 0 0 0 0 2 6020852D 18 8 67361740 19 8 6736174D 1 8 1 7347235D 18 5 20417040 18 1 12757030 17 1 3877788D 17 2 60208520 18 8 61361740 18 1 3877788D 17 b Output Continued Figure 6 Continued 89 PSEUDO TEST PROGRAM CASE 3 VASP PROGRAM C MATRIX 4 ROWS 2 COLUMNS lt lt _33 lt _ _ 23222 3 2232 2 2 2 32222 lt _ _ _ 3 50000000 02 0 0 2 5300000D 00 3 1000000D 01 2 _ T qc APSE MATRIX 2 ROWS 4 COLUMNS 20 a 2 u 13314210202 5 5012930D 03 3 0981008 E D5 ees m demone 0 0 2 55754170 01 2 8046074D 04 3 8798916D 06 cx AASA MATRIX 4 ROWS 2 COLUMNS 2 ec 0 0 0 0 1 08894560 17 3 91000000 00 2 53000000 00 3 10000000 01 MATRIX SREMS COLUMNS 0 0 0 0 1 08894560 17 4 4408921D 16 e ccu uu i E Lu uL eu cu LL uL C ES 2 Z204460D 16 0 0 0 0 3 13314710 02 5 50129305 03 3 95180810 01 0 0 2 55754170 01 2 8046074D 04 3 8798916 06 AERR MATRIX 2 ROWS 4 COLUMNS 0 0 0 0 0 0 0 0 2 71555769 1171 5 42101099 20 1 10335469 1
60. 0LTNESI0U 09400 JBLB VASP 09500LOAD BLKDTA u 09600 TNPUT 7s DDEF FTOSFO0UI INPUT DTSPLAY TNPUT FROM DATA SET 1NPUT gt 09700 IF INPUT s DISPLAY INPUT FROM TERMINAL 09800 IF OUTPUT 7z5 DDEF FTOGF001 0UTPUT DISPLAY OUTPUT PLACED IN DATA SET 0UTPUT 09900 ES OUTPUT ETS PLAT OUTPUT TO TERN rd Figure 12 List of data set VASPPROC Continued 691 10000 LIST 10100 END TODO ge ed 7 10500 VASP PROCEDURES NOW READY DO ABEND TO MAKE THEM AVAILABLE 10500 LOGOFF Figure 12 List of data set VASPPROC Concluded REFERENCE 1 Kalman R E and Engler T S A User s Manual for the Automatic Synthesis Program Program C NASA CR 475 1966 166 NASA Langley 1971 8 3882 NATIONAL AERONAUTICS AND SPACE ADMISTRATION WASHINGTON D C 20546 POSTAGE AND FEES PAID NATIONAL AERONAUTICS AND OFFICIAL BUSINESS FIRST CLASS MAIL SPACE ADMINISTRATION PENALTY FOR PRIVATE USE 300 0178 01 C2 UL 08 711008 50090355 720401 DEPT ATR FORCE AF WEAPONS LAB AFSC TECH LIBRARY WLOL ATTN E LOU BOWMAN CHIEF KIRTLAND NM 87117 sp If Undeliverable Section 15 POSTMASTER Postal Manual Do Nor gt __ 5 gt gt _ _ _ The aeronautical
61. 1 DO 20 I JP N GALL UNTETX CDUM y NT3 0013 II N I 1 4 0014 JJ N J 1 J T s n sN n 0016 JI N J 1 1 0017 CTH DCOS T IJ 06 18 SFH BS 0019 DUM 11 CTH 0020 DUM JJ CTH 6 3 ts s 0022 DUM JI 2 S TH 0023 CALL MULT DUM LDMI DUM NT DUM LDM2 NT 9024 ERtt EGUATE EDUMEH DM2 4NT 4 DUMCEDMETSNTI 0025 20 CONTINUE 0026 CALL TRANP DUM LDML NT T NT 88 33 fHd b Ut T CH 3NH rBUMHEBM NT BUMHH BM HHSNTH 0028 CALL PRNT T NT 4HT 1 0029 CALL DUM LDM2 4HT T 1 0030 0031 END b Subroutine ORTH Figure 5 Continued a LS TEST PROGRAM FOR DECGEN AND DECOM CASE 1 2X2 RANKI ZL 2 2 1 0 1 0 2 0 2 0 E T 2 2 E XR 1 1 EXC 1 1 TEST PROGRAM FOR DECGEN AND DECOM CASE 4 2X3 RANK2 P 1 220 T 2 2 ef EXR EXC 2 1 2 TEST PROGRAM FOR DECGEN AND DECOM CASE 7 4X4 RANK3 ZL 4 4 1 2 1 0 6 59 5 4 5 of EXR EXC I c Data Figure 5 Continued w 85 TEST PROGRAM AND DECOM CASE 1 2X2 RANK1 VASP PROGRAM Is nr ts
62. 10 NR NR 20 IR NR DO 30 I 1 IR LL I 1 IR 32 C gt lt 0000 60 DSUM KK K 1 IR DO 29 J 1 IR L Jj LL M J KK 29 DSUM DSUM 0 2 C ABOVE FORMAING T TRANSPOSE TIMES WHICH 15 APPROX OF SHARP KK IF OF 31 39 32 31 KK KD3 L D KK D L EE L DSUM 60 30 C 39 COMPUTE T TRANSPOSE T ONLY PROVIDES INVERSE SHARP 39 32 20 41 EE L DSUM GO TO 30 FE L D L KK KD3 L D KK DSUM CONTINUE 41 66 41 NV 1 IR NV 2 IR LET 1ST MATRIX EE Q 2D D KD3 CCC ROLES OF P AND O ARE GIVEN TO 2 MATRICES AT 31 SWITCHED AT 32 COMPUTE D K4 O P EEZD K4 Q EE EE D K3 Q P Q C FUNCTION IS NRM QxP Q Q NRM Q RESULT A SCALAR 9 CALL MULT D KD3 1 NV FE NV D KD44 1 NV NL IR IR CALL MULT D KD4 1 NV D KD3 1 NV NV KK I CC I D KK CALL NORM FE NV AN CALL NORM N KD341 NV RR N AS BSHRP APPROACH THAT FITTING 2 MOOR PENRSE AXIOM BDNRM FN FR RETURN 66 C66 IS A DUMMY REALLY WANT MATRX MULT ONLY C G0 TO 9 SIDE COMPUTATIONS J ANDREWS INF SYSTEMS MAY 1969 END vri SUBROUTINE JP C SUBROUTINE ANDRA DIAGONALTZES POS DEF SYMM J ANDREWS I S C SUBR
63. 133 06000 6 FORMAT 1X 1P6D13 6 06100 RETURN 06200 ENP 06300__END 06400 ERASE USERLIB VASPMN 065 99 F TN VAS PMNA 06600 XLIST VASPMN 06700 LOAD VASPMN Figure 12 List of data set VASPPROC Continued 06898 lt ee aa 06900 EDIT SOURCE MNPG 07000 5 1 LAST 07260 NSERT E i 07200__EXCERPT SOURCE VASPMN 100 1500 07300 DISPLAY FORTRAN STATEMENTS 02500 DEPAUL ee 07500 INSERT 100 07540_ PROCDEF RECMPT 07600 PROCDEF RECMPT 07700 DEFAULT LIMEN N 07800 CALE MNPG cm Ll v nIFT Y A U 07900 DEFAULT LIMEN W 08000 DISPLAY COMPUTING 0M 080 0 PROEREF REWRF ______ _ C _ FPVT 08080 5 1 LAST 08100 PROCDEF REWRT 08200 PARAM LINE 08500 BEFAULT EEMEN W F 08400 DEFAULT SYSINX E 08500 EDIT SOURCE MNPG 086000 EXCISE SEINE UNS E 08700 IF LINE 100 DISPLAY FORTRAN STATEMENTS 08800 IF LINE 7 100 LIST 100 LAST 08900 DpEFAULT Syst 6 v V F FT 09000 INSERT LINE 09040_ PROCDEF VASP 09080 EXCISE it NG sa 09100 PROCDEF VASP 09200 PARAM INPUT OUTPUT 9 A 8 5 0 7C W W X7X Y7Y 2 2
64. 300 FORMA TOP 310 JFN 1 GO TN 152 OO Chek ENGNT 2 WRITF 6 1000 NPHI NHNO NR 506 1000 FORMAT DIMENSION ERROR IN 5 7 735 2X4 4HNPHI 8X 2HNH 19X 2HNQ 9X 2HNR 9X 2HNP 5X 4HKDUM 3X 9HKD UM MIN 29X 56 3X 2 I 23X 145 5X T CALL ASPERR RETURN END IGER CHANGING X RTM SAMP vcl T 2 5 T1 T1 SUBROUTINE TRNST F NF G NG NJ R NR NK H NH IX NX T DUMMY KDUM DOUBLE PRECISION F G J K H X DUMMY DIMENSION F I NF V1 GC13 NGCI 4 JCL NJCI KC I NK 1 HC T 4 NHC3 1 1 NX 1 TUL R C1 NRCI s4 DUMMY I DIMENSION NNF 2 NNG 2 1ML 2 NNX 2 NY 2 DATA STAR x IF NF 2 NE NG 1 NJ 2 NE NR 1 ORe 2 1 ORe Y IC R NR x X e N e Re e IM L 2NF 2 1 GO TO 999 MAX NF 1 NF 2 NG 2 1 NH 1 9NK 1 IF LI GO TU 910 11 1 NSQ 11 NF 11 NX4 N EG IF KDUM LT GO TO 900 2 TRNST PROGRAM N2 Il NSO 2 NSO _ 2 2 999 __ N4 N3 NSQ L3 2 NF I1 NG 2 L4 L3 NJ I1 NR 2 L5 L4 NH 1 L L5 NJ T1 NR 2 1 T 1 11 LAST 16 11 T T gt T 4 i M 56 K ch CALL PRNT F NF F 1
65. 3052D 00 1 0339835D 00 2 21052460 00 1 00 2 5000002D 00 6 7846557D 00 1 50 4 2500002D 00 1 24040010 01 6 861206800 00 1 21859130 2 00 6 0000003D 00 1 9619383D 01 2 1972644D Ol 2 2361699D 2 50 2 82922420 01 5 9904692D 01 4 26147360 3 00 6 25844840 5 19287560 01 3 50 9 68767270 01 2 065301790 01 6 7 62400570 00 1 6201563D 02 5 7062075D 01 2 9312866D 02 1 3857167D 02 00 2 1972644D 2 57990800 01 4 26147360 01 3 49872850 01 5 1928756D 01 6 8828247D 01 2 6530179D 01 2 5687388D 02 2 2361699D 01 2 5799080D 01 6 0000003D 00 1 96193830 01 6 00000030 00 4 1981082D 01 7 6240057D 00 3 49872850 01 2 82922420 01 5 10620750 01 6 25844840 01 1 3857167D 02 9 68767270 01 Output Concluded Figure 2 Concluded 667528697 01 6 8828247D 01 2 19077700 02 2568173880 02 4 3170034D 02 6 8584482D 01 2 8292242D 01 6 8584482D 6 25844840 01 6 85844820 01 9 68767270 6 85844820 01 8 12198490 01 5 99046920 01 8 12198490 01 1 62015630 02 8 12198490 01 2 9312866 02 8 12198490 01 The user s main program A brief explanation of the statements using line numbers on figure 2 a as reference follows Lines 1 2 and 3 Lines 1 2 and 3 allocate storage same as lines 1 and 2 of example I Line 4 Common variables to be needed later are made available to the progr
66. 3204 CALL EQUATE W NN PHI NTH4 NN RETURN 0CI 600 CALL LNCNT 2 WRITE 6 1600 O04 1600 FORMAT DIMENSION ERROR RICAT T35 NPHI 7X4 NC 9X NPT 1 16X9 3X KDUM NIN 29X 3 3X 214 4X T4 5 14 CALL ASPERR RETURN END 11 SUBROUTINE 5 NPHI P MP OK SNK S RICONT DIMENSION NH NG INP CZ 4 TATA DOUBLE PRFCISTON PHI 1 HC1 0 1 RUIT 9PCL 4KC1 DUMCIT DOUBLE PRECISION SUM SUMN AL C DIMENSION OF DUM MUST AT LEAST 6 N N C CHECK FOR CONFORMABLE MATRICES N NPHI 1 M NH 1 NK 1 N NK 2 M NSQO N N ND 1 1 ND2 NSQ 1 02 50 50 6 6 5 0 IF 2 NH 2 NE 1 OR 2 OR 1 1 2 OR 1 NLOR NP 2 NE OR KDUM 1 I 5 06 260 0 900 NFIN NCONT 2 NPRINT NC ONT 1 ONT 3 C START OF MAIN COMPUTATION P O IS INPUT DATA IN P MATRIX KLN 0 JFN 0 1 0 1 100 CALL MULT H NH P NP DUM NZD DUM H P CALL TRANP ND2 NZ2D 23 C DUM2 HPRIME CALL MULT DUM NZD DUM ND 2 N7D 3 DUMCND3 NZD 5 C DUM3 CALL ADD DUM ND3 NZD 5 R NR DUM NZD C DUMZH P xHPRIME TR CALL PSEUDO DUM NZD DUM ND2 9NZDC3 S DU
67. 4 TT4DELTI DELT2 TEINAL 2 DET 0003 COMMON MAX MAXRC 66V _ _ _ _ _ __________ C _C_ _C_C_ _ _ lt 0005 CALL RDTITL 0006 MAXRC 36 DOT u u 0008 READ 5 100 TT DELT1y DELT2 TFINAL 0009 CALL READ 5 DOO A READ 2 2 20 0 2 2 00 11 100 FORMAT 8F 10 2 0012 101 FORMAT 1H 59X TIME RESPONSE 5 22 5 43 1 OUTPUTt 5Xt TT 14X4 XT 1 41 1X EX 2 LIEK 1 EX F433 1 21X 4t DEL TC 1 0013 102 10 2 6 3 16 7 10 2 16 7 10 16 7 0014 1 1 GOLS r I NCONT 2 TEINA DELTI he 0016 NCONT 3 4 0017 CALL EQUATE XO NXO XT NXT 0018 CALL EQUATELPO NPO PT NPT 0019 CALL FOUATECR RI NRI 0020 CALL NRI DET L CAL ANGLE GSNGSRTSNRIGH NH QNO C ANC alle 0022 CALL 2 NZ DELT1 NPHI 0023 CALL RICAT PHI NPHI PT NPT DUMMY KDIIM 0024 CALL MULTAG 0025 CALL 5 1 1 2 2 1 no 0026 CALL ADD F 2 NA2 FS TAR NFSTAR 0027 GALL PRNFLESTAR NAS FAR gL RS FR LY mou
68. 5870 01 6 2633587D 01 0 0 6 6370228 00 Output Concluded Figure 4 Concluded SE Example 5 Matrix Decomposition This example is a test program to check the operation of DECGEN IL first generates 2 matrix R to be decomposed then proceeds with the decomposition and checks the result printing all of the associated matrices The general procedure 15 to input diagonal matrix ZL and transform it into the matrix R 0 decomposed Figure 5 a 15 a listing of the main program figure 5 b is a listing of the subroutine figure 5 is the data deck and figures 5 d through 5 L are the output In the main program all matrices are dimensioned 100 although the actual matrix size used is 2 X 2 4 X4 Accordingly MAXRC is set to 100 The dummy matrix is dimensioned 700 since DECGEN requires that much The input matrices are read at line 8 Subroutine ORTH called at line 9 produces a n X n orthogonal matrix using the original T matrix and places the results back in T The procedure is as follows First generate an elementary rotation matrix This is a unity matrix with elements and ejj replaced by cos tjj and elements ejj ejj sin tjj Then T II Eg Lines 10 through 17 set up indices for referring to the seven dummy matrices The input matrix ZL is then transformed by the matrix T so that ZL T ZL T Note th
69. 7 105 105 II II 1 GO 106 XII 107 IF K 109 108 212 109 CALL LNCNT 1 WRITE 62110 110 FORMAT t ERROR IN ETPHI K IS NEGATIVE 112 IF 1 213 212 212 213 K 1 212 DO 111 J 1 K T 2 gt T CALL EQUATF R NB DUMMY 1 ND CALL FEOUATF DUMMY 1 ND DUMMY NDD ND 111 CALL MULT DUMMYCNDD ND DUMMYC1 NDB NR TT T 108 CONTINUE S 1 D0 SC CALL SCALF A NA Ay NA S RETURN 102 CALL LNCNT 1 WRITE 6 119 119 FORMAT ERROR IN K 1000 CALL ASPERR RETURN 998 CALL LNCNT 1 WRITE 6 50 NA KDIIM NDD 50 FORMAT DIMENSION ERROR IN FTPHTI NA 216 1 I1 1IM4 1 15 5 1 IKDUM MIN 15 CALL ASPFRR RETHRN SII SUBROUTINE AUG F NF G NG RI NRI H NH O0 NQ 7 7 TI DIMENSION F 1 GCV RE 1 HCV 32 QC1 72 1 2 CC1 DIMENSION NNP1 2 NNP2 2 NNP3 2 NNP4 2 NF 2 NG 2 NRI 2 1 2 7 2 2 2 0 2 DOUBLE PRECISION Gy 0 2 IF NF NE NF 2 OR NRI 1 NE NRI 2 f R 1NQ 1 NE NQ 2 GO TO 995 2 NF 1 I NG 1 NE NF 1 OR 2 NE NRI CI GO TO 995 IFCII EQ0 1 GO TO 206 IF NH 1 NE NQ 1 DR NH 2 NE NF 2 GO 995 206 TWO 2 N 1 50 NZ 1 NNZ NZ 2 MM7 1 1 2 1 M5 NP3 DR
70. 70 00 1 18438510 02 RANK MAT REX 1 ROMWS 3 COLUMNS 2 0000000D 00 f Case 7 Conciuded Figure 5 Concluded 5 65837100 01 finally figure 5 f case 7 is a 4 X 4 matrix of rank 3 with one very small eigenvalue equal to 10 The error matrices of the first two decompositions are extremely small but that from the third one has errors of the order of 10 These are caused by the built in pivot rejection device which rejects all pivots smaller than 2X10 times the largest of the diagonal elements see DECOM p 85 and PSEU p 70 This last matrix case 7 was also tried with an eigenvalue of 10 and the errors were then on the order of 1015 Example 6 Use of the Pseudoinverse Routine This program is designed to check the operation of PSEUDO The procedure is as follows First the input matrix A isread then B is computed The accuracy of the pseudoinverse is then checked by the first two Moore Penrose axioms BAB B A ABA A B All the various matrices are printed Figure 6 a is the program listing and figure 6 b the output Three cases are presented the first two are the examples presented in the ASP manual and the third one contains several zeros The first matrix printed for each case is the input matrix and each has a different label The other titles are abbreviations chosen to fit the allotted four character space as follows APSE gt
71. 78 0 8844306D 0 R672834n 0 8493472n 0 8306919N 48113857D 914 amp 954 0 1710861 0 75022120D 0 728 9623 0 7073689D 0 68549A87N 0 6634074D 0 6411488 0 61877430 0 5963334n 0 57387340 0 551 3960 0 5290751 0250682060 0448471510 06462 7957D 0 4410954n 0 4196482n 0 398R4B841n 0 3776314fh 0 3571164n 0 3369636n 0 3171955D 0 2978327 STATE XT 2 01 0 0 00 0 16739030 00 0 32773580 00 0 48093551 00 0 6269159D 00 0 18562920 0 8970518D 00 0 11380420 0 12476710 00 0 13501260 00 0 14454850 00 0 1533838D 00 0 1615292D 00 0 16899660 00 0 17579910 00 0 18195110 00 0 18746780 00 0 19236540 0 19666100 00 0 20037200 00 0 20351700 00 0 2061147 00 0 20818440 00 0 20974581 0 2108187 00 0 21142340 00 0 21158030 00 0 21130960 00 0 21063190 00 0 20956760 0 208137 1 00 0 20636060 0 2042582D 00 0 2018497D 00 0 19915470 0 1961926 0 19298220 00 0 18954240 00 0 185 891 3 ET TIME RESPONSE XT 3 0 0 0 49775060 0 1 0 98343630 0 14558150 0 2 3563540 0 2 TR2685D 21920330 0 35837720 0 2957379 0 4312427 0 464850 840 0 496560 4D 0 52633220 0 55416490 0 5800568D 0 60401260 0 62604350 0 64616620 0 66440255 0 68077930 0 6952277 0 1080831 0 7190843n 04 7283734n 0 7359
72. 8 0 Output Concluded Figure 6 Concluded APPENDIX A DESCRIPTION OF INTERNAL SUBROUTINES 25 READ DESCRIPTION This subroutine reads a single matrix from cards without a header card It is called by READ after the latter has read the header card The dimensions of the matrix to be read are in array NZ If this is zero no array will be read In any event the routine then prints either the array just read using NZ for dimensions or if NZ 0 the array already stored using NA for dimensions The subroutine reads the data from cards each row of the matrix starting on a new card using format 8F10 2 If the card data is in exponential form it must usea D exponent USAGE CALL READI A NA NZ NAM Input Arguments Matrix A if NZ 0 Dimension array NA NZ Constant NAM containing a four character or less name for the matrix which will be used by PRNT Output Arguments Matrix A if NZ 0 Dimension array NA 26 ASPERR DESCRIPTION This is an installation dependent subroutine It is called by the various subroutines when they detect an error It is intended to provide an error walkback so that the programmer can determine which call of a given subroutine is in error It also counts the number of errors and calls EXIT after ten entries into ASPERR 69 USAGE CALL ASPERR It has no arguments The user may if he wishes call
73. 9 CALL LNCNT 1 WRITE 6 50 50 FORMAT DIMENSION ERROR IN TRANP 216 CALL ASPERR RETURN F ND Y 101 SUBROUTINE DIMENSION 1 1 1 2 DOURLE PRECISION RIGA HOLD C OMM ON MAX MAXRC IF 1 2 GOTO 600 C SEARCH FOR LARGEST ELEMENT 1 N NA 1 50 IF N LT 1 0R NSQ GT MAXRC GU TO 600 NK 2 N DO 80 K 1 NK NK N L K K NPK N K K KK NK K BIGA A KK 20 J K N IZ N J 1 20 I K I IZ I 10 IF NABS RIGA DABS A TJ 15 20 20 15 RIGA A TJ L K 1 NPK N K L NPK J 20 CONTINUE C INTERCHANGE ROWS J L K IF J K 35 35 25 25 KI K N DO 30 I 1 KI KI N K J JI L Ip I J OLD INTERCHANGE COLUMNS 35 NPK N K LI L NPK I K 45 45 38 38 JP N 1 NO 40 J 1 N JK NK J JI JP J HOLD A JK A JK ACJT 40 ACJT HOLD C DIVIDE COLUMN RY MINUS PIVOT VALUE GF PIVOT ELFRENTS IS CAN TATARD C BIGA 45 IF BIGA 48 46 48 46 0 CALL LNCNT 1 KKK KK 1 WRITE 6 1046 KKK 1046 FORMAT INV ERROR OF 0 RANK OF A 14 CALL ASPFRR RETURN IF I K 50 55 50 50 IK NK I A IK BIGA 55 CONTINUE C REDUCE MATRIX DO 65 1 N IK NK J Eos 8 IJ 1 N N IF I ra K 604 65 60 60
74. 95Rn 174643340 0 7493510p 0 7508056D 75085275P 71495484N 0 7469508D 74311 78D 07381083P 72747961N 0 7166119n 0 7074875n 0 697658 70 c Output Concluded Figure 3 Concluded 01 00 00 00 00 VASP PROGRAM OUTPUT DELTC 0 6754622805 0 78731170 0 816825R8D 04 8432503D 0 9048R47D 0 9198623D 0 94208020 0 949459 0 95454 15D 0 95741470 0 95816290 0 95589360 0 95370790 0 948 7064 0 92479R4D 0 9125195 0 8999112 0 88R606650D 0 87107640 0 8550301D 0 8 3801440 0 82011440 0 7819890D 0 7619218D 04 7412862n 0 72015510 0 6985987D 04 6544T7R4f 0 63204200 0 609435 0 58671560 20903720 04 5411522Dp s The user s main program Some of the details of the main program discussed briefly The various matrices are first dimensioned and stated to be double precision The problem will be solved using basically 3 X 3 matrices but Z and PHI are 6 X 6 matrices so MAXRC is set to 36 line 6 double size dummy array is required in ETPHI so DUMMY is dimensioned at 72 and KDUM is set to 72 line 7 line 8 the timing information is read is the initial time DELTI is the time increment used in the computation of P DELT2 is the time increment 7 desired in the printout of the transient and TFINAL is the final time for the transient Lines 9 and 10 read data card
75. ATE MULT TRANP EQUATE ADD MULT INV EQUATE PRNT ADD SUBT MULT TRANP PSEUDO PSEU BDNRM ANDRA PRNT ADD EAT SUBT MULT PRNT EQUATE PSEU BDNRM ANDRA MULT TRANP INV NORM EQUATE DECOM PSEU BDNRM ANDRA PRNT None None MULT NORM BDNRM ANDRA MULT NORM LNCNT MULT NORM PSEU BDNRM ANDRA NCNT and ASPERR are additional external references 1 READ DESCRIPTION This subroutine reads 1 to 5 matrices from cards along with the names and dimensions prints the same information For each matrix the routine first reads a header card containing a four character title followed by two integers giving the row and column size of the matrix using format 4 4 214 Then the matrix data are read using 1 each row of the matrix starting on a new card using format 8F10 2 II the card data is in exponential form it must use a D exponent The matrix title and the matrix are then printed using subroutine PRNT If the header card contains no row and column size i e n 0 then the matrix in storage is unchanged no data cards are read for that matrix and the previously stored matrix 15 printed USAGE CALL READ I A NA B NB C NC D ND E NE Input Arguments Control constant I where I 15 integer from 1 to 5 and indicates the number of matrices to be read If I is less than 5 the extra matrices in the call list may be dummy variables or repeated refe
76. CALL UNITY A NA 15 Input Argument Dimension array NA Output Argument Matrix A 13 TRCE DESCRIPTION This subroutine computes the trace of the matrix A n TR 2 1 1 USAGE CALL Input Arguments Matrix A Dimension array NA Output Argument Scalar TR 14 EQUATE DESCRIPTION This subroutine copies the values stored in matrix A into matrix B as follows B A USAGE CALL EQUATE A NA B NB Input Arguments Matrix A Dimension array NA 16 Output Arguments Matrix B Dimension array NB 15 JUXTC DESCRIPTION This subroutine takes the m X n matrix A the m X p matrix B and forms the m X n p matrix C A USAGE CALL JUXTC A NA B NB C NC Input Arguments Matrices A B Dimension arrays NA NB Output Arguments Matrix C Dimension array NC 16 JUXTR DESCRIPTION This subroutine takes the m X n matrix A the p X n matrix and forms the X matrix USAGE CALL JUXTR A NA B NB C NC 17 Input Arguments Matrices A B Dimension Output Arguments Matrix Dimension array NC 17 EAT DESCRIPTION This subroutine computes eAt t c f eAT O For a linear time invariant system the system equation is x Gu Then t x t ety f eAT dr Gu 0 0L x t Bx CGu See ASP manual page 92 for reference USAGE CALL EAT A NA TT B NB C NC DUMMY KDUM Input Arguments Matrix A Dim
77. DCM IF ODCM LT 02 QITR ODCM 1 NR QR QIT 17 XC I TEST TO SEE IF SYMMETRIZATION IS NEEDED IF OPS 16 150 16 C TEST TO FIND SMALLER DIMENSION MATRIX 16 IF QR NCA 18 18 19 19 NR NCA QX 1 QR2 QR QLL OR 17 60 170 18 CONTINUE QX NR OR2 1 QLL NCA QTP 1 170 CONTINUE SET ROW COLUMN LIMIT TO APPROPRIATE CASE EITHER ROW OR COLM DIMENS 81 I NR DO 181 K i NR DSUM 0670 Q QN QX J QR L QN QR2 I M 2 183 DSUM DSUM A L A M C SUM A I LL LL RUNS 1 TO ROW OR TO COLM LIMIT AXA TRANS HAS COLM LIMIT B A TRANS A HAS ROW LIMIT L K 1 NR I 181 B L DSUM 60 188 C HERE MOVE TO P A IS ALREADY POSITIVE DEFINITE 150 DO 151 L 1 QNT 151 B L A L FORCE SYMMETRIZATION OF TO COMPENSATE FOR ROUND OFF MULTIPLIC 188 DO 189 I 1 NR DO 189 K 1 NR I K 1 xNR M K I 1 NR C 8 L 8 M 0 5D0 DSUM DSUM B L C HERE SET UP CALL INITIAL OF ANDRA ONLY COMES HERE ONCE PER MATRIX QNT NR NR KRC2 QNT KRC QNT KRC2 KRC4 QNT KRC3 C OMIT SAVING OF B IF RANK ONLY AND NO ITERATION IFC IP 2 NE IZ AND QITR EQ 12 GO TO 200 DO 1891 I 1 ONT 1891 D I BT 200 C ONT INUE C SEARCH DIAGONAL OF INPUT FOR LARGEST ELEMENT USE TO DEFINE FL QR1 NR 1 L 1 DFZD M Z DO 23 I 1
78. DIAGONALIZATION RETURN 691 CALL LNCNT 1 WRITE 6 1691 QR NCA 1691 FORMAT 1 DIMENSION ERROR IN PSEU NA 216 50 1700 692 CALL LNCNT 1 WRITE 6 1692 1692 FORMAT ERROR PSEU DIAGONAL ELEMENTS OF MATRIX 0 60 TO 1700 693 CALL LNCNT WRITE 6 1693 1693 FORMAT ERROR IN PSEU RANK HAS DECREASED COMPUTATION ENDED c ana Nun 1 38 dd dSV 1v9 00 1 13215 XId1VN d 1V d9 SI XXV 0354 NI Odd 1vWuOd 9091 9091 9 LIUM 895 01 09 dddSV 11V9 SdfHNEITNU _NUEIVIildMU d SV3dO NTI_SVH YNVY NASd NI i VMdU 769 7691 49 ALIN 1 LNON TIVI 769 0041 09 141 CVI FUNCTION BDNRM NR CT FE D KRV INTEGER 2 QF DOUBLE PRECISION AN BR DFZO DSIIM DIMENSION CT 400 EE 400 D 2000 NV 2 KRV C D HOLDS 5 MATRICES THE FIRST AND THE LAST 2 ARE USED HERE DIMENSION PPP 2 FQUIVALENCE AN FN BR FR DATA DFZO 0 D0 EQUIVALENCES BELOW JUST TO SAVE STORAGE EQUIVALENCE DFZ0 1Z DSUM CBR PPP Y 1 PPP 2 K 5 1 MV 1 L 2 CIR NL C TEST IF NR TRANSPOSE ROLES OF D AND CT TRANS CT QF NR 05 RRVU3 OOOO KD4 KRV 4 FC NR 0 10 20 ENTRY TTRM NR CT EE C TO DO T TR T ONLY ENTRY TTRM QF IZ GO TO 20
79. DUMMY 1 DUMMY NDD NA C NC ALL ADD DUMMY NDD NC Cy NC Di 111 CALL MULT DUMMY 1 DUMMY 1 NB B NB TT T 108 CONTINUE CALL EQUATE DUMMY NDD NC C NC S 1 N0 SC CALL SCALF A NA A NA S RETURN 102 CALL LNCNT 1 6 119 119 FORMAT ERROR IN EAT K 1000 CALL ASPERR RETURN 998 CALL LNCNT 1 WR 8 NA K DUM amp N DD 50 FORMAT DIMENSION ERROR 216 K 12X9 1 5 CALL ASPERR RETURN END NI CI SUBROUTINE BIMENSTON A TY BCL 5 DUMMY 1 NACZ DOUBLE PRECISION ANAA TMAX DUMMY S SC COMMON MAX MAXRC NR NA 1 NCC NA 2 NB 1 NR NB 2 NCC 1 NCC 1 OR LD 6T MAXRC GO TO 998 NDD 2 NA 1 1 IF KDUM LT NDD GO TO 998 NDD NA 1 NA 1 1 T TT CALL 10 0 1 ANAA 0 101 IF TMAX T 103 104 104 103 K K T TT 2 K IF K 1000 101 102 102 04 SC T CALL SCALE T CALL UNITY B NB LI 2 M 35 CALL 1 CALL EQUATE CA NA DUMMY NDD ND 106 CALL MULT As NA DUMMY NDD ND B NB S 1 DO 1 1 CALL SCALE R NB DUMMY NDD ND S CALL ADD DUMMY NDD ND DUMMY 1 ND NB CC CALL EQUATE CB NB DUMMY 1 ND N N 1 IF N 107 10
80. E C DUMMY DIMENSION A 1 B8 1 DUMMY 1 NAC 2 NBC2 NDC2 9CCT 9 NCC2 DOUBLE PRECISION A TMAX DUMMY 5 SC COMMON MAX MAXRC_ NR NA 1 NCC NA 2 NC C 1 NR NC 2 NCC NB 1 NR NBt 2 NCC LD NR NCC TF 1 eR LD GT MAXRC GO T 998 25 NA 1 1 IF KDUM LT NDD GO 998 NDD NA 1 XNA 1 1 T TT CALE ER 10 01 ANAA K 0 101 IF TMAX T 103 104 104 103 K K 1 ROK K IF K 1000 10 1 102 102 104 SC T CALL SCALE AyNA DO 401 I 1 2 821 NB T 1 CALL UNITY B NB CALL SCALE ByNByDUMMY 1 NB T S T 2 CALL SCALE A NA DUMMY NDID NA S N 35 1152 CALL ADD DUMMY 1 4 NA DUMMY NDD NA DUMMY NDD NA CALL ADD A NA B NB DUMMY 1 CALL EQUATF A NC 106 CALL MULT A NA C NC B NB CII Sz1 00 IT CALL SCALER NC S 6 141 CALL SCALF C NC B NB S CALL ADD amp NB DUMMY NDD DUMMY NDD NR CALL ADD C NC DUMMY 1 NC DUMM Y NC 1 107 107 105 105 11 11 1 GO TO 106 107 CALL EQUATE DUMMY 1 NB B NB IF K 109 108 212 109 CALL LNCNT 1 WRITE 6 110 110 FORMAT ERROR EAT IS NEGATIVE 112 IF K 1 213 212 212 213 K 1 212 111 J 1 K CALL EQUATF P DUMMY 1 CALL MULT
81. E OF 1ST MATRIX 1 ROW SIZE OF 2ND MATRIX 6 991 2 NG 1 WRITE 6 992 NJ 2 sNR 1 WRITE 6 993 2 1 WRITE 6 995 2 1 WRITE 6 994 1 1 2 2 9 99 FURMAT 1HO 992 FORMAT 1HO INTEGRAL EXP F T G 113 20 18 J R 117 115 20 18 993 FORMAT 1HO YK X 17 115 20 18 994 FORMAT 1HO 1 SIZE OF R 995 FURMAT 1HO 996 FURMAT 1HO GO TO 1000 52 FORMAT 1HO KX ROW SIZE OF J 15 I5 12X 0L K IS 15 23X COL IS 15 3 OF X IS I5 X 17X T115 20X 18 EXPCF T X 10 115 20 18 DUMMY MUST BE DIMENSIONED LEAST I4 X 1 BUT IS 1 DIMENSTONED ONLY I4 X 1 GO TO 1000 910 WRITE 6 52 MAX KDUM 1000 CALL ASPERR END CI SUBROUTINE PSEUDO A NA B NB DUM KDUM DIMENSION A 1 B 1 DUM 1 2 2 IP 4 DEP 2 DOUBLE PRECISION A B DUM DEP 1 2 10 I 1 2 1 0 0 10 0 20 IP 3 NA 1 IP 4 zNA 2 1 2 2 1 II KDUM LT NNW GOTO 999 NEEZNA 1 NA 2 1 ND 2 NEE 1 IS 1 1 MATRIX ROUTINE INVERTS 1 AND PUTS IT IM B CALL PSEU A DUM NEE DEP TP DUMEND GO TO 1000 x uL 60 TO 1000 999 WRITE 6 500 KDUM NNW 500 FORMAT DIMENSTON ERROR PSEUDO KD uMz T5 3X KDUM MIN Z T5 1000 RETURN END gcl C SUBROUTINE
82. FOR NIG ELEMENT MALREADY REDUCED 1 CODEMARKED ICC IF IDD ICC GOTO 30 IF FDI FMX 30 30 32 C TEST FOR NEGLIGIBLE FL THESE AND AS ZEROS 32 IF FDI EF GO TQ 30 SET NEW 2RLE PREC SAVE BEST ROW FOR PIVOT GO TO 30 40 CONTINUE 400 IF M 0 17 GO TD 505 DRS 1 DO DSQRT DMX C SET INDEX OF FIRST ROW OMR PIVOT COLUMN K 1 QS l OMR DO 41 1 1 OS DDM b L DRS T L T L DRS B L DDM C SYMMETRICALLY FORCE COLUMN TO SAME VALUE IN B ONLY B K DDM K K 1 61 L 05 L C FORCE PIVOT ELEMENT TO EXACT VALUE UNITY B M 1 06 NOW REDUCE ALL OTHER ROWS OF B T ELIMINATING COLUMN OF PIVOT VARIAB DO 460 I 1 QS C TEST FOR PIVOTAL ROW OTHER ROWS IF I EO OMR GO 460 COEFF TO ZERQED CAN NOT PREVIQUS J I OS K 052 J DRS B K TEST FOR ROW ALREADY REDUCED TO SKIP IF 11S EQ ICC GO TO 460 GET COEFF PIVOT COLUMN TO BE ELIMINATED K QMR QS J DMM B K L QMR K I DO 47 J 1 QS C L IS ROW USED TO REDUCE WITH PIVOT IS CURRENT ROW THAT PIVOT GETS ELIMINATED FROM B K B L DMM T K T K T L DMM L C5 L 47 K QS K 460 CONTINUE QMR 461 I 1 QS FORCE MOST OF PIVOT ROW TO ZERO COMPLETES REDUCTION WITH 1 B L DFZO 461
83. I a diagonal matrix with all 1 s moved to the uppermost left corner Given these matrices and the ability to find a pesudoinverse of A a decomposition of any matrix is possible PSEU and DECOM are called and the matrices then multiplied as described in the method to give a Kronecker decom position The routine calls PSEUP and INV to do most of the calculation Besides returning the matrices P and E it does nothing that could not be done by successive calls of other matrix routines It has parameters and error exits similar to that of PSEUP USAGE CALL DECOM A B C E J DCM KP D Arguments Description The symmetric non negative definite input B The output matrix E diagonal of and 1 with 1 s in the independent columns E J D and D1 are all of same size as A The output such that TAT diagonal of 0 and 17 87 1 2 DCM3 KP 1 2 4 88 Holds the inverse of B does not hold the inverse of A Not E of ASP A square integer matrix for housekeeping in INV and DECOM Parameters just as in subroutine PSEUP Multiplied times the largest magnitude of diagonal Of A to give a lower limit for an acceptable pivot PSEUP Set at 2 10 if zero is input Used only if the user elects to iterate in PSEUP Set at I DO no effect if zero is input Note DECGEN uses the default options for 1 and DCM2 The last pivot accepted by
84. IF J 62 65 62 62 KJ IJ I A 1J HOLD A KJ ATA 65 CONTINUE C DIVIDE ROW RY PIVAT KJ K N DO 75 J 1 N KJ KJ M IF K 70 15 10 TO A KJ RI GA 75 CONTINUE C PRODUCT OF PIVOTS DET DET BIGA C REPLACE PIVOT BY RECTPRRCAL A KK 1 RIGA RO CONTINUE C FINAL ROW AND COLUMN INTERCHANGE 100 K K 1 150 150 105 105 I L K IF T K 120 120 108 108 JO N K 1 JR N T 1 DO 110 J 1 N JK JQ J AK 7 JI JR J A JK AC JI 110 A JT HOLD 120 NPK N K J L NPK IF J K 100 100 125 125 KI K N 130 I 1 N KI KI N HOLD A KI JI KI K t J 130 A JI HOLD 100 vol 150 RETURN 600 CALL LNCNT 1 WR 6 1600 NA 1600 FORMAT 1 INV REOUIRES SQUARE MATRIX NA 216 CALL ASPERR RETURN END SOT SUBROUTINE NOBMLA NA ANORH K 3 DOURLE PRECISION A ANORM SUM RAWNAX CL i AX COMMON MAX MAXRC 1 NAC NA 2 L NAR IF 1 1 1 1 1 1 GA 999 ROWMAX 0 DO 300 I 1 NAC SUM 0 p coat 53 5 l eee K K 1 301 SUM SUM DABS A K IF CULMAX LT SUM COLMAX SUM 300 CONTINUE DO 302 I 1 NAR SUM 0 I 303 J 1 K NAR 303 SUM SUM DABS A K IF ROWMAX LT SUM ROWMAX SUM 302 CONTINUE ANOR
85. IX 1 ROWS 1 COLUMNS 200000000 00 a ois ms uuu ET e Case 4 Concluded Figure 5 Continued ON t2 TEST PROGRAM FOR DECGEN AND DECOM CASE 7 ILL COND 4X4 RANK3 VASP PROGRAM mm __ ___3 ______ __ _ ___ _ ___ __ ___ ZL MATRIX 4 ROWS 4 COLUMNS 1 00 00 0 0 5 0 0 0 0 0 2 00000000 00 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 8 0 1 9000900905 56 5 5 5 5 5 5 5 MATRIX 4 ROWS 4 COLUMNS 3 93291460 01 7 68929180 01 2 066722 30 01 2 59735070 01 2 64263040 01 5 48465710 01 6 97672110 01 3 11629410 01 1 4155739D 02 2 8968030D 01 64 3928160D 01 71 08452 15D 01 EXR MATRIX 1 ROWS 1 COLUMNS eS C XC MATRIX 1 ROWS 1 CULUMNS T MATRIX 4 ROWS 4 COLUMNS 3 8 a UU L DU J DU 9 6 68042710 02 8 98175610 01 1 37761430 01 4 12115950 01 1 16886670 01 1 21332550 02 9 4444935D 01 3 06905200 01 36803160 03 4 1690464D 01 H 93174226B 03 173110828 613 T T MATRIX 4 ROWS 4 COLUMNS 5 1 9090909009 00 2 gt 77555768 17 1 28777880 17
86. If no pivot is found skip to step 3 The square root of the pivot is taken The pivot row and pivot column are divided by the square root Thus the pivot at the intersection of the row and the column is reduced to unity The corresponding row of T is also divided by the square root Step 2 The new reduced pivot row is used to eliminate the elements of the pivot column Let K be the pivot row and column The pivot row is multiplied times the element in the j k position The resulting row vector is subtracted from the jth row This process is repeated for each row j that has not yet been a pivot row Exactly the same operations are carried out on the corresponding rows and columns of T except that of course the multiplier for a pivot row comes from B Then the pivot row of B except for the pivot is set to zero The pivot row and its corresponding row in T are never used again Step 3 If the rank is maximal exit If no pivot has been found a test is made to see if this should be an error exit or normal exit Otherwise repeat step 1 Computational In practice a number of modifications are made The actual calculations are rewritten to optimize speed and storage The reciprocal of the square root is used instead of a division For single precision straight division would probably be best In step 3 an artificial code is put into the pivot position This code is chosen as one that cannot be the result of floating point arithmetic
87. KRV 2 IR C WHY ONLY SNGLE PREC THIS IS ONLY A ROUGH TEST TO STOP ITERATION C 5 WHY SIMILARLY OTHER USES OF SINGLE PREC IF SOLD LT SS SOLD FZRO GO TO 650 C IF SUBSTANTIAL IMPROVEMENT TRY AGAIN OTHERWISE QUIT RETURN JHE A PSEUDO INVERSE EVEN IF OFF 220 CONTINUE QIT 01141 SOLD SS SAVE PREVIOUS ROW WHICH A PIVOT WAS FOUND M JP 5 IF QIT EN NR GO TO 700 PUT IN SMALLER TOLERANCE CASE DIAGONAL TOO SMALL QTHERWISE DPR 1 DPR 2 2 0 5 C TRY TO REDUCF 1 MORE ROW I IL NR 30 700 606 650 C ONT INUF RESTORE AND C TO THEIR PREVIOUS VALUES THE LAST PIVOT HAS BEEN CREJECTED BACK TRACK WHILE ITERATING LEI JP 3 JP 3 1 DO 652 I QNT J KRC J K KRC2 T B I D J 653 CCI D K 700 CONTINUE IR JP 3 M 7 HERE WISH REPLACE MARKERS IN DIAGONAL WITH LEGITIMATE 1 00 L 1 DO 704 I 1 NR DDI B L IF IDD 701 702 N TO 710 B L 1 00 60 10 704 7101 FORCE SMALL TRASH ZERO 7101 B L DFZO 702 Ms 704 L 061 L C IF ALREADY TRIED ANOTHER REDUCTION TO GET MATRIX IN UPPER DIAG COR OMIT PART OF CALCULATIONS IF ONLY RANK IS DESIREN IFC IP 2 IZ GO TO 877 M SUPPR PHAS E WAS CALLER IF M LT 1 OR QDCM LI OZ GO TO 80 C BELOW HAVE SING MATRIX THAT NEEDS FURTHER WORK C HAVE MATRIX DIAGONALIZED WITH 15 OS INTERSPERSED IS SINGULAR
88. L NASA TECHNICAL MEMORANDUM Gd gt lt RETUF L ru BE AFWL DOY L u lt KIRTLAND Hr USERS MANUAL FOR THE VARIABLE DIMENSION AUTOMATIC SYNTHESIS PROGRAM VASP by John S White and Homer Q Lee Ames Research Center ie Moffett Field Calif 94035 TECH LIBRARY NM INI 1 Report No ii 2 Government Accession No 3 Recipient s Catalog No NASA TM 241 74 Title and Subtitle 5 Report Date USERS MANUAL FOR THE VARIABLE DIMENSION AUTOMATIC October 1971 SYNTHESIS PROGRAM VASP 6 Performing Organization Code 7 Author s 8 Performing Organization Report No John S White and Homer Q Lee A 3882 10 Work Unit No 9 Performing Organization Name and Address 125 19 20 02 11 Contract Grant Moffett Field Calif 94035 13 Type of Report and Period Covered 12 Sponsoring Agency Name and Address Technical Memorandum National Aeronautics and Space Administration Washington D C 20546 14 Sponsoring Agency Code 15 Supplementary Notes Abstract VASP is a Variable dimension Fortran version of the Automatic Synthesis Program a computer implementation of the Kalman filtering and control theory It consists of 31 subprograms for analysis synthesis and optimization of complicated high order time variant problems associated with modern control theory
89. M DMINI COLMAX ROWMAX RETURN 999 CALL LNCNT 1 WRITE 6 50 50 FORMAT DIMENSION ERROR IN NORM NA 2 16 CALL ASPERR RETURN END 901 SUBROUTINE UNTTY A NA DIMENSION 1 2 DOUBLE PRFCISTON A 1 2 GO TO 999 CALL SCALF A NAS 0 DO J NA 1 1 T SD GOD J NA 8 J NAX J 1 300 A J 1 60 1000 999 CALL LNCNT 1 MRITE 6 5 1 1 2 50 FORMAT DIMENSION ERROR IN UNITY NAz 216 CALL ASPERR 1000 RETURN X v P END m LOT SUBROUTINE TRCF A NA IR DOUBLE PRFCISION A 1 TR DIMENSION 2 COMMON _ MAX MAXRC IF 1 2 GO TO 600 TR 0 DO NzNA 1 MM Nex N Na LT 1 HR NN GT MAXRC 50 600 DO 10 Iz1 N MzITNx I 1 10 TR TR A M RETURN 600 CALL LNCNT 1 WRITE 64 1600 NA 1600 FORMAT TRACE REQUIRES SQUARE MATRIX 1 216 CALL ASPERR RETURN END 801 SUBROUTINE FQUATE CA NA B NB DIMENSION 1 1 2 2 DOUBLE PRECISION A B COMMON MA X MAXRC NB 1 M 1 M8 2 2 L NA 1 NA 2 IF NA 1 LT 1 DR L LT 1 OR L GI MAXRC GO 999 00 300 I 1 L SA 1000 RETURN 999 CALL LNCNT 1 WRITE 6 50 NA 50 FORMAT 1 DIMENSION ERROR IN FQUATE 216 CALL ASPERR RETURN M h
90. M ND3 KDUM 3 NS 0 C NUM2 INVERSE CALL TRANP C DUM HPR IML CALL MULT DUM NZD yDUM ND2 NZD 3 DUM ND3 NZD 5 CCI CALL MULT M M M0 3 M7 5 DUM NZD DUM A 110 CALL MULT PHI NPHI DUM NZD Ky NK CALL MULT DUM NZD H NHsDUM ND2 NZD 3 C DUM2 CALL MULT 2 7 3 P NP DUM NZD DUMz AxHxP CALL SURT P NP DUM NZD DUM NZD C DUM P A H P CALL TRANP PHI NPHI DUM ND2 NZD 23 CALL MULT ND2 NZD 3 DUM ND3 NZD 5 CALL MULT PHI NPHI DUM ND3 NZD 5 DUM NZD CALL ADD DUM NZD4 0 NQ4 P NP DO 120 Mz2 N 1 1 DO 120 Lz1 NI1 Ll Nx L 1 M L2 N M 1 P L1 P L1 P L2 2 P P Li 130 IF 1 0 0 TO 150 SUM 0 SUMN 0 DO 140 IL 1 N I JUJ IL 1 N IL SUM SUM DARS P TJ NDL ND6 1L 140 SUMN SUMN DARSOP IJ DIM NDL AL SUMN SUM 150 DO 151 ILz1 N IJz IL 1 XM IL NDL ND6 151 DUM NDL P I ILST I J IF AL LFE 00001 GO TQ 300 IF I GE NF IN 60 TN 310 C INTERMEDIATE PRINT CCI JF I LT NPRINT GO TO 100 NPRINT NPRINT NCDNT 1 155 CALL LNCNT 2 WRITE 6 1152 ILST FORMA IOSTEP NUMBER 1 CALL PRINT 1 GO TO 170 170 170 156 7 156 CALL LNCNT 2 WRITE 6 1152 1 160 CALL PRNT 1 170 IF JFN E0 0 GO TO 100 RETURN 300 CALL LNCNT 2 WRITE 6 1300 1
91. NG ARRAY ONLY NEED FIRST COLUMN JLCT 17 K KD 1 39 D K DFZO 51 OL 1 P L 7803 7801 7803 7803 J JL QL CHECK FOR ROW OF DIAG THAT NEEDS 1 MDVED INTO II IF J EQ IZ GO TO 786 K 1 QS J KD C PUT 1 IN P TO MOVE KyK TO Jy K POSITION 2 PERMUTATIONS P P TRANS THE EFFECT IS TO MOVE K K TO J yK THENCE J J 1 1 00 OL QL 1 C MARK THIS 1 AS ZRO TO BE FILLEO IT IS MOVED UP AND OUTOF HERE 71801 JLIM K M M GO TO 780 i 786 J KD L C MAKE PART OF IDENTITY 786 DON T NEED TO MOVE 1 A HOLE D J 1 00 780 L C RETURN MATRICES COMPLETED E WITH 105 DELIBERATELY LEFT OUT RETURN gt 691 CALL LNCNT 1 WRITE 6 1691 05 0 OST 1691 FORMAT DIMENSION ERROR IN DECOM NCz I14 5X T4 RETURN 692 CALL LNCNT 1 NCC CM CNN CM CC C M Un 1692 FORMAT ERROR IN DECOM 0 7 0 KP 4 QS GO TO 38 END APPENDIX USE OF VASP ON AMES TSS NONCONVERSATIONAL BATCH RJE In using VASP on TSS the system must be told about the job library in which the VASP subroutines are located the source of input data and the location to send output data and the block data program must be loaded procedure has been written for doing this automatically The call to the procedure is VASP input data set output data set The procedure wil then perform th
92. OMPUTING DONE waw 29 recmpt 30 DATA _ 31 1 11 1 22 1 55 32 f h 22 11 2 22 2 55 33 f 7 3 11 3 22 5 55 _ ntas 7 35 M Figure 9 Example of conversational VASP D 36 F MATRIX 3 ROWS 5 COLUMNS 37 1 1100000D 00 2 1100000D 00 5 11000000 00 38 1 22000000 00 2 2200000D 00 3 2200000D 00 397 71 5300000D 00 2 55000000 00 3 3300000N 00 40 MATRIX 3 ROWS _ 41 1 503313D 02 2 68666110 02 2 8799809 02 42 1 5705153D 02 2 8346117D 02 1 07870809 02 bid 1 6476893D 02 2 9625623D 02 14 287435202 02 K COMPUTTN G DONE 46 FORTRAN STATEMENTS 47 0000100 call mult a na x nx y ny 48 0000200 call prnt Uy D C9 0000300 revise 200 50 0000200 call prnt y ny y 1 51 0000300_insert 159 2 0000150 print 6 f 93 1 54 w xMNPG NOW _ ee em 96 nx 3 1 x 1 0 0 0 57 1 110000D 00 1 2200000 00 1 530000D 00 2 110000D 00 2 220000D 00 2 330000D 00 58 371100000 00 3 220000D 00 3 330000D 00 3 200000D 00 3 300000D 00 3 500000D 00 99 1000000 00 4 200000D 00 4 300000D 00 4 400000D 00 9060 Y MATRIX 3 ROWS 1 COLUMNS 61 1 5033413D 02 DOLOR 63 1 6 76893D 02
93. RECISION COMMON X MAXRC NC 1 NA 1 NC 2 NB 2 IF NAC2 NE NB 1 GO TO 999 NAR NA 1 NAC NA 2 NBC 2 NAR NAC NBB NAR NRC IF NAR LT 1 0R NAA LT 1 0R NAA GI MAXRC UR NBB LT T1 0R 1 NBB GT MAXRC GO TO 999 0 DO 300 K 1 NBC IK IK NAC SNAR se gs es TR IR 1 TB IK LI J NAR C IR O DO 300 1 JI JI NAR IB TB 1 C IR C IR A JI B IB 300 CONTINUE 60 1000 999 CALL LNCNT 1 WRITE 6 500 1 1 1 2 1 1 1 2 500 FORMAT 1 DIMENSION ERROR IN MULT wNAsS 216 5X NBz7216 CALL ASPERR 1000 RETURN FND ________ a Y 66 SUBROUTINE SCALE A NA NB 5 DIMENSION A 1 B 1 2 NB 2 COMMON MAX MAXRC DOUBLE PRECISION A B 5 NB 1 NA 1 NB 2 NA 2 L 1 2 IF 1 17 1 11 1 1 1 300 300_B 1 A 1 S 1000 RETURN 999 CALL LNCNT 1 WRITE 6 50 GN TN 999 50 FORMAT DIMENSION ERROR IN SCALE NA 12I6 CALL ASPERR RETURN END 001 SUBROUTINE DIMENSION A 1 B 1 NA 2 NB 2 DOUBLE PRECISION COMMON NB 1 NA 2 NB 2 NA 1 NR NA 1 NC NA 2 NR NC IF NR 1 1 16 Gn TO 999 IRzO DO 300 1 1 j 1 NR DO 300 J 1 NC IJ IJ NR IR IR 1 300 B IR A TJ RETURN 99
94. T PROGRAM CASE 1 FROM ASP MANUAL PAGE 146 VASP PROGRAM B MATRIX 3 ROMS 4 COLUMNS 4 00009009 090 3 99900009 00 3 00900009 00 2 999090 99 2 00000000 00 5 00000000 00 1 00000000 3 00000000 00 2 00000000 1 30000000 01 9 00000000 00 5 00000000 00 APSE MATRIX 4 ROWS 3 COLUMNS 9 50296970 02 5 65801810 02 2 0318850D 02 232 049081720 62 2 gt 21352550 02 L 924 320 D 0 2 64 736480250 02 3 1884964D 02 3 9074711 02 5 06408250 02 3 61302280 02 8 90903410 03 F9 5 e _ 47 4 AASA MATRIX 3 ROWS 4 COLUMNS 4 0000000D 00 1 00000000 00 3 00000000 00 2 00000000 00 240000000D 00 5 0000000D 00 1 0000000 D 00 3 0000000D 00 2 00000000 00 1 3000000D 01 9 00000000 00 5 00000000 00 AERR MATRIX 3 ROWS 4 COLUMNS 64 6613381D 16 6 6613381D 16 2 22044600 16 6 6613381D 16 1 33226760 15 0 0 8 88178420 16 6 66133810 16 2 2204460D 15 3 996802 9D 15 4 2 1 88475D 15 2646613381016 tm ASAA MATRIX 4 ROWS 3 COLUMNS 9 50296970 02 5 6580181D 02 2 03188508 2 3 07908720 02 3 31353550 02 3 78243200 02 6 73648020 02 3 18849640 02 3 90747110 02 5 05408250 02 3 6730228D 02 8 9090341D 03 AERR MATRIX 4 ROWS 3 COLUMNS 3877788D 17 1 21430640 17 344694470D 18 2 3418767D 17 1 0408341D 17 1 64798730 17 0 0 8 67361740 19 1 47451500 17 1 8214596D0 17 2 8062556D 18 gt 4 3368087
95. These subprograms include operations of matrix algebra computation of the exponential of a matrix and its convolution integral solution of the matrix Ricatti equation and computation of dynamical response of a high order system Since VASP is programmed in Fortran the user has at his disposal not only the VASP subprograms but all Fortran built in functions and any other programs written in the Fortran language the storage allocation is controlled by the user so the largest system that the program can handle is limited only by the size of the computer the complexity of the problem and the ingenuity of the user No accuracy was lost in converting the original machine language program to Fortran The principal part of this report contains a VASP dictionary and some example problems The dictionary contains a description of each subroutine and instructions on its use The example problems give the user a better perspective on the use of VASP for solving problems in modern control theory These example problems include dynamical response optimal control gain solution of the sampled data matrix Ricatti equation matrix decomposition and pseudo inverse of a matrix Listings of all subroutines are also included The VASP program has been further adapted to run in the conversational mode on the Ames 360 67 computer The necessary procedures are given in appendix C 17 Key Words Suggested by Author s 18 Dist
96. UTINE LNCNT N C OMM DN L TNES NLP L IN 23 LINZLIN N IF LIN LE NLP Tn 20 WRITE 6 1010 TITLF I T21 23 1010 FORMAT 1H1 23A4 LTN 2 N _ 20 RETURN CN 96 SUBROUTINE ADD A 9 E eC 1 9 N 2 yl COMMON MAX MA DOUBLE PRECISION IF NA 1 1 2 2 GO 999 NC 1 NA 1 NC 2 zNA 2 1 2 IF NACL LT 1 DOR L LI eLeORebLeGleMAXRC GN 10 999 NO 300 I 1 L 30 GO 1000 999 CALL LNCNT 1 WRITE 64 50 NA NP 50 FORMAT DIMENSION ERROR IN NAZ 216 5X NR 216 CALL ASPERR 1000 RETURN 6 SUBROUTINE 5 DIMENSION 1 1 1 2 NBC2 4 NC 2 DOUBLE PRECISION COMMON MAX MAXRC IF M 1 NE NB 1 NAC2 NE NR 2 GO TO 999 1 1 NC 2 NA 2 L NA 1 NA 2 e e le 9 e UR e 7 7 299 DO 300 300 1 1 1 60 1000 999 CALL LNCNT 1 WRITE 6 50 50 FORMAT DIMENSTON ERROR IN SUBT 216 5 216 CALL ASPERR RF TURN CC 86 SUBROUTINE DIMENSION A 1 5 B 1 CU 3 NAL2 NBEU2 4 NC 2 DOUBLE P
97. X05 C FRR EXIT IF ARSURD DIMENSIONS IF OS LT 2 ONT GT MAXRC GO TO 69 C SET COLUMN SIZE RANK SIZE KP 4 QS OL KP 1 C SET SPECIAL PARAMS FOR PSEU CALL THESE ARE TO SUPPRESS THE WORK OF C RF INVERTING PSEUDO INVERSE IN THE CASE WHERE A SINGULAR KP 1 KP 1 1 C CALL PSEU P 10 GET MATRIX T IN C CNOTE THE LAST 3 MATRICES THE 5 IN D USED ONLY IF PSEUP ATTERATESA KP 2 CALL PSEU D KP 1 CL IF I IZ GO TO 38 C PLEASE DD NOT TRY 5 NAMES FOR MATRICES HERE SUCH MATRICES WERE RETURNED BY ASP NOR MY ROUTINE C 13 I 1 ONT 13 D I CCT NV 1 0S NV 2 0S ND 1 2 C ND IS PART OF FLAG PIV RETURNED BY INV CALL INV D NV PIV JL C TEST TO SEF IF PIV IS ZERO ERR MATRIX NOT INVERTIBLE ND1 IF ND 1 IZ 692 28 CONTINUE _ s F rn ms sOF TT UT P IS TO HOLD PERMUTATION MATRIX SUCH THAT THAT PTR ER C ER HAS ALL ONES MOVED TO EXTREME UPPER LEFT OF DIAGONAL Cx NOW SET UP TO MAKE P PERMUTATION MATRIX P 1 OR OS C ZERO OUT P WILL BE ZEROS AND ONES DO 39 I 1 ONT C ZERO HOUSEKEEPI
98. am Line b This statement reads the first card of the data deck see fig 2 b places its contents in TITLE array and prints the first line of the output see fig 2 c Lines 6 and 7 The matrices F G J R K H X and T are read in from data deck see fig 2 b and are printed Line 8 Line 8 calls the TRNSI subprogram performs the computation and prints outputs as explained in the example Example 3 An Optimum Control Problem Given a system x Fx Gu y Hx x o where x u and y are respectively the state control and observation vectors The system distribution and observation matrices are F G and H respectively We wish to define optimal control u t where u t Kx t so as to minimize the performance index J U RU dt The solution to this problem is K R G P where P is the solution of the matrix Ricatti equation The VASP program finds P by means of the subroutines AUG ETPHI and RICAT as follows First subroutine AUG is used to generate the matrices Z SE d C R G HQH Note This is the negative of the Z given on page 212 of the ASP manual Subroutine ETPHI is then used to compute the special transition matrix 40 _ to 02 022 Finally the matrix is computed by subroutine RICAT for a time increment of by repeated application of the formula P t 7 02 0 0 0 P t The computation is repeated fo
99. amp Data used by Conversational VASP REWRT N Edit SOURCE MNPG EXCISE from Statement N to last List program Figure 10 Flowchart VASP procedures 158 6S1 CHNGIN 0000000 PROCDEF CHNGIN 7 0 777777 CHNGIN 0000100 PARAM 1NPUT CHNGIN 0000200RFLEASE FT05F001 HNG1N 0000300 TE T T PUTY lt DDEF FTOSFOOI SINPUTSDISPLAY TNPUT FROM DATA SET SINPUT CHNGIN 0000500 IF INPUT 5 INPUT FROM TERMINAL M n CHNGOUT 0000000 PROCDEF CHNGOUT CHNGOUT 0000100 PARAM OUTPUT CHNGOUT 0000200 RELEASE FTOGFO0L CHNGOUT 0000300 IF OUTPUT DDEF FTOGFO01 OUTPUT DISPLAY OUTPUT PLACED IN DATA SET OUTPUT CHNGOUT 0000400 IF OUTPUT DISPLAY OUTPUT TO TERMINAL m 9 _ M KO VoO 0000900 LRVASP 0000000 PROCDEF CLRVASP 77 a CLRVASP 0000050 END CLRVASP 0000100 UNLOAD MNPG CLRVASP 0000200 UNLOAD VASPMN i CLRVASP 0000500 ERASE SOURCE VASPMN SOURCE MNPG USERLIB VASPMN CLRVASP 0000400 ERASE USERLIB MNPG 0000500 RELEASE VASPT 7 LOMA CLRVASP 0000600 DISPLAY
100. an easily determine the solution of the matrix Riccati equation for a 30 order system perhaps 40 depending on the size of other related matrices Another basic difference between these two Fortran versions is that VASP diagnostics are more self explanatory To recapitulate the objectives of VASP are flexibility and versatility so that it can serve the maximum number of users To achieve these goals FASP was revised extensively so as to have for example variable dimensioning more efficient storage and more self explanatory diagnostics In this report no attempt will be made to discuss any details of the theoretical background and the algorithms associated with the appropriate subprograms since they are well documented in the ASP manual an NASA contractor report NASA CR 475 1966 This report does not repeat information from the contractor report and the user is urged to consult that manual so that he may utilize VASP proficiently This program can be obtained from the centralized facility known as COSMIC located at the Computer Software Management and Information Center Barrow Hall University of Georgia Athens Georgia 30601 FEATURES OF PROGRAM The advantages of VASP over ASP are 1 more versatile programming language 2 more convenient input output format 3 some new programs and 4 variable dimensioning Since VASP is programmed in Fortran it can be used in a very large class of machines Moreover
101. arguments The initial value of P must be placed in PT before calling the subroutine The value of P is updated every iteration in the subroutine until the final P is reached This final P is of the outputs of the subroutine 21 SAMPL DESCRIPTION Subroutine SAMPL calculates the covariance and weighting matrices associated with the discrete case of either the control problem or the filter problem Consider the following filter problem Given the system where gaussian random sequence with variance Q and observations Hx v where v gaussian random variable with variance R The optimum estimate of the state is see p 234 in the ASP manual Ox Hx where K 9P H I HP HT P P P HT HP HT R T HP LI Pi Q pseudo inverse Here is the solution of the matrix Ricatti equation which is obtained by SAMPL The subrou tine has for inputs 0 H Q R and for output and K _ where Pus is written over Pi REMARKS 1 The routine will take n steps at a single call where n is an input parameter Further if P becomes constant then the routine will stop and exit before completing the n steps The actual test is as follows 24 If 1 Pa ae 105 then exit k k 2 The routine will print the value of Pi and or K every j steps and also when either exit occurs NCONT 3 controls which arrays are prin
102. as VASP is part of a larger main program all the Fortran built in functions are available to the main program Furthermore any subroutine available in the Fortran language may be used In other words the user has at his disposal the combined capabilities of VASP Fortran built in functions and all other subroutines written in Fortran The input output subroutines have been changed and now consist of READ RDTITL and PRNT In addition LNCNT has been added to control paging The VASP allows the user the option of using the existing standard VASP format or of supplying the output format of his own choice For a more detailed explanation of how to exercise this latter option see the dictionary entry under PRNT p 10 or Example 2 Our experience with ASP is that certain groups of statements are often repeated For the user s convenience these groups of statements are incorporated as new subroutines in the VASP They are AUG UNITY and SCALE Detailed explanations of them are available in the VASP dictionary in this report To utilize the storage space as efficiently as possible the subroutines are written with variable dimensioning with the storage allocation controlled by the user s calling program Consequently it is necessary to provide some dummy storage space as a part of the argument of the subroutine From the user s point of view the price for efficient storage is inconvenience All the subroutines are written in double precision
103. at ORTH leaves T in DUM3 Also if the T at input was the null matrix the rotation will be the identity matrix so that R ZL Lines 19 through 27 then juxtapose either the matrix EXR or the matrix EXC using JUXTR or JUXTC depending on the compatibility of the dimensions If both sets of dimensions are incompatible no juxtaposition is done In any case the result of this operation is placed R The decomposition routine is called next the original ZL matrix had zero in element 2 1 and no juxtaposition was done then R is assumed symmetric and the DECSYM entry is used If ZL was not symmetric the program will produce errors Otherwise the DECGEN entry is used lines 29 to 31 Finally the resulting matrices H and G are tested using R HG RE R R and all resulting matrices are printed In figure 5 c blank lines represent blank cards In the data cards for case 4 the header card for EXR has no dimension information and no associated data cards This indicates that the matrix EXR is to be left unchanged and that no data cards are to be read for EXR In case 7 EXR is again left unchanged A blank data card follows the EXC header card The output figs 5 d through 5 f contains the results of decomposing three different matrices Figure 5 d case 1 isa 2 X 2 rank 1 matrix figure 5 e case 4 is a 2 X 3 rank 2 matrix 53 vs C MAIN PROGRAM TO CHECK DECOM ET AL 0001 DIMENSION NZL 2 2 4
104. atical description and examples are given in NASA CR 475 Subroutine PSEU calls ANDRA to do each pivoting step after first forming a symmetric matrix B which is indeed forced to be perfectly symmetric The first call of ANDRA is an initialization call An identity matrix T is formed The rank counter is set to zero On an initialization call the routine proceeds to search the diagonal for pivots as always But after finding a pivot it always goes back and looks for another pivot regardless of the iteration option The process of searching for pivots continues until the number of tries is one greater than the row size no such test is made in the iteration case If the routine fails to find a single pivot in the initialization call it exits with an error message Pivots are accepted if and only if they are not less than a threshold input at every call Supposing that a pivot has been found in the diagonal the next step is always the same First the pivot is reduced to unity That is both the pivot row and column are divided by the square root of the pivot in B Only the row of T is so reduced The next step is to eliminate the pivot coefficient from all other rows not yet used as pivots This part is the same as in other inversion methods Both are treated exactly alike here Note that the actual algorithm checks the diagonal of a row to see if it is already marked asa pivot If so that entire row and the row in T are skipped The piv
105. e only a portion of the program from line 150 on Accordingly the REWRT command is issued with a parameter line 65 The system then erases SOURCE MNPG from line 150 inclusive to the end It then lists that portion of the program being used in this case line 100 only line 66 and prompts the user for additional lines with a line number line 67 The user then adds lines as desired lines 67 through 69 and requests a compile line 70 It can be seen that line 69 is missing a left parenthesis so the compiler prints a diagnostic and requests the line be corrected lines 72 and 73 The correction is entered line 74 after which the compilation is completed lines 75 through 77 No data are needed so the data request line 78 is answered with only line 79 The results are printed on lines 80 through 88 Since no more computations were desired a Logoff command was issued lines 89 and 90 Housecleaning A procedure called CLRVASP is available This procedure erases all data sets that have been set up by the various other procedures and allows the user to keep his storage low Use of the routine is not required since the other procedures have appropriate erase statements as needed 156 LISTINGS AND FLOWCHARTS Figure 10 shows all the procedures associated with VASP and indicates what each one does A complete listing of the procedures is given in figure 11 Figure 12 is a listing of data set VASPPROC If the user executes
106. e steps indicated above If the first parameter 15 omitted the data will be taken from SYSIN which is from cards in your data deck If an input data set is named then the data will be taken from the named data set which must have been stored previously Likewise if the second parameter is omitted the output will be placed in SYSOUT for printing on the high speed printer If an output data set is named the output will be placed in that data set If the name of the input or output data set must be changed use the procedure call CHNGIN new input data set name CHNGOUT new output data set name These two procedures will then change the DDEF to the new data set name If the parameter is omitted the new data set name will be SYSIN or SYSOUT A listing of these procedures is included in this appendix CONVERSATIONAL Provisions have also been made to allow conversational use of the VASP program so that the user can easily perform matrix operations The operations can be strung together in a sequence as desired with as much output as desired The user indicates the operations by use of Fortran statements and may not only call the VASP subroutines but also may execute any other Fortran statements that he wishes Data are requested for the program by means of subroutine INPUT allowing free form data from the typewriter If Fortran type input is used the data should also be obtained from the 151 typewriter If you try to use an input
107. e style as in line 17 with the elements of the matrices 152 om 1 1 vasp DDNAME JBLB0001 73 INPUT FROM TERMINAL 4 OUTPUT TO TERMINAL 5 convasp 16 f z fstar 8 5 AVAILABLE ALL DIMENSIONED 16 ARE A F C W X Y FSTAR FOR INPUT OR COMPUTATIONS D1 D2 03 D4 D5 D6 D7 __ FOR COMPUTATIONS ONLY _ _ 9 w FORTR N STATEMENTS 10 0000100 t 1 0 11 0000200 call etphi f nf t a na m TZ 0000300 call prnt E 13 0000500 call prnt a na a 14 0000500 1 T5 ww w wMNPG NOW COMPILING w O _ n ee U JU L ae bu aal q _ _ 16 DATA 17 6 1 1 1 2 1 35 1 2 1 2 2 2 5 2 3 1 3 2 35 3 2 1 2 6 5 18 F MATRIX 4 ROWS COLUMNS 19 1 10000000 00 2 1000000D 00 3 1000000D 00 10000000 00 20 1 2000000D 00 2 20000000 00 5 20000000 00 4 amp 2000000D 00 21 1 3000000D 00 2 3000000D 00 3 3000000D 00 4 3000000D 00 _22__1 10000000 00 2 40000000 00 3 4000000D CO 4 50000000 00 23 A MATRIX ROWS COLUMNS 24 7 72516 70 03 1 56901660 04 1 9656167D 0h _ 2 56221680 04 _ 25 8 0162773D 03 1 2088250 0h 2 0399372D 0h 2 665909190 26 8 3083898D 03 1 5725583D 221143577 Ok 2 7559671D Ok 27 8 6005025D 05 1 5243142D OF 2 18857820 Qh 2 85291220 Oh 798 C
108. e subroutine and are available to the calling program upon return When a dummy array is required it must be appropriately dimensioned in the calling program The required size is given in the appropriate dictionary entry Most of the routines check the array dimensions for compatibility and reasonableness and check for adequate storage available in the DUMMY array The reasonableness test is to see that all matrix dimensions are greater than zero and that the product of the matrix dimensions is less than the constant MAXRC In the program MAXRC has been set at 6400 It is recommended however that the user reset MAXRC to equal the size of his matrices and thus prevent accidental overflowing of the matrix dimensions If the matrices are incompatible or unreasonable or if a mathematical error is found a self explanatory error message is printed and no further computations are made in that subroutine However computation does go on to the subsequent steps which are likely to be wrong also After 10 such errors the program is terminated The VASP program uses double precision arrays so that the user s main program must define each matrix to be double precision and to contain a sufficient number of cells to accommodate the matrix The dimension statement may classify the array as one or two dimensional or even more For example to use the matrix A which is a 5 X 5 matrix any of the following dimension statements will be adequate
109. ension array NA Scalar TT where TT is the value of t used in equations 18 Output Arguments Matrices B C Dimension arrays Dummy Arguments Matrix DUMMY Constant KDUM Note KDUM contains the size of the DUMMY matrix which must be at least 2 NA 1 NA 2 18 ETPHI DESCRIPTION This subroutine computes the matrix exponential p e t See ASP manual page 92 and also EAT page 18 of this manual for reference USAGE CALL ETPHI A NA TT B NB DUMMY KDUM Input Arguments Matrix A Dimension array NA Scalar TT where TT is the final value of time Output Arguments Matrix B Dimension array NB Dummy Arguments Matrix DUMMY Constant KDUM Note KDUM contains the size of the DUMMY matrix which must be at least 2 NA 1 NA 2 19 19 AUG DESCRIPTION This subroutine computes C R G and F GR G HOH F The matrices C and eZt then used in RICAT to calculate the covariance and weighting matrices These matrices arise from a linear system of the form x Fx Gu with output equation y Hx and cost function J u Ru dt See ASP manual page 212 for reference In the special case where then F TrGR G Q F and the cost function is J I control index II is used to distinguish the two cases 20 REMARKS The inputs to this program are the matrices F G RI H Q F must be square b must
110. eparate statements as in the BLKDATA program If they are loaded in one statement they will probably load incorrectly because of the dimensionality of FMT2 Also NEPR must be consistent with the numbers in FMT1 and FMT2 USAGE CALL PRNT A NA NAM IP Input Arguments Matrix A Dimension array NA Matrix name NAM If NAM 0 a blank name will be printed NAM should contain four hollerith characters It can be written in the calling sequence as 4HAbbb If written the last three characters of the printed name will be garbage Control constant IP 1 heading same page 2 heading new page 3 no heading same page 4 no heading new page Output Arguments None 4 LNCNT DESCRIPTION This subroutine keeps track of the number of lines printed and automatically pages the output as required It is completely internal and the user need not refer to it unless he has WRITE state ments of his own In that case the user may should put the statement CALL LNCNT N before each WRITE statement where N is the number of lines to be printed Page length is controlled by the variable NLP set in the BLOCK DATA program to 45 This is an installation dependent variable and may be changed as necessary The subroutine provides one line of print at the top of each page This line contains 92 characters of which the first 72 are available for the programmers use and may be loaded by use of RDTITL The remainder contain
111. erating without having found at least one pivot In other words ANDRA always finds all the pivots it can before any side calculations are done If this first rank is maximal it never iterates The first pivots are not in doubt so these rules are more efficient The routine always uses a tolerance for pivot acceptance however it uses a new tolerance 50 000 times smaller than the last pivot found for each call to find one pivot in iterative mode The expensive test of matrix norms is avoided when no new pivot occurs The PSEU routine has only a finite number of tries to find a new pivot before it quits The exact number is the same as the maxi mum rank Since ANDRA has usually found several pivots initially this is ample 73 Iteration If DBP2 is larger than 1 it is raised to a power used as a factor and tends to make the routine stop with a smaller rank DEP2 of 1 actually works for most iterations Subroutine ANDRA The basic algorithm can be used as a separate routine by itself see ANDRA documentation The routine requires considerable setting and testing of parameters It has an escape exit for too many iterations calls to find only one pivot without finding any It returns a matrix T such that if X if pseudoinverse of positive definite matrix A then T T X Accuracy In double precision the accuracy has been very good Maximum accuracy can be obtained by using symmetric matrices and the PSEUP entry The test program included
112. error with all the 175 contiguous starting in the first diagonal If had been so permuted before diagonalizing then this different T resulting would be the one that gives an inverse that corresponds correctly to the old But since one is using a premultiplication and a postmultiplication simple substitution of a permuted matrix does not work It would if matrix multiplication were commutative Thus if it is necessary to transform the original starting B reinvert then transform back again The permuted form of T which does not actually occur has a nonsingular corner submatrix followed by the rest of diagonal set to I s These latter 175 correspond to the dependent equations of the original The rule for constructing the transforming matrix U is given below This matrix is made from T and put into the same storage T The explicit construction of U is more efficient in FORTRAN From here on the explanation concerns what is actually done rather than the mathematical reasons Let d denote the ith diagonal element of B In case the reader has forgotten this has been changed to a diagonal matrix of O s and 1 s Given T there are two cases Case One For Ui not on the diagonal Use tij if d 0 Use 0 if 1 Case Two For U on the diagonal Use the corresponding value of d Next using a copy of the original B form C U BU 77 The result is actually put in the same storage that held B originally The smallest a
113. ful 10 avoiding confusion between ranks determined from different calls Used also to output the effective rank PSEUDO sets IP1 and IP2 to zero The row size of the matrix input The column size of the matrix input Note IP4 need not be specified for PSEUP entry Holds the pseudoinverse output In rank only case holds a diagonal matrix with 0 s and 1 corresponding to pivots accepted or rejected In nonsingular case holds the matrix T of the diagonalization case In singular case holds that certain matrix U described in ASP manual Holds the pseudoinverse of the original B Note and B are the same size The other matrices are square of the size of C which is determined by the smaller dimension of A D is either five times the size of C if iterating or the same size as Matrix D In the nonsingular case D holds a copy of the B formed from A It equals A for a PSEUP entry In the singular case it holds a pseudo inverse for a permuted so that independent variables are all moved to the left most positions Note D has possibly four other matrices Let these be D1 D2 D3 and D4 in order They are used only if iterating D1 also used by DECOM D1 D2 hold old results D3 D4 holds intermedi ate values when doing the side calculations PSEUDO does not provide for D1 through D4 Notes on Usage Symmetry This method is well suited to symmetric non negative definite matrices The PSEUP en
114. in FAP Fortran Assembly Program and could be used only on the IBM 7090 7094 Second many complicated time variant analysis synthesis and optimization problems tax the capability of the ASP and other programs written in the Fortran language Fortran IV language makes the program adaptable to a much wider class of computers and expands its versatility The VASP is based extensively on a Fortran version of ASP called FASP Fortran ASP by its programmer Mr Don Kesler of Northrop Norair Two basic questions the user will inevitably ask are 1 How accurate is VASP compared with ASP 2 What is the highest order of system that VASP can handle The answer to these questions depends on the number of significant digits carried by the user s computer and the amount of available storage in the computer To answer the first question in a more concrete way the check cases given in the ASP manual were duplicated and the results were compared with those in the manual The accuracy of VASP was found to be the same as that of ASP The second question can best be answered by first noting some of the basic differences between FASP and VASP The pertinent difference between the two is that VASP has variable dimensioning and more efficient storage To allow a computer to handle the highest order system possible all matrix storage is assigned by the user s main program Consequently as an illustrative example a 125 000 byte version of the IBM 360 50 c
115. ine 12 Line 12 is the data format For this program the transient output was printed using the programmers write statement rather than PRNT The use of PRNT for this purpose is shown in the third example p 40 Line 13 Line 13 tells the line counter that the program will print 4 lines Line 14 Line 14 does the actual printing Lines 15 through 25 Lines 15 through 25 form a loop which increments TT line 23 and stops when is large enough line 24 Line 15 Line 15 computes the B and C matrices for time TT When C is computed the limits of the integral are and the present TT Note that W is specified for dummy storage and the 18 tells EAT the size of W 33 Line 16 Line 16 computes BXO stores the result in VI Array VI isset up for the programmers working storage Since W isalso available at this point in the program it could have been used instead of V1 if desired Line 17 Line 17 computes CG and stores the result m another working storage array Line 18 Line 18 computes CG U and stores the result in V2 still another working storage array Note that MULT obtains the product CG from 1 Line 19 Line 19 adds V1 and V2 to obtain XT Since the ADD subroutine allows the matrices to be repeated in the call the array V 1 could have been eliminated then line 16 would have stored its results in XT Line 19 then would have added XT and V2 to obtain the complete XT Line 20 Line 20 multiplies H times
116. initial condition for this differential equation is P t H S T H The optimal control in terms of the state variable x t 15 u t K t x t G P t x t AUG computes Z rather than Z so that the exponentiation for 0 uses positive time increments 22 and the optimal feedback gain K t is K t R 1G P t Letting C RIG then K t CP t REMARKS 1 This subroutine will be terminated when n n gt Pi lt where 10 1 1 1 1 or NC NT 2 steps have 1 2 Matrices P t and K t will be printed out every NCONT 1 steps as controlled by NCONT 3 3 Matrices 9 9 611 are submatrices of 0 Their dimensions are n X n where n is the order of the system i e the dimension of the F matrix They are partitioned from the 0 matrix as follows The Fortran symbol for 0 is PHI USAGE CALL RICAT PHI NPHI C NC NCONT K NK PT NPT DUM KDUM Input Arguments Matrices PHI C PT Dimension arrays NPHI NC NPT Control array NCONT 1 Number of steps per print NCONT 2 The maximum number of steps NCONT 3 Printout control 1 gt no P no K 2 P only 3 gt K only 4 P and K 23 Output Arguments Matrices K PT Dimension array NK Dummy Arguments Matrix DUM Constant KDUM Note KDUM contains the size of the DUMMY matrix which must be at least NPHI 1 2 Note PT is used for both input and output
117. inverse of the input matrix It sets up a standard set of options for use by PSEU which does the actual inversion For details of the method see PSEU p 70 26 USAGE CALL PSEUDO A NA B NB DUM KDUM Input Arguments Matrix A Dimension array NA Output Arguments Matrix B Dimension array NB Dummy Arguments Matrix DUM Constant KDUM Note KDUM contains the size of the dummy matrix which must be at least 3 NA 1 NA 2 24 DECGEN 24a DECSYM DESCRIPTION This subroutine decomposes a real matrix R with dimensions m X n and rank r lt min m n into two matrices H and G such that R HG Further both H and G are of maximal rank with dimensions m X r andr X n respectively It uses subroutine DECOM to provide matrices from which H and G can be computed The writeup of p 85 describes the method in detail Subroutine DECOM requires for input a matrix which is positive semidefinite symmetric Sub routine DECGEN computes this matrix by letting RRT or whichever is smaller and uses the former if R is square If the user knows that R is already positive semidefinite symmetric this step may be omitted by a call to DECSYM in which case A USAGE CALL DECGEN R NR G NG H NH DUM KDUM if R is general or CALL DECSYM R NR G NG H NH DUM KDUM if R is positive semidefinite symmetric 27 Input Arguments Matrix R Dimension array NR Output Arguments Matrices H G Dimension arrays
118. ired to generate the transient response of such a system in response to a given initial condition and fix control u In particular given I O ON w C I es find x t for 0 lt 1 lt 2 0 Also print x t and y t every 0 01 second The user s main program to solve this problem is shown in figure l a the corresponding data deck is shown in figure 1 b where each line represents one card and the beginning of the results is presented in figure 1 c 29 0 0001 DIMENSION 3 3 2 6 3 3 6 2 2 32 2 3 3 4 5 nd 1 2 0002 DOUBLE PRECISION F X0 U XT YT TT 11 9 0003 MAXRC O00424 COMMON LINES NLP LINyTITLE 23 0005 MAXRC 9 0006 CALL RDTITL 0 00 0008 READ 5 100 TTyDELTATyTFINAL 0009 100 FORMAT 8F10 2 OO1O ___ CALL READ 5 E NE G NG yH NUS XO NXO 5 Q 0011 101 FORMAT 1HO 59X TIME 5 5 22 5 54 1 OUTPUT TT 14X XT 1 11X XT 2 43 11X 3 21 9 1 1 a tS 0012 102 FORMAT 1HO 10 2 6 3 16 7 10 3 16 7 0013 CALL LNCNT 4 rigo NOD LA
119. is either obvious or is given the ASP manual NASA CR 475 Other routines were written by one of the authors These were quite simple and needed little description Subroutines ANDRA BDNRM DECOM and PSEU were written by John Andrews Information Systems Company Since no description of these subroutines is available elsewhere a detailed description of their theory and usage is included Because there were various programmers the nomenclature internal to the various subroutines is not completely consistent 6 gt to S 10 CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL READ LA NA B NB C NC D ND E NE RDTITL PRNT A NA NAM IP LNCNT N ADD A NA B NB C NC SUBT A NA B NB C NC MULT A NA B NB C NC SCALE A NA B NB S TRANP A NA B NB INV A NA DET DUM UNITY A NA TRCE A NA TR EQUATE A NA B NB JUXTC A NA B NB C NC TABLE 1 SUBROUTINE CALL STATEMENTS IN VASP CALL CALL CALL CALL CALL CALL CALL CALL CALL CALL JUXTR A NA B NB C NC EAT A NA TT B NB C NC DUMMY KDUM ETPHI A NA TT B NB DUMMY KDUM NF G NG RI NRI H NH Q NQ C NC Z NZ JI RICAT PHI NPHI C NC NCONT K NK PT NPT DUM KDUM SAMPL PHI NPHL H NH Q NQ R NR P NP K NK NCONT DUM KDUM TRNSKF NF G NG J NJ R NR K NK H NH X NX T DUMMY KDUM PSEUDO A NA B NB DUM KDUM DECGEN R NR G NG H NH DUM KDUM DECSYM R NR G NG H NH DUM
120. is not larger then a new pivot is searched for If not PSEU reverts to the previous values before the current pivot was used practice a number of modifications are made First the pivot used last js returned as DEP3 even if rejected so that the user can reconsider acceptance of it Second if maximum rank is achieved prior to iteration no side calculations are done Third the smallest possible pivot allow able is set to 0 00002 times the most recent pivot in order to reject many spurious pivots without doing the lengthy side calculations This modification is based on actual observation of pivot behavior The successive pivots of an ill conditioned matrix usually decrease fairly rapidly But there is usually a hugh jump in order of magnitude between the last good pivot and the first bad one Parts of the side calculations are actually done in single precision to save time Please note that a single iteration besides the ANDRA call makes ten subroutine calls and one library routine call Naturally this is slow 78 Matrix Storage Flow This section uses the same names as the Fortran IV routines It tells what is put into each matrix of PSEU at various times The call is CALL PSEU A B C EE DEP IP D The matrices and B are the same size possibly nonsquare Matrix C is square with dimension equal to the smaller dimension of A The other matrices are the same size as C Matrix D is divided into five matrices Let these be denoted a
121. llowable pivot for ANDRA is recalculated This result C is sent to ANDRA to do the diagonalization again The fact that C has rows and columns of 055 that ANDRA has to skip makes the diagonalization ineffi cient but this cannot be helped No iteration is done here Let T denote the result of this second ANDRA call Then the new pseudoinverse is X Tr Transform this back to get a correct answer The rest of the computation is as usual Note that if the rank changes the second ANDRA call error exits are taken Iteration The main method itself is purely algebraic The iteration option is a way of estimating the amount of error in the generalized inverse and using this to stop with a smaller effective rank Let B denote a matrix and X its pseudoinverse after taking so many pivot steps in ANDRA Then the two Moore Penrose axioms read BXB B XBX X If the iteration mode is selected ANDRA first finds all the pivots it can Then subroutine BDNRM is called twice Each call returns the value norm Q P Q Q norm Q The values of P and Q are B and X in one call X and B in the other The resulting two small scalars which would be zero if the axioms were perfectly satisfied are added together The result is taken as a factor times DEP2 raised to the current number of pivots From successive iterations one obtains a sequence of positive numbers decreasing as one approaches the largest possible rank As long as the new result
122. matrices and the sign controls multiplication procedure Output Arguments None This is a function subroutine Dummy Arguments Matrix D dummy array of size 5 NR Constant Array designates location of submatrices of D KRV1 NR KRV2 2 NR KRV3 3 NR KRV4 4 NR 30 ANDRA SUMMARY ANDRA is a Fortran routine to diagonalize a positive definite symmetric matrix The routine was originally designed to be used by subroutine PSEU The routine has a parameter to command it to initialize on the first call Two different modes can be used for pivoting steps In the first mode the routine does only one pivot to eliminate only one row at a time In the second mode as many pivots as possible are done in one call Pivots are chosen in order of decreasing magnitude They are rejected if smaller than a parameter threshold The original matrix input is destroyed and 82 replaced with artificial values However symmetry is after each pivot The answer matrix is such that if X isthe inverse of the input X T T The routine has error exits for matrices of the wrong size and for those that allow no pivot on the first try USAGE CALL ANDRA B T DPR JP Name B DPR 1 DPR2 JP JP1 JP2 JP3 JP4 5 Description Input symmetric matrix Destroyed Answer inverse of Parameter array of size 2 DPR1 is the tolerance for trial pivots Any less than this are rejected as zer
123. matrices must be used In the ordinary case matrix D may be omitted and also matric E is not used Naturally no pseudoinverse is returned 79 08 Set Values Pt Parameters Dri 002 Compute Qut 5 26 Input Array 2 691 Or Too Large Error P d No S Flags For Heration Set Count Of tterations 9 Yes Syouneine No 2 Yes Na MC 18 Start PSEUP QPS 0 Column Sizy Sae Farm Witr Column Sue Do B Matrix Computation Form With Row Suc Uso Symmigtric Matrix Already taput 198 Take Average To Make B Syri metre Set QNT NR x Nf 549 of B Sot Aftsots For D Matrices If flank Only Skip Save a Copy 6 tn D For Singular Qr tter E wt Up For Andra I Find Largest Qiagons Bement Sat Smallest Non Zero Pivot DEMI x Max Diag Get JP Purgineters AM Possible Pivots Andra Sord DEP 2 First Time Htecotion Flag Yes Save Old B And C Arrays 0102 Call ANDRAIB JP Sat JP To Value Heration Flag IR 7 Rank Returned Yes Omit Side Date Rastore Values Before Current Step Reduce Rank By One CaH ts There No Herati n Or Is Rank On Ist ANORA Sime Pivot Ar Before Accept Hank Gwen Yos 720 No Cali Twice Store Retuit SS
124. mmy storage and is dimensioned 18 as required by the EAT subroutine Lines 3 and 4 Common variables to be needed later are made available to the program Although the variables listed in line 4 are not needed in this program they are shown 101 reference Line 5 Since the basic matrices are 3 3 is set to 9 to prevent overfilling the matrices Note this will not protect from overfilling the arrays XO XT etc since they are expected to be 3 X 1 vectors and are dimensioned 3 Line 6 This statement reads the first card of the data deck see fig 1 b places its contents in the TITLE array and prints the first line of the output see fig 1 c Lines 8 and 9 The initial time the time increment and final time are read from the second data card Line 10 The arrays F G H U and XO are read from the remainder of the data deck and are printed fig 1 c Note that the dimensions used by the program are those given on the header card for each matrix If these were specified as 2 2 this same main program would solve a second order problem rather than the third order problem If the initial conditions were already stored in the XO array and you did not wish to disturb them then the header card for the XO array would contain only the matrix title no dimen sions and the associated data cards would be omitted The matrix XO would still be printed Line 11 Line 11 contains the information to head the main output L
125. mposition and pseudo inverse of a matrix Listings of all subroutines are also included The VASP program has been further adapted to run in the conversational mode on the Ames 360 67 computer The necessary procedures given in appendix INTRODUCTION VASP the Variable dimension Fortran version of the Automatic Synthesis Program is the new Fortran IV version of the ASP the Automatic Synthesis Program It is intended to implement the Kalman filtering and contro theory Basically it consists of 31 subprograms for solving most modern control problems in linear time variant or time invariant control systems These subpro grams include operations of matrix algebra computation of the exponential of a matrix and its con volution integral and the solution of the matrix Riccati equation The user calls these subprograms by means of a FORTRAN main program and so can easily obtain solutions to most general prob lems of extremization of a quadratic functional of the state of the linear dynamical system Particularly these problems include the synthesis of the Kalman filter gains and of the optimal feedback gains for minimization of a quadratic performance index The VASP is an outgrowth of ASP which was developed for NASA under contract with the Research Institute for Advanced Studies a division of the Martin Company There are two urgent reasons for reprogramming ASP into the present Fortran version First ASP was programmed
126. nary K is the weighting matrix corresponding to the beginning of the interval and P is the covariance matrix corresponding to the end of the interval This is appar ent in the output For example the first entry to SAMPL prints step number 0 and the K matrix followed by step number 1 and the P matrix On exit from SAMPL P and K contain the data corresponding to Pi and Ki which is the last interval If printing is requested the exit value of and K will always be printed and will be the last set of data 49 09 C CHECK PROCEDURE FOR SAMPLE SEE PAGES 234 AND 244 ASP MANUAL 0001 DIMENSION NPHI 2 NQ 2 NR 2 NP 2 NK 2 3 NH 2 2 1 IM L a L 1 DUM 54 0003 COMMON MAX MAXRC 0004 MAXRE 9 0005 NDUM 54 0006 RB 0008 10 CALL READ 4 PHI NPHI H NH 0 NQ R NR R NR 0009 NP 1 NPHI 1 90 9 NP 2 ENP H H 24 T 0011 CALL UNITY P NP 0012 NCONT 1 1 Q9 3 3 NC NF a 0014 NCONT 3 24 0015 CALL SAMPL PHI NPHI H NH Qs NQ R NR IP NP K NK NCONT DUM NDUM 0016 0017 GD TD 10 0018 END a Main program TEST PROGRAM FOR SAMPL CASE 1 FROM ASP MANUAL P234 AND P244 PHI 3 3 0 1 0 0 0 0 0 0 2 0 2 3 0 0 2 0 0 0 0 0
127. o DPR2 is the last pivot actually used Unchanged if no new pivot found Integer parameter array of size 5 Zero if all pivoting to be done on one call nonzero if only one pivot per call Zero if initialization call Subroutine sets to one when a pivot is found Holds the effective rank number of pivots found The integer giving the row and column size May range from one to a nominal figure The integer row where the last pivot was found The rows are left in the same positions as in the input matrix 83 v8 Start ANDRA EF Single precision Version of DPR1 FMX 0 Zero Maximum Diagonal Element And Index Of Max 1 0 Test Part Of DDI 1DD For Pivot Code IDD tec Skip Old Pivot No NOTE Do not use reciprocal for single precision Yes Not An Initialization JP2 0 QS JP4 Number of Rows QS x OS Size of Matrix INITIALIZE FORM IDENTITY MATRIX T Figure 8 Information Systems Co flow chart subroutine ANDRA B T DPR JP Error Initialize L Diagonal Index No I I 1 Update Index L To Current Diag DDI BiL Ool ES Test Diagonal FDI lt FMX Choose New Max DMX DDI Save Index MI L M 70 DRS 1 Initialize Indices To Start Of Pivot Row Pivot Column K Pivot Row L Pivot Column Do Loop 41 1 15 Index of Row DDM B L x DRS T L T L
128. oore Penrose generalized inverse of a non negative definite double precision matrix It has a separate entry PSEUP for input of a matrix that is already symmetric symmetric matrix is always used for the actual diagonalization process This 70 process is done a self contained subroutine ANDRA The routine never fails since it includes the singular case However it may fail to give the correct rank To control this an option to do side calculations is available After the first pivots have been found if the rank is not maximum the result of each pivot step is used in two axiomatic expressions subroutine BDNRM This side cal culation yields a measure of the worth of the pseudoinverse obtained so far This result is multi plied by a parameter factor raised to the power of the current rank nonlinear penalty function The routine can backtrack from the first bad step and stop with the previous rank IL hasan option to do the minimum calculations for getting a rank only The generalized inverse is useful for least squares solutions of Ax b it works when is singular This method is best suited to symmetric matrices The routine has suitable error exits USAGE CALL PSEU A B C EE DEP IP D or CALL PSEUP A B C EE DEP ID D Note PSEUDO uses PSEU entry Input Arguments Description Matrix A The array to be inverted left intact must be symmetric if PSEUP call is used Non negative definite or nearly so Control array
129. ot is then marked by an artifi cial code The routine always tests for this code and does not use this row again The code is put in the actual pivot position Thus the rows and columns are left in their starting places in the working matrix PSEU converts the result to a matrix of 15 and 075 that shows the independent and dependent variables The code is tested for an integer This is a considerable economy The resulting T is never singular If B were nonsingular and X the desired inverse of B T T X 76 This part is done by subroutine entry TTRM using coding shared with the iteration method The final answer is put back in matrix B PSEU always uses the original A again rather than the origi nal after this to give an answer for A Thus ANDRA is always supplied with a symmetric matrix B If B were singular at the start further reinversion would have to be done See the next section The Singular Case Suppose that the rank of in the diagonalization by ANDRA does not turn out to be maximal then PSEU must perform a number of matrix multiplications and call ANDRA TTRM to reinvert The accuracy is bound to suffer even though the reinversion is done on an exact copy of the original B A very short justification is given belov followed by a close description of how the work is actually done There exists a permutation matrix P such that is a matrix of 075 17 were it not for round off
130. ould arrange to turn off the underflow error messages THE VASP DICTIONARY A detailed description of all the subroutines is g ven in this dictionary Each entry is organized into subheadings that describe the subroutine and explain how to use it Other The storage in VASP is also compatible with the storage of general matrices in the IBM scientific subroutine package subheadings such as motivation and remarks are sometimes added to offer the user a better understanding of the theoretical background of the subroutine The dictionary proper lists only those routines that the user is expected to call directly Several additional subroutines used internally are also a part of VASP The user may however wish to call these routines himself since they are quite flexible These additional routines are described in appendix A and a complete listing of all programs is given in appendix B Several procedures have been written to facilitate the use of VASP on Ames time share system Their usage and listings are given in appendix C Table I lists all subroutines with their calling sequence and the TSS procedures for easy reference while table 2 lists the approximate storage used by each of the VASP subroutines when compiled on the NASA Ames 360 67 OS system Table 2 also lists the external references for each subroutine Some of the subroutines are almost direct copies from the Northrop FASP The detailed description of the theory
131. putations No flowchart is needed since there are no loops of any consequence Step 1 The matrix size is tested for reasonableness with an error exit if it is not KP 1 is set to special negative values to suppress reinversion by PSEUP and to change somewhat the matrix outputs This change is not discussed in PSEUP Step 2 Entry point PSEUP is used to diagonalize the input C holds a matrix T such that TAT B a matrix of 0 and 1 s If the rank only option is input the routine skips to step 4 Step 3 Subroutine INV puts an inverse of T into D The flag PIV is tested If zero INV failed the routine prints an error message INV uses matrix J Step 4 The matrix E which is matrix of 0 s and 1 s is scanned along its diagonal A matrix P of O s and 1 s is constructed such that I has all 1 s moved to extreme upper left corner A record of successive diagonal positions that are 0 is kept As each 1 is found in the diagonal in position IL the record is checked to see if there is an earlier 0 or 1 that needs to have a 1 permuted into its place j by permutation p If so a 1 is put into position j k of P Premultiplication by P will move position k k to j k Postmultiplication by pt will move jk to position j j Position is also marked as a hole that could be filled by a 1 lower on the diagonal since it vacates its old position The record in the first column of J has the structure of a q
132. r several steps until P t 7 gt P t which is the desired solution Subroutine RICAT will also stop after a specified number of steps if P has not con verged to a solution Finally having P and K we can compute the transient response of the system with optimum feedback from any desired initial condition The differential equation becomes x Fx GKx F GK x F x and the solution is x t ef ty The time history of the control is u t Kx t An alternate solution used in this example is to first calculate the transition matrix 2 where is the time increment at which the solution is desired then compute 2 x t x 0 7 xg The listing of a main program to solve this problem is given in figure 3 a the data for a particular case 15 given in figure 3 b and the first part of the results is given in figure 3 c In this problem I so the special case of AUG is used Asa result H is not used in AUG and need not have been used as an input 4 CV 0001 DIMENSION NF 2 NG 2 NFSTAR 2 NC 2 NCK 2 2 2 NRI 2 1 NZ 2 NPHI 2 NPO 2 NPT 2 MX0 2 SHXT 2 NAL 2 NA2 2 NH 2 uma co gt 353 LABUZ eset 0002 DOUBLE PRECISION 3 3 1 1 1 1 1 3 3 6 3 1 3 3 1 1 3 1 3 7 6 6 PHT 6 6 PO 3 3 PT 3 3 X0 3 XT 3 XX 3 ne ee 2 DELICLULEA4HUS 3 A143 3 A2032 3 4 DUMMYC72
133. rences to the same matrix for example CALL READ 1 A NA A NA A NA A NA A NA Output Arguments Matrices The first I of the matrices A B C D E Dimension arrays The first I of the arrays NA NB NC ND NE 2 RDTITL DESCRIPTION This subroutine reads a single card in hollerith format and loads it into the array TITLES It then calls LNCNT 100 This latter program in turn skips the printer to the next page prints the hollerith information in the array TITLES and positions the output to print next on line 3 USAGE CALL RDTITL It has no arguments 3 PRNT DESCRIPTION This subroutine prints a single matrix with or without a title line and either on the same page or a new page The matrix is printed using format 1P7D16 7 for the first line and 3x 1P7D16 7 for all subsequent lines The indents continuation lines for easier reading REMARKS The standard format is stored in arrays for the first line and FMT2 for subsequent lines both of which are stored in labeled COMMON as follows COMMON FORM NEPR FMT1 6 FMT2 6 where NEPR is the number of columns of data to be printed 7 in the standard case The standard format is loaded into COMMON in the BLKDATA program If other formats are desired they can be obtained either by changing the BLKDATA program or by having the users main program change the contents of COMMON CAUTION In writing a data statement for the formats put and FMT2 in s
134. ribution Statement Matrix computation Optimal control Kalman filtering 19 Security Classif of this report 20 Security Classif of this page 21 No of Pages 22 Price Unclassified Unclassified 167 3 00 For sale by the National Technical Information Service Springfield Virginia 22151 Unclassified Unlimited TABLE OF CONTENTS Page SU 1 INTRODUCTION 255 Xov 1 FEATURES OF THE PROGRAM 3 Universal Features 3 System Dependent Features 5 THE VASP DICTIONARY 5 EXAMPLE USES OF VASP 28 Example 1 Transient Response 28 Example 2 Transient Response Using TRNSI 34 Example 3 An Optimum Control Problem 4 Example 4 Sampled Data Ricatti Solution 49 Example 5 Matrix Decomposition 2 53 Example 6 Use of the Pseudoinverse Routine 64 APPENDIX A DESCRIPTION OF INTERNAL SUBROUTINES 69 APPENDIX B LISTINGS OF ALL VASP SUBROUTINES 90 APPENDIX C USE VASP ON AMES TSS _ 151 USERS MANUAL FOR THE VARIABLE DIMENSION AUTOMATIC SYNTHESIS PROGRAM VASP John S White and Homer Q Lee Ames Research Center
135. s DEP Values DEPI DEP2 DEP3 DEPI Default If zero user gets 2 D 6 used instead This number is multiplied times the largest mag nitude on the diagonal of B atstart If any trial pivots are found ess than this they are avoided as zero DEP2 Default If zero user gets 1 DO used instead Needed only if iteration The routine computes two numbers p q which would be zero if the first two Moore axioms were satisfied This num ber is raised to the power of the number of pivots found as a factor to use to make the product with the sum of p and q larger Mak ing this product larger tends to make the routine reject the current pivot Values between 1 and 2 work for ordinary purposes Note PSEUDO uses default values of and DEP2 7 IL 192 P3 IP4 Output Arguments Matrix B Matrix C Matrix EE 72 This is for output only It holds the last pivot actually accepted This gives the user or calling routine an estimate of the size of pivots found in case effective rank is not that desired operat ing with given value of If iterating this may be the last pivot rejected Parameter array of integers IP1 IP2 IP3 1 4 If zero do not iterate with side calculations If 1 iterate Note Other values should not be used since DECOM employs peculiar values If zero do all calculations otherwise do rank only Note Setting this to zero for each call is very use
136. s D D1 D2 D3 and D4 The last four are used only in iteration Maximal Rank Case A symmetric matrix from is placed B either directly as in PSEUP or indirectly from matrix multiplication A copy of B is put in D unless the rank only 70 iteration is used ANDRA is called to diagonalize B and place the result in C If the result is accepted TTRM puts the generalized inverse of B into EE Then the inverse of A is put into B The transpose may have to be used to get an answer for Singular Case The matrix U of the method is computed from C and put into C D holds original B EE Ct x D B EEXC ANDRA is called to diagonalize B Answer goes to LL TTRM puts pseudoinverse of B into D B CXD EE p x C The pseudoinverse is now in EE where the maximal rank case puts it Routine now forms pseudoinverse of teration Before each call of ANDRA the current values of B and are stored in D1 and D2 respectively B and C are changed when a new pivot is used in ANDRA BDNRM computes a number to decide if the pivot is to be rejected EE D3 and D4 are used as working storage in BDNRM EE actually has a matrix put in it that would be zero if the Moore Penrose axiom were perfectly satisfied If the pivot is rejected the old values from D1 and D2 are put back into B and C The work of the singular case is done next if the call was not made from DECOM Rank Only If iteration is used a full complement of
137. s to fill a total of 7 matrices Lines 14 15 and 16 set up the appropriate constants for RICAT specifying a print every step line 14 the maximum number of steps to be taken by RICAT line 15 and that both P and K should be printed line 16 Lines 17 18 store the initial values of and Pg in the running matrices and lines 19 through 23 do the necessary computations to obtain P and K called CK in program Then F and the transition matrix 2 lines 28 through 32 are computed and printed The transition matrix is labeled on the output lines 29 through 31 Lines 33 through 39 page the output print a heading for the transient response and print the first point Lines 40 through 47 then increment the solution and the time and print x t and u t called XT and DELTC in the program Example 4 Sampled Data Ricatti Solution This example is provided to show the general use of the subroutine SAMPL The theory of the example is given in the ASP manual page 222 and very briefly in the dictionary description of SAMPL page 24 in this manual A listing of the main program is shown in figure 4 a The data deck is shown in figure 4 b and the output in figure 4 c The main program is reasonably self explanatory The statement NCONT 2 4 line 13 indicates that SAMPL is to compute P for four successive time intervals and then stop Both P and K line 14 are to be printed at every step line 2 As mentioned in the dictio
138. t R be the matrix to decompose Matrices G and H are desired such that R HG H is to have only r nonzero columns G is to have only r nonzero rows Small r is the rank of R Case 1 Matrix IL is symmetric non negative definite Input R as the square input A to DECOM Then H and G are produced afterward from the matrices in the call statement as follows Parameter B is a diagonal matrix with r 175 H and G are computed by H DXB G DxB A original Case 2 R is nonsymmetric possibly not even positive definite Form RR subcase a or else form RtR subcase b subcases are chosen to give the smaller dimensions If R is square use RR to agree with both PSEU and DECOM Let this symmetric result be the input A to DECOM as in case 1 Obtain D and B as before and save them In subcase a X E but in subcase b X E X Rt Then for subcase a take H DB G XDB Similarly in subcase b take H X DB G DB Note The HandG matrices produced have the same dimensions as the smaller dimension of R If the rank of R is not maximal there will be zero rows or columns G If the matrix is used instead of B in the above calculations the zero rows or columns will be at the right or bot tom and the dimensions may be easily reduced This latter is the procedure used in subroutine DECGEN 89 Computation In practice the subroutine is very short it calls on PSEUP and INV to do the com
139. ted USAGE CALL SAMPL PHLNPHLH NH Q NQ R NR P NP K NK NCONT DUM KDUM Input Arguments Matrices PHLH Q R P Dimension arrays NPHI NH NQ NR NP Control arrays NCONT NCONT 1 j number of steps per print NCONT 2 n maximum number of steps NCONT 3 print control 1 no print 2 print P only 3 print K only 4 print both P and Output Arguments Matrices Dimension arrays Dummy Arguments Matrix DUM Constant KDUM Note KDUM contains the size of the DUM matrix which must be at least 6 NPHI 1 NPHI 2 22 TRNSI This subroutine computes t x t eFtx t_ f eF7 dr Gu to where u t JR Kx t it 25 and u is held constant for any interval specified by itj t to G t1 t i 0 1 2 The system output y t is given by y t Hx t The state vector x and system outputs y are printed every t intervals Also t must be positive integral multiple of t The program terminates at t gt See ASP manual pages 120 121 for reference USAGE CALL TRNSI F NF G NG J NJ R NR K NK H NH X NX T DUMMY KDUM Input Arguments Matrices F G J R K H X T Dimension arrays NF NG NJ NR NK NH NX Note Dimension of T is 4 where t T2 t T3 te T4 to Dummy Argument Matrix DUMMY Constant KDUM Note KDUM contains the size of the dummy matrix which must be at least 4 NF 1 NF 2 23 PSEUDO DESCRIPTION This subroutine computes the Moore Penrose generalized
140. this data set it will generate all the procedures and place them in the user s USERLIB TSS ACCESS For access to the VASP program an Ames TSS user should issue the following statements SHARE VASP FSTJSW VASP which allows access to the VASP subroutines SHARE VASPPROC FSTJSW VASPPROC EXECUTE VASPPROC which first allows access to a data set containing the various procedures and then enters these procedures in the user s USERLIB Note that the EXECUTE command sets up a batch job and that the procedures will not be available until that batch job is completed and the user has issued either a LOGOFF or ABEND command After once issuing these commands the user need only call the procedure as discussed earlier Further for conversational use issue the command SHARE VASPI FSTJSW VASP which allows access to the proper version of subroutine INPUT 157 VASPS JBLB VASP LOAD BLKDTA Input Output DDEF Default Options CONVASP JBLB VASP1 DISPLAY Matrix Names Edit SOURCE VASPMN Compile VASPMN Load VASPMN Edit SOURCE MNPG EXCERPT beginning of Fortran Programs Display FORTRAN STATEMENTS CMPL Add end of Fortran Program Display MNPG Now Compiling Compile MNPG Call MNPG Display COMPUTING DONE RECMPT Call MINPG REWRT EDIT SOURCE MNPG EXCISE Program Display FORTRAN STATEMENTS CHNGIN Change Input DDEF CHNGOUT Change Output DDEF CLRVASP Erase all Programs
141. this program to help him track down errors Subroutine ASPERR calls in turn a system program which provides the actual walkback In Ames OS this system routine is called ERRTRA while in Ames TSS it is called TRACE The calling statement should be changed to match the user s operating system or else deleted altogether 27 BLKDATA DESCRIPTION This is an installation dependent subroutine It loads certain common areas used by VASP with appropriate constants as follows 1 COMMON FQRM NEPR FMT1 6 FMT2 6 These three variables control the printing procedure and are set to 7 1P7D16 7 and 3x 1P7D16 7 respectively They assume a line length of at least 115 characters 2 COMMON LINES NLP LIN TITLE 23 NLP controls the number of lines per page and is set at 45 to agree with the NASA Ames system It should be changed to match each installation LIN is a counter which keeps track of the number of lines printed on each page It is incremented and used only in LNCNT TITLE contains 72 blank characters which can be loaded as desired by use of RDTITL plus 20 more characters containing VASP PROGRAM Subroutines LNCNT prints TITLE at the head of every page 3 COMMON MAX MAXRC MAXRC is used by most subroutines to check the reasonableness of the matrix dimensions The user should set MAXRC to match the storage available for each matrix It is preset to 6400 28 PSEU SUMMARY PSEU is a FORTRAN routine to find the M
142. ting the test constant MAXRC to 25 The error test in the READ subroutine will note that the product of NC 1 and NC 2 is greater than MAXRC and will return an error message This selection of MAXRC will limit all other VASP arrays to 25 so it is frequently desirable to dimension all arrays the same Occasionally the user may wish to refer to a single element of a matrix Since FORTRAN statements use the FORTRAN dimension statement a reference to B 4 4 in the previous example will select the 19 element in the B storage However VASP using the VASP dimension has stored B 4 4 in the 16 element of B which is not the same To actually select a specific element say 80 1 one must refer to B j 1 NB 2 1 1 In the above example references to A i j will be correct since the FORTRAN and VASP dimensions are the same System Dependent Features Two subroutines in the VASP package are system dependent The firstis BLKDTA Data statements in this subroutine control the printing They require a printer with at least 115 charac ters per line and place 45 lines on each page These requirements may be changed as needed The second is ASPERR which calls a system subroutine for error tracing The description of ASPERR indicates any necessary changes to match the system The VASP programs frequently generate very small numbers The computer operating system may detect these small numbers as underflows and print error messages If so the user sh
143. tions include Tech Briefs technical information generated under a NASA Technology Utilization Reports and contract or grant and considered an important ee ee Technology Surveys contribution to existing knowledge zd d _ Details the availability of th se publications may be obtained from SCIENTIFIC AND TECHNICAL INFORMATION OFFICE NATIONAL AERONAUTICS AND SPACE ADMINISTRATION Washinston D C 90546
144. try assumes this Matrices formed by computer arithmetic will not always be symmetric Hence the routine always forces the symmetric matrix B to actually be symmetric by taking the average of the element and its transpose The nonsymmetric entry unfortunately approximately squares the ratio between largest and smallest eigenvalue There is a nonsymmetric feature The routine choses AAT instead of the other way around if A is a square matrix This arbitrary choice agrees with the DECOM routine and the ASP routines As a result in the singular case multiplying A by its pseudo inverse from the left is more likely to give a diagonal matrix of 175 0 than multiplying from the right side of A Pivot Size is used to compute a smallest allowable pivot In no case is it reasonable or desirable to worry about exact equality in the use of such tolerances Fortunately work with ill conditioned systems shows a series of pivots that decrease steadily in magnitude Furthermore the first bad erroneous pivot is at most 10 to 1000 times smaller than its predecessors Since ANDRA is choos ing largest pivots first the user has considerable latitude in actual choice All positive elements can be accepted if the matrix is known to be nonsingular by choosing DEP very small By choosing DEP very large say nearly 1 0 the routine can be forced to reject pivots after the first At present there is no way of making it start it
145. ueue Matrix P isin DI the second matrix of dummy array D Step 5 Return NOMENCLATURE The nomenclature used in DECOM is compatible with that used in PSEU but differs from that used in the ASP manual description of the decomposition routine p 154 Also since DECGEN requires dummy storage the nomenclature there is different again The following table lists the correlations 90 D1 DECGEN DUMI DUM N7 DUM N4 DUM N5 DUM N2 DUM N3 DUM N6 ASP AAT 91 92 APPENDIX B LISTINGS OF ALL VASP SUBROUTINES q C6 SUBROUTINE READ I By DIMENS ION A 1 Bll Cll D 1 DOUBLE PRECISION READ 5 1100 2 1 CALL READ1 A NA NZ LAB lt E0 1 9959 5 100 LAB NZ 1 CALL READ1 B NB NZ LAB EO 2 999 READ 5 1 NZ 1 CALL RFAD1 C NCX NZ IF 0 3 GO TO 999 READ 5 100 LAR NZ 1 CALL READ1 D ND NZ LAB IC I 4 GO TO 999 100 LAR 1 RFAD 5 CALL RFAD1 F LAP 00 4 4 214 999 RFTURN FND SUBROUTINE RDTITL COMM OM REAN 100 FORMAT LINES NLP LIN TITLE 23 5 100 1 1 1 18 18 4 CALL LNCNT 100 RETURN END 1 DIMENSION NA 2 NB 2 NCX 2 ND 2 gt 2 2 NZ 2 M7 2 M7 2 M7 NCX

NASA Technical Memorandum X-2417

Contents

Download Pdf Manuals

Related Search

Related Contents