0 引言

线性变换方程包含最基本的线性方程组、线性矩阵方程, 如Lyapunov矩阵方程, Sylvester矩阵方程等.这些方程(组) 的求解问题广泛来源于控制理论^[1]、图像恢复^[2]、迁移理论中Riccati微分方程的数值解^[3]、模型降阶^[4-6]、线性系统的稳定性分析^[7]等许多科学领域.

目前, 求解大型线性变换方程通常采用迭代法, 如ADI方法^[8](Alternating Direction Implicit) 及Krylov子空间法^[9-13]等.ADI方法的不便之处在于确定最优参数; 共轭梯度法^[14](conjugate gradient method, 简称CG法) 具有存储量小、适合并行计算等优点, 但它仅适用于系数矩阵为对称正定的情形; 变形共轭梯度法^[15](modified conjugate gradient method, 简称MCG法) 的收敛速度与系数矩阵的条件数密切相关, 条件数越小, 收敛速度越快, 因此MCG算法不适用于条件数较大的情形; SYMMLQ算法^[16]适合求解系数矩阵为对称这一情形, 其易于实现并行计算且收敛速度快.但是该算法的缺点是在并行实现过程中, 需要两次全归约, 由此引起的大规模全局通讯使得算法的可扩放性比较差^[17-18].

基于对上述算法的研究, 本文首先对SYMMLQ算法实现并行运算, 然后将该并行算法进行改进, 以减少并行计算时间, 提高并行效率.文中涉及到的所有矩阵都是实矩阵, 定义同阶矩阵A与B的内积为[A, B]=tr (A^TB), 由此导出矩阵的Frobenius范数${{\left\| \mathit{\boldsymbol{A}} \right\|}_{F}}=\sqrt{\left[\mathit{\boldsymbol{A}}, \mathit{\boldsymbol{A}} \right]} $.

1 SYMMLQ算法的推广 1.1 求解线性对称变换方程的SYMMLQ算法

文献[16]给出了求解对称情形下线性方程组的SYMMLQ算法, 文中将该算法推广到求解线性对称变换方程, 以拓宽其使用范围.

设线性变换T是n维欧式空间V上的对称变换, 即 X, Y∈V, 有[T(X), Y]=[X, T(Y)].线性对称变换方程为

$ T\left( \mathit{\boldsymbol{X}} \right) = \mathit{\boldsymbol{F}}, $

(1)

其中X∈V是未知的元素, F∈V是已知元素.

根据文献[12]的Lanczos过程, 推广到线性对称变换方程，可得到如下关系式

$ T\left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}},{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}, \cdots ,{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right) = \left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}},{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}, \cdots ,{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right){\mathit{\boldsymbol{T}}_k} + {b_k}{\mathit{\boldsymbol{Q}}^{\left( {k + 1} \right)}}\mathit{\boldsymbol{e}}_k^{\rm{T}}, $

其中Q⁽¹⁾, Q⁽²⁾, …, Q^(k)表示由Lanczos正交过程得到的n维欧式空间V中相互正交的单位元素(i, j=1, …, k).e_k^T表示k维空间第k个向量, 即e_k^T=(0, 0, …, 0, 1).T_k是一个k×k的对称三对角矩阵, 即

$ {\mathit{\boldsymbol{T}}_k} = \left( {\begin{array}{*{20}{c}} {{a_1}}&{{b_1}}&{}&{}&{}\\ {{b_1}}&{{a_2}}&{{b_2}}&{}&{}\\ {}&{{b_2}}& \ddots & \ddots &{}\\ {}&{}& \ddots & \ddots &{{b_{k - 1}}}\\ {}&{}&{}&{{b_{k - 1}}}&{{a_k}} \end{array}} \right). $

考虑方程(1) 的求解, 令

$ {\mathit{\boldsymbol{X}}^{\left( k \right)}} = {\mathit{\boldsymbol{X}}^{\left( 0 \right)}} + \left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}},{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}, \cdots ,{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right){\mathit{\boldsymbol{f}}^{\left( k \right)}}, $

其中f^(k)是k×1列向量.取f^(k)使得R^(k)=F-T(X^(k)) 正交于Q⁽¹⁾, Q⁽²⁾, …, Q^(k), 即

$ \left[ {{\mathit{\boldsymbol{Q}}^{\left( i \right)}},{\mathit{\boldsymbol{R}}^{\left( k \right)}}} \right] = 0\left( {i = 1, \cdots ,k} \right). $

(2)

