%!PS %%Version: 3.3.1 %%DocumentFonts: (atend) %%Pages: (atend) %%EndComments % % Version 3.3.1 prologue for troff files. % /#copies 1 store /aspectratio 1 def /formsperpage 1 def /landscape false def /linewidth .3 def /magnification 1 def /margin 0 def /orientation 0 def /resolution 720 def /rotation 1 def /xoffset 0 def /yoffset 0 def /roundpage true def /useclippath true def /pagebbox [0 0 612 792] def /R /Times-Roman def /I /Times-Italic def /B /Times-Bold def /BI /Times-BoldItalic def /H /Helvetica def /HI /Helvetica-Oblique def /HB /Helvetica-Bold def /HX /Helvetica-BoldOblique def /CW /Courier def /CO /Courier def /CI /Courier-Oblique def /CB /Courier-Bold def /CX /Courier-BoldOblique def /PA /Palatino-Roman def /PI /Palatino-Italic def /PB /Palatino-Bold def /PX /Palatino-BoldItalic def /Hr /Helvetica-Narrow def /Hi /Helvetica-Narrow-Oblique def /Hb /Helvetica-Narrow-Bold def /Hx /Helvetica-Narrow-BoldOblique def /KR /Bookman-Light def /KI /Bookman-LightItalic def /KB /Bookman-Demi def /KX /Bookman-DemiItalic def /AR /AvantGarde-Book def /AI /AvantGarde-BookOblique def /AB /AvantGarde-Demi def /AX /AvantGarde-DemiOblique def /NR /NewCenturySchlbk-Roman def /NI /NewCenturySchlbk-Italic def /NB /NewCenturySchlbk-Bold def /NX /NewCenturySchlbk-BoldItalic def /ZD /ZapfDingbats def /ZI /ZapfChancery-MediumItalic def /S /S def /S1 /S1 def /GR /Symbol def /inch {72 mul} bind def /min {2 copy gt {exch} if pop} bind def /setup { counttomark 2 idiv {def} repeat pop landscape {/orientation 90 orientation add def} if /scaling 72 resolution div def linewidth setlinewidth 1 setlinecap pagedimensions xcenter ycenter translate orientation rotation mul rotate width 2 div neg height 2 div translate xoffset inch yoffset inch neg translate margin 2 div dup neg translate magnification dup aspectratio mul scale scaling scaling scale addmetrics 0 0 moveto } def /pagedimensions { useclippath userdict /gotpagebbox known not and { /pagebbox [clippath pathbbox newpath] def roundpage currentdict /roundpagebbox known and {roundpagebbox} if } if pagebbox aload pop 4 -1 roll exch 4 1 roll 4 copy landscape {4 2 roll} if sub /width exch def sub /height exch def add 2 div /xcenter exch def add 2 div /ycenter exch def userdict /gotpagebbox true put } def /addmetrics { /Symbol /S null Sdefs cf /Times-Roman /S1 StandardEncoding dup length array copy S1defs cf } def /pagesetup { /page exch def currentdict /pagedict known currentdict page known and { page load pagedict exch get cvx exec } if } def /decodingdefs [ {counttomark 2 idiv {y moveto show} repeat} {neg /y exch def counttomark 2 idiv {y moveto show} repeat} {neg moveto {2 index stringwidth pop sub exch div 0 32 4 -1 roll widthshow} repeat} {neg moveto {spacewidth sub 0.0 32 4 -1 roll widthshow} repeat} {counttomark 2 idiv {y moveto show} repeat} {neg setfunnytext} ] def /setdecoding {/t decodingdefs 3 -1 roll get bind def} bind def /w {neg moveto show} bind def /m {neg dup /y exch def moveto} bind def /done {/lastpage where {pop lastpage} if} def /f { dup /font exch def findfont exch dup /ptsize exch def scaling div dup /size exch def scalefont setfont linewidth ptsize mul scaling 10 mul div setlinewidth /spacewidth ( ) stringwidth pop def } bind def /changefont { /fontheight exch def /fontslant exch def currentfont [ 1 0 fontheight ptsize div fontslant sin mul fontslant cos div fontheight ptsize div 0 0 ] makefont setfont } bind def /sf {f} bind def /cf { dup length 2 idiv /entries exch def /chtab exch def /newencoding exch def /newfont exch def findfont dup length 1 add dict /newdict exch def {1 index /FID ne {newdict 3 1 roll put}{pop pop} ifelse} forall newencoding type /arraytype eq {newdict /Encoding newencoding put} if newdict /Metrics entries dict put newdict /Metrics get begin chtab aload pop 1 1 entries {pop def} for newfont newdict definefont pop end } bind def % % A few arrays used to adjust reference points and character widths in some % of the printer resident fonts. If square roots are too high try changing % the lines describing /radical and /radicalex to, % % /radical [0 -75 550 0] % /radicalex [-50 -75 500 0] % % Move braceleftbt a bit - default PostScript character is off a bit. % /Sdefs [ /bracketlefttp [201 500] /bracketleftbt [201 500] /bracketrighttp [-81 380] /bracketrightbt [-83 380] /braceleftbt [203 490] /bracketrightex [220 -125 500 0] /radical [0 0 550 0] /radicalex [-50 0 500 0] /parenleftex [-20 -170 0 0] /integral [100 -50 500 0] /infinity [10 -75 730 0] ] def /S1defs [ /underscore [0 80 500 0] /endash [7 90 650 0] ] def % % Tries to round clipping path dimensions, as stored in array pagebbox, so they % match one of the known sizes in the papersizes array. Lower left coordinates % are always set to 0. % /roundpagebbox { 7 dict begin /papersizes [8.5 inch 11 inch 14 inch 17 inch] def /mappapersize { /val exch def /slop .5 inch def /diff slop def /j 0 def 0 1 papersizes length 1 sub { /i exch def papersizes i get val sub abs dup diff le {/diff exch def /j i def} {pop} ifelse } for diff slop lt {papersizes j get} {val} ifelse } def pagebbox 0 0 put pagebbox 1 0 put pagebbox dup 2 get mappapersize 2 exch put pagebbox dup 3 get mappapersize 3 exch put end } bind def %%EndProlog %%BeginSetup mark /resolution 720 def setup 2 setdecoding %%EndSetup %%Page: 1 1 /saveobj save def mark 1 pagesetup 10 R f (Appendix 2)1 469 1 2825 120 t (SYMMETRIC MATRICES)1 1126 1 2497 360 t ( Estimation)1 459( Condition)1 576(SYCE -)1 389 3 720 720 t ( DeComposition)1 809(SYDC -)1 400 2 720 840 t ( and Back Solve)3 652( Forward)1 513(SYFBS -)1 440 3 720 960 t ( Equation solution)2 734( Linear)1 435(SYLE -)1 383 3 720 1080 t (SYMD -)1 422 1 720 1200 t 10 I f (MDM)1317 1200 w 7 I f (T)1566 1160 w 10 R f (decomposition)1638 1200 w ( MuLtiplication)1 781(SYML -)1 411 2 720 1320 t ( NorM)1 419(SYNM -)1 422 2 720 1440 t ( Solution)1 365( System)1 470(SYSS -)1 373 3 720 1560 t cleartomark showpage saveobj restore %%EndPage: 1 1 %%Page: 2 2 /saveobj save def mark 2 pagesetup 10 R f (SYCE \320 decomposition of a symmetric matrix with condition estimation)9 2953 1 1403 120 t 9 B f (Purpose:)720 480 w 10 R f ( lower bound for the condition num-)6 1460(SYCE \(SYmmetric matrix Condition Estimation\) gives a)6 2284 2 1296 480 t ( also supplies the MDM)4 961( It)1 111(ber of a symmetric matrix A, which need not be positive de\256nite.)11 2616 3 1296 600 t 7 R f (T)4989 560 w 10 R f (decomposition of A and may be used in a linear equation package.)11 2656 1 1296 720 t 9 B f (Usage:)720 1080 w 10 R f (CALL SYCE \(N, C, INTER, COND\))5 1499 1 1296 1080 t (N)1440 1320 w 14 S f (\256)1980 1340 w 10 R f (the number of rows in A)5 979 1 2196 1320 t (C)1440 1560 w 14 S f (\256)1980 1580 w 10 R f ( lower tri-)2 417(a one-dimensional array of length N\(N+1\)/2 into which the)8 2427 2 2196 1560 t ( A is packed by columns as illustrated in the)9 1804(angular part of the matrix)4 1040 2 2196 1680 t (following 4)1 464 1 2196 1800 t 10 S f (\264)2668 1800 w 10 R f (4 example:)1 441 1 2731 1800 t 10 I f (a)2232 2040 w 7 R f (11)2293 2060 w 10 I f (c)3312 2040 w 7 R f (1)3367 2060 w 10 I f (a)2232 2160 w 7 R f (21)2293 2180 w 10 I f (a)2412 2160 w 7 R f (22)2473 2180 w 14 S f (\256)2952 2160 w 10 I f (c)3312 2160 w 7 R f (2)3367 2180 w 10 I f (c)3451 2160 w 7 R f (5)3506 2180 w 10 I f (a)2232 2280 w 7 R f (31)2293 2300 w 10 I f (a)2412 2280 w 7 R f (32)2473 2300 w 10 I f (a)2592 2280 w 7 R f (33)2653 2300 w 10 I f (c)3312 2280 w 7 R f (3)3367 2300 w 10 I f (c)3451 2280 w 7 R f (6)3506 2300 w 10 I f (c)3590 2280 w 7 R f (8)3645 2300 w 10 I f (a)2232 2400 w 7 R f (41)2293 2420 w 10 I f (a)2412 2400 w 7 R f (42)2473 2420 w 10 I f (a)2592 2400 w 7 R f (43)2653 2420 w 10 I f (a)2772 2400 w 7 R f (44)2833 2420 w 10 I f (c)3312 2400 w 7 R f (4)3367 2420 w 10 I f (c)3451 2400 w 7 R f (7)3506 2420 w 10 I f (c)3590 2400 w 7 R f (9)3645 2420 w 10 I f (c)3729 2400 w 7 R f (10)3784 2420 w 14 S f (\254)1980 2660 w 10 R f (the D and M matrices of the decomposition for SYFBS)9 2207 1 2196 2640 t (\(see)2196 2760 w 8 B f (Note 2)1 219 1 2381 2760 t 10 R f (\))2600 2760 w (INTER)1440 3000 w 14 S f (\254)1980 3020 w 10 R f ( the interchanges,)2 719(an integer vector of length N containing a record of)9 2125 2 2196 3000 t (i.e. the matrix P, described in)5 1171 1 2196 3120 t 8 B f (Note 2)1 219 1 3392 3120 t 10 R f (below.)3636 3120 w (COND)1440 3360 w 14 S f (\254)1980 3380 w 10 R f (an estimate of the condition number of A \(see)8 1830 1 2196 3360 t 8 B f (Note 1)1 219 1 4051 3360 t 10 R f (\))4270 3360 w 9 B f (Note 1:)1 278 1 720 3720 t 10 R f ( to errors in)3 481(The condition number measures the sensitivity of the solution of a linear system)12 3263 2 1296 3720 t ( the elements of the matrix and the right-hand side\(s\))9 2134( If)1 118( and in the right-hand side.)5 1081(the matrix)1 411 4 1296 3840 t (of your linear system have)4 1091 1 1296 3960 t 10 B f (d)2420 3960 w 10 R f ( solution might have as few as)6 1264(decimal digits of precision, the)4 1267 2 2509 3960 t 10 B f (d)1296 4080 w 10 S f (-)1377 4080 w 10 R f (log)1457 4080 w 7 R f (10)1596 4100 w 10 R f ( if COND is greater than 10)6 1121( Thus)1 252(\(COND\) correct decimal digits.)3 1270 3 1674 4080 t 7 R f (Bd)4327 4040 w 7 I f (P)4419 4040 w 10 R f ( be)1 120(, there may)2 450 2 4470 4080 t (no correct digits.)2 674 1 1296 4200 t 9 B f (Note 2:)1 278 1 720 4560 t 10 R f (The MDM)1 435 1 1296 4560 t 7 R f (T)1736 4520 w 10 R f ( satis\256es P)2 435(decomposition of a symmetric matrix A)5 1621 2 1817 4560 t 7 R f (T)3878 4520 w 10 R f (AP = MDM)2 496 1 3929 4560 t 7 R f (T)4430 4520 w 10 R f (, where P is a)4 559 1 4481 4560 t ( D is a block diagonal matrix)6 1232(permutation matrix, M is a unit lower triangular matrix, and)9 2512 2 1296 4680 t (with blocks of order 1)4 922 1 1296 4800 t 10 S f (\264)2226 4800 w 10 R f (1 and blocks of order 2)5 973 1 2289 4800 t 10 S f (\264)3270 4800 w 10 R f (2. Whenever)1 544 1 3333 4800 t 10 I f (d)3912 4800 w 7 I f (i)3973 4820 w 7 S f (+)4009 4820 w 7 R f (1 ,)1 58 1 4059 4820 t 7 I f (i)4122 4820 w 10 R f ( 2)1 86(is nonzero \(in a)3 648 2 4185 4800 t 10 S f (\264)4927 4800 w 10 R f (2)4990 4800 w (block of D\),)2 491 1 1296 4920 t 10 I f (m)1815 4920 w 7 I f (i)1898 4940 w 7 S f (+)1934 4940 w 7 R f (1 ,)1 58 1 1984 4940 t 7 I f (i)2047 4940 w 10 R f ( return from SYCE,)3 797( On)1 175(is zero.)1 291 3 2103 4920 t 10 I f (d)3394 4920 w 7 I f (ii)3455 4940 w 10 R f ( posi-)1 227(, the diagonal of D, occupies the)6 1310 2 3503 4920 t (tion of C which contained)4 1062 1 1296 5040 t 10 I f (a)2389 5040 w 7 I f (ii)2450 5060 w 10 R f ( the elements of the strictly lower portions of M)9 1975(on entry, and)2 536 2 2529 5040 t ( the diagonal elements of M)5 1158( Since)1 280(and D appear permuted in the remaining positions of C.)9 2306 3 1296 5160 t ( INTER contain information for con-)5 1513( positive elements of)3 848( The)1 211(are all 1, they are not stored.)6 1172 4 1296 5280 t ( if any, in-)3 419(structing P \(see the introduction to this chapter\). The negative elements of INTER,)12 3325 2 1296 5400 t (dicate the presence of 2)4 989 1 1296 5520 t 10 S f (\264)2326 5520 w 10 R f ( is negative, then D contains a 2)7 1367(2 blocks in D. If INTER\(I\))5 1130 2 2422 5520 t 10 S f (\264)4927 5520 w 10 R f (2)4990 5520 w (block beginning at row I)4 982 1 1296 5640 t 10 S f (-)2278 5640 w 10 R f ( this case,)2 391(1. In)1 208 2 2333 5640 t 10 I f (d)2957 5640 w 7 I f (i)3018 5660 w 7 R f (,)3043 5660 w 7 I f (i)3066 5660 w 7 S f (-)3102 5660 w 7 R f (1)3152 5660 w 10 R f (directly follows)1 630 1 3220 5640 t 10 I f (d)3875 5640 w 7 I f (i)3936 5660 w 7 S f (-)3972 5660 w 7 R f (1 ,)1 58 1 4022 5660 t 7 I f (i)4085 5660 w 7 S f (-)4121 5660 w 7 R f (1)4171 5660 w 10 R f (in C.)1 195 1 4239 5640 t 9 B f (Note 3:)1 278 1 720 6000 t 10 R f ( the)1 160(For complex Hermitian matrices \(A = A*, where A* is the conjugate transpose of A\),)14 3584 2 1296 6000 t ( MDM* decomposition and re-)4 1287(complex Hermitian version of this subroutine computes the)7 2457 2 1296 6120 t (turns the conjugate of M rather than M in C.)9 1770 1 1296 6240 t 9 B f (Error situations:)1 653 1 720 6600 t 10 R f ( see)1 183( \320)1 125( those errors marked with an asterisk)6 1505(*\(The user can elect to `recover' from)6 1546 4 1517 6600 t 10 I f (Er-)4907 6600 w (ror Handling)1 531 1 1512 6720 t 10 R f (, Framework Chapter\))2 884 1 2043 6720 t (Linear Algebra)1 606 1 360 7500 t (\320 2 \320)2 300 1 2730 7620 t (SYCE)360 7740 w cleartomark showpage saveobj restore %%EndPage: 2 2 %%Page: 3 3 /saveobj save def mark 3 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYCE)1 4305(February 11, 1993)2 735 2 360 840 t (Number Error)1 1506 1 1800 1320 t (1 N)1 1008 1 1944 1560 t 10 S f (<)2977 1560 w 10 R f (1)3057 1560 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 1800 t 9 B f (Double-precision version:)1 988 1 720 2160 t 10 R f (DSYCE with C declared double precision.)5 1698 1 1800 2160 t 9 B f (Complex symmetric version:)2 1106 1 720 2520 t 10 R f (CSYCE with C declared complex)4 1349 1 1901 2520 t 9 B f (Complex Hermitian version:)2 1101 1 720 2880 t 10 R f (CHECE with C declared complex \(see)5 1539 1 1896 2880 t 8 B f (Note 3)1 219 1 3460 2880 t 10 R f (\).)3679 2880 w 9 B f (Storage:)720 3240 w 10 R f ( storage in the)3 566(N real \(double precision for DSYCE, complex for CSYCE\) locations of scratch)11 3178 2 1296 3240 t (dynamic storage stack)2 887 1 1296 3360 t 9 B f (Time:)720 3770 w 10 R f (6)1356 3840 w 10 I f (N)1321 3710 w 7 R f (3)1399 3670 w 10 S1 f (_ ___)1 151 1 1306 3740 t 10 S f (+)1507 3770 w 10 R f (4)1652 3840 w (19)1627 3710 w 10 S1 f (_ __)1 130 1 1612 3740 t 10 I f (N)1760 3770 w 7 R f (2)1838 3730 w 10 S f (+)1921 3770 w 10 R f (6)2066 3840 w (25)2041 3710 w 10 S1 f (_ __)1 130 1 2026 3740 t 10 I f (N)2174 3770 w 10 R f (additions)2266 3770 w (6)1356 4180 w 10 I f (N)1321 4050 w 7 R f (3)1399 4010 w 10 S1 f (_ ___)1 151 1 1306 4080 t 10 S f (+)1507 4110 w 10 R f (4)1652 4180 w (13)1627 4050 w 10 S1 f (_ __)1 130 1 1612 4080 t 10 I f (N)1760 4110 w 7 R f (2)1838 4070 w 10 S f (+)1921 4110 w 10 R f (6)2041 4180 w (5)2041 4050 w 10 S1 f (_ _)1 80 1 2026 4080 t 10 I f (N)2124 4110 w 10 R f (multiplications)2216 4110 w (2)1356 4520 w 10 I f (N)1321 4390 w 7 R f (2)1399 4350 w 10 S1 f (_ ___)1 151 1 1306 4420 t 10 S f (+)1507 4450 w 10 R f (2)1627 4520 w (5)1627 4390 w 10 S1 f (_ _)1 80 1 1612 4420 t 10 I f (N)1710 4450 w 10 R f (divisions)1802 4450 w (at most)1 292 1 1296 4680 t 10 I f (N)1613 4680 w 7 R f (2)1691 4640 w 10 S f (+)1750 4680 w 10 I f (N)1821 4680 w 10 R f (comparisons)1913 4680 w 9 B f (Method:)720 5040 w 10 R f ( - Kaufman algorithm, described in reference [1] below, is used to determine the)13 3297(The Bunch)1 447 2 1296 5040 t ( algorithm in [2] is used to get the condition estimate.)10 2141(decomposition. The)1 819 2 1296 5160 t (SYCE calls SYMD after setting EPS to 0.0)7 1728 1 1296 5400 t 9 B f (See also:)1 333 1 720 5760 t 10 R f (SYLE, SYFBS, SYMD, SYDC, SYSS)4 1553 1 1296 5760 t 9 B f (Author:)720 6120 w 10 R f (Linda Kaufman)1 629 1 1296 6120 t (Linear Algebra)1 606 1 4794 7500 t (\320 3 \320)2 300 1 2730 7620 t (SYCE)5144 7740 w cleartomark showpage saveobj restore %%EndPage: 3 3 %%Page: 4 4 /saveobj save def mark 4 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYCE February)1 4665 2 360 840 t 9 B f (References:)720 1320 w 10 R f ( Kaufman, L., and Parlett, B., Decomposition of a symmetric matrix,)10 2790( J. R.,)2 237([1] Bunch,)1 538 3 1296 1320 t 10 I f (Nu-)4890 1320 w (mer. Math 27)2 541 1 1548 1440 t 10 R f (\(1976\), 95-109.)1 624 1 2114 1440 t ( H., An estimate for the)5 980( A. K., Moler, C. B., Stewart, G. W., and Wilkinson, J.)11 2270([2] Cline,)1 494 3 1296 1680 t (condition number,)1 733 1 1548 1800 t 10 I f (SIAM J. Numer. Anal. 16)4 1007 1 2331 1800 t 10 R f (\(1979\), 368-375.)1 674 1 3363 1800 t 9 B f (Example:)720 2160 w 10 R f ( is an encoding of the iterative re\256nement algorithm which may be used to ob-)14 3199(This example)1 545 2 1296 2160 t ( with an ill-conditioned)3 1017(tain a highly accurate solution to a system of linear equations)10 2727 2 1296 2280 t ( is not excessively high, the program usually re-)8 1981( the condition number)3 901( If)1 123(coef\256cient matrix.)1 739 4 1296 2400 t (turns a solution that is accurate to the working precision of the machine.)12 2882 1 1296 2520 t (The iterative re\256nement algorithm is essentially:)5 1940 1 1296 2760 t ( Ax = b)3 303(\(1\) Solve)1 408 2 1980 3000 t ( tol =)2 212(\(2\) Set)1 308 2 1980 3120 t 10 S f (e)2525 3120 w 15 S f (S)2626 3150 w 10 S f (\357)2764 3137 w 10 R f (x)2812 3120 w 7 I f (i)2873 3140 w 10 S f (\357)2901 3137 w 10 R f (where)2160 3250 w 10 S f (e)2428 3250 w 10 R f (is the precision of the machine)5 1223 1 2497 3250 t ( in double precision the residual)5 1279(\(3\) Compute)1 547 2 1980 3370 t (r = Ax)2 261 1 2310 3490 t 10 S f (-)2596 3490 w 10 R f (b)2676 3490 w ( A)1 97(\(4\) Solve)1 408 2 1980 3610 t 10 S f (d)2510 3610 w 10 R f (x = r)2 189 1 2559 3610 t ( norm =)2 317(\(5\) Compute)1 547 2 1980 3730 t 15 S f (S)2885 3760 w 10 S f (\357)3023 3747 w (d)3071 3730 w 10 R f (x)3128 3730 w 7 I f (i)3189 3750 w 10 S f (\357)3217 3747 w 10 R f ( x to x +)4 334(\(6\) Set)1 308 2 1980 3860 t 10 S f (d)2647 3860 w 10 R f (x)2721 3860 w ( norm)1 236(\(7\) If)1 246 2 1980 3980 t 10 S f (\243)2487 3980 w 10 R f (tol stop, else return to step 3)6 1130 1 2567 3980 t ( subroutines SYCE)2 777(In the program below, step \(1\) is accomplished using the two lower-level)11 2967 2 1296 4220 t ( the system using)3 707( decomposes A into several factors and SYFBS solves)8 2218( SYCE)1 312(and SYFBS.)1 507 4 1296 4340 t ( A is destroyed by SYCE and needed in step \(3\) of the algorithm, a copy)15 2941( Since)1 275(these factors.)1 528 3 1296 4460 t ( step \(4\) the decomposition created earlier in SYCE is reused and)11 2620( In)1 134( is saved.)2 371(of the A matrix)3 619 4 1296 4580 t ( it is possible that the matrix is so)8 1348( Since)1 273( required.)1 383(only the forward and back solver SYFBS is)7 1740 4 1296 4700 t ( steps \(3\) through \(7\) in)5 1001(ill-conditioned that the iterative re\256nement algorithm will diverge,)7 2743 2 1296 4820 t ( number is chosen to be an up-)7 1230( This)1 229( only a \256nite number of times.)6 1214(the program are performed)3 1071 4 1296 4940 t ( number of bits in the mantissa of the \257oating-point number supported by)12 3040(per bound on the)3 704 2 1296 5060 t (the machine.)1 510 1 1296 5180 t (This algorithm is not included in)5 1394 1 1296 5420 t 8 R f (PORT)2733 5420 w 10 R f (because for double-precision matrices part of the)6 2058 1 2982 5420 t (computation would have to be done in extended precision.)8 2333 1 1296 5540 t (Linear Algebra)1 606 1 360 7500 t (\320 4 \320)2 300 1 2730 7620 t (SYCE)360 7740 w cleartomark showpage saveobj restore %%EndPage: 4 4 %%Page: 5 5 /saveobj save def mark 5 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYCE)1 4305(February 11, 1993)2 735 2 360 840 t 8 CW f (INTEGER N, JEND, IREAD, I1MACH, I, JBEGIN, J, IWRITE)8 2496 1 1992 1300 t (INTEGER INTER\(6\), IEND, ITER, L, IFIX)5 1776 1 1992 1400 t (REAL C\(20\), SAVEC\(36\), B\(6\), SAVEB\(6\), R\(6\))5 2064 1 1992 1500 t (REAL COND, R1MACH, BNORM, RNORM, ABS, ALOG10)6 2112 1 1992 1600 t (DOUBLE PRECISION D\(6\))2 1008 1 1992 1700 t (N=5)1992 1800 w (C)1656 1900 w (C READ IN A SYMMETRIC MATRIX WHOSE UPPER TRIANGULAR)8 2448 1 1656 2000 t (C PORTION IS STORED ONE ROW PER CARD. MAKE A)9 2112 1 1656 2100 t (C COPY OF THE MATRIX SO THAT IT CAN BE USED LATER)11 2352 1 1656 2200 t (C)1656 2300 w (JEND=0)2040 2400 w (IREAD=I1MACH\(1\))2040 2500 w (DO 20 I=1,N)2 528 1 2040 2600 t (JBEGIN=JEND+1)2184 2700 w (JEND=JBEGIN+N - I)2 816 1 2184 2800 t (READ\(IREAD,1\)\(C\(J\),J=JBEGIN,JEND\))2184 2900 w (1 FORMAT\(5F8.0\))1 1056 1 1752 3000 t (DO 10 J=JBEGIN,JEND)2 912 1 2184 3100 t (SAVEC\(J\)=C\(J\))2424 3200 w (10 CONTINUE)1 816 1 1752 3300 t (20 CONTINUE)1 672 1 1752 3400 t (C READ IN RIGHT HAND SIDE AND SAVE IT)8 1776 1 1656 3500 t (DO 30 I=1,N)2 528 1 2040 3600 t (READ\(IREAD,1\)B\(I\))2184 3700 w (SAVEB\(I\)=B\(I\))2184 3800 w (30 CONTINUE)1 672 1 1752 3900 t (C)1656 4000 w ( AX = B USING SEPARATE CALLS TO SYCE AND SYFBS)10 2208(C SOLVE)1 384 2 1656 4100 t (C)1656 4200 w (CALL SYCE\(N,C,INTER,COND\))1 1200 1 1992 4300 t (CALL SYFBS\(N,C,B,6,1,INTER\))1 1296 1 1992 4400 t (IWRITE=I1MACH\(2\))1992 4500 w (IF\(COND.GE.1.0/R1MACH\(4\)\)WRITE\(IWRITE,31\))1992 4600 w ( CONDITION NUMBER HIGH,ACCURATE SOLUTION UNLIKELY\))5 2400(31 FORMAT\(49H)1 720 2 1752 4700 t (WRITE\(IWRITE,32\) COND)1 1008 1 1992 4800 t ( CONDITION NUMBER IS ,1PE16.8\))4 1440(32 FORMAT\(21H)1 720 2 1752 4900 t ( NORM OF SOLUTION)3 816(C COMPUTE)1 672 2 1656 5000 t (BNORM=0.0)1992 5100 w (WRITE\(IWRITE,33\))1992 5200 w ( THE FIRST SOLUTION X, FROM SYCE AND SYFBS=\))8 2112(33 FORMAT\(43H)1 720 2 1752 5300 t (DO 40 I=1,N)2 528 1 1992 5400 t (BNORM=BNORM+ABS\(B\(I\)\))2136 5500 w (40 WRITE\(IWRITE,41\)B\(I\))1 1344 1 1752 5600 t ( ,F20.7\))1 384(41 FORMAT\(1H)1 672 2 1752 5700 t (C)1656 5800 w (C IEND IS THE UPPER BOUND ON THE NUMBER OF BITS PER WORD)12 2688 1 1656 5900 t (C)1656 6000 w (IEND=I1MACH\(11\)*IFIX\(R1MACH\(5\)/ALOG10\(2.0\)+1.0\))1992 6100 w (C)1656 6200 w (C REFINE SOLUTION)2 816 1 1656 6300 t (C)1656 6400 w (DO 90 ITER=1,IEND)2 816 1 1992 6500 t (C)1656 6600 w (C COMPUTE RESIDUAL R = B - AX, IN DOUBLE PRECISION)10 2400 1 1656 6700 t (C)1656 6800 w (DO 50 I=1,N)2 528 1 2136 6900 t 10 R f (Linear Algebra)1 606 1 4794 7560 t (\320 5 \320)2 300 1 2730 7680 t (SYCE)5144 7800 w cleartomark showpage saveobj restore %%EndPage: 5 5 %%Page: 6 6 /saveobj save def mark 6 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYCE February)1 4665 2 360 840 t 8 CW f (50 D\(I\)=DBLE\(SAVEB\(I\)\))1 1440 1 1752 1300 t (L=1)2136 1400 w (DO 70 I=1,N)2 528 1 2136 1500 t (DO 60 J=I,N)2 528 1 2280 1600 t (IF \(I.NE.J\) D\(J\)=D\(J\) - DBLE\(SAVEC\(L\)\)*B\(I\))4 2064 1 2424 1700 t (D\(I\) = D\(I\) - DBLE\(SAVEC\(L\)\)*B\(J\))4 1584 1 2424 1800 t (60 L=L+1)1 768 1 1752 1900 t (R\(I\) = D\(I\))2 528 1 2280 2000 t (70 CONTINUE)1 768 1 1752 2100 t (C)1656 2200 w (C SOLVE A\(DELTAX\) =R)3 960 1 1656 2300 t (C)1656 2400 w (CALL SYFBS\(N,C,R,8,1,INTER\))1 1296 1 2136 2500 t (C)1656 2600 w (C DETERMINE NORM OF CORRECTION AND ADD IN CORRECTION)8 2496 1 1656 2700 t (C)1656 2800 w (WRITE\(IWRITE,71\)ITER)2136 2900 w ( THE SOLUTION AFTER ITERATION ,I5\))5 1632(71 FORMAT\(30H)1 864 2 1752 3000 t (RNORM=0.0)2136 3100 w (DO 80 I=1,N)2 528 1 2136 3200 t (B\(I\) = B\(I\) + R\(I\))4 864 1 2424 3300 t (RNORM=RNORM+ABS\(R\(I\)\))2424 3400 w (WRITE\(IWRITE,41\)B\(I\))2424 3500 w (80 CONTINUE)1 768 1 1752 3600 t (IF\(RNORM.LT.R1MACH\(4\)*BNORM\) STOP)1 1584 1 1992 3700 t (90 CONTINUE)1 624 1 1752 3800 t (WRITE\(IWRITE,91\))1992 3900 w ( ITERATIVE IMPROVEMENT FAILED\))3 1440(91 FORMAT\(29H)1 720 2 1752 4000 t (STOP)1992 4100 w (END)1992 4200 w 10 R f ( the)1 160(The above program was applied to a problem in which the upper triangular portion of)14 3584 2 1296 4560 t (symmetric matrix A was given by)5 1357 1 1296 4680 t 10 S f (-)2370 4920 w 10 R f (4.0 0.0)1 400 1 2425 4920 t 10 S f (-)2975 4920 w 10 R f (16.)3030 4920 w 10 S f (-)3355 4920 w 10 R f (32. 28.0)1 555 1 3410 4920 t ( 10.0)1 380(1.0 5.0)1 505 2 2700 5040 t 10 S f (-)3785 5040 w 10 R f (6.0)3840 5040 w 10 S f (-)2975 5160 w 10 R f (37.0)3030 5160 w 10 S f (-)3355 5160 w 10 R f (66.0 64.0)1 555 1 3410 5160 t 10 S f (-)3355 5280 w 10 R f (85.0 53.0)1 555 1 3410 5280 t 10 S f (-)3735 5400 w 10 R f (15.0)3790 5400 w ( conform to FORTRAN conventions)4 1478(Of course when the matrix was read in to the C array, to)12 2266 2 1296 5640 t (each of the above lines had to be left justi\256ed.)9 1836 1 1296 5760 t (When the following right hand side was read in)8 1897 1 1296 6000 t (448.)3108 6240 w 10 S f (-)3053 6360 w 10 R f (111.)3108 6360 w (1029.)3058 6480 w (1207.)3058 6600 w 10 S f (-)3053 6720 w 10 R f (719.)3108 6720 w (Linear Algebra)1 606 1 360 7500 t (\320 6 \320)2 300 1 2730 7620 t (SYCE)360 7740 w cleartomark showpage saveobj restore %%EndPage: 6 6 %%Page: 7 7 /saveobj save def mark 7 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYCE)1 4305(February 11, 1993)2 735 2 360 840 t (the following results were obtained on the Honeywell 6000 computer at Bell Labs:)12 3307 1 1296 1320 t 8 CW f ( 05)1 144( 6.71523875E)1 672(CONDITION NUMBER IS)2 912 3 1704 1660 t (THE FIRST SOLUTION X, FROM SYCE AND SYFBS=)7 2016 1 1704 1760 t (-8.0002890)2112 1860 w (-3.0003689)2112 1960 w (-1.9998967)2112 2060 w (-5.0000131)2112 2160 w (8.0000029)2160 2260 w ( 1)1 288(THE SOLUTION AFTER ITERATION)3 1344 2 1704 2360 t (-8.0000000)2112 2460 w (-3.0000000)2112 2560 w (-2.0000000)2112 2660 w (-5.0000000)2112 2760 w (8.0000000)2160 2860 w ( 2)1 288(THE SOLUTION AFTER ITERATION)3 1344 2 1704 2960 t (-8.0000000)2112 3060 w (-3.0000000)2112 3160 w (-2.0000000)2112 3260 w (-5.0000000)2112 3360 w (8.0000000)2160 3460 w 10 R f ( change in)2 424(As in most iterative algorithms, the algorithm implemented above stops when the)11 3320 2 1296 3800 t ( of iteration 1 is correct,)5 1004( the solution at the end)5 958( Although)1 438(the solution is suf\256ciently small.)4 1344 4 1296 3920 t ( hence the program decided to take one)7 1628(the change from the original solution was large and)8 2116 2 1296 4040 t (more step.)1 416 1 1296 4160 t ( the estimate of the)4 784(The \256rst solution above is inaccurate, as would have been expected from)11 2960 2 1296 4400 t ( iterative re\256nement algorithm successfully improved)5 2184( The)1 214( matrix.)1 320(condition number for the)3 1026 4 1296 4520 t ( be exactly rep-)3 632(the solution to this problem because the matrix and the right-hand side could)12 3112 2 1296 4640 t ( the input matrix)3 682( Often)1 283( the condition number was not high.\))6 1510( \(Also)1 278(resented in the machine.)3 991 5 1296 4760 t (cannot be represented exactly and the iterative re\256nement algorithm produces a very accu-)12 3744 1 1296 4880 t (rate, but worthless, solution to a slightly incorrect problem.)8 2368 1 1296 5000 t (Linear Algebra)1 606 1 4794 7500 t (\320 7 \320)2 300 1 2730 7620 t (SYCE)5144 7740 w cleartomark showpage saveobj restore %%EndPage: 7 7 %%Page: 8 8 /saveobj save def mark 8 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYCE February)1 4665 2 360 840 t (SYDC \320 decomposition of a symmetric matrix)6 1916 1 2210 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( the MDM)2 424(SYDC \(SYmmetric matrix DeComposition\) forms)4 2022 2 1296 1680 t 7 R f (T)3747 1640 w 10 R f (decomposition of a symmetric)3 1216 1 3824 1680 t ( as the \256rst step of the so-)7 1038(matrix A, which need not be positive de\256nite. It is called by SYLE)12 2706 2 1296 1800 t (lution of a symmetric linear system.)5 1438 1 1296 1920 t 9 B f (Usage:)720 2280 w 10 R f (CALL SYDC \(N, C, INTER\))4 1177 1 1296 2280 t (N)1440 2520 w 14 S f (\256)1980 2540 w 10 R f (the order of the matrix A)5 995 1 2196 2520 t (C)1440 2760 w 14 S f (\256)1980 2780 w 10 R f ( lower tri-)2 417(a one-dimensional array of length N\(N+1\)/2 into which the)8 2427 2 2196 2760 t ( A is packed by columns as illustrated in the)9 1804(angular part of the matrix)4 1040 2 2196 2880 t (following 4)1 464 1 2196 3000 t 10 S f (\264)2668 3000 w 10 R f (4 example:)1 441 1 2731 3000 t 10 I f (a)2232 3240 w 7 R f (11)2293 3260 w 10 I f (c)3312 3240 w 7 R f (1)3367 3260 w 10 I f (a)2232 3360 w 7 R f (21)2293 3380 w 10 I f (a)2412 3360 w 7 R f (22)2473 3380 w 14 S f (\256)2952 3360 w 10 I f (c)3312 3360 w 7 R f (2)3367 3380 w 10 I f (c)3451 3360 w 7 R f (5)3506 3380 w 10 I f (a)2232 3480 w 7 R f (31)2293 3500 w 10 I f (a)2412 3480 w 7 R f (32)2473 3500 w 10 I f (a)2592 3480 w 7 R f (33)2653 3500 w 10 I f (c)3312 3480 w 7 R f (3)3367 3500 w 10 I f (c)3451 3480 w 7 R f (6)3506 3500 w 10 I f (c)3590 3480 w 7 R f (8)3645 3500 w 10 I f (a)2232 3600 w 7 R f (41)2293 3620 w 10 I f (a)2412 3600 w 7 R f (42)2473 3620 w 10 I f (a)2592 3600 w 7 R f (43)2653 3620 w 10 I f (a)2772 3600 w 7 R f (44)2833 3620 w 10 I f (c)3312 3600 w 7 R f (4)3367 3620 w 10 I f (c)3451 3600 w 7 R f (7)3506 3620 w 10 I f (c)3590 3600 w 7 R f (9)3645 3620 w 10 I f (c)3729 3600 w 7 R f (10)3784 3620 w 14 S f (\254)1980 3860 w 10 R f (the D and M matrices of the decomposition for SYFBS \(see)10 2392 1 2196 3840 t 8 B f (Note 1)1 219 1 4613 3840 t 10 R f (\))4832 3840 w (INTER)1440 4080 w 14 S f (\254)1980 4100 w 10 R f ( the interchanges,)2 719(an integer vector of length N containing a record of)9 2125 2 2196 4080 t (i.e. the matrix P, described in)5 1171 1 2196 4200 t 8 B f (Note 1)1 219 1 3392 4200 t 10 R f (below.)3636 4200 w 9 B f (Note 1:)1 278 1 720 4560 t 10 R f (The MDM)1 437 1 1296 4560 t 7 R f (T)1738 4520 w 10 R f (decomposition of a symmetric matrix A satis\256es P)7 2068 1 1821 4560 t 7 R f (T)3894 4520 w 10 R f (AP = MDM)2 498 1 3945 4560 t 7 R f (T)4448 4520 w 10 R f (where P is a)3 509 1 4531 4560 t ( D is a block diagonal matrix)6 1232(permutation matrix, M is a unit lower triangular matrix, and)9 2512 2 1296 4680 t (with blocks of order 1)4 922 1 1296 4800 t 10 S f (\264)2226 4800 w 10 R f (1 and blocks of order 2)5 973 1 2289 4800 t 10 S f (\264)3270 4800 w 10 R f (2. Whenever)1 544 1 3333 4800 t 10 I f (d)3912 4800 w 7 I f (i)3973 4820 w 7 S f (+)4009 4820 w 7 R f (1 ,)1 58 1 4059 4820 t 7 I f (i)4122 4820 w 10 R f ( 2)1 86(is nonzero \(in a)3 648 2 4185 4800 t 10 S f (\264)4927 4800 w 10 R f (2)4990 4800 w (block of D\),)2 491 1 1296 4920 t 10 I f (m)1815 4920 w 7 I f (i)1898 4940 w 7 S f (+)1934 4940 w 7 R f (1 ,)1 58 1 1984 4940 t 7 I f (i)2047 4940 w 10 R f ( return from SYDC,)3 805( On)1 174(is zero.)1 291 3 2103 4920 t 10 I f (d)3400 4920 w 7 I f (ii)3461 4940 w 10 R f (, the diagonal of D, occupies the posi-)7 1531 1 3509 4920 t (tion of C which contained)4 1062 1 1296 5040 t 10 I f (a)2389 5040 w 7 I f (ii)2450 5060 w 10 R f ( the elements of the strictly lower portions of M)9 1975(on entry, and)2 536 2 2529 5040 t ( C. Since the diagonal elements of M are)8 1628(and D appear permuted in the remaining positions of)8 2116 2 1296 5160 t ( INTER contain information for construct-)5 1710( positive elements of)3 833( The)1 206(all 1, they are not stored.)5 995 4 1296 5280 t ( indicate)1 343(ing P \(see the introduction to this chapter\). The negative elements of INTER, if any,)14 3401 2 1296 5400 t (the presence of 2)3 687 1 1296 5520 t 10 S f (\264)1991 5520 w 10 R f ( contains a 2)3 514(2 blocks in D. If INTER\(I\) is negative, D)8 1671 2 2054 5520 t 10 S f (\264)4247 5520 w 10 R f (2 block beginning)2 730 1 4310 5520 t (at row I)2 310 1 1296 5640 t 10 S f (-)1606 5640 w 10 R f ( this case,)2 391(1. In)1 208 2 1661 5640 t 10 I f (d)2285 5640 w 7 I f (i)2346 5660 w 7 R f (,)2371 5660 w 7 I f (i)2394 5660 w 7 S f (-)2430 5660 w 7 R f (1)2480 5660 w 10 R f (directly follows)1 630 1 2548 5640 t 10 I f (d)3203 5640 w 7 I f (i)3264 5660 w 7 S f (-)3300 5660 w 7 R f (1 ,)1 58 1 3350 5660 t 7 I f (i)3413 5660 w 7 S f (-)3449 5660 w 7 R f (1)3499 5660 w 10 R f (in C.)1 195 1 3567 5640 t 9 B f (Note 2:)1 278 1 720 6000 t 10 R f ( is the conjugate transpose of A\), the)7 1541( *)1 58( where A)2 387( *,)1 83(For complex Hermitian matrices \(A = A)6 1675 5 1296 6000 t ( re-)1 145( decomposition and)2 805( *)1 58(complex Hermitian version of this subroutine computes the MDM)8 2736 4 1296 6120 t (turns the conjugate of M rather than M in C.)9 1770 1 1296 6240 t (Linear Algebra)1 606 1 360 7500 t (\320 8 \320)2 300 1 2730 7620 t (SYDC)360 7740 w cleartomark showpage saveobj restore %%EndPage: 8 8 %%Page: 9 9 /saveobj save def mark 9 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYDC)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Error situations:)1 653 1 720 1320 t 10 R f ( see)1 183( \320)1 125( those errors marked with an asterisk)6 1505(*\(The user can elect to `recover' from)6 1546 4 1517 1320 t 10 I f (Er-)4907 1320 w (ror Handling)1 531 1 1512 1440 t 10 R f (, Framework Chapter\))2 884 1 2043 1440 t (Number Error)1 1506 1 1800 1680 t (1 N)1 1008 1 1944 1920 t 10 S f (<)2977 1920 w 10 R f (1)3057 1920 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2160 t 9 B f (Double-precision version:)1 988 1 720 2520 t 10 R f (DSYDC with C declared double precision.)5 1709 1 1800 2520 t 9 B f (Complex symmetric version:)2 1106 1 720 2880 t 10 R f (CSYDC with C declared complex)4 1360 1 1901 2880 t 9 B f (Complex Hermitian version:)2 1101 1 720 3240 t 10 R f (CHEDC with C declared complex \(see)5 1550 1 1896 3240 t 8 B f (Note 2)1 219 1 3471 3240 t 10 R f (above\).)3715 3240 w 9 B f (Storage:)720 3600 w 10 R f (None)1296 3600 w 9 B f (Time:)720 4010 w 10 R f (6)1356 4080 w (N)1321 3950 w 7 R f (3)1398 3910 w 10 S1 f (___)1306 3980 w 10 S f (+)1515 4010 w 10 R f (4)1644 4080 w (3)1644 3950 w 10 S1 f (_ _)1 80 1 1629 3980 t 10 R f (N)1727 4010 w 7 R f (2)1804 3970 w 10 S f (+)1896 4010 w 10 R f (6)2025 4080 w (7)2025 3950 w 10 S1 f (_ _)1 80 1 2010 3980 t 10 R f (N additions)1 464 1 2108 4010 t (6)1356 4300 w (N)1321 4170 w 7 R f (3)1398 4130 w 10 S1 f (___)1306 4200 w 10 S f (+)1515 4230 w 10 R f (4)1679 4300 w (N)1644 4170 w 7 R f (2)1721 4130 w 10 S1 f (___)1629 4200 w 10 S f (+)1838 4230 w 10 R f (6)1992 4300 w (11)1967 4170 w 10 S1 f (_ __)1 130 1 1952 4200 t 10 R f (N multiplications)1 698 1 2100 4230 t (2)1356 4520 w (N)1321 4390 w 7 R f (2)1398 4350 w 10 S1 f (___)1306 4420 w 10 S f (+)1515 4450 w 10 R f (2)1655 4520 w (N)1644 4390 w 10 S1 f (_ __)1 102 1 1629 4420 t 10 R f (divisions)1766 4450 w (at most N)2 389 1 1296 4740 t 7 R f (2)1690 4700 w 10 S f (-)1749 4740 w 10 R f (1 comparisons)1 580 1 1820 4740 t 9 B f (Method:)720 5100 w 10 R f (The Bunch - Kaufman algorithm, described in reference [1] below, is used.)11 2998 1 1296 5100 t ( =)1 81( after setting EPS)3 730(SYDC calls SYMD)2 809 3 1296 5340 t 10 S f (\357\357)2941 5340 w 10 R f (A)3039 5340 w 10 S f (\357\357 e)1 167 1 3111 5340 t 10 R f (, where)1 304 1 3278 5340 t 10 S f (e)3618 5340 w 10 R f (is the machine precision, i.e. the)5 1342 1 3698 5340 t (value returned by R1MACH\(4\) \(or, for double precision, by D1MACH\(4\)\).)9 3022 1 1296 5460 t 9 B f (See also:)1 333 1 720 5820 t 10 R f (SYLE, SYFBS, SYMD, SYCE, SYSS)4 1542 1 1296 5820 t 9 B f (Author:)720 6180 w 10 R f (Linda Kaufman)1 629 1 1296 6180 t 9 B f (Reference:)720 6540 w 10 R f ( a symmetric matrix,)3 857(Bunch, J. R., Kaufman, L., and Parlett, B., Decomposition of)9 2521 2 1330 6540 t 10 I f (Numer.)4743 6540 w (Math 27)1 336 1 1296 6660 t 10 R f (\(1976\), 95-109.)1 624 1 1657 6660 t (Linear Algebra)1 606 1 4794 7500 t (\320 9 \320)2 300 1 2730 7620 t (SYDC)5133 7740 w cleartomark showpage saveobj restore %%EndPage: 9 9 %%Page: 10 10 /saveobj save def mark 10 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYDC February)1 4665 2 360 840 t 9 B f (Example:)720 1320 w 10 R f (The program fragment below determines how many positive eigenvalues a symmetric matrix,)11 3744 1 1296 1320 t (A, has.)1 280 1 1296 1440 t ( the eigenvalues of A correspond to the signs of)9 1923(According to the reference above the signs of)7 1821 2 1296 1680 t (the eigenvalues of D in the MDM)6 1420 1 1296 1800 t 7 R f (T)2721 1760 w 10 R f ( each 2)2 306( Moreover,)1 480(decomposition of A.)2 843 3 2809 1800 t 10 S f (\264)4446 1800 w 10 R f ( in D)2 222(2 block)1 309 2 4509 1800 t ( the number of positive eigenvalues)5 1443( Thus)1 254( positive-negative eigenvalue pair.)3 1387(corresponds to a)2 660 4 1296 1920 t (of D is equal to the number of 2)8 1276 1 1296 2040 t 10 S f (\264)2580 2040 w 10 R f (2 blocks of D plus the number of positive 1)9 1735 1 2643 2040 t 10 S f (\264)4386 2040 w 10 R f (1 blocks.)1 361 1 4449 2040 t (The program \256rst counts the number of 2)7 1659 1 1296 2280 t 10 S f (\264)2963 2280 w 10 R f ( D by counting the number of negative)7 1566(2 blocks of)2 448 2 3026 2280 t ( the array INTER, since each negative element \(see)8 2062(elements in)1 461 2 1296 2400 t 8 B f (Note 1)1 221 1 3846 2400 t 10 R f (above\) signals the exis-)3 946 1 4094 2400 t (tence of a 2)3 477 1 1296 2520 t 10 S f (\264)1781 2520 w 10 R f ( adds to this count the number of positive 1)9 1775( It)1 116(2 block.)1 327 3 1844 2520 t 10 S f (\264)4070 2520 w 10 R f ( D in order)3 453(1 blocks of)2 454 2 4133 2520 t (to \256nd the total number of positive eigenvalues of A.)9 2115 1 1296 2640 t 8 CW f (CALL SYDC\(N,C,INTER\))1 960 1 2040 2960 t (NPOS=0)2040 3060 w (C)1656 3160 w (C COUNT THE NUMBER OF POSITIVE EIGENVALUES OF D AND HENCE)10 2736 1 1656 3260 t (C OF THE MATRIX WHICH HAD ORIGINALLY BEEN PACKED INTO C)10 2640 1 1656 3360 t (C)1656 3460 w (C THE INDEX K PICKS OUT THE DIAGONAL OF THE MATRIX D OF)12 2640 1 1656 3560 t (C THE DECOMPOSITION)2 912 1 1656 3660 t (C)1656 3760 w (K=1)2040 3860 w (I=1)2040 3960 w ( \(I-N\) 20,30,40)2 720(10 IF)1 384 2 1752 4060 t ( .GT. 0\) GO TO 30)5 816(20 \(INTER\(I+1\))1 816 2 1752 4160 t (C)1656 4260 w (C OTHERWISE WE HAVE A 2)5 1104 1 1656 4360 t 8 S f (\264)2760 4360 w 8 CW f (2 BLOCK)1 336 1 2804 4360 t (C)1656 4460 w (NPOS=NPOS+1)2136 4560 w (K=K+2*\(N-I\)+1)2136 4660 w (I=I+2)2136 4760 w (GO TO 10)2 384 1 2136 4860 t (C)1656 4960 w (C WE HAVE A 1)4 624 1 1656 5060 t 8 S f (\264)2280 5060 w 8 CW f (1 BLOCK)1 336 1 2324 5060 t (C)1656 5160 w ( .GT. 0.0\) NPOS=NPOS+1)3 1056(30 IF\(C\(K\))1 720 2 1752 5260 t (K = K + N - I)6 624 1 2136 5360 t (I=I+1)2136 5460 w (GO TO 10)2 384 1 2136 5560 t (40 IWRITE=I1MACH\(2\))1 1056 1 1752 5660 t (WRITE\(IWRITE,41\)NPOS)2040 5760 w ( THE NUMBER OF POSITIVE EIGENVALUES IS,I5\))6 2016(41 FORMAT\(38H)1 768 2 1752 5860 t 10 R f (Linear Algebra)1 606 1 360 7500 t (\320 10 \320)2 350 1 2705 7620 t (SYDC)360 7740 w cleartomark showpage saveobj restore %%EndPage: 10 10 %%Page: 11 11 /saveobj save def mark 11 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYDC)1 4305(February 11, 1993)2 735 2 360 840 t (SYFBS \320 forward and back solve for symmetric matrices)8 2341 1 1997 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( matrix Forward and Back Solution\) solves AX = B where A is a sym-)14 2912(SYFBS \(SYmmetric)1 832 2 1296 1680 t ( is)1 104( It)1 123( or SYMD.)2 471(metric matrix using the decomposition of A computed by SYCE, SYDC,)10 3046 4 1296 1800 t (called by both SYSS and SYLE to solve symmetric linear systems.)10 2680 1 1296 1920 t 9 B f (Usage:)720 2280 w 10 R f (CALL SYFBS \(N, C, B, IB, NB, INTER\))7 1673 1 1296 2280 t (N)1440 2520 w 14 S f (\256)1980 2540 w 10 R f (the number of equations)3 968 1 2196 2520 t (C)1440 2760 w 14 S f (\256)1980 2780 w 10 R f ( of length N\(N+1\)/2 containing the matrices M)7 1880(a one-dimensional array)2 964 2 2196 2760 t (and D computed by SYDC, SYCE, or SYMD)7 1830 1 2196 2880 t (B)1440 3120 w 14 S f (\256)1980 3140 w 10 R f ( right-hand sides, dimensioned \(IB,KB\) in the calling pro-)8 2328(the matrix of)2 516 2 2196 3120 t (gram, where IB)2 623 1 2196 3240 t 10 S f (\263)2819 3240 w 10 R f (N and KB)2 405 1 2874 3240 t 10 S f (\263)3279 3240 w 10 R f (NB)3334 3240 w 14 S f (\254)1980 3500 w 10 R f (the solution X)2 567 1 2196 3480 t (IB)1440 3720 w 14 S f (\256)1980 3740 w 10 R f ( dimension of B, as dimensioned in the calling pro-)9 2139(the row \(leading\))2 705 2 2196 3720 t (gram)2196 3840 w (NB)1440 4080 w 14 S f (\256)1980 4100 w 10 R f (the number of right-hand sides)4 1226 1 2196 4080 t (INTER)1440 4320 w 14 S f (\256)1980 4340 w 10 R f ( a record of the interchanges)5 1187(an integer vector of length N containing)6 1657 2 2196 4320 t (performed by SYDC, SYCE, or SYMD)5 1585 1 2196 4440 t 9 B f (Error situations:)1 653 1 720 4800 t 10 R f ( see)1 183( \320)1 125( those errors marked with an asterisk)6 1505(*\(The user can elect to `recover' from)6 1546 4 1517 4800 t 10 I f (Er-)4907 4800 w (ror Handling)1 531 1 1512 4920 t 10 R f (, Framework Chapter\))2 884 1 2043 4920 t (Number Error)1 1506 1 1800 5160 t (1 N)1 1008 1 1944 5400 t 10 S f (<)2977 5400 w 10 R f (1)3057 5400 w (2 IB)1 1036 1 1944 5640 t 10 S f (<)3005 5640 w 10 R f (N)3085 5640 w (3 NB)1 1075 1 1944 5880 t 10 S f (<)3044 5880 w 10 R f (1)3124 5880 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 6120 t 9 B f (Double-precision version:)1 988 1 720 6480 t 10 R f (DSYFBS with C and B declared double precision.)7 2010 1 1800 6480 t (Linear Algebra)1 606 1 4794 7500 t (\320 11 \320)2 350 1 2705 7620 t (SYFBS)5093 7740 w cleartomark showpage saveobj restore %%EndPage: 11 11 %%Page: 12 12 /saveobj save def mark 12 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYFBS February)1 4665 2 360 840 t 9 B f (Complex version:)1 678 1 720 1320 t 10 R f (CSYFBS with C and B declared complex)6 1661 1 1473 1320 t 9 B f (Complex Hermitian version:)2 1101 1 720 1680 t 10 R f (CHEFBS with C and B declared complex)6 1666 1 1896 1680 t 9 B f (Storage:)720 2040 w 10 R f (None)1296 2040 w 9 B f (Time:)720 2400 w 10 R f (NB)1296 2400 w 10 S f (\264)1443 2400 w 10 R f (\( N)1 113 1 1506 2400 t 7 R f (2)1624 2360 w 10 S f (-)1683 2400 w 10 R f ( additions)1 392(N \))1 113 2 1754 2400 t (NB)1296 2520 w 10 S f (\264)1443 2520 w 10 R f (\( N)1 113 1 1506 2520 t 7 R f (2)1624 2480 w 10 S f (+)1683 2520 w 10 R f ( multiplications)1 626(N \))1 113 2 1754 2520 t (NB)1296 2640 w 10 S f (\264)1435 2640 w 10 R f (N divisions)1 459 1 1490 2640 t 9 B f (Method:)720 3000 w 10 R f (The Bunch - Kaufman algorithm, described in the reference below, is used.)11 3004 1 1296 3000 t 9 B f (See also:)1 333 1 720 3360 t 10 R f (SYLE, SYDC, SYMD, SYSS, SYCE)4 1502 1 1296 3360 t 9 B f (Author:)720 3720 w 10 R f (Linda Kaufman)1 629 1 1296 3720 t 9 B f (Reference:)720 4080 w 10 R f ( a symmetric matrix,)3 857(Bunch, J. R., Kaufman, L., and Parlett, B., Decomposition of)9 2521 2 1330 4080 t 10 I f (Numer.)4743 4080 w (Math 27)1 336 1 1296 4200 t 10 R f (\(1976\), 95-109.)1 624 1 1657 4200 t (Linear Algebra)1 606 1 360 7500 t (\320 12 \320)2 350 1 2705 7620 t (SYFBS)360 7740 w cleartomark showpage saveobj restore %%EndPage: 12 12 %%Page: 13 13 /saveobj save def mark 13 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYFBS)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Example:)720 1320 w 10 R f ( below replaces the K)4 872(The program fragment)2 905 2 1296 1320 t 10 S f (\264)3073 1320 w 10 R f (N matrix B with BA)4 825 1 3128 1320 t 7 S f (-)3964 1280 w 7 R f (1)4014 1280 w 10 R f (where A is a symmetric)4 956 1 4084 1320 t ( Note)1 255( the parameter list of SYDC.)5 1198(matrix packed into C according to the scheme given in)9 2291 3 1296 1440 t (that A)1 247 1 1296 1560 t 7 S f (-)1554 1520 w 7 R f (1)1604 1520 w 10 R f (is not formed explicitly since forming BA)6 1677 1 1672 1560 t 7 S f (-)3360 1520 w 7 R f (1)3410 1520 w 10 R f ( the)1 148(is equivalent to solving XA = B for)7 1414 2 3478 1560 t ( A is symmetric, solving XA = B is, in turn, equivalent to solving)13 2930( Because)1 406(matrix X.)1 408 3 1296 1680 t (AX)1296 1800 w 7 R f (T)1445 1760 w 10 S f (=)1536 1800 w 10 R f (B)1631 1800 w 7 R f (T)1703 1760 w 10 R f ( solving a linear system with `K' right-hand sides,)8 2049( the problem reduces to)4 957(. Thus)1 280 3 1754 1800 t ( which unfortunately resides in a `row' of the array B, rather than a column of an ar-)17 3448(each of)1 296 2 1296 1920 t ( B but to invoke SYFBS K)6 1080( the program fragment we chose not to transpose the matrix)10 2379(ray. In)1 285 3 1296 2040 t (times and store each row of B temporarily in the one-dimensional array Y.)12 2975 1 1296 2160 t 8 CW f (CALL SYDC\(N,C,INTER\))1 960 1 1992 2480 t (DO 30 I=1,K)2 528 1 1992 2580 t (DO 10 J=1,N)2 528 1 2136 2680 t (Y\(J\)=B\(I,J\))2280 2780 w (10 CONTINUE)1 864 1 1656 2880 t (CALL SYFBS\(N,C,Y,N,1,INTER\))1 1296 1 2136 2980 t (DO 20 J=1,N)2 528 1 2136 3080 t (B\(I,J\)=Y\(J\))2280 3180 w (20 CONTINUE)1 864 1 1656 3280 t (30 CONTINUE)1 720 1 1656 3380 t 10 R f (Linear Algebra)1 606 1 4794 7500 t (\320 13 \320)2 350 1 2705 7620 t (SYFBS)5093 7740 w cleartomark showpage saveobj restore %%EndPage: 13 13 %%Page: 14 14 /saveobj save def mark 14 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYFBS February)1 4665 2 360 840 t (SYLE \320 symmetric linear system solution)5 1725 1 2305 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( where A is a symmetric matrix.)6 1296( = B)2 173( Equation solution\) solves AX)4 1219(SYLE \(SYmmetric Linear)2 1056 4 1296 1680 t (A does not have to be positive de\256nite.)7 1560 1 1296 1800 t 9 B f (Usage:)720 2160 w 10 R f (CALL SYLE \(N, C, B, IB, NB\))6 1272 1 1296 2160 t (N)1440 2400 w 14 S f (\256)1980 2420 w 10 R f (the number of equations)3 968 1 2196 2400 t (C)1440 2640 w 14 S f (\256)1980 2660 w 10 R f ( lower tri-)2 417(a one-dimensional array of length N\(N+1\)/2 into which the)8 2427 2 2196 2640 t ( A is packed by columns as illustrated in the)9 1804(angular part of the matrix)4 1040 2 2196 2760 t (following 4)1 464 1 2196 2880 t 10 S f (\264)2668 2880 w 10 R f (4 example:)1 441 1 2731 2880 t 10 I f (a)2232 3120 w 7 R f (11)2293 3140 w 10 I f (c)3312 3120 w 7 R f (1)3367 3140 w 10 I f (a)2232 3240 w 7 R f (21)2293 3260 w 10 I f (a)2412 3240 w 7 R f (22)2473 3260 w 14 S f (\256)2952 3240 w 10 I f (c)3312 3240 w 7 R f (2)3367 3260 w 10 I f (c)3451 3240 w 7 R f (5)3506 3260 w 10 I f (a)2232 3360 w 7 R f (31)2293 3380 w 10 I f (a)2412 3360 w 7 R f (32)2473 3380 w 10 I f (a)2592 3360 w 7 R f (33)2653 3380 w 10 I f (c)3312 3360 w 7 R f (3)3367 3380 w 10 I f (c)3451 3360 w 7 R f (6)3506 3380 w 10 I f (c)3590 3360 w 7 R f (8)3645 3380 w 10 I f (a)2232 3480 w 7 R f (41)2293 3500 w 10 I f (a)2412 3480 w 7 R f (42)2473 3500 w 10 I f (a)2592 3480 w 7 R f (43)2653 3500 w 10 I f (a)2772 3480 w 7 R f (44)2833 3500 w 10 I f (c)3312 3480 w 7 R f (4)3367 3500 w 10 I f (c)3451 3480 w 7 R f (7)3506 3500 w 10 I f (c)3590 3480 w 7 R f (9)3645 3500 w 10 I f (c)3729 3480 w 7 R f (10)3784 3500 w 10 R f (C is overwritten during the solution.)5 1450 1 2196 3720 t (B)1440 3960 w 14 S f (\256)1980 3980 w 10 R f ( right-hand sides, dimensioned \(IB,KB\) in the calling pro-)8 2328(the matrix of)2 516 2 2196 3960 t (gram, where IB)2 623 1 2196 4080 t 10 S f (\263)2819 4080 w 10 R f (N and KB)2 405 1 2874 4080 t 10 S f (\263)3279 4080 w 10 R f (NB)3334 4080 w 14 S f (\254)1980 4340 w 10 R f (the solution X)2 567 1 2196 4320 t (IB)1440 4560 w 14 S f (\256)1980 4580 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 4560 t (calling program)1 635 1 2196 4680 t (NB)1440 4920 w 14 S f (\256)1980 4940 w 10 R f (the number of right-hand sides)4 1226 1 2196 4920 t 9 B f (Note 1:)1 278 1 720 5280 t 10 R f (Unless the given matrix A, is known in advance to be well-conditioned, the user should use)15 3744 1 1296 5280 t (the routine SYSS instead of SYLE.)5 1411 1 1296 5400 t 9 B f (Note 2:)1 278 1 720 5760 t 10 R f ( coef\256cient matrix, but differ-)4 1201(Users who wish to solve a sequence of problems with the same)11 2543 2 1296 5760 t (ent right-hand sides)2 796 1 1296 5880 t 10 I f (not all known in advance,)4 1046 1 2121 5880 t 10 R f ( use SYLE, but should call subpro-)6 1420(should not)1 424 2 3196 5880 t ( is called once to get the)6 1055( SYDC)1 332( the example in SYCE.\))4 1012( \(See)1 241(grams SYDC and SYFBS.)3 1104 5 1296 6000 t (MDM)1296 6120 w 7 R f (T)1551 6080 w 10 R f ( to this chapter\) and then SYFBS is called for)9 1926(decomposition \(see the introduction)3 1474 2 1640 6120 t (each new right-hand side.)3 1025 1 1296 6240 t (Linear Algebra)1 606 1 360 7500 t (\320 14 \320)2 350 1 2705 7620 t (SYLE)360 7740 w cleartomark showpage saveobj restore %%EndPage: 14 14 %%Page: 15 15 /saveobj save def mark 15 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYLE)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Error situations:)1 653 1 720 1320 t 10 R f ( see)1 183( \320)1 125( those errors marked with an asterisk)6 1505(*\(The user can elect to `recover' from)6 1546 4 1517 1320 t 10 I f (Er-)4907 1320 w (ror Handling)1 531 1 1512 1440 t 10 R f (, Framework Chapter\))2 884 1 2043 1440 t (Number Error)1 1506 1 1800 1680 t (1 N)1 1008 1 1944 1920 t 10 S f (<)2977 1920 w 10 R f (1)3057 1920 w (2 IB)1 1036 1 1944 2160 t 10 S f (<)3005 2160 w 10 R f (N)3085 2160 w (3 NB)1 1075 1 1944 2400 t 10 S f (<)3044 2400 w 10 R f (1)3124 2400 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2640 t 9 B f (Double-precision version:)1 988 1 720 3000 t 10 R f (DSYLE with C and B declared double precision.)7 1953 1 1783 3000 t 9 B f (Complex symmetric version:)2 1106 1 720 3360 t 10 R f (CSYLE with C and B declared complex)6 1604 1 1901 3360 t 9 B f (Complex Hermitian version:)2 1101 1 720 3720 t 10 R f (CHELE with C and B declared complex)6 1609 1 1896 3720 t 9 B f (Storage:)720 4080 w 10 R f (N integer locations of scratch storage in the dynamic storage stack)10 2650 1 1296 4080 t 9 B f (Time:)720 4490 w 10 R f (6)1356 4560 w (N)1321 4430 w 7 R f (3)1398 4390 w 10 S1 f (___)1306 4460 w 10 S f (+)1515 4490 w 10 R f (\()1619 4490 w (4)1685 4560 w (3)1685 4430 w 10 S1 f (_ _)1 80 1 1670 4460 t 10 S f (+)1776 4490 w 10 R f (NB \))1 180 1 1847 4490 t 10 S f (\264)2043 4490 w 10 R f (N)2106 4490 w 7 R f (2)2183 4450 w 10 S f (+)2275 4490 w 10 R f (\()2379 4490 w (6)2445 4560 w (7)2445 4430 w 10 S1 f (_ _)1 80 1 2430 4460 t 10 S f (+)2536 4490 w 10 R f ( additions)1 392( N)1 88(NB \))1 180 3 2607 4490 t (6)1356 4780 w (N)1321 4650 w 7 R f (3)1398 4610 w 10 S1 f (___)1306 4680 w 10 S f (+)1515 4710 w 10 R f (\()1619 4710 w (4)1685 4780 w (1)1685 4650 w 10 S1 f (_ _)1 80 1 1670 4680 t 10 S f (+)1776 4710 w 10 R f (NB \))1 180 1 1847 4710 t 10 S f (\264)2043 4710 w 10 R f (N)2106 4710 w 7 R f (2)2183 4670 w 10 S f (+)2275 4710 w 10 R f (\()2379 4710 w (6)2470 4780 w (11)2445 4650 w 10 S1 f (_ __)1 130 1 2430 4680 t 10 S f (+)2586 4710 w 10 R f (NB \))1 180 1 2657 4710 t 10 S f (\264)2853 4710 w 10 R f (N multiplications)1 698 1 2916 4710 t (2)1356 5000 w (N)1321 4870 w 7 R f (2)1398 4830 w 10 S1 f (___)1306 4900 w 10 S f (+)1515 4930 w 10 R f (\()1619 4930 w (2)1685 5000 w (1)1685 4870 w 10 S1 f (_ _)1 80 1 1670 4900 t 10 S f (+)1776 4930 w 10 R f (NB \))1 180 1 1847 4930 t 10 S f (\264)2043 4930 w 10 R f (N divisions)1 459 1 2106 4930 t (at most N)2 389 1 1296 5220 t 7 R f (2)1690 5180 w 10 S f (-)1749 5220 w 10 R f (1 comparisons)1 580 1 1820 5220 t 9 B f (Method:)720 5580 w 10 R f (The Bunch - Kaufman algorithm, described in the reference below, is used.)11 3004 1 1296 5580 t 9 B f (See also:)1 333 1 720 5940 t 10 R f (SYDC, SYFBS, SYMD, SYSS, SYCE)4 1559 1 1296 5940 t 9 B f (Author:)720 6300 w 10 R f (Linda Kaufman)1 629 1 1296 6300 t 9 B f (Reference:)720 6660 w 10 R f ( a symmetric matrix,)3 857(Bunch, J. R., Kaufman, L., and Parlett, B., Decomposition of)9 2521 2 1330 6660 t 10 I f (Numer.)4743 6660 w (Math 27)1 336 1 1296 6780 t 10 R f (\(1976\), 95-109.)1 624 1 1657 6780 t (Linear Algebra)1 606 1 4794 7500 t (\320 15 \320)2 350 1 2705 7620 t (SYLE)5150 7740 w cleartomark showpage saveobj restore %%EndPage: 15 15 %%Page: 16 16 /saveobj save def mark 16 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYLE February)1 4665 2 360 840 t 9 B f (Examples:)720 1320 w 10 R f ( following call to SYLE replaces the N)7 1556(1. The)1 335 2 1116 1440 t 10 S f (\264)3032 1440 w 10 R f ( B with A)3 395(K matrix)1 358 2 3112 1440 t 7 S f (-)3876 1400 w 7 R f (1)3926 1400 w 10 R f (B where A is a symmetric)5 1045 1 3995 1440 t ( that A)2 272( Note)1 244(matrix packed into a vector C.)5 1209 3 1296 1560 t 7 S f (-)3032 1520 w 7 R f (1)3082 1520 w 10 R f (is not computed)2 639 1 3150 1560 t 8 CW f (CALL SYLE \(N, C, B, N, K\))6 1200 1 2304 1900 t 10 R f ( the other hand, to compute the inverse A)8 1666(2. On)1 302 2 1116 2260 t 7 S f (-)3095 2220 w 7 R f (1)3145 2220 w 10 R f ( set B to the identity matrix and call)8 1457(, one may)2 395 2 3188 2260 t ( following program fragment leaves A)5 1543(SYLE. The)1 482 2 1296 2380 t 7 S f (-)3332 2340 w 7 R f (1)3382 2340 w 10 R f ( does not use the fact)5 845( It)1 112(in the matrix B.)3 631 3 3452 2380 t (that)1296 2500 w 10 I f (A)1471 2500 w 7 S f (-)1543 2460 w 7 R f (1)1593 2460 w 10 R f (is symmetric.)1 539 1 1661 2500 t 8 CW f (DO 20 I=1,N)2 528 1 2232 2840 t (DO 10 J=1,N)2 528 1 2376 2940 t (B\(I,J\)=0.0)2520 3040 w (10 CONTINUE)1 816 1 1944 3140 t (B\(I,I\)=1.0)2376 3240 w (20 CONTINUE)1 672 1 1944 3340 t (CALL SYLE\(N, C, B, N, N\))5 1152 1 2232 3440 t 10 R f (Linear Algebra)1 606 1 360 7500 t (\320 16 \320)2 350 1 2705 7620 t (SYLE)360 7740 w cleartomark showpage saveobj restore %%EndPage: 16 16 %%Page: 17 17 /saveobj save def mark 17 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYLE)1 4305(February 11, 1993)2 735 2 360 840 t (SYMD \320 MDM)2 689 1 2033 1320 t 7 R f (T)2727 1280 w 10 R f (decomposition of a symmetric matrix)4 1499 1 2803 1320 t 9 B f (Purpose:)720 1680 w 10 R f (SYMD \(SYmmetric MDM)2 1089 1 1296 1680 t 7 R f (T)2390 1640 w 10 R f (decomposition\) forms the decomposition PMDM)4 1988 1 2470 1680 t 7 R f (T)4463 1640 w 10 R f (P)4514 1680 w 7 R f (T)4575 1640 w 10 R f (of a sym-)2 385 1 4655 1680 t ( and D is)3 361(metric matrix A, where P is a permutation matrix, M is unit lower triangular matrix,)14 3383 2 1296 1800 t ( subroutine allows the user)4 1099( This)1 235(block diagonal. The matrix A need not be positive de\256nite.)9 2410 3 1296 1920 t ( is called by the decomposition)5 1276( It)1 118(to specify a threshold for considering the matrix singular.)8 2350 3 1296 2040 t (routines SYDC and SYCE.)3 1089 1 1296 2160 t 9 B f (Usage:)720 2520 w 10 R f (CALL SYMD \(N, C, INTER, EPS\))5 1422 1 1296 2520 t (N)1440 2760 w 14 S f (\256)1980 2780 w 10 R f (the order of the matrix A)5 995 1 2196 2760 t (C)1440 3000 w 14 S f (\256)1980 3020 w 10 R f ( lower tri-)2 417(a one-dimensional array of length N\(N+1\)/2 into which the)8 2427 2 2196 3000 t ( A is packed by columns as illustrated in the)9 1804(angular part of the matrix)4 1040 2 2196 3120 t (following 4)1 464 1 2196 3240 t 10 S f (\264)2660 3240 w 10 R f (4 example:)1 441 1 2715 3240 t 10 I f (a)2232 3480 w 7 R f (11)2293 3500 w 10 I f (c)3312 3480 w 7 R f (1)3367 3500 w 10 I f (a)2232 3600 w 7 R f (21)2293 3620 w 10 I f (a)2412 3600 w 7 R f (22)2473 3620 w 14 S f (\256)2952 3600 w 10 I f (c)3312 3600 w 7 R f (2)3367 3620 w 10 I f (c)3451 3600 w 7 R f (5)3506 3620 w 10 I f (a)2232 3720 w 7 R f (31)2293 3740 w 10 I f (a)2412 3720 w 7 R f (32)2473 3740 w 10 I f (a)2592 3720 w 7 R f (33)2653 3740 w 10 I f (c)3312 3720 w 7 R f (3)3367 3740 w 10 I f (c)3451 3720 w 7 R f (6)3506 3740 w 10 I f (c)3590 3720 w 7 R f (8)3645 3740 w 10 I f (a)2232 3840 w 7 R f (41)2293 3860 w 10 I f (a)2412 3840 w 7 R f (42)2473 3860 w 10 I f (a)2592 3840 w 7 R f (43)2653 3860 w 10 I f (a)2772 3840 w 7 R f (44)2833 3860 w 10 I f (c)3312 3840 w 7 R f (4)3367 3860 w 10 I f (c)3451 3840 w 7 R f (7)3506 3860 w 10 I f (c)3590 3840 w 7 R f (9)3645 3860 w 10 I f (c)3729 3840 w 7 R f (10)3784 3860 w 14 S f (\254)1980 4100 w 10 R f (the D and M matrices of the decomposition for SYFBS \(see)10 2392 1 2196 4080 t 8 B f (Note 1)1 219 1 4613 4080 t 10 R f (\))4832 4080 w (INTER)1440 4320 w 14 S f (\254)1980 4340 w 10 R f (a vector of length N containing a record of the interchanges,)10 2402 1 2196 4320 t (i.e. the matrix P, described in)5 1171 1 2196 4440 t 8 B f (Note 1)1 219 1 3392 4440 t 10 R f (below.)3636 4440 w (EPS)1440 4680 w 14 S f (\256)1980 4700 w 10 R f (if)2196 4680 w 10 S f (\357)2291 4680 w 10 I f (d)2375 4680 w 7 I f (kk)2436 4700 w 10 S f (\357 \243)1 139 1 2541 4680 t 10 R f (EPS and)1 352 1 2705 4680 t 10 I f (d)3092 4680 w 7 I f (kk)3153 4700 w 10 R f (corresponds to a 1)3 759 1 3258 4680 t 10 S f (\264)4017 4680 w 10 R f (1 block of D, then C is)6 968 1 4072 4680 t (considered a singular matrix whose rank is at least k)9 2088 1 2196 4800 t 10 S f (-)4284 4800 w 10 R f (1)4339 4800 w 9 B f (Note 1:)1 278 1 720 5160 t 10 R f (The MDM)1 437 1 1296 5160 t 7 R f (T)1738 5120 w 10 R f (decomposition of a symmetric matrix A satis\256es P)7 2068 1 1821 5160 t 7 R f (T)3894 5120 w 10 R f (AP = MDM)2 498 1 3945 5160 t 7 R f (T)4448 5120 w 10 R f (where P is a)3 509 1 4531 5160 t ( D is a block diagonal matrix)6 1232(permutation matrix, M is a unit lower triangular matrix, and)9 2512 2 1296 5280 t (with blocks of order 1)4 922 1 1296 5400 t 10 S f (\264)2226 5400 w 10 R f (1 and blocks of order 2)5 973 1 2289 5400 t 10 S f (\264)3270 5400 w 10 R f (2. Whenever)1 544 1 3333 5400 t 10 I f (d)3912 5400 w 7 I f (i)3973 5420 w 7 S f (+)4009 5420 w 7 R f (1 ,)1 58 1 4059 5420 t 7 I f (i)4122 5420 w 10 R f ( 2)1 86(is nonzero \(in a)3 648 2 4185 5400 t 10 S f (\264)4927 5400 w 10 R f (2)4990 5400 w (block of D\),)2 487 1 1296 5520 t 10 I f (m)1809 5520 w 7 I f (i)1892 5540 w 7 S f (+)1928 5540 w 7 R f (1 ,)1 58 1 1978 5540 t 7 I f (i)2041 5540 w 10 R f ( return from SYMD,)3 824( On)1 173(is zero.)1 289 3 2095 5520 t 10 I f (d)3407 5520 w 7 I f (ii)3468 5540 w 10 R f (, the diagonal of D, occupies the posi-)7 1524 1 3516 5520 t (tion of C which contained)4 1062 1 1296 5640 t 10 I f (a)2389 5640 w 7 I f (ii)2450 5660 w 10 R f ( the elements of the strictly lower portions of M)9 1975(on entry, and)2 536 2 2529 5640 t ( C. Since the diagonal elements of M are)8 1628(and D appear permuted in the remaining positions of)8 2116 2 1296 5760 t ( INTER contain information for construct-)5 1710( positive elements of)3 833( The)1 206(all 1, they are not stored.)5 995 4 1296 5880 t ( indicate)1 343(ing P \(see the introduction to this chapter\). The negative elements of INTER, if any,)14 3401 2 1296 6000 t (the presence of 2)3 687 1 1296 6120 t 10 S f (\264)1991 6120 w 10 R f ( contains a 2)3 514(2 blocks in D. If INTER\(I\) is negative, D)8 1671 2 2054 6120 t 10 S f (\264)4247 6120 w 10 R f (2 block beginning)2 730 1 4310 6120 t (at row I)2 310 1 1296 6240 t 10 S f (-)1606 6240 w 10 R f ( this case,)2 391(1. In)1 208 2 1661 6240 t 10 I f (d)2285 6240 w 7 I f (i)2346 6260 w 7 R f (,)2371 6260 w 7 I f (i)2394 6260 w 7 S f (-)2430 6260 w 7 R f (1)2480 6260 w 10 R f (directly follows)1 630 1 2548 6240 t 10 I f (d)3203 6240 w 7 I f (i)3264 6260 w 7 S f (-)3300 6260 w 7 R f (1 ,)1 58 1 3350 6260 t 7 I f (i)3413 6260 w 7 S f (-)3449 6260 w 7 R f (1)3499 6260 w 10 R f (in C.)1 195 1 3567 6240 t 9 B f (Note 2:)1 278 1 720 6600 t 10 R f ( is the conjugate transpose of A\), the)7 1562( *)1 58( A)1 111( where)1 281( *,)1 83(For complex Hermitian matrices\(A = A)5 1649 6 1296 6600 t ( re-)1 145( decomposition and)2 805( *)1 58(complex Hermitian version of this subroutine computes the MDM)8 2736 4 1296 6720 t (turns the conjugate of M rather than M in C.)9 1770 1 1296 6840 t (Linear Algebra)1 606 1 4794 7500 t (\320 17 \320)2 350 1 2705 7620 t (SYMD)5111 7740 w cleartomark showpage saveobj restore %%EndPage: 17 17 %%Page: 18 18 /saveobj save def mark 18 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYMD February)1 4665 2 360 840 t 9 B f (Error situations:)1 653 1 720 1320 t 10 R f ( see)1 183( \320)1 125( those errors marked with an asterisk)6 1505(*\(The user can elect to `recover' from)6 1546 4 1517 1320 t 10 I f (Er-)4907 1320 w (ror Handling)1 531 1 1512 1440 t 10 R f (, Framework Chapter\))2 884 1 2043 1440 t (Number Error)1 1506 1 1800 1680 t (1 N)1 1008 1 1944 1920 t 10 S f (<)2977 1920 w 10 R f (1)3057 1920 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2160 t 9 B f (Double-precision version:)1 988 1 720 2520 t 10 R f (DSYMD with C and EPS declared double precision.)7 2098 1 1800 2520 t 9 B f (Complex symmetric version:)2 1106 1 720 2880 t 10 R f (CSYMD with C declared complex)4 1382 1 1901 2880 t 9 B f (Complex Hermitian version:)2 1101 1 720 3240 t 10 R f (CHEMD with C declared complex \(see)5 1572 1 1896 3240 t 8 B f (Note 2)1 219 1 3493 3240 t 10 R f (\).)3712 3240 w 9 B f (Storage:)720 3600 w 10 R f (None)1296 3600 w 9 B f (Time:)720 4010 w 10 R f (6)1356 4080 w (N)1321 3950 w 7 R f (3)1398 3910 w 10 S1 f (___)1306 3980 w 10 S f (+)1515 4010 w 10 R f (4)1679 4080 w (N)1644 3950 w 7 R f (2)1721 3910 w 10 S1 f (___)1629 3980 w 10 S f (+)1838 4010 w 10 R f (6)1967 4080 w (7)1967 3950 w 10 S1 f (_ _)1 80 1 1952 3980 t 10 R f (N additions)1 464 1 2050 4010 t (6)1356 4360 w (N)1321 4230 w 7 R f (3)1398 4190 w 10 S1 f (___)1306 4260 w 10 S f (+)1515 4290 w 10 R f (4)1679 4360 w (N)1644 4230 w 7 R f (2)1721 4190 w 10 S1 f (___)1629 4260 w 10 S f (+)1838 4290 w 10 R f (6)1992 4360 w (11)1967 4230 w 10 S1 f (_ __)1 130 1 1952 4260 t 10 R f (N additions)1 464 1 2100 4290 t (2)1356 4640 w (N)1321 4510 w 7 R f (2)1398 4470 w 10 S1 f (___)1306 4540 w 10 S f (+)1515 4570 w 10 R f (2)1655 4640 w (N)1644 4510 w 10 S1 f (_ __)1 102 1 1629 4540 t 10 R f (divisions)1766 4570 w (at most N)2 389 1 1296 4776 t 7 R f (2)1690 4736 w 10 S f (-)1749 4776 w 10 R f (1 comparisons)1 580 1 1820 4776 t 9 B f (Method:)720 5136 w 10 R f (The Bunch - Kaufman algorithm, described in the reference below, is used.)11 3004 1 1296 5136 t 9 B f (See also:)1 333 1 720 5496 t 10 R f (SYLE, SYDC, SYFBS, SYSS, SYCE)4 1520 1 1296 5496 t 9 B f (Author:)720 5856 w 10 R f (Linda Kaufman)1 629 1 1296 5856 t (Linear Algebra)1 606 1 360 7500 t (\320 18 \320)2 350 1 2705 7620 t (SYMD)360 7740 w cleartomark showpage saveobj restore %%EndPage: 18 18 %%Page: 19 19 /saveobj save def mark 19 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYMD)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Reference:)720 1320 w 10 R f ( a symmetric matrix,)3 857(Bunch, J. R., Kaufman, L., and Parlett, B., Decomposition of)9 2521 2 1330 1320 t 10 I f (Numer.)4743 1320 w (Math 27)1 336 1 1296 1440 t 10 R f (\(1976\), 95-109.)1 624 1 1657 1440 t 9 B f (Example:)720 1800 w 10 R f (The following program fragment determines whether a matrix is positive de\256nite. According)11 3744 1 1296 1800 t ( is positive de\256nite only if D in)7 1248(to the theory given in the reference above, a symmetric matrix)10 2496 2 1296 1920 t ( computed by SYMD is diagonal with positive diagonal elements. If D is)12 3002(the decomposition)1 742 2 1296 2040 t ( of D are packed into C in)7 1062(diagonal, all the elements of INTER are positive and the elements)10 2682 2 1296 2160 t (the same positions that the diagonal of A had originally occupied.)10 2628 1 1296 2280 t 8 CW f (CALL SYMD\(N,C,INTER,0.0\))1 1152 1 1992 2600 t (IWRITE=I1MACH\(2\))1992 2700 w (C)1656 2800 w (C DETERMINE IF THE MATRIX PACKED INTO C IS POSITIVE DEFINITE.)10 2928 1 1656 2900 t (C THE INDEX K PICKS OUT THE DIAGONAL OF THE MATRIX D)11 2496 1 1656 3000 t (C OF THE DECOMPOSITION)3 1056 1 1656 3100 t (K=1)2040 3200 w (DO 10 I=1,N)2 528 1 2040 3300 t (IF\(INTER\(I\).LT.0.OR.C\(K\).LE.0.0\) GO TO 20)3 1968 1 2184 3400 t (C FIND NEXT DIAGONAL ELEMENT)4 1344 1 1656 3500 t (K = K + N - I)6 624 1 2184 3600 t (10 CONTINUE)1 768 1 1656 3700 t (WRITE\(IWRITE,11\))2040 3800 w ( THE MATRIX IS POSITIVE DEFINITE\))5 1584(11 FORMAT\(32H)1 864 2 1656 3900 t (GO TO 30)2 384 1 2040 4000 t (20 WRITE\(IWRITE,21\))1 1152 1 1656 4100 t ( THE MATRIX IS NOT POSITIVE DEFINITE\))6 1776(21 FORMAT\(36H)1 864 2 1656 4200 t (30 CONTINUE)1 768 1 1656 4300 t 10 R f (Linear Algebra)1 606 1 4794 7500 t (\320 19 \320)2 350 1 2705 7620 t (SYMD)5111 7740 w cleartomark showpage saveobj restore %%EndPage: 19 19 %%Page: 20 20 /saveobj save def mark 20 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYMD February)1 4665 2 360 840 t (SYML \320 symmetric matrix - vector multiplication)6 2055 1 2140 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( the product)2 489(SYML \(SYmmetric matrix MuLtiplication\) forms)4 2025 2 1296 1680 t 10 I f (Ax)3841 1680 w 10 R f (where)3977 1680 w 10 I f (A)4251 1680 w 10 R f (is a general sym-)3 697 1 4343 1680 t (metric matrix stored in packed form.)5 1464 1 1296 1800 t 9 B f (Usage:)720 2160 w 10 R f (CALL SYML\(N, C, X, B\))4 1058 1 1296 2160 t (N)1440 2400 w 14 S f (\256)1980 2420 w 10 R f (the length of)2 505 1 2196 2400 t 10 I f (x)2726 2400 w 10 R f (C)1440 2640 w 14 S f (\256)1980 2660 w 10 R f ( lower tri-)2 417(a one-dimensional array of length N\(N+1\)/2 into which the)8 2427 2 2196 2640 t ( A is packed by columns as illustrated in the)9 1804(angular part of the matrix)4 1040 2 2196 2760 t (following 4)1 464 1 2196 2880 t 10 S f (\264)2660 2880 w 10 R f (4 example:)1 441 1 2715 2880 t 10 I f (a)2232 3120 w 7 R f (11)2293 3140 w 10 I f (c)3312 3120 w 7 R f (1)3367 3140 w 10 I f (a)2232 3240 w 7 R f (21)2293 3260 w 10 I f (a)2412 3240 w 7 R f (22)2473 3260 w 14 S f (\256)2952 3240 w 10 I f (c)3312 3240 w 7 R f (2)3367 3260 w 10 I f (c)3451 3240 w 7 R f (5)3506 3260 w 10 I f (a)2232 3360 w 7 R f (31)2293 3380 w 10 I f (a)2412 3360 w 7 R f (32)2473 3380 w 10 I f (a)2592 3360 w 7 R f (33)2653 3380 w 10 I f (c)3312 3360 w 7 R f (3)3367 3380 w 10 I f (c)3451 3360 w 7 R f (6)3506 3380 w 10 I f (c)3590 3360 w 7 R f (8)3645 3380 w 10 I f (a)2232 3480 w 7 R f (41)2293 3500 w 10 I f (a)2412 3480 w 7 R f (42)2473 3500 w 10 I f (a)2592 3480 w 7 R f (43)2653 3500 w 10 I f (a)2772 3480 w 7 R f (44)2833 3500 w 10 I f (c)3312 3480 w 7 R f (4)3367 3500 w 10 I f (c)3451 3480 w 7 R f (7)3506 3500 w 10 I f (c)3590 3480 w 7 R f (9)3645 3500 w 10 I f (c)3729 3480 w 7 R f (10)3784 3500 w 10 R f (X)1440 3720 w 14 S f (\256)1980 3740 w 10 R f (the vector)1 396 1 2196 3720 t 10 I f (x)2617 3720 w 10 R f (to be multiplied)2 634 1 2686 3720 t (B)1440 3960 w 14 S f (\254)1980 3980 w 10 R f (the vector A)2 493 1 2196 3960 t 10 I f (x)2689 3960 w 9 B f (Error situations:)1 648 1 720 4320 t 10 R f (\(All errors in this subprogram are fatal \320)7 1666 1 1512 4320 t (see)1512 4440 w 10 I f (Error Handling)1 631 1 1664 4440 t 10 R f (, Framework Chapter\))2 884 1 2295 4440 t (Number Error)1 1506 1 1800 4680 t (1 N)1 1108 1 1944 4920 t 10 S f (<)3077 4920 w 10 R f (1)3157 4920 w 9 B f (Double-precision version:)1 988 1 720 5280 t 10 R f (DSYML with C, X, and B declared double precision.)8 2128 1 1800 5280 t 9 B f (Complex version:)1 678 1 720 5640 t 10 R f (CSYML with C, X, and B declared complex)7 1779 1 1473 5640 t 9 B f (Complex Hermitian version:)2 1101 1 720 6000 t 10 R f (CHEML with C, X, and B declared complex)7 1784 1 1896 6000 t 9 B f (Time:)720 6360 w 10 R f (N)1296 6360 w 7 R f (2)1373 6320 w 10 R f (additions)1441 6360 w (N)1296 6480 w 7 R f (2)1373 6440 w 10 R f (multiplications)1441 6480 w (Linear Algebra)1 606 1 360 7500 t (\320 20 \320)2 350 1 2705 7620 t (SYML)360 7740 w cleartomark showpage saveobj restore %%EndPage: 20 20 %%Page: 21 21 /saveobj save def mark 21 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYML)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (See also:)1 333 1 720 1320 t 10 R f (SYFBS, SYCE, SYDC, SYMD, SYLE, SYSS)5 1859 1 1296 1320 t 9 B f (Author:)720 1680 w 10 R f (Linda Kaufman)1 629 1 1296 1680 t 9 B f (Example:)720 2040 w 10 R f ( and SYLE, the symmetric linear equation)6 1774(This example checks the consistency of SYML)6 1970 2 1296 2040 t (solver.)1296 2160 w (First the example uses SYML to compute for a given vector)10 2397 1 1296 2400 t 10 I f (x)3718 2400 w 10 R f (and matrix)1 430 1 3787 2400 t 10 I f (A)4242 2400 w 10 R f (, the vector)2 446 1 4303 2400 t 10 I f (b)2736 2640 w 10 S f (=)2835 2640 w 10 I f (Ax .)1 138 1 2939 2640 t 10 R f (Then the problem is inverted, i.e., SYLE is used to \256nd the vector)12 2631 1 1296 2880 t 10 I f (x)3952 2880 w 10 R f (which satis\256es)1 586 1 4021 2880 t 10 I f (Ax)2736 3120 w 10 S f (=)2890 3120 w 10 I f (b)2994 3120 w 10 R f (This)1296 3360 w 10 I f (x)1506 3360 w 10 R f ( 10)1 132( The)1 212(is then compared with the original vector.)6 1709 3 1582 3360 t 10 S f (\264)3635 3360 w 10 R f (10 symmetric matrix)2 847 1 3690 3360 t 10 I f (A)4569 3360 w 10 R f (is chosen)1 377 1 4663 3360 t (so that)1 264 1 1296 3480 t 10 I f (a)2592 3720 w 7 I f (i j)1 45 1 2653 3740 t 10 S f (= \357)1 153 1 2755 3720 t 10 I f (i)2916 3720 w 10 S f (-)2968 3720 w 10 I f (j)3039 3720 w 10 S f (\357)3075 3720 w 10 I f (.)3132 3720 w 10 R f (The vector)1 429 1 1296 3960 t 10 I f (x)1750 3960 w 10 R f (is chosen randomly.)2 802 1 1819 3960 t 8 CW f (INTEGER N, L, I, J, IWRITE, I1MACH)6 1632 1 2040 4300 t (REAL C\(55\), X\(10\), B\(10\))3 1152 1 2040 4400 t (REAL UNI, ERR, SASUM, ABS)4 1200 1 2040 4500 t (N=10)2040 4600 w (C)1656 4700 w (C CONSTRUCT THE MATRIX A\(I,J\)=ABS\(J-I\) AND PACK INTO C)8 2592 1 1656 4800 t (C)1656 4900 w (L=0)2040 5000 w (DO 20 I=1,N)2 528 1 2040 5100 t (DO 10 J=I,N)2 528 1 2184 5200 t (L=L+1)2328 5300 w (C\(L\)=J-I)2328 5400 w (10 CONTINUE)1 816 1 1752 5500 t (20 CONTINUE)1 672 1 1752 5600 t (C)1656 5700 w (C CONSTRUCT A RANDOM VECTOR X)5 1392 1 1656 5800 t (C)1656 5900 w (DO 30 I=1,N)2 528 1 2040 6000 t (X\(I\)=UNI\(0\))2184 6100 w (30 CONTINUE)1 672 1 1752 6200 t (C)1656 6300 w (C FIND THE VECTOR B=AX)4 1056 1 1656 6400 t (C)1656 6500 w (CALL SYML\(N,C,X,B\))1 864 1 1992 6600 t (C)1656 6700 w (C SOLVE THE SYSTEM AX=B)4 1104 1 1656 6800 t (C)1656 6900 w 10 R f (Linear Algebra)1 606 1 4794 7560 t (\320 21 \320)2 350 1 2705 7680 t (SYML)5122 7800 w cleartomark showpage saveobj restore %%EndPage: 21 21 %%Page: 22 22 /saveobj save def mark 22 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYML February)1 4665 2 360 840 t 8 CW f (CALL SYLE\(N,C,B,N,1\))1 960 1 1992 1300 t (C)1656 1400 w (C PRINT THE COMPUTED AND TRUE SOLUTION)6 1824 1 1656 1500 t (C)1656 1600 w (IWRITE=I1MACH\(2\))1992 1700 w (WRITE\(IWRITE,31\))1992 1800 w ( SOLUTION\))1 480( COMPUTED)1 528( TRUE SOLUTION)2 672(31 FORMAT\(34H)1 720 4 1752 1900 t (WRITE\(IWRITE,32\)\(X\(I\),B\(I\),I=1,N\))1992 2000 w ( ,2E17.8\))1 432(32 FORMAT\(1H)1 672 2 1752 2100 t (C)1656 2200 w (C COMPUTE THE RELATIVE ERROR)4 1344 1 1656 2300 t (C)1656 2400 w (ERR=0.0)1992 2500 w (DO 40 I=1,N)2 528 1 1992 2600 t (ERR=ERR+ABS\(B\(I\)-X\(I\)\))2136 2700 w (40 CONTINUE)1 624 1 1752 2800 t (ERR=ERR/SASUM\(N,X,1\))1992 2900 w (WRITE\(IWRITE,41\)ERR)1992 3000 w ( RELATIVE ERROR IS ,1PE15.7\))4 1344(41 FORMAT\(19H)1 720 2 1752 3100 t (STOP)1992 3200 w (END)1992 3300 w 10 R f ( on the Honeywell 6000 machine at Bell Laboratories,)8 2174(When the above program was executed)5 1570 2 1296 3640 t (the following was printed:)3 1052 1 1296 3760 t 8 CW f ( SOLUTION)1 432( COMPUTED)1 816(TRUE SOLUTION)1 624 3 1848 4100 t ( 00)1 144( 0.22925607E)1 1008(0.22925607E 00)1 672 3 1848 4200 t ( 00)1 144( 0.76687498E)1 1008(0.76687502E 00)1 672 3 1848 4300 t ( 00)1 144( 0.68317696E)1 1008(0.68317685E 00)1 672 3 1848 4400 t ( 00)1 144( 0.50919100E)1 1008(0.50919111E 00)1 672 3 1848 4500 t ( 00)1 144( 0.87455969E)1 1008(0.87455959E 00)1 672 3 1848 4600 t ( 00)1 144( 0.64464090E)1 1008(0.64464101E 00)1 672 3 1848 4700 t ( 00)1 144( 0.84746833E)1 1008(0.84746840E 00)1 672 3 1848 4800 t ( 00)1 144( 0.35396342E)1 1008(0.35396343E 00)1 672 3 1848 4900 t ( 00)1 144( 0.39889174E)1 1008(0.39889160E 00)1 672 3 1848 5000 t ( 00)1 144( 0.45709418E)1 1008(0.45709422E 00)1 672 3 1848 5100 t ( 1.2794319E-07)1 768(RELATIVE ERROR IS)2 816 2 1704 5200 t 10 R f ( and the machine)3 691(The condition number of the matrix \(see the example in SYSS\) is about 54.)13 3053 2 1296 5540 t ( is about 10)3 482(precision, on the Honeywell computer)4 1556 2 1296 5660 t 7 S f (-)3345 5620 w 7 R f (8)3395 5620 w 10 R f ( even in the absence of roundoff)6 1321(. Thus)1 281 2 3438 5660 t (error in SYML, a relative error of 5)7 1509 1 1296 5780 t 10 S f (\264)2813 5780 w 10 R f (10)2876 5780 w 7 S f (-)2987 5740 w 7 R f (7)3037 5740 w 10 R f ( value computed)2 688( The)1 219(would not be surprising.)3 1014 3 3119 5780 t (above is quite reasonable.)3 1031 1 1296 5900 t (Linear Algebra)1 606 1 360 7500 t (\320 22 \320)2 350 1 2705 7620 t (SYML)360 7740 w cleartomark showpage saveobj restore %%EndPage: 22 22 %%Page: 23 23 /saveobj save def mark 23 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYML)1 4305(February 11, 1993)2 735 2 360 840 t ( of a symmetric matrix)4 910( norm)1 261(SYNM \320)1 414 3 2375 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( the norm of a symmetric matrix A stored in)9 1888(SYNM \(SYmmetric matrix NorM\) computes)4 1856 2 1296 1680 t (packed form. The in\256nity norm is de\256ned as)7 1776 1 1296 1850 t 7 R f (1)3097 1920 w 7 S f (\243)3143 1920 w 7 I f (i)3198 1920 w 7 S f (\243)3234 1920 w 7 I f (n)3284 1920 w 10 R f (max)3122 1850 w 7 I f (j)3329 1950 w 7 S f (=)3360 1950 w 7 R f (1)3410 1950 w 15 S f (S)3343 1880 w 7 I f (n)3370 1750 w 10 S f (\357)3456 1850 w 10 I f (a)3513 1850 w 7 I f (i j)1 45 1 3574 1870 t 10 S f (\357)3659 1850 w 9 B f (Type:)720 2290 w 10 R f (Real function)1 541 1 1296 2290 t 9 B f (Usage:)720 2650 w 10 S f (<)1296 2650 w 10 R f (answer)1351 2650 w 10 S f (>)1633 2650 w 10 R f (= SYNM \(N, C\))3 650 1 1713 2650 t (N)1440 2890 w 14 S f (\256)1980 2910 w 10 R f (the number of rows in A)5 979 1 2196 2890 t (C)1440 3130 w 14 S f (\256)1980 3150 w 10 R f ( lower tri-)2 417(a one-dimensional array of length N\(N+1\)/2 into which the)8 2427 2 2196 3130 t ( A is packed by columns as illustrated in the)9 1804(angular part of the matrix)4 1040 2 2196 3250 t (following 4)1 464 1 2196 3370 t 10 S f (\264)2668 3370 w 10 R f (4 example:)1 441 1 2731 3370 t 10 I f (a)2232 3610 w 7 R f (11)2293 3630 w 10 I f (c)3312 3610 w 7 R f (1)3367 3630 w 10 I f (a)2232 3730 w 7 R f (21)2293 3750 w 10 I f (a)2412 3730 w 7 R f (22)2473 3750 w 14 S f (\256)2952 3730 w 10 I f (c)3312 3730 w 7 R f (2)3367 3750 w 10 I f (c)3451 3730 w 7 R f (5)3506 3750 w 10 I f (a)2232 3850 w 7 R f (31)2293 3870 w 10 I f (a)2412 3850 w 7 R f (32)2473 3870 w 10 I f (a)2592 3850 w 7 R f (33)2653 3870 w 10 I f (c)3312 3850 w 7 R f (3)3367 3870 w 10 I f (c)3451 3850 w 7 R f (6)3506 3870 w 10 I f (c)3590 3850 w 7 R f (8)3645 3870 w 10 I f (a)2232 3970 w 7 R f (41)2293 3990 w 10 I f (a)2412 3970 w 7 R f (42)2473 3990 w 10 I f (a)2592 3970 w 7 R f (43)2653 3990 w 10 I f (a)2772 3970 w 7 R f (44)2833 3990 w 10 I f (c)3312 3970 w 7 R f (4)3367 3990 w 10 I f (c)3451 3970 w 7 R f (7)3506 3990 w 10 I f (c)3590 3970 w 7 R f (9)3645 3990 w 10 I f (c)3729 3970 w 7 R f (10)3784 3990 w 10 S f (<)1440 4260 w 10 R f (answer)1495 4260 w 10 S f (>)1777 4260 w 14 S f (\254)1980 4280 w 7 R f (1)2196 4330 w 7 S f (\243)2242 4330 w 7 I f (i)2292 4330 w 7 S f (\243)2317 4330 w 7 I f (n)2361 4330 w 10 R f (max)2210 4260 w 7 I f (j)2406 4360 w 7 S f (=)2437 4360 w 7 R f (1)2487 4360 w 15 S f (S)2420 4290 w 7 I f (n)2447 4160 w 10 S f (\357)2533 4260 w 10 I f (a)2590 4260 w 7 I f (i j)1 45 1 2651 4280 t 10 S f (\357)2736 4260 w 9 B f (Error situations:)1 648 1 720 4700 t 10 R f (\(All errors in this subprogram are fatal \320)7 1666 1 1512 4700 t (see)1512 4820 w 10 I f (Error Handling)1 631 1 1664 4820 t 10 R f (, Framework Chapter\))2 884 1 2295 4820 t (Number Error)1 1506 1 1800 5060 t (1 N)1 1108 1 1944 5300 t 10 S f (<)3077 5300 w 10 R f (1)3157 5300 w 9 B f (Double precision version:)2 981 1 720 5660 t 10 R f (DSYNM with C and DSYNM declared double precision)7 2261 1 1776 5660 t 9 B f (Complex version:)1 678 1 720 6020 t 10 R f (CSYNM with C declared complex)4 1382 1 1448 6020 t 9 B f (Storage:)720 6380 w 10 R f (None)1296 6380 w (Linear Algebra)1 606 1 4794 7500 t (\320 23 \320)2 350 1 2705 7620 t (SYNM)5111 7740 w cleartomark showpage saveobj restore %%EndPage: 23 23 %%Page: 24 24 /saveobj save def mark 24 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYNM February)1 4665 2 360 840 t 9 B f (Time:)720 1320 w 10 R f (N)1296 1320 w 7 R f (2)1373 1280 w 10 R f (additions)1441 1320 w (N comparisons)1 602 1 1296 1440 t 9 B f (See also:)1 333 1 720 1800 t 10 R f (SYDC, SYMD, SYLE, SYSS, SYCE)4 1502 1 1296 1800 t 9 B f (Author:)720 2160 w 10 R f (Linda Kaufman)1 629 1 1296 2160 t 9 B f (Example:)720 2520 w 10 R f ( for solving)2 467( library)1 355(The subroutines in the)3 900 3 1424 2520 t 10 I f (Ax)3174 2520 w 10 S f (=)3328 2520 w 10 I f (b)3432 2520 w 10 R f ( to return computed solu-)4 1026(are designed)1 504 2 3510 2520 t (tions)1296 2640 w 10 I f (x)1516 2640 w 10 R f (such that the residual)3 846 1 1585 2640 t 10 I f (r)2456 2640 w 10 S f (=)2544 2640 w 10 I f (Ax)2648 2640 w 10 S f (-)2802 2640 w 10 I f (b)2906 2640 w 10 R f (sati\256es)2981 2640 w 10 S f (\357 \357)1 106 1 2523 2990 t 10 I f (A)2637 2990 w 10 S f ( \357)1 57( \357)1 90(\357 \357)1 106 3 2706 2990 t 10 I f (x)2967 2990 w 10 S f (\357 \357)1 106 1 3019 2990 t (\357 \357)1 106 1 2691 2860 t 10 I f (r)2805 2860 w 10 S f (\357 \357)1 106 1 2852 2860 t 10 S1 f (_ ____________)1 632 1 2509 2890 t 10 S f (\243 e)1 107 1 3159 2920 t 10 R f (\(1.1\))4016 2920 w (where)1296 3210 w 10 S f (e)1571 3210 w 10 R f ( then)1 205( this example we show that if A is ill-condtioned,)9 2034( In)1 140(is the machine precision.)3 1014 4 1647 3210 t ( calose to the true solution even though equation \(1.1\) is)10 2325(the computed solution need not be)5 1419 2 1296 3330 t ( matrix in)2 393( The)1 206(satis\256ed. The subroutine SYNM is used to compute the left-hand side of \(1.1\).)12 3145 3 1296 3450 t (this example is given by)4 972 1 1296 3570 t 10 I f (a)2448 3810 w 7 I f (i j)1 45 1 2509 3830 t 10 S f (= \357)1 120 1 2578 3810 t 10 I f (i)2706 3810 w 10 S f (-)2758 3810 w 10 I f (j)2829 3810 w 10 S f (\357)2889 3810 w 10 R f (and the true solution is)4 911 1 1296 4050 t 10 I f (x)2232 4050 w 7 I f (i)2287 4070 w 10 S f (=)2331 4050 w 10 I f (i)2402 4050 w 10 R f ( the computed)2 568(. The right hand side is generated using SYML and)9 2042 2 2430 4050 t ( SYLE. The subroutine SAMAX is used to compute the 1-norm of)11 2703(solution is obtained using)3 1041 2 1296 4170 t (a vector.i.e.)1 465 1 1296 4290 t 7 R f (1)1811 4360 w 7 S f (\243)1851 4360 w 7 I f (i)1895 4360 w 7 S f (\243)1920 4360 w 7 I f (n)1964 4360 w 10 R f (max)1819 4290 w 10 S f (\357)2007 4290 w 10 I f (x)2064 4290 w 7 I f (i)2119 4310 w 10 S f (\357)2179 4290 w 8 CW f (INTEGER I, J, L, N, I1MACH, IWRITE)6 1632 1 2040 4680 t (REAL C\(1300\), CC\(1300\), B\(50\), X\(50\))4 1728 1 2040 4780 t (REAL RELERR, RELRES, XNORM, RNORM, ERR, R\(50\))6 2160 1 2040 4880 t (REAL SYNM, SAMAX)2 768 1 2040 4980 t (L=0)2040 5080 w (C)1656 5180 w (C GENERATE MATRIX)2 816 1 1656 5280 t (C)1656 5380 w (N=50)2040 5480 w (DO 20 I=1,N)2 528 1 2040 5580 t (DO 10 J=I,N)2 528 1 2184 5680 t (L=L+1)2328 5780 w (C\(L\)=J-I)2328 5880 w (CC\(L\)=C\(L\))2328 5980 w (10 CONTINUE)1 864 1 1704 6080 t (B\(I\)=I)2184 6180 w (20 CONTINUE)1 720 1 1704 6280 t (C)1656 6380 w (C GENERATE RIGHT HAND SIDE)4 1248 1 1656 6480 t (C)1656 6580 w (CALL SYML\(N,C,B,X\))1 864 1 2040 6680 t (C)1656 6780 w (C MAKE COPY OF RIGHT HAND SIDE)6 1440 1 1656 6880 t 10 R f (Linear Algebra)1 606 1 360 7540 t (\320 24 \320)2 350 1 2705 7660 t (SYNM)360 7780 w cleartomark showpage saveobj restore %%EndPage: 24 24 %%Page: 25 25 /saveobj save def mark 25 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYNM)1 4305(February 11, 1993)2 735 2 360 840 t 8 CW f (C)1656 1300 w (CALL MOVEFR\(N,X,B\))1 864 1 2040 1400 t (C)1656 1500 w (C SOLVE THE SYSTEM)3 864 1 1656 1600 t (C)1656 1700 w (CALL SYLE\(N,C,B,N,1\))1 960 1 2040 1800 t (C)1656 1900 w (C COMPUTE THE RELATIVE ERROR AND THE RELATIVE RESIDUAL)8 2592 1 1656 2000 t (C)1656 2100 w (CALL SYML\(N,CC,B,R\))1 912 1 2040 2200 t (ERR=0.0)2040 2300 w (DO 30 I=1,N)2 528 1 2040 2400 t (ERR=AMAX1\(ERR,ABS\(B\(I\)-FLOAT\(I\)\)\))2184 2500 w (R\(I\)=R\(I\)-X\(I\))2184 2600 w (30 CONTINUE)1 720 1 1704 2700 t (XNORM=SAMAX\(N,X,1\))2040 2800 w (RNORM=SAMAX\(N,R,1\))2040 2900 w (RELERR=ERR/XNORM)2040 3000 w (RELRES=RNORM/\(XNORM*SYNM\(N,CC\)\))2040 3100 w (IWRITE=I1MACH\(2\))2040 3200 w (WRITE\(IWRITE,31\)RELERR,RELRES)2040 3300 w ( RELATIVE ERROR=,E15.5,19H RELATIVE RESIDUAL=,)4 2208(31 FORMAT\(16H)1 816 2 1704 3400 t (1 E15.5\))1 480 1 1896 3500 t (STOP)2040 3600 w (END)2040 3700 w 10 R f ( on the Honeywell 6000 machine at Bell Laboratories,)8 2174(When the above program was executed)5 1570 2 1296 4040 t (the following was printed:)3 1052 1 1296 4160 t 8 CW f ( 0.95702E-11)1 624( RELATIVE RESIDUAL=)2 912( 0.10544E-07)1 624(RELATIVE ERROR=)1 720 4 1704 4500 t 10 R f ( is about 1300, and the ma-)6 1139(The condition number of the matrix \(see the example in SYSS\))10 2605 2 1296 4840 t ( is about 10)3 521(chine precision on the Honeywell computer)5 1844 2 1296 4960 t 7 S f (-)3672 4920 w 7 R f (8)3722 4920 w 10 R f (. Thus even in the absence of)6 1275 1 3765 4960 t ( of 1.3)2 264(roundoff error in SYML, a relative error)6 1621 2 1296 5080 t 10 S f (\264)3209 5080 w 10 R f (10)3292 5080 w 7 S f (-)3403 5040 w 7 R f (5)3453 5040 w 10 R f ( relative)1 327( The)1 208(would not be surprising.)3 981 3 3524 5080 t ( The relative residual, as promised, satis\256es \(1.1\))7 2049(error given above is quite within reason.)6 1695 2 1296 5200 t (even though the problem is slightly ill-conditioned.)6 2053 1 1296 5320 t (Linear Algebra)1 606 1 4794 7500 t (\320 25 \320)2 350 1 2705 7620 t (SYNM)5111 7740 w cleartomark showpage saveobj restore %%EndPage: 25 25 %%Page: 26 26 /saveobj save def mark 26 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYNM February)1 4665 2 360 840 t (SYSS \320 symmetric linear system solution with condition estimation)8 2763 1 1786 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( where A is a symmetric ma-)6 1177( = B)2 173( system AX)2 480(SYSS \(SYmmetric System Solution\) solves the)5 1914 4 1296 1680 t ( does not have to be posi-)6 1045( A)1 126( estimate of the condition number of A.)7 1604( also provides an)3 689(trix. It)1 280 5 1296 1800 t (tive de\256nite.)1 500 1 1296 1920 t 9 B f (Usage:)720 2280 w 10 R f (CALL SYSS \(N, C, B, IB, NB, COND\))7 1595 1 1296 2280 t (N)1440 2520 w 14 S f (\256)1980 2540 w 10 R f (the number of equations)3 968 1 2196 2520 t (C)1440 2760 w 14 S f (\256)1980 2780 w 10 R f ( lower tri-)2 417(a one-dimensional array of length N\(N+1\)/2 into which the)8 2427 2 2196 2760 t ( A is packed by columns as illustrated in the)9 1804(angular part of the matrix)4 1040 2 2196 2880 t (following 4)1 464 1 2196 3000 t 10 S f (\264)2660 3000 w 10 R f (4 example:)1 441 1 2715 3000 t 10 I f (a)2232 3240 w 7 R f (11)2293 3260 w 10 I f (c)3312 3240 w 7 R f (1)3367 3260 w 10 I f (a)2232 3360 w 7 R f (21)2293 3380 w 10 I f (a)2412 3360 w 7 R f (22)2473 3380 w 14 S f (\256)2952 3360 w 10 I f (c)3312 3360 w 7 R f (2)3367 3380 w 10 I f (c)3451 3360 w 7 R f (5)3506 3380 w 10 I f (a)2232 3480 w 7 R f (31)2293 3500 w 10 I f (a)2412 3480 w 7 R f (32)2473 3500 w 10 I f (a)2592 3480 w 7 R f (33)2653 3500 w 10 I f (c)3312 3480 w 7 R f (3)3367 3500 w 10 I f (c)3451 3480 w 7 R f (6)3506 3500 w 10 I f (c)3590 3480 w 7 R f (8)3645 3500 w 10 I f (a)2232 3600 w 7 R f (41)2293 3620 w 10 I f (a)2412 3600 w 7 R f (42)2473 3620 w 10 I f (a)2592 3600 w 7 R f (43)2653 3620 w 10 I f (a)2772 3600 w 7 R f (44)2833 3620 w 10 I f (c)3312 3600 w 7 R f (4)3367 3620 w 10 I f (c)3451 3600 w 7 R f (7)3506 3620 w 10 I f (c)3590 3600 w 7 R f (9)3645 3620 w 10 I f (c)3729 3600 w 7 R f (10)3784 3620 w 10 R f (C is overwritten during the solution.)5 1450 1 2196 3840 t (B)1440 4080 w 14 S f (\256)1980 4100 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 4080 t (the calling program, where IB)4 1200 1 2196 4200 t 10 S f (\263)3396 4200 w 10 R f (N and KB)2 405 1 3451 4200 t 10 S f (\263)3856 4200 w 10 R f (NB)3911 4200 w 14 S f (\254)1980 4460 w 10 R f (the solution X)2 567 1 2196 4440 t (IB)1440 4680 w 14 S f (\256)1980 4700 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 4680 t (calling program)1 635 1 2196 4800 t (NB)1440 5040 w 14 S f (\256)1980 5060 w 10 R f (the number of right-hand sides)4 1226 1 2196 5040 t (COND)1440 5280 w 14 S f (\254)1980 5300 w 10 R f ( \(see)1 210(an estimate of the condition number of A)7 1645 2 2196 5280 t 8 B f (Note 1)1 219 1 4076 5280 t 10 R f (\))4295 5280 w 9 B f (Note 1:)1 278 1 720 5640 t 10 R f ( to errors in)3 481(The condition number measures the sensitivity of the solution of a linear system)12 3263 2 1296 5640 t ( the elements of the matrix and the right-hand side\(s\))9 2134( If)1 118( and in the right-hand side.)5 1081(the matrix)1 411 4 1296 5760 t (of your linear system have)4 1091 1 1296 5880 t 10 B f (d)2420 5880 w 10 R f ( solution might have as few as)6 1264(decimal digits of precision, the)4 1267 2 2509 5880 t 10 B f (d)1296 6000 w 10 S f (-)1377 6000 w 10 R f (log)1457 6000 w 7 R f (10)1596 6020 w 10 R f ( if COND is greater than 10)6 1121( Thus)1 252(\(COND\) correct decimal digits.)3 1270 3 1674 6000 t 7 R f (Bd)4327 5960 w 7 I f (P)4419 5960 w 10 R f ( be)1 120(, there may)2 450 2 4470 6000 t (no correct digits.)2 674 1 1296 6120 t ( advance to be well-conditioned, then the user may wish to)10 2359(If the given matrix, A, is known in)7 1385 2 1296 6360 t ( however, the user is)4 867( Ordinarily,)1 502( a little faster than SYSS.)5 1066(use the routine SYLE, which is)5 1309 4 1296 6480 t (strongly urged to choose SYSS, and to follow it by a test of the condition estimate.)15 3313 1 1296 6600 t (Linear Algebra)1 606 1 360 7500 t (\320 26 \320)2 350 1 2705 7620 t (SYSS)360 7740 w cleartomark showpage saveobj restore %%EndPage: 26 26 %%Page: 27 27 /saveobj save def mark 27 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYSS)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Note 2:)1 278 1 720 1320 t 10 R f ( coef\256cient matrix, but differ-)4 1201(Users who wish to solve a sequence of problems with the same)11 2543 2 1296 1320 t (ent right-hand sides)2 798 1 1296 1440 t 10 I f ( known in advance,)3 783(not all)1 264 2 2124 1440 t 10 R f (should not use SYSS, but should call subpro-)7 1840 1 3200 1440 t ( is called once to get the)6 1061( SYCE)1 322( the example of SYCE.\))4 1021( \(See)1 243( SYFBS.)1 373(grams SYCE and)2 724 6 1296 1560 t (MDM)1296 1680 w 7 R f (T)1551 1640 w 10 R f ( to this chapter\) and then SYFBS is called for)9 1926(decomposition \(see the introduction)3 1474 2 1640 1680 t (each new right-hand side.)3 1025 1 1296 1800 t 9 B f (Error situations:)1 653 1 720 2160 t 10 R f ( see)1 183( \320)1 125( those errors marked with an asterisk)6 1505(*\(The user can elect to `recover' from)6 1546 4 1517 2160 t 10 I f (Er-)4907 2160 w (ror Handling)1 531 1 1512 2280 t 10 R f (, Framework Chapter\))2 884 1 2043 2280 t (Number Error)1 1506 1 1800 2520 t (1 N)1 1008 1 1944 2760 t 10 S f (<)2977 2760 w 10 R f (1)3057 2760 w (2 IB)1 1036 1 1944 3000 t 10 S f (<)3005 3000 w 10 R f (N)3085 3000 w (3 NB)1 1075 1 1944 3240 t 10 S f (<)3044 3240 w 10 R f (1)3124 3240 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 3480 t 9 B f (Double-precision version:)1 988 1 720 3840 t 10 R f (DSYSS with C, B, and COND declared double precision.)8 2301 1 1800 3840 t 9 B f (Complex symmetric version:)2 1106 1 720 4200 t 10 R f (CSYSS with C and B declared complex)6 1594 1 1901 4200 t 9 B f (Complex Hermitian version:)2 1101 1 720 4560 t 10 R f (CHESS with C and B declared complex)6 1599 1 1896 4560 t 9 B f (Storage:)720 4920 w 10 R f (N integer locations and)3 929 1 1296 4920 t ( scratch)1 315(N real \(double precision for DSYSS, complex for CSYSS and CHESS\) locations of)12 3429 2 1296 5040 t (storage in the dynamic storage stack)5 1450 1 1296 5160 t 9 B f (Time:)720 5570 w 10 R f (6)1356 5640 w (N)1321 5510 w 7 R f (3)1398 5470 w 10 S1 f (___)1306 5540 w 10 S f (+)1515 5570 w 10 R f (\()1619 5570 w (4)1710 5640 w (19)1685 5510 w 10 S1 f (_ __)1 130 1 1670 5540 t 10 S f (+)1826 5570 w 10 R f (NB \))1 180 1 1897 5570 t 10 S f (\264)2126 5570 w 10 R f (N)2222 5570 w 7 R f (2)2299 5530 w 10 S f (+)2391 5570 w 10 R f (\()2495 5570 w (6)2586 5640 w (25)2561 5510 w 10 S1 f (_ __)1 130 1 2546 5540 t 10 S f (+)2702 5570 w 10 R f (NB \))1 180 1 2773 5570 t 10 S f (\264)3002 5570 w 10 R f (N additions)1 464 1 3098 5570 t (6)1356 5920 w (N)1321 5790 w 7 R f (3)1398 5750 w 10 S1 f (___)1306 5820 w 10 S f (+)1515 5850 w 10 R f (\()1619 5850 w (4)1710 5920 w (13)1685 5790 w 10 S1 f (_ __)1 130 1 1670 5820 t 10 S f (+)1826 5850 w 10 R f (NB \))1 180 1 1897 5850 t 10 S f (\264)2126 5850 w 10 R f (N)2222 5850 w 7 R f (2)2299 5810 w 10 S f (+)2391 5850 w 10 R f (\()2495 5850 w (6)2561 5920 w (5)2561 5790 w 10 S1 f (_ _)1 80 1 2546 5820 t 10 S f (+)2652 5850 w 10 R f (NB \))1 180 1 2723 5850 t 10 S f (\264)2952 5850 w 10 R f (N multiplications)1 698 1 3048 5850 t (2)1356 6200 w (N)1321 6070 w 7 R f (2)1398 6030 w 10 S1 f (___)1306 6100 w 10 S f (+)1515 6130 w 10 R f (\()1619 6130 w (2)1685 6200 w (5)1685 6070 w 10 S1 f (_ _)1 80 1 1670 6100 t 10 S f (+)1776 6130 w 10 R f (NB \))1 180 1 1847 6130 t 10 S f (\264)2076 6130 w 10 R f (N divisions)1 459 1 2172 6130 t (at most N)2 389 1 1296 6420 t 7 R f (2)1690 6380 w 10 S f (+)1749 6420 w 10 R f (N comparisons)1 602 1 1820 6420 t (Linear Algebra)1 606 1 4794 7500 t (\320 27 \320)2 350 1 2705 7620 t (SYSS)5160 7740 w cleartomark showpage saveobj restore %%EndPage: 27 27 %%Page: 28 28 /saveobj save def mark 28 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYSS February)1 4665 2 360 840 t 9 B f (Method:)720 1320 w 10 R f ( reference [2])2 541( See)1 197(The Bunch - Kaufman algorithm described in reference [1] below, is used.)11 3006 3 1296 1320 t (below for the method used to estimate the condition number.)9 2431 1 1296 1440 t (SYSS calls SYCE and SYFBS.)4 1255 1 1296 1560 t 9 B f (See also:)1 333 1 720 1920 t 10 R f (SYDC, SYFBS, SYMD, SYLE, SYCE)4 1569 1 1296 1920 t 9 B f (Author:)720 2280 w 10 R f (Linda Kaufman)1 629 1 1296 2280 t 9 B f (References:)720 2640 w 10 R f ( Kaufman, L., and Parlett, B., Decomposition of a symmetric matrix,)10 2790( J. R.,)2 237([1] Bunch,)1 538 3 1296 2640 t 10 I f (Nu-)4890 2640 w (mer. Math 27)2 541 1 1548 2760 t 10 R f (\(1976\), 95-109.)1 624 1 2114 2760 t ( H., An estimate for the)5 980( A. K., Moler, C. B., Stewart, G. W., and Wilkinson, J.)11 2270([2] Cline,)1 494 3 1296 3000 t (condition number,)1 733 1 1548 3120 t 10 I f (SIAM J. Numer. Anal. 16)4 1007 1 2331 3120 t 10 R f (\(1979\), 368-375.)1 674 1 3363 3120 t 9 B f (Example:)720 3480 w 10 R f (The following program packs the matrix)5 1617 1 1296 3480 t 10 I f (a)2592 3720 w 7 I f (i)2653 3740 w 7 R f (,)2678 3740 w 7 I f (j)2707 3740 w 10 S f (= \357)1 120 1 2751 3720 t 10 I f (i)2879 3720 w 10 S f (-)2931 3720 w 10 I f (j)3002 3720 w 10 S f (\357)3038 3720 w 10 R f ( right-hand side so that the solution will be all ones. It then)12 2424(into the vector C and sets up the)7 1320 2 1296 3960 t ( solution. This example was)4 1151(calls SYSS to solve the problem and calculates the error in the)11 2593 2 1296 4080 t ( looking problems may be ill-conditioned and to)7 2021(included to show that seemingly innocent)5 1723 2 1296 4200 t ( effect of ill-conditioning on a solution. It also demonstrates how to pack a symmet-)14 3383(show the)1 361 2 1296 4320 t (ric matrix into a vector.)4 940 1 1296 4440 t 8 CW f (INTEGER N, L, I, IWRITE, I1MACH)5 1488 1 2088 4760 t (REAL C\(5000\), B\(100\))2 960 1 2088 4860 t (REAL SUM, FLOAT, ABS, ERR, COND)5 1488 1 2088 4960 t ( N=10,90,40)1 576(DO 40)1 240 2 2088 5060 t (C)1656 5160 w (C CREATE THE MATRIX A\(I,J\)=ABS\(I)4 1536 1 1656 5260 t 8 S f (-)3192 5260 w 8 CW f (J\), PACK IT INTO)3 768 1 3236 5260 t (C THE VECTOR C AND FORM THE RIGHT-HAND SIDE SO THE)10 2400 1 1656 5360 t (C SOLUTION HAS ALL ONES.)4 1152 1 1656 5460 t (C)1656 5560 w (L=1)2184 5660 w (SUM=\(N*\(N-1\)\)/2)2184 5760 w (DO 20 I=1,N)2 528 1 2184 5860 t (DO 10 J=I,N)2 528 1 2328 5960 t (C\(L\)=J-I)2472 6060 w (L=L+1)2472 6160 w (10 CONTINUE)1 960 1 1752 6260 t (B\(I\)=SUM)2328 6360 w (SUM=SUM+FLOAT\(I-\(N-I\)\))2328 6460 w (20 CONTINUE)1 816 1 1752 6560 t (C)1656 6660 w (C SOLVE THE SYSTEM AND GET THE CONDITION NUMBER OF THE MATRIX)11 2928 1 1656 6760 t (CALL SYSS\(N,C,B,100,1,COND\))1 1296 1 2184 6860 t 10 R f (Linear Algebra)1 606 1 360 7520 t (\320 28 \320)2 350 1 2705 7640 t (SYSS)360 7760 w cleartomark showpage saveobj restore %%EndPage: 28 28 %%Page: 29 29 /saveobj save def mark 29 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( SYSS)1 4305(February 11, 1993)2 735 2 360 840 t 8 CW f (C)1656 1300 w (C COMPUTE THE ERROR IN THE SOLUTION)6 1680 1 1656 1400 t (ERR=0.0)2184 1500 w (DO 30 I=1,N)2 528 1 2184 1600 t (30 ERR=ERR+ABS\(B\(I\)-1.0\))1 1632 1 1752 1700 t (ERR=ERR/FLOAT\(N\))2184 1800 w (IWRITE=I1MACH\(2\))2184 1900 w (WRITE\(IWRITE,31\)N)2184 2000 w ( FOR N= ,I5\))3 576(31 FORMAT\(/8H)1 912 2 1752 2100 t (WRITE\(IWRITE,32\)COND)2184 2200 w ( CONDITION ESTIMATE IS 1PE15.7\))4 1488(32 FORMAT\(23H)1 912 2 1752 2300 t (WRITE\(IWRITE,33\)ERR)2184 2400 w ( RELATIVE ERROR IN SOLUTION IS,1PE15.7\))5 1872(33 FORMAT\(30H)1 912 2 1752 2500 t (40 CONTINUE)1 720 1 1752 2600 t (STOP)2088 2700 w (END)2088 2800 w 10 R f (The output from the above program run on the Honeywell 6000 at Bell Labs was)14 3234 1 1296 3140 t 8 CW f ( 10)1 288(FOR N=)1 288 2 1704 3560 t ( 01)1 144( 5.4831256E)1 624(CONDITION ESTIMATE IS)2 1008 3 1704 3660 t ( 7.3015689E-08)1 720(RELATIVE ERROR IN SOLUTION IS)4 1392 2 1704 3760 t ( 50)1 288(FOR N=)1 288 2 1704 3960 t ( 03)1 144( 1.3551582E)1 624(CONDITION ESTIMATE IS)2 1008 3 1704 4060 t ( 4.9673021E-06)1 720(RELATIVE ERROR IN SOLUTION IS)4 1392 2 1704 4160 t ( 90)1 288(FOR N=)1 288 2 1704 4360 t ( 03)1 144( 4.4305656E)1 624(CONDITION ESTIMATE IS)2 1008 3 1704 4460 t ( 1.2567225E-05)1 720(RELATIVE ERROR IN SOLUTION IS)4 1392 2 1704 4560 t 10 R f ( solution)1 350(As the output indicates, for small values of N, the matrix is well-conditioned and the)14 3394 2 1296 4920 t ( as N increases, the matrix becomes more ill-conditioned and the error in)12 2932(is very accurate, but)3 812 2 1296 5040 t ( is not unrea-)3 548( the relative error for N = 90 appears large, it)10 1862( Although)1 435(the solution increases.)2 899 4 1296 5160 t (sonable as the following analysis indicates:)5 1729 1 1296 5280 t (Linear Algebra)1 606 1 4794 7500 t (\320 29 \320)2 350 1 2705 7620 t (SYSS)5160 7740 w cleartomark showpage saveobj restore %%EndPage: 29 29 %%Page: 30 30 /saveobj save def mark 30 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(SYSS February)1 4665 2 360 840 t (Let)1296 1320 w 10 S f (D)1479 1320 w 10 I f (b)1548 1320 w 10 R f (represent a perturbation in the right-hand side of a linear system.)10 2581 1 1623 1320 t (If)1296 1440 w 10 I f (Ax)1387 1440 w 10 S f (=)1541 1440 w 10 I f (b)1645 1440 w 10 R f (then)1720 1440 w 10 I f (A)2448 1680 w 10 R f (\()2517 1680 w 10 I f (x)2558 1680 w 10 S f (+ D)1 132 1 2626 1680 t 10 I f (x)2766 1680 w 10 R f (\))2818 1680 w 10 S f (=)2908 1680 w 10 I f (b)3012 1680 w 10 S f (+ D)1 132 1 3086 1680 t 10 I f (b)3226 1680 w 10 R f (where)1296 1920 w 10 S f (\357 \357)1 106 1 2363 2315 t 10 I f (x)2477 2315 w 10 S f (\357 \357)1 106 1 2529 2315 t (\357 \357 D)2 175 1 2329 2185 t 10 I f (x)2512 2185 w 10 S f (\357 \357)1 106 1 2564 2185 t 10 S1 f (_ _______)1 371 1 2314 2215 t 10 S f (\243)2736 2245 w 10 I f (K)2832 2245 w 10 R f (\()2907 2245 w 10 I f (A)2948 2245 w 10 R f (\))3017 2245 w 10 S f (\354)3074 2158 w (\357)3074 2258 w (\356)3074 2358 w (\357 \357)1 106 1 3182 2315 t 10 I f (b)3296 2315 w 10 S f (\357 \357)1 106 1 3354 2315 t (\357 \357 D)2 175 1 3148 2185 t 10 I f (b)3331 2185 w 10 S f (\357 \357)1 106 1 3389 2185 t 10 S1 f (_ _______)1 377 1 3133 2215 t 10 S f (\374)3520 2158 w (\357)3520 2258 w (\376)3520 2358 w 10 R f (\(1.1\))4044 2245 w (where)1296 2590 w 10 I f (K)1573 2590 w 10 R f (\()1648 2590 w 10 I f (A)1689 2590 w 10 R f ( condition number of)3 871(\) is the)2 290 2 1758 2590 t 10 I f (A)2954 2590 w 10 R f (,)3015 2590 w 10 I f (K)3075 2590 w 10 R f (\()3150 2590 w 10 I f (A)3191 2590 w 10 R f (\))3260 2590 w 10 S f ( \357)1 57(= \357)1 153 2 3350 2590 t 10 I f (A)3568 2590 w 10 S f ( \357)1 57( \357)1 90(\357 \357)1 106 3 3637 2590 t 10 I f (A)3898 2590 w 7 S f (-)3970 2550 w 7 R f (1)4020 2550 w 10 S f (\357 \357)1 106 1 4095 2590 t 10 R f (and)4236 2590 w 10 S f (\357 \357)1 106 1 4415 2590 t 10 R f (.)4529 2560 w 10 S f (\357 \357)1 106 1 4586 2590 t 10 R f (is some)1 313 1 4727 2590 t (norm e.g.,)1 405 1 1296 2760 t 10 S f (\357 \357)1 106 1 1726 2760 t 10 I f (x)1840 2760 w 10 S f (\357 \357)1 106 1 1892 2760 t 7 R f (1)2009 2780 w 10 S f (=)2101 2760 w 7 I f (i)2205 2860 w 7 S f (=)2241 2860 w 7 R f (1)2291 2860 w 15 S f (S)2221 2790 w 7 I f (n)2248 2660 w 10 S f (|)2334 2760 w 10 I f (x)2362 2760 w 7 I f (i)2417 2780 w 10 S f (|)2453 2760 w 10 R f (if)2498 2760 w 10 I f (x)2584 2760 w 10 R f (is a vector.)2 435 1 2653 2760 t ( provide an accurate an-)4 993(The methods used in our linear equation package are guaranteed to)10 2751 2 1296 3080 t ( we assume that our method produces the correct an-)9 2130( If)1 119( problem.)1 386(swer to a slightly perturbed)4 1109 4 1296 3200 t (swer to a problem where)4 1006 1 1296 3320 t 10 S f (\357 \357 D)2 175 1 2332 3320 t 10 I f (b)2515 3320 w 10 S f ( \357 \357)2 114( \243 e)2 181(\357 \357)1 106 3 2573 3320 t 10 I f (b)2982 3320 w 10 S f (\357 \357)1 106 1 3040 3320 t 10 R f (, where)1 298 1 3146 3320 t 10 S f (e)3474 3320 w 10 R f (is the machine precision, then on the)6 1492 1 3548 3320 t (Honeywell 6000 where)2 955 1 1296 3440 t 10 S f (e)2288 3440 w 10 R f (is about 10)2 463 1 2369 3440 t 7 S f (-)2843 3400 w 7 R f (8)2893 3400 w 10 R f ( error of 4)3 434(, a relative)2 442 2 2936 3440 t 10 S f (\264)3820 3440 w 10 R f (10)3883 3440 w 7 S f (-)3994 3400 w 7 R f (5)4044 3400 w 10 R f (for the above problem)3 917 1 4123 3440 t (with)1296 3560 w 10 I f (N)1499 3560 w 10 R f (of 90 would not be surprising.)5 1205 1 1591 3560 t (Linear Algebra)1 606 1 360 7500 t (\320 30 \320)2 350 1 2705 7620 t (SYSS)360 7740 w (-- --)1 5544 1 0 8030 t cleartomark showpage saveobj restore %%EndPage: 30 30 %%Trailer done %%Pages: 30 %%DocumentFonts: Courier Times-Bold Times-Italic Times-Roman Times-Roman Symbol