图3.1 GPS Rayshooting模式的相关树结构片段
3.3 自动生成测试程序
基于IO相关分析的结果，DFT系统自动生成微分测试代码，分别对切线性模式的可靠性和运行代价做统计评价测试。特别地，DFT系统还可将任何模式参数都视为输入输出参数，生成在数学意义上等价的测试代码，这样处理的不利之处在于往往需要极高的存储开销。
3.4 基于语句级的代码优化
目前，DFT系统仅仅具备局地优化能力。在语句级微分实现上采用二元归约的方法对微分代码进行优化是DFT系统的一个重要特色。根据右端表达式的乘法复杂性及含变元数目的不同，DFT系统采取不同的分解策略。二元归约的方法避免了微分计算中的许多冗余计算，在一些复杂的非线性表达式的微分计算中具有最小的计算代价，同时也非常适合于微分系统的软件实现。同时，对于某些特殊的运算操作（除法、乘方）和特殊函数（如sqrt、exp），DFT系统较好地利用了基态值计算得到的中间结果，避免了微分实现中的冗余计算。
4.系统应用
运用自动微分工具得到的切线性模式，可以在无截断误差意义下求解函数的数值微分和导数、稀疏雅可比矩阵。同时这些结果在数值参数敏感性分析、非线性最优化以及其它数值理论分析中有着非常重要的应用。这里简单介绍切线性模式的几个基本应用。
4.1 符号导数和微分
如果输入为数学关系式，DFT系统可以自动生成对应的微分表达式和梯度，而与数学关系式的复杂程度无关。例如我们输入关系式：
， （1）
DFT系统将自动生成其符号微分形式及其梯度形式分别为
， （2）
4.2 数值导数和微分
切线性模式最基本的应用就是在一定扰动输入下求解输出变量的扰动（响应）。表4.1给出了DFT系统在对IAP 9L模式、GPS Rayshooting模式和GPS Raytrace模式三个数值模式做切线性化的具体应用中，一些不同计算粒度、不同引用深度和不同程序风格的核心子过程，以及它们的切线性模式在SGI 2000上运行的统计评价测试结果，其中切线性模式的可靠性指标都准确到六个有效数字以上，在运行时间、存储开销和代码复杂性方面分别是原模式的两倍左右，比较接近于理想的微分代价结果（1.5倍）。除了IAP 9L模式由于过于复杂仅做粗略统计外，其余模式都用非注释语句行数来表示各自的代码复杂性。
表4.1 DFT系统在三个数值模式中的统计评价测试结果
性能指标
对象模式 运行时间（10-3秒） 存储开销（字节数） 代码复杂性
原模式 切线性
模式
原模式 切线性
模式
原模式 切线性
模式
Xyz2g 2.530 6.160 5524 11048 55 89
IntCIRA 1.560 2.750 1334 2661 41 65
Dabel 0.035 0.072 60 120 27 49
LSS 8.300 17.50 669 1338 79 143
RP 42.40 85.10 3605 7210 22 38
Vgrad1 0.100 0.212 18564 36828 24 54
RefGr 43.00 86.00 718654 1437308 35 78
LL2JK 0.626 1.350 2622 5244 22 32
RayFind4 62.70
×103 125.4
×103 9856 18212 111 179
EPSIMP 1.760 11.50 4455 8910 13 27
Hlimits 0.830 1.880 242577 484254 37 74
Int3sL 26.90 51.20 820029 1639458 46 90
MAKE
NCEP 1340 3920 722925 1445850 45 84
Curvcent 0.013 0.0385 27 54 27 54
DYFRM 3.800
×103 7.250
×103 5000* 9500* 161 279
PHYSIC 2.750
×103 5.385
×103 3000 5000* 1399*
(含注释行) 2826*
(含注释行)
适当设置输入扰动的初值，运用切线性模式可以简单求解输出变量对输入的偏导数。例如，对于一个含有 个输入参数的实型函数
（3）
这里设 ， 。运用DFT系统，可以得到对应的切线性模式
（4）
其中 ， 为切线性模式的扰动输入参数。可以通过以下办法来求得偏导数：
（5）
其中 。如果对于某个 既是输入参数又是输出参数，可以类似以下过程引用的办法来处理。对于过程引用的情形，例如一个含有 个输入参数的子过程
（6）
其中 ， 为输入参数； ， 为输出参数； ， 既为输入参数又为输出参数。运用DFT系统，可以得到对应的切线性模式为
（7）
其中 ， ， ， 分别为切线性模式的微分扰动输入、输出和输入输出参数。可以通过以下输入扰动设置并引用切线性模式（7）来求得偏导数：
a) 设置 ； （ ， ）； （ ）可以同时求得 （ ）和 （ ），其中 。
b) 设置 （ ）； ； （ ， ）可以同时求得 （ ）和 （ ），其中 。
4.3 稀疏雅可比矩阵
运用上节讨论的方法来求解稀疏雅可比矩阵,具有极高的计算代价。例如，一个含 个独立和 个依赖参数的子过程，为求解整个雅可比矩阵就需要反复调用 次切线性模式，当 相当大时，这对许多实际的数值计算问题是不能接受的。事实上，如果雅可比矩阵的任意两列（行）相互正交，那么可以通过适当设置扰动输入值，这两列（行）的元素就可以通过一次引用切线性模式（伴随模式）完全得到。设 和 分别为雅可比矩阵的行宽度和列宽度，即各行和各列非零元素数目的最大值，显然有 ， 。这里介绍几种常用的求解方法。
正向积分 当 时，通常采用切线性模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择右乘初始输入矩阵，可以获得接近 的计算时间代价。DFT系统采用一种逐列（行）求解的方法，来有效求解右（左）乘初始输入矩阵。其基本思路是：按照某种列次序考察雅可比矩阵的各列；考察当前列中所有非零元素，并对这些非零元素所在行的行向量做类似模二和累加运算（即将非零元素视为逻辑“1”， 零元素视为逻辑“0”），从而得到一个描述当前列与各行存在“某种”相关的标志向量（其元素都是“1”或“0”）；依据此标志向量，就很容易得到一个与之正交的列初始向量 ，其中与当前列序号对应的元素设置为“1”，而与标志向量中非零元素序号对应的元素设置为“0”， 与标志向量中非零元素序号对应的元素设置为“-1”，显然，该列初始向量是唯一的，并且对应着当前右乘初始输入矩阵的最后一列；逐一考察已求解得到的列初始向量，如果某列初始向量与当前求解得到的列初始向量按下面定义的乘法（见过程4）正交，那么这两列就可以合并，即将当前列初始向量中非“-1”的元素按照对应关系分别赋值给该初始向量，并从记录中删除当前列初始向量；重复以上过程，继续按照给定列次序考察雅可比矩阵的“下一列”。不难说明，按照不同列次序求解得到的右乘初始输入矩阵可能不同。其中逐列求解右乘初始输入矩阵的过程可以简单叙述为：
1）将右乘初始输入矩阵 所有元素的初值均设置为 ， ， 。 。
2）如果 ，转6）。否则，如果雅可比矩阵 的第 列中的所有元素均为 ， ，重复2）的判断。否则转3）。
3）计算标志向量 。令 ，做如下计算：
， ;
4）设 为 的列向量。在 上定义乘法 ，对任意的 ，我们有：a） ；b)如果 ，必有 和 。然后，做如下计算：
, ;
， 6);
2);
5）令 ，并做如下计算：
， ;
令 , 。如果 ，转6）；否则，重复2）的判断。
6）对 ， ，如果 ，则 。取 的前 列，这样，我们就得到了一个 维右乘初始输入矩阵。
这里需要说明的是，运用上面的方法求得的右乘初始输入矩阵不仅与求解雅可比矩阵的列序有关，而且与过程4）中的合并顺序也有关系。至于如何最优求解右乘初始输入矩阵，目前还很难讨论清楚。但是，大量模拟试验结果表明，运用上面自然次序求得的右乘初始输入矩阵宽度 已经非常接近于其下界值 。
反向积分 当 和 时，通常采用伴随模式来计算雅可比矩阵。根据雅可比矩阵的稀疏结构，适当选择左乘初始输入矩阵，可以获得接近 的计算时间代价。其中左乘初始输入矩阵的求解过程完全可以按照上面的方法进行，但是在处理前必须先将雅可比矩阵转置，最后还需将得到的初始输入矩阵转
