%!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 4)1 469 1 2825 120 t (BANDED)1746 360 w 10 S f (<)2162 360 w 10 R f (SYMMETRIC)2242 360 w 10 S f (<)2837 360 w 10 R f (POSITIVE-DEFINITE MATRICES)1 1457 1 2917 360 t ( Solve)1 253( Back)1 380(BPBS -)1 379 3 720 720 t ( Estimation)1 459( Condition)1 576(BPCE -)1 384 3 720 840 t ( DeComposition)1 809(BPDC -)1 395 2 720 960 t ( Solve)1 253( Forward)1 513(BPFS -)1 368 3 720 1080 t ( Equation solution)2 734( Linear)1 435(BPLE -)1 378 3 720 1200 t (BPLD -)1 389 1 720 1320 t 10 I f (LDL)1284 1320 w 7 I f (T)1479 1280 w 10 R f (decomposition)1551 1320 w ( MuLtiplication)1 781(BPML -)1 406 2 720 1440 t ( NorM)1 419(BPNM -)1 417 2 720 1560 t ( Solution)1 365( System)1 470(BPSS -)1 368 3 720 1680 t cleartomark showpage saveobj restore %%EndPage: 1 1 %%Page: 2 2 /saveobj save def mark 2 pagesetup 10 R f (BPBS \320 band positive de\256nite upper triangular linear system solution)9 2825 1 1467 120 t 9 B f (Purpose:)720 480 w 10 R f ( where D is)3 475( = B)2 173(BPBS \(Banded symmetric Positive de\256nite matrix Back Solve\) solves DRX)9 3096 3 1296 480 t ( \(1's on the diagonal and 0's below the diago-)9 1860(diagonal and R is banded unit upper triangular)7 1884 2 1296 600 t ( is used)2 300( \(It)1 144( can be used for the back solution phase of a banded linear system solution.)14 3009(nal\). It)1 291 4 1296 720 t (in this way by the routines BPSS and BPLE.\))8 1815 1 1296 840 t 9 B f (Usage:)720 1200 w 10 R f (CALL BPBS \(N, ML, G, IG, B, IB, NB\))8 1628 1 1296 1200 t (N)1440 1440 w 14 S f (\256)1980 1460 w 10 R f (the number of equations)3 968 1 2196 1440 t (ML)1440 1680 w 14 S f (\256)1980 1700 w 10 R f (the number of nonzero diagonals of R \(including the)8 2097 1 2196 1680 t (unit diagonal\))1 558 1 2196 1800 t (G)1440 2040 w 14 S f (\256)1980 2060 w 10 R f (a matrix \(which may contain results obtained by the)8 2075 1 2196 2040 t ( or BPDC\) into which D and R are packed as)10 1897(routines BPLD, BPCE,)2 947 2 2196 2160 t (follows:)2196 2280 w (G\(1, i\) =)2 347 1 2916 2460 t 10 I f (d)3288 2460 w 7 I f (i)3349 2480 w 10 R f (G\(j)2916 2580 w 10 S f (-)3049 2580 w 10 R f (i + 1, i\) =)4 376 1 3104 2580 t 10 I f (r)3505 2580 w 7 I f (i j)1 45 1 3555 2600 t 10 R f (for j)1 169 1 3633 2580 t 10 S f (>)3827 2580 w 10 R f (i)3907 2580 w (\(See the introduction to this chapter.\))5 1487 1 2196 2760 t ( in the calling program, where)5 1348(G should be dimensioned \(IG,KG\))4 1496 2 2196 2880 t (IG)2196 3000 w 10 S f (\263)2301 3000 w 10 R f (ML and KG)2 488 1 2356 3000 t 10 S f (\263)2844 3000 w 10 R f (N.)2899 3000 w (IG)1440 3240 w 14 S f (\256)1980 3260 w 10 R f (the row \(leading\) dimension of G, as dimensioned in the)9 2253 1 2196 3240 t (calling program)1 635 1 2196 3360 t (B)1440 3600 w 14 S f (\256)1980 3620 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 3600 t (the calling program, where IB)4 1200 1 2196 3720 t 10 S f (\263)3396 3720 w 10 R f (N and KB)2 405 1 3451 3720 t 10 S f (\263)3856 3720 w 10 R f (NB.)3911 3720 w 14 S f (\254)1980 3980 w 10 R f (the solution X)2 567 1 2196 3960 t (IB)1440 4200 w 14 S f (\256)1980 4220 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 4200 t (calling program)1 635 1 2196 4320 t (NB)1440 4560 w 14 S f (\256)1980 4580 w 10 R f (the number of right-hand sides)4 1226 1 2196 4560 t 9 B f (Note:)720 4920 w 10 R f ( can be used directly on the output matrix produced by BPDC, BPLD, or)13 3049(BPFS and BPBS)2 695 2 1296 4920 t (BPCE to solve a banded symmetric positive de\256nite linear system.)9 2666 1 1296 5040 t (Linear Algebra)1 606 1 360 7500 t (\320 2 \320)2 300 1 2730 7620 t (BPBS)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 ( BPBS)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 ML)1 1086 1 1944 2160 t 10 S f (<)3055 2160 w 10 R f (1)3135 2160 w (3 IG)1 1041 1 1944 2400 t 10 S f (<)3010 2400 w 10 R f (ML)3090 2400 w (4 IB)1 1036 1 1944 2640 t 10 S f (<)3005 2640 w 10 R f (N)3085 2640 w (5 NB)1 1075 1 1944 2880 t 10 S f (<)3044 2880 w 10 R f (1)3124 2880 w ( D with kth diagonal element 0.0)6 1313( singular)1 952(10 + k*)2 306 3 1944 3120 t 9 B f (Double-precision version:)1 988 1 720 3480 t 10 R f (DBPBS, with G and B declared double precision)7 1954 1 1758 3480 t 9 B f (Complex Hermitian version:)2 1101 1 720 3840 t 10 R f (CBPBS with G and B declared complex)6 1605 1 1896 3840 t 9 B f (Storage:)720 4200 w 10 R f (None)1296 4200 w 9 B f (Time:)720 4560 w 10 R f (NB)1296 4560 w 10 S f (\264)1435 4560 w 10 R f (\(ML)1490 4560 w 10 S f (-)1673 4560 w 10 R f (1\))1728 4560 w 10 S f (\264)1811 4560 w 10 R f (N additions)1 464 1 1866 4560 t (NB)1296 4680 w 10 S f (\264)1435 4680 w 10 R f (\(ML)1490 4680 w 10 S f (-)1673 4680 w 10 R f (1\))1728 4680 w 10 S f (\264)1811 4680 w 10 R f (N multiplications)1 698 1 1866 4680 t (NB)1296 4800 w 10 S f (\264)1435 4800 w 10 R f (N divisions)1 459 1 1490 4800 t 9 B f (See also:)1 333 1 720 5160 t 10 R f (BPCE, BPDC, BPFS, BPLD, BPLE, BPSS)5 1734 1 1296 5160 t 9 B f (Author:)720 5520 w 10 R f (Linda Kaufman)1 629 1 1296 5520 t (Linear Algebra)1 606 1 4794 7500 t (\320 3 \320)2 300 1 2730 7620 t (BPBS)5154 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(BPBS February)1 4665 2 360 840 t 9 B f (Example:)720 1320 w 10 R f ( a linear system AX = B, where A is a symmetric posi-)12 2273(The program fragment below solves)4 1471 2 1296 1320 t ( A matrix has been packed into G according)8 1797( is assumed that the)4 807( It)1 117(tive de\256nite band matrix.)3 1023 4 1296 1440 t (to the scheme G\(j)3 734 1 1296 1560 t 10 S f (-)2030 1560 w 10 R f (i+1, i\)=)1 310 1 2085 1560 t 10 I f (a)2395 1560 w 7 I f (i j)1 45 1 2456 1580 t 10 R f (. The subroutine BPCE factors A into LDL)7 1784 1 2509 1560 t 7 R f (T)4298 1520 w 10 R f ( L is unit)3 389(, where)1 302 2 4349 1560 t ( factors are returned in G so that BPFS can forward)10 2111( The)1 211(lower triangular and D is diagonal.)5 1422 3 1296 1680 t (solve \(solve LY=B\) and BPBS can back solve \(solve DL)9 2273 1 1296 1800 t 7 R f (T)3574 1760 w 10 R f (X=Y\).)3625 1800 w ( also provides an estimate of the condition number of A. In the code)13 2848(The subroutine BPCE)2 896 2 1296 2040 t ( is larger than the reciprocal of the machine precision \(given by)11 2530(below if the condition number)4 1214 2 1296 2160 t (R1MACH\(4\)\), the matrix is considered too ill-conditioned and the system is not solved.)12 3515 1 1296 2280 t 8 CW f (IWRITE=I1MACH\(2\))2040 2620 w (CALL BPCE\(N,ML,G,IG,COND\))1 1200 1 2040 2720 t (IF \(COND .GT. 1.0/R1MACH\(4\)\) GO TO 10)6 1776 1 2040 2820 t (CALL BPFS\(N,ML,G,IG,B,IB,NB\))1 1344 1 2040 2920 t (CALL BPBS\(N,ML,G,IG,B,IB,NB\))1 1344 1 2040 3020 t (GO TO 20)2 384 1 2040 3120 t (10 WRITE\(IWRITE,11\))1 1152 1 1656 3220 t ( MATRIX TOO ILL-CONDITIONED\))3 1344(11 FORMAT\(26H)1 864 2 1656 3320 t (20 CONTINUE)1 768 1 1656 3420 t 10 R f (Linear Algebra)1 606 1 360 7500 t (\320 4 \320)2 300 1 2730 7620 t (BPBS)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 ( BPBS)1 4305(February 11, 1993)2 735 2 360 840 t (BPCE \320 LDL)2 595 1 2011 1320 t 7 R f (T)2611 1280 w 10 R f (decomposition with condition estimation)3 1637 1 2687 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( symmetric Positive de\256nite matrix Condition Estimation\) gives a lower)9 3106(BPCE \(Banded)1 638 2 1296 1680 t ( re-)1 136( also)1 188( It)1 113(bound for the condition number of a banded symmetric positive de\256nite matrix A.)12 3307 4 1296 1800 t (turns the LDL)2 566 1 1296 1920 t 7 R f (T)1867 1880 w 10 R f (decomposition of A and may be used in a linear equation package.)11 2656 1 1943 1920 t 9 B f (Usage:)720 2280 w 10 R f (CALL BPCE \(N, MU, G, IG, COND\))6 1521 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 (MU)1440 2760 w 14 S f (\256)1980 2780 w 10 R f (the number of nonzero bands on and above the diagonal of A)11 2442 1 2196 2760 t (G)1440 3000 w 14 S f (\256)1980 3020 w 10 R f ( which the upper triangular portion of the matrix A has)10 2311(a matrix into)2 533 2 2196 3000 t (been packed as follows:)3 956 1 2196 3120 t (G\(j)2916 3300 w 10 S f (-)3049 3300 w 10 R f (i + 1, i\) =)4 376 1 3104 3300 t 10 I f (a)3505 3300 w 7 I f (i j)1 45 1 3566 3320 t 10 R f (for j)1 169 1 3644 3300 t 10 S f (\263)3813 3300 w 10 R f (i)3868 3300 w (\(See the introduction to this chapter.\))5 1487 1 2196 3480 t ( in the calling program, where)5 1348(G should be dimensioned \(IG,KG\))4 1496 2 2196 3600 t (IG)2196 3720 w 10 S f (\263)2301 3720 w 10 R f (MU and KG)2 499 1 2356 3720 t 10 S f (\263)2855 3720 w 10 R f (N.)2910 3720 w 14 S f (\254)1980 3980 w 10 R f (LDL)2196 3960 w 7 R f (T)2395 3920 w 10 R f ( \(see)1 188(decomposition suitable for input into BPFS and BPBS)7 2192 2 2473 3960 t 8 B f (Note)4881 3960 w (2)2196 4080 w 10 R f (\))2236 4080 w (IG)1440 4320 w 14 S f (\256)1980 4340 w 10 R f (the row \(leading\) dimension of G, as dimensioned in the)9 2253 1 2196 4320 t (calling program)1 635 1 2196 4440 t (COND)1440 4680 w 14 S f (\254)1980 4700 w 10 R f (an estimate of the condition number of A \(see)8 1830 1 2196 4680 t 8 B f (Note 1)1 219 1 4051 4680 t 10 R f (\))4270 4680 w 9 B f (Note 1:)1 278 1 720 5040 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 5040 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 5160 t (of your linear system have)4 1091 1 1296 5280 t 10 B f (d)2420 5280 w 10 R f ( solution might have as few as)6 1264(decimal digits of precision, the)4 1267 2 2509 5280 t 10 B f (d)1296 5400 w 10 S f (-)1377 5400 w 10 R f (log)1457 5400 w 7 R f (10)1596 5420 w 10 R f ( if COND is greater than 10)6 1109( Thus)1 250(\(COND\) correct decimal digits.)3 1264 3 1699 5400 t 7 R f (Bd)4332 5360 w 7 I f (P)4424 5360 w 10 R f (, there may be)3 565 1 4475 5400 t (no correct digits.)2 674 1 1296 5520 t 9 B f (Note 2:)1 278 1 720 5880 t 10 R f (The LDL)1 378 1 1296 5880 t 7 R f (T)1679 5840 w 10 R f (decomposition of A satis\256es the equation A = LDL)8 2081 1 1759 5880 t 7 R f (T)3845 5840 w 10 R f (where L is lower unit trian-)5 1115 1 3925 5880 t ( 0's above the diagonal\) and D is diagonal. On return from BPCE,)12 2677(gular \(1's on the diagonal,)4 1067 2 1296 6000 t (the diagonal of D occupies the \256rst row of G and)10 1952 1 1296 6120 t (G\(i)1296 6240 w 10 S f (-)1454 6240 w 10 R f (j+1,i\))1534 6240 w 10 S f (=)1795 6240 w 10 I f (l)1899 6240 w 7 I f (i j)1 45 1 1938 6260 t 10 R f (for i)1 169 1 2016 6240 t 10 S f (>)2185 6240 w 10 R f (j.)2240 6240 w 9 B f (Note 3:)1 278 1 720 6600 t 10 R f ( represents the conjugate transpose of A\))6 1629( *)1 58( where A)2 365( *,)1 83(For complex Hermitian matrices \(A = A)6 1609 5 1296 6600 t ( decomposition and returns the)4 1260( *)1 58( this subroutine computes the LDL)5 1426(, the complex version of)4 1000 4 1296 6720 t (conjugate of L rather than L in G.)7 1347 1 1296 6840 t (Linear Algebra)1 606 1 4794 7500 t (\320 5 \320)2 300 1 2730 7620 t (BPCE)5149 7740 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(BPCE 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 (2 MU)1 1097 1 1944 2160 t 10 S f (<)3066 2160 w 10 R f (1)3146 2160 w (3 IG)1 1041 1 1944 2400 t 10 S f (<)3010 2400 w 10 R f (MU)3090 2400 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2640 t ( k)1 75( the)1 574(10 + N + k*)4 484 3 1944 2880 t 7 R f (th)3082 2840 w 10 R f (principal minor is not positive de\256nite)5 1531 1 3170 2880 t 9 B f (Double-precision version:)1 988 1 720 3240 t 10 R f (DBPCE with G and COND declared double precision)7 2150 1 1783 3240 t 9 B f (Complex Hermitian version:)2 1101 1 720 3600 t 10 R f (CBPCE with G declared complex \(see)5 1534 1 1896 3600 t 8 B f (Note 3)1 219 1 3455 3600 t 10 R f (above\).)3699 3600 w 9 B f (Storage:)720 3960 w 10 R f (N real \(double precision for DBPCE, complex for CBPCE\) locations of scratch)11 3168 1 1296 3960 t (storage in the dynamic storage stack)5 1450 1 1296 4080 t 9 B f (Time:)720 4440 w 10 R f (N)1296 4440 w 10 S f (\264)1393 4440 w 10 R f (\(\(MU)1473 4440 w 10 S f (-)1700 4440 w 10 R f (1\))1755 4440 w 10 S f (\264)1863 4440 w 10 R f (\(10+MU/2\)+8\) additions)1 992 1 1943 4440 t (N)1296 4560 w 10 S f (\264)1393 4560 w 10 R f (\(\(MU)1473 4560 w 10 S f (-)1700 4560 w 10 R f (1\))1755 4560 w 10 S f (\264)1863 4560 w 10 R f (\(6+MU/2\)+4\) multiplications)1 1176 1 1943 4560 t (N)1296 4680 w 10 S f (\264)1393 4680 w 10 R f (\(MU+1\) divisions)1 720 1 1473 4680 t 9 B f (Method:)720 5040 w 10 R f (Gaussian elimination without pivoting)3 1537 1 1296 5040 t (See the reference below for the method used to estimate the condition number.)12 3141 1 1296 5160 t (BPCE calls BPLD with EPS=0.0)4 1322 1 1296 5280 t 9 B f (See also:)1 333 1 720 5640 t 10 R f (BPBS, BPDC, BPFS, BPLD, BPLE, BPSS)5 1729 1 1296 5640 t 9 B f (Authors:)720 6000 w 10 R f (Doris Ryan and Linda Kaufman)4 1281 1 1296 6000 t 9 B f (Reference:)720 6360 w 10 R f ( W., and Wilkinson, J. H., An estimate for the condi-)10 2154(Cline, A. K., Moler, C. B., Stewart, G.)7 1562 2 1324 6360 t (tion number,)1 511 1 1296 6480 t 10 I f (SIAM J. Numer. Anal. 16)4 1007 1 1857 6480 t 10 R f (\(1979\), 368-375.)1 674 1 2889 6480 t (Linear Algebra)1 606 1 360 7500 t (\320 6 \320)2 300 1 2730 7620 t (BPCE)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 ( BPCE)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Example:)720 1320 w 10 R f (In the following example we obtain estimates of the condition numbers of the matrix)13 3388 1 1296 1320 t (1+x)2097 1680 w 10 S f (-)2403 1680 w 10 R f (1)2458 1680 w 10 S f (-)2042 1800 w 10 R f (1 2)1 411 1 2097 1800 t 10 S f (-)2658 1800 w 10 R f (1)2713 1800 w 10 S f (-)2403 1920 w 10 R f (1 2)1 305 1 2458 1920 t 10 S f (-)2913 1920 w 10 R f (1)2968 1920 w 10 S f (-)2658 2040 w 10 R f (1 2)1 305 1 2713 2040 t 10 S f (-)3168 2040 w 10 R f (1)3223 2040 w (.)2953 2160 w (.)3208 2280 w 10 S f (-)3423 2400 w 10 R f (1 2)1 305 1 3478 2400 t 10 S f (-)3933 2400 w 10 R f (1)3988 2400 w 10 S f (-)3678 2520 w 10 R f (1 1+x)1 411 1 3733 2520 t ( becomes more ill-)3 821( matrix is singular when x is 0, and)8 1594( The)1 229( and .0001.)2 492( .01,)1 175(with x=1.,)1 433 6 1296 2760 t (conditioned as x approaches 0.)4 1228 1 1296 2880 t ( condition number)2 737(In this example we compare the)5 1277 2 1296 3120 t 10 I f (K)3337 3120 w 10 S f ( \357)1 57(= \357)1 153 2 3453 3120 t 10 R f (A)3671 3120 w 10 S f ( \357)1 57( \357)1 90(\357 \357)1 106 3 3751 3120 t 10 R f (A)4012 3120 w 7 S f (-)4095 3080 w 7 R f (1)4145 3080 w 10 S f (\357 \357)1 106 1 4220 3120 t 10 R f (with the estimate)2 687 1 4353 3120 t ( the code below A)4 770( In)1 145(obtained from BPCE.)2 888 3 1296 3240 t 7 S f (-)3110 3200 w 7 R f (1)3160 3200 w 10 R f (is computed one column at a time and)7 1595 1 3239 3240 t 10 I f (K)4870 3240 w 10 R f (is)4973 3240 w ( the 1-norm,)2 495( In)1 135(computed using the 1-norm.)3 1133 3 1296 3360 t 10 S f (\357 \357)1 106 1 3086 3360 t 10 I f (A)3200 3360 w 10 S f (\357 \357)1 106 1 3269 3360 t 10 R f ( sum in absolute)3 662(is the maximum column)3 976 2 3402 3360 t (value, which is obviously 4.)4 1122 1 1296 3480 t 8 CW f (INTEGER N, MU, IG, K, I, IWRITE, I1MACH, J)8 2016 1 2040 3820 t (REAL G\(2,100\), B\(200\))2 1008 1 2040 3920 t (REAL X, COND, AINVNO, AINORM, ABS)5 1584 1 2040 4020 t (N=100)2040 4120 w (X=1.0)2040 4220 w (MU=2)2040 4320 w (IG=2)2040 4420 w (DO 50 K=1,3)2 528 1 2040 4520 t (C CONSTRUCT MATRIX)2 864 1 1656 4620 t (DO 10 I=1,N)2 528 1 2232 4720 t (G\(1,I\)=2.0)2376 4820 w (G\(2,I\)=-1.0)2376 4920 w (10 CONTINUE)1 864 1 1752 5020 t (G\(1,1\)=1.0+X)2232 5120 w (G\(1,N\)=1.0+X)2232 5220 w (C GET ESTIMATE OF CONDITION NUMBER FROM BPCE)7 2112 1 1656 5320 t (CALL BPCE\(N,MU,G,IG,COND\))1 1200 1 2232 5420 t (IWRITE=I1MACH\(2\))2232 5520 w (WRITE\(IWRITE,11\)X)2232 5620 w ( WHEN X IS,E14.6\))3 816(11 FORMAT\(/10H)1 1008 2 1752 5720 t (WRITE\(IWRITE,12\)COND)2232 5820 w ( ,E15.8\))1 480( CONDITION ESTIMATE IS)3 1056(12 FORMAT\(25H)1 960 3 1752 5920 t (C SINCE CONDITION NUMBER IS NORM\(A\)*NORM\(INVERSE\(A\)\),)5 2544 1 1656 6020 t (C FIND THE NORM OF EACH COLUMN OF INVERSE\(A\). GENERATE)9 2592 1 1656 6120 t (C THE COLUMNS ONE AT A TIME AND REUSE SPACE)9 2064 1 1656 6220 t (AINVNO=0.0)2232 6320 w (DO 40 I=1,N)2 528 1 2232 6420 t (C GENERATE ITH COLUMN OF IDENTITY MATRIX AS RIGHT HAND SIDE)10 2832 1 1656 6520 t (DO 20 J=1,N)2 528 1 2424 6620 t (B\(J\)=0.0)2568 6720 w (20 CONTINUE)1 1056 1 1752 6820 t (B\(I\)=1.0)2424 6920 w 10 R f (Linear Algebra)1 606 1 4794 7580 t (\320 7 \320)2 300 1 2730 7700 t (BPCE)5149 7820 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(BPCE February)1 4665 2 360 840 t 8 CW f (C SOLVE AX=B TO GET ITH COLUMN OF A\(INVERSE\))8 2112 1 1656 1300 t (CALL BPFS\(N,MU,G,IG,B,N,1\))1 1248 1 2424 1400 t (CALL BPBS\(N,MU,G,IG,B,N,1\))1 1248 1 2424 1500 t (C FIND NORM OF COLUMN)4 1008 1 1656 1600 t (AINORM=0.0)2424 1700 w (DO 30 J=1,N)2 528 1 2424 1800 t (AINORM=AINORM+ABS\(B\(J\)\))2568 1900 w (30 CONTINUE)1 1056 1 1752 2000 t (IF\(AINVNO.LT.AINORM\)AINVNO=AINORM)2424 2100 w (40 CONTINUE)1 864 1 1752 2200 t (COND=4.0*AINVNO)2232 2300 w (WRITE\(IWRITE,41\)COND)2232 2400 w ( TRUE CONDITION NUMBER IS,E15.8\))4 1536(41 FORMAT\(25H)1 960 2 1752 2500 t (X=X/100.0)2232 2600 w (50 CONTINUE)1 720 1 1752 2700 t (STOP)2088 2800 w (END)2088 2900 w 10 R f ( was executed on the Honeywell 6000 machine at Bell Laboratories, the)11 2883(When the above code)3 861 2 1296 3260 t (following was printed:)2 905 1 1296 3380 t 8 CW f ( 01)1 144( 0.100000E)1 528(WHEN X IS)2 432 3 1704 3720 t ( 04)1 144( 0.40807862E)1 720(CONDITION ESTIMATE IS)2 1008 3 1704 3820 t (TRUE CONDITION NUMBER IS 0.50999824E 04)5 1872 1 1704 3920 t ( 0.100000E-01)1 672(WHEN X IS)2 432 2 1704 4120 t ( 05)1 144( 0.23329148E)1 720(CONDITION ESTIMATE IS)2 1008 3 1704 4220 t (TRUE CONDITION NUMBER IS 0.24899523E 05)5 1872 1 1704 4320 t ( 0.100000E-03)1 672(WHEN X IS)2 432 2 1704 4520 t ( 07)1 144( 0.19933923E)1 720(CONDITION ESTIMATE IS)2 1008 3 1704 4620 t (TRUE CONDITION NUMBER IS 0.19950487E 07)5 1872 1 1704 4720 t 10 R f ( condition)1 412(The comparison above of the condition number estimated by BPCE with the true)12 3332 2 1296 5080 t ( that the order of magnitude \(which is all one usually is interested in\) of the)15 3055(number indicates)1 689 2 1296 5200 t ( the inverse of a band matrix is usually a full)10 1791( that)1 175( Note)1 244(estimated condition number is correct.)4 1534 4 1296 5320 t 10 I f (n)1296 5440 w 10 S f (\264)1354 5440 w 10 I f (n)1417 5440 w 10 R f (matrix and should rarely be calculated.)5 1552 1 1492 5440 t (Linear Algebra)1 606 1 360 7500 t (\320 8 \320)2 300 1 2730 7620 t (BPCE)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 ( BPCE)1 4305(February 11, 1993)2 735 2 360 840 t (BPDC \320 LDL)2 606 1 1583 1320 t 7 R f (T)2194 1280 w 10 R f (decomposition of a band symmetric positive de\256nite matrix A)8 2482 1 2270 1320 t 9 B f (Purpose:)720 1680 w 10 R f (BPDC \(Banded symmetric Positive de\256nite matrix DeComposition\) decomposes a banded)9 3744 1 1296 1680 t (symmetric positive de\256nite matrix A into LDL)6 1902 1 1296 1800 t 7 R f (T)3203 1760 w 10 R f ( \(1's on the)3 464(where L is lower unit triangular)5 1292 2 3284 1800 t ( step)1 189( is called by BPLE as the \256rst)7 1200( It)1 113(diagonal and 0's above the diagonal\) and D is diagonal.)9 2242 4 1296 1920 t (of the solution of a banded symmetric positive de\256nite linear system.)10 2762 1 1296 2040 t 9 B f (Usage:)720 2400 w 10 R f (CALL BPDC \(N, MU, G, IG\))5 1199 1 1296 2400 t (N)1440 2640 w 14 S f (\256)1980 2660 w 10 R f (the number of equations)3 968 1 2196 2640 t (MU)1440 2880 w 14 S f (\256)1980 2900 w 10 R f (the number of nonzero bands on and above the diagonal of A)11 2442 1 2196 2880 t (G)1440 3120 w 14 S f (\256)1980 3140 w 10 R f ( which the upper triangular portion of the matrix A has)10 2311(a matrix into)2 533 2 2196 3120 t (been packed as follows:)3 956 1 2196 3240 t (G\(j)2916 3420 w 10 S f (-)3049 3420 w 10 R f (i + 1, i\) =)4 376 1 3104 3420 t 10 I f (a)3505 3420 w 7 I f (i j)1 45 1 3566 3440 t 10 R f (for j)1 169 1 3644 3420 t 10 S f (\263)3813 3420 w 10 R f (i)3868 3420 w (\(See the introduction to the chapter.\))5 1464 1 2196 3600 t ( in the calling program, where)5 1348(G should be dimensioned \(IG,KG\))4 1496 2 2196 3720 t (IG)2196 3840 w 10 S f (\263)2301 3840 w 10 R f (MU and KG)2 499 1 2356 3840 t 10 S f (\263)2855 3840 w 10 R f (N.)2910 3840 w 14 S f (\254)1980 4100 w 10 R f (LDL)2196 4080 w 7 R f (T)2395 4040 w 10 R f ( \(see)1 188(decomposition suitable for input into BPFS and BPBS)7 2192 2 2473 4080 t 8 B f (Note)4881 4080 w (1)2196 4200 w 10 R f (\))2236 4200 w (IG)1440 4440 w 14 S f (\256)1980 4460 w 10 R f (the row \(leading\) dimension of G, as dimensioned in the)9 2253 1 2196 4440 t (calling program)1 635 1 2196 4560 t 9 B f (Note 1:)1 278 1 720 4920 t 10 R f (The LDL)1 378 1 1296 4920 t 7 R f (T)1679 4880 w 10 R f (decomposition of A satis\256es the equation A = LDL)8 2081 1 1759 4920 t 7 R f (T)3845 4880 w 10 R f (where L is lower unit trian-)5 1115 1 3925 4920 t ( D is diagonal. On return from BPDC,)7 1538(gular \(1's on the diagonal, 0's above the diagonal\) and)9 2206 2 1296 5040 t (the diagonal of D occupies the \256rst row of G and)10 1952 1 1296 5160 t 10 I f (G)3273 5160 w 10 R f (\()3353 5160 w 10 I f (i)3394 5160 w 10 S f (-)3430 5160 w 10 I f (j)3501 5160 w 10 S f (+)3545 5160 w 10 R f (1 ,)1 83 1 3616 5160 t 10 I f (i)3707 5160 w 10 R f (\))3743 5160 w 10 S f (=)3833 5160 w 10 I f (l)3937 5160 w 7 I f (i j)1 45 1 3976 5180 t 10 R f (for i)1 169 1 4054 5160 t 10 S f (>)4223 5160 w 10 R f (j.)4278 5160 w 9 B f (Note 2:)1 278 1 720 5520 t 10 R f ( represents the conjugate transpose of)5 1549( *)1 58( where A)2 385( *,)1 83(For complex Hermitian matrices \(A = A)6 1669 5 1296 5520 t ( decomposition and returns)3 1118( *)1 58( subroutine computes the LDL)4 1265(A\), the complex version of this)5 1303 4 1296 5640 t (the conjugate of L, rather than L, in G.)8 1544 1 1296 5760 t (Linear Algebra)1 606 1 4794 7500 t (\320 9 \320)2 300 1 2730 7620 t (BPDC)5138 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(BPDC 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 (2 MU)1 1097 1 1944 2160 t 10 S f (<)3066 2160 w 10 R f (1)3146 2160 w (3 IG)1 1041 1 1944 2400 t 10 S f (<)3010 2400 w 10 R f (MU)3090 2400 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2640 t ( k)1 75( the)1 574(10 + N + k*)4 484 3 1944 2880 t 7 R f (th)3082 2840 w 10 R f (principal minor is not positive de\256nite)5 1531 1 3170 2880 t 9 B f (Double-precision version:)1 988 1 720 3240 t 10 R f (DBPDC, with G declared double precision.)5 1734 1 1783 3240 t 9 B f (Complex Hermitian version:)2 1101 1 720 3600 t 10 R f (CBPDC with G declared complex \(see)5 1545 1 1896 3600 t 8 B f (Note 2)1 219 1 3466 3600 t 10 R f (\).)3685 3600 w 9 B f (Storage:)720 3960 w 10 R f (None)1296 3960 w 9 B f (Time:)720 4320 w 10 R f (N)1296 4320 w 10 S f (\264)1368 4320 w 10 R f (\(MU)1423 4320 w 10 S f (-)1617 4320 w 10 R f (1\))1672 4320 w 10 S f (\264)1755 4320 w 10 R f (MU/2 + N)2 417 1 1835 4320 t 10 S f (\264)2252 4320 w 10 R f (\(2MU)2307 4320 w 10 S f (-)2551 4320 w 10 R f (1\) additions)1 475 1 2606 4320 t (N)1296 4440 w 10 S f (\264)1368 4440 w 10 R f (\(MU)1423 4440 w 10 S f (-)1617 4440 w 10 R f (1\))1672 4440 w 10 S f (\264)1755 4440 w 10 R f (MU/2 multiplications)1 865 1 1835 4440 t (N)1296 4560 w 10 S f (\264)1368 4560 w 10 R f (\(MU)1423 4560 w 10 S f (-)1617 4560 w 10 R f (1\) divisions)1 470 1 1672 4560 t 9 B f (Method:)720 4920 w 10 R f ( =)1 81(BPDC calls BPLD after setting EPS)5 1488 2 1296 4920 t 10 S f (\357\357)2890 4920 w 10 R f (A)2988 4920 w 10 S f (\357\357e)3060 4920 w 10 R f (, where)1 301 1 3202 4920 t 10 S f (e)3536 4920 w 10 R f (is machine precision, i.e. the value)5 1426 1 3614 4920 t (returned by R1MACH\(4\) \(or, for double precision, by D1MACH\(4\)\).)8 2781 1 1296 5040 t 9 B f (See also:)1 333 1 720 5400 t 10 R f (BPBS, BPCE, BPFS, BPLD, BPLE, BPSS)5 1718 1 1296 5400 t 9 B f (Author:)720 5760 w 10 R f (Linda Kaufman)1 629 1 1296 5760 t 9 B f (Example:)720 6120 w 10 R f ( BPLD, and BPCE as a)5 948(This example is designed to indicate the relative ef\256ciency of BPDC,)10 2796 2 1296 6120 t ( the band. All three subroutines compute the same factorization, but)10 2767(function of the width of)4 977 2 1296 6240 t ( all the subroutines)3 785( In)1 142( the three.)2 414(the criterion for singularity is treated differently in each of)9 2403 4 1296 6360 t ( singular if for some i,)5 903(the matrix is considered)3 969 2 1296 6480 t 10 I f (d)3196 6480 w 7 I f (i)3257 6500 w 10 S f (\243)3313 6480 w 10 R f (EPS. In BPLD the user provides EPS; in)7 1644 1 3396 6480 t (BPDC the subroutine computes EPS \(see)5 1677 1 1296 6600 t 8 B f (Method)3004 6600 w 10 R f ( \(However,)1 480( as EPS.)2 345(\); in BPCE, 0.0 is used)5 945 3 3270 6600 t (BPCE also provides an estimate of the condition number of the matrix.\))11 2870 1 1296 6720 t (Linear Algebra)1 606 1 360 7500 t (\320 10 \320)2 350 1 2705 7620 t (BPDC)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 ( BPDC)1 4305(February 11, 1993)2 735 2 360 840 t ( example is derived from the traditional 5-point approximation to the La-)11 3026(The matrix in the)3 718 2 1296 1320 t ( N)1 97( The)1 205(place operator on the unit square.)5 1330 3 1296 1440 t 10 S f (\264)2928 1440 w 10 R f (N matrix A has the form)5 979 1 2983 1440 t (C)2408 1560 w 10 S f (-)2636 1560 w 10 R f (I)2691 1560 w 10 S f (-)2398 1680 w 10 R f (I C)1 260 1 2453 1680 t 10 S f (-)2930 1680 w 10 R f (I)2985 1680 w 10 S f (-)2636 1800 w 10 R f (I C)1 316 1 2691 1800 t 10 S f (-)3224 1800 w 10 R f (I)3279 1800 w (........)2874 1920 w 10 S f (-)3224 2040 w 10 R f (I C)1 260 1 3279 2040 t 10 S f (-)3700 2040 w 10 R f (I)3755 2040 w 10 S f (-)3462 2160 w 10 R f (I C)1 260 1 3517 2160 t (where I is the identity matrix and C is the matrix)10 1943 1 1296 2400 t (4)2458 2640 w 10 S f (-)2658 2640 w 10 R f (1)2713 2640 w 10 S f (-)2403 2760 w 10 R f (1 4)1 305 1 2458 2760 t 10 S f (-)2913 2760 w 10 R f (1)2968 2760 w 10 S f (-)2658 2880 w 10 R f (1 4)1 305 1 2713 2880 t 10 S f (-)3168 2880 w 10 R f (1)3223 2880 w (...)3183 3000 w 10 S f (-)3168 3120 w 10 R f (1 4)1 305 1 3223 3120 t 10 S f (-)3678 3120 w 10 R f (1)3733 3120 w 10 S f (-)3423 3240 w 10 R f (1 4)1 305 1 3478 3240 t ( a)1 73( For)1 193(In this problem we vary the sizes of the blocks C and I which changes the bandwidth.)16 3478 3 1296 3480 t (given blocksize, the number of blocks is varied which changes N.)10 2623 1 1296 3600 t ( counts)1 287( It)1 112( is a timer on the Honeywell 6000 with about 1% accuracy.)11 2379(The subroutine ILAPSZ)2 966 4 1296 3840 t (in 1/64 milliseconds.)2 837 1 1296 3960 t 8 CW f (INTEGER IG, MLM1, IWRITE, I1MACH, K, N, MU)7 2016 1 1992 4420 t (INTEGER NBLOK, KBLOK, KK, I, J, IT, ILAPSZ, IT2)8 2256 1 1992 4520 t (REAL G\(17, 100\), G2\(17, 100\), G3\(17, 100\))6 1968 1 1992 4620 t (REAL COND)1 432 1 1992 4720 t (IG=17)1992 4820 w (MLM1=4)1992 4920 w (IWRITE=I1MACH\(2\))1992 5020 w ( K=1,3)1 336(DO 70)1 240 2 1992 5120 t (DO 60 N=48,96,48)2 768 1 2136 5220 t (MU=MLM1+1)2280 5320 w (I=0)2280 5420 w (NBLOK=N/MLM1)2280 5520 w (C)1656 5620 w (C SET UP THREE MATRICES FOR ELLIPTIC PDE IN 2 DIMENSION)10 2640 1 1656 5720 t (C)1656 5820 w (DO 30 KBLOK=1,NBLOK)2 912 1 2280 5920 t (DO 20 KK=1,MLM1)2 720 1 2424 6020 t (I=I+1)2568 6120 w (G\(1,I\)=4.0)2568 6220 w (G\(2,I\)=-1.0)2568 6320 w (G\(MU,I\)=-1.0)2568 6420 w (DO 10 J=3,MLM1)2 672 1 2568 6520 t (G\(J,I\)=0.0)2712 6620 w (10 CONTINUE)1 1200 1 1752 6720 t (20 CONTINUE)1 1056 1 1752 6820 t (G\(2,I\)=0.0)2424 6920 w 10 R f (Linear Algebra)1 606 1 4794 7580 t (\320 11 \320)2 350 1 2705 7700 t (BPDC)5138 7820 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(BPDC February)1 4665 2 360 840 t 8 CW f (30 CONTINUE)1 912 1 1752 1300 t (DO 50 I=1,N)2 528 1 2280 1400 t (DO 40 J=1,MU)2 576 1 2424 1500 t (G2\(J,I\)=G\(J,I\))2568 1600 w (G3\(J,I\)=G\(J,I\))2568 1700 w (40 CONTINUE)1 1056 1 1752 1800 t (50 CONTINUE)1 912 1 1752 1900 t (WRITE\(IWRITE,51\)N,MU)2280 2000 w ( N IS ,I4,30H ,NUMBER OF UPPER DIAGONALS IS,I3\))8 2256(51 FORMAT\(/6H)1 1008 2 1752 2100 t (C TIME DECOMPOSITION BY BPLD)4 1344 1 1656 2200 t (IT=ILAPSZ\(0\))2280 2300 w (CALL BPLD\(N,MU,G,IG,0.0\))1 1152 1 2280 2400 t (IT=ILAPSZ\(0\)-IT)2280 2500 w (WRITE\(IWRITE,52\)IT)2280 2600 w ( TIME FOR BPLD,I7\))3 864(52 FORMAT\(14H)1 1008 2 1752 2700 t (C TIME DECOMPOSITION BY BPDC)4 1344 1 1656 2800 t (IT2=ILAPSZ\(0\))2280 2900 w (CALL BPDC\(N,MU,G2,IG\))1 1008 1 2280 3000 t (IT2=ILAPSZ\(0\)-IT2)2280 3100 w (WRITE\(IWRITE,53\)IT2)2280 3200 w ( TIME FOR BPDC,I7\))3 864(53 FORMAT\(14H)1 1008 2 1752 3300 t (C TIME DECOMPOSITION BY BPCE)4 1344 1 1656 3400 t (IT3=ILAPSZ\(0\))2280 3500 w (CALL BPCE\(N,MU,G3,IG,COND\))1 1248 1 2280 3600 t (IT3=ILAPSZ\(0\)-IT3)2280 3700 w (WRITE\(IWRITE,54\)IT3)2280 3800 w ( TIME FOR BPCE,I7\))3 864(54 FORMAT\(14H)1 1008 2 1752 3900 t (60 CONTINUE)1 720 1 1752 4000 t (MLM1=MLM1*2)2088 4100 w (70 CONTINUE)1 576 1 1752 4200 t (STOP)1944 4300 w (END)1944 4400 w 10 R f ( at Bell Laboratories with an optimiz-)6 1540(When the above code was run on the Honeywell 6000)9 2204 2 1296 4760 t (ing compiler the following was printed)5 1557 1 1296 4880 t 8 CW f ( 5)1 144( ,NUMBER OF UPPER DIAGONALS IS)5 1440( 48)1 240(N IS)1 192 4 1704 5300 t ( 959)1 336(TIME FOR BPLD)2 624 2 1704 5400 t ( 1413)1 336(TIME FOR BPDC)2 624 2 1704 5500 t ( 3670)1 336(TIME FOR BPCE)2 624 2 1704 5600 t ( 5)1 144( ,NUMBER OF UPPER DIAGONALS IS)5 1440( 96)1 240(N IS)1 192 4 1704 5800 t ( 1868)1 336(TIME FOR BPLD)2 624 2 1704 5900 t ( 2673)1 336(TIME FOR BPDC)2 624 2 1704 6000 t ( 7965)1 336(TIME FOR BPCE)2 624 2 1704 6100 t ( 9)1 144( ,NUMBER OF UPPER DIAGONALS IS)5 1440( 48)1 240(N IS)1 192 4 1704 6300 t ( 2395)1 336(TIME FOR BPLD)2 624 2 1704 6400 t ( 3181)1 336(TIME FOR BPDC)2 624 2 1704 6500 t ( 6377)1 336(TIME FOR BPCE)2 624 2 1704 6600 t ( 9)1 144( ,NUMBER OF UPPER DIAGONALS IS)5 1440( 96)1 240(N IS)1 192 4 1704 6800 t ( 4854)1 336(TIME FOR BPLD)2 624 2 1704 6900 t 10 R f (Linear Algebra)1 606 1 360 7560 t (\320 12 \320)2 350 1 2705 7680 t (BPDC)360 7800 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 ( BPDC)1 4305(February 11, 1993)2 735 2 360 840 t 8 CW f ( 6450)1 336(TIME FOR BPDC)2 624 2 1704 1300 t ( 13170)1 336(TIME FOR BPCE)2 624 2 1704 1400 t ( ,NUMBER OF UPPER DIAGONALS IS 17)6 1584( 48)1 240(N IS)1 192 3 1704 1600 t ( 7143)1 336(TIME FOR BPLD)2 624 2 1704 1700 t ( 8227)1 336(TIME FOR BPDC)2 624 2 1704 1800 t ( 12722)1 336(TIME FOR BPCE)2 624 2 1704 1900 t ( ,NUMBER OF UPPER DIAGONALS IS 17)6 1584( 96)1 240(N IS)1 192 3 1704 2100 t ( 15001)1 336(TIME FOR BPLD)2 624 2 1704 2200 t ( 17160)1 336(TIME FOR BPDC)2 624 2 1704 2300 t ( 49809)1 336(TIME FOR BPCE)2 624 2 1704 2400 t 10 R f ( example indicates, the cost of computing the condition estimate can be substantially)12 3479(As the)1 265 2 1296 2740 t ( of)1 108( ratio)1 209( The)1 206(greater than the cost of factoring the matrix when the width of the band is small.)15 3221 4 1296 2860 t ( general lin-)2 486(BPCE to BPDC does not always decrease as N increases, as it does in the case of)16 3258 2 1296 2980 t ( BPCE to prevent)3 730(ear systems, but depends rather on the scaling that is sometimes done by)12 3014 2 1296 3100 t ( no scaling is done, the ratio of)7 1281( If)1 123( the calculation of the condition estimate.)6 1693(over\257ow during)1 647 4 1296 3220 t ( to BPDC remains constant, but for some examples most of the time is)13 2933(the times for BPCE)3 811 2 1296 3340 t (spent scaling and the time ratio increases as N increases.)9 2256 1 1296 3460 t (Linear Algebra)1 606 1 4794 7500 t (\320 13 \320)2 350 1 2705 7620 t (BPDC)5138 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(BPDC February)1 4665 2 360 840 t (BPFS \320 band symmetric lower \(unit\) triangular linear system solution)9 2841 1 1747 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( where L)2 358( = B)2 173( de\256nite matrix Forward Solution\) solves LX)6 1817(BPFS \(Banded symmetric Positive)3 1396 4 1296 1680 t ( unit lower triangular matrix, \(i.e. 1's on the diagonal, 0's above the diagonal\).)13 3251( banded)1 347(is a)1 146 3 1296 1800 t ( be used for the forward solution phase of a band symmetric positive de\256nite linear)14 3345(BPFS can)1 399 2 1296 1920 t ( is used in this way by the routines BPSS and BPLE.\))11 2140(system. \(It)1 447 2 1296 2040 t 9 B f (Usage:)720 2400 w 10 R f (CALL BPFS \(N, ML, G, IG, B, IB, NB\))8 1617 1 1296 2400 t (N)1440 2640 w 14 S f (\256)1980 2660 w 10 R f (the number of equations)3 968 1 2196 2640 t (ML)1440 2880 w 14 S f (\256)1980 2900 w 10 R f (the number of nonzero diagonals on and below)7 1877 1 2196 2880 t (the diagonal of L)3 685 1 2196 3000 t (G)1440 3240 w 14 S f (\256)1980 3260 w 10 R f ( results obtained by the routines)5 1414(a matrix \(which may contain the)5 1430 2 2196 3240 t (BPCE, BPDC, or BPLD\) into which L is packed as follows:)10 2406 1 2196 3360 t (G\(1+i)2916 3540 w 10 S f (-)3155 3540 w 10 R f (j, i\) =)2 220 1 3210 3540 t 10 I f (l)3455 3540 w 7 I f (i j)1 45 1 3494 3560 t 10 R f (for i)1 169 1 3572 3540 t 10 S f (>)3741 3540 w 10 R f (j)3796 3540 w (\(See the introduction to this chapter.\))5 1487 1 2196 3720 t ( in the calling program, where)5 1348(G should be dimensioned \(IG,KG\))4 1496 2 2196 3840 t (IG)2196 3960 w 10 S f (\263)2301 3960 w 10 R f (ML and KG)2 488 1 2356 3960 t 10 S f (\263)2844 3960 w 10 R f (N.)2899 3960 w (IG)1440 4200 w 14 S f (\256)1980 4220 w 10 R f (the row \(leading\) dimension of G, as dimensioned in the)9 2253 1 2196 4200 t (calling program)1 635 1 2196 4320 t (B)1440 4560 w 14 S f (\256)1980 4580 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 4560 t (the calling program, where IB)4 1200 1 2196 4680 t 10 S f (\263)3396 4680 w 10 R f (N and KB)2 405 1 3451 4680 t 10 S f (\263)3856 4680 w 10 R f (NB.)3911 4680 w 14 S f (\254)1980 4940 w 10 R f (the solution X)2 567 1 2196 4920 t (IB)1440 5160 w 14 S f (\256)1980 5180 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 5160 t (calling program)1 635 1 2196 5280 t (NB)1440 5520 w 14 S f (\256)1980 5540 w 10 R f (the number of right-hand sides)4 1226 1 2196 5520 t 9 B f (Note 1:)1 278 1 720 5880 t 10 R f ( can be used directly on the output matrix produced by BPDC, BPLD, or)13 3049(BPFS and BPBS)2 695 2 1296 5880 t (BPCE to solve a general linear system.)6 1557 1 1296 6000 t (Linear Algebra)1 606 1 360 7500 t (\320 14 \320)2 350 1 2705 7620 t (BPFS)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 ( BPFS)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 ( same coef\256cient matrix, but differ-)5 1434(Users who have to solve a sequence of problems with the)10 2310 2 1296 1320 t (ent right-hand sides,)2 815 1 1296 1440 t 10 I f ( advance,)1 382(not all known in)3 651 2 2137 1440 t 10 R f (should not call BPSS or BPLE repeatedly, but)7 1845 1 3195 1440 t (should use BPDC to \256nd the LDL)6 1404 1 1296 1560 t 7 R f (T)2705 1520 w 10 R f ( then)1 205(decompostion of the coef\256cient matrix, which can)6 2047 2 2788 1560 t ( solutions \(using BPFS\) and)4 1148(be used repeatedly \(and ef\256ciently\) for the sequence of forward)9 2596 2 1296 1680 t (back solutions \(using BPBS\) for each set of right-hand sides.)9 2437 1 1296 1800 t 9 B f (Error situations:)1 648 1 720 2160 t 10 R f (\(All errors in this subprogram are fatal \320)7 1666 1 1512 2160 t (see)1512 2280 w 10 I f (Error Handling)1 631 1 1664 2280 t 10 R f (, Framework Chapter\))2 884 1 2295 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 ML)1 1086 1 1944 3000 t 10 S f (<)3055 3000 w 10 R f (1)3135 3000 w (3 IG)1 1041 1 1944 3240 t 10 S f (<)3010 3240 w 10 R f (ML)3090 3240 w (4 IB)1 1036 1 1944 3480 t 10 S f (<)3005 3480 w 10 R f (N)3085 3480 w (5 NB)1 1075 1 1944 3720 t 10 S f (<)3044 3720 w 10 R f (1)3124 3720 w 9 B f (Double-precision version:)1 988 1 720 4080 t 10 R f (DBPFS, with G and B declared double precision)7 1943 1 1800 4080 t 9 B f (Complex Hermitian version:)2 1101 1 720 4440 t 10 R f (CBPFS with G and B declared complex)6 1594 1 1896 4440 t (G should contain the conjugate of L)6 1437 1 2041 4560 t 9 B f (Storage:)720 4920 w 10 R f (None)1296 4920 w 9 B f (Time:)720 5280 w 10 R f (N)1296 5280 w 10 S f (\264)1368 5280 w 10 R f (\(ML)1423 5280 w 10 S f (-)1606 5280 w 10 R f (1\))1661 5280 w 10 S f (\264)1744 5280 w 10 R f (NB additions)1 531 1 1799 5280 t (N)1296 5400 w 10 S f (\264)1368 5400 w 10 R f (\(ML)1423 5400 w 10 S f (-)1606 5400 w 10 R f (1\))1661 5400 w 10 S f (\264)1744 5400 w 10 R f (NB multiplications)1 765 1 1799 5400 t 9 B f (See also:)1 333 1 720 5760 t 10 R f (BPBS, BPCE, BPDC, BPLD, BPLE, BPSS)5 1745 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 15 \320)2 350 1 2705 7620 t (BPFS)5165 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(BPFS February)1 4665 2 360 840 t 9 B f (Example:)720 1320 w 10 R f ( a discretization of a)4 829(In the following example we solve three problems that might arise from)11 2915 2 1296 1320 t (1-dimensional differential equation. The coef\256cient matrix in each problem has the form)11 3538 1 1296 1440 t (1+x)2097 1680 w 10 S f (-)2403 1680 w 10 R f (1)2458 1680 w 10 S f (-)2042 1800 w 10 R f (1 2)1 411 1 2097 1800 t 10 S f (-)2658 1800 w 10 R f (1)2713 1800 w 10 S f (-)2403 1920 w 10 R f (1 2)1 305 1 2458 1920 t 10 S f (-)2913 1920 w 10 R f (1)2968 1920 w 10 S f (-)2658 2040 w 10 R f (1 2)1 305 1 2713 2040 t 10 S f (-)3168 2040 w 10 R f (1)3223 2040 w (.)2953 2160 w (.)3208 2280 w 10 S f (-)3423 2400 w 10 R f (1 2)1 305 1 3478 2400 t 10 S f (-)3933 2400 w 10 R f (1)3988 2400 w 10 S f (-)3678 2520 w 10 R f (1 1+x)1 411 1 3733 2520 t ( errors in the)3 516( make it easy to detect)5 900( To)1 163(with x=1., .01, and .0001 de\256ning the three problems.)8 2165 4 1296 2760 t (solution, the right-hand side has been chosen to make the solution a vector of all 1's.)15 3382 1 1296 2880 t ( below is an encoding of an iterative re\256nement algorithm which may be used to)14 3226(The program)1 518 2 1296 3120 t ( an ill-conditioned ma-)3 928(obtain a highly accurate solution to a system of linear equations with)11 2816 2 1296 3240 t ( high, the iterative re\256nement algorithm)5 1633( the condition number is not excessively)6 1652(trix. When)1 459 3 1296 3360 t ( itera-)1 237( The)1 207( a solution that is accurate to the working precision of the machine.)12 2706(usually returns)1 594 4 1296 3480 t (tive re\256nement algorithm is essentially:)4 1583 1 1296 3600 t ( Ax = b)3 303(\(1\) Solve)1 408 2 1980 3840 t ( tol =)2 212(\(2\) Set)1 308 2 1980 3960 t 10 S f (e)2525 3960 w 15 S f (S)2626 3990 w 10 S f (\357)2764 3977 w 10 R f (x)2812 3960 w 7 I f (i)2873 3980 w 10 S f (\357)2901 3977 w 10 R f (where)2160 4090 w 10 S f (e)2428 4090 w 10 R f (is the precision of the machine)5 1223 1 2497 4090 t ( in double precision the residual, Ax)6 1451(\(3\) Compute)1 547 2 1980 4210 t 10 S f (-)3978 4210 w 10 R f (b,)4033 4210 w (and put it in the real vector r.)7 1159 1 2160 4330 t ( A)1 97(\(4\) Solve)1 408 2 1980 4450 t 10 S f (d)2510 4450 w 10 R f (x = r)2 189 1 2559 4450 t ( norm =)2 317(\(5\) Compute)1 547 2 1980 4570 t 15 S f (S)2885 4600 w 10 S f (\357)3023 4587 w (d)3071 4570 w 10 R f (x)3128 4570 w 7 I f (i)3189 4590 w 10 S f (\357)3217 4587 w 10 R f ( x to x +)4 334(\(6\) Set)1 308 2 1980 4700 t 10 S f (d)2647 4700 w 10 R f (x)2721 4700 w ( norm)1 236(\(7\) If)1 246 2 1980 4820 t 10 S f (\243)2487 4820 w 10 R f ( else return to step 3)5 807(tol stop,)1 348 2 2567 4820 t ( using the the linear equation solver for band positive)9 2218(In our code, step \(1\) is accomplished)6 1526 2 1296 5060 t ( leaves in the G matrix a decomposition which)8 1859(de\256nite matrices, BPLE. The subroutine BPLE)5 1885 2 1296 5180 t ( but dif-)2 330(may be used by BPFS and BPBS to solve problems with the same coef\256cient matrix)14 3414 2 1296 5300 t ( the matrix will be so ill-)6 1115( it is possible that)4 789( Since)1 294(ferent right-hand sides as in step\(4\).)5 1546 4 1296 5420 t ( diverge, steps \(3\) through \(7\) in our)7 1509(conditioned that the iterative re\256nement algorithm will)6 2235 2 1296 5540 t ( number is chosen to be an upper)7 1408( This)1 242(code are performed only a \256nite number of times.)8 2094 3 1296 5660 t ( the ma-)2 329(bound on number of of bits in the mantissa of the \257oating-point number supported by)14 3415 2 1296 5780 t (chine.)1296 5900 w (This algorithm is not included in)5 1394 1 1296 6140 t 8 R f (PORT)2733 6140 w 10 R f (because for double-precision matrices part of the)6 2058 1 2982 6140 t (computation would have to be done in extended precision.)8 2333 1 1296 6260 t (Linear Algebra)1 606 1 360 7500 t (\320 16 \320)2 350 1 2705 7620 t (BPFS)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 ( BPFS)1 4305(February 11, 1993)2 735 2 360 840 t (The following program has been speci\256cally tailored to the three problems.)10 3009 1 1296 1320 t 8 CW f (INTEGER N, ML, IG, NM1, K, I, IWRITE, I1MACH, IT, IEND, ITER)11 2880 1 1992 1660 t (REAL G\(2,100\), B\(200\), R\(200\))3 1392 1 1992 1760 t (REAL X, ERR, AMAX1, RNORM, BNORM, R1MACH, ABS)7 2160 1 1992 1860 t (DOUBLE PRECISION DBLE)2 1008 1 1992 1960 t (C CONSTRUCT MATRIX AND RIGHT HAND SIDE SO TRUE SOLUTION IS)10 2784 1 1656 2060 t (C COMPOSED ENTIRELY OF 1S)4 1200 1 1656 2160 t (N=100)1992 2260 w (X=1)1992 2360 w (ML=2)1992 2460 w (IG=2)1992 2560 w (NM1=N-1)1992 2660 w (DO 90 K=1,3)2 528 1 1992 2760 t (DO 10 I=1,N)2 528 1 2136 2860 t (G\(1,I\)=2.0)2280 2960 w (G\(2,I\)=-1.0)2280 3060 w (B\(I\)=0.0)2280 3160 w (10 CONTINUE)1 768 1 1752 3260 t (G\(1,1\)=1.0+X)2136 3360 w (G\(1,N\)=1.0+X)2136 3460 w (B\(1\)=X)2136 3560 w (B\(N\)=X)2136 3660 w (C SOLVE THE SYSTEM)3 864 1 1656 3760 t (CALL BPLE\(N,ML,G,IG,B,N,1\))1 1248 1 2136 3860 t (IWRITE=I1MACH\(2\))2136 3960 w (WRITE\(IWRITE,11\)X)2136 4060 w ( X IS,F16.8\))2 576(11 FORMAT\(/5H)1 864 2 1752 4160 t (C COMPUTE THE ERROR)3 912 1 1656 4260 t (ERR=0.0)2136 4360 w (DO 20 I=1,N)2 528 1 2136 4460 t (ERR=AMAX1\(ERR,ABS\(B\(I\)-1.0\)\))2280 4560 w (20 CONTINUE)1 768 1 1752 4660 t (WRITE\(IWRITE,21\)ERR)2136 4760 w ( FOR BPLE THE ERROR IS,F16.8\))5 1392(21 FORMAT\(22H)1 864 2 1752 4860 t (IEND=I1MACH\(11\)*IFIX\(R1MACH\(5\)/ALOG10\(2.0\)+1.0\))2136 4960 w (C FIND THE NORM OF THE SOLUTION)6 1488 1 1656 5060 t (BNORM=0.0)2136 5160 w (DO 30 I=1,N)2 528 1 2136 5260 t (BNORM=AMAX1\(BNORM,ABS\(B\(I\)\)\))2328 5360 w (30 CONTINUE)1 768 1 1752 5460 t (C REFINE THE SOLUTION)3 1008 1 1656 5560 t (DO 60 ITER=1,IEND)2 816 1 2136 5660 t (IT=ITER)2280 5760 w (C COMPUTE THE RESIDUAL R=B-AX, IN DOUBLE PRECISION)7 2400 1 1656 5860 t (DO 40 I=2,NM1)2 624 1 2280 5960 t (R\(I\)=DBLE\(B\(I-1\)\)+DBLE\(B\(I+1\)\)-2.D0*DBLE\(B\(I\)\))2424 6060 w (40 CONTINUE)1 912 1 1752 6160 t (R\(1\)=X+B\(2\)-DBLE\(1.0+X\)*DBLE\(B\(1\)\))2280 6260 w (R\(N\)=X+B\(N-1\)-DBLE\(1.+X\)*DBLE\(B\(N\)\))2280 6360 w (C SOLVE A\(DELTAX\)=R)2 912 1 1656 6460 t (CALL BPFS\(N,ML,G,IG,R,N,1\))1 1248 1 2280 6560 t (CALL BPBS\(N,ML,G,IG,R,N,1\))1 1248 1 2280 6660 t (C DETERMINE NORM OF CORRECTION AND ADD IN CORRECTION)8 2496 1 1656 6760 t (RNORM=0.0)2280 6860 w 10 R f (Linear Algebra)1 606 1 4794 7520 t (\320 17 \320)2 350 1 2705 7640 t (BPFS)5165 7760 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(BPFS February)1 4665 2 360 840 t 8 CW f (DO 50 I=1,N)2 528 1 2280 1300 t (B\(I\)=B\(I\)+R\(I\))2424 1400 w (RNORM=RNORM+ABS\(R\(I\)\))2424 1500 w (50 CONTINUE)1 912 1 1752 1600 t (IF\(RNORM.LT.R1MACH\(4\)*BNORM\) GO TO 70)3 1776 1 2280 1700 t (60 CONTINUE)1 768 1 1752 1800 t (WRITE\(IWRITE,61\))2136 1900 w ( REFINEMENT FAILED\))2 912(61 FORMAT\(18H)1 864 2 1752 2000 t (C COMPUTE NEW ERROR)3 912 1 1656 2100 t (70 ERR=0.0)1 720 1 1752 2200 t (DO 80 I=1,N)2 528 1 2136 2300 t (ERR=AMAX1\(ERR,ABS\(B\(I\)-1.0\)\))2280 2400 w (80 CONTINUE)1 768 1 1752 2500 t (WRITE\(IWRITE,81\)IT,ERR)2136 2600 w ( ERROR AFTER REFINEMENT ,I4,3H IS,E14.7\))5 1920(81 FORMAT\(24H)1 864 2 1752 2700 t (X=X/100.0)2136 2800 w (90 CONTINUE)1 624 1 1752 2900 t (STOP)1992 3000 w (END)1992 3100 w 10 R f ( above program was executed on the Honeywell 6000 at Bell Laboratories, the fol-)13 3356(When the)1 388 2 1296 3440 t (lowing was printed)2 766 1 1296 3560 t 8 CW f ( 1.00000000)1 768(X IS)1 192 2 1704 3880 t ( 0.00000431)1 768(FOR BPLE THE ERROR IS)4 1008 2 1704 3980 t ( IS 0.)2 288( 2)1 240(ERROR AFTER REFINEMENT)2 1056 3 1704 4080 t ( 0.01000000)1 768(X IS)1 192 2 1704 4280 t ( 0.00002055)1 768(FOR BPLE THE ERROR IS)4 1008 2 1704 4380 t ( IS 0.)2 288( 3)1 240(ERROR AFTER REFINEMENT)2 1056 3 1704 4480 t ( 0.00010000)1 768(X IS)1 192 2 1704 4680 t ( 0.00491761)1 768(FOR BPLE THE ERROR IS)4 1008 2 1704 4780 t ( IS 0.)2 288( 4)1 240(ERROR AFTER REFINEMENT)2 1056 3 1704 4880 t 10 R f ( as x approaches 0, the matrix becomes more ill-)9 2250(This problem was chosen because)4 1494 2 1296 5220 t ( grows, and more steps of iterative re\256nement are re-)9 2219(conditioned, the error in the solution)5 1525 2 1296 5340 t ( represent the matrix)3 835(quired. The reader should be aware that in this example we were able to)13 2909 2 1296 5460 t ( case. Often)2 479(and the right-hand side precisely, but due to roundoff error this is not always the)14 3265 2 1296 5580 t ( useless, solution to a slightly incor-)6 1466(the iterative re\256nement algorithm produces an exact, but)7 2278 2 1296 5700 t (rect problem.)1 532 1 1296 5820 t (Linear Algebra)1 606 1 360 7500 t (\320 18 \320)2 350 1 2705 7620 t (BPFS)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 ( BPFS)1 4305(February 11, 1993)2 735 2 360 840 t (BPLD \320 LDL)2 600 1 1635 1320 t 7 R f (T)2240 1280 w 10 R f (decomposition of a band symmetric positive de\256nite matrix)7 2385 1 2316 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( LDL)1 225(BPLD \(Band Positive de\256nite)3 1213 2 1296 1680 t 7 R f (T)2739 1640 w 10 R f (decomposition\) decomposes a banded symmetric posi-)5 2219 1 2821 1680 t (tive de\256nite matrix A into LDL)5 1268 1 1296 1800 t 7 R f (T)2569 1760 w 10 R f ( allows the)2 435( It)1 112(where L is lower triangular and D is diagonal.)8 1846 3 2647 1800 t ( a threshold for considering a matrix singular BPLD is called by the decompo-)13 3145(user to provide)2 599 2 1296 1920 t (sition routines BPCE and BPDC.)4 1327 1 1296 2040 t 9 B f (Usage:)720 2400 w 10 R f (CALL BPLD \(N, MU, G, IG, EPS\))6 1416 1 1296 2400 t (N)1440 2640 w 14 S f (\256)1980 2660 w 10 R f (the number of equations)3 968 1 2196 2640 t (MU)1440 2880 w 14 S f (\256)1980 2900 w 10 R f (the number of nonzero bands on and above the diagonal of A)11 2442 1 2196 2880 t (G)1440 3120 w 14 S f (\256)1980 3140 w 10 R f ( which the upper triangular portion of the matrix A has)10 2311(a matrix into)2 533 2 2196 3120 t (been packed as follows:)3 956 1 2196 3240 t (G\(j)2916 3420 w 10 S f (-)3049 3420 w 10 R f (i + 1, i\) =)4 376 1 3104 3420 t 10 I f (a)3505 3420 w 7 I f (i j)1 45 1 3566 3440 t 10 R f (for j)1 169 1 3644 3420 t 10 S f (\263)3813 3420 w 10 R f (i)3868 3420 w (\(See the introduction to this chapter.\))5 1487 1 2196 3600 t ( in the calling program, where)5 1348(G should be dimensioned \(IG,KG\))4 1496 2 2196 3720 t (IG)2196 3840 w 10 S f (\263)2301 3840 w 10 R f (MU and KG)2 499 1 2356 3840 t 10 S f (\263)2855 3840 w 10 R f (N.)2910 3840 w 14 S f (\254)1980 4100 w 10 R f (the LDL)1 346 1 2196 4080 t 7 R f (T)2547 4040 w 10 R f (decomposition suitable for input into BPFS and BPBS \(see)8 2411 1 2629 4080 t 8 B f (Note 2)1 219 1 2196 4200 t 10 R f (\))2415 4200 w (IG)1440 4440 w 14 S f (\256)1980 4460 w 10 R f (the row \(leading\) dimension of G, as dimensioned in the)9 2253 1 2196 4440 t (calling program)1 635 1 2196 4560 t (EPS)1440 4800 w 14 S f (\256)1980 4820 w 10 R f (if there exists an index)4 920 1 2196 4800 t 10 I f (k)3145 4800 w 10 R f (such that)1 362 1 3218 4800 t 10 S f (\357)3610 4800 w 10 I f (d)3689 4800 w 7 I f (kk)3750 4820 w 10 S f (\357 \243)1 129 1 3850 4800 t 10 R f (EPS then A is considered)4 1036 1 4004 4800 t (singular)2196 4920 w 9 B f (Note 1:)1 278 1 720 5280 t 10 R f ( singular\), the determinant of)4 1182(After the execution of BPLD, \(if the matrix has not been found)11 2562 2 1296 5280 t (A is the product of the elements of the \256rst row of G.)12 2122 1 1296 5400 t 9 B f (Note 2:)1 278 1 720 5760 t 10 R f (The LDL)1 378 1 1296 5760 t 7 R f (T)1679 5720 w 10 R f (decomposition of A satis\256es the equation A = LDL)8 2081 1 1759 5760 t 7 R f (T)3845 5720 w 10 R f (where L is lower unit trian-)5 1115 1 3925 5760 t ( BPLD,)1 308(gular \(1's on the diagonal, 0's above the diagonal\) and D is diagonal. On return from)15 3436 2 1296 5880 t (the diagonal of D occupies the \256rst row of G and G\(i)11 2110 1 1296 6000 t 10 S f (-)3406 6000 w 10 R f (j+1,i\) =)1 301 1 3461 6000 t 10 I f (l)3787 6000 w 7 I f (i j)1 45 1 3826 6020 t 10 R f (for i)1 169 1 3904 6000 t 10 S f (>)4073 6000 w 10 R f (j.)4128 6000 w 9 B f (Note 3:)1 278 1 720 6360 t 10 R f ( the conjugate transpose of A)5 1182( represents)1 431( *)1 58( where A)2 369( *,)1 83(For complex Hermitian matrices \(A = A)6 1621 6 1296 6360 t ( and returns the)3 630( decomposition)1 619( *)1 58(\), the complex version of this subroutine computes the LDL)9 2437 4 1296 6480 t (conjugate of L rather than L in G.)7 1347 1 1296 6600 t (Linear Algebra)1 606 1 4794 7500 t (\320 19 \320)2 350 1 2705 7620 t (BPLD)5144 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(BPLD 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 (2 MU)1 1097 1 1944 2160 t 10 S f (<)3066 2160 w 10 R f (1)3146 2160 w (3 IG)1 1041 1 1944 2400 t 10 S f (<)3010 2400 w 10 R f (MU)3090 2400 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2640 t ( the)1 574(10 + N + k*)4 484 2 1944 2880 t 10 I f (k)3027 2880 w 7 I f (th)3082 2840 w 10 R f (principal minor is not positive de\256nite)5 1531 1 3170 2880 t 9 B f (Double-precision version:)1 988 1 720 3240 t 10 R f (DBPLD with G and EPS declared double precision)7 2045 1 1783 3240 t 9 B f (Complex Hermitian version:)2 1101 1 720 3600 t 10 R f (CBPLD with G declared complex \(see)5 1539 1 1896 3600 t 8 B f (Note 3)1 219 1 3460 3600 t 10 R f (above\))3704 3600 w 9 B f (Storage:)720 3960 w 10 R f (None)1296 3960 w 9 B f (Timing:)720 4320 w 10 R f (N)1296 4320 w 10 S f (\264)1368 4320 w 10 R f (\(MU)1423 4320 w 10 S f (-)1617 4320 w 10 R f (1\))1672 4320 w 10 S f (\264)1755 4320 w 10 R f (MU/2 additions)1 631 1 1810 4320 t (N)1296 4440 w 10 S f (\264)1368 4440 w 10 R f (\(MU)1423 4440 w 10 S f (-)1617 4440 w 10 R f (1\))1672 4440 w 10 S f (\264)1755 4440 w 10 R f (MU/2 multiplications)1 865 1 1810 4440 t (\(MU)1296 4560 w 10 S f (-)1490 4560 w 10 R f (1\))1545 4560 w 10 S f (\264)1628 4560 w 10 R f (N divisions)1 459 1 1683 4560 t 9 B f (Method:)720 4920 w 10 R f (Gaussian elimination without pivoting)3 1537 1 1296 4920 t 9 B f (See also:)1 333 1 720 5280 t 10 R f (BPBS, BPCE, BPDC, BPFS, BPLE, BPSS)5 1724 1 1296 5280 t 9 B f (Author:)720 5640 w 10 R f (Linda Kaufman)1 629 1 1296 5640 t (Linear Algebra)1 606 1 360 7500 t (\320 20 \320)2 350 1 2705 7620 t (BPLD)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 ( BPLD)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Example:)720 1320 w 10 R f ( matrix stored in G is the)6 1060(As noted above, after execution of BPLD, the determinant of the)10 2684 2 1296 1320 t ( G. The subroutine computes this product)6 1717(product of the elements stored in the \256rst row of)9 2027 2 1296 1440 t ( subroutine UMKFL is used to decompose)6 1723( The)1 209(taking care to avoid under\257ow and over\257ow.)6 1812 3 1296 1560 t ( mantissa, M, and an exponent E such that F)9 1786(a \257oating-point number, F, into a)5 1340 2 1296 1680 t 10 S f (=)4471 1680 w 10 R f (Mb)4575 1680 w 7 R f (E)4719 1640 w 10 R f (where)4797 1680 w (b is the base of the machine and 1/b)8 1431 1 1296 1800 t 10 S f (\243 \357)1 129 1 2752 1800 t 10 R f (M)2881 1800 w 10 S f (\357 \243)1 129 1 2970 1800 t 10 R f (1)3124 1800 w 8 CW f (SUBROUTINE BPDET\(N,MU,G,IG,DETMAN,IDETEX\))1 1968 1 1992 2140 t (C)1656 2240 w (C THIS SUBROUTINE COMPUTES THE DETERMINANT OF A)7 2256 1 1656 2340 t (C BAND SYMMETRIC POSITIVE DEFINITE MATRIX STORED IN G.)8 2592 1 1656 2440 t (C IT IS GIVEN BY DETMAN*BETA**IDETEX)5 1728 1 1656 2540 t (C WHERE BETA IS THE BASE OF THE MACHINE)8 1872 1 1656 2640 t (C AND DETMAN IS BETWEEN 1/BETA AND 1)7 1728 1 1656 2740 t (C)1656 2840 w (REAL G\(IG,N\),DETMAN)1 912 1 1992 2940 t (REAL ONOVBE,M)1 624 1 1992 3040 t (INTEGER E)1 432 1 1992 3140 t (INTEGER IDETEX)1 672 1 1992 3240 t (CALL BPLD\(N,MU,G,IG,0.0\))1 1152 1 1992 3340 t (C)1656 3440 w (C THE DETERMINANT IS THE PRODUCT OF THE ELEMENTS OF ROW 1 OF G)13 2976 1 1656 3540 t (C WE TRY TO COMPUTE THIS PRODUCT IN A WAY THAT WILL)11 2448 1 1656 3640 t (C AVOID UNDERFLOW AND OVERFLOW)4 1440 1 1656 3740 t (C)1656 3840 w (ONOVBE=1.0/FLOAT\(I1MACH\(10\)\))1992 3940 w (DETMAN=ONOVBE)1992 4040 w (BETA=FLOAT\(I1MACH\(10\)\))1992 4140 w (IDETEX=1)1992 4240 w (DO 10 I=1,N)2 528 1 1992 4340 t (CALL UMKFL\(G\(1,I\),E,M\))1 1056 1 2184 4440 t (DETMAN=DETMAN*M)2184 4540 w (IDETEX=IDETEX+E)2184 4640 w (IF\(DETMAN.GE.ONOVBE\) GO TO 10)3 1392 1 2184 4740 t (IDETEX=IDETEX-1)2328 4840 w (DETMAN=DETMAN*BETA)2328 4940 w (10 CONTINUE)1 624 1 1752 5040 t (RETURN)1992 5140 w (END)1992 5240 w 10 R f (Linear Algebra)1 606 1 4794 7500 t (\320 21 \320)2 350 1 2705 7620 t (BPLD)5144 7740 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(BPLD February)1 4665 2 360 840 t (BPLE \320 band symmetric positive de\256nite linear system solution)8 2606 1 1865 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( symmetric Positive de\256nite Linear Equation solution\) solves the system)9 3113(BPLE \(Banded)1 631 2 1296 1680 t (AX = B where A is a banded symmetric positive de\256nite matrix)11 2556 1 1296 1800 t 9 B f (Usage:)720 2280 w 10 R f (CALL BPLE \(N, MU, G, IG, B, IB, NB\))8 1638 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 (MU)1440 2760 w 14 S f (\256)1980 2780 w 10 R f (the number of nonzero bands on and below the diagonal of A)11 2448 1 2196 2760 t (G)1440 3000 w 14 S f (\256)1980 3020 w 10 R f ( which the upper triangular portion of the matrix A has)10 2311(a matrix into)2 533 2 2196 3000 t (been packed as follows:)3 956 1 2196 3120 t (G \( j)2 149 1 2916 3300 t 10 S f (-)3073 3300 w 10 R f (i)3136 3300 w 10 S f (+)3213 3300 w 10 R f (1 ,)1 83 1 3317 3300 t 10 I f (i)3408 3300 w 10 R f (\))3444 3300 w 10 S f (=)3534 3300 w 10 I f (a)3638 3300 w 7 I f (i j)1 45 1 3699 3320 t 10 R f (for)3793 3300 w 10 I f (j)3925 3300 w 10 S f (\263)3994 3300 w 10 I f (i)4090 3300 w 10 R f (\(See the introduction to this chapter.\))5 1487 1 2196 3480 t ( in the calling program, where)5 1348(G should be dimensioned \(IG,KG\))4 1496 2 2196 3600 t (IG)2196 3720 w 10 S f (\263)2301 3720 w 10 R f (MU and KG)2 499 1 2356 3720 t 10 S f (\263)2855 3720 w 10 R f (N.)2910 3720 w 14 S f (\254)1980 3980 w 10 R f (L and D from the factorization of A into LDL)9 1827 1 2196 3960 t 7 R f (T)4028 3920 w 10 R f (\(see)4104 3960 w 8 B f (Note 3)1 219 1 4289 3960 t 10 R f (\))4508 3960 w (IG)1440 4200 w 14 S f (\256)1980 4220 w 10 R f (the row \(leading\) dimension of G, as dimensioned in the)9 2253 1 2196 4200 t (calling program)1 635 1 2196 4320 t (B)1440 4560 w 14 S f (\256)1980 4580 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 4560 t (the calling program, where IB)4 1200 1 2196 4680 t 10 S f (\263)3396 4680 w 10 R f (N and KB)2 405 1 3451 4680 t 10 S f (\263)3856 4680 w 10 R f (NB.)3911 4680 w 14 S f (\254)1980 4940 w 10 R f (the solution X)2 567 1 2196 4920 t (IB)1440 5160 w 14 S f (\256)1980 5180 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 5160 t (calling program)1 635 1 2196 5280 t (NB)1440 5520 w 14 S f (\256)1980 5540 w 10 R f (the number of right-hand sides)4 1226 1 2196 5520 t 9 B f (Note 1:)1 278 1 720 5880 t 10 R f ( known in advance to be well-conditioned, the user should use)10 2568(Unless the given matrix A is)5 1176 2 1296 5880 t (BPSS instead of BPLE.)3 946 1 1296 6000 t 9 B f (Note 2:)1 278 1 720 6360 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 6360 t (ent right-hand sides)2 816 1 1296 6480 t 10 I f ( advance,)1 395(not all known in)3 690 2 2151 6480 t 10 R f (should call subprograms BPLE, BPFS and)5 1766 1 3274 6480 t ( is called once to get the LDL)7 1268( BPLE)1 307( example in BPFS.\))3 820( the)1 158(BPBS. \(See)1 509 5 1296 6600 t 7 R f (T)4363 6560 w 10 R f (decomposition)4451 6600 w (\(see the introduction to this chapter\) and to solve for the \256rst solution and then the pair,)16 3744 1 1296 6720 t (BPFS \(forward solve\) and BPBS \(back solve\), is called for each additional right-hand side.)13 3637 1 1296 6840 t (Linear Algebra)1 606 1 360 7500 t (\320 22 \320)2 350 1 2705 7620 t (BPLE)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 ( BPLE)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Note 3:)1 278 1 720 1320 t 10 R f (The LDL)1 378 1 1296 1320 t 7 R f (T)1679 1280 w 10 R f (decomposition of A satis\256es the equation A = LDL)8 2081 1 1759 1320 t 7 R f (T)3845 1280 w 10 R f (where L is lower unit trian-)5 1115 1 3925 1320 t ( is diagonal. On return from BPLE,)6 1428(gular \(1's on the diagonal, 0's above the diagonal\) and D)10 2316 2 1296 1440 t (the diagonal of D occupies the \256rst row of G and)10 1952 1 1296 1560 t 10 I f (G)3273 1560 w 10 R f (\()3353 1560 w 10 I f (i)3394 1560 w 10 S f (-)3430 1560 w 10 I f (j)3501 1560 w 10 S f (+)3545 1560 w 10 R f (1 ,)1 83 1 3616 1560 t 10 I f (i)3707 1560 w 10 R f (\))3743 1560 w 10 S f (=)3833 1560 w 10 I f (l)3937 1560 w 7 I f (i j)1 45 1 3976 1580 t 10 R f (for i)1 169 1 4054 1560 t 10 S f (>)4223 1560 w 10 R f (j.)4278 1560 w 9 B f (Note 4:)1 278 1 720 1920 t 10 R f ( represents the conjugate transpose of)5 1549( *)1 58( where A)2 385( *,)1 83(For complex Hermitian matrices \(A = A)6 1669 5 1296 1920 t ( decomposition and returns)3 1118( *)1 58( subroutine computes the LDL)4 1265(A\), the complex version of this)5 1303 4 1296 2040 t (the conjugate of L rather than L in G.)8 1494 1 1296 2160 t 9 B f (Error situations:)1 653 1 720 2520 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 2520 t 10 I f (Er-)4907 2520 w (ror Handling)1 531 1 1512 2640 t 10 R f (, Framework Chapter\))2 884 1 2043 2640 t (Number Error)1 1506 1 1800 2880 t (1 N)1 1008 1 1944 3120 t 10 S f (<)2977 3120 w 10 R f (1)3057 3120 w (2 MU)1 1097 1 1944 3360 t 10 S f (<)3066 3360 w 10 R f (1)3146 3360 w (3 IG)1 1041 1 1944 3600 t 10 S f (<)3010 3600 w 10 R f (MU)3090 3600 w (4 IB)1 1036 1 1944 3840 t 10 S f (<)3005 3840 w 10 R f (N)3085 3840 w (5 NB)1 1075 1 1944 4080 t 10 S f (<)3044 4080 w 10 R f (1)3124 4080 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 4320 t (10 + N + k*)4 484 1 1944 4560 t 10 I f (k)2880 4560 w 7 I f (th)2935 4520 w 10 R f (principal minor is not positive de\256nite)5 1531 1 3023 4560 t 9 B f (Double-precision version:)1 988 1 720 4920 t 10 R f (DBPLE with G and B declared double precision)7 1928 1 1783 4920 t 9 B f (Complex Hermitian version:)2 1101 1 720 5280 t 10 R f (CBPLE with G and B declared complex \(see)7 1789 1 1896 5280 t 8 B f (Note 4)1 219 1 3710 5280 t 10 R f (\).)3929 5280 w 9 B f (Storage:)720 5640 w 10 R f (None)1296 5640 w 9 B f (Time:)720 6000 w 10 R f (N)1296 6000 w 10 S f (\264)1368 6000 w 10 R f (\(\(MU)1423 6000 w 10 S f (-)1650 6000 w 10 R f (1\))1705 6000 w 10 S f (\264)1788 6000 w 10 R f (\(MU/2+2)1843 6000 w 10 S f (\264)2221 6000 w 10 R f (\(NB+1\)\)+1\) additions)1 875 1 2276 6000 t (N)1296 6120 w 10 S f (\264)1368 6120 w 10 R f (\(MU)1423 6120 w 10 S f (-)1617 6120 w 10 R f (1\))1672 6120 w 10 S f (\264)1755 6120 w 10 R f (\(MU/2+2)1810 6120 w 10 S f (\264)2188 6120 w 10 R f (NB\) multiplications)1 798 1 2243 6120 t (N)1296 6240 w 10 S f (\264)1368 6240 w 10 R f (\(MU)1423 6240 w 10 S f (-)1617 6240 w 10 R f (1+NB\) divisions)1 665 1 1672 6240 t 9 B f (Method:)720 6600 w 10 R f ( form the LDL)3 588(BPLE calls BPDC to)3 843 2 1296 6600 t 7 R f (T)2732 6560 w 10 R f (decomposition, and then calls BPFS for the forward so-)8 2231 1 2809 6600 t (lution, and BPBS for the back solution.)6 1573 1 1296 6720 t (Linear Algebra)1 606 1 4794 7500 t (\320 23 \320)2 350 1 2705 7620 t (BPLE)5155 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(BPLE February)1 4665 2 360 840 t 9 B f (See also:)1 333 1 720 1320 t 10 R f (BPBS, BPCE, BPDC, BPFS, BPLD, BPSS)5 1735 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 ( derived from the usual \256ve-point approximation to the Laplace)9 2565(The matrix in this example is)5 1179 2 1296 2040 t (operator on the unit square with an 11)7 1517 1 1296 2160 t 10 S f (\264)2813 2160 w 10 R f (11 mesh. The 100)3 716 1 2868 2160 t 10 S f (\264)3584 2160 w 10 R f (100 matrix A has the form)5 1057 1 3639 2160 t (C)2408 2280 w 10 S f (-)2636 2280 w 10 R f (I)2691 2280 w 10 S f (-)2398 2400 w 10 R f (I C)1 260 1 2453 2400 t 10 S f (-)2930 2400 w 10 R f (I)2985 2400 w 10 S f (-)2636 2520 w 10 R f (I C)1 316 1 2691 2520 t 10 S f (-)3224 2520 w 10 R f (I)3279 2520 w (........)2874 2640 w 10 S f (-)3224 2760 w 10 R f (I C)1 260 1 3279 2760 t 10 S f (-)3700 2760 w 10 R f (I)3755 2760 w 10 S f (-)3462 2880 w 10 R f (I C)1 260 1 3517 2880 t (where I is the identity matrix of order 10 and C is the matrix)13 2411 1 1296 3120 t (4)2458 3360 w 10 S f (-)2658 3360 w 10 R f (1)2713 3360 w 10 S f (-)2403 3480 w 10 R f (1 4)1 305 1 2458 3480 t 10 S f (-)2913 3480 w 10 R f (1)2968 3480 w 10 S f (-)2658 3600 w 10 R f (1 4)1 305 1 2713 3600 t 10 S f (-)3168 3600 w 10 R f (1)3223 3600 w (...)3183 3720 w 10 S f (-)3168 3840 w 10 R f (1 4)1 305 1 3223 3840 t 10 S f (-)3678 3840 w 10 R f (1)3733 3840 w 10 S f (-)3423 3960 w 10 R f (1 4)1 305 1 3478 3960 t ( side has been chosen to make)6 1227(To make it easy to detect errors in the example, the right-hand)11 2517 2 1296 4200 t ( the right-hand side the subroutine BPML is used)8 1966( construct)1 392( To)1 162(the solution a vector of all 1's.)6 1224 4 1296 4320 t ( where A is a band symmetric positive de\256nite matrix packed into the)12 2852(which produces b=Ax)2 892 2 1296 4440 t (matrix G.)1 383 1 1296 4560 t 8 CW f (INTEGER IG, N, MU, MLM1, I, KBLOK, KK, J)8 1920 1 2184 5020 t (INTEGER IWRITE, I1MACH)2 1056 1 2184 5120 t (REAL G\(11,100\), B\(100\), X\(100\))3 1440 1 2184 5220 t (REAL ERR, AMAX1)2 720 1 2184 5320 t (IG=11)2184 5420 w (N=100)2184 5520 w (MU=11)2184 5620 w (C)1656 5720 w (C SET UP MATRIX FOR ELLIPTIC PDE IN 2 DIMENSIONS)9 2304 1 1656 5820 t (C)1656 5920 w (MLM1=MU-1)2184 6020 w (I=0)2184 6120 w (DO 30 KBLOK=1,MLM1)2 864 1 2184 6220 t (DO 20 KK=1,MLM1)2 720 1 2328 6320 t (I=I+1)2472 6420 w (G\(1,I\)=4.0)2472 6520 w (G\(2,I\)=-1.0)2472 6620 w (DO 10 J=3,MLM1)2 672 1 2472 6720 t (G\(J,I\)=0.0)2616 6820 w (10 CONTINUE)1 1056 1 1800 6920 t 10 R f (Linear Algebra)1 606 1 360 7580 t (\320 24 \320)2 350 1 2705 7700 t (BPLE)360 7820 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 ( BPLE)1 4305(February 11, 1993)2 735 2 360 840 t 8 CW f (G\(MU,I\)=-1.0)2472 1300 w (20 CONTINUE)1 912 1 1800 1400 t (G\(2,I\)=0.0)2328 1500 w (30 CONTINUE)1 768 1 1800 1600 t (C)1656 1700 w (C SET UP RIGHT HAND SIDE SO SOLUTION IS ALL 1'S)10 2256 1 1656 1800 t (C)1656 1900 w (DO 40 I=1,N)2 528 1 2184 2000 t (X\(I\)=1.0)2328 2100 w (40 CONTINUE)1 816 1 1752 2200 t (CALL BPML\(N,MU,G,IG,X,B\))1 1152 1 2184 2300 t (C)1656 2400 w (C SOLVE THE SYSTEM)3 864 1 1656 2500 t (C)1656 2600 w (CALL BPLE\(N,MU,G,IG,B,100,1\))1 1344 1 2184 2700 t (C)1656 2800 w (C COMPUTE THE ERROR)3 912 1 1656 2900 t (C)1656 3000 w (ERR=0.0)2184 3100 w (DO 50 I=1,N)2 528 1 2184 3200 t (ERR=AMAX1\(ERR,ABS\(B\(I\)-1.0\)\))2328 3300 w (50 CONTINUE)1 816 1 1752 3400 t (IWRITE=I1MACH\(2\))2184 3500 w (WRITE\(IWRITE,51\)ERR)2184 3600 w ( ERROR IN SOLUTION FROM BPLE IS,F15.8\))6 1824(51 FORMAT\(31H)1 912 2 1752 3700 t (STOP)2184 3800 w (END)2184 3900 w 10 R f ( Laboratories, the fol-)3 885(When the above code was run on the Honeywell 6000 machine at Bell)12 2859 2 1296 4240 t (lowing was printed:)2 794 1 1296 4360 t ( 0.00000012)1 600(ERROR IN SOLUTION FROM BPLE IS)5 1681 2 1321 4600 t (Linear Algebra)1 606 1 4794 7500 t (\320 25 \320)2 350 1 2705 7620 t (BPLE)5155 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(BPLE February)1 4665 2 360 840 t (BPML \320 banded positive de\256nite matrix - vector multiplication)8 2583 1 1876 1320 t 9 B f (Purpose:)720 1680 w 10 R f (BPML \(Banded Positive de\256nite matrix MuLtiplication\) forms the product)8 3058 1 1296 1680 t 10 I f (Ax)4387 1680 w 10 R f (where)4525 1680 w 10 I f (A)4801 1680 w 10 R f (is a)1 145 1 4895 1680 t (symmetric banded positive matrix stored in packed form.)7 2286 1 1296 1800 t 9 B f (Usage:)720 2160 w 10 R f ( IG, X, B\))3 402( G,)1 147(CALL BPML \(N, MU,)3 925 3 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 (MU)1440 2640 w 14 S f (\256)1980 2660 w 10 R f (the number of nonzero bands on and above the diagonal of A)11 2442 1 2196 2640 t (G)1440 2880 w 14 S f (\256)1980 2900 w 10 R f ( which the upper triangular portion of the matrix A has)10 2311(a matrix into)2 533 2 2196 2880 t (been packed as follows:)3 956 1 2196 3000 t (G\(1 + j)2 289 1 2916 3180 t 10 S f (-)3205 3180 w 10 R f (i,i\) = a)2 264 1 3260 3180 t 7 R f (ij)3535 3200 w 10 R f (for j)1 169 1 3608 3180 t 10 S f (\263)3802 3180 w 10 R f (i.)3882 3180 w ( should be dimensioned)3 1043( G)1 154(\(See the introduction to this chapter.\))5 1647 3 2196 3360 t (\(IG, KG\) in the calling program, where IG)7 1698 1 2196 3480 t 10 S f (\263)3894 3480 w 10 R f (MU and KG)2 499 1 3949 3480 t 10 S f (\263)4448 3480 w 10 R f (N.)4503 3480 w (IG)1440 3720 w 14 S f (\256)1980 3740 w 10 R f ( in the calling pro-)4 778(the row \(leading\) dimension of G, as dimensioned)7 2066 2 2196 3720 t (gram)2196 3840 w (X)1440 4080 w 14 S f (\256)1980 4100 w 10 R f (the vector)1 396 1 2196 4080 t 10 I f (x)2617 4080 w 10 R f (to be multiplied)2 634 1 2686 4080 t (B)1440 4320 w 14 S f (\254)1980 4340 w 10 R f (the vector A)2 493 1 2196 4320 t 10 I f (x)2689 4320 w 9 B f (Error situations:)1 648 1 720 4680 t 10 R f (\(All errors in this subprogram are fatal \320)7 1666 1 1512 4680 t (see)1512 4800 w 10 I f (Error Handling)1 631 1 1664 4800 t 10 R f (, Framework Chapter\))2 884 1 2295 4800 t (Number Error)1 1506 1 1800 5040 t (1 N)1 1008 1 1944 5280 t 10 S f (<)2977 5280 w 10 R f (1)3057 5280 w (2 MU)1 1097 1 1944 5520 t 10 S f (<)3066 5520 w 10 R f (1)3146 5520 w (3 IG)1 1041 1 1944 5760 t 10 S f (<)3010 5760 w 10 R f (MU)3090 5760 w 9 B f (Double-precision version:)1 988 1 720 6120 t 10 R f (DBPML with G, X, and B declared double precision.)8 2128 1 1800 6120 t 9 B f (Complex Hermitian version:)2 1101 1 720 6480 t 10 R f (CBPML with G, X, and B declared complex)7 1779 1 1896 6480 t (Linear Algebra)1 606 1 360 7500 t (\320 26 \320)2 350 1 2705 7620 t (BPML)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 ( BPML)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Storage:)720 1320 w 10 R f (None)1296 1320 w 9 B f (Time:)720 1680 w 10 R f (\(2MU)1296 1680 w 10 S f (-)1540 1680 w 10 R f (1\))1595 1680 w 10 S f (\264)1678 1680 w 10 R f (N additions)1 464 1 1733 1680 t (\(2MU)1296 1800 w 10 S f (-)1540 1800 w 10 R f (1\))1595 1800 w 10 S f (\264)1678 1800 w 10 R f (N multiplications)1 698 1 1733 1800 t 9 B f (See also:)1 333 1 720 2160 t 10 R f (BPBS, BPCE, BPDC, BPLU, BPLE, BPSS)5 1745 1 1296 2160 t 9 B f (Author:)720 2520 w 10 R f (Linda Kaufman)1 629 1 1296 2520 t 9 B f (Example:)720 2880 w 10 R f ( First)1 245(This example checks the consistency of BPML and BPSS the banded system solver.)12 3499 2 1296 2880 t (the example uses BPML to compute for a given vector)9 2300 1 1296 3000 t 10 I f (x)3634 3000 w 10 R f (and a given matrix)3 782 1 3715 3000 t 10 I f (A)4534 3000 w 10 R f (the vector)1 408 1 4632 3000 t 10 I f (b)1296 3120 w 10 S f (=)1395 3120 w 10 I f (Ax.)1499 3120 w 10 R f (Then the problem is inverted, i.e., BPSS is used to \256nd the vector)12 2664 1 1683 3120 t 10 I f (x)4376 3120 w 10 R f (which satis\256es)1 591 1 4449 3120 t 10 I f (Ax)1296 3240 w 10 S f (=)1450 3240 w 10 I f (b)1554 3240 w 10 R f (. This)1 266 1 1604 3240 t 10 I f (x)1908 3240 w 10 R f ( vector)1 286( The)1 217( with the original vector.)4 1033(is then compared)2 708 4 1990 3240 t 10 I f (x)4271 3240 w 10 R f (is generated ran-)2 688 1 4352 3240 t (domly and the 10)3 697 1 1296 3360 t 10 S f (\264)1993 3360 w 10 R f (10 matrix)1 386 1 2048 3360 t 10 I f (A)2459 3360 w 10 R f (is given by)2 439 1 2545 3360 t (4)2278 3600 w 10 S f (-)2478 3600 w 10 R f (1)2533 3600 w 10 S f (-)2733 3600 w 10 R f (1)2788 3600 w 10 S f (-)2223 3720 w 10 R f (1 4)1 305 1 2278 3720 t 10 S f (-)2733 3720 w 10 R f (1)2788 3720 w 10 S f (-)2988 3720 w 10 R f (1)3043 3720 w 10 S f (-)2223 3840 w 10 R f (1)2278 3840 w 10 S f (-)2478 3840 w 10 R f (1 4)1 305 1 2533 3840 t 10 S f (-)2988 3840 w 10 R f (1)3043 3840 w 10 S f (-)3243 3840 w 10 R f (1)3298 3840 w 10 S f (-)2478 3960 w 10 R f (1)2533 3960 w 10 S f (-)2733 3960 w 10 R f (1 4)1 305 1 2788 3960 t 10 S f (-)3243 3960 w 10 R f (1)3298 3960 w 10 S f (-)3498 3960 w 10 R f (1)3553 3960 w (. . .)2 535 1 2773 4080 t 10 S f (-)2988 4200 w 10 R f (1)3043 4200 w 10 S f (-)3243 4200 w 10 R f (1 4)1 305 1 3298 4200 t 10 S f (-)3753 4200 w 10 R f (1)3808 4200 w 10 S f (-)4008 4200 w 10 R f (1)4063 4200 w 10 S f (-)3243 4320 w 10 R f (1)3298 4320 w 10 S f (-)3498 4320 w 10 R f (1 4)1 305 1 3553 4320 t 10 S f (-)4008 4320 w 10 R f (1)4063 4320 w 10 S f (-)3498 4440 w 10 R f (1)3553 4440 w 10 S f (-)3753 4440 w 10 R f (1 4)1 305 1 3808 4440 t 8 CW f (INTEGER IG, N, MU, I, IWRITE, I1MACH)6 1728 1 2088 4780 t (REAL G\(3,20\), X\(20\), B\(20\))3 1248 1 2088 4880 t (REAL UNI, ERR, COND, SASUM, ABS)5 1488 1 2088 4980 t (IG=3)2088 5080 w (N=10)2088 5180 w (MU=3)2088 5280 w (C)1656 5380 w (C CONSTRUCT MATRIX A AND PACK IT INTO G)8 1872 1 1656 5480 t (C)1656 5580 w (DO 10 I=1,N)2 528 1 2136 5680 t (G\(1,I\)=4.0)2280 5780 w (G\(2,I\)=-1.0)2280 5880 w (G\(3,I\)=-1.0)2280 5980 w (10 CONTINUE)1 768 1 1752 6080 t (C)1656 6180 w (C CONSTRUCT A RANDOM VECTOR)4 1296 1 1656 6280 t (C)1656 6380 w (DO 20 I=1,N)2 528 1 2136 6480 t (X\(I\)=UNI\(0\))2280 6580 w (20 CONTINUE)1 768 1 1752 6680 t (C)1656 6780 w (C CONSTRUCT B=AX)2 768 1 1656 6880 t 10 R f (Linear Algebra)1 606 1 4794 7540 t (\320 27 \320)2 350 1 2705 7660 t (BPML)5127 7780 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(BPML February)1 4665 2 360 840 t 8 CW f (C)1656 1300 w (CALL BPML\(N,MU,G,IG,X,B\))1 1152 1 2136 1400 t (C)1656 1500 w (C SOLVE THE SYSTEM AX=B)4 1104 1 1656 1600 t (C)1656 1700 w (CALL BPSS\(N,MU,G,IG,B,N,1,COND\))1 1488 1 2136 1800 t (C)1656 1900 w (C PRINT OUT THE TRUE SOLUTION AND THE COMPUTED SOLUTION)9 2640 1 1656 2000 t (C)1656 2100 w (IWRITE=I1MACH\(2\))2136 2200 w (WRITE\(IWRITE,21\))2136 2300 w ( SOLUTION\))1 480( COMPUTED)1 528( TRUE SOLUTION)2 672(21 FORMAT\(34H)1 864 4 1752 2400 t (WRITE\(IWRITE,22\)\(X\(I\),B\(I\),I=1,N\))2136 2500 w ( ,2E16.8\))1 432(22 FORMAT\(1H)1 816 2 1752 2600 t (ERR=0.0)2136 2700 w (DO 30 I=1,N)2 528 1 2136 2800 t (ERR=ERR+ABS\(B\(I\)-X\(I\)\))2328 2900 w (30 CONTINUE)1 768 1 1752 3000 t (ERR=ERR/SASUM\(N,X,1\))2136 3100 w (WRITE\(IWRITE,31\)ERR)2136 3200 w ( RELATIVE ERROR IS ,1PE15.7\))4 1344(31 FORMAT\(19H)1 864 2 1752 3300 t (WRITE\(IWRITE,32\)COND)2136 3400 w ( CONDITION NUMBER IS,1PE15.7\))3 1392(32 FORMAT\(20H)1 864 2 1752 3500 t (STOP)2136 3600 w (END)2136 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 (which has a machine precision of 1.)6 1433 1 1296 4160 t 10 S f (\264)2737 4160 w 10 R f (10)2800 4160 w 7 S f (-)2905 4120 w 7 R f (8)2949 4120 w 10 R f (, the following was printed:)4 1102 1 2992 4160 t 8 CW f ( SOLUTION)1 432( COMPUTED)1 528(TRUE SOLUTION)1 624 3 1704 4500 t ( 00)1 144( 0.22925608E)1 624(0.22925607E 00)1 672 3 1800 4600 t ( 00)1 144( 0.76687504E)1 624(0.76687502E 00)1 672 3 1800 4700 t ( 00)1 144( 0.68317687E)1 624(0.68317685E 00)1 672 3 1800 4800 t ( 00)1 144( 0.50919112E)1 624(0.50919111E 00)1 672 3 1800 4900 t ( 00)1 144( 0.87455962E)1 624(0.87455959E 00)1 672 3 1800 5000 t ( 00)1 144( 0.64464103E)1 624(0.64464101E 00)1 672 3 1800 5100 t ( 00)1 144( 0.84746842E)1 624(0.84746840E 00)1 672 3 1800 5200 t ( 00)1 144( 0.35396345E)1 624(0.35396343E 00)1 672 3 1800 5300 t ( 00)1 144( 0.39889160E)1 624(0.39889160E 00)1 672 3 1800 5400 t ( 00)1 144( 0.45709422E)1 624(0.45709422E 00)1 672 3 1800 5500 t ( 3.1985797E)1 624(RELATIVE ERROR IS)2 816 2 1704 5600 t 8 S f (-)3144 5600 w 8 CW f (08)3188 5600 w ( 01)1 144( 2.1901962E)1 576(CONDITION NUMBER IS)2 912 3 1704 5700 t 10 R f ( number of the matrix and the precision of the Honeywell suggest that even in)14 3182(The condition)1 562 2 1296 6060 t (the absence of roundoff error in BPML, a relative error of 2.2)11 2490 1 1296 6180 t 10 S f (\264)3786 6180 w 10 R f (10)3841 6180 w 7 S f (-)3946 6140 w 7 R f (7)3990 6140 w 10 R f ( be surprising.)2 573(would not)1 406 2 4061 6180 t (The value computed above is quite reasonable.)6 1871 1 1296 6300 t (Linear Algebra)1 606 1 360 7500 t (\320 28 \320)2 350 1 2705 7620 t (BPML)360 7740 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 ( BPML)1 4305(February 11, 1993)2 735 2 360 840 t ( of a banded symmetric positive de\256nite matrix)7 1890( norm)1 261(BPNM \320)1 409 3 1888 1320 t 9 B f (Purpose:)720 1680 w 10 R f (BPNM \(Banded Positive de\256nite matrix NorM\) computes the norm of a banded symmetric)12 3744 1 1296 1680 t (positive de\256nite matrix A stored in packed form. The in\256nity norm is de\256ned as)13 3744 1 1296 1800 t 7 R f (1)1296 2040 w 7 S f (\243)1342 2040 w 7 I f (i)1397 2040 w 7 S f (\243)1433 2040 w 7 I f (n)1483 2040 w 10 R f (max)1321 1970 w 7 I f (j)1528 2070 w 7 S f (=)1559 2070 w 7 R f (1)1609 2070 w 15 S f (S)1542 2000 w 7 I f (n)1569 1870 w 10 S f (\357)1655 1970 w 10 I f (a)1712 1970 w 7 I f (i j)1 45 1 1773 1990 t 10 S f (\357)1858 1970 w 9 B f (Type:)720 2410 w 10 R f (Real function)1 541 1 1296 2410 t 9 B f (Usage:)720 2770 w 10 S f (<)1296 2770 w 10 R f (answer)1351 2770 w 10 S f (>)1633 2770 w 10 R f (= BPNM \(N, MU, G, IG\))5 1016 1 1713 2770 t (N)1440 3010 w 14 S f (\256)1980 3030 w 10 R f (the number of rows in A)5 979 1 2196 3010 t (MU)1440 3250 w 14 S f (\256)1980 3270 w 10 R f (the number of nonzero bands in A on and above the diagonal)11 2437 1 2196 3250 t (G)1440 3490 w 14 S f (\256)1980 3510 w 10 R f ( which the upper triangular portion of the matrix A has)10 2311(a matrix into)2 533 2 2196 3490 t (been packed as follows:)3 956 1 2196 3610 t (G \(1 + j)3 314 1 2916 3790 t 10 S f (-)3230 3790 w 10 R f (i, i\) =)2 220 1 3285 3790 t 10 I f (a)3530 3790 w 7 I f (i j)1 45 1 3591 3810 t 10 R f (for j)1 169 1 3669 3790 t 10 S f (\263)3838 3790 w 10 R f (i)3893 3790 w (i.e. the main diagonal of A is in the \256rst row of G)12 1976 1 2196 3970 t (\(See the introduction to this chapter.\))5 1487 1 2196 4090 t ( in the calling program, where)5 1348(G should be dimensioned \(IG,KG\))4 1496 2 2196 4210 t (IG)2196 4330 w 10 S f (\263)2301 4330 w 10 R f (MU and KG)2 499 1 2356 4330 t 10 S f (\263)2855 4330 w 10 R f (N.)2910 4330 w (IG)1440 4570 w 14 S f (\256)1980 4590 w 10 R f (the row \(leading\) dimension of G, as dimensioned in the)9 2253 1 2196 4570 t (calling program)1 635 1 2196 4690 t 10 S f (<)1440 4980 w 10 R f (answer)1495 4980 w 10 S f (>)1777 4980 w 14 S f (\254)1980 5000 w 7 R f (1)2196 5050 w 7 S f (\243)2242 5050 w 7 I f (i)2292 5050 w 7 S f (\243)2317 5050 w 7 I f (n)2361 5050 w 10 R f (max)2210 4980 w 7 I f (j)2406 5080 w 7 S f (=)2437 5080 w 7 R f (1)2487 5080 w 15 S f (S)2420 5010 w 7 I f (n)2447 4880 w 10 S f (\357)2533 4980 w 10 I f (a)2590 4980 w 7 I f (i j)1 45 1 2651 5000 t 10 S f (\357)2736 4980 w 9 B f (Error situations:)1 648 1 720 5420 t 10 R f (\(All errors in this subprogram are fatal \320)7 1666 1 1512 5420 t (see)1512 5540 w 10 I f (Error Handling)1 631 1 1664 5540 t 10 R f (, Framework Chapter\))2 884 1 2295 5540 t (Number Error)1 1506 1 1800 5780 t (1 N)1 1008 1 1944 6020 t 10 S f (<)2977 6020 w 10 R f (1)3057 6020 w (2 MU)1 1097 1 1944 6260 t 10 S f (<)3066 6260 w 10 R f (1)3146 6260 w (3 IG)1 1041 1 1944 6500 t 10 S f (<)3010 6500 w 10 R f (MU)3090 6500 w (Linear Algebra)1 606 1 4794 7500 t (\320 29 \320)2 350 1 2705 7620 t (BPNM)5116 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 ( 28, 1979)2 375(BPNM March)1 4665 2 360 840 t 9 B f (Double-precision version:)1 988 1 720 1320 t 10 R f (DBPNM with G and DBPNM declared double precision)7 2256 1 1783 1320 t 9 B f (Complex version:)1 678 1 720 1680 t 10 R f (CBPNM with G declared complex)4 1382 1 1448 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 (N)1296 2400 w 10 S f (\264)1393 2400 w 10 R f (\(2)1473 2400 w 10 S f (\264)1581 2400 w 10 R f (MU)1661 2400 w 10 S f (-)1847 2400 w 10 R f (1\) additions)1 475 1 1927 2400 t 9 B f (See also:)1 333 1 720 2760 t 10 R f (BPDC, BPLD, BPLE, BPSS, BPCE)4 1449 1 1296 2760 t 9 B f (Author:)720 3120 w 10 R f (Linda Kaufman)1 629 1 1296 3120 t 9 B f (Example:)720 3480 w 10 R f ( a matrix is singular. One)5 1020(Because of roundoff error it is often very dif\256cult to decide whether)11 2724 2 1296 3480 t (criteria often used for symmetric positive de\256nite matrices is to compute the LDL)12 3268 1 1296 3600 t 7 I f (T)4569 3560 w 10 R f (decompo-)4641 3600 w ( element of the diagonal matrix D is)7 1454(sition of the matrix and declare the matrix singular if any)10 2290 2 1296 3720 t (less than)1 347 1 1296 3840 t 10 S f (e \357 \357)2 158 1 1668 3840 t 10 I f (A)1834 3840 w 10 S f (\357 \357)1 106 1 1903 3840 t 10 R f (where)2034 3840 w 10 S f (e)2302 3840 w 10 R f (is the machine precision and is computed by R1MACH\(4\).)8 2356 1 2371 3840 t ( used to indicate whether the symmetric positive)7 2004(The following program fragment might be)5 1740 2 1296 4080 t ( uses the fact that the)5 875( It)1 118( G is singular or nearly singular.)6 1332(de\256nite banded matrix packed into)4 1419 4 1296 4200 t (subroutine BPLD, which computes the LDL)5 1796 1 1296 4320 t 7 I f (T)3097 4280 w 10 R f ( a banded symmetric matrix,)4 1164(decomposition of)1 702 2 3174 4320 t ( of D less than EPS, an input parameter)8 1605(issues a recoverable error when it detects an element)8 2139 2 1296 4440 t (to the subroutine.)2 697 1 1296 4560 t 8 CW f (IWRITE=I1MACH\(2\))2136 4900 w (CALL ENTSRC\(IROLD,1\))1 960 1 2136 5000 t (EPS=BPNM\(N,MU,G,IG\)*R1MACH\(4\))2136 5100 w (CALL BPLD\(N,MU,G,IG,EPS\))1 1152 1 2136 5200 t (IF \(NERROR\(IERR\).EQ.0\) GO TO 10)4 1488 1 2136 5300 t (CALL ERROFF)1 528 1 2280 5400 t (WRITE\(IWRITE,1\))2280 5500 w ( SINGULAR MATRIX\))2 816(1 FORMAT\(16H)1 1056 2 1704 5600 t (10 CONTINUE)1 816 1 1704 5700 t 10 R f (Linear Algebra)1 606 1 360 7500 t (\320 30 \320)2 350 1 2705 7620 t (BPNM)360 7740 w cleartomark showpage saveobj restore %%EndPage: 30 30 %%Page: 31 31 /saveobj save def mark 31 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 ( BPNM)1 4405(March 28, 1979)2 635 2 360 840 t (BPSS \320 band positive de\256nite linear system solution with condition estimation)10 3197 1 1569 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( where A is a band)5 760( = B)2 173( solves the system AX)4 906(BPSS \(Band Positive de\256nite System Solution\))5 1905 4 1296 1680 t ( estimate of the condition number of A is provided.)9 2048( An)1 172(symmetric positive de\256nite matrix.)3 1400 3 1296 1800 t 9 B f (Usage:)720 2160 w 10 R f (CALL BPSS \(N, MU, G, IG, B, IB, NB, COND\))9 1961 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 (MU)1440 2640 w 14 S f (\256)1980 2660 w 10 R f (the number of nonzero bands on and above the diagonal of A)11 2442 1 2196 2640 t (G)1440 2880 w 14 S f (\256)1980 2900 w 10 R f ( which the upper triangular portion of the matrix A has)10 2311(a matrix into)2 533 2 2196 2880 t (been placed as follows:)3 934 1 2196 3000 t (G\(j)2916 3180 w 10 S f (-)3049 3180 w 10 R f (i+1, i\) =)2 326 1 3104 3180 t 10 I f (a)3455 3180 w 7 I f (i j)1 45 1 3516 3200 t 10 R f (for j)1 169 1 3594 3180 t 10 S f (\263)3763 3180 w 10 R f (i)3818 3180 w (\(See the introduction to this chapter.\))5 1487 1 2196 3360 t ( in the calling program, where)5 1348(G should be dimensioned \(IG,KG\))4 1496 2 2196 3480 t (IG)2196 3600 w 10 S f (\263)2301 3600 w 10 R f (MU and KG)2 499 1 2356 3600 t 10 S f (\263)2855 3600 w 10 R f (N.)2910 3600 w 14 S f (\254)1980 3860 w 10 R f (L and D from the factorization of A into LDL)9 1827 1 2196 3840 t 7 R f (T)4028 3800 w 10 R f (\(see)2196 3960 w 8 B f (Note 3)1 219 1 2381 3960 t 10 R f (\))2600 3960 w (IG)1440 4200 w 14 S f (\256)1980 4220 w 10 R f (the row \(leading\) dimension of G, as dimensioned in the)9 2253 1 2196 4200 t (calling program)1 635 1 2196 4320 t (B)1440 4560 w 14 S f (\256)1980 4580 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 4560 t (the calling program, where IB)4 1200 1 2196 4680 t 10 S f (\263)3396 4680 w 10 R f (N and KB)2 405 1 3451 4680 t 10 S f (\263)3856 4680 w 10 R f (NB.)3911 4680 w 14 S f (\254)1980 4940 w 10 R f (the solution X)2 567 1 2196 4920 t (IB)1440 5160 w 14 S f (\256)1980 5180 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 5160 t (calling program)1 635 1 2196 5280 t (NB)1440 5520 w 14 S f (\256)1980 5540 w 10 R f (the number of right-hand sides)4 1226 1 2196 5520 t (COND)1440 5760 w 14 S f (\254)1980 5780 w 10 R f ( \(See)1 227(an estimate of the condition number of A)7 1645 2 2196 5760 t 8 B f (Note 1)1 219 1 4093 5760 t 10 R f (\))4312 5760 w 9 B f (Note 1:)1 278 1 720 6120 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 6120 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 6240 t (of your linear system have)4 1091 1 1296 6360 t 10 B f (d)2420 6360 w 10 R f ( solution might have as few as)6 1264(decimal digits of precision, the)4 1267 2 2509 6360 t 10 B f (d)1296 6480 w 10 S f (-)1377 6480 w 10 R f (log)1457 6480 w 7 R f (10)1596 6500 w 10 R f ( if COND is greater than 10)6 1121( Thus)1 252(\(COND\) correct decimal digits.)3 1270 3 1674 6480 t 7 R f (Bd)4327 6440 w 7 I f (P)4419 6440 w 10 R f ( be)1 120(, there may)2 450 2 4470 6480 t (no correct digits.)2 674 1 1296 6600 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 6840 t (Linear Algebra)1 606 1 4794 7500 t (\320 31 \320)2 350 1 2705 7620 t (BPSS)5165 7740 w cleartomark showpage saveobj restore %%EndPage: 31 31 %%Page: 32 32 /saveobj save def mark 32 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 ( 28, 1979)2 375(BPSS March)1 4665 2 360 840 t ( however, the user is)4 871( Ordinarily,)1 503(use the routine BPLE, which is a little faster than BPSS.)10 2370 3 1296 1320 t (strongly urged to choose BPSS, and to follow it by a test of the condition estimate.)15 3308 1 1296 1440 t 9 B f (Note 2:)1 278 1 720 1800 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 1800 t (ent right-hand sides)2 798 1 1296 1920 t 10 I f (not all known in advance,)4 1050 1 2124 1920 t 10 R f ( use BPSS, but should call subpro-)6 1411(should not)1 425 2 3204 1920 t ( is called once to get the)6 983( BPCE)1 304( the example of BPFS.\))4 948( \(See)1 230( BPBS.)1 299(grams BPCE, BPFS and)3 980 6 1296 2040 t (LDL)1296 2160 w 7 R f (T)1495 2120 w 10 R f ( and then the pair, BPFS \(forward)6 1375(decomposition \(see the introduction to this chapter\))6 2089 2 1576 2160 t (solve\) and BPBS \(back solve\), is called for each new right-hand side.)11 2770 1 1296 2280 t 9 B f (Note 3:)1 278 1 720 2640 t 10 R f (The LDL)1 378 1 1296 2640 t 7 R f (T)1679 2600 w 10 R f (decomposition of A satis\256es the equation A = LDL)8 2081 1 1759 2640 t 7 R f (T)3845 2600 w 10 R f (where L is lower unit trian-)5 1115 1 3925 2640 t ( 0's above the diagonal\) and D is diagonal. On return from BPSS,)12 2673(gular \(1's on the diagonal,)4 1071 2 1296 2760 t (the diagonal of D occupies the \256rst row of G and G\(i)11 2110 1 1296 2880 t 10 S f (-)3406 2880 w 10 R f (j+1,i\))3461 2880 w 10 S f (=)3722 2880 w 10 I f (l)3826 2880 w 7 I f (i j)1 45 1 3865 2900 t 10 R f (for i)1 169 1 3943 2880 t 10 S f (>)4112 2880 w 10 R f (j.)4167 2880 w 9 B f (Note 4:)1 278 1 720 3240 t 10 R f ( represents the conjugate transpose of)5 1549( *)1 58( where A)2 385( *,)1 83(For complex Hermitian matrices \(A = A)6 1669 5 1296 3240 t ( decomposition and returns)3 1118( *)1 58( subroutine computes the LDL)4 1265(A\), the complex version of this)5 1303 4 1296 3360 t (the conjugate of L rather than L in G.)8 1494 1 1296 3480 t 9 B f (Error situations:)1 653 1 720 3840 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 3840 t 10 I f (Er-)4907 3840 w (ror Handling)1 531 1 1512 3960 t 10 R f (, Framework Chapter\))2 884 1 2043 3960 t (Number Error)1 1506 1 1800 4200 t (1 N)1 1008 1 1944 4440 t 10 S f (<)2977 4440 w 10 R f (1)3057 4440 w (2 MU)1 1097 1 1944 4680 t 10 S f (<)3066 4680 w 10 R f (1)3146 4680 w (3 IG)1 1041 1 1944 4920 t 10 S f (<)3010 4920 w 10 R f (MU)3090 4920 w (4 IB)1 1036 1 1944 5160 t 10 S f (<)3005 5160 w 10 R f (N)3085 5160 w (5 NB)1 1075 1 1944 5400 t 10 S f (<)3044 5400 w 10 R f (1)3124 5400 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 5640 t (10 + N + k*)4 484 1 1944 5880 t 10 I f (k)2880 5880 w 7 I f (th)2935 5840 w 10 R f (principal minor is not positive de\256nite)5 1531 1 3023 5880 t 9 B f (Double-precision version:)1 988 1 720 6240 t 10 R f (DBPSS, with G, B, and COND declared double precision.)8 2326 1 1783 6240 t 9 B f (Complex Hermitian version:)2 1101 1 720 6600 t 10 R f (CBPSS with G and B declared complex \(see)7 1779 1 1896 6600 t 8 B f (Note 4)1 219 1 3700 6600 t 10 R f (above\))3944 6600 w (Linear Algebra)1 606 1 360 7500 t (\320 32 \320)2 350 1 2705 7620 t (BPSS)360 7740 w cleartomark showpage saveobj restore %%EndPage: 32 32 %%Page: 33 33 /saveobj save def mark 33 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 ( BPSS)1 4405(March 28, 1979)2 635 2 360 840 t 9 B f (Storage:)720 1320 w 10 R f ( storage in the)3 575(N real \(double precision for DBPSS, complex for CBPSS\) locations of scratch)11 3169 2 1296 1320 t (dynamic storage stack)2 887 1 1296 1440 t 9 B f (Time:)720 1800 w 10 R f (N)1296 1800 w 10 S f (\264)1368 1800 w 10 R f (\(\(MU)1423 1800 w 10 S f (-)1650 1800 w 10 R f (1\))1705 1800 w 10 S f (\264)1788 1800 w 10 R f (\(MU/2+2)1843 1800 w 10 S f (\264)2221 1800 w 10 R f (NB+8\)+8\) additions)1 809 1 2276 1800 t (N)1296 1920 w 10 S f (\264)1368 1920 w 10 R f (\(\(MU)1423 1920 w 10 S f (-)1650 1920 w 10 R f (1\))1705 1920 w 10 S f (\264)1788 1920 w 10 R f (\(MU/2+2)1843 1920 w 10 S f (\264)2221 1920 w 10 R f (NB+6\)+4\) multiplications)1 1043 1 2276 1920 t (N)1296 2040 w 10 S f (\264)1368 2040 w 10 R f (\(MU+1+NB\) divisions)1 915 1 1423 2040 t 9 B f (Method:)720 2400 w 10 R f (BPSS calls BPCE to form the LDL)6 1419 1 1296 2400 t 7 R f (T)2720 2360 w 10 R f ( the forward so-)3 643(decomposition, and then calls BPFS for)5 1599 2 2798 2400 t ( to esti-)2 304( the reference below for the method used)7 1652( See)1 197(lution, and BPBS for the back solution.)6 1591 4 1296 2520 t (mate the condition number.)3 1099 1 1296 2640 t 9 B f (See also:)1 333 1 720 3000 t 10 R f (BPBS, BPCE, BPDC, BPFS, BPLD, BPLE)5 1745 1 1296 3000 t 9 B f (Author:)720 3360 w 10 R f (Linda Kaufman)1 629 1 1296 3360 t 9 B f (Reference:)720 3720 w 10 R f ( W., and Wilkinson, J. H., An estimate for the condi-)10 2154(Cline, A. K., Moler, C. B., Stewart, G.)7 1562 2 1324 3720 t (tion number,)1 511 1 1296 3840 t 10 I f (SIAM J. Numer. Anal. 16)4 1007 1 1857 3840 t 10 R f (\(1979\), 368-375.)1 674 1 2889 3840 t 9 B f (Example:)720 4200 w 10 R f ( a discretization of a)4 829(In the following example we solve three problems that might arise from)11 2915 2 1296 4200 t (1-dimensional differential equation. The coef\256cient matrix in each problem has the form)11 3538 1 1296 4320 t (1+x)2097 4560 w 10 S f (-)2403 4560 w 10 R f (1)2458 4560 w 10 S f (-)2042 4680 w 10 R f (1 2)1 411 1 2097 4680 t 10 S f (-)2658 4680 w 10 R f (1)2713 4680 w 10 S f (-)2403 4800 w 10 R f (1 2)1 305 1 2458 4800 t 10 S f (-)2913 4800 w 10 R f (1)2968 4800 w 10 S f (-)2658 4920 w 10 R f (1 2)1 305 1 2713 4920 t 10 S f (-)3168 4920 w 10 R f (1)3223 4920 w (.)2953 5040 w (.)3208 5160 w 10 S f (-)3423 5280 w 10 R f (1 2)1 305 1 3478 5280 t 10 S f (-)3933 5280 w 10 R f (1)3988 5280 w 10 S f (-)3678 5400 w 10 R f (1 1+x)1 411 1 3733 5400 t ( matrix is singular when x is 0,)7 1282( The)1 212( the three problems.)3 814(with x=1.0, .01, and .0001 de\256ning)5 1436 4 1296 5640 t ( To)1 166( our output suggests, the matrix becomes more ill-conditioned as x approaches 0.)12 3295(and, as)1 283 3 1296 5760 t ( solution, the right-hand side has been chosen to make the)10 2351(make it easy to detect errors in the)7 1393 2 1296 5880 t ( should notice in the output the correlation between)8 2125( reader)1 283( The)1 215(solution a vector of all 1's.)5 1121 4 1296 6000 t ( with most problems with nearly constant diagonal, it)8 2170(the condition number and the error. As)6 1574 2 1296 6120 t (was very easy to write the code to set up the problem for BPSS.)13 2547 1 1296 6240 t 8 CW f (INTEGER N, K, I, IWRITE, I1MACH, MU)6 1680 1 1776 6580 t (REAL G\(2,100\), B\(200\))2 1008 1 1776 6680 t (REAL X, COND, ERR, AMAX1)4 1152 1 1776 6780 t (C CONSTRUCT MATRIX AND RIGHT-HAND SIDE SO TRUE SOLUTION IS)9 2784 1 1440 6880 t 10 R f (Linear Algebra)1 606 1 4794 7540 t (\320 33 \320)2 350 1 2705 7660 t (BPSS)5165 7780 w cleartomark showpage saveobj restore %%EndPage: 33 33 %%Page: 34 34 /saveobj save def mark 34 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 ( 28, 1979)2 375(BPSS March)1 4665 2 360 840 t 8 CW f (C COMPOSED ENTIRELY OF ONES)4 1296 1 1440 1300 t (N=100)1776 1400 w (X=1)1776 1500 w ( K=1,3)1 336(DO 30)1 240 2 1776 1600 t (DO 10 I=1,N)2 528 1 1920 1700 t (G\(1,I\)=2.0)2064 1800 w (G\(2,I\)=-1.0)2064 1900 w (B\(I\)=0.0)2064 2000 w (10 CONTINUE)1 768 1 1536 2100 t (G\(1,1\)=1.0+X)1920 2200 w (G\(1,N\)=1.0+X)1920 2300 w (B\(1\)=X)1920 2400 w (B\(N\)=X)1920 2500 w (C SOLVE THE SYSTEM)3 864 1 1440 2600 t (MU=2)1920 2700 w (CALL BPSS\(N,MU,G,2,B,N,1,COND\))1 1440 1 1920 2800 t (IWRITE=I1MACH\(2\))1920 2900 w (WRITE\(IWRITE,11\)X)1920 3000 w ( X IS,F15.7\))2 576(11 FORMAT\(/5H)1 864 2 1536 3100 t (WRITE\(IWRITE,12\)COND)1920 3200 w ( CONDITION NUMBER IS,1PE15.7\))3 1392(12 FORMAT\(20H)1 864 2 1536 3300 t (C COMPUTE THE ERROR)3 912 1 1440 3400 t (ERR=0.0)1920 3500 w (DO 20 I=1,N)2 528 1 1920 3600 t (ERR=AMAX1\(ERR,ABS\(B\(I\)-1.0\)\))2064 3700 w (20 CONTINUE)1 768 1 1536 3800 t (WRITE\(IWRITE,21\)ERR)1920 3900 w ( FOR BPSS THE ERROR IS,F16.8\))5 1392(21 FORMAT\(22H)1 864 2 1536 4000 t (X=X/100.)1920 4100 w (30 CONTINUE)1 624 1 1536 4200 t (STOP)1776 4300 w (END)1776 4400 w 10 R f ( on the Honeywell 6000 machine at Bell Laboratories,)8 2174(When the above program was executed)5 1570 2 1296 4740 t (the following was printed)3 1024 1 1296 4860 t 8 CW f ( 1.0000000)1 720(X IS)1 192 2 1704 5280 t ( 03)1 144( 4.0807862E)1 576(CONDITION NUMBER IS)2 912 3 1704 5380 t ( 0.00000431)1 768(FOR BPSS THE ERROR IS)4 1008 2 1704 5480 t ( 0.0100000)1 720(X IS)1 192 2 1704 5680 t ( 04)1 144( 2.3329148E)1 576(CONDITION NUMBER IS)2 912 3 1704 5780 t ( 0.00002055)1 768(FOR BPSS THE ERROR IS)4 1008 2 1704 5880 t ( 0.0001000)1 720(X IS)1 192 2 1704 6080 t ( 06)1 144( 1.9933923E)1 576(CONDITION NUMBER IS)2 912 3 1704 6180 t ( 0.00491761)1 768(FOR BPSS THE ERROR IS)4 1008 2 1704 6280 t 10 R f (Linear Algebra)1 606 1 360 7500 t (\320 34 \320)2 350 1 2705 7620 t (BPSS)360 7740 w (-- --)1 5544 1 0 8030 t cleartomark showpage saveobj restore %%EndPage: 34 34 %%Trailer done %%Pages: 34 %%DocumentFonts: Courier Times-Bold Times-Italic Times-Roman Symbol