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)令 ，并做如下计算：