Higham - Accuracy and Stability of Numerical Algorithms (523152), страница 91
Текст из файла (страница 91)
For example, SLAS2 computes the singular values of a 2 x 2 triangular matrix, SLAEV2 computes the eigenvalues andeigenvectors of a 2 x 2 symmetric matrix, and SLASY2 solves a Sylvester equation AX – XB = C where A and B have order 1 or 2. These routines needto be reliable and robust because they are called repeatedly by higher-levelsolvers. Because LAPACK makes minimal assumptions about the underlyingfloating point arithmetic, the 2 x 2 routines are nontrivial to write and aresurprisingly long: the counts of executable statements are 39 for SLAS2, 61 for25.7 P ORTABILITY501SLAEV2, and 207 for SLASY2.
If the availability of extended precision arithmetic (possibly simulated using the working precision) or IEEE arithmeticcan be assumed, the codes can be simplified significantly. Complicated andless efficient code for these 2 x 2 problem solvers is a price to be paid forportability across a wide range of floating point arithmetics.25.7.3. Numerical ConstantsHow should numerical constants be specified, e.g., for a Gauss quadrature ruleor a Runge–Kutta method? Some compilers limit the number of significantdigits that can be specified in a DATA statement or an assignment statement.One possibility is to use rational approximations, but again the number ofdigits required to cater for all precision may be too large for some compilers.Another scheme is to store a carefully designed table of leading and trailingparts of constants. A constant such as π is best computed via a statement suchas pi = 4.0*atan(1.0), if you trust your language’s mathematics library.25.7.4.
Models of Floating Point ArithmeticThe model (2.4) of floating point arithmetic contains only one parameter, theunit roundoff. Other more detailed models of floating point arithmetic havebeen proposed, with the intention of promoting the development of portablemathematical software.
A program developed under a model and provedcorrect for that model has the attraction that it necessarily performs correctlyon any machine for which the model is valid.The first detailed model was proposed by van Wijngaarden [1046, 1966];it contained 32 axioms and was unsuccessful because it was “mathematicallyintimidating” and did not cover the CDC 600, an important high-performancecomputer of that time [635, 1991]. A more recent model is that of Brown [150,1981]. Brown’s model contains four parameters: the base β, the precision t,and the minimum and maximum exponents emin and emax together with anumber of axioms describing the behaviour of the floating point arithmetic.Aberrant arithmetics are accommodated in the model by penalizing theirparameters (for example, reducing t by 1 from its actual machine value ifthere is no guard digit).
Brown builds a substantial theory on the modeland gives an example program to compute the 2-norm of a vector, which isaccompanied by a correctness proof for the model.Brown’s model is intended to cater to diverse computer arithmetics. Infact, “any behaviour permitted by the axioms of the model is actually exhibited by at least one commercially important computer” [150, 1981, p.
457].This broadness is a weakness. It means that there are important properties shared by many (but not all) machines that the model cannot reflect,such as the exactness of fl(x – y) described in Theorem 2.4. As Kahan [630,502S OFTWARE I SSUESINFLOATING P OINT A RITHMETIC1981] notes, “Programming for the [IEEE] standard is like programming forone of a family of well-known machines, whereas programming for a model islike programming for a horde of obscure and ill-understood machines all atonce.” Although Brown’s model was used in the Ada language to provide adetailed specification of floating point arithmetic, the model is still somewhatcontroversial.Wichmann [1077, 1989] gives a formal specification of floating point arithmetic in the VDM specification language based on a modification of Brown’smodel.The most recent model is the Language Independent Arithmetic (LIA1) [702, 1993].
The LIA-1 specifies floating point arithmetic far more tightlythan Brown’s model. It, too, is controversial; an explanation of flaws in anearlier version (know then as the Language Compatible Arithmetic Standard)was published by Kahan [635, 1991].For a more detailed critique of models of floating point arithmetic seePriest [844, 1992].25.8.
Avoiding Underflow and OverflowA classic example showing how care is needed to avoid underflow and overflowis the evaluation of a vector 2-norm,For about half ofall machine numbers x, x2 either underflows or overflows! Overflow is avoidedby the following algorithm:s =0for i = l:ns = s + (x(i)/t)2endThe trouble with this algorithm is that it requires n divisions and two passesover the data (the first to evaluateso it is slow. (It also involvesmore rounding errors than the unscaled evaluation, which could be obviatedby scaling by a power of the machine base.) Blue [128, 1978] develops a onepass algorithm that avoids both overflow and underflow and requires betweenO and n divisions, depending on the data, and he gives an error analysis toshow that the algorithm produces an accurate answer. The idea behind thealgorithm is simple.
The sum of squares is collected in three accumulators,one each for small, medium, and large numbers. After the summation, variouslogical tests are used to decide how to combine the accumulators to obtainthe final answer.The original, portable implementation of the level-l BLAS routine xNRM2(listed in [307, 1979]) was written by C. L. Lawson in 1978 and, according to25.8 A VOIDING U NDERFLOWANDO VERFLOW503Lawson, Hanson, Kincaid, and Krogh [694, 1979], makes use of Blue’s ideas.However, xNRM2 is not a direct coding of Blue’s algorithm and is extremelydifficult to understand—a classic piece of “Fortrans spaghetti”! Nevertheless,the routine clearly works and is reliable, because it has been very widely usedwithout any reported problems.
Lawson’s version of xNRM2 has now been superseded in the LAPACK distribution by a concise and elegant version by S.J. Harnmarling, which implements a one-pass algorithm; see Problem 25.5.A special case of the vector 2-norm is the Pythagorean sum,which occurs in many numerical computations. One way to compute itis by a beautiful cubically convergent, square-root-free iteration devised byMoler and Morrison [773, 1983 ; see Problem 25.6. LAPACK has a routinexLAPY2 that computesand another routine xLAPY3 that computesboth routines avoid overflow by using the algorithm listed atthe start of this section with n = 2 or 3.Pythagorean sums arise in computing the l-norm of a complex vector:The level-l BLAS routine xCASUM does not compute the l-norm, but the moreeasily computed quantityThe reason given by the BLAS developers is that it was assumed that userswould expect xCASUM to compute a less expensive measure of size than the 2norm [694, 1979].
This reasoning is sound, but many users have been confusednot to receive the true l-norm. See Problem 6.16 for more on this pseudo 1norm.Another example where underflow and overflow can create problems is incomplex division. If we use the formulathen overflow will occur whenever cord exceeds the square root of the overflowthreshold, even if the quotient (a + ib)/(c + id) does not overflow. Certain Crayand NEC machines implement complex division in this way [290, 1992]; onthese machines the exponent range is effectively halved from the point of viewof the writer of robust software.
Smith [929, 1962] suggests a simple way toavoid overflow: if |c| > |d| use the following formula, obtained by multiplyingthe numerators and denominators by c–1 ,(25.1)504S OFTWARE I SSUESINFLOATING P OINT A RITHMETICand if |d| > |c| use the analogous formula involving d–1. Stewart [947, 1985]points out that underflow is still possible in these formulae, and suggests amore complicated algorithm that avoids both underflow and overflow.Demmel [280, 1984] discusses in detail the effects of underflow on numericalsoftware and analyses different ways of implementing underflows, the main twoof which are flush to zero and gradual underflow (as used in IEEE standardarithmetic).
Cody [222, 1988] makes the following observations:The variety in the treatment of underflow is wondrous. There exist machines in which the exponent wraps around and underflowresults in a large floating-point number, perhaps with an accompanying sign change. There are also machines on which underflow isreplaced by the smallest in magnitude nonzero number, machineson which small positive numbers vanish when added to themselves,and machines on which such small numbers vanish when multiplied by 1.0. Finally, on IEEE-style machines, underflow may begraceful.25.9.
Multiple Precision ArithmeticIf we wish to work at a precision beyond that provided by a Fortran compiler’sDOUBLE PRECISION (or the quadruple precision supported as an extension bysome compilers) there are two basic approaches. The first is to implementhigher precision arithmetic using arrays of double precision variables, wherethe unevaluated sum of the elements of an array represents a high precisionvariable. Arithmetic can be implemented in this system using techniquessuch as those of Dekker [275, 1971] (cf. compensated summation, described in§4.3), as refined by Linnainmaa [706, 1981] and Priest [843, 1991], [844, 1992].These techniques require a “well-behaved” floating point arithmetic, and arenot easy to apply.The second approach is to use arrays to represent numbers in “standard”floating point form with a large base and a long mantissa spread across theelements of an array.
All the software described in this section is of the secondtype.The first major piece of multiple precision Fortran software was Brent’spackage [145, 1978], [146, 1978]. This is a portable collection of Fortran 66subroutines with wide capabilities, including the evaluation of special functions. Brent’s package represents multiple precision numbers as arrays ofintegers and operates on them with integer arithmetic.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
