3 数值健壮性
3 数值健壮性 |
|
从很多方面来说, 数值计算都是计算科学的核心, 因为很多计算科学牵涉到数学模型的模拟. 因此, 好的数值计算能力就成为关键所在. 传统可获得的数值能力含有基本的单精度和双精度实型数据, 复杂数据类型(只有单精度), 和丰富的库函数用于提供很广泛的数值计算. 偶尔的情况下会遇到附加的能力, 例如双精度复合型或者四倍精度实型数据. 这些传统的数值能力能够满足大量的科学计算应用的要求.
然而，在有些情况下，需要操作数的数值精度更优，或需要关于当前数值执行环境信息比从传统工具获得的可用信息更多，用以保证收敛，产生更精确的结果，或提供数值计算的其它形式的健壮性。第一点即数值精度的更优选择，是由Fortran 90的"type kind"结构和使用户和实现定义的函数通用化于不同类型参数的能力来提供，第二点即关于数值环境的信息，是由数值近拟模型和相应的返回与该模型相联系的环境信息的查寻本征函数提供，以下五小节中前三节描述了数值健壮性两个增强功能中的第一个功能，后两节描述了第二个功能。
Fortran 90 的 KIND 结构是对Fortran77 "* size"扩展的一个标准化、规则化、概括化的结构。例如，通常认为REAL*4与单精度REAL等价，REAL*8 与双精度DOUBLE PRECISION 等价，"*size"是不同实数类型的参数化形式。Fortran90将一个数据类型种类的概念形式化并把一个整型种类值与每个本征数据类型联系起来。这些kind值是依赖于具体实现的, 但是如果一个实现选择4作为单精度REAL，8作为DOUBLE PRECISION，那么定义REAL的两种方式是REAL(4)和REAL(KIND＝4)，定义DOUBLE PRECISION的两种方式是REAL(8)和REAL(KIND)=8。因此，在语句说明中定义REAL类型的完全形式化结构为：
REAL[<[KIND]=kind-valne>]
虽然有些实现可以选择参数值4作为单精度实数，参数值8作为双精度实数，以使得这些数与通常的“*4”和“*8”扩展类似或代表的种类型的对象的字节数(4和8的最初目的)，但是这样并不是对所有的体系结构都合适，例如，有些机器(象Cray向量超级计算机)，使用不同(于4)的字节数表示单精度实型对象，也许更重要的是，一个实现可以支持不只一种形式的单精度数：一种是缺省的实型REAL，另一种(一些)使用相同数量的存贮单元但是不同的表示方法。例如，在一台本地算法不是IEEE的机器上，一个Fortran实现可以支持基于IEEE算法另外形式的(单精度)实数类型。在Fortran90术语中，这仅仅是另外一种具有不同类型种类值的实数类型。因此，对每一个实现而言，没有明显的最优kind值，这就是kind值依赖实现原因。
幸运的是，有一种方法既可以将kind值实现的依赖性与应程序代码分离，同时又能使代码具有更强的可读性.这种技术就是使用KIND内部函数来为(适当命名的)整型常数建立正确的值。(KIND内部函数返回针对那个实现的参数类型的整型kind值)，例如，假设SINGLE和DOUBLE是整型常数的名字，分别具有(可能要求)单精度和双精度实数类型的合适的kind值，那么可用下列语句来定义单精度和双精度实型变量
REAL(SINGLE)…single precision variables…
REAL(DOUBLE)…double precision variables…
另一种"full blown"的定义方式为：
REAL(KIND=SINGLE)…single precision variables…
REAL(KIND=DOUBLE)…double precision variables…
在使用这种方式定义之前，SINGLE和DOUBLE必须具有合适的值，这通过以下语句来完成：
INTEGER， PARAMETER；；SINGLE=KIND(1.0),
DOUBLE=KIND(1D0)
(在本例中，采用了Fortran90的“面向实体”的定义方式，将INTEGER定义和PARAMETER定义合并到一个定义语句中)。在前面的例子中，使用的第一个KIND内部函数具有一个单精度实型参数(1.0)，这样它返回的值是单精度实型的整型种类值. 第二个例子中KIND内部函数具有一个双精度实型参数(1D0)，所以它返回双精度实型的整型种类值。使用这种技术的程序比较灵活，不考虑由实现选择的种类值。
复数COMPLEX类型使用与实数类型REAL一样的KIND值——复数种类实际上是实/虚两部分组成——, 对每个实型种类的支持, 都有一种一修复数种类实现的对应。这样, 定义复数对象的最佳方式为：
COMPLEX(SINGLE)...sigle precision complex variables...
COMPLEX(DOUBLE)...double precision complex variables...
进一步地，如果一个实现的本地算法不是IEEE，但它支持IEEE数据类型，则由IEEE算法的实现所指定的kind值可以赋给整型常量IEEE. 这时IEEE变量可以定义如下：
REAL(IEEE) ...real variables of type IEEE...
COMPLEX(IEEE) ...complex variables of type IEEE...
这使得对Fortran90的数值有用特征的讨论容易些，SELECTED_REAL_KIND内部函数允许程序员指定最小的十进制精度和/或指数范围特征. SELECTED_REAL_KIND函数有两个可选择的整型参数，至少要提供其中的一个参数，一个用来指定需要的十进制精度，另一个指定的(十进制)指数范围；然后SELECTED_REAL_KIND返回一个种类值，该种类值是由实现支持的符合或超过所指定条件的实数类型的最小种类值，如果没有这样的数据类型就会存在一个错误条件。因此，例如语句
INTEGER， PARAMETER::P9=SELECTED_REAL_KIND(9)
定义了一个整型常量P9，这个常量的实型种类值同有至少9位十进制精度的实型对象相适合。符合这种要求的变量可用如下语句定义：
REAL(P9) ... 9 drgit (at least) precision real variable...
在Sun工作站上，P9的值与DOUBLE(如前定义)一样；在克雷(Cray)巨型机上，P9的值与SINGLE一样。也就是说，REAL(P9)在Sun工作站上定义了双精度变量而在克雷(Cray)机上则定义了单精度变量。采用这种技术，可将整个程序写成基于变量所需精度／范围特性，而不是基于实现时的单精度和双精度实型数的缺省值，这本身对于数值健壮性来说就是一个有力的工具。
任何种类的实型常量可通过在由语言提供的缺省单精度实型常量后下划线和种类的值来形成，这样，用上面定义的kind值:
1.41_SINGLE 和 1.41 等价，
1.41_DOUBLE 和 1.41D0 等价，
1.41_IEEEE 是 1.41 的IEEE版本，
1.41_P9 是 1.41 的9位精度表示形式，典型地，它将与1.41_SINGLE或1.41_DOUBLE等价。
所有计算的内部函数继承实现提供的所有类型种类。例如，根据X的种类是SINGLE,DOUBLE,IEEE或P9， 函数COS(X)的返回结果分别是INGLE， DOUBLE， IEEE或P9的合适种类值. 这种继承特性在具有灵活健壮性的应用代码开发上起了重要作用。在Fortran 77 中内部函数同样具有继承性，但是它在健壮性方面的一个缺陷是用户和实现(和第三方软件卖主)提供的过程不能继承参数类型。Fortran90 弥补了这一缺陷，(在这一章“多样性”可假定意味着我们熟悉的Fortran 77 内部函数的继承特性可以被指定，因此扩展了用户定义的过程)
接口块可被用来为用户提供的过程集合指明一个类属名，或为已经存在的类属名增加过程。在下面的两个例子中，第一个接口块定义了一个新的与四个具体过程相联的类属名(SMOOTH)， 第二个接口块将COS内部函数扩展成类型RATIONAL的参数。
INTERFACE SMOOTH | SMOOTH 是类属名
INTERGER FUNCTION SMOOTH_INT(AA) | for procedures SMOOTH_INT
INTEGER :: AA (:,:) | SMOOTH_SINGLE
END FUNCTION SMOOTH_INT | SMOOTH_DOUBLE
| SMOOTH_RATIONAL
INTEGER FUNCTION SMOOTH_SINGLE(AA) REAL(SINGLE) :: AA(:,:) | AA is an assumed shape two- END FUNTION SMOOTH_SINGLE | dimensional array in each case.
INTEGER FUNCTION SMOOTH_DOUBLE(AA) REAL(DOUBLE) :: AA(:,:) END FUNCTION SMOOTH_DOUBLE
INTEGER FUNCTION SMOOTH_RATIONAL(AA) TYPE(RATIONAL) :: AA(:,:) END FUNCTION SMOOTH_RATIONAL
END INTERFACE
INTERFACE COS | Extends the generic properties FUNCTION RATIONAL_COS(X) | of COS to return results of TYPE(RATIONAL) :: RATIONAL_COM | type RATIONAL, assuming the TYPE(RATIONAL) :: X | argument is of type RATIONAL END FUNCTION RATIONAL_COS END INTERFACE
注意，在例子SMOOTH中，每种情况下参数(AA)具有不同的类型，但函数结果是相同类型(INTEGER). 函数结果也可能是某些情况下或在所有情况下都不相同，这种类属的唯一要求是在每个具体情况下哑元集合的 类型/种类/秩(type/kind/rank)和每个指定情况一致. 在例子COS中注意它与第一节第一页用户定义类型和操作符中提到的操作符"+"的扩展的概念上的相似性。
这些类属定义能力使得程序员可以使用SELECTED_REAL_KIND函数指定精度。因而，程序能配合应用的最优精度而不是集中在由实现提供的具体精度上。这不仅仅是开发精值软件的更自然的方法，而且也有利于数值健壮性. 在2.1节给出了例程SMOOTH的Fortran 77，C, 和 Fortran 90 版本。
对一个实现的数值特征的动态访问使得能够开发具有灵活数值健壮性的软件。Fortran 90 提供了这样的访问，它是通过将近似实数值的方法模型与相应的提取模型值的内部函数集合起来实现的，而这些内部函数(共有16个)被称为环境查询(9个)和数值操作(7个)内部函数
每种类型的实数被模型化表示为
x=sbe∑fib-i i=1,...,P (1)
其中:
这里x是实型值，
S是值的符号，
b是基数(通常是2)；b是对一个给定的实型kind值是个常数
p是基数b的精度；p对一个给定实型kind值是一个常数；
e是基数b的指数值；
fi是第i位数的值; 0<fi<b；
f1可以为0仅当所有fi都是0。
一个给定实型数的基本特性是它的基数b，基数精度p的值以及指数e的范围。
IEEE运算是基于二进制(b=2)表示，在这个运算中，p=24(单精度)，p=25(双精度)，-127<e<127；IEEE使用指数-127来表示 0 和 NaNs (非法值或溢出值)。非二进制的例子是 IBM370 的实数运算，在该运算中，b=16，p=6(单精度)，p=14(双精度)，127<e<127。在大多数实现中，不同实数类型kinds表示的主要差别在于p的值，尽管可能(或偶而发生)b的值或e的范围在不同类型间也有变化。
对一个给定的实型kind值，有9个重要的模型相关的固定特征值。这些值构成这种kind的“指纹”，它们在各种计算上下文中非常有用，在(a)数据对象定义的数值表达式中(b)任何计算中环境查寻内部函数允许程序员在任意时刻访问这些值。这些环境查询内部函数中每个都带有一个单一参数，这个参数可以是常量，标量变量，或数组变量，如果是变量，就不必定义它的值，因为这些函数使用的是只是参数的类型kind，而不是参数的值。
这9个特征值和与之相联的环境内部函数如表2所示：
| 实型特征值 | 内部函数名 |
| 十进制精度 | PRECISION |
| 十进制指数范围 | RANGE |
| 最大值 | HUGE |
| 最小值 | TINY |
| 1和b1-p中的较小值 | EPSILON |
| 基数 b | RADIX |
| p的值 | DIGITS |
| e的最小值 | MINEXPONENT |
| e的最大值 | MAXEXPONENT |
表 2
这9个值中，HUGE，TINY和EPSILON对编写健壮灵活的数值软件相当有用，虽然在特殊上下文中，其它值也同样必不可少。
七个数值操作函数允许程序员访问与给定那种类型具体kind值的相联系的重要模块值，给定值是第一个参数的值，它可以是具有那种kind的任何表达式，(这七个函数中有三个还带有第二个参数)，这七个有用的值和相应的数值操作本征函数如表3所示：
| 与参数值相关的值 | 内部函数名 |
| 指数值 e | EXPONENT |
| 小数部分 ，S∑fib-i | FRACTION |
| 最近值第(二个参数指定方向) | NEAREST |
| 接近参数的绝对空间的倒数 | RRSPACING |
| 通过第二个参数值改变的值 | SCALE |
| 将e设为第二个参数的值 | SET_EXPONENT |
| 接近参数的绝对空间 | SPACING |
表 3
一个(理想化的)例子将说明如何使用这些函数来编写健壮灵活的代码，在这个例子中使用牛顿迭代法用最小迭代次数寻找函数F的最精确的(对于所使用的实型kind)根。这个例子假设函数F和它的派生函数DF是可用的，X的值在牛顿迭代收敛至根的范围内已经建立。
...
do DX=F(X)/OF(X) | Compute the next delta-X X = X-DX if (DX<2*spacing(X)) exit | Stop if near the spacing end do | limits of that kind ... | in this region
| Copyright: NPACT |
|