因为

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}^{\left( k \right)}} = {\mathit{\boldsymbol{R}}^{\left( 0 \right)}} - T\left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}},{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}, \cdots ,{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right){\mathit{\boldsymbol{f}}^{\left( k \right)}}}\\ { = {\mathit{\boldsymbol{R}}^{\left( 0 \right)}} - \left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}},{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}, \cdots ,{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right){\mathit{\boldsymbol{T}}_k}{\mathit{\boldsymbol{f}}^{\left( k \right)}} - {b_k}{\mathit{\boldsymbol{Q}}^{\left( {k + 1} \right)}}\mathit{\boldsymbol{f}}_k^{\left( k \right)},} \end{array} $

(3)

将式(3) 左右两端分别与Q⁽ⁱ⁾(i=1, …, k) 做内积, 由式(2) 和(3), 可得

$ {\left\| {{\mathit{\boldsymbol{R}}^{\left( 0 \right)}}} \right\|_\mathit{\boldsymbol{F}}}{\mathit{\boldsymbol{e}}_1} - {\mathit{\boldsymbol{T}}_k}{\mathit{\boldsymbol{f}}^{\left( k \right)}} = 0. $

这样f^(k)是方程${{\mathit{\boldsymbol{T}}}_{k}}f_{k}^{\left( k \right)}={{\left\| {{\mathit{\boldsymbol{R}}}^{\left( 0 \right)}} \right\|}_{F}}{{\mathit{\boldsymbol{e}}}_{1}}$的解.

因为对称三对角矩阵T_k可能奇异, 因而对它采用LQ分解, 该过程类似于文献[12].于是, 得到求解线性对称变换方程(1) 的SYMMLQ算法(算法1):

(1) 任意给定初始元素X⁽⁰⁾∈V, 置k:=1, 计算

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}^{\left( 0 \right)}} = \mathit{\boldsymbol{F}} - T\left( {{\mathit{\boldsymbol{X}}^{\left( 0 \right)}}} \right),{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}} = \frac{{{\mathit{\boldsymbol{R}}^{\left( 0 \right)}}}}{{\left\| {{\mathit{\boldsymbol{R}}^{\left( 0 \right)}}} \right\|}},{a_1} = \left[ {T\left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}}} \right),{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}}} \right],{{\tilde r}_1} = {a_1},}\\ {{b_1} = \left\| {T\left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}}} \right) - {a_1}{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}}} \right\|,{{\tilde \xi }_1} = \frac{{\left\| {{\mathit{\boldsymbol{R}}^{\left( 0 \right)}}} \right\|}}{{{{\tilde r}_1}}},{{\mathit{\boldsymbol{\tilde W}}}^{\left( 1 \right)}} = {\mathit{\boldsymbol{Q}}^{\left( 1 \right)}},{\mathit{\boldsymbol{X}}^{\left( 1 \right)}} = {\mathit{\boldsymbol{X}}^{\left( 0 \right)}} + {{\tilde \xi }_1}{{\mathit{\boldsymbol{\tilde W}}}^{\left( 1 \right)}};} \end{array} $

(2) 置k:=2, 计算

$ \begin{array}{l} {\mathit{\boldsymbol{Q}}^{\left( 2 \right)}} = \frac{{\left( {T\left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}}} \right) - {a_1}{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}}} \right)}}{{{b_1}}},{a_2} = \left[ {T\left( {{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}} \right),{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}} \right],\\ {b_2} = \left\| {T\left( {{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}} \right) - {a_2}{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}} - {b_1}{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}}} \right\|,{c_1} = \frac{{{a_1}}}{{{{\left( {a_1^2 + b_1^2} \right)}^{1/2}}}},{s_1} = \frac{{{b_1}}}{{{{\left( {a_1^2 + b_1^2} \right)}^{1/2}}}},\\ {\delta _2} = {b_1}{c_1} + {a_2}{s_1},{\varepsilon _3} = {b_2}{s_1},{r_1} = {\left( {a_1^2 + b_1^2} \right)^{1/2}},{\xi _1} = \left\| {{\mathit{\boldsymbol{r}}_0}} \right\|/{r_1},\\ {{\tilde r}_2} = {b_1}{s_1} - {a_2}{c_1}{r_2} = {\left( {\tilde r_2^2 + b_2^2} \right)^{1/2}},{{\tilde \delta }_3} = - {b_2}{c_1},{{\tilde \xi }_2} = - {\delta _2}{\xi _1}/{{\tilde r}_2},{\xi _2} = - \frac{{{\delta _2}{\xi _1}}}{{{r_2}}},\\ {\mathit{\boldsymbol{W}}^{\left( 1 \right)}} = {c_1}{{\mathit{\boldsymbol{\tilde W}}}^{\left( 1 \right)}} + {s_1}{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}},{{\mathit{\boldsymbol{\tilde W}}}^{\left( 2 \right)}} = {s_1}{{\mathit{\boldsymbol{\tilde W}}}^{\left( 1 \right)}} - {c_1}{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}},\\ {\mathit{\boldsymbol{Y}}^{\left( 1 \right)}} = {\mathit{\boldsymbol{X}}^{\left( 0 \right)}} + {\xi _1}{\mathit{\boldsymbol{W}}^{\left( 1 \right)}},{\mathit{\boldsymbol{X}}^{\left( 2 \right)}} = {\mathit{\boldsymbol{Y}}^{\left( 1 \right)}} + {{\tilde \xi }_2}{\mathit{\boldsymbol{W}}^{\left( 2 \right)}}; \end{array} $

