%!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: 0 1
/saveobj save def
mark
1 pagesetup
12 B f
(Programs for Solving Linear Equations in the PORT Library)8 3151 1 1304 1230 t
10 I f
(Linda Kaufman)1 626 1 2567 1470 t
10 R f
(AT&T Bell Laboratories)2 993 1 2383 1650 t
(Murray Hill, New Jersey 07974)4 1267 1 2246 1770 t
10 I f
(ABSTRACT)2643 2270 w
10 R f
(This paper describes the subroutines that have recently been inserted into the PORT)12 3350 1 1330 2566 t
( of the subroutines are high-level drivers which)7 1937( Some)1 285(library for solving linear systems.)4 1378 3 1080 2686 t
(solve)1080 2806 w
10 I f
(AX)2532 2986 w
10 S f
(=)2703 2986 w
10 I f
(B)2807 2986 w
10 R f
( are)1 157( Others)1 327( the solution to perturbations in the problems.)7 1899(and indicate the sensitivity of)4 1217 4 1080 3166 t
( for complicated problems such as solving a sequence of)9 2335(low level subroutines designed)3 1265 2 1080 3286 t
( depend on pre-)3 638(problems with the same matrix but with different right-hand sides, which)10 2962 2 1080 3406 t
( subroutines are classified on the basis of the structure of the)11 2624( The)1 224(vious solutions.)1 648 3 1080 3526 t
10 I f
(A)4619 3526 w
10 R f
( it is symmetric, banded, sparse, etc.)6 1448( whether)1 371(matrix, e.g.)1 455 3 1080 3646 t
(February 11, 1993)2 735 1 720 4126 t
cleartomark
showpage
saveobj restore
%%EndPage: 0 1
%%Page: 1 2
/saveobj save def
mark
2 pagesetup
12 B f
(Programs for Solving Linear Equations in the PORT Library)8 3151 1 1304 1230 t
10 I f
(Linda Kaufman)1 626 1 2567 1470 t
10 R f
(AT&T Bell Laboratories)2 993 1 2383 1650 t
(Murray Hill, New Jersey 07974)4 1267 1 2246 1770 t
10 B f
(1. Introduction)1 670 1 720 2130 t
10 R f
( many routines for solving linear equations)6 1775(The linear algebra chapter in the PORT library contains)8 2295 2 970 2286 t
( are high-level drivers \(see section 2.1\) designed for those users)10 2641( Some)1 288( least squares problems.)3 984(and linear)1 407 4 720 2406 t
(who simply want to solve a system of equations)8 1916 1 720 2526 t
10 I f
(AX)2712 2706 w
10 S f
(=)2883 2706 w
10 I f
(B)2987 2706 w
10 R f
(Others are lower-level routines \(see section 2.2\) designed for users with slightly more complicated problem.)14 4320 1 720 2886 t
( same matrix but different right-hand)5 1528(e.g., users who wish to solve a sequence of linear systems with the)12 2792 2 720 3006 t
(sides.)720 3126 w
( the following categories, with the first two letters in the name)11 2552(Presently, the package is divided into)5 1518 2 970 3282 t
(of each subroutine specifying the type of matrix.)7 1941 1 720 3402 t
10 B f
(Table 1)1 320 1 2720 3642 t
10 R f
( Meaning)1 1407( letter prefix)2 493(Category 2)1 923 3 1195 3942 t
10 S f
(_ _________________________________________________________________________________)1 4083 1 838 3952 t
(_ _________________________________________________________________________________)1 4083 1 838 3972 t
10 R f
( known specific properties)3 1056( no)1 455(General GE)1 1568 3 838 4082 t
10 S f
(_ _________________________________________________________________________________)1 4083 1 838 4102 t
10 R f
(Symmetric SY)1 1565 1 838 4222 t
10 I f
(A)2761 4222 w
10 R f
(\()2830 4222 w
10 I f
(i)2871 4222 w
10 R f
(,)2907 4222 w
10 I f
(j)2948 4222 w
10 R f
(\))2984 4222 w
10 S f
(=)3074 4222 w
10 I f
(A)3178 4222 w
10 R f
(\()3247 4222 w
10 I f
(j)3304 4222 w
10 R f
(,)3340 4222 w
10 I f
(i)3373 4222 w
10 R f
(\))3409 4222 w
10 S f
(_ _________________________________________________________________________________)1 4083 1 838 4242 t
10 R f
( non-zero elements are near the main diagonal)7 1842( all)1 452(Banded BA)1 1571 3 838 4362 t
10 S f
(_ _________________________________________________________________________________)1 4083 1 838 4382 t
10 R f
( posi-)1 252( symmetric,)1 498(BP Banded,)1 -1110 3 2278 4502 t
(tive definite)1 480 1 838 4622 t
( and all eigenvalues)3 811(banded and symmetric, as above,)4 1349 2 2761 4502 t
(positive)2761 4622 w
10 S f
(_ _________________________________________________________________________________)1 4083 1 838 4642 t
10 R f
( known)1 321( placements of all nonzero elements are)6 1717( the)1 488(Sparse SP)1 1557 4 838 4762 t
( 1% of the ma-)4 625(and they occupy not more than about)6 1535 2 2761 4882 t
(trix)2761 5002 w
10 S f
(\347)1993 5002 w
(\347)1993 4922 w
(\347)1993 4822 w
(\347)1993 4722 w
(\347)1993 4622 w
(\347)1993 4522 w
(\347)1993 4422 w
(\347)1993 4322 w
(\347)1993 4222 w
(\347)1993 4122 w
(\347)1993 4022 w
(\347)1993 3922 w
(\347)2686 5002 w
(\347)2686 4922 w
(\347)2686 4822 w
(\347)2686 4722 w
(\347)2686 4622 w
(\347)2686 4522 w
(\347)2686 4422 w
(\347)2686 4322 w
(\347)2686 4222 w
(\347)2686 4122 w
(\347)2686 4022 w
(\347)2686 3922 w
10 R f
( example, for symmetric matrices, the)5 1541( For)1 195( each category is a distinct data structure.)7 1688(Associated with)1 646 4 970 5338 t
( structures are discussed in)4 1090( section 3 the data)4 745( In)1 139(upper triangular portion only is stored in one long vector.)9 2346 4 720 5458 t
(detail.)720 5578 w
( implemented a subroutine for computing)5 1707(Note that we have not)4 912 2 970 5734 t
10 I f
(A)3624 5734 w
7 S f
(-)3696 5694 w
7 R f
(1)3746 5694 w
10 R f
( necessary, the matrix)3 900(. If)1 151 2 3789 5734 t
10 I f
(A)4875 5734 w
7 S f
(-)4947 5694 w
7 R f
(1)4997 5694 w
10 R f
(can be found by solving)4 960 1 720 5854 t
10 I f
(AX)1705 5854 w
10 S f
(=)1876 5854 w
10 I f
(I)1980 5854 w
10 R f
(.)2013 5854 w
( this package have been implemented in single precision, double precision, and)11 3302(The subroutines in)2 768 2 970 6010 t
(complex arithmetic.)1 799 1 720 6130 t
10 B f
( of the Package)3 646(2. Structure)1 535 2 720 6370 t
( Drivers)1 346(2.1. Higher-Level)1 766 2 720 6610 t
10 R f
( high-level subroutines \320 the Linear Equation)6 1967(Associated with each category in Table 1 are two)8 2103 2 970 6766 t
(solver \(LE\) and the System Solver \(SS\) which call lower-level subroutines.)10 3010 1 720 6886 t
( LE subroutines \(GELE, SYLE, BALE,)5 1617(Both the SS subroutines \(GESS, SYSS, BASS, etc.\) and the)9 2453 2 970 7042 t
( systems)1 342(etc.\) solve)1 435 2 720 7162 t
cleartomark
showpage
saveobj restore
%%EndPage: 1 2
%%Page: 2 3
/saveobj save def
mark
3 pagesetup
10 R f
(- 2 -)2 166 1 2797 480 t
10 I f
(AX)2712 840 w
10 S f
(=)2883 840 w
10 I f
(B)2987 840 w
10 R f
(\(2.1\))4849 840 w
( example, partial pivot-)3 938( For)1 191( matrix[8].)1 429(using variants of Gaussian elimination adapted to the structure of the)10 2762 4 720 1020 t
( The)1 208( for band symmetric positive definite matrices.)6 1885(ing is used for general matrices but no pivoting is used)10 2227 3 720 1140 t
(matrix)720 1260 w
10 I f
(B)1018 1260 w
10 R f
( all the subroutines the matrix)5 1251( In)1 145(in \(2.1\) may have several columns.)5 1460 3 1116 1260 t
10 I f
(A)4009 1260 w
10 R f
(is overwritten, and the)3 932 1 4108 1260 t
(solution replaces the right-hand side matrix B.)6 1851 1 720 1380 t
( the SS subroutines returns an estimate of the condition number)10 2787(Besides solving \(2.2\), each of)4 1283 2 970 1536 t
(\(COND\) of the matrix)3 914 1 720 1656 t
10 I f
(A)1667 1656 w
10 R f
( condition number provides a measure of the sensitivity of)9 2400(. The)1 238 2 1728 1656 t
10 I f
(X)4399 1656 w
10 R f
(to errors in)2 454 1 4493 1656 t
10 I f
(A)4979 1656 w
10 R f
(and)720 1776 w
10 I f
(B)892 1776 w
10 R f
( a machine having)3 741( one is on)3 395( If)1 119( larger it is, the more ill conditioned the matrix.)9 1923(. The)1 233 5 953 1776 t
10 I f
(d)4393 1776 w
10 R f
(decimal digits)1 568 1 4472 1776 t
(of precision, then normally one may expect the elements of)9 2467 1 720 1896 t
10 I f
(X)3223 1896 w
10 R f
( at most)2 337(to have)1 302 2 3320 1896 t
10 I f
(d)3994 1896 w
10 S f
(-)4068 1896 w
10 R f
(log)4139 1896 w
7 R f
(10)4278 1916 w
10 R f
( correct)1 311(\( COND \))2 365 2 4364 1896 t
( section 2.2 for further)4 951( \(See)1 241( the accuracy also depends on the right-hand sides.)8 2138(decimal digits, although)2 990 4 720 2016 t
(details.\))720 2136 w
( the matrix is)3 546( If)1 123( generally urged to use the SS subroutines rather than the LE subroutines.)12 3022(Users are)1 379 4 970 2292 t
( However,)1 442( some computation time can be saved by using the LE subroutine.)11 2653(known to be well conditioned,)4 1225 3 720 2412 t
( In)1 139( with an SS subroutine should rarely be a deterrent.)9 2101(the added cost for computing the condition number)7 2080 3 720 2532 t
( cost of obtaining the condition estimate will be relatively inconse-)10 2680(general, if solving \(2.1\) is expensive, the)6 1640 2 720 2652 t
( the condition estimate will indeed be)6 1541( solving \(2.1\) is rather inexpensive, the cost of obtaining)9 2310(quential. If)1 469 3 720 2772 t
( the time required to compute the condition estimate is roughly)10 2665( Since)1 287( rarely worrisome.)2 764(noticeable, but)1 604 4 720 2892 t
(equivalent to the additional time that would be required by the LE subroutines if)13 3315 1 720 3012 t
10 I f
(B)4068 3012 w
10 R f
(were augmented with)2 877 1 4163 3012 t
( course the ratio is)4 744( Of)1 158(two more columns, the overhead decreases as the number of columns in B increases.)13 3418 3 720 3132 t
( decreases as)2 528( general and symmetric systems, the added cost)7 1939( For)1 195(dependent on the structure of the matrix.)6 1658 4 720 3252 t
( 10)1 127( small systems \(say)3 784( With)1 253(the size of the system increases.)5 1290 4 720 3372 t
10 S f
(\264)3182 3372 w
10 R f
(10\) with one right-hand side, the LE subrou-)7 1795 1 3245 3372 t
( 100)1 208( \(say)1 199(tine might execute in half the time for SS, while for larger systems)12 2745 3 720 3492 t
10 S f
(\264)3880 3492 w
10 R f
(100\) there might be only a)5 1097 1 3943 3492 t
( banded and sparse systems the added cost of SS over LE depends on the sparsity of the)17 3622( For)1 196(10% saving.)1 502 3 720 3612 t
( in)1 113(matrix, i.e. the ratio of zero to nonzero elements)8 1993 2 720 3732 t
10 I f
(A)2861 3732 w
10 R f
( operation counts given in the example in the)8 1878(. The)1 240 2 2922 3732 t
( for BASS for narrow banded matrices with one right)9 2235(Port user sheet for BALE suggest that the penalty)8 2085 2 720 3852 t
(hand side is quite noticeable.)4 1157 1 720 3972 t
10 B f
( Subroutines)1 543(2.2. Lower-Level)1 743 2 720 4212 t
10 R f
( users)1 243( Casual)1 335( the higher level subroutines.)4 1209(The lower-level subroutines are the building blocks for)7 2283 4 970 4368 t
( some people have compli-)4 1100( However)1 421( safely skip this section.)4 982(will probably not use them directly and may)7 1817 4 720 4488 t
( which cannot be handled by the high-level drivers, while others may wish to have more)15 3702(cated problems)1 618 2 720 4608 t
(control over the process.)3 981 1 720 4728 t
( of two phases:)3 602( composed)1 455(The solving of a linear system is)6 1299 3 970 4884 t
( thereof into the product of)5 1128( decomposition \(DC\) of the coefficient matrix or a permutation)9 2611(\(1\) The)1 331 3 970 5040 t
( example, one might factor)4 1073( For)1 189(two or three matrices which have simple structures.)7 2057 3 720 5160 t
10 I f
(A)4064 5160 w
10 R f
(as)4150 5160 w
10 I f
(PA)2678 5340 w
10 S f
(=)2849 5340 w
10 I f
(LU)2953 5340 w
10 R f
(\(2.2\))4849 5340 w
(where)720 5520 w
10 I f
(P)988 5520 w
10 R f
(is a permutation matrix,)3 955 1 1074 5520 t
10 I f
(L)2054 5520 w
10 R f
(is a lower triangular matrix and)5 1256 1 2135 5520 t
10 I f
(U)3416 5520 w
10 R f
(is an upper triangular matrix.)4 1162 1 3513 5520 t
( have been decomposed into the)5 1343( the matrices)2 534( Usually)1 373( solving of the decomposed system.)5 1487(\(2\) The)1 333 5 970 5676 t
( triangular matrices and one usually thinks of forward solving \(FS\) with the first matrix and)15 3810(product of 2)2 510 2 720 5796 t
( if)1 86( Thus)1 250(back solving \(BS\) with the second matrix.)6 1685 3 720 5916 t
10 I f
(PA)2766 5916 w
10 S f
(=)2937 5916 w
10 I f
(LU)3041 5916 w
10 R f
(and we wish to solve)4 838 1 3194 5916 t
10 I f
(Ax)4057 5916 w
10 S f
(=)4211 5916 w
10 I f
(b)4315 5916 w
10 R f
(we will forward solve)3 873 1 970 6072 t
10 I f
(Ly)2698 6252 w
10 S f
(=)2847 6252 w
10 I f
(Pb)2951 6252 w
10 R f
(\(2.3\))4849 6252 w
(and then back solve)3 790 1 970 6432 t
10 I f
(Ux)2723 6612 w
10 S f
(=)2888 6612 w
10 I f
(y)2992 6612 w
10 R f
(\(2.4\))4849 6612 w
( exam-)1 280( For)1 195( subroutine, an FS subroutine, and then a BS subroutine.)9 2309(A typical LE routine would call a DC)7 1536 4 720 6828 t
(ple, GELE calls first GEDC, then GEFS, and finally GEBS.)9 2396 1 720 6948 t
( phase, the matrix)3 728(During the decomposition)2 1054 2 970 7104 t
10 I f
(A)2783 7104 w
10 R f
( permuta-)1 391( The)1 211(is often permuted to promote stability.)5 1563 3 2875 7104 t
(tion matrix)1 472 1 720 7224 t
10 I f
(P)1247 7224 w
10 R f
( an output parameter of the DC)6 1416(can be obtained by unraveling the n-vector INTER,)7 2261 2 1363 7224 t
cleartomark
showpage
saveobj restore
%%EndPage: 2 3
%%Page: 3 4
/saveobj save def
mark
4 pagesetup
10 R f
(- 3 -)2 166 1 2797 480 t
( an)1 119(subroutines. For)1 675 2 720 840 t
10 I f
(n)1539 840 w
10 S f
(\264)1597 840 w
10 I f
(n)1660 840 w
10 R f
(matrix Gaussian elimination proceeds in)4 1615 1 1735 840 t
10 I f
(n)3375 840 w
10 S f
(-)3449 840 w
10 R f
(1 stages. At stage)3 699 1 3520 840 t
10 I f
(k)4244 840 w
10 R f
(a particular row is)3 726 1 4314 840 t
(interchanged with row)2 906 1 720 960 t
10 I f
(k)1655 960 w
10 R f
(and elements below the diagonal in the)6 1583 1 1728 960 t
10 I f
(k)3340 960 w
7 I f
(th)3395 920 w
10 R f
( our subrou-)2 494( In)1 136(column are zeroed out.)3 923 3 3487 960 t
(tines at the)2 437 1 720 1080 t
10 I f
(k)1184 1080 w
7 I f
(th)1239 1040 w
10 R f
( rows)1 222(stage we interchange)2 840 2 1329 1080 t
10 I f
(k)2419 1080 w
10 R f
(and)2491 1080 w
10 I f
(i)2663 1080 w
10 R f
(for some)1 355 1 2719 1080 t
10 I f
(i)3102 1080 w
10 S f
(\263)3171 1080 w
10 I f
(k)3267 1080 w
10 R f
(and set INTER\()2 638 1 3339 1080 t
10 I f
(k)3977 1080 w
10 R f
(\) to)1 139 1 4021 1080 t
10 I f
(i)4188 1080 w
10 R f
( if pivoting is)3 546(. Thus)1 278 2 4216 1080 t
( INTER\()1 359(never done,)1 473 2 720 1200 t
10 I f
(k)1552 1200 w
10 R f
(\)=)1596 1200 w
10 I f
(k)1685 1200 w
10 R f
(for 1)1 198 1 1761 1200 t
10 S f
(\243)1975 1200 w
10 I f
(k)2046 1200 w
10 S f
(\243)2114 1200 w
10 I f
(n)2185 1200 w
10 S f
(-)2259 1200 w
10 R f
( interchanges after the)3 904(1. Because)1 464 2 2330 1200 t
10 I f
(k)3730 1200 w
7 I f
(th)3785 1160 w
10 R f
(stage of the elimination pro-)4 1160 1 3880 1200 t
(cess are not applied to the first)6 1220 1 720 1320 t
10 I f
(k)1965 1320 w
10 R f
(columns of the)2 594 1 2034 1320 t
10 I f
(L)2653 1320 w
10 R f
(matrix, in the FS subroutines the exact order in which the)10 2305 1 2735 1320 t
( do not use the common alternative)6 1440( that we)2 330( Note)1 251(interchanges were performed must be readily accessible.)6 2299 4 720 1440 t
(of initializing INTER\()2 912 1 720 1560 t
10 I f
(i)1632 1560 w
10 R f
(\) to)1 146 1 1660 1560 t
10 I f
(i)1841 1560 w
10 R f
( INTER\()1 362( In)1 143(and interchanging the elements of that vector.)6 1887 3 1904 1560 t
10 I f
(n)4296 1560 w
10 R f
(\) we return +1 if)4 694 1 4346 1560 t
( of interchanges and)3 806(there has been an even number)5 1237 2 720 1680 t
10 S f
(-)2788 1680 w
10 R f
( information is)2 589( This)1 228(1 if there has been an odd number.)7 1380 3 2843 1680 t
(useful in determining the sign of the determinant \(see the Port user sheet for GELU\).)14 3386 1 720 1800 t
( In)1 146( whether the matrix is singular \(or nearly singular\).)8 2142(During the decomposition phase we decide)5 1782 3 970 1956 t
(exact arithmetic if a matrix were singular the DC subroutines would generate a)12 3257 1 720 2076 t
10 I f
(U)4011 2076 w
10 R f
(matrix with some zero)3 923 1 4117 2076 t
( error,when)1 462( to roundoff)2 481( However, due)2 638(diagonal elements.)1 751 4 720 2196 t
10 I f
(A)3080 2196 w
10 R f
(is singular generally no diagonal element of)6 1771 1 3169 2196 t
10 I f
(U)4968 2196 w
10 R f
(will be exactly zero, but some diagonal element will be very small relative to the elements of)16 3714 1 720 2316 t
10 I f
(A)4459 2316 w
10 R f
(.)4520 2316 w
( still lower level subroutine \(GELU)5 1428(Since choosing a cutoff criterion for singularity is very difficult, a)10 2642 2 970 2472 t
( These)1 294( which allows the user to specify his own criterion for singularity.)11 2696(for general matrices\) is provided)4 1330 3 720 2592 t
( declare singular any matrix which pro-)6 1592(lowest level subroutines have an additional parameter EPS and they)9 2728 2 720 2712 t
(duces a diagonal element of)4 1142 1 720 2832 t
10 I f
(U)1894 2832 w
10 R f
( the choice of EPS)4 762( Because)1 388( than or equal to EPS.)5 902(whose magnitude is less)3 990 4 1998 2832 t
( numerical analysis than can reasonably be expected of a user, the DC sub-)13 3013(often demands more expertise in)4 1307 2 720 2952 t
(routines compute a default value for EPS, namely)7 1986 1 720 3072 t
(EPS)2390 3287 w
10 S f
(=)2612 3287 w
7 I f
(j)2792 3357 w
10 R f
(max)2716 3287 w
10 S f
(\354)2896 3250 w
(\356)2896 3350 w
7 I f
(i)2995 3387 w
15 S f
(S)2961 3317 w
10 S f
(\357)3066 3304 w
10 I f
(a)3114 3287 w
7 I f
(i j)1 45 1 3175 3307 t
10 S f
(\357)3228 3304 w
(\374)3244 3250 w
(\376)3244 3350 w
(e)3325 3287 w
10 R f
(where)720 3522 w
10 S f
(e)988 3522 w
10 R f
(is the machine precision, and then call the lowest level subroutines.)10 2691 1 1057 3522 t
(In the lowest level subroutine, if for some index)8 1917 1 970 3678 t
10 I f
(i)2912 3678 w
10 S f
(\357)2626 3875 w
10 I f
(u)2674 3858 w
7 I f
(ii)2735 3878 w
10 S f
(\357)2783 3875 w
(\243)2864 3858 w
10 R f
(EPS)2960 3858 w
(then we set)2 451 1 720 4038 t
10 I f
(u)1197 4038 w
7 I f
(ii)1258 4058 w
10 S f
(=)1355 4038 w
10 I f
(sign)1459 4038 w
10 R f
(\()1667 4038 w
10 I f
(u)1708 4038 w
7 I f
(ii)1769 4058 w
10 R f
(\))1825 4038 w
10 S f
(\264)1874 4038 w
10 R f
( EPS is 0, the)4 545( If)1 118( computation continues.)2 962(EPS and)1 343 4 1937 4038 t
10 I f
(i)3932 4038 w
7 I f
(th)3971 3998 w
10 R f
(column of)1 410 1 4061 4038 t
10 I f
(L)4498 4038 w
10 R f
(is set to the)3 459 1 4581 4038 t
10 I f
(i)720 4158 w
7 I f
(th)759 4118 w
10 R f
( general, users who legitimately)4 1301( In)1 140( are avoided.)2 526(column of the identity matrix and numerical problems)7 2218 4 855 4158 t
( eigenvalue inverse iteration)3 1164(wish to solve singular or near singular problems \(e.g. those arising from the)12 3156 2 720 4278 t
( if the matrix is)4 639( However,)1 447(algorithm\) may do so without encountering avoidable overflow and underflow.)9 3234 3 720 4398 t
( naively, users may encounter an overflow condi-)7 2029(complex and the compiler implements complex division)6 2291 2 720 4518 t
(tion if the modulus of a diagonal element of)8 1754 1 720 4638 t
10 I f
(U)2499 4638 w
10 R f
(underflows.)2596 4638 w
( mentioned in section)3 864( As)1 163( level subroutines are the condition estimators \(CE\).)7 2095(Another group of lower)3 948 4 970 4794 t
( the condition)2 554( If)1 118( measures the sensitivity of the solution to errors in the system.)11 2543(\(2.1\), the condition number)3 1105 4 720 4914 t
( people may be under)4 871( Some)1 281(number is rather large, one should be wary of the solution to the linear system.)14 3168 3 720 5034 t
( some very ill con-)4 757( However,)1 442(the impression that the determinant is a good, cheap measure of conditioning.)11 3121 3 720 5154 t
( example the determinant of the)5 1267( For)1 189(ditioned problems have very reasonable determinants.)5 2157 3 720 5274 t
10 I f
(n)4358 5274 w
10 R f
(by)4433 5274 w
10 I f
(n)4558 5274 w
10 R f
(matrix)4633 5274 w
(1)2495 5574 w
10 S f
(-)2595 5574 w
10 R f
(1)2650 5574 w
( .)1 100(. .)1 180 2 2830 5544 t
10 S f
(-)3160 5574 w
10 R f
(1)3215 5574 w
(1)2650 5694 w
10 S f
(-)2750 5694 w
10 R f
(1)2805 5694 w
(. .)1 125 1 2985 5664 t
10 S f
(-)3160 5694 w
10 R f
(1)3215 5694 w
(1)2805 5814 w
10 S f
(-)2905 5814 w
10 R f
(1)2960 5814 w
(.)3085 5784 w
10 S f
(-)3160 5814 w
10 R f
(1)3215 5814 w
(.)3085 5904 w
(1)3060 6054 w
10 S f
(-)3160 6054 w
10 R f
(1)3215 6054 w
(1)3215 6174 w
(is 1, but its condition number is approximately 2)8 1958 1 720 6474 t
7 I f
(n)2683 6434 w
10 R f
( condition number is defined as)5 1267(. The)1 232 2 2726 6474 t
10 S f
(| |)1 48 1 4253 6474 t
10 I f
(A)4309 6474 w
10 S f
( |)1 28( |)1 61(| |)1 48 3 4378 6474 t
10 I f
(A)4523 6474 w
7 S f
(-)4595 6434 w
7 R f
(1)4645 6434 w
10 S f
(| |)1 48 1 4696 6474 t
10 R f
(, where)1 296 1 4744 6474 t
10 S f
(| |)1 48 1 720 6594 t
10 I f
(A)776 6594 w
10 S f
( =)1 104(| |)1 48 2 845 6594 t
7 S f
(| |)1 33 1 1046 6664 t
7 I f
(x)1084 6664 w
7 S f
( =)1 50(| |)1 33 2 1120 6664 t
7 R f
(1)1214 6664 w
10 R f
(max)1061 6594 w
10 S f
(| |)1 48 1 1256 6594 t
10 I f
(Ax)1312 6594 w
10 S f
(| |)1 48 1 1425 6594 t
10 R f
( condition estimator does not compute)5 1560(. Our)1 236 2 1473 6594 t
10 I f
(A)3300 6594 w
7 S f
(-)3372 6554 w
7 R f
(1)3422 6554 w
10 R f
(explicitly, but estimates the size of its)6 1544 1 3496 6594 t
( uses the algorithm defined in [2] which solves)8 1870( It)1 111(largest elements.)1 671 3 720 6764 t
10 I f
(A)2687 6944 w
7 I f
(T)2759 6904 w
10 B f
(x)2814 6944 w
10 S f
(=)2913 6944 w
10 B f
(b)3017 6944 w
10 R f
(\(2.4\))4849 6944 w
(and)720 7124 w
10 I f
(A)2719 7304 w
10 B f
(y)2788 7304 w
10 S f
(=)2887 7304 w
10 B f
(x)2991 7304 w
cleartomark
showpage
saveobj restore
%%EndPage: 3 4
%%Page: 4 5
/saveobj save def
mark
5 pagesetup
10 R f
(- 4 -)2 166 1 2797 480 t
(with the vector)2 611 1 720 840 t
10 B f
(b)1362 840 w
10 R f
( that)1 182(suitably chosen so)2 745 2 1449 840 t
10 B f
(y)2408 840 w
10 R f
(looks like the column of)4 1000 1 2490 840 t
10 I f
(A)3522 840 w
7 S f
(-)3594 800 w
7 R f
(1)3644 800 w
10 R f
( Since)1 279(with the largest elements.)3 1042 2 3719 840 t
( decomposition of)2 724(this algorithm requires a)3 980 2 720 960 t
10 I f
(A)2450 960 w
10 R f
(, each CE subroutine returns that decomposition along with the)9 2529 1 2511 960 t
(condition estimate.)1 761 1 720 1080 t
( NM \(NorM\) sub-)3 733( The)1 210( package.)1 381(Several additional families of subroutines have been included in the)9 2746 4 970 1236 t
( ML \(MuLtiply\) subroutines)3 1196( The)1 225( a matrix.)2 420(routines \(GENM, BANM, BPNM, etc\) return the norm of)8 2479 4 720 1356 t
(\(GEML, BAML, BPML, etc.\) form the product)6 1940 1 720 1476 t
10 I f
(y)2692 1476 w
10 S f
(=)2760 1476 w
10 I f
(Ax)2831 1476 w
10 R f
(, where)1 300 1 2936 1476 t
10 I f
(A)3268 1476 w
10 R f
(is a matrix having the structure indicated)6 1679 1 3361 1476 t
(by the first 2 letters of the subroutine.)7 1504 1 720 1596 t
(The following figure shows the hierarchy of subroutines for general matrices.)10 3096 1 970 1752 t
(GESS GELE)1 1950 1 1520 2052 t
( GEFS GEBS)2 1001( GEDC)1 822(GECE GEFS GEBS)2 1262 3 1070 2532 t
( GELU)1 566( GENM)1 1244(GENM GELU)1 860 3 820 3012 t
cleartomark
showpage
saveobj restore
%%EndPage: 4 5
%%Page: 5 6
/saveobj save def
mark
6 pagesetup
10 R f
(- 5 -)2 166 1 2797 480 t
10 B f
( Structures for Special Cases)4 1226(3. Data)1 330 2 720 840 t
( Matrices)1 401(3.1. Symmetric)1 665 2 720 1080 t
10 R f
(A real matrix)2 538 1 970 1236 t
10 I f
(A)1536 1236 w
10 R f
( if)1 90(is symmetric)1 517 2 1625 1236 t
10 I f
(a)2261 1236 w
7 I f
(i j)1 45 1 2322 1256 t
10 S f
(=)2424 1236 w
10 I f
(a)2528 1236 w
7 I f
(j i)1 45 1 2589 1256 t
10 R f
( subroutines for symmetric matrices are included)6 1976(. Separate)1 422 2 2642 1236 t
( is possible to decompose a symmetric matrix with ele-)9 2282(in PORT because a recent paper [1] has shown it)9 2038 2 720 1356 t
( nonsymmetric matrix.)2 908( a)1 94(mentary transformations twice as efficiently as)5 1871 3 720 1476 t
(The algorithm in [1] for symmetric matrices also requires the storage of only the upper triangular por-)16 4070 1 970 1632 t
( store the upper triangle portion of a symmetric matrix)9 2199( We)1 191(tion of the matrix.)3 731 3 720 1752 t
10 I f
(A)3869 1752 w
10 R f
(by rows in a vector)4 773 1 3957 1752 t
10 I f
(C)4757 1752 w
10 R f
(. For)1 216 1 4824 1752 t
(example the 4)2 560 1 720 1872 t
10 S f
(\264)1288 1872 w
10 R f
(4 matrix)1 336 1 1351 1872 t
( 5 7)2 400(1 3)1 200 2 2580 2172 t
( 4 6)2 400(3 2)1 200 2 2580 2292 t
( 10)1 200( 8)1 200(5 4)1 200 3 2580 2412 t
( 9)1 200(7 6 10)2 400 2 2580 2532 t
(is stored as the vector)4 865 1 720 2832 t
( 1 3 5 7 2 4 6 8 10 9 \))11 858(C = \()2 256 2 995 3072 t
(In FORTRAN, one can pack an)5 1264 1 720 3312 t
10 I f
(n)2009 3312 w
10 S f
(\264)2067 3312 w
10 I f
(n)2130 3312 w
10 R f
(symmetric matrix)1 708 1 2205 3312 t
10 I f
(A)2938 3312 w
10 R f
(into a vector)2 499 1 3024 3312 t
10 I f
(C)3548 3312 w
10 R f
(using the following scheme:)3 1130 1 3640 3312 t
(L=1)1620 3552 w
(DO 10 I=1,N)2 530 1 1620 3672 t
(DO 5 J=I,N)2 469 1 1800 3792 t
(C\(L\)=A\(I,J\))1980 3912 w
(L=L+1)1980 4032 w
(5 CONTINUE)1 820 1 1490 4152 t
(10 CONTINUE)1 690 1 1440 4272 t
( course the exact code depends)5 1242( Of)1 157( complicated.)1 540(Reading a matrix into a packed array is slightly more)9 2131 4 970 4548 t
( upper)1 253( each record of the input file contains one row of the)11 2114( If)1 118(on the information contained in the input file.)7 1835 4 720 4668 t
(triangular portion of the matrix, the matrix can be read in as follows.)12 2740 1 720 4788 t
(LBEGIN=1)1620 5028 w
(DO 10 I=1,N)2 530 1 1620 5148 t
(LEND=LBEGIN+N)1800 5268 w
10 S f
(-)2616 5268 w
10 R f
(I)2671 5268 w
(READ\(5,1\)\(C\(L\),L=LBEGIN,LEND\))1800 5388 w
(1)1440 5508 w
10 S f
(<)1800 5508 w
10 R f
(appropriate FORMAT statement)2 1304 1 1855 5508 t
10 S f
(>)3159 5508 w
10 R f
(LBEGIN=LEND+1)1800 5628 w
(10 CONTINUE)1 690 1 1440 5748 t
( upper and lower triangular por-)5 1317(If each record of the file contains one row of the original matrix \(both its)14 3003 2 720 5988 t
(tion\), the following FORTRAN program fragment which uses an N-vector TEMP, can be used.)13 3811 1 720 6108 t
cleartomark
showpage
saveobj restore
%%EndPage: 5 6
%%Page: 6 7
/saveobj save def
mark
7 pagesetup
10 R f
(- 6 -)2 166 1 2797 480 t
(L=1)1620 960 w
(DO 20 I=1,N)2 530 1 1620 1080 t
(READ\(5,1\)\(TEMP\(J\),J=1,N\))1800 1200 w
(1)1440 1320 w
10 S f
(<)1800 1320 w
10 R f
(appropriate FORMAT statement)2 1304 1 1855 1320 t
10 S f
(>)3159 1320 w
10 R f
(DO 10 J=I,N)2 519 1 1800 1440 t
(C\(L\)=TEMP\(J\))1980 1560 w
(L=L+1)1980 1680 w
(10 CONTINUE)1 870 1 1440 1800 t
(20 CONTINUE)1 690 1 1440 1920 t
( the following FOR-)3 838(If each record of the file contains one row of only the lower triangular portion,)14 3232 2 970 2196 t
(TRAN program can be used.)4 1150 1 720 2316 t
(DO 20 I=1,N)2 580 1 1570 2556 t
(READ\(5,1\)\(TEMP\(J\),J=1,I\))1720 2676 w
(1)1370 2796 w
10 S f
(<)1670 2796 w
10 R f
(appropriate FORMAT statement\))2 1387 1 1775 2796 t
( NOW CONTAINS A\(J,I\))3 1128(C TEMP\(J\))1 689 2 1270 2916 t
(L=I)1720 3036 w
( WILL CONTAIN A\(1,I\))3 1094(C C\(I\))1 483 2 1270 3156 t
(DO 10 J=1,I)2 547 1 1720 3276 t
(C\(L\)=TEMP\(J\))1870 3396 w
(L=L+N)1870 3516 w
10 S f
(-)2176 3516 w
10 R f
(J)2231 3516 w
(10 CONTINUE)1 810 1 1420 3636 t
(20 CONTINUE)1 660 1 1420 3756 t
( matrix is considered complex symmetric if)6 1941(A complex)1 474 2 970 4032 t
10 I f
(a)3444 4032 w
7 I f
(i j)1 45 1 3505 4052 t
10 S f
(=)3607 4032 w
10 I f
(a)3711 4032 w
7 I f
(j i)1 45 1 3772 4052 t
10 R f
(, and complex Hermitian if)4 1215 1 3825 4032 t
10 I f
(a)720 4157 w
7 I f
(i j)1 45 1 781 4177 t
10 S f
(=)883 4157 w
10 I f
(a)987 4157 w
10 S1 f
(_)989 4089 w
7 I f
(j i)1 45 1 1048 4177 t
10 R f
( beginning with)2 656( Subroutines)1 542(, i.e. the imaginary parts have opposite signs across the diagonal.)10 2741 3 1101 4157 t
( are designed for com-)4 933(CSY are designed for complex symmetric matrices and those beginning with CHE)11 3387 2 720 4277 t
(plex Hermitian matrices.)2 990 1 720 4397 t
(In the symmetric decomposition subroutines one decomposes)6 2459 1 970 4553 t
10 I f
(A)3454 4553 w
10 R f
(as)3540 4553 w
10 I f
(PAP)2535 4733 w
7 I f
(T)2729 4693 w
10 S f
(=)2825 4733 w
10 I f
(MDM)2929 4733 w
7 I f
(T)3178 4693 w
10 R f
(where)720 4913 w
10 I f
(P)991 4913 w
10 R f
(is a permutation matrix,)3 967 1 1080 4913 t
10 I f
(M)2076 4913 w
10 R f
(is a unit lower triangular matrix and)6 1461 1 2188 4913 t
10 I f
(D)3678 4913 w
10 R f
(is a block diagonal matrix with)5 1261 1 3779 4913 t
( vector INTER, an output parameter of SYCE, SYDC, and SYMD, indicates)11 3126( The)1 211(blocks of order 1 and 2.)5 983 3 720 5033 t
(not only the permutation array as in the GE decomposition subroutines but also the block structure of)16 4078 1 720 5153 t
10 I f
(D)4825 5153 w
10 R f
(. If)1 143 1 4897 5153 t
10 I f
(D)720 5273 w
10 R f
(contains a 2)2 499 1 828 5273 t
10 S f
(\264)1335 5273 w
10 R f
( beginning at row I, INTER\(I+1\) will be set to)9 1938(2 block)1 308 2 1398 5273 t
10 S f
(-)3679 5273 w
10 R f
(1 and during the decomposition)4 1306 1 3734 5273 t
( D contains a 1)4 611( If)1 119( row I+1.)2 375(phase row INTER\(I\) was interchanged with)5 1758 4 720 5393 t
10 S f
(\264)3591 5393 w
10 R f
(1 block beginning at row I, during)6 1386 1 3654 5393 t
(the decomposition phase row INTER\(I\) was interchanged with row I.)9 2772 1 720 5513 t
10 B f
( Matrices)1 401(3.2. Banded)1 529 2 720 5753 t
10 R f
(A matrix)1 358 1 970 5909 t
10 I f
(A)1353 5909 w
10 R f
(is banded if all its nonzero elements lie on diagonals close to the main diagonal as in)16 3380 1 1439 5909 t
( x)1 100( x)1 850(x x)1 150 3 2180 6209 t
( x x)2 200( x)1 750(x x x)2 250 3 2180 6329 t
( x x x)3 300( x)1 650(x x x)2 250 3 2280 6449 t
( x x x)3 300( x)1 650(x x x)2 250 3 2380 6569 t
( x x)2 200( x)1 750(x x)1 150 3 2480 6689 t
10 B f
( 2)1 75( Fig.)1 875(Fig. 1)1 239 3 2285 6929 t
10 R f
( in PORT to cover general banded matrices because most)9 2458(Separate programs have been included)4 1612 2 970 7265 t
cleartomark
showpage
saveobj restore
%%EndPage: 6 7
%%Page: 7 8
/saveobj save def
mark
8 pagesetup
10 R f
(- 7 -)2 166 1 2797 480 t
( space by taking advantage)4 1083(problems with banded matrices are very large and it is possible to save time and)14 3237 2 720 840 t
(of the structure.)2 629 1 720 960 t
(In our subroutines the user must specify two quantities:)8 2215 1 970 1116 t
10 I f
(m)1220 1296 w
7 I f
(l)1303 1316 w
10 R f
( number of diagonals in the lower triangular)7 1873(, the)1 221 2 1339 1296 t
(portion of)1 413 1 1438 1476 t
10 I f
(A)1892 1476 w
10 R f
(\(3.1\))4849 1476 w
10 I f
(m)1259 1656 w
10 R f
( total number of diagonals.)4 1138(, the)1 221 2 1339 1656 t
(In Figure 1)2 444 1 720 1872 t
10 I f
(m)1189 1872 w
7 I f
(l)1272 1892 w
10 S f
(=)1325 1872 w
10 R f
(2 and)1 219 1 1429 1872 t
10 I f
(m)1673 1872 w
10 S f
(=)1769 1872 w
10 R f
(3, and in Figure 2,)4 733 1 1840 1872 t
10 I f
(m)2598 1872 w
7 I f
(l)2681 1892 w
10 S f
(=)2734 1872 w
10 R f
(3 and)1 219 1 2838 1872 t
10 I f
(m)3082 1872 w
10 S f
(=)3178 1872 w
10 R f
(4.)3249 1872 w
( subrou-)1 334( those)1 239( In)1 136(Our subroutines and data structures are based on those in Martin and Wilkinson[6].)12 3361 4 970 2028 t
(tines each diagonal of)3 879 1 720 2148 t
10 I f
(A)1626 2148 w
10 R f
( in a)2 174(occupies a column of an array; in ours each diagonal is stored)11 2486 2 1714 2148 t
10 I f
(row)4400 2148 w
10 R f
(of an array,)2 458 1 4582 2148 t
( a banded matrix)3 698( Basically)1 432( leftmost \(lowest\) diagonal occupying the first row.)7 2123(with the)1 334 4 720 2268 t
10 I f
(A)4342 2268 w
10 R f
(will be packed)2 602 1 4438 2268 t
(into an)1 275 1 720 2388 t
10 I f
(m)1020 2388 w
10 S f
(\264)1100 2388 w
10 I f
(n)1163 2388 w
10 R f
(array)1238 2388 w
10 I f
(G)1467 2388 w
10 R f
(according to the formula)3 984 1 1564 2388 t
10 I f
(G)1220 2568 w
10 R f
(\()1300 2568 w
10 I f
(m)1341 2568 w
7 I f
(l)1424 2588 w
10 S f
(+)1468 2568 w
10 I f
(j)1539 2568 w
10 S f
(-)1583 2568 w
10 I f
(i)1654 2568 w
10 R f
(,)1690 2568 w
10 I f
(i)1723 2568 w
10 R f
(\))1759 2568 w
10 S f
(=)1849 2568 w
10 I f
(a)1953 2568 w
7 I f
(i j)1 45 1 2014 2588 t
10 I f
(i)1772 2748 w
10 S f
(=)1849 2748 w
10 R f
(1 , 2 ,... ,)4 282 1 1953 2748 t
10 I f
(n)2243 2748 w
10 R f
(,)2301 2748 w
10 I f
(j)1772 2928 w
10 S f
(=)1849 2928 w
10 I f
(i)1953 2928 w
10 S f
(-)2005 2928 w
10 I f
(m)2076 2928 w
7 I f
(l)2159 2948 w
10 S f
(+)2203 2928 w
10 R f
(1 ,... ,)2 191 1 2274 2928 t
10 I f
(i)2473 2928 w
10 S f
(+)2525 2928 w
10 I f
(m)2596 2928 w
10 S f
(-)2692 2928 w
10 I f
(m)2763 2928 w
7 I f
(l)2846 2948 w
10 R f
(\(3.2\))4849 2928 w
(Since the expression)2 839 1 970 3144 t
10 I f
(i)1843 3144 w
10 S f
(-)1895 3144 w
10 I f
(j)1966 3144 w
10 R f
( on any diagonal, equation \(3.2\) turns a diagonal of)9 2134(is a constant)2 512 2 2028 3144 t
10 I f
(A)4709 3144 w
10 R f
(into a)1 235 1 4805 3144 t
(row of)1 269 1 720 3264 t
10 I f
(G)1020 3264 w
10 R f
( the nonzero portions of rows of)6 1317(. Also)1 270 2 1092 3264 t
10 I f
(A)2710 3264 w
10 R f
(are turned into columns of)4 1078 1 2802 3264 t
10 I f
(G)3911 3264 w
10 R f
( example the banded)3 838(. For)1 219 2 3983 3264 t
(matrix)720 3384 w
(11 12 13)2 500 1 2330 3684 t
(21 22 23 24)3 700 1 2330 3804 t
(32 33 34 35)3 700 1 2530 3924 t
(43 44 45 46)3 700 1 2730 4044 t
(54 55 56)2 500 1 2930 4164 t
(65 66)1 300 1 3130 4284 t
(will be stored as)3 652 1 720 4584 t
(21 32 43 54 65)4 900 1 2530 4884 t
(11 22 33 44 55 66)5 1100 1 2330 5004 t
(12 23 34 45 56)4 900 1 2330 5124 t
(13 24 35 46)3 700 1 2330 5244 t
( of our algorithm adds a multiple)6 1328( the basic step)3 566( Since)1 273(This storage scheme is admittedly complicated.)5 1903 4 970 5580 t
( matrix)1 291(of one row of the)4 711 2 720 5700 t
10 I f
(A)1752 5700 w
10 R f
(\(i.e. a column of)3 672 1 1843 5700 t
10 I f
(G)2545 5700 w
10 R f
(\) to another row and since FORTRAN stores 2 dimensional)9 2423 1 2617 5700 t
( the elements)2 531(arrays by columns, this arrangement should prevent excessive paging because in any such step)13 3789 2 720 5820 t
(referenced in)1 522 1 720 5940 t
10 I f
(G)1267 5940 w
10 R f
(will be close together.)3 882 1 1364 5940 t
(The user does not have to `zero' the unused corners in the)11 2396 1 970 6096 t
10 I f
(G)3399 6096 w
10 R f
( will not)2 352(matrix, since our subroutines)3 1184 2 3504 6096 t
(use them.)1 383 1 720 6216 t
( system)1 328(To solve a linear)3 668 2 970 6372 t
10 I f
(AX)1991 6372 w
10 S f
(=)2137 6372 w
10 I f
(B)2208 6372 w
10 R f
(with an)1 297 1 2294 6372 t
10 I f
(n)2616 6372 w
10 S f
(\264)2674 6372 w
10 I f
(n)2737 6372 w
10 R f
(banded matrix)1 574 1 2812 6372 t
10 I f
(A)3411 6372 w
10 R f
(one first decomposes)2 843 1 3497 6372 t
10 I f
(A)4365 6372 w
10 R f
(as)4451 6372 w
10 I f
(PA)2678 6552 w
10 S f
(=)2849 6552 w
10 I f
(LU)2953 6552 w
10 R f
(\(3.3\))4849 6552 w
(where)720 6732 w
10 I f
(L)996 6732 w
10 R f
(is lower triangular and)3 925 1 1085 6732 t
10 I f
(U)2043 6732 w
10 R f
(is upper triangular as in \(2.2\) and)6 1382 1 2149 6732 t
10 I f
(P)3565 6732 w
10 R f
(is a permutation matrix chosen to)5 1380 1 3660 6732 t
( matrix)1 289( The)1 208(promote stability.)1 709 3 720 6852 t
10 I f
(U)1954 6852 w
10 R f
(will have at most)3 695 1 2054 6852 t
10 I f
(m)2777 6852 w
10 R f
(diagonals \(see \(3.1\)\) and)3 992 1 2876 6852 t
10 I f
(L)3895 6852 w
10 R f
(will have at most)3 692 1 3978 6852 t
10 I f
(m)4697 6852 w
7 I f
(l)4780 6872 w
10 R f
(diag-)4835 6852 w
(onals. Thus)1 488 1 720 6972 t
10 I f
(U)1235 6972 w
10 R f
(can be stored in the space which)6 1303 1 1334 6972 t
10 I f
(A)2664 6972 w
10 R f
( space is needed for)4 798(once occupied but extra)3 956 2 2752 6972 t
10 I f
(L)4534 6972 w
10 R f
( the)1 150(. Since)1 300 2 4590 6972 t
(diagonal of)1 461 1 720 7092 t
10 I f
(L)1215 7092 w
10 R f
(contains all 1's, the actual storage space reserved for)8 2174 1 1305 7092 t
10 I f
(L)3513 7092 w
10 R f
(need only be \()3 592 1 3603 7092 t
10 I f
(m)4203 7092 w
7 I f
(l)4286 7112 w
10 S f
(-)4330 7092 w
10 R f
(1 \))1 91 1 4401 7092 t
10 S f
(\264)4508 7092 w
10 I f
(n)4571 7092 w
10 R f
(locations.)4654 7092 w
(The information needed to form)4 1281 1 720 7212 t
10 I f
(P)2026 7212 w
10 R f
(requires only)1 524 1 2112 7212 t
10 I f
(n)2661 7212 w
10 S f
(-)2735 7212 w
10 R f
(1 integer locations.)2 763 1 2806 7212 t
cleartomark
showpage
saveobj restore
%%EndPage: 7 8
%%Page: 8 9
/saveobj save def
mark
9 pagesetup
10 R f
(- 8 -)2 166 1 2797 480 t
(Typically after forming \(3.3\) one forward solves)6 1937 1 970 840 t
10 I f
(LY)2686 1020 w
10 S f
(=)2847 1020 w
10 I f
(PB)2951 1020 w
10 R f
(\(3.4\))4849 1020 w
(and then back solves)3 829 1 720 1200 t
10 I f
(UX)2692 1380 w
10 S f
(=)2874 1380 w
10 I f
(Y .)1 89 1 2978 1380 t
10 R f
(\(3.5\))4849 1380 w
( are combined and there is no need to store either)10 2082(In our BALE subroutine, steps \(3.3\) and \(3.4\))7 1902 2 720 1560 t
10 I f
(P)4741 1560 w
10 R f
(or)4839 1560 w
10 I f
(L)4959 1560 w
10 R f
(.)5015 1560 w
( \(3.4\) separately for those users whose)6 1604(However lower level subroutines are provided for steps \(3.3\) and)9 2716 2 720 1680 t
( the driver BASS, which solves)5 1290( Moreover,)1 475(problems cannot be solved by BALE.)5 1539 3 720 1800 t
10 I f
(AX)4056 1800 w
10 S f
(=)4227 1800 w
10 I f
(B)4331 1800 w
10 R f
(and determines)1 615 1 4425 1800 t
(the condition number of)3 966 1 720 1920 t
10 I f
(A)1712 1920 w
10 R f
( BASS requires more storage than)5 1367( Thus)1 251(demands the separation of the two steps.)6 1623 3 1799 1920 t
(BALE.)720 2040 w
(Our subroutine for computing)3 1211 1 970 2196 t
10 I f
(U)2211 2196 w
10 R f
( in)1 109(in \(3.3\) also keeps track of the number of diagonals)9 2102 2 2313 2196 t
10 I f
(U)4555 2196 w
10 R f
(. Because)1 413 1 4627 2196 t
(of pivoting it can have as many as)7 1397 1 720 2316 t
10 I f
(m)2147 2316 w
10 R f
( have only)2 424(diagonals; but in the absence of pivoting it will)8 1923 2 2249 2316 t
10 I f
(m)4625 2316 w
10 S f
(-)4721 2316 w
10 I f
(m)4792 2316 w
7 I f
(l)4875 2336 w
10 S f
(+)4919 2316 w
10 R f
(1)4990 2316 w
( there is no need during backsolving to worry about diagonals which are all zero, the vari-)16 3638(diagonals. Since)1 682 2 720 2436 t
( diagonals)1 409(able MU determined in the decomposition subroutines tells BABS the number of nonzero)12 3610 2 720 2556 t
10 I f
(U)4765 2556 w
10 R f
(con-)4863 2556 w
(tains.)720 2676 w
( sides,)1 262(Users with problems involving matrices with narrow bands and using only a few right-hand)13 3808 2 970 2832 t
( example in the Port user sheet)6 1232( The)1 206(should acquaint themselves with the relative costs of BASS and BALE.)10 2882 3 720 2952 t
(for BALE gives some computational times which may be used as a guide.)12 2955 1 720 3072 t
10 B f
( Symmetric Positive Definite Matrices)4 1612(3.3. Banded)1 529 2 720 3312 t
10 R f
( beginning with the letters BP are designed for matrices)9 2267(The subroutines)1 645 2 970 3468 t
10 I f
(A)3912 3468 w
10 R f
(with three special proper-)3 1037 1 4003 3468 t
(ties:)720 3588 w
( are banded, i.e. their nonzero elements lie close to the main diagonal.)12 2789(\(1\) They)1 371 2 970 3744 t
( are symmetric, i.e.)3 765(\(2\) They)1 371 2 970 3900 t
10 I f
(a)2131 3900 w
7 I f
(i j)1 45 1 2192 3920 t
10 S f
(=)2294 3900 w
10 I f
(a)2398 3900 w
7 I f
(j i)1 45 1 2459 3920 t
10 R f
(.)2512 3900 w
( are positive definite, i.e. all their eigenvalues are positive.)9 2332(\(3\) They)1 371 2 970 4056 t
( property that the sum of the absolute values of the off diagonal elements in)14 3130(A matrix satisfying the)3 940 2 970 4212 t
(any row is less than the diagonal, is certain to be positive definite.)12 2636 1 720 4332 t
(The matrix)1 441 1 970 4488 t
(4 1)1 200 1 2630 4788 t
(1 4 1)2 350 1 2630 4908 t
(1 4 1)2 350 1 2780 5028 t
(1 4)1 200 1 2930 5148 t
( are often encountered in physical)5 1389( single out these matrices because they)6 1589( We)1 195(satisfies all these properties.)3 1147 4 720 5448 t
( one third the storage and one third the time for the decomposition)12 2647(problems and the BP subroutines requires)5 1673 2 720 5568 t
(phase of the BA subroutines..)4 1182 1 720 5688 t
(Our subroutines assume that)3 1156 1 970 5844 t
10 I f
(A)2158 5844 w
10 R f
( into the)2 344(has been packed)2 667 2 2251 5844 t
10 I f
(m)3295 5844 w
7 I f
(l)3378 5864 w
10 S f
(\264)3414 5844 w
10 I f
(n)3477 5844 w
10 R f
(matrix)3560 5844 w
10 I f
(G)3854 5844 w
10 R f
(according to the following)3 1081 1 3959 5844 t
(rule:)720 5964 w
10 I f
(G)1220 6144 w
10 R f
(\()1300 6144 w
10 I f
(j)1357 6144 w
10 S f
(-)1401 6144 w
10 I f
(i)1472 6144 w
10 S f
(+)1524 6144 w
10 R f
(1 ,)1 83 1 1595 6144 t
10 I f
(i)1686 6144 w
10 R f
(\))1722 6144 w
10 S f
(=)1812 6144 w
10 I f
(a)1916 6144 w
7 I f
(i j)1 45 1 1977 6164 t
10 R f
(for)2104 6144 w
10 I f
(i)2261 6144 w
10 S f
(=)2313 6144 w
10 R f
(1 ,... ,)2 191 1 2384 6144 t
10 I f
(n)2583 6144 w
10 R f
(;)2665 6144 w
10 I f
(j)2709 6144 w
10 S f
(=)2753 6144 w
10 I f
(i)2824 6144 w
10 R f
(, ,.. ,)2 141 1 2860 6144 t
10 I f
(i)3009 6144 w
10 S f
(+)3061 6144 w
10 I f
(m)3132 6144 w
7 I f
(l)3215 6164 w
10 S f
(-)3259 6144 w
10 R f
(1.)3330 6144 w
(where)720 6324 w
10 I f
(m)997 6324 w
7 I f
(l)1080 6344 w
10 R f
( the diagonal of)3 654( Hence)1 314( number of nonzero diagonals on and below the diagonal of A.)11 2611(is the)1 223 4 1142 6324 t
10 I f
(A)4979 6324 w
10 R f
(becomes the first row of)4 970 1 720 6444 t
10 I f
(G)1715 6444 w
10 R f
(, the first subdiagonal becomes the second row, etc.)8 2060 1 1787 6444 t
(Thus the matrix)2 633 1 970 6600 t
cleartomark
showpage
saveobj restore
%%EndPage: 8 9
%%Page: 9 10
/saveobj save def
mark
10 pagesetup
10 R f
(- 9 -)2 166 1 2797 480 t
(8 2 1)2 350 1 2555 1020 t
(2 7 2 1)3 500 1 2555 1140 t
(1 2 7 2 1)4 650 1 2555 1260 t
(1 2 7 2)3 500 1 2705 1380 t
(1 2 8)2 350 1 2855 1500 t
(will be stored as)3 652 1 720 1800 t
(8 7 7 7 8)4 650 1 2555 2100 t
(2 2 2 2)3 500 1 2555 2220 t
(1 1 1)2 350 1 2555 2340 t
(The BP subroutines do not touch the lower right-hand corner.)9 2458 1 970 2676 t
(Our subroutines for banded symmetric positive definite matrices use Gaussian elimination without)11 4070 1 970 2832 t
( of the matrix, is not needed for stability)8 1735( pivoting, which ruins the band structure)6 1720(pivoting. Symmetric)1 865 3 720 2952 t
(because the matrix is positive definite.)5 1537 1 720 3072 t
10 B f
( Matrices)1 401(4. Sparse)1 414 2 720 3312 t
10 R f
(A matrix is considered sparse if not more than about 1% of its elements are nonzero.)15 3377 1 970 3468 t
10 B f
( to PFORT, PORT, and LINPACK)5 1498(5. Relation)1 486 2 720 3708 t
10 R f
( achieves portability by using the PFORT subset of ANSI FORTRAN, as)11 2982(The linear algebra package)3 1088 2 970 3864 t
( [8], and by using the machine constants module from the)10 2440(defined and enforced by the PFORT Verifier)6 1880 2 720 3984 t
( package also depends)3 913( This)1 236( Subroutine Library[4] to characterize the host computer.)7 2334(PORT Mathematical)1 837 4 720 4104 t
( space, and on PORT's error)5 1161(on PORT's dynamic storage allocation module to provide temporary working)9 3159 2 720 4224 t
( three modules consti-)3 900( These)1 293( which may be "fatal" or "recoverable".)6 1603(handling module to respond to errors,)5 1524 4 720 4344 t
(tute the PORT kernel[5], which is in the public domain.)9 2229 1 720 4464 t
( tested the)2 421(Concurrently with the development of the linear algebra subroutines in PORT, the author)12 3649 2 970 4620 t
( fact, the idea of the condition esti-)7 1392( In)1 133( and greatly benefited from the experience.)6 1714(programs in LINPACK[3],)2 1081 4 720 4740 t
( from the)2 368(mator can be attributed directly to that effort. The major differences between the two libraries stem)15 3952 2 720 4860 t
( a result, the PORT library includes high-)7 1731( As)1 172( the philosophy of the PORT library.)6 1535(author's adherence to)2 882 4 720 4980 t
( some cases)2 481( In)1 138( and uses the PORT kernel when applicable.)7 1804(level drivers, checks arguments when possible,)5 1897 4 720 5100 t
( use essentially the same algorithms, which are ele-)8 2113(the two libraries use different data structures but both)8 2207 2 720 5220 t
(gantly explained in the LINPACK User's Guide[3].)6 2066 1 720 5340 t
10 B f
(6. References)1 589 1 720 5580 t
10 R f
( R. Bunch and L. Kaufman,)5 1149(1. J.)1 314 2 720 5736 t
10 I f
(Some Stable Methods for Calculating Inertia and Solving Symmetric)8 2823 1 2217 5736 t
(Linear Systems,)1 633 1 970 5856 t
10 R f
(Math. Comp. 31, January 1977, pp. 163-179.)6 1799 1 1628 5856 t
( K. Cline, C. B. Moler, G. W. Stewart and J. H. Wilkinson,)12 2384(2. A.)1 347 2 720 6012 t
10 I f
( for the Condition Number)4 1074(An Estimate)1 488 2 3478 6012 t
(of a Matrix,)2 475 1 970 6132 t
10 R f
(SIAM J. Numer. Anal.16 April 1979, 368-375.)6 1879 1 1470 6132 t
( G. W. Stewart,)3 636( Dongarra, J. R. Bunch, C. B. Moler,)7 1499(3. J.J.)1 378 3 720 6288 t
10 I f
(Linpack Users' Guide,)2 917 1 3263 6288 t
10 R f
(SIAM, Philadelphia.)1 830 1 4210 6288 t
(1979.)970 6408 w
( Fox, A.D. Hall, and N.L. Schryer,)6 1396(4. P.A.)1 428 2 720 6564 t
10 I f
(The PORT Mathematical Subroutine Library,)4 1837 1 2571 6564 t
10 R f
(ACM Transac-)1 604 1 4436 6564 t
(tions on Mathematical Software 4, pp. 104-126, June 1978.)8 2370 1 970 6684 t
( Fox, A.D. Hall, and N.L. Schryer,)6 1474(5. P.A.)1 428 2 720 6840 t
10 I f
(Algorithm 528:Framework for a Portable Library,)5 2108 1 2663 6840 t
10 R f
(ACM)4812 6840 w
(Transactions on Mathematical Software 4, pp. 177-188, June 1978.)8 2685 1 970 6960 t
( Lawson, R.J. Hanson, D.R. Kincaid, and F.T. Krough,)8 2294(6. C.L.)1 428 2 720 7116 t
10 I f
(Basic Linear Algebra Modules,)3 1294 1 3480 7116 t
10 R f
(ACM)4812 7116 w
(Transactions on Mathematical Software 5, pp. 308-325, Sept. 1979.)8 2705 1 970 7236 t
cleartomark
showpage
saveobj restore
%%EndPage: 9 10
%%Page: 10 11
/saveobj save def
mark
11 pagesetup
10 R f
(- 10 -)2 216 1 2772 480 t
( Martin and J.H. Wilkinson,)4 1151(7. R.S.)1 423 2 720 840 t
10 I f
( and Unsymmetric Band Equations and the)6 1775(Solutions of Symmetric)2 938 2 2327 840 t
(Calculations of Eigenvectors of Band Matrices,)5 1905 1 970 960 t
10 R f
(Numer. Math. 9, pp. 279-301, 1967.)5 1446 1 2900 960 t
( Ryder,)1 294(8. B.G.)1 439 2 720 1116 t
10 I f
(The PFORT Verifier,)2 847 1 1478 1116 t
10 R f
(Software Practice and Experience 4, pp.359-377, October 1974.)7 2552 1 2350 1116 t
( H. Wilkinson,)2 589(9. J.)1 314 2 720 1272 t
10 I f
(The Algebraic Eigenvalue Problem,)3 1437 1 1648 1272 t
10 R f
(Clarendon Press, Oxford, 1965.)3 1265 1 3110 1272 t
cleartomark
showpage
saveobj restore
%%EndPage: 10 11
%%Trailer
done
%%Pages: 11
%%DocumentFonts: Times-Bold Times-Italic Times-Roman Times-Roman Symbol