说明书对带假结的RNA二级结构预测进行硬件加速的方法
技术领域
本发明涉及一种对基于四维动态规划方法的带假结的RNA二级结构预测进 行加速的方法，目的是加快带假结的RNA二级结构预测的速度。
背景技术
RNA二级结构是识别ncRNA的重要依据，是研究RNA功能的基础和前提。实 验手段是获取RNA二级结构的最可靠方法，目前主要的RNA结构测定方法有X射 线衍射和核磁共振，虽然采用实验方法获取的结果精确可靠但是其过程非常耗 时，并且代价高昂，因此研究RNA结构预测计算方法就显得特别重要，近年来 采用计算机和数学模型预测RNA序列二级结构的方法被广泛采用，成为RNA研究 领域的热点问题。RNA二级结构预测方法一般包含以下三个部分：
(一)几何表示方法
由于构成配对的碱基对RNA二级结构的稳定性起促进作用，而由未配对碱 基构成的各种环都是破坏RNA结构稳定性的，所以RNA二级结构预测的核心是寻 找序列中的配对碱基，通常采用插入“·”的嵌套括号图“·(·(··)··(·)·)·” 表示碱基配对情况，其中左括号“(”、右括号“)”和“·”都对应序列中的碱 基，配对的左右括号表示对应位置上的碱基构成互补配对，而“·”表示序列 中对应位置的碱基构成环结构。图1为RNA序列二级结构示意图。
(二)打分函数
通过试验测定和对已知结构RNA序列的统计分析，对相邻碱基对和碱基对 与独立碱基之间相互作用因素实现参数化，并采用打分函数为各种可能的RNA 二级结构进行打分，从而评价预测结果的好坏。
(三)搜索策略
RNA二级结构预测不是一个对各种可能的结构进行穷举的过程，需要选择 一个优化方法对结构空间进行快速搜索，迅速找到与某一全局极大(小)值相 对应的二级结构。
定义1：假设R为长度为n的RNA序列，R＝r1r2r3...rn，(i·j)表示RNA序列R中的 碱基ri和rj构成互补配对，i，j，k，l分别表示碱基ri、rj、rk、rl在RNA 序列中的序号，且1≤i≤j≤n，1≤k≤l≤n，则RNA二级结构预测问题实际上 是寻找当打分函数y＝f(g(x1),g(x2),...,g(xi))取全局极大(小)值情况下，序列R 中碱基对的集合S，其中f为组合函数，xi表示子序列r1r2...ri,(1≤i≤n)。
目前存在两类主要的RNA二级结构预测方法：第一种方法是从头预测法， 该方法以单条RNA序列作为输入。Nussinov算法是最早提出的基于单序列的RNA 结构预测算法，该算法通过寻找具有最大碱基配对个数的结构来实现结构预 测，因此也称为最大碱基配对算法。由于该方法仅考虑了配对碱基对RNA二级 结构稳定性所起的作用，算法的预测精度较差。
由于碱基配对可以是RNA分子的能量降低，结构趋于稳定，因此最小自由 能算法(Minimum Free Energy，简称MFE)认为在一定的温度下，RNA分子通 过构象调整达成某种热力学平衡，使自由能最小，从而形成最稳定的状态，此 时的二级结构即被认为是RNA的真实二级结构。最小自由能算法由M.Zuker于 1981年提出，又被称为Zuker算法。该算法的计算对象不是简单的碱基配对数 量，而是子序列的自由能。算法基本思想基于RNA二级结构中各种子结构的自 由能具有独立性和可加性假设，采用试验方法测定的各种子结构自由能参数 表，将序列所有可能形成的子结构的自由能相加，整条RNA序列的最小自由能 等于所有可能的子结构能量之和的最小值。Zuker算法是目前最好的针对单条 RNA序列的结构预测算法，尤其是针对小分子RNA的结构预测取得了很好的预测 结果，该方法的缺点是不支持包含假结的RNA二级结构的预测。
基于随机上下文无关文法(SCFG)模型的预测方法也是针对单条RNA序列 的典型结构预测方法，是目前最适于描述和建模RNA二级结构的全概率模型， 在RNA二级结构预测研究领域占有重要地位。目前基于SCFG理论模型的标准比 对算法为Coche-Younger-Kasami，简称CYK算法[18][19]。CYK算法用于实现单 条序列与单个RNA家族的共变模型(covariance model，简称CM模型)进行比 对，从而判断该RNA序列是否属于该家族并且进一步地得到该序列的二级结构。 CYK算法虽然以单序列作为输入，但家族的共变模型的建立需要大量的RNA序列 来进行参数估计。
以上方法都属于单序列从头预测方法，随着基因组测序技术的发展，已知 的RNA序列也越来越多，这使得利用比较基因组方法预测RNA二级结构成为可 能。这类方法以多条同源RNA序列或由它们构成的比对作为输入，其理论基础 是生物序列的结构保守性大于序列保守性。同源比较方法基于多序列比对，首 先利用多序列比对工具如ClustalW程序构造RNA多序列比对，然后通过突变 检测获取该组序列的保守结构。同源预测法的典型代表是RNAalifold，该算 法是MFE方法在RNA多序列比对中的扩展，它在计算该组序列平均最小自由能 的同时考虑了协变信息(covariance information)，通过能量计算和协变分 值计算相结合的方法预测该组同源序列的公共二级结构。
上述算法无论是基于单序列的还是基于多序列同源比对的都限制了碱基 对(i·j)和(k·l)相互间的位置关系，即满足i<k<j<l或者k<i<l<j，而没有考 虑碱基构成的交错结构，因此都不能够预测假结。而假结对病毒基因组的复制 和蛋白质合成调控起到非常重要的作用[34]，是RNA三级结构的重要构成要 素，因此对假结的预测成为目前RNA二级结构预测领域的热点问题。由于基于 MFE模型的带假结的RNA二级结构预测已被证明是NP-完全问题。为了提高算法 的实用性，研究者根据实验观测数据加强假结预测的约束条件，降低了算法复 杂度，使得通过计算方法预测假结成为可能。目前，已有几种支持假结预测的 近似算法。1999年Rivas和Eddy首次采用动态规划算法实现了对RNA假结的预 测，算法的时间复杂度为O(n6)，空间复杂度为O(n4)。文献通过约束假结的类 型将计算复杂度降低为O(n5)，进一步将计算复杂度降低为O(n4)，但只能预测 最简单的假结，实用性不高。由于Rivas和Eddy提出的算法能够很好的预测平 面假结和受限制的非平面假结，是目前公认的最完整、最权威的支持假结的RNA 二级结构预测算法。
上述算法尽管采用了不同的RNA二级结构几何表示方法和打分函数，但 是都采用了相同的搜索策略，即采用动态规划思想将整个序列的结构预测问题 分解为一系列子序列的结构预测问题，从最短子序列出发逐步得到整个序列结 构的最优解。按照绪论中提出的动态规划问题分类标准，常规的结构预测属于 三维动态规划问题，而带假结的RNA结构预测属于四维动态规划问题。
前面讨论的是不包含假结的RNA二级结构预测，而实际上假结也是一种常 见的RNA二级结构类型，它对RNA三级结构的构成以及在影响RNA的功能活性方 面起到非常重要的作用。
定义2：在定义1的条件(3)中，如果(k,l)∈S,若满足i<k<j<l或者 k<i<l<j，则碱基对(i·j)和碱基对(k·l)构成的交错结构就称为假结 (Pseudoknots)。图2为RNA二级结构中的假结及对应的圆顶图表示。
因此，在定义2的基础上放宽条件(3)中对碱基序号的限制，即对序列 R＝r1r2r3...rn中的碱基ri、rj、rk、rl，允许出现当i<k<j<l或者k<i<l<j时，(i·j) 和(k·l)的交错配对，则普通RNA二级结构预测问题即拓展为带假结的RNA 二级结构预测问题。虽然相对于常规的碱基配对数量来说，假结在总体数量上 很少，但假结是RNA三级结构的重要构成要素，近年来假结预测能力逐渐成 为衡量RNA二级结构预测算法性能的重要指标。
在所有的支持假结预测的算法中，Rivas和Eddy提出的PKNOTS算法能 够很好的预测平面假结和受限制的非平面假结，是目前公认的最好的支持假结 预测的RNA二级结构预测算法，其预测结果已得到实验检验。
PKNOTS算法是基于MFE模型的单序列从头预测方法，也采用了动态规 划思想，将序列结构预测问题分解为较短的子序列的结构预测问题，通过计算 子序列的结构得到序列本身的结构。由于每一个假结需要ijkl四个参数来描 述，所以动态规划的递归关系中需要四重循环来寻找最优，因此称为四维动态 规划问题。PKNOTS算法的输入是单条RNA序列，输出是包含假结(如果存 在)的碱基配对结果。
由于包含假结的结构预测问题具有极高的计算复杂性，不适合处理规模较 长的序列。测试表明，在AMD Phenom 9650Quad CPU上使用PKNOTS-1.05 程序对长度为64bps的序列实现预测需要140s，而长度为128bps的序列结构 预测的时间超过9000s，约2.5小时。目前虽有不少基于PKNOTS的改良算法 能够支持假结预测，但是这些算法均已牺牲精度和正确性换取执行速度，预测 效果不佳；此外还有其他一些假结预测方法，如基于堆叠稳定的权重匹配算法， 虽然具有理想的计算复杂度，但只对特定类型的假结具有较好的预测结果。比 较而言，PKNOTS算法的预测效果明显优于其他算法，但是高昂的时空复杂 度限制了PKNOTS算法的实用性，目前只能对包含几十个碱基的短序列结构 实现预测。
2010年，Krishnan等人基于IBM Cell多核处理器，以PKNOTS算法为对 象首次开展了针对四维动态规划算法的并行化研究，与运行在通用微处理器上 的标准PKNOTS软件相比，在Sony Play Station3平台上运行的并行版本获得 了大约3倍的加速效果，但只能支持长度小于100bps的RNA序列。本节基于 FPGA平台，对带假结的RNA二级结构预测领域涉及的四维动态规划问题展 开研究，针对高维动态规划问题复杂的计算特征提出一种数据相关分析方法， 并在此基础上实现存储优化和细粒度并行，相对于现有串行算法实现获得了 3～5倍的全局加速效果。
PKNOTS算法的计算过程涉及三个二维矩阵VX、WX、WBX和四个四维 矩阵VHX、ZHX、YHX和WHX，算法的时间复杂度为O(n6)，空间复杂度为 O(n4)，n为序列长度。其中四维矩阵VHX、ZHX和YHX的迭代公式如下：
VHX ( i , j : k , l ) = min EIS 2 ( i , j : k , l ) EIS 2 ( i , j : r , s ) + VHX ( r , s : k , l ) EIS 2 ( r , s : k , l ) + VHX ( i , j : r , s ) WHX ( i + 1 , j - 1 : k - 1 , l + 1 ) + M - - - ( 1 ) ]]>
ZHX ( i , j : k , l ) = min VHX ( i , j : k , l ) + P ZHX ( i , j : k - 1 , l ) + Q ZHX ( i , j : k , l + 1 ) + Q ZHX ( i , j : r , l ) + WX ( r + 1 , k ) ZHX ( i , j : k , s ) + WX ( l , s - 1 ) EIS 2 ( i , j : r , s ) + ZHX ( r , s : k , l ) WHX ( i + 1 , j - 1 : k , l ) + P + M - - - ( 2 ) ]]>
YHX ( i , j : k , l ) = min VHX ( i , j : k , l ) + P YHX ( i + 1 , j : k , l ) + Q YHX ( i , j - 1 : k , l ) + Q YHX ( r + 1 , j : k , l ) + WX ( i , r ) YHX ( i , s : k , l ) + WX ( s + 1 , j ) EIS 2 ( r , s : k , l ) + YHX ( i , j : r , s ) WHX ( i , j : k - 1 , l + 1 ) + P + M - - - ( 3 ) ]]>
上述公式中的变量i,j,k,l,r,s表示碱基在RNA序列中的序号，同时也是元 素在矩阵中的坐标，满足i≤r≤k≤l≤s≤j的关系。VHX、ZHX和YHX 为四维动态规划矩阵，P、Q、M和EIS2为能量参数。和Zuker算法计算公式 相比，包含假结的PKNOTS算法由于考虑了假结这种特殊子结构，因此增加 候选结构的搜索空间，从公式形式上看便是每个元素的能量计算都相应地增加 了候选分支，但是基本思想仍基于最小能量模型，计算过程中的每一步都是从 所有可能的子结构对应的能量值中选取最小值作为局部结构，整条RNA序列 的最小能量等于所有可能的子结构能量之和的最小值这一基本原理依然未变。
图3为四维动态规划算法中的计算空间图。如果采用二维平面图表示四维 矩阵，根据元素下标间的关系可以发现PKNOTS算法中，每一个四维矩阵的 高两维i和j构成一个二维上三角矩阵，并且其中的每一个单元(称为Cell) 都是一个二维三角矩阵。
通过对迭代公式进一步分析可以发现PKNOTS算法中的三个二维矩阵 VX、WX、WBX和四个四维矩阵VHX、ZHX、YHX和WHX之间存在远比 三位动态规划问题复杂的自我调用和相互调用关系。图4是矩阵间的数据依赖 关系图，图中的圆圈代表矩阵，箭头代表存在数据依赖关系，自圈代表存在自 身数据依赖，WHX矩阵也有自圈，为了简洁起见在图中没有标出。
由于带假结结构预测算法是在基于最小自由能模型的常规结构预测算法 上的扩展，并且都采用了动态规划方法实现由区列局部到全局能量的计算，因 此与前文所述的RNA二级结构预测算法具有类似的计算特征：(1)基本的 计算区域为三角形；(2)计算过程存在大量参数表查询操作；(3)存在行列交 替和不规则存储访问特征，数据相关距离随着计算位置的移动而变化。此外， 由于带假结的结构预测算法放松了对配对碱基的位置限制，搜索空间的增加导 致存储调度和IO带宽成为设计瓶颈，而计算复杂度的增加给数据相关的分析 和并行策略的制定带来挑战。
上述算法尽管采用了不同的RNA二级结构几何表示方法和打分函数，但是 都采用了相同的搜索策略，即采用动态规划思想将整个序列的结构预测问题分 解为一系列子序列的结构预测问题，从最短子序列出发逐步得到整个序列结构 的最优解。
调研结果表明，目前在国内外对序列分析领域的硬件加速研究大多限于序 列一级结构层次，目前没有对高维尤其是四维结构预测算法实现硬件加速的报 道。
发明内容
本发明的目的就是针对现有方法的缺陷，首次提出了一种对基于四维动态 规划方法的带假结的RNA二级结构预测进行加速的方法，目的是加快带假结的 RNA二级结构预测的速度。
本发明的技术方案是这样实现的：它是先构建由主机和可重构算法加速器 组成的异构计算系统，接着主机将格式化后的热力学模型参数和编码后的 RNA序列发送至可重构算法加速器，算法加速器的七个计算模块采用MPMD 方式执行无回溯的PKNOTS算法计算；计算中采用矩阵降维方法将四维矩阵 分解为N个三维矩阵，然后采用逐层按区域轮转分割和区域内按列并行处理 的任务划分策略实现细粒度并行计算，每一个计算模块内部的n个PE采用 SPMD方式同时计算位于区域内不同列的n个数据，n为自然数。
PKNOTS算法计算通过三个二维矩阵和四个四维矩阵计算模块实现，所 述的三个二维矩阵计算模块为PE_VX、PE_WX和PE_WBX；所述的四个四 维矩阵计算模块为PE_WHX、PE_VHX、PE_ZHX和PE_YHX。
所述的二维矩阵计算模块PE_VX、PE_WX和PE_WBX的结构完全相同， 它们内部均包含一个子PE控制器、一个子PE计算单元、一个局部存储器和 一个数据传递寄存器，其中子PE控制器用于实现对计算和数据访存时序的控 制；子PE计算单元的核心是一个32bit加法器，用于实现对两个输入操作数 的加法运算，其计算结果同时写入局部存储器和数据传递寄存器，局部存储器 用于缓存二维矩阵一整列元素的计算结果，数据传递寄存器仅存储当前元素的 计算结果，并供下一个计算模块立即使用。
所述的四维矩阵计算模块PE_WHX、PE_VHX、PE_ZHX和PE_YHX的 结构完全相同，四维矩阵计算模块的功能是实现对二维Cell的并行计算，每 个Cell都是一个二维三角矩阵。
本发明以RNA序列分析领域中的复杂结构预测算法对高性能计算的需求 为背景，基于通用微处理器结合FPGA硬件算法加速器的异构体系结构，从提 取典型方法的动态计算特征入手，研究对复杂数据相关和不规则存储访问的优 化方法，对典型算法实现细粒度并行，达到高效加速计算的目的；并在此基础 上提出了一种基于可重构硬件的对带假结的RNA二级结构实现预测的细粒度 并行算法，为特定领域的算法群提供一种基本的硬件结构模板和并行程序设计 框架，为有效降低算法加速器设计复杂度、实现加速器快速生成奠定基础。
本发明公开的四维动态规划矩阵细粒度并行计算方法，不仅能够加快带假 结的RNA二级结构预测的速度，而且该方法以及基于异构体系结构下的并行程 序设计模板和设计框架不仅能指导硬件加速器快速生成，还能够为解决其他领 域的高维动态规划矩阵的计算问题提供技术参考。
附图说明
图1为RNA序列二级结构示意图
图2为RNA二级结构中的假结及对应的圆顶图
图3为四维动态规划算法中的计算空间图
图4为矩阵间的数据依赖关系图
图5为时空域重叠数据相关性分析方法流程图
图6为二维矩阵计算模块内部结构图
图7为四维矩阵的计算过程图
图8为三维矩阵的计算过程图
图9为四维矩阵计算模块内部的线性阵列图
图10为四维动态规划算法并行计算结构图
具体实施方式
以下结合实施例对本发明作进一步描述：
本发明先构建由主机和可重构算法加速器组成的异构计算系统，接着主机 将格式化后的热力学模型参数和编码后的RNA序列发送至可重构算法加速 器，算法加速器的七个计算模块采用MPMD方式执行无回溯的PKNOTS算法 计算；计算中采用矩阵降维方法将四维矩阵分解为N个三维矩阵，然后采用 逐层按区域轮转分割和区域内按列并行处理的任务划分策略实现细粒度并行 计算，每一个计算模块内部的n个PE采用SPMD方式同时计算位于区域内不 同列的n个数据，n为自然数。
PKNOTS算法计算通过三个二维矩阵和四个四维矩阵计算模块实现，所 述的三个二维矩阵计算模块为PE_VX、PE_WX和PE_WBX；所述的四个四 维矩阵计算模块为PE_WHX、PE_VHX、PE_ZHX和PE_YHX。七个模块(PE) 的并行计算结构，如图10所示。每个PE负责一个矩阵的计算。模块的命名 采用PE_后接所计算矩阵名称的方式进行命名(例如：模块PE_VX表示当前 PE负责完成矩阵VX的计算)。
图10中的七个模块(PE)构成一个PE阵列，PE阵列每次同时计算七个 矩阵中具有相同下标的单元(i,j)，对二维矩阵VX、WX和WBX而言是一个 元素，而对四维矩阵WHX、VHX、YHX、ZHX而言则是图7中的一层(即 一个三角形区域)，称为一个“Cell”。
前三个二维矩阵计算模块PE_VX、PE_WX和PE_WBX的结构完全相同， 其内部结构如图6所示。二维矩阵计算模块PE_VX、PE_WX和PE_WBX的 内部均包含一个子PE控制器(Sub PE Controller)、一个子PE计算单元 (Sub_PE)、一个局部存储器(Mem)和一个数据传递寄存器(Trans Regs)。 各组件之间的连接关系如图6所示，箭头表示数据传递方向。其中子PE控制 器(Sub PE Controller)用于实现对计算和数据访存时序的控制；子PE计算单 元的核心是一个32bit加法器，用于实现对两个输入操作数的加法运算，其计 算结果同时写入局部存储器(Mem)和数据传递寄存器。其中局部存储器 (Mem)用于缓存二维矩阵一整列元素的计算结果，数据传递寄存器仅存储 当前元素的计算结果，并供下一个计算模块立即使用。
后四个四维矩阵计算模块PE_WHX、PE_VHX、PE_ZHX和PE_YHX的 结构完全相同，其内部结构如图9所示。每个计算模块的内部都包含一个线性 PE阵列。由于四维矩阵计算模块的功能是实现对二维“Cell”的并行计算，而每 个Cell都是一个二维三角矩阵，因此使用了图9所示的多PE线性阵列结构， 采用“按列轮转划分”策略实现对多列元素的并行计算。图9中构成PE阵列的 所有子处理单元的结构完全相同，其内部结构与图9所示的子模块结构相同。 子PE控制模块(Sub PE Controller)实现任务分配，每次将三角矩阵中一列元 素的计算任务加载至对应的子处理单元(Sub PE)，并控制阵列的同步。计算 启动后，阵列中的每个子PE单元每次计算当前列中的一个元素，这样将整个 阵列作为一个整体来看，就实现了对图8(b)子图某个区域中的一条对角线 的同步计算。随着每个子PE当前所计算单元的向上位移，整个PE阵列当前 所计算的对角线也就随之向上移动，这样也就逐步实现了对二维“Cell”的并行 计算。
本发明采用了一种称为“时空域重叠”的数据相关性分析方法，通过程序特 征分析生成操作执行顺序表，并从中提取数据相关性；通过对多单元执行顺序 表中的项目在时间和空间域上合并建立数据源到目的地的映射关系(通过时间 域重叠发现不相关的操作，实现并行计算；通过空间域重叠发现相同数据源， 实现数据重用)，并构建访存调度矩阵，指导数据调度优化和并行策略的制定。
图5为“时空域重叠”数据相关性分析方法主要流程，它显示了从源代码 到数据载入与传递示意图的生成过程。该过程包括运算类型和数据源统计、数 据相关和操作数来源分析、多单元执行顺序表合并和生成访存调度矩阵四个步 骤。
1.运算类型和数据源统计
对每个矩阵对应的处理单元，按照软件中代码的执行顺序列出操作类型和 数据来源；如果遇到循环语句则对循环进行展开，以四维矩阵中的“Cell”为 基本数据块，对循环变量的变化规律和数据相关区域进行统计和分析，画出当 前元素和计算所依赖元素的移动轨迹，提取数据相关性；按照数据相关性对代 码中的操作执行按顺序进行编号，列出操作类型、数据来源和相关序号，生成 如表1所示的源代码中的操作执行顺序表。表1中操作类型一栏既可以表示一 个操作也可以表示循环体中的一组类型相同的操作。操作数来源一项中的 YHX(1,1)表示当前操作依赖于YHX矩阵Cell(1,1)中的数据，Para(1,1)表示当 前操作依赖于参数表(1,1)区域中的数据。相关序号一栏表示当前操作的计算结 果将会被序号对应的后续操作使用，即存在序号之间存在先写后读相关。如果 本栏为空则表示该操作不依赖于之前操作的运算结果。
表1操作执行顺序表
2.数据相关和操作数来源分析
表1反映了代码的串行执行过程，通过对表1进行分析可以获得操作之间 的数据依赖关系。接下来按照满足数据真相关条件下的执行顺序对操作进行重 新标注，对存在数据相关的操作分配顺序的编号，对不存在数据相关的操作分 配相同的编号，然后对编号相同的操作进行合并。由于编号相同的操作意味着 可以并行执行，也就是说该组操作的真实执行时间是相同的，因此将该步骤称 为时间重叠。接下来对执行时间重叠操作的数据来源进行分析，对操作数来源 一栏中相同的项进行合并。同一行中不同操作的数据来源相同意味着真实的访 存地址相同或者相近，访问的数据属于同一个Cell，因此可以看作是地址空间 重叠的存储访问。经过上述两个处理步骤得到如表2所示的时空重叠处理后的 执行顺序表。最后对顺序表操作数来源一栏相邻两行之间相同的项进行标注。
表2时空重叠处理后的执行顺序表
比较表1和表2可以发现，由于表1中的前两项不存在数据相关，因此将 其合并；并且两组加法操作的操作数来源部分重叠，都依赖于YHX矩阵的 Cell(1,1)因此将其合并。观察表2操作数来源一栏中的相邻表项可以发现， YHX(1,1)和WHX(1,1)都被相邻的操作使用，因此对其进行标注。同一数据源 如果被相邻操作所使用意味着在存储调度时可以考虑对其进行片内缓存，通过 数据重用减少片外存储访问开销。
3.多单元执行顺序表合并
按照步骤1和2对每一个处理单元建立操作执行顺序表，然后对多个表格 进行合并，生成表多单元执行顺序表。表3并列列出了三个计算模块的执行顺 序表，每个模块都包括数据来源和表示数据相关的序号。
表3多单元执行顺序表
由于不同计算模块同时计算不同矩阵中坐标相同的元素，并将结果存储在 FPGA片内，如果这些元素之间存在数据依赖关系则可通过片内数据传递网络 实现数据重用，而不会涉及片外存储调度问题，因此本步骤不考虑计算模块之 间的数据相关。因此，表2中同一计算模块的纵向操作之间存在先后关系，而 同一序号的横向操作可以并行执行。
接下来在保证每个模块纵向相对执行顺序不变的前提下，再次执行时间重 叠操作，对计算模块的执行顺序进行上下调整，使得处于不同计算模块但数据 来源相同的操作尽量位于表格的同一行。
表4时间重叠后的多单元执行顺序表
比较表3和表4可以发现，由于计算模块2的第一个操作和计算模块1 的第三个操作都要使用ZHX(1,1)，所以将表3中计算模块2的第一行移至表4 的第三行，因此计算模块2在前两个执行时间段处于空闲状态。基于同样的原 因，计算模块3的第二个操作被移至表4的第四行，与计算模块2的第二个操 作对齐。由于计算模块3的第一个操作和计算模块1的第一个操作都要使用 YHX(1,1),Para(1,2)，所以计算模块3第一行的位置保持不变。
4.生成访存调度矩阵
在表4的基础上再次进行数据源合并，首先按照执行序号将每一行所有计 算模块需要使用的数据源进行整理，将相同的数据来源进行横向合并，按照数 据源的载入顺序生成访存调度矩阵；其次，考虑执行时间上相邻的访存操作， 如果使用了相同的数据源则进行纵向合并，通过FPGA内部的数据缓存实现数 据重用，避免重复载入，最终生成最终的访存调度矩阵。表5为访存调度矩阵 的示意图，表格纵向为按载入顺序排列的数据源地址，横向为数据传输的目的 地，“1”表示左侧数据会被对应的计算模块使用，“0”表示左侧数据不会被对 应模块使用。“●”表示对应的数据块已经载入FPGA，并处于有效状态，不再 需要从片外载入。
表5访存调度矩阵
数据来源 计算模块1 计算模块2 计算模块3
YHX(1,1) 1 0 1
Para(1,2) 1 0 1
Para(1,1) 1 0 0
VHX(1,1) 0 0 1
YHX(1,1)● 1 0 0
WHX(1,1) 1 0 0
WHX(1,1)● 1 0 0
ZHX(1,1) 1 1 0
Para(1,1)● 0 1 0
Para(1,2)● 0 1 0
ZHX(1,1)● 0 1 1
VHX(1,1)● 0 1 1
YHX(1,1)● 0 1 0
ZHX(1,1)● 0 1 1
WHX(1,1)● 0 0 1
… … … …
生成访存调度矩阵时需要考虑以下几个因素：(1)数据相关性，如果当前 计算所依赖的不同数据源存储在不同的存储模块中，则从不同通道同时载入； (2)如果不同数据源存储在同一个存储模块中，则按使用先后顺序载入，尽 快启动流水线；(3)如果数据源的使用不存在相关性，则先载入大块数据；(4) 如果IO通道空闲并且FPGA内部有空闲的缓冲区则立即预取下一数据块。
实验结果表明，使用最终生成的访存调度矩阵指导片外存储访问调度、数 据分配和重用能够减少50％左右的访存请求，从而有效降低存储访问开销。
由于二维矩阵VX、WX和WBX的存储需求小、计算过程简单，因此将 其存放在FPGA片内，设计时不用过多考虑它们的计算和存储问题，而四维矩 阵WHX、VHX、YHX、ZHX的计算是PKNOTS算法的核心。按照图4所示 的矩阵间数据依赖关系，WHX矩阵处于数据依赖图的核心，因此本节以WHX 矩阵为例说明四维矩阵的填充过程。
四维三角矩阵WHX(i,j,k,l)可以被分解为N个三维三角矩阵WHXi(j,k,l) (1≤i≤N)，每个三维矩阵WHXi(j,k,l)由N个边长为N的二维三角矩阵构 成，每个二维上三角矩阵对应于图3中的一个Cell。
如图7所示，计算过程以Cell为基本单位，从WHX1的第1层(Cell1) 开始：当矩阵WHX1(j,k,l)的第1层WHX1(1,k,l)算完后再计算第2层WHX1(2,k,l)，直到WHX1的最后一层算完；然后再计算第2个三维矩阵WHX2的 第1层，第2层，…,第N层；接下来再计算WHX3…，直到最后一个矩阵 WHXN的第N层算完。图中的虚线和二维三角矩阵右上角标识的数值代表Cell 的计算顺序。
对四维矩阵中的每个三维矩阵WHXi(j,k,l)的计算，采用图8所示的区 域划分计算策略，将每一层(Cell)按列分割为若干区域，然后逐个区域进行 计算。图8(a)子图中的三维矩阵矩阵每层都被分为三个区域，区域的编号 代表对应区域的计算顺序，带箭头的虚线表示每个区域中元素的计算顺序。对 每个区域的计算，按照图8中(b)子图所示的按列轮转划分的任务分配策略， 使用多个处理单元按照由下至上的顺序沿矩阵对角线实现对当前Cell中当前 区域的并行计算。
对图8(b)子图所示的每个区域，每个PE负责计算当前区域中的一列， PE计算的元素的列号与PE在阵列中的序号一一对应。图中填充阴影的p列 元素表示当前的计算区域，它们被同时分配至p个PE上并行计算。每个PE 都从各自列的底部位置开始，按照由下至上的顺序计算。当前区域的计算启动 时，所有PE计算的元素都位于三角矩阵的主对角线上(图中标星号的单元表 示初始计算位置)：PE_1计算元素(k,l)，第2个PE计算元素(k+1,l+1)，…， 第p个PE计算元素(k+p-1,l+p-1)。根据算法的数据相关性，对角线上的元素 不存在数据相关，因此处于不同PE上的p个元素可以并行计算。而且由于位 于同一条对角线上元素的计算量相等，因此所有PE都可以同步向上推进，在 任意时刻，PE阵列当前计算的元素总是处于矩阵同一条对角线上。由于三角 矩阵每一列元素的个数不等，当PE计算的元素行坐标k＝1时，PE计算暂停， 进入等待状态(如果计算结果需要写回片外存储器，PE将在等待状态下发出 写回请求)。阵列中的所有PE将按编号顺序依次进入同步等待状态，并发出 写回请求。
对其他三个四维矩阵VHX、YHX和ZHX采用和WHX相同的计算顺序。 而对于三个二维矩阵VX、WX和WBX，则采用从左至右逐列计算、每列由 下至上的顺序实现填充。为了实现并行计算，本文设计了七个计算模块(PE)， 每个PE负责一个矩阵的计算。PE阵列每次同时计算七个矩阵中具有相同下 标的单元(i,j)，对二维矩阵VX、WX和WBX而言是一个元素，而对四维矩 阵WHX、VHX、YHX、ZHX而言则是图7中的一层，即一个“Cell”。
图10为基于异构多PE线性阵列的四维动态规划算法并行计算结构，主 要由阵列控制模块、多PE计算阵列、存储模块和阵列同步和写回控制模块构 成。其中阵列控制模块负责计算阵列的初始化、任务分配和控制计算区域的切 换。
计算阵列由七个PE模块构成，分别实现对七个自由能矩阵的计算，所有 PE的地位是对等的，都和数据总线相连。每个计算模块都拥有独立的输入和 输出数据缓存(Data Buf和Cache)，其中Data Buf用于缓存从片外载入的数 据，Cache用于存储本模块的计算结果，所有PE的计算逻辑和输出数据缓存 都通过数据传递网络连接，通过数据共享方式实现重用。整个计算阵列的数据 缓存和PE独有的数据Cache都使用FPGA片内多端口BlockRAM存储块实现。 为了避免访问冲突，每个计算模块都保存一份RNA序列和自由能参数表的拷 贝，使用FPGA片内分布式存储资源实现。此外，PE之间还设计了数据传递 寄存器组，实现PE计算结果的快速传递。阵列同步和写回控制模块与每个PE 的输出缓存和运算逻辑相连，用于控制PE阵列的同步和将Cache中保存的计 算结果依次写回片外存储器。
对比试验：我们在测试平台上实现了硬件PKNOTS算法加速器。测试平 台由一台通用计算机和一个算法加速器构成。主机配置为Intel Core2四核 Q94002.66GHz处理器，4.0GB主存。算法加速器硬件主要包括1片Xilinx Virtex7系列FPGA芯片(XC7VX485T)，三条容量为8GB的DDR3-1600DRAM 存储条，加速器通过SFP+光纤数据通道与主机相连(采用XC7VX485T芯片 内部集成的GTX Transceiver实现)，有效数据传输带宽可达10Gb/s。算法加 速器支持动态重构，可在60ms内完成不同规模的CM模型间的快速切换，与 配置时间为秒级的常规配置方法如JTAG或并行SlectMAP相比，FPGA的配 置效率提高了2～3个数量级。RNA二级结构预测软件版本为PKNOTS-1.08， 由美国华盛顿大学医学院Elena Rivas开发，分别在Intel Core2四核Q9400、 Intel Xeon(R)X5670CPU以及FPGA算法加速器三种不同平台上运行。
实验结果(表6)表明，在XC7VX485T FPGA平台上只能实现一个PKNOTS 算法加速引擎，主要原因是缓存四维矩阵的“Cell”数据块占用了过多的存储 容量，存储资源的使用率达到了82％，而逻辑资源使用率仅为28％。由于运 算类型主要是乘加操作，设计中不存在大规模多路选择器和集中式的存储访问 端口，系统时钟频率可以达到210MHz，可见存储资源不足是系统实现的主要 瓶颈。如果使用目前最大规模的商用FPGA器件XC6VSX1140T可以实现至少 2个PKNOTS算法加速引擎，可以同时实现2条RNA序列的结构预测，并且 能够支持更长的序列。
表6四维动态规划算法在FPGA平台上的实现结果
并行效果
表7PKNOTS算法加速效果(时间单位：秒)
实验选择了4组长度在30～176bps之间的RNA序列，测试了PKNOTS-1.08 程序在Inter Q9400和Intel Xeon(R)X5670CPU平台下的平均执行时间，并与 硬件加速器进行了比较。从表7可以看到，在算法加速器上执行包含30个碱 基的序列结构预测可获得2倍的加速比，而当测试序列长度为176bps时，可 获得51.8倍的加速效果。与Intel Xeon(R)X5670相比也可获得超过25倍的加 速效果。受限于XC7VX485T FPGA器件的逻辑和存储容量，目前还只能对长 度小于256bps的序列实现带假结的结构预测。使用Xilinx EDA工具的综合结 果表明，在XC7VX1140T芯片上可以实现2个PKNOTS加速引擎，同时对2 条序列实现结构预测，相对于目前主流CPU平台，可获得超过60倍的加速效 果。