(3) 计算ΔX^(k-1)=X^(k-1)-X^(k-2), 若‖ΔX^(k-1)‖ < ε停止计算; 否则, 计算

$ \begin{array}{l} {\mathit{\boldsymbol{Q}}^{\left( k \right)}} = \frac{{\left( {T\left( {{\mathit{\boldsymbol{Q}}^{\left( {k - 1} \right)}}} \right) - {a_{k - 1}}{\mathit{\boldsymbol{Q}}^{\left( {k - 1} \right)}} - {b_{k - 2}}{\mathit{\boldsymbol{Q}}^{\left( {k - 2} \right)}}} \right)}}{{{b_{k - 1}}}},{a_k} = \left[ {T\left( {{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right),{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right],\\ {b_k} = \left\| {T\left( {{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right) - {a_k}{\mathit{\boldsymbol{Q}}^{\left( k \right)}} - {b_{k - 1}}{\mathit{\boldsymbol{Q}}^{\left( {k - 1} \right)}}} \right\|,{c_{k - 1}} = \frac{{{{\tilde r}_{k - 1}}}}{{{{\left( {\tilde r_{k - 1}^2 + b_{k - 1}^2} \right)}^{1/2}}}},{s_{k - 1}} = \frac{{{b_{k - 1}}}}{{{{\left( {\tilde r_{k - 1}^2 + b_{k - 1}^2} \right)}^{1/2}}}},\\ {\delta _k} = {{\tilde \delta }_k}{c_{k - 1}} + {a_k}{s_{k - 1}},{\varepsilon _{k + 1}} = {b_k}{s_{k - 1}},{{\tilde r}_k} = {{\tilde \delta }_k}{s_{k - 1}} - {a_k}{c_{k - 1}},{r_k} = {\left( {\tilde r_k^2 + b_k^2} \right)^{1/2}},\\ {{\tilde \delta }_{k + 1}} = - {b_k}{c_{k - 1}},{{\tilde \xi }_k} = - \left( {{\varepsilon _k}{\xi _{k - 2}} + {\delta _k}{\xi _{k - 1}}} \right)/{{\tilde r}_k},{\xi _k} = - \left( {{\varepsilon _k}{\xi _{k - 2}} + {\delta _k}{\xi _{k - 1}}} \right)/{r_k},\\ {\mathit{\boldsymbol{W}}^{\left( {k - 1} \right)}} = {c_{k - 1}}{{\mathit{\boldsymbol{\tilde W}}}^{\left( {k - 1} \right)}} + {s_{k - 1}}{\mathit{\boldsymbol{Q}}^{\left( k \right)}},{{\mathit{\boldsymbol{\tilde W}}}^{\left( k \right)}} = {s_{k - 1}}{{\mathit{\boldsymbol{\tilde W}}}^{\left( {k - 1} \right)}} - {c_{k - 1}}{\mathit{\boldsymbol{Q}}^{\left( k \right)}},\\ {\mathit{\boldsymbol{Y}}^{\left( {k - 1} \right)}} = {\mathit{\boldsymbol{Y}}^{\left( {k - 2} \right)}} + {\xi _{k - 1}}{\mathit{\boldsymbol{W}}^{\left( {k - 1} \right)}},{\mathit{\boldsymbol{X}}^{\left( k \right)}} = {\mathit{\boldsymbol{Y}}^{\left( {k - 1} \right)}} + {{\tilde \xi }_k}{{\mathit{\boldsymbol{\tilde W}}}^{\left( k \right)}}; \end{array} $

(4) 置k:=k+1, 转(3).

1.2 SYMMLQ算法的收敛性分析

定理1  按照1.1中推导过程所求得的近似解X^(k)有‖R^(k)‖_F=|b_kf_k^(k)|, 其中f_k^(k)是f^(k)的第k个分量.

证明   由于

$ \begin{array}{l} {\mathit{\boldsymbol{R}}^{\left( k \right)}} = {\mathit{\boldsymbol{R}}^{\left( 0 \right)}} - T\left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}},{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}, \cdots ,{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right)\mathit{\boldsymbol{f}}\left( k \right) = \\ \;\;\;\;\;\;\;\;\;\left( {{\mathit{\boldsymbol{Q}}^{\left( 1 \right)}},{\mathit{\boldsymbol{Q}}^{\left( 2 \right)}}, \cdots ,{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right)\left( {{{\left\| {{\mathit{\boldsymbol{R}}^{\left( 0 \right)}}} \right\|}_F}{\mathit{\boldsymbol{e}}_1} - {\mathit{\boldsymbol{T}}_k}{\mathit{\boldsymbol{f}}^{\left( k \right)}}} \right) - {b_k}{\mathit{\boldsymbol{Q}}^{\left( {k + 1} \right)}}\mathit{\boldsymbol{f}}_k^{\left( k \right)}, \end{array} $

因为Q^(k+1)正交于Q⁽¹⁾, Q⁽²⁾, …, Q^(k), 所以

$ \left\| {{\mathit{\boldsymbol{R}}^{\left( k \right)}}} \right\|_F^2 = {\left\| {{{\left\| {{\mathit{\boldsymbol{R}}^{\left( 0 \right)}}} \right\|}_F}{\mathit{\boldsymbol{e}}_1} - {\mathit{\boldsymbol{T}}_k}{\mathit{\boldsymbol{f}}^{\left( k \right)}}} \right\|^2} + b_k^2\mathit{\boldsymbol{f}}_k^{\left( k \right)2}. $

而第一项为0, 故‖R^(k)‖_F=|b_kf_k^(k)|.证毕.

因此当|b_kf_k^(k)|很小时, 近似解较好地接近真解.

定理2  残量‖R^(k)‖满足[R^(k), R^(j)]=0, 其中j=1, 2, …, k-1.

证明   由上述定理1的证明过程可知

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{R}}^{\left( k \right)}} = - {b_k}{\mathit{\boldsymbol{Q}}^{\left( {k + 1} \right)}}\mathit{\boldsymbol{f}}_k^{\left( k \right)},{\mathit{\boldsymbol{R}}^{\left( j \right)}} = - {b_{jk}}{\mathit{\boldsymbol{Q}}^{\left( {j + 1} \right)}}\mathit{\boldsymbol{f}}_j^{\left( j \right)},}\\ {\left[ {{\mathit{\boldsymbol{Q}}^{\left( {k + 1} \right)}},{\mathit{\boldsymbol{Q}}^{\left( {j + 1} \right)}}} \right] = 0\left( {j = 1,2, \cdots ,k - 1} \right),} \end{array} $

所以[R^(k), R^(j)]=0.证毕.

定理3  设方程(1) 相容, 则对任意初始元素X⁽⁰⁾∈V, 算法1至多需n步有R⁽ⁿ⁾=0, 即至多需n步算法收敛(n为欧式空间的维数).

证明  n维欧式空间V中两两正交的非零元素至多为n个, 又由定理2的结论可知对任意初始元素X⁽⁰⁾∈V所得到的残量R⁽⁰⁾与R⁽¹⁾, …, R^(n-1), R⁽ⁿ⁾相互正交, 因此算法1至多需n步有R⁽ⁿ⁾=0, 即结论成立.证毕.

2 SYMMLQ算法的并行实现

针对线性空间V=R^n×n上的对称变换T(X)=AX+XB, 即方程(1) 形式如下:

$ \mathit{\boldsymbol{AX}} + \mathit{\boldsymbol{XB}} = \mathit{\boldsymbol{F}}. $

(4)

其中A, B, F∈R^n×n, 并且A=A^T, B=B^T.下面给出求解方程(4) 的SYMMLQ算法的具体并行实现过程.设处理机台数为p, p整除n, 即n=pl, p_i(i=1, 2, …, p) 表示第i台处理机.记

$ \begin{align} & \mathit{\boldsymbol{A}}={{\left( \mathit{\boldsymbol{A}}_{1}^{T},\mathit{\boldsymbol{A}}_{2}^{T},\cdots ,\mathit{\boldsymbol{A}}_{p}^{T} \right)}^{\rm{T}}},\mathit{\boldsymbol{B}}={{\left( \mathit{\boldsymbol{B}}_{1}^{T},\mathit{\boldsymbol{B}}_{2}^{T},\cdots ,\mathit{\boldsymbol{B}}_{p}^{T} \right)}^{\rm{T}}},\mathit{\boldsymbol{F}}={{\left( \mathit{\boldsymbol{F}}_{1}^{T},\mathit{\boldsymbol{F}}_{2}^{T},\cdots ,\mathit{\boldsymbol{F}}_{p}^{T} \right)}^{\rm{T}}}, \\ & {{\mathit{\boldsymbol{X}}}^{\left( k \right)}}={{\left( {{\left( \mathit{\boldsymbol{X}}_{1}^{\left( k \right)} \right)}^{\rm{T}}},{{\left( \mathit{\boldsymbol{X}}_{2}^{\left( k \right)} \right)}^{\rm{T}}},\cdots ,{{\left( \mathit{\boldsymbol{X}}_{p}^{\left( k \right)} \right)}^{\rm{T}}} \right)}^{\rm{T}}},{{\mathit{\boldsymbol{R}}}^{\left( k \right)}}= \\ & {{\left( {{\left( \mathit{\boldsymbol{R}}_{1}^{\left( k \right)} \right)}^{\rm{T}}},{{\left( \mathit{\boldsymbol{R}}_{2}^{\left( k \right)} \right)}^{\rm{T}}},\cdots ,{{\left( \mathit{\boldsymbol{R}}_{p}^{\left( k \right)} \right)}^{\rm{T}}} \right)}^{\rm{T}}}. \\ \end{align} $

类似地, 矩阵变量Q^(k), W^(k), ${{\mathit{\boldsymbol{\tilde{W}}}}^{\left( k \right)}} $, Y^(k), G^(k), H^(k), L^(k)均按照上述方式进行行行划分^[17-18].其中, A_i, B_i, F_i, X_i^(k), R_i^(k)是l×n的行块子矩阵(i=1, 2, …, p).

注   算法涉及到的所有矩阵乘法皆采用行行划分矩阵乘法的并行算法^[17-18], 文中所涉及的矩阵乘法并行计算不再描述.

2.1 未改进的SYMMLQ算法的并行实现

(1) 存储方式:将A_i, B_i, F_i, X_i⁽⁰⁾均存储在p_i(i=1, 2, …, p) 处理机中.

(2) 计算过程:

①在处理机p_i中, 计算R_i⁽⁰⁾=A_iX⁽⁰⁾+X_i⁽⁰⁾B; 再计算内积[R_i⁽⁰⁾, R_i⁽⁰⁾], 全归约后得到[R⁽⁰⁾, R⁽⁰⁾], 然后计算Q_i⁽¹⁾=R_i⁽⁰⁾/‖R⁽⁰⁾‖;

②在处理机p_i中, 计算G_i⁽¹⁾=A_iQ⁽¹⁾+Q_i⁽¹⁾B, 再计算内积[G_i⁽¹⁾, Q_i⁽¹⁾], 全归约后得到a₁=[G⁽¹⁾, Q⁽¹⁾]，计算H_i⁽¹⁾=G_i⁽¹⁾-a₁Q_i⁽¹⁾, 计算内积[H_i⁽¹⁾, H_i⁽¹⁾], 全归约后得到[H⁽¹⁾, H⁽¹⁾], 计算${{b}_{1}}=\sqrt{\left[{{\mathit{\boldsymbol{H}}}^{\left( 1 \right)}}, {{\mathit{\boldsymbol{H}}}^{\left( 1 \right)}} \right]} $, 在各进程计算变量${{{\tilde{r}}}_{1}}, {{{\tilde{\xi }}}_{1}}, \mathit{\boldsymbol{\tilde{W}}}_{i}^{\left( 1 \right)} $以及$\mathit{\boldsymbol{X}}_{i}^{\left( 1 \right)}=\mathit{\boldsymbol{X}}_{i}^{\left( 0 \right)}+{{{\tilde{\xi }}}_{1}}\mathit{\boldsymbol{\tilde{W}}}_{i}^{\left( 1 \right)} $;

③在p_i中, 计算$\mathit{Q}_{i}^{\left( 2 \right)}=\frac{\mathit{\boldsymbol{H}}_{i}^{\left( 1 \right)}}{{{b}_{1}}} $, 然后计算G_i⁽²⁾=A_iQ⁽²⁾+Q_i⁽²⁾B, 计算内积[G_i⁽²⁾, Q_i⁽²⁾], 全归约后得到a₂=[G⁽²⁾, Q⁽²⁾]，计算H_i⁽²⁾=G_i⁽²⁾-a₂Q_i⁽²⁾-b₁Q_i⁽¹⁾, 计算内积[H_i⁽²⁾, H_i⁽²⁾], 全归约后得到[H⁽²⁾, H⁽²⁾], 计算$ {{b}_{2}}=\sqrt{\left[{{\mathit{\boldsymbol{H}}}^{\left( 2 \right)}}, {{\mathit{\boldsymbol{H}}}^{\left( 2 \right)}} \right]}$;

④在各进程中计算变量c₁, s₁, δ₂, ε₃, r₁, ξ₁, ${{{\tilde{r}}}_{2}} $, r₂, ${{{\tilde{\delta }}}_{3}} $, 若$ {{{\tilde{r}}}_{2}} < 1.0\times {{10}^{-15}}$, 停止计算; 否则计算${{{\tilde{\xi }}}_{2}}=-{{\delta }_{2}}{{\xi }_{1}}/{{{\tilde{r}}}_{2}} $;

⑤在I=(I₁^T, I₂^T, …, I_p^T)^T中, 计算变量ξ₂以及$\mathit{\boldsymbol{W}}_{i}^{\left( 1 \right)}, \mathit{\boldsymbol{\tilde{W}}}_{i}^{\left( 2 \right)}, \mathit{\boldsymbol{Y}}_{i}^{\left( 1 \right)}, \mathit{\boldsymbol{X}}_{i}^{\left( 2 \right)} $.

(3) 循环过程:

步骤1   计算$\left\| {{\mathit{\boldsymbol{R}}}^{\left( k-1 \right)}} \right\|=\left| {{b}_{k-1}}\left( {{s}_{k-2}}{{\xi }_{k-2}}-{{c}_{k-2}}{{{\tilde{\xi }}}_{k-1}} \right) \right| $, 若‖R^(k-1)‖ < ε停止计算, 否则转步骤2;

步骤2  在p_i中, 计算$\mathit{\boldsymbol{Q}}_{i}^{\left( k \right)}=\frac{H_{i}^{\left( k-1 \right)}}{{{b}_{k-1}}}, \mathit{\boldsymbol{G}}_{i}^{\left( k \right)}={{\mathit{\boldsymbol{A}}}_{i}}{{\mathit{\boldsymbol{Q}}}^{\left( k \right)}}+\mathit{\boldsymbol{Q}}_{i}^{\left( k \right)}\mathit{\boldsymbol{B}} $, 然后计算内积$\left[\mathit{\boldsymbol{G}}_{i}^{\left( k \right)}, \mathit{\boldsymbol{Q}}_{i}^{\left( k \right)} \right] $, 全归约后得到a_k=[G^(k),Q^(k)], 计算$ \mathit{\boldsymbol{H}}_{i}^{\left( k \right)}=\mathit{\boldsymbol{G}}_{i}^{\left( k \right)}-{{a}_{k}}\mathit{\boldsymbol{Q}}_{i}^{\left( k \right)}-{{b}_{k-1}}\mathit{\boldsymbol{Q}}_{i}^{\left( k-1 \right)}$,计算内积$\left[\mathit{\boldsymbol{H}}_{i}^{\left( k \right)}, \mathit{\boldsymbol{H}}_{i}^{\left( k \right)} \right] $, 全归约后得到[H^(k),H^(k)],计算 $ {{b}_{k}}=\sqrt{\left[{{\mathit{\boldsymbol{H}}}^{\left( k \right)}}, {{\mathit{\boldsymbol{H}}}^{\left( k \right)}} \right]}$;

步骤3   在各进程计算变量$ {{c}_{k-1}}, {{s}_{k-1}}, {{\delta }_{k}}, {{\varepsilon }_{k+1}}, {{r}_{k}}, {{{\tilde{r}}}_{k}}, {{{\tilde{\delta }}}_{k+1}}$, 若$ {{{\tilde{r}}}_{k}} < 1.0\times {{10}^{-15}}$则停止计算, 否则计算$ {{{\tilde{\xi }}}_{k}}=-\left( {{\varepsilon }_{k}}{{\xi }_{k-2}}+{{\delta }_{k}}{{\xi }_{k-1}} \right)/{{{\tilde{r}}}_{k}}$;

步骤4   在p_i中, 计算变量${{\xi }_{k}}=-\left( {{\varepsilon }_{k}}{{\xi }_{k-2}}+{{\delta }_{k}}{{\xi }_{k-1}} \right)/{{r}_{k}} $以及

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{W}}_i^{\left( {k - 1} \right)} = {c_{k - 1}}\mathit{\boldsymbol{\tilde W}}_i^{\left( {k - 1} \right)} + {s_{k - 1}}\mathit{\boldsymbol{Q}}_i^{\left( {k} \right)},\mathit{\boldsymbol{\tilde W}}_i^{\left( k \right)} = {s_{k - 1}}\mathit{\boldsymbol{\tilde W}}_i^{\left( {k - 1} \right)} - {c_{k - 1}}\mathit{\boldsymbol{Q}}_i^{\left( {k} \right)},}\\ {\mathit{\boldsymbol{Y}}_i^{\left( {k - 1} \right)} = \mathit{\boldsymbol{Y}}_i^{\left( {k - 2} \right)} + {\xi _{k - 1}}\mathit{\boldsymbol{W}}_i^{\left( {k - 1} \right)},\mathit{\boldsymbol{X}}_i^{\left( k \right)} = \mathit{\boldsymbol{Y}}_i^{\left( {k - 1} \right)} + {{\tilde \xi }_k}\mathit{\boldsymbol{\tilde W}}_i^{\left( k \right)};} \end{array} $

步骤5   置k:=k+1, 返回步骤1.

2.2 改进的SYMMLQ算法的并行实现

在SYMMLQ算法的并行实现过程中, 每迭代一次需要计算

$ {a_k} = \left[ {T\left( {{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right),{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right].{b_k} = \left\| {T\left( {{\mathit{\boldsymbol{Q}}^{\left( k \right)}}} \right) - {a_k}{\mathit{\boldsymbol{Q}}^{\left( k \right)}} - {b_{k - 1}}{\mathit{\boldsymbol{Q}}^{\left( {k - 1} \right)}}} \right\|, $

显然需要2次全归约.为了保证改进的SYMMLQ算法的结构稳定性, 且基本不影响原算法的计算精度, 在经过多次尝试后, 对计算步骤进行重新安排, 将上述计算步骤调整为如下形式:

$ \begin{align} & {{a}_{k}}=\left[ T\left( {{\mathit{\boldsymbol{Q}}}^{\left( k \right)}} \right),{{\mathit{\boldsymbol{Q}}}^{\left( k \right)}} \right],{{\mathit{\boldsymbol{e}}}_{k}}={{\left\| T\left( {{\mathit{\boldsymbol{Q}}}^{\left( k \right)}} \right)-{{a}_{k-1}}{{\mathit{\boldsymbol{Q}}}^{\left( k \right)}}-{{b}_{k-1}}{{\mathit{\boldsymbol{Q}}}^{\left( k-1 \right)}} \right\|}^{2}},{{b}_{k}}= \\ & \sqrt{{{\mathit{\boldsymbol{e}}}_{k}}-{{\left( {{a}_{k}}-{{a}_{k-1}} \right)}^{2}}}. \\ \end{align} $

这样, a_k, e_k的计算只需要1次全归约, 从而达到减少通讯的目的.

调整步骤后, 得到了改进的SYMMLQ算法.虽然这个过程需增加一些计算量, 但相对很小.同时, 改进的算法结构保持不变, 所以编写程序代码时只需调整原算法相对应的循环过程步骤2即可, 其余步骤保持不变.具体实现过程如下:

步骤2  在p_i中, 计算$\mathit{\boldsymbol{Q}}_{i}^{\left( k \right)}=\frac{\left( H_{i}^{\left( k-1 \right)}-{{a}_{k-1}}\mathit{\boldsymbol{Q}}_{i}^{\left( k-1 \right)} \right)}{{{b}_{k-1}}}, \mathit{\boldsymbol{G}}_{i}^{\left( k \right)}={{\mathit{\boldsymbol{A}}}_{i}}{{\mathit{\boldsymbol{Q}}}^{\left( k \right)}}+\mathit{\boldsymbol{Q}}_{i}^{\left( k \right)}\mathit{\boldsymbol{B}}, \mathit{\boldsymbol{H}}_{i}^{\left( k \right)}=\mathit{\boldsymbol{G}}_{i}^{\left( k \right)}-{{b}_{k-1}}\mathit{\boldsymbol{Q}}_{i}^{\left( k-1 \right)} $和$\mathit{\boldsymbol{L}}_{i}^{\left( k \right)}=\mathit{\boldsymbol{H}}_{i}^{\left( k \right)}-{{a}_{k-1}}\mathit{\boldsymbol{Q}}_{i}^{\left( k \right)} $, 然后计算内积[G_i^(k), Q_i^(k)]和[L_i^(k), L_i^(k)], 进行一次全归约后得到a_k=[G^(k), Q^(k)]和e_k=[L^(k), L^(k)], 计算$ {{b}_{k}}=\sqrt{{{\mathit{\boldsymbol{e}}}_{k}}-{{\left( {{a}_{k}}-{{a}_{k-1}} \right)}^{2}}}$.

3 数值算例与结果分析 3.1 算例

在HPZ80工作站集群并行机上求解3个算例, 分别采用SYMMLQ算法与本文提出的改进的SYMMLQ算法求解, 并对计算结果进行比较.在迭代计算过程中, 所有的初始矩阵X⁽⁰⁾为零矩阵, 终止条件ε=10^-6.

例1    Poisson方程第一边值问题

$ \left\{ \begin{array}{l} \Delta u = \frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}} = - f\left( {x,y} \right),\;\;\;\;0 < x,y < 1\\ u\left| {_{\mathit{\Gamma }{\rm{ = }}0}} \right. \end{array} \right.,f\left( {x,y} \right) = x + y. $

其中, Γ为单位正方形区域的边界, 取步长h=1/1201.采用中心差分格式, 可将本题转化为求解矩阵方程AX+XB=F的问题.设矩阵A=(a_ij)_n×n, B=A, F=(f_ij)_n×n,

$ {a_{ij}} = \left\{ \begin{array}{l} 2,i = j,\\ - 1,\left| {i - j} \right| = 1,{f_{ij}} = {h^3}\left( {i + j} \right).\\ 0,\;\;其他. \end{array} \right. $

当n=1 200时, 计算结果见表 1.

例2   在挡土墙的弹性力学分析中, 由于其结构的形状和受力、约束情况具有一定的特点, 所以将应力分量之间的关系近似表示为平衡偏微分方程，即

$ \frac{{{\partial ^2}{u_1}}}{{\partial {x^2}}} + \frac{{{\partial ^2}{u_2}}}{{\partial {y^2}}} + 2\frac{{{\partial ^2}{u_3}}}{{\partial x\partial y}} = f\left( {x,y} \right). $

其中, 假设正应力u₁，u₂与切应力u₃平衡, 即u₁=u₂=u₃=u, 并且f(x, y)=x+y, 应力边界条件为u|_Γ=0, Γ为受力边界.于是, 平衡偏微分方程转化为

$ \frac{{{\partial ^2}u}}{{\partial {x^2}}} + \frac{{{\partial ^2}u}}{{\partial {y^2}}} + 2\frac{{{\partial ^2}u}}{{\partial x\partial y}} = x + y. $

取步长h=1/1 001, 采用中心差分格式以及一些近似处理, 可将本题转化为求解Sylvester矩阵方程AX+XB=F的问题.设矩阵A=(a_ij)_n×n, B=(b_ij)_n×n, F=(f_ij)_n×n,

$ {a_{ij}} = \left\{ \begin{array}{l} - 4,\;\;\;\;\;\;i = j,\\ 4,\;\;\;\;\;\;\;\;\left| {i = j} \right| = 1,\\ 1,\;\;\;\;\;\;\;\;\left| {i = j} \right| = 2,\\ 0,\;\;\;\;\;\;\;\;其他. \end{array} \right.\;\;{b_{ij}} = \left\{ \begin{array}{l} - 8,\;\;\;\;\;\;i = j,\\ 2,\;\;\;\;\;\;\;\;\left| {i = j} \right| = 1,\\ - 1,\;\;\;\;\;\;\left| {i = j} \right| = 2,\\ 0,\;\;\;\;\;\;\;\;其他. \end{array} \right.{f_{ij}} = 3{h^3}\left( {i + j} \right), $

计算结果见表 2.

例3   设矩阵A=(a_ij)_n×n, B=(b_ij)_n×n, 其中

$ {a_{ij}} = \left\{ \begin{array}{l} 1,9,\;\;\;\;\;i = j,\\ - 1,\;\;\;\;\;\;\left| {i = j} \right| = 1,\\ 0,\;\;\;\;\;\;\;\;其他. \end{array} \right.\;\;\;{b_{ij}} = \left\{ \begin{array}{l} 1.8,\;\;\;\;\;i = j,\\ - 1,\;\;\;\;\;\left| {i = j} \right| = 1,\\ 0,\;\;\;\;\;\;\;其他. \end{array} \right. $

求解Sylvester矩阵方程AX+XB=F.在此取n=1 500, F为任意给定的矩阵.计算结果见表 3.

表 1~3中p表示处理机台数, T/s表示时间, K表示迭代次数, E表示相对并行效率, E₁表示绝对并行效率, Δ表示误差‖R^(k)‖.由于一台处理机的计算时间较长, 因此采用多台处理机的计算时间, 并且使用多台处理机中最小的计算时间作为基准来计算加速比和并行效率.

3.2 结果分析

(1) 改进的SYMMLQ算法比未改进的算法减少了并行计算时间, 这是由于改进后的算法在循环迭代中减少了全归约的次数, 从而缩短了全局通讯所引起的时间消耗.

(2) 改进的SYMMLQ算法的迭代次数与原SYMMLQ算法的迭代次数基本相同, 说明改进的算法没有破坏原有SYMMLQ算法的结构与数值稳定性, 而且由表中误差数据可以看出本文改进的算法对计算的精度影响不大.

(3) 表 3数据表明当计算规模一定时, 处理机增加到30台以上, 并行效率会急剧下降, 这是由于整体内积的计算与全收集而引起的大规模通讯使得消耗增加, 进一步说明本文算法与原算法的可扩放性比较差, 因此为保证较高的并行效率应采用适当数量的处理机.

致谢: 感谢西北工业大学高性能计算研究与发展中心的大力支持!