%!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 1)1 469 1 2825 120 t (GENERAL MATRICES)1 997 1 2561 360 t ( Solve)1 253( Back)1 355(GEBS -)1 477 3 720 720 t ( Estimation)1 459( Condition)1 551(GECE -)1 477 3 720 840 t ( DeComposition)1 784(GEDC -)1 477 2 720 960 t ( Solve)1 253( Forward)1 488(GEFS -)1 477 3 720 1080 t ( Equation solution)2 734( Linear)1 410(GELE -)1 477 3 720 1200 t ( decomposition)1 614( LU)1 283(GELU -)1 477 3 720 1320 t ( MuLtiplication)1 756(GEML -)1 477 2 720 1440 t (GENM - NorM)2 871 1 720 1560 t ( Solution)1 365( System)1 445(GESS -)1 477 3 720 1680 t cleartomark showpage saveobj restore %%EndPage: 1 1 %%Page: 2 2 /saveobj save def mark 2 pagesetup 10 R f (GEBS \320 upper triangular linear system solution)6 1949 1 1905 120 t 9 B f (Purpose:)720 480 w 10 R f ( It)1 112( where A is an upper triangular matrix.)7 1559(GEBS \(GEneral matrix Back-Solve\) solves AX = B)7 2073 3 1296 480 t ( general linear system solution. \(It is used in this)9 1965(can be used for the back solution phase of a)9 1779 2 1296 600 t (way by the routines GESS and GELE.\))6 1562 1 1296 720 t 9 B f (Usage:)720 1080 w 10 R f (CALL GEBS \(N, A, IA, B, IB, NB\))7 1438 1 1296 1080 t (N)1440 1320 w 14 S f (\256)1980 1340 w 10 R f (the number of equations)3 968 1 2196 1320 t (A)1440 1560 w 14 S f (\256)1980 1580 w 10 R f (the array, dimensioned \(IA, KA\) in the calling program,)8 2237 1 2196 1560 t (where IA)1 376 1 2196 1680 t 10 S f (\263)2600 1680 w 10 R f (N and KA)2 416 1 2683 1680 t 10 S f (\263)3127 1680 w 10 R f (N, containing the N)3 797 1 3210 1680 t 10 S f (\264)4032 1680 w 10 R f ( triangular ma-)2 601(N upper)1 327 2 4112 1680 t (trix)2196 1800 w (The strictly lower triangular portion of A is not)8 1887 1 2196 1920 t (used or changed.)2 673 1 2196 2040 t (IA)1440 2280 w 14 S f (\256)1980 2300 w 10 R f (the row \(leading\) dimension of A, as dimensioned in the)9 2253 1 2196 2280 t (calling program)1 635 1 2196 2400 t (B)1440 2640 w 14 S f (\256)1980 2660 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 2640 t (the calling program, where IB)4 1200 1 2196 2760 t 10 S f (\263)3421 2760 w 10 R f (N and KB)2 405 1 3501 2760 t 10 S f (\263)3931 2760 w 10 R f (NB)4011 2760 w 14 S f (\254)1980 3020 w 10 R f (the solution X)2 567 1 2196 3000 t (IB)1440 3240 w 14 S f (\256)1980 3260 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 3240 t (calling program)1 635 1 2196 3360 t (NB)1440 3600 w 14 S f (\256)1980 3620 w 10 R f (the number of right-hand sides)4 1226 1 2196 3600 t 9 B f (Note:)720 3960 w 10 R f ( GEDC, GELU, or)3 773(GEFS and GEBS can be used directly on the output matrix produced by)12 2971 2 1296 3960 t (GECE to solve a general linear system.)6 1567 1 1296 4080 t 9 B f (Error situations:)1 653 1 720 4440 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 4440 t 10 I f (Er-)4907 4440 w (ror Handling)1 531 1 1512 4560 t 10 R f (, Framework Chapter\))2 884 1 2043 4560 t (Number Error)1 1506 1 1800 4800 t (1 N)1 1008 1 1944 5040 t 10 S f (<)2977 5040 w 10 R f (1)3057 5040 w (2 IA)1 1041 1 1944 5280 t 10 S f (<)3010 5280 w 10 R f (N)3090 5280 w (3 IB)1 1036 1 1944 5520 t 10 S f (<)3005 5520 w 10 R f (N)3085 5520 w (4 NB)1 1075 1 1944 5760 t 10 S f (<)3044 5760 w 10 R f (1)3124 5760 w ( matrix with)2 489( singular)1 952(10 + k*)2 306 3 1944 6000 t 10 I f (k)3716 6000 w 7 I f (th)3771 5960 w 10 R f (diagonal element 0.0)2 835 1 3859 6000 t 9 B f (Double-precision version:)1 988 1 720 6360 t 10 R f (DGEBS with A and B declared double precision)7 1939 1 1783 6360 t (Linear Algebra)1 606 1 360 7500 t (\320 2 \320)2 300 1 2730 7620 t (GEBS)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 ( GEBS)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Complex version:)1 678 1 720 1320 t 10 R f (CGEBS with A and B declared complex)6 1615 1 1473 1320 t 9 B f (Storage:)720 1680 w 10 R f (None)1296 1680 w 9 B f (Time:)720 2040 w 10 R f (\( N)1 113 1 1296 2040 t 7 R f (2)1414 2000 w 10 R f (/ 2)1 86 1 1465 2040 t 10 S f (-)1567 2040 w 10 R f (N / 2 \))3 207 1 1638 2040 t 10 S f (\264)1894 2040 w 10 R f (NB additions)1 531 1 1990 2040 t (\( N)1 113 1 1296 2160 t 7 R f (2)1414 2120 w 10 R f (/ 2)1 86 1 1465 2160 t 10 S f (-)1567 2160 w 10 R f (N / 2 \))3 207 1 1638 2160 t 10 S f (\264)1894 2160 w 10 R f (NB multiplications)1 765 1 1990 2160 t (N)1296 2280 w 10 S f (\264)1409 2280 w 10 R f (NB divisions)1 526 1 1505 2280 t 9 B f (See also:)1 333 1 720 2640 t 10 R f (GECE, GEDC, GEFS, GELE, GELU, GESS)5 1794 1 1296 2640 t 9 B f (Author:)720 3000 w 10 R f (Linda Kaufman)1 629 1 1296 3000 t 9 B f (Example:)720 3360 w 10 R f ( is used as part of a package designed for general matrices as in)13 2555(Usually the subroutine GEBS)3 1189 2 1296 3360 t ( triangu-)1 343(the example for GECE. However it may also be used to solve a linear system with a)16 3401 2 1296 3480 t ( this example the last column of the inverse of the)10 2032( In)1 137(lar matrix as in the following example.)6 1575 3 1296 3600 t (triangular matrix)1 674 1 1296 3720 t (1)2485 3960 w 10 S f (-)2685 3960 w 10 R f (1)2740 3960 w 10 S f (-)2940 3960 w 10 R f ( .)1 187(1 . .)2 400 2 2995 3960 t 10 S f (-)3745 3960 w 10 R f (1)3800 3960 w (1)2740 4080 w 10 S f (-)2940 4080 w 10 R f ( .)1 187(1 . .)2 400 2 2995 4080 t 10 S f (-)3745 4080 w 10 R f (1)3800 4080 w ( .)1 187(1 . .)2 400 2 2995 4200 t 10 S f (-)3745 4200 w 10 R f (1)3800 4200 w ( .)1 228( .)1 187(. .)1 200 3 3195 4320 t ( .)1 228(. .)1 212 2 3370 4440 t (1)3545 4560 w 10 S f (-)3745 4560 w 10 R f (1)3800 4560 w (1)3800 4680 w (is computed by setting the right-hand side to the vector \(0, 0, . . . , 0, 1\).)17 2867 1 1296 4920 t ( elements.)1 406(Although the determinant of the matrix is 1, the inverse can have large off-diagonal)13 3338 2 1296 5160 t ( of this matrix increases, the off-diagonal elements of the inverse increase)11 3020(In fact as the size)4 724 2 1296 5280 t (and the matrix becomes more ill-conditioned.)5 1820 1 1296 5400 t (Linear Algebra)1 606 1 4794 7500 t (\320 3 \320)2 300 1 2730 7620 t (GEBS)5144 7740 w cleartomark showpage saveobj restore %%EndPage: 3 3 %%Page: 4 4 /saveobj save def mark 4 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(GEBS February)1 4665 2 360 840 t 8 CW f (INTEGER N, I, J, IWRITE, I1MACH)5 1488 1 2040 1300 t (REAL A\(15,15\), B\(15\))2 960 1 2040 1400 t (N=15)2040 1500 w (C)1656 1600 w (C FORM THE MATRIX AND SET THE RIGHT-HAND SIDE)8 2160 1 1656 1700 t (C TO THE LAST COLUMN OF THE IDENTITY MATRIX)8 2064 1 1656 1800 t (DO 20 I=1,N)2 528 1 2040 1900 t (DO 10 J=I,N)2 528 1 2184 2000 t (A\(I,J\) = -1.0)2 624 1 2328 2100 t (10 CONTINUE)1 816 1 1752 2200 t (A\(I,I\) = 1.0)2 576 1 2184 2300 t ( 0.0)1 192(B\(I\) =)1 384 2 2184 2400 t (20 CONTINUE)1 672 1 1752 2500 t (B\(N\)=1.0)2040 2600 w (C FIND THE LAST COLUMN OF THE INVERSE MATRIX)8 2112 1 1656 2700 t (CALL GEBS\(N,A,15,B,N,1\))1 1104 1 2040 2800 t (IWRITE=I1MACH\(2\))2040 2900 w (WRITE\(IWRITE,21\)\(I,B\(I\),I=1,N\))2040 3000 w ( B\(,I3,3H \)=,F15.4\))2 912(21 FORMAT\(3H)1 720 2 1752 3100 t (STOP)2040 3200 w (END)2040 3300 w 10 R f (Execution on the Honeywell 6000 computer at Bell Labs yields the result:)11 2961 1 1296 3640 t ( 8192.0000)1 575( \)=)1 114(B\( 1)1 200 3 1321 3880 t ( 4096.0000)1 575( \)=)1 114(B\( 2)1 200 3 1321 4000 t ( 2048.0000)1 575( \)=)1 114(B\( 3)1 200 3 1321 4120 t ( 1024.0000)1 575( \)=)1 114(B\( 4)1 200 3 1321 4240 t ( 512.0000)1 550( \)=)1 114(B\( 5)1 200 3 1321 4360 t ( 256.0000)1 550( \)=)1 114(B\( 6)1 200 3 1321 4480 t ( 128.0000)1 550( \)=)1 114(B\( 7)1 200 3 1321 4600 t ( 64.0000)1 525( \)=)1 114(B\( 8)1 200 3 1321 4720 t ( 32.0000)1 525( \)=)1 114(B\( 9)1 200 3 1321 4840 t ( 16.0000)1 525(B\( 10 \)=)2 339 2 1321 4960 t ( 8.0000)1 500(B\( 11 \)=)2 339 2 1321 5080 t ( 4.0000)1 500(B\( 12 \)=)2 339 2 1321 5200 t ( 2.0000)1 500(B\( 13 \)=)2 339 2 1321 5320 t ( 1.0000)1 500(B\( 14 \)=)2 339 2 1321 5440 t ( 1.0000)1 500(B\( 15 \)=)2 339 2 1321 5560 t (Linear Algebra)1 606 1 360 7500 t (\320 4 \320)2 300 1 2730 7620 t (GEBS)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 ( GEBS)1 4305(February 11, 1993)2 735 2 360 840 t (GECE \320 LU decomposition of a general matrix with condition estimation)10 2987 1 1674 1320 t 9 B f (Purpose:)720 1680 w 10 R f (GECE \(GEneral matrix Condition Estimation\) gives a lower bound for the condition number)12 3744 1 1296 1680 t ( using par-)2 435(of a real general matrix A. It also supplies the LU decomposition of the matrix A)15 3309 2 1296 1800 t (tial pivoting and may be used to replace GEDC or GELU in a linear equation package.)15 3460 1 1296 1920 t 9 B f (Usage:)720 2280 w 10 R f (CALL GECE \(N, A, IA, INTER, COND\))6 1664 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 (A)1440 2760 w 14 S f (\256)1980 2780 w 10 R f (the array, dimensioned \(IA, KA\) in the calling program,)8 2237 1 2196 2760 t (where IA)1 373 1 2196 2880 t 10 S f (\263)2594 2880 w 10 R f (N and KA)2 410 1 2674 2880 t 10 S f (\263)3109 2880 w 10 R f (N, containing the N)3 788 1 3189 2880 t 10 S f (\264)4002 2880 w 10 R f (N coef\256cient matrix)2 804 1 4082 2880 t 14 S f (\254)1980 3140 w 10 R f (the LU decomposition of A \(see)5 1284 1 2196 3120 t 8 B f (Note 2)1 219 1 3505 3120 t 10 R f (\))3724 3120 w (IA)1440 3360 w 14 S f (\256)1980 3380 w 10 R f (the row \(leading\) dimension of A, as dimensioned in)8 2106 1 2196 3360 t (the calling program)2 782 1 2196 3480 t (INTER)1440 3720 w 14 S f (\254)1980 3740 w 10 R f ( row interchanges performed)3 1197(an integer vector of length N recording)6 1647 2 2196 3720 t (during the decomposition \(see)3 1207 1 2196 3840 t 8 B f (Note 2)1 219 1 3428 3840 t 10 R f (\))3647 3840 w (COND)1440 4080 w 14 S f (\254)1980 4100 w 10 R f (an estimate of the condition number of A \(see)8 1830 1 2196 4080 t 8 B f (Note 1)1 219 1 4051 4080 t 10 R f (\))4270 4080 w 9 B f (Note 1:)1 278 1 720 4440 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 4440 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 4560 t (of your linear system have)4 1091 1 1296 4680 t 10 B f (d)2420 4680 w 10 R f ( solution might have as few as)6 1264(decimal digits of precision, the)4 1267 2 2509 4680 t 10 B f (d)1296 4800 w 10 S f (-)1377 4800 w 10 R f (log)1457 4800 w 7 R f (10)1596 4820 w 10 R f ( if COND is greater than 10)6 1121( Thus)1 252(\(COND\) correct decimal digits.)3 1270 3 1674 4800 t 7 R f (Bd)4327 4760 w 7 I f (P)4419 4760 w 10 R f ( be)1 120(, there may)2 450 2 4470 4800 t (no correct digits.)2 674 1 1296 4920 t 9 B f (Note 2:)1 278 1 720 5280 t 10 R f ( in A are suitable for input into GEFS and GEBS.)10 2000(INTER and the LU decomposition returned)5 1744 2 1296 5280 t ( matrix,)1 316(The LU decomposition of A satis\256es the equation PA=LU where P is a permutation)13 3428 2 1296 5400 t ( matrix, and U is an upper triangular matrix. On return from)11 2591(L is a unit lower triangular)5 1153 2 1296 5520 t ( INTER \(see the)3 654(GECE, U occupies the upper triangular portion of A, P can be obtained from)13 3090 2 1296 5640 t ( and the elements of L appear permuted in the strictly lower tri-)12 2594(introduction to this chapter\),)3 1150 2 1296 5760 t (angular portion of A. Since the diagonal elements of L are all 1, they are not stored.)16 3341 1 1296 5880 t (Linear Algebra)1 606 1 4794 7500 t (\320 5 \320)2 300 1 2730 7620 t (GECE)5139 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(GECE 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 IA)1 1041 1 1944 2160 t 10 S f (<)3010 2160 w 10 R f (N)3090 2160 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2400 t 9 B f (Double-precision version:)1 988 1 720 2760 t 10 R f (DGECE with A and COND declared double precision)7 2160 1 1783 2760 t 9 B f (Complex version:)1 678 1 720 3120 t 10 R f (CGECE with A declared complex)4 1359 1 1498 3120 t 9 B f (Storage:)720 3480 w 10 R f ( for DGECE, complex for CGECE\) locations of scratch storage in)10 2747(N real \(double precision)3 997 2 1296 3480 t (the dynamic storage stack)3 1034 1 1296 3600 t 9 B f (Time:)720 4010 w 10 R f (3)1356 4080 w (N)1321 3950 w 7 R f (3)1398 3910 w 10 S1 f (___)1306 3980 w 10 S f (+)1506 4010 w 10 R f (2)1626 4080 w (9)1626 3950 w 10 S1 f (_ _)1 80 1 1611 3980 t 10 R f (N)1709 4010 w 7 R f (2)1786 3970 w 10 S f (+)1869 4010 w 10 R f (6)2014 4080 w (19)1989 3950 w 10 S1 f (_ __)1 130 1 1974 3980 t 10 R f (N additions)1 464 1 2122 4010 t (3)1356 4420 w (N)1321 4290 w 7 R f (3)1398 4250 w 10 S1 f (___)1306 4320 w 10 S f (+)1506 4350 w 10 R f (2)1626 4420 w (5)1626 4290 w 10 S1 f (_ _)1 80 1 1611 4320 t 10 R f (N)1709 4350 w 7 R f (2)1786 4310 w 10 S f (+)1869 4350 w 10 R f (6)1989 4420 w (7)1989 4290 w 10 S1 f (_ _)1 80 1 1974 4320 t 10 R f (N multiplications)1 698 1 2072 4350 t (2)1356 4760 w (N)1321 4630 w 7 R f (2)1398 4590 w 10 S1 f (___)1306 4660 w 10 S f (+)1506 4690 w 10 R f (2)1626 4760 w (3)1626 4630 w 10 S1 f (_ _)1 80 1 1611 4660 t 10 R f (N divisions)1 459 1 1709 4690 t 9 B f (Method:)720 5100 w 10 R f (Gaussian elimination with partial pivoting.)4 1714 1 1296 5100 t (See the reference below for the method used to estimate the condition number.)12 3141 1 1296 5220 t (GECE calls GELU after setting EPS to 0.)7 1660 1 1296 5340 t 9 B f (See also:)1 333 1 720 5700 t 10 R f (GEBS, GEDC, GEFS, GELE, GELU, GESS)5 1789 1 1296 5700 t 9 B f (Authors:)720 6060 w 10 R f (Doris Ryan and Linda Kaufman)4 1281 1 1296 6060 t 9 B f (Reference:)720 6420 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 6420 t (tion number,)1 511 1 1296 6540 t 10 I f (SIAM J. Numer. Anal. 16)4 1007 1 1857 6540 t 10 R f (\(1979\), 368-375.)1 674 1 2889 6540 t (Linear Algebra)1 606 1 360 7500 t (\320 6 \320)2 300 1 2730 7620 t (GECE)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 ( GECE)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Example:)720 1320 w 10 R f ( which can be used to)5 872(The example below is an encoding of the iterative re\256nement algorithm)10 2872 2 1296 1320 t ( highly accurate solution to a system of linear equations with an ill-conditioned)12 3404(obtain a)1 340 2 1296 1440 t ( not excessively high, the program usually re-)7 1895(coef\256cient matrix. If the condition number is)6 1849 2 1296 1560 t (turns a solution that is accurate to the working precision of the machine.)12 2882 1 1296 1680 t (The iterative re\256nement algorithm is essentially:)5 1940 1 1296 1920 t ( Ax = b)3 303(\(1\) Solve)1 408 2 1980 2160 t ( tol =)2 212(\(2\) Set)1 308 2 1980 2280 t 10 S f (e)2525 2280 w 15 S f (S)2626 2310 w 10 S f (\357)2764 2297 w 10 R f (x)2812 2280 w 7 I f (i)2873 2300 w 10 S f (\357)2901 2297 w 10 R f (where)2160 2410 w 10 S f (e)2428 2410 w 10 R f (is the precision of the machine)5 1223 1 2497 2410 t ( in double precision the residual)5 1279(\(3\) Compute)1 547 2 1980 2530 t (r = Ax)2 261 1 2310 2650 t 10 S f (-)2596 2650 w 10 R f (b)2676 2650 w ( A)1 97(\(4\) Solve)1 408 2 1980 2770 t 10 S f (d)2510 2770 w 10 R f (x = r)2 189 1 2559 2770 t ( norm =)2 317(\(5\) Compute)1 547 2 1980 2890 t 15 S f (S)2885 2920 w 10 S f (\357)3023 2907 w (d)3071 2890 w 10 R f (x)3128 2890 w 7 I f (i)3189 2910 w 10 S f (\357)3217 2907 w 10 R f ( x to x +)4 334(\(6\) Set)1 308 2 1980 3020 t 10 S f (d)2647 3020 w 10 R f (x)2721 3020 w ( norm)1 236(\(7\) If)1 246 2 1980 3140 t 10 S f (\243)2487 3140 w 10 R f ( else return to step 3)5 807(tol stop,)1 348 2 2567 3140 t ( \(1\) is accomplished using the three lower-level subroutines GECE, GEFS,)10 3061(In our code, step)3 683 2 1296 3380 t ( lower triangular and U is)5 1068( subroutine GECE factors A into LU where L is)9 2001( The)1 215(and GEBS.)1 460 4 1296 3500 t ( A)1 99( Since)1 274( U.)1 124( GEFS forward solves with L and GEBS back solves with)10 2325( Then)1 256(upper triangular.)1 666 6 1296 3620 t ( by GECE and needed in step \(3\) of the algorithm, a copy of the A matrix is)17 3182(is overwritten)1 562 2 1296 3740 t ( \(4\) the decomposition created earlier in GECE is reused and only GEFS and)13 3161( step)1 192(saved. In)1 391 3 1296 3860 t ( is possible that the matrix is so ill-conditioned that the iterative)11 2657( it)1 92( Since)1 283(GEBS are called.)2 712 4 1296 3980 t ( only a)2 302(re\256nement algorithm will diverge, steps \(3\) through \(7\) in our code are performed)12 3442 2 1296 4100 t ( on the number of bits in)6 989( number is chosen to be an upper bound)8 1608( This)1 230(\256nite number of times.)3 917 4 1296 4220 t (the mantissa of the \257oating-point number supported by the machine.)9 2731 1 1296 4340 t ( yet included in)3 643(This algorithm is not)3 858 2 1296 4580 t 8 R f (PORT)2830 4580 w 10 R f (because the double-precision version of the pro-)6 1971 1 3069 4580 t (gram would require the residuals to be computed in extended precision.)10 2859 1 1296 4700 t 8 CW f (INTEGER N, IA, IB, NB, INTER\(5\), IREAD, I1MACH)7 2208 1 1704 5160 t (INTEGER I, J, IWRITE, ITER, IEND)5 1536 1 1704 5260 t (REAL A\(5, 5\), SAVEA\(5, 5\), B\(5\), SAVEB\(5\), R\(5\))7 2256 1 1704 5360 t (REAL COND, BNORM, R1MACH, ABS, RNORM)5 1728 1 1704 5460 t (DOUBLE PRECISION DSDOT)2 1056 1 1704 5560 t (C)1368 5660 w (N=5)1704 5760 w (IA=5)1704 5860 w (IB=5)1704 5960 w (NB=1)1704 6060 w (IREAD=I1MACH\(1\))1704 6160 w (C)1368 6260 w (DO 10 I=1,N)2 528 1 1704 6360 t ( \(A\(I,J\),J=1,N\))1 720(10 READ\(IREAD,11\))1 1008 2 1512 6460 t (11 FORMAT\(1X,5F8.0\))1 1104 1 1512 6560 t (DO 20 I=1,IB)2 576 1 1704 6660 t ( B\(I\))1 240(20 READ\(IREAD,21\))1 1008 2 1512 6760 t (21 FORMAT\(F8.0\))1 912 1 1512 6860 t 10 R f (Linear Algebra)1 606 1 4794 7520 t (\320 7 \320)2 300 1 2730 7640 t (GECE)5139 7760 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(GECE February)1 4665 2 360 840 t 8 CW f (C)1368 1300 w (C SAVE THE MATRIX AND RIGHT-HAND SIDE \(WHICH WILL BE OVERWRITTEN\))10 3120 1 1368 1400 t (C)1368 1500 w (DO 40 I=1,N)2 528 1 1704 1600 t (SAVEB\(I\)=B\(I\))1800 1700 w (DO 30 J=1,N)2 528 1 1800 1800 t (30 SAVEA\(I,J\)=A\(I,J\))1 1248 1 1512 1900 t (40 CONTINUE)1 576 1 1512 2000 t (C)1368 2100 w (C SOLVE AX = B USING SEPARATE CALLS TO GECE, GEFS, GEBS)11 2640 1 1368 2200 t (C)1368 2300 w (CALL GECE\(N,A,IA,INTER,COND\))1 1344 1 1704 2400 t (IWRITE=I1MACH\(2\))1704 2500 w (IF \(COND.GE.1.0/R1MACH\(4\)\) WRITE\(IWRITE,41\))2 2064 1 1704 2600 t ( CONDITION NUMBER HIGH,ACCURATE SOLUTION UNLIKELY\))5 2400(41 FORMAT\(49H)1 672 2 1512 2700 t (C)1368 2800 w (CALL GEFS\(N,A,IA,B,IB,NB,INTER\))1 1488 1 1704 2900 t (C)1368 3000 w (CALL GEBS\(N,A,IA,B,IB,NB\))1 1200 1 1704 3100 t (WRITE\(IWRITE,42\))1704 3200 w ( ESTIMATED CONDITION NUMBER OF THE MATRIX A,\))7 2160(42 FORMAT\(44H)1 672 2 1512 3300 t (WRITE\(IWRITE,43\) COND)1 1008 1 1704 3400 t ( ONE CALL TO GECE = ,E15.7\))6 1296(43 FORMAT\(27H USING)2 1008 2 1512 3500 t (BNORM=0.0)1704 3600 w (WRITE\(IWRITE,44\))1704 3700 w ( THE FIRST SOLUTION X,\))4 1104(44 FORMAT\(/22H)1 720 2 1512 3800 t (WRITE\(IWRITE,45\))1704 3900 w ( \(USING CALLS TO GECE, GEFS, AND GEBS\) = \))9 2016(45 FORMAT\(41H)1 672 2 1512 4000 t (C)1368 4100 w (C COMPUTE NORM OF SOLUTION)4 1248 1 1368 4200 t (C)1368 4300 w (DO 50 I=1,N)2 528 1 1704 4400 t (BNORM=BNORM + ABS\(B\(I\)\))2 1104 1 1848 4500 t ( B\(I\))1 240(50 WRITE\(IWRITE,51\))1 1104 2 1512 4600 t ( 5F20.7\))1 384(51 FORMAT\(1X,)1 816 2 1512 4700 t (C)1368 4800 w (C REFINE THE SOLUTION DEPENDING ON THE LENGTH OF THE MANTISSA)10 2928 1 1368 4900 t (C)1368 5000 w (IEND=I1MACH\(11\))2040 5100 w 8 S f (*)2760 5100 w 8 CW f (IFIX\(R1MACH\(5\)/ALOG10\(2.0\) + 1.0\))2 1584 1 2800 5100 t (DO 90 ITER=1,IEND)2 816 1 1704 5200 t ( RESIDUAL R = B - AX, IN DOUBLE PRECISION)9 1968(C COMPUTE)1 480 2 1368 5300 t (C)1368 5400 w (WRITE\(IWRITE,52\))1848 5500 w ( THE RESIDUAL R = B - AX = \))9 1344(52 FORMAT\(/27H)1 864 2 1512 5600 t (DO 70 I=1,IA)2 576 1 1848 5700 t (DSDOT=0.0)2136 5800 w (DO 60 J=1,N)2 528 1 2136 5900 t ( + DBLE\(SAVEA\(I,J\)\))2 1008( = DSDOT)2 384(60 DSDOT)1 1008 3 1512 6000 t 8 S f (*)3912 6000 w 8 CW f (B\(J\))3952 6000 w (R\(I\) = SAVEB\(I\) - DSDOT)4 1104 1 2136 6100 t ( R\(I\))1 240(70 WRITE\(IWRITE,51\))1 1392 2 1512 6200 t (C)1368 6300 w ( LU)1 144(C SOLVE)1 384 2 1368 6400 t 8 S f (*)1896 6400 w 8 CW f (\(DELTA X\) = R USING SEPARATE CALLS TO GEFS AND GEBS)10 2448 1 1936 6400 t (C)1368 6500 w (CALL GEFS\(N,A,IA,R,IB,NB,INTER\))1 1488 1 1848 6600 t (CALL GEBS\(N,A,IA,R,IB,NB\))1 1200 1 1848 6700 t (C)1368 6800 w ( X)1 96( + DELTA)2 480( NEW SOLUTION X = X)5 912(C THE)1 288 4 1368 6900 t 10 R f (Linear Algebra)1 606 1 360 7560 t (\320 8 \320)2 300 1 2730 7680 t (GECE)360 7800 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 ( GECE)1 4305(February 11, 1993)2 735 2 360 840 t 8 CW f (C)1368 1300 w (WRITE\(IWRITE,71\))1848 1400 w ( X = \))3 288( + DELTA)2 480( THE NEW SOLUTION X = X)6 1104(71 FORMAT\(/36H)1 864 4 1512 1500 t (C)1368 1600 w (C DETERMINE NORM OF CORRECTION AND ADD IN CORRECTION)8 2496 1 1368 1700 t (C)1368 1800 w (RNORM=0.0)1848 1900 w (DO 80 I=1,N)2 528 1 1848 2000 t ( + R\(I\))2 432(B\(I\) = B\(I\))2 528 2 2136 2100 t (RNORM=RNORM + ABS\(R\(I\)\))2 1104 1 2136 2200 t ( B\(I\))1 240(80 WRITE\(IWRITE,51\))1 1392 2 1512 2300 t (C)1368 2400 w (C TEST FOR CONVERGENCE)3 1056 1 1368 2500 t (C)1368 2600 w (IF\(RNORM.LT.R1MACH\(4\))1848 2700 w 8 S f (*)2856 2700 w 8 CW f (BNORM\) GO TO 100)3 768 1 2896 2700 t (90 CONTINUE)1 576 1 1512 2800 t (WRITE\(IWRITE,91\))1704 2900 w ( ITERATIVE IMPROVEMENT FAILED\))3 1440(91 FORMAT\(/29H)1 720 2 1512 3000 t (100 CONTINUE)1 624 1 1464 3100 t (STOP)1704 3200 w (END)1704 3300 w 10 R f (For the input matrix)3 803 1 1008 4500 t (1.)2356 4620 w 10 S f (-)2581 4620 w 10 R f (2. 3. 7.)2 735 1 2636 4620 t 10 S f (-)3616 4620 w 10 R f (9.)3671 4620 w 10 S f (-)2301 4740 w 10 R f (2. 8.)1 355 1 2356 4740 t 10 S f (-)2911 4740 w 10 R f ( 50.)1 375(6. 2.)1 405 2 2966 4740 t (3.)2356 4860 w 10 S f (-)2581 4860 w 10 R f (6. 18.)1 405 1 2636 4860 t 10 S f (-)3191 4860 w 10 R f (15.)3246 4860 w 10 S f (-)3566 4860 w 10 R f (18.)3621 4860 w (7. 2.)1 355 1 2356 4980 t 10 S f (-)2861 4980 w 10 R f ( 174.)1 375(15. 273.)1 455 2 2916 4980 t 10 S f (-)2301 5100 w 10 R f (9. 50.)1 355 1 2356 5100 t 10 S f (-)2861 5100 w 10 R f ( 1667.)1 375(18. 173.)1 455 2 2916 5100 t (with the following right-hand side:)4 1394 1 1008 5460 t (78.)2678 5700 w 10 S f (-)2573 5820 w 10 R f (320.)2628 5820 w 10 S f (-)2623 5940 w 10 R f (81.)2678 5940 w (215.)2628 6060 w 10 S f (-)2473 6180 w 10 R f (10856.)2528 6180 w (the following results were obtained on the Honeywell 6000 computer at Bell Labs:)12 3307 1 1008 6420 t (Linear Algebra)1 606 1 4794 7500 t (\320 9 \320)2 300 1 2730 7620 t (GECE)5139 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(GECE February)1 4665 2 360 840 t 8 CW f (ESTIMATED CONDITION NUMBER OF THE MATRIX A,)6 2064 1 1416 1300 t ( 07)1 144( 0.7263499E)1 624(USING ONE CALL TO GECE =)5 1152 3 1464 1400 t (THE FIRST SOLUTION X,)3 1008 1 1416 1600 t (\(USING CALLS TO GECE, GEFS, AND GEBS\) =)7 1872 1 1416 1700 t 10 S f (-)2388 1940 w 10 R f (5.9872369)2443 1940 w 10 S f (-)2388 2060 w 10 R f (4.9984942)2443 2060 w 10 S f (-)2388 2180 w 10 R f (8.0019742)2443 2180 w (4.9995200)2443 2300 w 10 S f (-)2388 2420 w 10 R f (6.9999477)2443 2420 w 8 CW f (THE RESIDUAL R = B)4 864 1 1416 2760 t 8 S f (-)2328 2760 w 8 CW f (AX =)1 192 1 2420 2760 t 10 R f (0.0000020)2443 3000 w 10 S f (-)2388 3120 w 10 R f (0.0000190)2443 3120 w (0.0000215)2443 3240 w 10 S f (-)2388 3360 w 10 R f (0.0000073)2443 3360 w 10 S f (-)2388 3480 w 10 R f (0.0000473)2443 3480 w 8 CW f (THE NEW SOLUTION X = X + DELTA X =)9 1632 1 1416 3820 t 10 S f (-)2388 4060 w 10 R f (6.0000002)2443 4060 w 10 S f (-)2388 4180 w 10 R f (5.0000000)2443 4180 w 10 S f (-)2388 4300 w 10 R f (7.9999999)2443 4300 w (5.0000000)2443 4420 w 10 S f (-)2388 4540 w 10 R f (7.0000000)2443 4540 w 8 CW f (THE RESIDUAL R = B)4 864 1 1416 4880 t 8 S f (-)2328 4880 w 8 CW f (AX =)1 192 1 2420 4880 t 10 R f (0.0000001)2443 5120 w 10 S f (-)2388 5240 w 10 R f (0.0000001)2443 5240 w 10 S f (-)2388 5360 w 10 R f (0.0000004)2443 5360 w (0.0000026)2443 5480 w 10 S f (-)2388 5600 w 10 R f (0.0000011)2443 5600 w 8 CW f (THE NEW SOLUTION X = X + DELTA X =)9 1632 1 1416 5940 t 10 S f (-)2388 6180 w 10 R f (6.0000000)2443 6180 w 10 S f (-)2388 6300 w 10 R f (5.0000000)2443 6300 w 10 S f (-)2388 6420 w 10 R f (8.0000000)2443 6420 w (5.0000000)2443 6540 w 10 S f (-)2388 6660 w 10 R f (7.0000000)2443 6660 w (Linear Algebra)1 606 1 360 7500 t (\320 10 \320)2 350 1 2705 7620 t (GECE)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 ( GECE)1 4305(February 11, 1993)2 735 2 360 840 t ( expected from the estimate of the condi-)7 1680(The \256rst solution above is inaccurate, as would have been)9 2352 2 1008 1320 t ( iterative re\256nement algorithm successfully improved the solution to)8 2737( The)1 206(tion number for the matrix.)4 1089 3 1008 1440 t ( matrix and the right-hand side could be represented exactly in the machine.)12 3042(this problem because the)3 990 2 1008 1560 t ( the input matrix cannot be represented exactly and)8 2048( Often)1 278( not high.\))2 416(\(Also the condition number was)4 1290 4 1008 1680 t ( a slightly in-)3 552(the iterative re\256nement algorithm produces a very accurate, but worthless, solution to)11 3480 2 1008 1800 t (correct problem.)1 659 1 1008 1920 t (Linear Algebra)1 606 1 4794 7500 t (\320 11 \320)2 350 1 2705 7620 t (GECE)5139 7740 w cleartomark showpage saveobj restore %%EndPage: 11 11 %%Page: 12 12 /saveobj save def mark 12 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(GECE February)1 4665 2 360 840 t (GEDC \320 LU decomposition of a general matrix)7 1950 1 2049 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( computes the LU decomposition of a dense general)8 2098(GEDC \(GEneral matrix DeComposition\))3 1646 2 1296 1680 t ( step of the solution of a general)7 1284(matrix using partial pivoting. It is called by GELE as the \256rst)11 2460 2 1296 1800 t (linear system.)1 555 1 1296 1920 t 9 B f (Usage:)720 2280 w 10 R f (CALL GEDC \(N, A, IA, INTER\))5 1342 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 (A)1440 2760 w 14 S f (\256)1980 2780 w 10 R f ( the calling program, where IA)5 1250( in)1 106( KA\))1 202(the array, dimensioned \(IA,)3 1109 4 2196 2760 t 10 S f (\263)4888 2760 w 10 R f (N)4968 2760 w (and KA)1 313 1 2196 2880 t 10 S f (\263)2534 2880 w 10 R f (N, containing the N)3 788 1 2614 2880 t 10 S f (\264)3427 2880 w 10 R f (N coef\256cient matrix)2 804 1 3507 2880 t 14 S f (\254)1980 3140 w 10 R f (the LU decomposition of A \(see)5 1284 1 2196 3120 t 8 B f (Note)3505 3120 w 10 R f (\))3664 3120 w (IA)1440 3360 w 14 S f (\256)1980 3380 w 10 R f (the row \(leading\) dimension of A, as dimensioned in)8 2106 1 2196 3360 t (the calling program)2 782 1 2196 3480 t (INTER)1440 3720 w 14 S f (\254)1980 3740 w 10 R f ( row interchanges performed)3 1197(an integer vector of length N recording)6 1647 2 2196 3720 t (during the decomposition \(see)3 1207 1 2196 3840 t 8 B f (Note)3428 3840 w 10 R f (\))3587 3840 w 9 B f (Note:)720 4200 w 10 R f ( in A are suitable for input into GEFS and GEBS.)10 2000(INTER and the LU decomposition returned)5 1744 2 1296 4200 t ( matrix,)1 316(The LU decomposition of A satis\256es the equation PA=LU where P is a permutation)13 3428 2 1296 4320 t ( return from)2 516( On)1 188(L is a unit lower triangular matrix, and U is an upper triangular matrix.)13 3040 3 1296 4440 t ( the upper triangular portion of A, P can be obtained from INTER \(see the)14 2972(GEDC, U occupies)2 772 2 1296 4560 t ( and the elements of L appear permuted in the strictly lower tri-)12 2594(introduction to this chapter\),)3 1150 2 1296 4680 t (angular portion of A. Since the diagonal elements of L are all 1, they are not stored.)16 3341 1 1296 4800 t 9 B f (Error situations:)1 653 1 720 5160 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 5160 t 10 I f (Er-)4907 5160 w (ror Handling)1 531 1 1512 5280 t 10 R f (, Framework Chapter\))2 884 1 2043 5280 t (Number Error)1 1506 1 1800 5520 t (1 N)1 1008 1 1944 5760 t 10 S f (<)2977 5760 w 10 R f (1)3057 5760 w (2 IA)1 1041 1 1944 6000 t 10 S f (<)3010 6000 w 10 R f (N)3090 6000 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 6240 t 9 B f (Double-precision version:)1 988 1 720 6600 t 10 R f (DGEDC with A declared double precision)5 1694 1 1800 6600 t (Linear Algebra)1 606 1 360 7500 t (\320 12 \320)2 350 1 2705 7620 t (GEDC)360 7740 w cleartomark showpage saveobj restore %%EndPage: 12 12 %%Page: 13 13 /saveobj save def mark 13 pagesetup 10 R f (-- --)1 5544 1 0 120 t 8 R f (PORT)360 600 w 10 R f ( Algebra)1 346(library Linear)1 4463 2 591 600 t ( GEDC)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Complex version:)1 678 1 720 1320 t 10 R f (CGEDC with A declared complex)4 1370 1 1473 1320 t 9 B f (Storage:)720 1680 w 10 R f (None)1296 1680 w 9 B f (Time:)720 2090 w 10 R f (3)1356 2160 w (N)1321 2030 w 7 R f (3)1398 1990 w 10 S1 f (___)1306 2060 w 10 S f (+)1515 2090 w 10 R f (2)1679 2160 w (N)1644 2030 w 7 R f (2)1721 1990 w 10 S1 f (___)1629 2060 w 10 S f (+)1838 2090 w 10 R f (6)1978 2160 w (N)1967 2030 w 10 S1 f (_ __)1 102 1 1952 2060 t 10 R f (additions)2089 2090 w (3)1356 2440 w (N)1321 2310 w 7 R f (3)1398 2270 w 10 S1 f (___)1306 2340 w 10 S f (-)1482 2370 w 10 R f (2)1613 2440 w (N)1578 2310 w 7 R f (2)1655 2270 w 10 S1 f (___)1563 2340 w 10 S f (+)1772 2370 w 10 R f (6)1912 2440 w (N)1901 2310 w 10 S1 f (_ __)1 102 1 1886 2340 t 10 R f (multiplications)2023 2370 w (2)1476 2720 w (\( N)1 113 1 1321 2590 t 7 R f (2)1439 2550 w 10 S f (-)1498 2590 w 10 R f (N \))1 113 1 1569 2590 t 10 S1 f (_ _______)1 391 1 1306 2620 t 10 R f (divisions)1732 2650 w 9 B f (Method:)720 3060 w 10 R f ( =)1 81( calls GELU after setting EPS)5 1231( GEDC)1 329( with partial pivoting.)3 888(Gaussian elimination)1 853 5 1296 3060 t 10 S f (\357\357)4703 3060 w 10 R f (A)4801 3060 w 10 S f (\357\357e)4873 3060 w 10 R f (,)5015 3060 w (where)1296 3180 w 10 S f (e)1573 3180 w 10 R f ( preci-)1 265(is machine precision, i.e. the value returned by R1MACH\(4\) \(or, for double)11 3124 2 1651 3180 t (sion, by D1MACH\(4\)\).)2 938 1 1296 3300 t 9 B f (See also:)1 333 1 720 3660 t 10 R f (GEBS, GECE, GEFS, GELE, GELU, GESS)5 1778 1 1296 3660 t 9 B f (Author:)720 4020 w 10 R f (Linda Kaufman)1 629 1 1296 4020 t 9 B f (Example:)720 4380 w 10 R f ( special case, how the building blocks, GEDC, GEFS and)9 2334(In this example, we illustrate, for a)6 1410 2 1296 4380 t (GEBS of our linear equation solver can be used to circumvent a limitation in memory space.)15 3695 1 1296 4500 t (We consider a matrix)3 856 1 1296 4740 t 10 I f (A)2177 4740 w 10 R f (which has the special form)4 1070 1 2263 4740 t 10 I f (A)2016 5065 w 10 S f (=)2126 5065 w (\351)2238 4978 w (\357)2238 5078 w (\353)2238 5178 w 10 I f (E)2291 5145 w (C)2288 5005 w (F)2410 5145 w (D)2405 5005 w 10 S f (\371)2477 4978 w (\357)2477 5078 w (\373)2477 5178 w 10 R f (where)1296 5410 w 10 I f (C)1566 5410 w 10 R f (and)1660 5410 w 10 I f (F)1831 5410 w 10 R f (are dense)1 375 1 1919 5410 t 10 I f (n)2321 5410 w 10 S f (\264)2412 5410 w 10 I f (n)2508 5410 w 10 R f (matrices, and)1 534 1 2585 5410 t 10 I f (D)3146 5410 w 10 R f (and)3245 5410 w 10 I f (E)3416 5410 w 10 R f (are)3505 5410 w 10 I f (n)3654 5410 w 10 S f (\264)3745 5410 w 10 I f (n)3841 5410 w 10 R f ( \(Thus)1 286(diagonal matrices.)1 735 2 3919 5410 t 10 I f (D)4968 5410 w 10 R f (and)1296 5530 w 10 I f (E)1469 5530 w 10 R f (can be stored in vectors of length)6 1349 1 1559 5530 t 10 I f (n)2937 5530 w 10 R f (.\) If)1 178 1 2987 5530 t 10 I f (n)3194 5530 w 10 R f ( the problem is to be solved on a)8 1323(is 200, and)2 444 2 3273 5530 t ( be solved by)3 533(computer with only 100K words of data space, the set of linear equations cannot)13 3211 2 1296 5650 t ( because too much addi-)4 1010(either the general subprogram GELE or the band package, BALE,)9 2734 2 1296 5770 t ( sparse matrix package is also eliminated because those)8 2220( The)1 206( required.)1 383(tional storage would be)3 935 4 1296 5890 t ( column indices for each nonzero ele-)6 1573(subroutines demand additional storage for additional)5 2171 2 1296 6010 t (ment.)1296 6130 w (Linear Algebra)1 606 1 4794 7500 t (\320 13 \320)2 350 1 2705 7620 t (GEDC)5128 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(GEDC February)1 4665 2 360 840 t (However, if)1 476 1 1296 1320 t 10 I f (C)1797 1320 w 7 S f (-)1875 1280 w 7 R f (1)1925 1280 w 10 R f (exists then one can solve)4 993 1 1993 1320 t 10 I f (AX)3011 1320 w 10 S f (=)3157 1320 w 10 I f (B)3228 1320 w 10 R f (using the fact that)3 713 1 3314 1320 t 10 S f (\351)2024 1558 w (\357)2024 1658 w (\353)2024 1758 w 10 I f (E)2077 1725 w (C)2074 1585 w (F)2196 1725 w (D)2191 1585 w 10 S f (\371)2263 1558 w (\357)2263 1658 w (\373)2263 1758 w (=)2383 1645 w (\351)2528 1558 w (\357)2528 1658 w (\353)2528 1758 w 10 I f (EC)2578 1730 w 7 S f (-)2717 1690 w 7 R f (1)2767 1690 w 10 I f (I)2677 1580 w (I)2868 1730 w 10 R f (0)2860 1580 w 10 S f (\371)2910 1558 w (\357)2910 1658 w (\373)2910 1758 w (\351)3030 1558 w (\357)3030 1658 w (\353)3030 1758 w 10 R f (0)3088 1730 w 10 I f (C)3080 1580 w (F)3197 1730 w 10 S f (-)3282 1730 w 10 I f (EC)3353 1730 w 7 S f (-)3492 1690 w 7 R f (1)3542 1690 w 10 I f (D)3593 1730 w (D)3395 1580 w 10 S f (\371)3665 1558 w (\357)3665 1658 w (\373)3665 1758 w 10 I f (.)3810 1645 w 10 R f (If)1296 1990 w 10 I f (B)1387 1990 w 10 R f (is partitioned into)2 706 1 1473 1990 t 10 I f (B)2520 2195 w 10 S f (=)2630 2195 w (\351)2742 2108 w (\357)2742 2208 w (\353)2742 2308 w 10 I f (B)2822 2265 w 7 I f (L)2894 2285 w 10 I f (B)2817 2115 w 7 I f (U)2889 2135 w 10 S1 f (_ ___)1 160 1 2802 2165 t 10 S f (\371)2972 2108 w (\357)2972 2208 w (\373)2972 2308 w 10 R f (where)1296 2540 w 10 I f (B)1564 2540 w 7 I f (U)1636 2560 w 10 R f (has length)1 408 1 1719 2540 t 10 I f (n)2152 2540 w 10 R f (, then, applying the inverse,)4 1107 1 2202 2540 t 10 S f (\351)2168 2778 w (\357)2168 2878 w (\353)2168 2978 w 10 I f (EC)2218 2950 w 7 S f (-)2357 2910 w 7 R f (1)2401 2910 w 10 I f (I)2314 2800 w (I)2502 2950 w 10 R f (0)2494 2800 w 10 S f (\371)2544 2778 w (\357)2544 2878 w (\373)2544 2978 w 7 S f (-)2587 2716 w 7 R f (1)2631 2716 w 10 S f (=)2723 2865 w (\351)2843 2778 w (\357)2843 2878 w (\353)2843 2978 w (-)2893 2950 w 10 I f (EC)2964 2950 w 7 S f (-)3103 2910 w 7 R f (1)3147 2910 w 10 I f (I)3025 2800 w (I)3248 2950 w 10 R f (0)3240 2800 w 10 S f (\371)3290 2778 w (\357)3290 2878 w (\373)3290 2978 w 10 R f (,)3468 2865 w (to)1296 3210 w 10 I f (AX)1399 3210 w 10 S f (=)1545 3210 w 10 I f (B)1616 3210 w 10 R f (, we see that X is the solution to)8 1280 1 1677 3210 t 10 S f (\351)1808 3448 w (\357)1808 3548 w (\353)1808 3648 w 10 R f (0)1965 3620 w 10 I f (C)1957 3470 w (F)2173 3620 w 10 S f (-)2283 3620 w 10 I f (E C)1 169 1 2387 3620 t 7 S f (-)2567 3580 w 7 R f (1)2611 3580 w 10 I f (D)2662 3620 w (D)2417 3470 w 10 S f (\371)2734 3448 w (\357)2734 3548 w (\373)2734 3648 w 10 I f (X)2813 3535 w 10 S f (=)2923 3535 w (\351)3035 3448 w (\357)3035 3548 w (\353)3035 3648 w 10 I f (B)3085 3620 w 7 I f (L)3157 3640 w 10 S f (-)3253 3620 w 10 I f (E C)1 169 1 3357 3620 t 7 S f (-)3537 3580 w 7 R f (1)3581 3580 w 10 I f (B)3632 3620 w 7 I f (U)3704 3640 w 10 I f (B)3358 3450 w 7 I f (U)3430 3470 w 10 S f (\371)3762 3448 w (\357)3762 3548 w (\373)3762 3648 w 10 R f (Finally, if)1 395 1 1296 3880 t 10 I f (X)1716 3880 w 10 R f (is partitioned into)2 706 1 1802 3880 t 10 I f (X)2520 4205 w 10 S f (=)2630 4205 w (\351)2742 4118 w (\357)2742 4218 w (\353)2742 4318 w 10 I f (X)2822 4275 w 7 I f (L)2894 4295 w 10 I f (X)2817 4125 w 7 I f (U)2889 4145 w 10 S1 f (_ ___)1 160 1 2802 4175 t 10 S f (\371)2972 4118 w (\357)2972 4218 w (\373)2972 4318 w 10 R f (where)1296 4550 w 10 I f (X)1564 4550 w 7 I f (U)1636 4570 w 10 R f (has length)1 408 1 1719 4550 t 10 I f (n)2152 4550 w 10 R f (, then)1 222 1 2202 4550 t 10 I f (X)2449 4550 w 10 R f (can be found by the following algorithm:)6 1643 1 2535 4550 t (\(1\) replace)1 428 1 1596 4790 t 10 I f (F)2049 4790 w 10 R f (by)2135 4790 w 10 I f (F)2260 4790 w 10 S f (-)2345 4790 w 10 I f (EC)2416 4790 w 7 S f (-)2555 4750 w 7 R f (1)2599 4750 w 10 I f (D)2650 4790 w 10 R f (\(2\) replace)1 428 1 1596 4910 t 10 I f (B)2049 4910 w 7 I f (L)2121 4930 w 10 R f (by)2193 4910 w 10 I f (B)2318 4910 w 7 I f (L)2390 4930 w 10 S f (-)2453 4910 w 10 I f (EC)2524 4910 w 7 S f (-)2663 4870 w 7 R f (1)2707 4870 w 10 I f (B)2758 4910 w 7 I f (U)2830 4930 w 10 R f (\(3\) solve)1 352 1 1596 5030 t 10 I f (FX)1973 5030 w 7 I f (L)2106 5050 w 10 S f (=)2202 5030 w 10 I f (B)2306 5030 w 7 I f (L)2378 5050 w 10 R f (\(4\) replace)1 428 1 1596 5150 t 10 I f (B)2049 5150 w 7 I f (U)2121 5170 w 10 R f (by)2204 5150 w 10 I f (B)2329 5150 w 7 I f (U)2401 5170 w 10 S f (-)2508 5150 w 10 I f (DX)2612 5150 w 7 I f (L)2756 5170 w 10 R f (\(5\) solve)1 352 1 1596 5270 t 10 I f (CX)1973 5270 w 7 I f (U)2112 5290 w 10 S f (=)2219 5270 w 10 I f (B)2323 5270 w 7 I f (U)2395 5290 w 10 R f (Of course in steps \(1\) and \(2\) we do not form)10 1827 1 1296 5510 t 10 I f (C)3150 5510 w 7 S f (-)3228 5470 w 7 R f (1)3272 5470 w 10 R f (explicitly, but solve a system of equations)6 1698 1 3342 5510 t (with)1296 5630 w 10 I f (C)1514 5630 w 10 R f ( steps \(1\),\(2\), and \(5\) solve systems with the same)9 2126( Thus)1 264(as the coef\256cient matrix.)3 1029 3 1621 5630 t ( the right-hand side for step \(5\) is)7 1348( Moreover,)1 470( right-hand sides.)2 695(coef\256cient matrix but different)3 1231 4 1296 5750 t ( the \256rst two systems have been solved so that GELE, the general linear)13 2894(not known until after)3 850 2 1296 5870 t ( correct approach is to)4 902( The)1 209( used to solve all three simultaneously.)6 1576(equation solver, cannot be)3 1057 4 1296 5990 t (compute the LU factorization of)4 1330 1 1296 6110 t 10 I f (C)2661 6110 w 10 R f ( subroutine GEFS)2 737( Then)1 265(once, and to use it three times.)6 1275 3 2763 6110 t (forward solves with the L portion and GEBS back solves with the U portion.)13 3064 1 1296 6230 t (Linear Algebra)1 606 1 360 7500 t (\320 14 \320)2 350 1 2705 7620 t (GEDC)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 ( GEDC)1 4305(February 11, 1993)2 735 2 360 840 t (In the code below,)3 737 1 1296 1320 t 10 I f (X)2058 1320 w 7 I f (U)2130 1340 w 10 R f (is left in the array)4 704 1 2213 1320 t 10 I f (B)2942 1320 w 7 I f (U)3014 1340 w 10 R f (and)3097 1320 w 10 I f (X)3266 1320 w 7 I f (L)3338 1340 w 10 R f (is left in the array)4 704 1 3410 1320 t 10 I f (B)4139 1320 w 7 I f (L)4211 1340 w 10 R f (.)4258 1320 w 8 CW f (INTEGER INTER\(200\), I, J, N)4 1296 1 2040 1660 t (REAL C\(200, 200\), D\(200\), E\(200\), F\(200, 200\))6 2160 1 2040 1760 t (REAL TEMP\(200\), BL\(200\), BU\(200\))3 1536 1 2040 1860 t (N=200)2040 1960 w (.)2040 2060 w (code for filling in C,D,E,F,BL, and BU belongs here)8 2448 1 2040 2160 t (.)2040 2260 w (C DO AN LU DECOMPOSITION OF C)6 1392 1 1656 2360 t (CALL GEDC\(N,C,200,INTER\))1 1152 1 2040 2460 t (C)1656 2560 w (C FORM F - EC\(INVERSE\)D IN F)6 1344 1 1656 2660 t (C)1656 2760 w (DO 30 J=1,N)2 528 1 2040 2860 t (DO 10 I=1,N)2 528 1 2184 2960 t (TEMP\(I\)=0)2328 3060 w (10 CONTINUE)1 816 1 1752 3160 t (TEMP\(J\)=D\(J\))2184 3260 w (CALL GEFS\(N,C,200,TEMP,200,1,INTER\))1 1680 1 2184 3360 t (CALL GEBS\(N,C,200,TEMP,200,1\))1 1392 1 2184 3460 t (C TEMP CONTAINS THE JTH COLUMN OF)6 1584 1 1656 3560 t (C C\(INVERSE\)D)1 624 1 1656 3660 t (DO 20 I=1,N)2 528 1 2184 3760 t (F\(I,J\)=F\(I,J\) - E\(I\)*TEMP\(I\))2 1344 1 2328 3860 t (20 CONTINUE)1 816 1 1752 3960 t (30 CONTINUE)1 672 1 1752 4060 t (C)1656 4160 w (C FORM BL - EC\(INVERSE\)BU)4 1200 1 1656 4260 t (C)1656 4360 w (DO 40 I=1,N)2 528 1 2040 4460 t (TEMP\(I\)=BU\(I\))2184 4560 w (40 CONTINUE)1 672 1 1752 4660 t (CALL GEFS\(N,C,200,TEMP,200,1,INTER\))1 1680 1 2040 4760 t (CALL GEBS\(N,C,200,TEMP,200,1\))1 1392 1 2040 4860 t (DO 50 I=1,N)2 528 1 2040 4960 t (BL\(I\)=BL\(I\) - E\(I\)*TEMP\(I\))2 1248 1 2184 5060 t (50 CONTINUE)1 672 1 1752 5160 t (C)1656 5260 w (C SOLVE FOR LOWER PART OF X)6 1296 1 1656 5360 t (C)1656 5460 w (CALL GELE\(N,F,200,BL,200,1\))1 1296 1 2040 5560 t (C)1656 5660 w (C FORM RIGHT HAND SIDE TO SOLVE FOR UPPER PART OF X)11 2448 1 1656 5760 t (C)1656 5860 w (DO 60 I=1,N)2 528 1 2040 5960 t (BU\(I\)=BU\(I\) - D\(I\)*BL\(I\))2 1152 1 2184 6060 t (60 CONTINUE)1 672 1 1752 6160 t (C)1656 6260 w (C SOLVE FOR UPPER PART OF X)6 1296 1 1656 6360 t (C)1656 6460 w (CALL GEFS\(N,C,200,BU,200,1,INTER\))1 1584 1 2040 6560 t (CALL GEBS\(N,C,200,BU,200,1\))1 1296 1 2040 6660 t (.)2040 6760 w (.)2040 6860 w 10 R f (Linear Algebra)1 606 1 4794 7520 t (\320 15 \320)2 350 1 2705 7640 t (GEDC)5128 7760 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(GEDC February)1 4665 2 360 840 t (GEFS \320 lower \(unit\) triangular linear system solution)7 2185 1 2075 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( triangular)1 424( where A is a unit lower)6 1019( = PB)2 229(GEFS \(GEneral matrix Forward-Solve\) solves AX)5 2072 4 1296 1680 t ( be used for the for-)5 824( can)1 171( It)1 119(matrix, \(i.e. 1's on the diagonal\), and P is a permutation matrix.)11 2630 4 1296 1800 t ( in this way by the routines)6 1113( is used)2 308( \(It)1 148(ward solution phase of a general linear system solver.)8 2175 4 1296 1920 t (GESS and GELE.\))2 752 1 1296 2040 t 9 B f (Usage:)720 2400 w 10 R f (CALL GEFS \(N, A, IA, B, IB, NB, INTER\))8 1771 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 (A)1440 2880 w 14 S f (\256)1980 2900 w 10 R f (the array, dimensioned \(IA, KA\) in the calling program,)8 2237 1 2196 2880 t (where IA)1 373 1 2196 3000 t 10 S f (\263)2594 3000 w 10 R f (N and KA)2 410 1 2674 3000 t 10 S f (\263)3109 3000 w 10 R f (N, containing the N)3 788 1 3189 3000 t 10 S f (\264)4002 3000 w 10 R f (N coef\256cient matrix)2 804 1 4082 3000 t ( the main diagonal\) is not)5 1059(The upper triangular portion of A\(including)5 1785 2 2196 3120 t (used or changed.)2 673 1 2196 3240 t (IA)1440 3480 w 14 S f (\256)1980 3500 w 10 R f (the row \(leading\) dimension of A, as dimensioned in the)9 2253 1 2196 3480 t (calling program)1 635 1 2196 3600 t (B)1440 3840 w 14 S f (\256)1980 3860 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 3840 t (the calling program, where IB)4 1200 1 2196 3960 t 10 S f (\263)3421 3960 w 10 R f (N and KB)2 405 1 3501 3960 t 10 S f (\263)3931 3960 w 10 R f (NB)4011 3960 w 14 S f (\254)1980 4220 w 10 R f (the solution X)2 567 1 2196 4200 t (IB)1440 4440 w 14 S f (\256)1980 4460 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 4440 t (calling program)1 635 1 2196 4560 t (NB)1440 4800 w 14 S f (\256)1980 4820 w 10 R f (the number of right-hand sides)4 1226 1 2196 4800 t (INTER)1440 5040 w 14 S f (\256)1980 5060 w 10 R f ( of length N recording interchanges performed in)7 2106(the integer vector)2 738 2 2196 5040 t ( solve a unit lower triangular system, set)7 1615( To)1 161( or GECE.)2 419(GELU, GEDC,)1 638 4 2196 5160 t (INTER\(J\) = J, J = 1,...,N.)5 1011 1 2196 5340 t 9 B f (Note 1:)1 278 1 720 5820 t 10 R f ( GEDC, GELU, or)3 773(GEFS and GEBS can be used directly on the output matrix produced by)12 2971 2 1296 5820 t (GECE to solve a general linear system.)6 1567 1 1296 5940 t 9 B f (Note 2:)1 278 1 720 6420 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 6420 t (ent right-hand sides,)2 833 1 1296 6540 t 10 I f (not all known in advance,)4 1070 1 2164 6540 t 10 R f ( repeatedly,)1 474(should not call GESS or GELE)5 1297 2 3269 6540 t (but should use the sequence shown in the example on page 3.)11 2452 1 1296 6660 t (Linear Algebra)1 606 1 360 7500 t (\320 16 \320)2 350 1 2705 7620 t (GEFS)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 ( GEFS)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Error situations:)1 648 1 720 1320 t 10 R f (\(All errors in this subprogram are fatal \320)7 1666 1 1512 1320 t (see)1512 1440 w 10 I f (Error Handling)1 631 1 1664 1440 t 10 R f (, Framework Chapter\))2 884 1 2295 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 IA)1 1041 1 1944 2160 t 10 S f (<)3010 2160 w 10 R f (N)3090 2160 w (3 IB)1 1036 1 1944 2400 t 10 S f (<)3005 2400 w 10 R f (N)3085 2400 w (4 NB)1 1075 1 1944 2640 t 10 S f (<)3044 2640 w 10 R f (1)3124 2640 w ( of INTER out of range)5 934(5 elements)1 1291 2 1944 2880 t 9 B f (Double-precision version:)1 988 1 720 3240 t 10 R f (DGEFS with A and B declared double precision)7 1928 1 1800 3240 t 9 B f (Complex version:)1 678 1 720 3600 t 10 R f (CGEFS with A and B declared complex)6 1604 1 1473 3600 t 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)1 113 1 1296 4320 t 7 R f (2)1414 4280 w 10 S f (-)1473 4320 w 10 R f (N \))1 113 1 1544 4320 t 10 S f (\264)1706 4320 w 10 R f ( additions)1 392(NB / 2)2 233 2 1802 4320 t (\( N)1 113 1 1296 4440 t 7 R f (2)1414 4400 w 10 S f (-)1473 4440 w 10 R f (N \))1 113 1 1544 4440 t 10 S f (\264)1706 4440 w 10 R f ( multiplications)1 626(NB / 2)2 233 2 1802 4440 t 9 B f (See also:)1 333 1 720 4800 t 10 R f (GEBS, GECE, GEDC, GELE, GELU, GESS)5 1805 1 1296 4800 t 9 B f (Author:)720 5160 w 10 R f (Linda Kaufman)1 629 1 1296 5160 t (Linear Algebra)1 606 1 4794 7500 t (\320 17 \320)2 350 1 2705 7620 t (GEFS)5155 7740 w cleartomark showpage saveobj restore %%EndPage: 17 17 %%Page: 18 18 /saveobj save def mark 18 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(GEFS February)1 4665 2 360 840 t 9 B f (Example:)720 1320 w 10 R f ( of problems with)3 720(In this example we give a general outline of a program to solve a sequence)14 3024 2 1296 1320 t (the same coef\256cient matrix but different right-hand sides.)7 2296 1 1296 1440 t ( decomposition can)2 781( This)1 230( computes the LU decomposition of the matrix A.)8 2006(The call to GEDC)3 727 4 1296 1680 t ( solutions \(using GEFS\))3 974(be then used repeatedly \(and ef\256ciently\) for the sequence of forward)10 2770 2 1296 1800 t (and back solutions \(using GEBS\) for each set of right-hand sides.)10 2616 1 1296 1920 t 8 CW f (declare matrix with leading dimension IA)5 1920 1 2184 2260 t (declare right-hand side vector)3 1440 1 2184 2360 t (declare integer vector INTER)3 1344 1 2184 2460 t (assign appropriate values to N and IA)6 1776 1 2184 2560 t (.)2184 2660 w (compute or read in the coefficient matrix)6 1968 1 2184 2760 t (.)2184 2860 w (CALL GEDC\(N,A,IA,INTER\))1 1104 1 2184 2960 t ( a right-hand side)3 864(10 compute)1 576 2 1944 3160 t (CALL GEFS\(N,A,IA,B,N,1,INTER\))1 1392 1 2184 3360 t (CALL GEBS\(N,A,IA,B,N,1\))1 1104 1 2184 3460 t (if sequence of problems is not finished go to 10)9 2304 1 2184 3660 t (.)2184 3760 w (.)2184 3860 w (.)2184 3960 w 10 R f (Linear Algebra)1 606 1 360 7500 t (\320 18 \320)2 350 1 2705 7620 t (GEFS)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 ( GEFS)1 4305(February 11, 1993)2 735 2 360 840 t (GELE \320 general linear system solution)5 1601 1 2367 1320 t 9 B f (Purpose:)720 1680 w 10 R f (GELE \(GEneral Linear Equation solution\) solves the system AX = B where A is a dense gen-)16 3744 1 1296 1680 t (eral matrix.)1 460 1 1296 1800 t 9 B f (Usage:)720 2160 w 10 R f (CALL GELE \(N, A, IA, B, IB, NB\))7 1437 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 (A)1440 2640 w 14 S f (\256)1980 2660 w 10 R f (the array, dimensioned \(IA, KA\) in the calling program,)8 2237 1 2196 2640 t (where IA)1 373 1 2196 2760 t 10 S f (\263)2594 2760 w 10 R f (N and KA)2 410 1 2674 2760 t 10 S f (\263)3109 2760 w 10 R f (N, containing the N)3 788 1 3189 2760 t 10 S f (\264)4002 2760 w 10 R f (N coef\256cient matrix)2 804 1 4082 2760 t (A is overwritten during the solution.)5 1455 1 2196 2880 t (IA)1440 3120 w 14 S f (\256)1980 3140 w 10 R f (the row \(leading\) dimension of A, as dimensioned in the)9 2253 1 2196 3120 t (calling program)1 635 1 2196 3240 t (B)1440 3480 w 14 S f (\256)1980 3500 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 3480 t (the calling program, where IB)4 1200 1 2196 3600 t 10 S f (\263)3421 3600 w 10 R f (N and KB)2 405 1 3501 3600 t 10 S f (\263)3931 3600 w 10 R f (NB)4011 3600 w 14 S f (\254)1980 3860 w 10 R f (the solution X)2 567 1 2196 3840 t (IB)1440 4080 w 14 S f (\256)1980 4100 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 4080 t (calling program)1 635 1 2196 4200 t (NB)1440 4440 w 14 S f (\256)1980 4460 w 10 R f (the number of right-hand sides)4 1226 1 2196 4440 t 9 B f (Note 1:)1 278 1 720 4800 t 10 R f ( well-conditioned, the user should use)5 1540(Unless the given matrix, A, is known in advance to be)10 2204 2 1296 4800 t (GESS instead of GELE.)3 966 1 1296 4920 t 9 B f (Note 2:)1 278 1 720 5280 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 5280 t (ent right-hand sides)2 796 1 1296 5400 t 10 I f ( known in advance,)3 780(not all)1 263 2 2121 5400 t 10 R f (should not use GELE, but should call subpro-)7 1848 1 3192 5400 t ( is called once to get)5 853( GEDC)1 329( in GEDC.\))2 472( the example)2 522( \(See)1 233(grams GEDC, GEFS and GEBS.)4 1335 6 1296 5520 t ( to this chapter\) and then the pair, GEFS \(forward)9 1994(the LU decomposition \(see the introduction)5 1750 2 1296 5640 t (solve\) and GEBS \(back solve\), is called for each new right-hand side.)11 2780 1 1296 5760 t (Linear Algebra)1 606 1 4794 7500 t (\320 19 \320)2 350 1 2705 7620 t (GELE)5145 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(GELE 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 IA)1 1041 1 1944 2160 t 10 S f (<)3010 2160 w 10 R f (N)3090 2160 w (3 IB)1 1036 1 1944 2400 t 10 S f (<)3005 2400 w 10 R f (N)3085 2400 w (4 NB)1 1075 1 1944 2640 t 10 S f (<)3044 2640 w 10 R f (1)3124 2640 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2880 t 9 B f (Double-precision version:)1 988 1 720 3240 t 10 R f (DGELE with A and B declared double precision)7 1938 1 1783 3240 t 9 B f (Complex version:)1 678 1 720 3600 t 10 R f (CGELE with A and B declared complex)6 1614 1 1473 3600 t 9 B f (Storage:)720 3960 w 10 R f (N integer locations of scratch storage in the dynamic storage stack)10 2650 1 1296 3960 t 9 B f (Time:)720 4320 w 10 R f (N)1296 4320 w 7 R f (3)1373 4280 w 10 R f (/3 + N)2 256 1 1416 4320 t 7 R f (2)1677 4280 w 10 R f (/2 + N/6 + NB)4 579 1 1720 4320 t 10 S f (\264)2299 4320 w 10 R f (\(N)2354 4320 w 7 R f (2)2464 4280 w 10 S f (-)2507 4320 w 10 R f (N\) additions)1 497 1 2562 4320 t (N)1296 4440 w 7 R f (3)1373 4400 w 10 R f (/3)1416 4440 w 10 S f (-)1494 4440 w 10 R f (N)1549 4440 w 7 R f (2)1626 4400 w 10 R f (/2 + N/6 + NB)4 579 1 1669 4440 t 10 S f (\264)2248 4440 w 10 R f (\(N)2303 4440 w 7 R f (2)2413 4400 w 10 S f (-)2456 4440 w 10 R f (N\) multiplications)1 731 1 2511 4440 t (\(N)1296 4560 w 7 R f (2)1406 4520 w 10 S f (-)1449 4560 w 10 R f (N\)/2 + N)2 361 1 1504 4560 t 10 S f (\264)1865 4560 w 10 R f (NB divisions)1 526 1 1920 4560 t 9 B f (Method:)720 4920 w 10 R f (Gaussian elimination with partial pivoting.)4 1714 1 1296 4920 t (GELE calls GEDC, GEFS, and GEBS.)5 1555 1 1296 5040 t 9 B f (See also:)1 333 1 720 5400 t 10 R f (GEBS, GECE, GEDC, GEFS, GELU, GESS)5 1795 1 1296 5400 t 9 B f (Authors:)720 5760 w 10 R f (Linda Kaufman and Doris Ryan)4 1281 1 1296 5760 t 9 B f (Example:)720 6120 w 10 R f ( the relative ef\256ciencies of GELE)5 1368( First,)1 265(Two things are illustrated in the following example.)7 2111 3 1296 6120 t ( as a function of the size of the system. The example indicates that)13 2742(and GESS are compared)3 1002 2 1296 6240 t ( cost associated with computing the condition number \(GESS\) decreases as the size)12 3394(the extra)1 350 2 1296 6360 t ( small)1 242( For)1 189( N = 90, GELE is only about 10% faster that GESS.)11 2074( When)1 288(of the system increases.)3 951 5 1296 6480 t ( required by)2 506(systems with one right-hand side, however, GESS takes almost twice the time)11 3238 2 1296 6600 t (GELE.)1296 6720 w (Linear Algebra)1 606 1 360 7500 t (\320 20 \320)2 350 1 2705 7620 t (GELE)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 ( GELE)1 4305(February 11, 1993)2 735 2 360 840 t (Secondly, the example illustrates the increase of the error in the solution as the condition)14 3744 1 1296 1320 t ( 3)1 58( number increases from 1.)4 1079( that as the condition)4 877( Notice)1 327(number grows.)1 610 5 1296 1440 t 10 S f (\264)4288 1440 w 10 R f (10)4384 1440 w 7 R f (3)4489 1400 w 10 R f (to 1)1 163 1 4567 1440 t 10 S f (\264)4771 1440 w 10 R f (10)4867 1440 w 7 R f (6)4972 1400 w 10 R f (,)5015 1440 w (the error in the solution grows from 3)7 1501 1 1296 1560 t 10 S f (\264)2838 1560 w 10 R f (10)2934 1560 w 7 S f (-)3045 1520 w 7 R f (7)3095 1520 w 10 R f (to 6)1 153 1 3163 1560 t 10 S f (\264)3357 1560 w 10 R f (10)3453 1560 w 7 S f (-)3564 1520 w 7 R f (4)3614 1520 w 10 R f (.)3657 1560 w (The N)1 252 1 1296 1800 t 10 S f (\264)1573 1800 w 10 R f (N matrix used in the example)5 1179 1 1653 1800 t 10 I f (a)1626 2125 w 7 I f (i j)1 45 1 1687 2145 t 10 S f (=)1789 2125 w (\354)1901 2038 w (\355)1901 2138 w (\356)1901 2238 w 10 I f (i)1991 2205 w 10 S f (-)2043 2205 w 10 I f (j)2114 2205 w 10 S f (+)2191 2205 w 10 R f (1 for)1 207 1 2295 2205 t 10 I f (i)2543 2205 w 10 S f (\263)2579 2205 w 10 I f (j)2650 2205 w (j)1991 2065 w 10 S f (-)2035 2065 w 10 I f (i)2106 2065 w 10 R f (for)2241 2065 w 10 I f (i)2398 2065 w 10 S f (<)2450 2065 w 10 I f (j)2521 2065 w 10 R f ( right-hand side was chosen to make the)7 1670( The)1 215( increases.)1 425(becomes more ill-conditioned as N)4 1434 4 1296 2470 t ( 1's, and the maximum error is computed as)8 1796(solution vector all)2 734 2 1296 2590 t 7 I f (I)3930 2660 w 10 R f (max)3856 2590 w 10 S f (\357)4028 2607 w 10 R f (X \()1 113 1 4076 2590 t 10 I f (I)4197 2590 w 10 R f (\))4238 2590 w 10 S f (-)4320 2590 w 10 R f (1. 0)1 133 1 4416 2590 t 10 S f (\357)4549 2607 w 10 R f (, for the so-)3 475 1 4565 2590 t (lution, X.)1 381 1 1296 2760 t (The timing subroutine, ILAPSZ, on the Honeywell 6000 system has about 1% accuracy.)12 3532 1 1296 3000 t 8 CW f (INTEGER IA, IB, I1MACH, N, I, J, IT, ILAPSZ, IWRITE)9 2448 1 2040 3340 t (REAL A\(100, 100\), AA\(100, 100\), B\(100\), BB\(100\))6 2256 1 2040 3440 t (REAL SUM, ERR, COND, ABS, TIME, TIMES, AMAX1)7 2112 1 2040 3540 t (IA=100)2040 3640 w (IB =100)1 336 1 2040 3740 t (C)1656 3840 w (C GENERATE THE MATRIX AND RIGHT-HAND SIDE)6 1968 1 1656 3940 t (C)1656 4040 w (DO 40 N=10,90,40)2 768 1 2040 4140 t (DO 20 I=1,N)2 528 1 2184 4240 t (SUM=0.0)2328 4340 w (DO 10 J=1,N)2 528 1 2328 4440 t (A\(I,J\)=ABS\(I-J\))2472 4540 w (IF \(I.GE.J\) A\(I,J\)=A\(I,J\) + 1.0)4 1488 1 2472 4640 t (AA\(I,J\)=A\(I,J\))2472 4740 w (SUM=SUM + AA\(I,J\))2 816 1 2472 4840 t (10 CONTINUE)1 960 1 1752 4940 t (B\(I\)=SUM)2328 5040 w (BB\(I\)=SUM)2328 5140 w (20 CONTINUE)1 816 1 1752 5240 t (C)1656 5340 w ( AND TIME IT)3 576( GELE)1 288(C CALL)1 288 3 1656 5440 t (IT =ILAPSZ\(0\))1 624 1 2184 5540 t (CALL GELE\(N,A,IA,B,IB,1\))1 1152 1 2184 5640 t (TIME=FLOAT\(ILAPSZ\(0\)-IT\)/64.0)2184 5740 w (C)1656 5840 w (C COMPUTE THE MAXIMUM ERROR)4 1296 1 1656 5940 t (C)1656 6040 w (ERR=0.0)2184 6140 w (DO 30 I=1,N)2 528 1 2184 6240 t (ERR=AMAX1\(ERR, ABS\(B\(I\)-1.0\)\))1 1392 1 2328 6340 t (30 CONTINUE)1 816 1 1752 6440 t (C)1656 6540 w (C CALL GESS)2 528 1 1656 6640 t (C)1656 6740 w (IT =ILAPSZ\(0\))1 624 1 2184 6840 t 10 R f (Linear Algebra)1 606 1 4794 7500 t (\320 21 \320)2 350 1 2705 7620 t (GELE)5145 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(GELE February)1 4665 2 360 840 t 8 CW f (CALL GESS\(N,AA,IA,BB,IB,1,COND\))1 1488 1 2184 1300 t (TIMES=FLOAT\(ILAPSZ\(0\)-IT\)/64.0)2184 1400 w (IWRITE=I1MACH\(2\))2184 1500 w (WRITE\(IWRITE,31\)N,COND)2184 1600 w ( FOR N= ,I4,20H CONDITION NUMBER = ,E15.7\))7 2016(31 FORMAT\(8H)1 864 2 1752 1700 t (WRITE\(IWRITE,32\)ERR)2184 1800 w ( MAXIMUM ERROR IN SOLUTION IS ,F15.7\))6 1776(32 FORMAT\(30H)1 912 2 1752 1900 t (WRITE\(IWRITE,33\)TIME)2184 2000 w ( TIME IN MILLISECONDS FOR GELE IS ,F10.2\))7 1968(33 FORMAT\(34H)1 912 2 1752 2100 t (WRITE\(IWRITE,34\)TIMES)2184 2200 w ( TIME IN MILLISECONDS FOR GESS IS ,F10.2\))7 1968(34 FORMAT\(34H)1 912 2 1752 2300 t (40 CONTINUE)1 672 1 1752 2400 t (STOP)2040 2500 w (END)2040 2600 w 10 R f ( above was run on the Honeywell 6000 machine at Bell Laboratories, the)12 2986(When the program)2 758 2 1296 2960 t (following was printed.)2 902 1 1296 3080 t 8 CW f ( 04)1 144( 0.1349031E)1 624( CONDITION NUMBER =)3 912( 10)1 240(FOR N=)1 288 5 1704 3420 t ( 0.0000003)1 768(MAXIMUM ERROR IN SOLUTION IS)4 1344 2 1704 3520 t ( 13.27)1 528(TIME IN MILLISECONDS FOR GELE IS)5 1536 2 1704 3620 t ( 22.70)1 528(TIME IN MILLISECONDS FOR GESS IS)5 1536 2 1704 3720 t ( 06)1 144( 0.1674143E)1 624( CONDITION NUMBER =)3 912( 50)1 240(FOR N=)1 288 5 1704 3920 t ( 0.0000704)1 768(MAXIMUM ERROR IN SOLUTION IS)4 1344 2 1704 4020 t ( 499.30)1 528(TIME IN MILLISECONDS FOR GELE IS)5 1536 2 1704 4120 t ( 597.75)1 528(TIME IN MILLISECONDS FOR GESS IS)5 1536 2 1704 4220 t ( 06)1 144( 0.9745692E)1 624( CONDITION NUMBER =)3 912( 90)1 240(FOR N=)1 288 5 1704 4420 t ( 0.0005805)1 768(MAXIMUM ERROR IN SOLUTION IS)4 1344 2 1704 4520 t ( 2602.00)1 528(TIME IN MILLISECONDS FOR GELE IS)5 1536 2 1704 4620 t ( 2919.06)1 528(TIME IN MILLISECONDS FOR GESS IS)5 1536 2 1704 4720 t 10 R f (Linear Algebra)1 606 1 360 7500 t (\320 22 \320)2 350 1 2705 7620 t (GELE)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 ( GELE)1 4305(February 11, 1993)2 735 2 360 840 t (GELU \320 LU decomposition of a general matrix)7 1944 1 2196 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( a dense general)3 669(GELU \(GEneral matrix LU decomposition\) \256nds the LU decomposition of)9 3075 2 1296 1680 t ( allows the user to specify a threshold for considering a ma-)11 2442( It)1 117(matrix using partial pivoting.)3 1185 3 1296 1800 t ( is called by the LU decomposition routines GECE and GEDC.)10 2523( GELU)1 316(trix singular.)1 511 3 1296 1920 t 9 B f (Usage:)720 2280 w 10 R f (CALL GELU \(N, A, IA, INTER, EPS\))6 1559 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 (A)1440 2760 w 14 S f (\256)1980 2780 w 10 R f (the array, dimensioned \(IA, KA\) in the calling program,)8 2237 1 2196 2760 t (where IA)1 373 1 2196 2880 t 10 S f (\263)2594 2880 w 10 R f (N and KA)2 410 1 2674 2880 t 10 S f (\263)3109 2880 w 10 R f (N, containing the N)3 788 1 3189 2880 t 10 S f (\264)4002 2880 w 10 R f (N coef\256cient matrix)2 804 1 4082 2880 t 14 S f (\254)1980 3140 w 10 R f (the LU decomposition of A \(see)5 1284 1 2196 3120 t 8 B f (Note 2)1 219 1 3505 3120 t 10 R f (\))3724 3120 w (IA)1440 3360 w 14 S f (\256)1980 3380 w 10 R f (the row \(leading\) dimension of A, as dimensioned in)8 2106 1 2196 3360 t (the calling program)2 782 1 2196 3480 t (INTER)1440 3720 w 14 S f (\254)1980 3740 w 10 R f (an integer vector of length N recording row interchanges)8 2266 1 2196 3720 t (performed during the decomposition \(see)4 1647 1 2196 3840 t 8 B f (Note 2)1 219 1 3868 3840 t 10 R f (\))4087 3840 w (EPS)1440 4080 w 14 S f (\256)1980 4100 w 10 R f (if A = LU and)4 566 1 2196 4080 t 10 S f (\357)2787 4080 w 10 I f (u)2836 4080 w 7 R f (kk)2897 4100 w 10 S f (\357 <)1 129 1 2999 4080 t 10 R f (EPS, for some 1)3 650 1 3153 4080 t 10 S f (\243)3803 4080 w 10 R f (k)3866 4080 w 10 S f (\243)3924 4080 w 10 R f (N,)3979 4080 w (the matrix is considered singular.)4 1329 1 2196 4200 t 9 B f (Note 1:)1 278 1 720 4560 t 10 R f ( execution of GELU, \(if the matrix has not been found singular\), the value of the de-)16 3387(After the)1 357 2 1296 4560 t (terminant is INTER\(N\))2 946 1 1296 4680 t 10 S f (\264)2307 4680 w 10 R f (A\(1,1\))2427 4680 w 10 S f (\264)2755 4680 w 10 R f (A\(2,2\))2874 4680 w 10 S f (\264)3201 4680 w 10 R f (. . .)2 125 1 3322 4655 t 10 S f (\264)3513 4680 w 10 R f (A\(N,N\) where INTER\(N\) contains)3 1408 1 3632 4680 t (the sign of the permutation \(the number of row interchanges\).)9 2462 1 1296 4800 t 9 B f (Note 2:)1 278 1 720 5160 t 10 R f ( in A are suitable for input into GEFS and GEBS.)10 2000(INTER and the LU decomposition returned)5 1744 2 1296 5160 t ( matrix,)1 316(The LU decomposition of A satis\256es the equation PA=LU where P is a permutation)13 3428 2 1296 5280 t ( matrix, and U is an upper triangular matrix. On return from)11 2591(L is a unit lower triangular)5 1153 2 1296 5400 t ( P can be obtained from INTER \(see the)8 1610(GELU, U occupies the upper triangular portion of A,)8 2134 2 1296 5520 t ( and the elements of L appear permuted in the strictly lower tri-)12 2594(introduction to this chapter\),)3 1150 2 1296 5640 t (angular portion of A. Since the diagonal elements of L are all 1, they are not stored.)16 3341 1 1296 5760 t (Linear Algebra)1 606 1 4794 7500 t (\320 23 \320)2 350 1 2705 7620 t (GELU)5134 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(GELU 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 IA)1 1041 1 1944 2160 t 10 S f (<)3010 2160 w 10 R f (N)3090 2160 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 2400 t 9 B f (Double-precision version:)1 988 1 720 2760 t 10 R f (DGELU with A and EPS declared double precision)7 2055 1 1783 2760 t 9 B f (Complex version:)1 678 1 720 3120 t 10 R f (CGELU with A declared complex)4 1364 1 1473 3120 t 9 B f (Storage:)720 3480 w 10 R f (None)1296 3480 w 9 B f (Time:)720 3890 w 10 R f (3)1356 3960 w (N)1321 3830 w 7 R f (3)1398 3790 w 10 S1 f (___)1306 3860 w 10 S f (-)1515 3890 w 10 R f (2)1679 3960 w (N)1644 3830 w 7 R f (2)1721 3790 w 10 S1 f (___)1629 3860 w 10 S f (+)1838 3890 w 10 R f (6)1978 3960 w (N)1967 3830 w 10 S1 f (_ __)1 102 1 1952 3860 t 10 R f (additions)2089 3890 w (3)1356 4240 w (N)1321 4110 w 7 R f (3)1398 4070 w 10 S1 f (___)1306 4140 w 10 S f (-)1515 4170 w 10 R f (2)1679 4240 w (N)1644 4110 w 7 R f (2)1721 4070 w 10 S1 f (___)1629 4140 w 10 S f (+)1838 4170 w 10 R f (6)1978 4240 w (N)1967 4110 w 10 S1 f (_ __)1 102 1 1952 4140 t 10 R f (multiplications)2089 4170 w (2)1476 4520 w (\( N)1 113 1 1321 4390 t 7 R f (2)1439 4350 w 10 S f (-)1498 4390 w 10 R f (N \))1 113 1 1569 4390 t 10 S1 f (_ _______)1 391 1 1306 4420 t 10 R f (divisions)1732 4450 w 9 B f (Method:)720 4860 w 10 R f (Gaussian elimination with partial pivoting)4 1689 1 1296 4860 t 9 B f (See also:)1 333 1 720 5220 t 10 R f (GEBS, GECE, GEDC, GEFS, GELE, GESS)5 1784 1 1296 5220 t 9 B f (Author:)720 5580 w 10 R f (Linda Kaufman)1 629 1 1296 5580 t 9 B f (Example:)720 5940 w 10 R f ( subroutine uses the)3 836( The)1 220( to compute a determinant.)4 1128(The following subroutine uses GELU)4 1560 4 1296 5940 t ( INTER. Care is taken to avoid over\257ow and)8 1925(stack to obtain space for the integer vector)7 1819 2 1296 6060 t ( is used to decompose a \257oating)6 1337( subroutine UMKFL)2 840( The)1 214(under\257ow during the calculation.)3 1353 4 1296 6180 t (point number, F, into a mantissa, M, and an exponent E such that)12 2604 1 1296 6300 t (F)2592 6540 w 10 S f (=)2697 6540 w 10 R f (Mb)2801 6540 w 7 R f (E)2945 6500 w 10 R f (where b is the base of the machine and 1/b)9 1699 1 1296 6780 t 10 S f (\243)3020 6780 w 10 R f (M)3100 6780 w 10 S f (<)3214 6780 w 10 R f (1.)3294 6780 w (Linear Algebra)1 606 1 360 7500 t (\320 24 \320)2 350 1 2705 7620 t (GELU)360 7740 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 ( GELU)1 4305(February 11, 1993)2 735 2 360 840 t 8 CW f (SUBROUTINE DET\(N,A,IA,DETMAN,IDETEX\))1 1728 1 1992 1300 t (C)1656 1400 w (C THIS SUBROUTINE COMPUTES THE DETERMINANT OF A)7 2256 1 1656 1500 t (C THE RESULT IS GIVEN BY DETMAN*BETA**IDETEX)6 2112 1 1656 1600 t (C WHERE BETA IS THE BASE OF THE MACHINE)8 1872 1 1656 1700 t (C AND DETMAN IS BETWEEN 1/BETA AND 1)7 1728 1 1656 1800 t (C)1656 1900 w (INTEGER N, IA, IDETEX)3 1008 1 1992 2000 t (INTEGER E, IPOINT, ISTKGT, I1MACH, ISIGN, I)6 2064 1 1992 2100 t (INTEGER IN\(1000\))1 768 1 1992 2200 t (REAL A\(IA, N\), DETMAN, BETA, FLOAT, ONOVBE, M, ABS)8 2400 1 1992 2300 t (DOUBLE PRECISION D\(500\))2 1104 1 1992 2400 t (COMMON /CSTAK/ D)2 768 1 1992 2500 t (EQUIVALENCE\(D\(1\),IN\(1\)\))1992 2600 w (C)1656 2700 w (C ALLOCATE SPACE FROM THE STACK FOR THE PIVOT ARRAY)9 2448 1 1656 2800 t (C)1656 2900 w (IPOINT=ISTKGT\(N,2\))1992 3000 w (CALL GELU\(N,A,IA,IN\(IPOINT\),0.0\))1 1536 1 1992 3100 t (C)1656 3200 w (C THE DETERMINANT IS THE PRODUCT OF THE DIAGONAL ELEMENTS)9 2736 1 1656 3300 t (C AND THE LAST ELEMENT OF THE INTERCHANGE ARRAY)8 2256 1 1656 3400 t (C WE TRY TO COMPUTE THIS PRODUCT IN A WAY THAT WILL)11 2448 1 1656 3500 t (C AVOID UNDERFLOW AND OVERFLOW)4 1440 1 1656 3600 t (C)1656 3700 w (BETA=FLOAT\(I1MACH\(10\)\))1992 3800 w (ONOVBE=1.0/BETA)1992 3900 w (ISIGN=IPOINT + N-1)2 864 1 1992 4000 t (DETMAN=IN\(ISIGN\)*ONOVBE)1992 4100 w (IDETEX=1)1992 4200 w (DO 10 I=1,N)2 528 1 1992 4300 t (CALL UMKFL\(A\(I,I\),E,M\))1 1056 1 2184 4400 t (DETMAN=DETMAN*M)2184 4500 w (IDETEX=IDETEX+E)2184 4600 w (IF\(ABS\(DETMAN\).GE.ONOVBE\) GO TO 10)3 1632 1 2184 4700 t (IDETEX=IDETEX-1)2328 4800 w (DETMAN=DETMAN*BETA)2328 4900 w (10 CONTINUE)1 624 1 1752 5000 t (RETURN)1992 5100 w (END)1992 5200 w 10 R f (Linear Algebra)1 606 1 4794 7500 t (\320 25 \320)2 350 1 2705 7620 t (GELU)5134 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(GELU February)1 4665 2 360 840 t ( - vector multiplication)3 919( matrix)1 311(GEML \320)1 408 3 2349 1320 t 9 B f (Purpose:)720 1680 w 10 R f (GEML \(GEneral matrix MuLtiplication\) forms the product)6 2358 1 1296 1680 t 10 I f (Ax)3679 1680 w 10 R f (where)3809 1680 w 10 I f (A)4077 1680 w 10 R f (is a general matrix.)3 765 1 4163 1680 t 9 B f (Usage:)720 2040 w 10 R f (CALL GEML\(N, A, IA, X, B\))5 1223 1 1296 2040 t (N)1440 2280 w 14 S f (\256)1980 2300 w 10 R f (the length of)2 505 1 2196 2280 t 10 I f (x)2726 2280 w 10 R f (A)1440 2520 w 14 S f (\256)1980 2540 w 10 R f (the array, dimensioned \(IA, KA\) in the calling program,)8 2237 1 2196 2520 t (where IA)1 373 1 2196 2640 t 10 S f (\263)2594 2640 w 10 R f (N and KA)2 410 1 2674 2640 t 10 S f (\263)3109 2640 w 10 R f (N, containing the N)3 788 1 3189 2640 t 10 S f (\264)4002 2640 w 10 R f (N coef\256cient matrix)2 804 1 4082 2640 t (IA)1440 2880 w 14 S f (\256)1980 2900 w 10 R f (the row \(leading\) dimension of A, as dimensioned in)8 2106 1 2196 2880 t (the calling program)2 782 1 2196 3000 t (X)1440 3240 w 14 S f (\256)1980 3260 w 10 R f (the vector)1 396 1 2196 3240 t 10 I f (x)2617 3240 w 10 R f (to be multiplied)2 634 1 2686 3240 t (B)1440 3480 w 14 S f (\254)1980 3500 w 10 R f (the vector A)2 493 1 2196 3480 t 10 I f (x)2689 3480 w 9 B f (Error situations:)1 648 1 720 3840 t 10 R f (\(All errors in this subprogram are fatal \320)7 1666 1 1512 3840 t (see)1512 3960 w 10 I f (Error Handling)1 631 1 1664 3960 t 10 R f (, Framework Chapter\))2 884 1 2295 3960 t (Number Error)1 1506 1 1800 4200 t (1 N)1 1108 1 1944 4440 t 10 S f (<)3077 4440 w 10 R f (1)3157 4440 w (2 IA)1 1141 1 1944 4680 t 10 S f (<)3110 4680 w 10 R f (N)3190 4680 w 9 B f (Double-precision version:)1 988 1 720 5040 t 10 R f (DGEML with A, X, and B declared double precision.)8 2138 1 1800 5040 t 9 B f (Complex version:)1 678 1 720 5400 t 10 R f (CGEML with A, X, and B declared complex)7 1789 1 1473 5400 t 9 B f (Time:)720 5760 w 10 R f (N)1296 5760 w 7 R f (2)1373 5720 w 10 R f (additions)1441 5760 w (N)1296 5880 w 7 R f (2)1373 5840 w 10 R f (multiplications)1441 5880 w 9 B f (See also:)1 333 1 720 6240 t 10 R f (GECE, GEDC, GELU, GELE, GESS)4 1499 1 1296 6240 t 9 B f (Author:)720 6600 w 10 R f (Linda Kaufman)1 629 1 1296 6600 t (Linear Algebra)1 606 1 360 7500 t (\320 26 \320)2 350 1 2705 7620 t (GEML)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 ( GEML)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Example:)720 1320 w 10 R f (This example checks the consistency of GEML and GESS, the linear system solver.)12 3351 1 1296 1320 t (First the example uses GEML to compute for a given vector)10 2402 1 1296 1560 t 10 I f (x)3723 1560 w 10 R f (and matrix)1 430 1 3792 1560 t 10 I f (A)4247 1560 w 10 R f (, the vector)2 446 1 4308 1560 t 10 I f (b)2997 1680 w 10 S f (=)3096 1680 w 10 I f (Ax .)1 138 1 3200 1680 t 10 R f (Then the problem is inverted, i.e., GESS is used to \256nd the vector)12 2626 1 1296 1920 t 10 I f (x)3947 1920 w 10 R f (which satis\256es)1 586 1 4016 1920 t 10 I f (Ax)3014 2160 w 10 S f (=)3168 2160 w 10 I f (b)3272 2160 w 10 R f (This)1296 2400 w 10 I f (x)1499 2400 w 10 R f ( 10)1 125( The)1 205(is then compared with the original vector.)6 1667 3 1568 2400 t 10 S f (\264)3565 2400 w 10 R f (10 matrix)1 386 1 3620 2400 t 10 I f (A)4031 2400 w 10 R f (is chosen so that)3 658 1 4117 2400 t 10 I f (a)2778 2725 w 7 I f (i j)1 45 1 2839 2745 t 10 S f (=)2941 2725 w (\354)3053 2638 w (\355)3053 2738 w (\356)3053 2838 w 10 I f (i)3143 2805 w 10 S f (-)3195 2805 w 10 I f (j)3266 2805 w 10 S f (+)3343 2805 w 10 R f (1 for)1 207 1 3447 2805 t 10 I f (i)3695 2805 w 10 S f (\263)3731 2805 w 10 I f (j)3802 2805 w (j)3143 2665 w 10 S f (-)3187 2665 w 10 I f (i)3258 2665 w 10 R f (for)3393 2665 w 10 I f (i)3550 2665 w 10 S f (<)3602 2665 w 10 I f (j)3673 2665 w 10 R f (The vector)1 429 1 1296 3070 t 10 I f (x)1750 3070 w 10 R f (is chosen randomly.)2 802 1 1819 3070 t 8 CW f (INTEGER I, J, IWRITE, I1MACH, N)5 1488 1 2040 3410 t (REAL A\(10, 10\), X\(10\), B\(10\))4 1344 1 2040 3510 t (REAL ERR, SASUM, UNI, COND)4 1248 1 2040 3610 t (N=10)2040 3710 w (C)1656 3810 w (C CONSTRUCT A MATRIX)3 960 1 1656 3910 t (C)1656 4010 w (DO 20 I=1,N)2 528 1 2040 4110 t (DO 10 J=I,N)2 528 1 2184 4210 t (A\(I,J\)=J-I)2328 4310 w (A\(J,I\)=J-I + 1)2 672 1 2328 4410 t (10 CONTINUE)1 816 1 1752 4510 t (20 CONTINUE)1 672 1 1752 4610 t (C)1656 4710 w (C CONSTRUCT A RANDOM VECTOR X)5 1392 1 1656 4810 t (C)1656 4910 w (DO 30 I=1,N)2 528 1 2040 5010 t (X\(I\)=UNI\(0\))2184 5110 w (30 CONTINUE)1 672 1 1752 5210 t (C)1656 5310 w (C FIND THE VECTOR B=AX)4 1056 1 1656 5410 t (C)1656 5510 w (CALL GEML\(N,A,10,X,B\))1 1008 1 1992 5610 t (C)1656 5710 w (C SOLVE THE SYSTEM AX=B)4 1104 1 1656 5810 t (C)1656 5910 w (CALL GESS\(N,A,10,B,N,1,COND\))1 1344 1 1992 6010 t (C)1656 6110 w (C PRINT THE COMPUTED AND TRUE SOLUTION)6 1824 1 1656 6210 t (C)1656 6310 w (IWRITE=I1MACH\(2\))1992 6410 w (WRITE\(IWRITE,31\))1992 6510 w ( SOLUTION\))1 480( COMPUTED)1 528( TRUE SOLUTION)2 672(31 FORMAT\(34H)1 720 4 1752 6610 t (WRITE\(IWRITE,32\)\(X\(I\),B\(I\),I=1,N\))1992 6710 w ( ,2E17.8\))1 432(32 FORMAT\(1H)1 672 2 1752 6810 t (C)1656 6910 w 10 R f (Linear Algebra)1 606 1 4794 7570 t (\320 27 \320)2 350 1 2705 7690 t (GEML)5117 7810 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(GEML February)1 4665 2 360 840 t 8 CW f (C COMPUTE THE RELATIVE ERROR)4 1344 1 1656 1300 t (C)1656 1400 w (ERR=0.0)1992 1500 w (DO 40 I=1,N)2 528 1 1992 1600 t (ERR=ERR + ABS\(B\(I\)-X\(I\)\))2 1152 1 2136 1700 t (40 CONTINUE)1 624 1 1752 1800 t (ERR=ERR/SASUM\(N,X,1\))1992 1900 w (WRITE\(IWRITE,41\)ERR)1992 2000 w ( RELATIVE ERROR IS ,1PE15.7\))4 1344(41 FORMAT\(19H)1 720 2 1752 2100 t (WRITE\(6,42\)COND)1992 2200 w ( CONDITION NUMBER IS ,1PE15.7\))4 1440(42 FORMAT\(21H)1 720 2 1752 2300 t (STOP)1992 2400 w (END)1992 2500 w 10 R f ( on the Honeywell 6000 machine at Bell Laboratories,)8 2174(When the above program was executed)5 1570 2 1296 2840 t (the following was printed:)3 1052 1 1296 2960 t 8 CW f (TRUE SOLUTION COMPUTED SOLUTION)3 1488 1 2064 3300 t (0.22925607E 00 0.22925687E 00)3 1392 1 2064 3400 t (0.76687502E 00 0.76687336E 00)3 1392 1 2064 3500 t (0.68317685E 00 0.68317838E 00)3 1392 1 2064 3600 t (0.50919111E 00 0.50918986E 00)3 1392 1 2064 3700 t (0.87455959E 00 0.87456071E 00)3 1392 1 2064 3800 t (0.64464101E 00 0.64463982E 00)3 1392 1 2064 3900 t (0.84746840E 00 0.84746962E 00)3 1392 1 2064 4000 t (0.35396343E 00 0.35396226E 00)3 1392 1 2064 4100 t (0.39889160E 00 0.39889258E 00)3 1392 1 2064 4200 t (0.45709422E 00 0.45709377E 00)3 1392 1 2064 4300 t (RELATIVE ERROR IS 1.9705190E-06)3 1488 1 2064 4400 t (CONDITION NUMBER IS 1.3490306E 03)4 1584 1 2064 4500 t 10 R f ( the precision of the Honeywell computer suggest)7 2074(The condition number of the matrix and)6 1670 2 1296 4840 t ( in GEML, a relative error of 10)7 1343(that even in the absence of roundoff error)7 1723 2 1296 4960 t 7 S f (-)4373 4920 w 7 R f (5)4423 4920 w 10 R f (would not be)2 540 1 4500 4960 t ( value computed above is quite reasonable.)6 1716(surprising. The)1 630 2 1296 5080 t (Linear Algebra)1 606 1 360 7500 t (\320 28 \320)2 350 1 2705 7620 t (GEML)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 ( GEML)1 4305(February 11, 1993)2 735 2 360 840 t ( of a general matrix)4 781( norm)1 261(GENM \320)1 419 3 2437 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( norm is)2 364( The)1 223( a general matrix A.)4 867(GENM \(GEneral matrix NorM\) computes the norm of)7 2290 4 1296 1680 t (de\256ned as)1 402 1 1296 1850 t 7 R f (1)1723 1920 w 7 S f (\243)1769 1920 w 7 I f (j)1830 1920 w 7 S f (\243)1861 1920 w 7 I f (n)1911 1920 w 10 R f (max)1748 1850 w 7 I f (i)1953 1950 w 7 S f (=)1989 1950 w 7 R f (1)2039 1950 w 15 S f (S)1969 1880 w 7 I f (n)1996 1750 w 10 S f (\357)2082 1850 w 10 I f (a)2139 1850 w 7 I f (i j)1 45 1 2200 1870 t 10 S f (\357)2285 1850 w 9 B f (Type:)720 2290 w 10 R f (Real function)1 541 1 1296 2290 t 9 B f (Usage:)720 2650 w 10 S f (<)1296 2650 w 10 R f (answer)1351 2650 w 10 S f (>)1633 2650 w 10 R f (= GENM \(N, A, IA\))4 815 1 1713 2650 t (N)1440 2890 w 14 S f (\256)1980 2910 w 10 R f (the number of rows in A)5 979 1 2196 2890 t (A)1440 3130 w 14 S f (\256)1980 3150 w 10 R f ( in the calling program, where IA)6 1351( KA\))1 202( dimensioned \(IA,)2 730(the array,)1 378 4 2196 3130 t 10 S f (\263)4885 3130 w 10 R f (N)4968 3130 w (and KA)1 313 1 2196 3250 t 10 S f (\263)2534 3250 w 10 R f (N, containing the N)3 788 1 2614 3250 t 10 S f (\264)3427 3250 w 10 R f (N coef\256cient matrix)2 804 1 3507 3250 t 14 S f (\254)1980 3510 w 10 R f (the LU decomposition of A \(see)5 1284 1 2196 3490 t 8 B f (Note 2)1 219 1 3505 3490 t 10 R f (\))3724 3490 w (IA)1440 3730 w 14 S f (\256)1980 3750 w 10 R f (the row \(leading\) dimension of A, as dimensioned in)8 2106 1 2196 3730 t (the calling program)2 782 1 2196 3850 t 10 S f (<)1440 4140 w 10 R f (answer)1495 4140 w 10 S f (>)1777 4140 w 14 S f (\254)1980 4160 w 7 R f (1)2196 4210 w 7 S f (\243)2242 4210 w 7 I f (j)2292 4210 w 7 S f (\243)2317 4210 w 7 I f (n)2361 4210 w 10 R f (max)2210 4140 w 7 I f (i)2404 4240 w 7 S f (=)2440 4240 w 7 R f (1)2490 4240 w 15 S f (S)2420 4170 w 7 I f (n)2447 4040 w 10 S f (\357)2533 4140 w 10 I f (a)2590 4140 w 7 I f (i j)1 45 1 2651 4160 t 10 S f (\357)2736 4140 w 9 B f (Error situations:)1 648 1 720 4580 t 10 R f (\(All errors in this subprogram are fatal \320)7 1666 1 1512 4580 t (see)1512 4700 w 10 I f (Error Handling)1 631 1 1664 4700 t 10 R f (, Framework Chapter\))2 884 1 2295 4700 t (Number Error)1 1506 1 1800 4940 t (1 N)1 1108 1 1944 5180 t 10 S f (<)3077 5180 w 10 R f (1)3157 5180 w (2 IA)1 1141 1 1944 5420 t 10 S f (<)3110 5420 w 10 R f (N)3190 5420 w 9 B f (Double precision version:)2 981 1 720 5780 t 10 R f (DGENM with A and DGENM declared double precision)7 2276 1 1776 5780 t 9 B f (Complex version:)1 678 1 720 6140 t 10 R f (CGENM with A declared complex)4 1392 1 1448 6140 t 9 B f (Storage:)720 6500 w 10 R f (None)1296 6500 w (Linear Algebra)1 606 1 4794 7500 t (\320 29 \320)2 350 1 2705 7620 t (GENM)5106 7740 w cleartomark showpage saveobj restore %%EndPage: 29 29 %%Page: 30 30 /saveobj save def mark 30 pagesetup 10 R f (-- --)1 5544 1 0 120 t (Linear Algebra)1 606 1 360 600 t 8 R f (PORT)4903 600 w 10 R f (library)5134 600 w ( 11, 1993)2 375(GENM February)1 4665 2 360 840 t 9 B f (Time:)720 1320 w 10 R f (N)1296 1320 w 7 R f (2)1373 1280 w 10 R f (additions)1441 1320 w (N comparisons)1 602 1 1296 1440 t 9 B f (See also:)1 333 1 720 1800 t 10 R f (GEDC, GELU, GELE, GESS, GECE)4 1499 1 1296 1800 t 9 B f (Author:)720 2160 w 10 R f (Linda Kaufman)1 629 1 1296 2160 t 9 B f (Example:)720 2520 w 10 R f (The subroutines in the PORT library for solving)7 1938 1 1296 2520 t 10 I f (Ax)3261 2520 w 10 S f (=)3415 2520 w 10 I f (b)3519 2520 w 10 R f ( so-)1 150(are designed to return computed)4 1294 2 3596 2520 t (lutions)1296 2640 w 10 I f (x)1594 2640 w 10 R f (such that the residual)3 846 1 1663 2640 t 10 I f (r)2534 2640 w 10 S f (=)2622 2640 w 10 I f (Ax)2726 2640 w 10 S f (-)2880 2640 w 10 I f (b)2984 2640 w 10 R f (satis\256es)3059 2640 w 10 S f (\357 \357)1 106 1 2809 2990 t 10 I f (A)2923 2990 w 10 S f ( \357)1 57( \357)1 90(\357 \357)1 106 3 2992 2990 t 10 I f (x)3253 2990 w 10 S f (\357 \357)1 106 1 3305 2990 t (\357 \357)1 106 1 2977 2860 t 10 I f (r)3091 2860 w 10 S f (\357 \357)1 106 1 3138 2860 t 10 S1 f (_ ____________)1 632 1 2795 2890 t 10 S f (\243 e)1 107 1 3445 2920 t 10 R f (where)1296 3210 w 10 S f (e)1569 3210 w 10 R f ( if A is ill-conditioned, then)5 1141( this example we show that)5 1110( In)1 138(is the machine precision.)3 1008 4 1643 3210 t ( \(1.1\))1 219(the computed solution need not be very close to the true solution even though equation)14 3525 2 1296 3330 t ( matrix)1 288( The)1 207( compute the left-hand side of \(1.1\).)6 1448(is satis\256ed. The subroutine GENM is used to)7 1801 4 1296 3450 t (in this example is given by)5 1075 1 1296 3570 t 10 I f (a)2807 3895 w 7 I f (i j)1 45 1 2868 3915 t 10 S f (=)2970 3895 w (\354)3082 3808 w (\355)3082 3908 w (\356)3082 4008 w 10 I f (i)3172 3975 w 10 S f (-)3224 3975 w 10 I f (j)3295 3975 w 10 S f (+)3372 3975 w 10 R f (1 for)1 207 1 3476 3975 t 10 I f (i)3724 3975 w 10 S f (\263)3760 3975 w 10 I f (j)3831 3975 w (j)3172 3835 w 10 S f (-)3216 3835 w 10 I f (i)3287 3835 w 10 R f (for)3422 3835 w 10 I f (i)3579 3835 w 10 S f (<)3631 3835 w 10 I f (j)3702 3835 w 10 R f (and the true solution is)4 955 1 1296 4240 t 10 I f (x)2287 4240 w 7 I f (i)2342 4260 w 10 S f (=)2386 4240 w 10 I f (i)2457 4240 w 10 R f ( right hand side is generated using GEML and the com-)10 2339(. The)1 216 2 2485 4240 t ( GELE. The function SAMAX is used to compute the 1-norm)10 2467(puted solution is obtained using)4 1277 2 1296 4360 t (of a vector; i.e.)3 601 1 1296 4480 t 7 R f (1)1922 4550 w 7 S f (\243)1962 4550 w 7 I f (i)2006 4550 w 7 S f (\243)2031 4550 w 7 I f (n)2075 4550 w 10 R f (max)1930 4480 w 10 S f (\357)2118 4480 w 10 I f (x)2175 4480 w 7 I f (i)2230 4500 w 10 S f (\357)2290 4480 w 8 CW f (INTEGER I, J, L, N, IA, IWRITE, I1MACH)7 1824 1 2040 4870 t (REAL A\(50, 50\), AA\(50, 50\), B\(50\), X\(50\))6 1920 1 2040 4970 t (REAL RELERR, RELRES, XNORM, RNORM, ERR, R\(50\))6 2160 1 2040 5070 t (REAL GENM, SAMAX)2 768 1 2040 5170 t (IA = 50)2 336 1 2040 5270 t (C)1656 5370 w (C GENERATE MATRIX)2 816 1 1656 5470 t (C)1656 5570 w (N=50)2040 5670 w (DO 20 I=1,N)2 528 1 2040 5770 t (DO 10 J=I,N)2 528 1 2184 5870 t (A\(I,J\)=J-I)2328 5970 w (A\(J,I\)=J-I + 1)2 672 1 2328 6070 t (AA\(I,J\)=A\(I,J\))2328 6170 w (AA\(J,I\)=A\(J,I\))2328 6270 w (10 CONTINUE)1 864 1 1704 6370 t (B\(I\)=I)2184 6470 w (20 CONTINUE)1 720 1 1704 6570 t (C)1656 6670 w (C GENERATE RIGHT HAND SIDE)4 1248 1 1656 6770 t (C)1656 6870 w 10 R f (Linear Algebra)1 606 1 360 7530 t (\320 30 \320)2 350 1 2705 7650 t (GENM)360 7770 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 ( GENM)1 4305(February 11, 1993)2 735 2 360 840 t 8 CW f (CALL GEML\(N,A,IA,B,X\))1 1008 1 2040 1300 t (C)1656 1400 w (C MAKE COPY OF RIGHT HAND SIDE)6 1440 1 1656 1500 t (C)1656 1600 w (CALL MOVEFR\(N,X,B\))1 864 1 2040 1700 t (C)1656 1800 w (C SOLVE THE SYSTEM)3 864 1 1656 1900 t (C)1656 2000 w (CALL GELE\(N,A,IA,B,N,1\))1 1104 1 2040 2100 t (C)1656 2200 w (C COMPUTE THE RELATIVE ERROR AND THE RELATIVE RESIDUAL)8 2592 1 1656 2300 t (C)1656 2400 w (CALL GEML\(N,AA,IA,B,R\))1 1056 1 2040 2500 t (ERR=0.0)2040 2600 w (DO 30 I=1,N)2 528 1 2040 2700 t (ERR=AMAX1\(ERR,ABS\(B\(I\)-FLOAT\(I\)\)\))2184 2800 w (R\(I\)=R\(I\)-X\(I\))2184 2900 w (30 CONTINUE)1 720 1 1704 3000 t (XNORM=SAMAX\(N,X,1\))2040 3100 w (RNORM=SAMAX\(N,R,1\))2040 3200 w (RELERR=ERR/XNORM)2040 3300 w (RELRES=RNORM/\(XNORM*GENM\(N,AA,IA\)\))2040 3400 w (IWRITE=I1MACH\(2\))2040 3500 w (WRITE\(IWRITE,31\)RELERR,RELRES)2040 3600 w ( RELATIVE ERROR=,E15.5,19H RELATIVE RESIDUAL=,)4 2208(31 FORMAT\(16H)1 816 2 1704 3700 t (1 E15.5\))1 480 1 1896 3800 t (STOP)2040 3900 w (END)2040 4000 w 10 R f ( on the Honeywell 6000 machine at Bell Laboratories,)8 2174(When the above program was executed)5 1570 2 1296 4340 t (the following was printed:)3 1052 1 1296 4460 t 8 CW f ( 0.22987E-10)1 720( RELATIVE RESIDUAL=)2 912( 0.13554E-06)1 720(RELATIVE ERROR=)1 720 4 1488 4900 t 10 R f ( of the matrix\(see the example in GELE\) is about 10)10 2101(The condition number)2 888 2 1296 5240 t 7 R f (5)4290 5200 w 10 R f (, and the machine)3 707 1 4333 5240 t ( computer is about 10)4 866(precision on the Honeywell)3 1104 2 1296 5360 t 7 S f (-)3277 5320 w 7 R f (8)3327 5320 w 10 R f (. Thus even in the absence of roundoff er-)8 1670 1 3370 5360 t ( a relative error of 10)5 929(ror in GEML,)2 584 2 1296 5480 t 7 S f (-)2820 5440 w 7 R f (3)2870 5440 w 10 R f ( relative error given)3 840( The)1 222(would not be surprising.)3 1023 3 2955 5480 t ( as promised, satis\256es \(1.1\) even though)6 1658(above is quite within reason. The relative residual,)7 2086 2 1296 5600 t (the problem is ill-conditioned.)3 1211 1 1296 5720 t (Linear Algebra)1 606 1 4794 7500 t (\320 31 \320)2 350 1 2705 7620 t (GENM)5106 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 ( 11, 1993)2 375(GENM February)1 4665 2 360 840 t (GESS \320 general linear system solution with condition estimation)8 2639 1 1848 1320 t 9 B f (Purpose:)720 1680 w 10 R f ( where A is a general matrix. An)7 1316( = B)2 173( AX)1 171(GESS \(GEneral System Solution\) solves the system)6 2084 4 1296 1680 t (estimate of the condition of A is provided.)7 1693 1 1296 1800 t 9 B f (Usage:)720 2160 w 10 R f (CALL GESS \(N, A, IA, B, IB, NB, COND\))8 1760 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 (A)1440 2640 w 14 S f (\256)1980 2660 w 10 R f (the array, dimensioned \(IA, KA\) in the calling program,)8 2237 1 2196 2640 t (where IA)1 373 1 2196 2760 t 10 S f (\263)2594 2760 w 10 R f (N and KA)2 410 1 2674 2760 t 10 S f (\263)3109 2760 w 10 R f (N, containing the N)3 788 1 3189 2760 t 10 S f (\264)4002 2760 w 10 R f (N coef\256cient matrix)2 804 1 4082 2760 t (A is overwritten during the solution.)5 1455 1 2196 2880 t (IA)1440 3120 w 14 S f (\256)1980 3140 w 10 R f (the row \(leading\) dimension of A, as dimensioned in the)9 2253 1 2196 3120 t (calling program)1 635 1 2196 3240 t (B)1440 3480 w 14 S f (\256)1980 3500 w 10 R f (the matrix of right-hand sides, dimensioned \(IB, KB\) in)8 2226 1 2196 3480 t (the calling program, where IB)4 1200 1 2196 3600 t 10 S f (\263)3421 3600 w 10 R f (N and KB)2 405 1 3501 3600 t 10 S f (\263)3931 3600 w 10 R f (NB)4011 3600 w 14 S f (\254)1980 3860 w 10 R f (the solution X)2 567 1 2196 3840 t (IB)1440 4080 w 14 S f (\256)1980 4100 w 10 R f (the row \(leading\) dimension of B, as dimensioned in the)9 2248 1 2196 4080 t (calling program)1 635 1 2196 4200 t (NB)1440 4440 w 14 S f (\256)1980 4460 w 10 R f (the number of right-hand sides)4 1226 1 2196 4440 t (COND)1440 4680 w 14 S f (\254)1980 4700 w 10 R f ( \(see)1 210(an estimate of the condition number of A)7 1645 2 2196 4680 t 8 B f (Note 1)1 219 1 4076 4680 t 10 R f (\))4295 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 ( 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 5760 t ( however, the user is)4 863( Ordinarily,)1 501(use the routine GELE, which is a little faster than GESS.)10 2380 3 1296 5880 t (strongly urged to choose GESS, and to follow it by a test of the condition estimate.)15 3318 1 1296 6000 t (Linear Algebra)1 606 1 360 7500 t (\320 32 \320)2 350 1 2705 7620 t (GESS)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 ( GESS)1 4305(February 11, 1993)2 735 2 360 840 t 9 B f (Note 2:)1 278 1 720 1320 t 10 R f ( coef\256cient matrix, but differ-)4 1201(Users who wish to solve a sequence of problems with the same)11 2543 2 1296 1320 t (ent right-hand sides)2 796 1 1296 1440 t 10 I f (not all known in advance,)4 1046 1 2121 1440 t 10 R f ( subpro-)1 333(should not use GESS, but should call)6 1511 2 3196 1440 t ( is called once to get)5 858( GECE)1 319( the example of GEDC.\))4 1005( \(See)1 234(grams GECE, GEFS and GEBS.)4 1328 5 1296 1560 t ( to this chapter\) and then the pair, GEFS \(forward)9 1994(the LU decomposition \(see the introduction)5 1750 2 1296 1680 t (solve\) and GEBS \(back solve\), is called for each new right-hand side.)11 2780 1 1296 1800 t 9 B f (Error situations:)1 653 1 720 2160 t 10 R f ( see)1 183( \320)1 125( those errors marked with an asterisk)6 1505(*\(The user can elect to `recover' from)6 1546 4 1517 2160 t 10 I f (Er-)4907 2160 w (ror Handling)1 531 1 1512 2280 t 10 R f (, Framework Chapter\))2 884 1 2043 2280 t (Number Error)1 1506 1 1800 2520 t (1 N)1 1008 1 1944 2760 t 10 S f (<)2977 2760 w 10 R f (1)3057 2760 w (2 IA)1 1041 1 1944 3000 t 10 S f (<)3010 3000 w 10 R f (N)3090 3000 w (3 IB)1 1036 1 1944 3240 t 10 S f (<)3005 3240 w 10 R f (N)3085 3240 w (4 NB)1 1075 1 1944 3480 t 10 S f (<)3044 3480 w 10 R f (1)3124 3480 w ( matrix whose rank is at least k)7 1240( singular)1 952(10 + k*)2 306 3 1944 3720 t 9 B f (Double-precision version:)1 988 1 720 4080 t 10 R f (DGESS with A, B, and COND declared double precision)8 2286 1 1783 4080 t 9 B f (Complex version:)1 678 1 720 4440 t 10 R f (CGESS with A and B declared complex)6 1604 1 1498 4440 t 9 B f (Storage:)720 4800 w 10 R f (N integer locations and)3 929 1 1296 4800 t ( for DGESS, complex for CGESS\) locations of scratch storage in the)11 2774(N real \(double precision)3 970 2 1296 4920 t (dynamic storage stack)2 887 1 1296 5040 t 9 B f (Time:)720 5450 w 10 R f (3)1356 5520 w (N)1321 5390 w 7 R f (3)1398 5350 w 10 S1 f (___)1306 5420 w 10 S f (+)1515 5450 w 10 R f (N)1619 5450 w 7 R f (2)1696 5410 w 10 S f (\264)1780 5450 w 10 R f (\()1876 5450 w (2)1942 5520 w (9)1942 5390 w 10 S1 f (_ _)1 80 1 1927 5420 t 10 S f (+)2066 5450 w 10 R f (NB \))1 180 1 2170 5450 t 10 S f (+)2407 5450 w 10 R f (N)2511 5450 w 10 S f (\264)2624 5450 w 10 R f (\()2720 5450 w (6)2811 5520 w (19)2786 5390 w 10 S1 f (_ __)1 130 1 2771 5420 t 10 S f (+)2960 5450 w 10 R f ( additions)1 392(NB \))1 180 2 3064 5450 t (3)1356 5800 w (N)1321 5670 w 7 R f (3)1398 5630 w 10 S1 f (___)1306 5700 w 10 S f (+)1515 5730 w 10 R f (N)1619 5730 w 7 R f (2)1696 5690 w 10 S f (\264)1780 5730 w 10 R f (\()1876 5730 w (2)1942 5800 w (5)1942 5670 w 10 S1 f (_ _)1 80 1 1927 5700 t 10 S f (+)2033 5730 w 10 R f (NB \))1 180 1 2104 5730 t 10 S f (+)2341 5730 w 10 R f (N)2445 5730 w 10 S f (\264)2558 5730 w 10 R f (\()2654 5730 w (6)2720 5800 w (7)2720 5670 w 10 S1 f (_ _)1 80 1 2705 5700 t 10 S f (+)2811 5730 w 10 R f ( multiplications)1 626(NB \))1 180 2 2882 5730 t (2)1356 6080 w (N)1321 5950 w 7 R f (2)1398 5910 w 10 S1 f (___)1306 5980 w 10 S f (+)1515 6010 w 10 R f (N)1619 6010 w 10 S f (\264)1732 6010 w 10 R f (\()1828 6010 w (2)1894 6080 w (3)1894 5950 w 10 S1 f (_ _)1 80 1 1879 5980 t 10 S f (+)2018 6010 w 10 R f ( divisions)1 387(NB \))1 180 2 2122 6010 t 9 B f (Method:)720 6420 w 10 R f (Gaussian elimination with partial pivoting.)4 1714 1 1296 6420 t (See the reference below for the method used to estimate the condition number.)12 3141 1 1296 6540 t (GESS calls GECE, GEFS, and GEBS.)5 1534 1 1296 6660 t (Linear Algebra)1 606 1 4794 7500 t (\320 33 \320)2 350 1 2705 7620 t (GESS)5155 7740 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 ( 11, 1993)2 375(GESS February)1 4665 2 360 840 t 9 B f (See also:)1 333 1 720 1320 t 10 R f (GEBS, GECE, GEDC, GEFS, GELE, GELU)5 1805 1 1296 1320 t 9 B f (Authors:)720 1680 w 10 R f (Linda Kaufman and Doris Ryan)4 1281 1 1296 1680 t 9 B f (Reference:)720 2040 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 2040 t (tion number,)1 511 1 1296 2160 t 10 I f (SIAM J. Numer. Anal. 16)4 1007 1 1857 2160 t 10 R f (\(1979\), 368-375.)1 674 1 2889 2160 t 9 B f (Example:)720 2520 w 10 R f (The following program solves a 5)5 1351 1 1296 2520 t 10 S f (\264)2672 2520 w 10 R f (5 system with two right-hand sides)5 1397 1 2752 2520 t 8 CW f (INTEGER N, IREAD, I1MACH, I, NB, IWRITE, J)7 2016 1 1416 2860 t (REAL A\(5,5\), B\(5,2\), COND)3 1200 1 1416 2960 t (N=5)1416 3060 w (IREAD=I1MACH\(1\))1416 3160 w (C)1080 3260 w (DO 10 I=1,N)2 528 1 1416 3360 t (READ\(IREAD,1\) \(A\(I,J\),J=1,N\))1 1344 1 1608 3460 t (1 FORMAT\(1X,5F10.0\))1 1152 1 1272 3560 t (10 CONTINUE)1 576 1 1224 3660 t (C)1080 3760 w (NB=2)1416 3860 w (DO 20 I=1,N)2 528 1 1416 3960 t (READ\(IREAD,11\) \(B\(I,J\),J=1,NB\))1 1440 1 1608 4060 t (11 FORMAT\(1X,2F10.3\))1 1200 1 1224 4160 t (20 CONTINUE)1 576 1 1224 4260 t (C)1080 4360 w ( CALLING GESS)2 624( BY)1 192(C SOLVE AX = B)4 672 3 1080 4460 t (C)1080 4560 w (CALL GESS\(N,A,N,B,N,NB,COND\))1 1344 1 1416 4660 t (IWRITE=I1MACH\(2\))1416 4760 w (WRITE\(IWRITE,21\) COND)1 1008 1 1416 4860 t ( AN ESTIMATE OF THE CONDITION NUMBER OF THE MATRIX =,)10 2544(21 FORMAT\(52H)1 672 2 1224 4960 t (1 E14.7\))1 816 1 1320 5060 t (C)1080 5160 w (WRITE\(IWRITE,22\))1416 5260 w ( THE COMPUTED SOLUTION X IS,//\))5 1488(22 FORMAT\(27H)1 672 2 1224 5360 t (DO 30 I=1,N)2 528 1 1416 5460 t (WRITE\(IWRITE,23\) \(B\(I,J\),J=1,NB\))1 1536 1 1608 5560 t (23 FORMAT\(1H,5F20.7\))1 1200 1 1224 5660 t (30 CONTINUE)1 576 1 1224 5760 t (C)1080 5860 w (STOP)1416 5960 w (END)1416 6060 w 10 R f (Linear Algebra)1 606 1 360 7500 t (\320 34 \320)2 350 1 2705 7620 t (GESS)360 7740 w cleartomark showpage saveobj restore %%EndPage: 34 34 %%Page: 35 35 /saveobj save def mark 35 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 ( GESS)1 4305(February 11, 1993)2 735 2 360 840 t (For the input matrix given by:)5 1203 1 720 1320 t (1.)2212 1560 w 10 S f (-)2437 1560 w 10 R f (2. 3. 7.)2 735 1 2492 1560 t 10 S f (-)3472 1560 w 10 R f (9.)3527 1560 w 10 S f (-)2157 1680 w 10 R f (2. 8.)1 355 1 2212 1680 t 10 S f (-)2767 1680 w 10 R f ( 50.)1 375(6. 9.)1 405 2 2822 1680 t (11.)2162 1800 w 10 S f (-)2437 1800 w 10 R f (6. 18.)1 405 1 2492 1800 t 10 S f (-)3047 1800 w 10 R f (15.)3102 1800 w 10 S f (-)3422 1800 w 10 R f (18.)3477 1800 w (7. 2.)1 355 1 2212 1920 t 10 S f (-)2717 1920 w 10 R f ( 173.)1 375(15. 273.)1 455 2 2772 1920 t 10 S f (-)2157 2040 w 10 R f (9. 50.)1 355 1 2212 2040 t 10 S f (-)2717 2040 w 10 R f ( 1667.)1 375(18. 6.)1 455 2 2772 2040 t (and the following right-hand sides:)4 1399 1 720 2280 t (30. 29.419)1 705 1 2605 2520 t 10 S f (-)2500 2640 w 10 R f (191.)2555 2640 w 10 S f (-)2930 2640 w 10 R f (190.994)2985 2640 w (133. 133.072)1 755 1 2555 2760 t 10 S f (-)2500 2880 w 10 R f (986.)2555 2880 w 10 S f (-)2930 2880 w 10 R f (985.775)2985 2880 w 10 S f (-)2450 3000 w 10 R f (6496.)2505 3000 w 10 S f (-)2880 3000 w 10 R f (6495.553)2935 3000 w (the following results were obtained on the Honeywell 6000 computer at Bell Labs:)12 3307 1 720 3240 t 8 CW f (AN ESTIMATE OF THE CONDITION NUMBER OF THE MATRIX = 0.2759414E 04)11 3120 1 1128 3580 t (THE COMPUTED SOLUTION X IS)4 1248 1 1128 3680 t 10 R f (2.0000004 2.4800003)1 1055 1 2380 4040 t (4.9999970 4.8709986)1 1055 1 2380 4160 t (2.9999988 2.6439993)1 1055 1 2380 4280 t 10 S f (-)2325 4400 w 10 R f (1.0000001)2380 4400 w 10 S f (-)2955 4400 w 10 R f (1.0320001)3010 4400 w 10 S f (-)2325 4520 w 10 R f (3.9999999)2380 4520 w 10 S f (-)2955 4520 w 10 R f (3.9970000)3010 4520 w (The true solution to this problem is:)6 1434 1 720 4880 t (2. 2.48)1 455 1 2655 5120 t (5. 4.871)1 505 1 2655 5240 t (3. 2.644)1 505 1 2655 5360 t 10 S f (-)2600 5480 w 10 R f (1.)2655 5480 w 10 S f (-)2880 5480 w 10 R f (1.032)2935 5480 w 10 S f (-)2600 5600 w 10 R f (4.)2655 5600 w 10 S f (-)2880 5600 w 10 R f (3.997)2935 5600 w ( Fur-)1 226(Notice that a seemingly slight change in the right-hand side causes the solution to change noticeably.)15 4094 2 720 5840 t ( about 2)2 352(thermore, the relative error in the solution is)7 1876 2 720 5960 t 10 S f (\264)2989 5960 w 10 R f (10)3085 5960 w 7 S f (-)3196 5920 w 7 R f (7)3246 5920 w 10 R f ( the Honeywell computer, which has)5 1539(. On)1 212 2 3289 5960 t ( A)1 127( 1.5 decimal digits.)3 779(about 8 decimal digits for single-precision numbers, this represents the loss of about)12 3414 3 720 6080 t (loss of up to 2)4 567 1 720 6200 t 10 S f (\264)1328 6200 w 10 R f (10)1424 6200 w 7 S f (-)1535 6160 w 7 R f (5)1585 6160 w 10 R f (could be expected in light of the analysis given below.)9 2175 1 1653 6200 t (Linear Algebra)1 606 1 4794 7500 t (\320 35 \320)2 350 1 2705 7620 t (GESS)5155 7740 w cleartomark showpage saveobj restore %%EndPage: 35 35 %%Page: 36 36 /saveobj save def mark 36 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(GESS February)1 4665 2 360 840 t (Let)720 1320 w 10 S f (D)903 1320 w 10 I f (b)972 1320 w 10 R f (represent a perturbation in the right-hand side of a linear system.)10 2581 1 1047 1320 t (If)720 1440 w 10 I f (Ax)811 1440 w 10 S f (=)965 1440 w 10 I f (b)1069 1440 w 10 R f (then)1144 1440 w 10 I f (A)1944 1680 w 10 R f (\()2013 1680 w 10 I f (x)2054 1680 w 10 S f (+ D)1 165 1 2147 1680 t 10 I f (x)2320 1680 w 10 R f (\))2372 1680 w 10 S f (=)2462 1680 w 10 I f (b)2566 1680 w 10 S f (+ D)1 165 1 2665 1680 t 10 I f (b)2838 1680 w 10 R f (where)720 1920 w 10 S f (| |)1 48 1 2003 2315 t 10 I f (x)2059 2315 w 10 S f (| |)1 48 1 2111 2315 t (| | D)2 117 1 1969 2185 t 10 I f (x)2094 2185 w 10 S f (| |)1 48 1 2146 2185 t 10 S1 f (_ _____)1 255 1 1954 2215 t 10 S f (\243)2260 2245 w 10 I f (K)2356 2245 w 10 R f (\()2431 2245 w 10 I f (A)2472 2245 w 10 R f (\))2541 2245 w 10 S f (\354)2598 2158 w (\357)2598 2258 w (\356)2598 2358 w (| |)1 48 1 2706 2315 t 10 I f (b)2762 2315 w 10 S f (| |)1 48 1 2820 2315 t (| | D)2 117 1 2672 2185 t 10 I f (b)2797 2185 w 10 S f (| |)1 48 1 2855 2185 t 10 S1 f (_ _____)1 261 1 2657 2215 t 10 S f (\374)2928 2158 w (\357)2928 2258 w (\376)2928 2358 w 10 R f (where)720 2590 w 10 I f (K)1031 2590 w 10 R f (\()1106 2590 w 10 I f (A)1147 2590 w 10 R f (\) is the condition number of)5 1333 1 1216 2590 t 10 I f (A)2618 2590 w 10 R f (,)2679 2590 w 10 I f (K)2773 2590 w 10 R f (\()2848 2590 w 10 I f (A)2889 2590 w 10 R f (\))2958 2590 w 10 S f ( |)1 28(= |)1 124 2 3048 2590 t 10 I f (A)3208 2590 w 10 S f ( |)1 28( |)1 61(| |)1 48 3 3277 2590 t 10 I f (A)3422 2590 w 7 S f (-)3494 2550 w 7 R f (1)3538 2550 w 10 S f (| |)1 48 1 3589 2590 t 10 R f (and)3706 2590 w 10 S f (| |)1 48 1 3919 2590 t 10 R f (.)3975 2560 w 10 S f (| |)1 48 1 4008 2590 t 10 R f (, is some norm, e.g.,)4 984 1 4056 2590 t 10 S f (| |)1 48 1 720 2760 t 10 I f (x)776 2760 w 10 S f (| |)1 48 1 828 2760 t 7 R f (1)887 2780 w 10 S f (=)979 2760 w 7 I f (i)1083 2860 w 7 S f (=)1119 2860 w 7 R f (1)1169 2860 w 15 S f (S)1099 2790 w 7 I f (n)1126 2660 w 10 S f (|)1212 2760 w 10 I f (x)1240 2760 w 7 I f (i)1295 2780 w 10 S f (|)1331 2760 w 10 R f (if)1376 2760 w 10 I f (x)1462 2760 w 10 R f (is a vector.)2 435 1 1531 2760 t ( accurate answer to a slightly)5 1176(The methods used in our linear equation package are guaranteed to provide an)12 3144 2 720 3080 t ( the correct answer to a problem where)7 1728( we assume that our method produces)6 1659( If)1 142(perturbed problem.)1 791 4 720 3200 t 10 S f (| | D)2 117 1 720 3320 t 10 I f (b)845 3320 w 10 S f ( | |)2 56( \243 e)2 181(| |)1 48 3 903 3320 t 10 I f (b)1196 3320 w 10 S f (| |)1 48 1 1254 3320 t 10 R f (, where)1 294 1 1328 3320 t 10 S f (e)1648 3320 w 10 R f ( where)1 270(is the machine precision, then on the Honeywell 6000)8 2158 2 1718 3320 t 10 S f (e)4173 3320 w 10 R f (is about 10)2 443 1 4244 3320 t 7 S f (-)4692 3280 w 7 R f (8)4736 3280 w 10 R f (, a rel-)2 261 1 4779 3320 t (ative error for the above example of 2)7 1509 1 720 3440 t 10 S f (\264)2270 3440 w 10 R f (10)2366 3440 w 7 S f (-)2471 3400 w 7 R f (5)2515 3400 w 10 R f (would not be surprising.)3 972 1 2583 3440 t (In our example one may consider the \256rst column of)9 2247 1 720 3680 t 10 I f (B)3009 3680 w 10 R f (as)3113 3680 w 10 I f (b)3239 3680 w 10 R f (in \(1.1\), and the second column of)6 1478 1 3332 3680 t 10 I f (B)4853 3680 w 10 R f (as)4957 3680 w 10 I f (b)720 3800 w 10 S f (+ D)1 165 1 819 3800 t 10 I f (b)992 3800 w 10 R f (, so that)2 316 1 1042 3800 t 10 S f (| | D)2 117 1 1384 3800 t 10 I f (b)1509 3800 w 10 S f (| |)1 48 1 1567 3800 t 10 I f (/)1623 3800 w 10 S f (| |)1 48 1 1659 3800 t 10 I f (b)1715 3800 w 10 S f (| |)1 48 1 1773 3800 t 10 R f (is approximately .00015 using the)4 1362 1 1847 3800 t 10 S f (| |)1 48 1 3234 3800 t 10 R f (.)3290 3770 w 10 S f (| |)1 48 1 3323 3800 t 7 R f (1)3382 3820 w 10 R f ( we look at the second solution)6 1238(norm. If)1 352 2 3450 3800 t (as)720 3920 w 10 I f (x)830 3920 w 10 S f (+ D)1 165 1 923 3920 t 10 I f (x)1096 3920 w 10 R f ( \256rst solution as)3 646(in \(1.1\) and the)3 616 2 1167 3920 t 10 I f (x)2457 3920 w 10 R f (, then)1 225 1 2501 3920 t 10 S f (| | D)2 117 1 2754 3920 t 10 I f (x)2879 3920 w 10 S f (| |)1 48 1 2931 3920 t 10 I f (/)2987 3920 w 10 S f (| |)1 48 1 3023 3920 t 10 I f (x)3079 3920 w 10 S f (| |)1 48 1 3131 3920 t 10 R f ( equation \(1.1\) in-)3 730( Thus)1 253(is approximately .07.)2 850 3 3207 3920 t ( condition estimate, which at \256rst appeared to be)8 1995(dicates that the condition number is at least 400, and the)10 2325 2 720 4040 t (conservative, was in fact quite realistic.)5 1577 1 720 4160 t (Linear Algebra)1 606 1 360 7500 t (\320 36 \320)2 350 1 2705 7620 t (GESS)360 7740 w (-- --)1 5544 1 0 8030 t cleartomark showpage saveobj restore %%EndPage: 36 36 %%Trailer done %%Pages: 36 %%DocumentFonts: Courier Times-Bold Times-Italic Times-Roman Times-Roman Symbol