%!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