CN103631761B

CN103631761B - 并行处理架构进行矩阵运算并用于严格波耦合分析的方法

Info

Publication number: CN103631761B
Application number: CN201210313665.9A
Authority: CN
Inventors: 刘志钧; 徐益平; 施耀明
Original assignee: Raintree Scientific Instruments Shanghai Corp
Current assignee: Raintree Scientific Instruments Shanghai Corp
Priority date: 2012-08-29
Filing date: 2012-08-29
Publication date: 2018-02-27
Anticipated expiration: 2032-08-29
Also published as: CN103631761A

Abstract

为提供速度较快并且成本较低的矩阵运算和严格波耦合分析的技术，本发明提供了一种并行处理架构进行矩阵运算并用于严格波耦合分析的方法，该并行处理架构包括多个处理器模块，各处理器模块分别与独立的共享存储器关联并执行独立的线程块。矩阵运算的方法包括如下步骤：i.各处理器模块分别执行并行处理架构调用中的运算指令，其中，各运算指令一一对应于矩阵运算中的各运算部分，各运算部分能够并行执行并且互不相关；ii.将该运算部分所用的数据分别读入相应处理器模块的共享存储器中；iii.各处理器模块基于相应的运算指令，读取共享存储器中的相应数据，并行地执行线程块来完成该矩阵运算中的该运算部分。

Description

并行处理架构进行矩阵运算并用于严格波耦合分析的方法

技术领域

本发明涉及并行计算，特别涉及基于并行处理架构进行矩阵运算。

背景技术

随着高端计算机图形显示卡的发展，多核图形处理单元(Graphic ProcessingUnit，简称GPU)越来越强大，GPU不仅为显示图像做了优化，还有天生的并行性。经过对硬件和软件的改进，GPU的可编程力不断提高，在计算上已经超越了通用的CPU。为了充分利用芯片的强大的计算功能，显卡厂商NVIDIA推出了新的运算平台-统一计算设备架构(ComputeUnified Device Architecture，简称CUDA)。CUDA是一种通用并行处理架构，该架构使GPU能够解决复杂的计算问题。它包含了CUDA指令集架构(ISA)以及GPU内部的并行计算引擎。CUDA采用C语言作为编程语言，开发人员可以使用C语言来为CUDA架构编写程序，能够在GPU的强大计算能力的基础上建立起一种效率更高的密集数据计算解决方案。

CUDA支持大量的线程并行，并在硬件中动态地创建，调度和执行这些线程。CUDA编程模型将CPU作为主机，而将GPU作为协处理器，以CPU来控制程序整体的窜行逻辑和任务调度，而让GPU来运行一些能够被高度线性化的数据并行部分。

能够使用GPU计算的程序必需有以下特点：需要处理的数据量比较大，数据以数组或矩阵形式有序存储，并且对这些数据要进行的处理方式基本相同，各个数据间的依赖性或耦合很小。

CUDA提供了一种编译工具(nvcc)。开发人员只要将含CUDA指令的原程序，以.cu作文件后缀。nvcc会将.cu文档拆解出在GPU上执行的部份，及在CPU上执行的部份，并调用适当的程序进行编译动作。在GPU执行的部份会透过NVIDIA提供的编译器编译成中介码，而主机执行的部份则会透过系统上的C++编译器来编译(在Windows上使用Visual C++而在Linux上使用gcc)。

所有目前支持CUDA的NVIDIA显示芯片，都是由多个多核处理器(或称multiprocessors)组成，如图1所示。每个多核处理器里包含了八个流处理器，其组成是四个一组，也就是说实际上可以看成是有两组、每组四个处理器。此外，每个多核处理器通常具有8192个寄存器，16KB～64KB的共享存储器，以及纹理缓存和固定缓存。此外，显示芯片还有显卡内存。每个流处理器都可以读写显卡内存，但只能读写它所在的多核处理器中的寄存器和共享存储器。流处理器读写寄存器和共享存储器很快，但读写显卡内存较慢。所以，编程时，应尽量多用共享存储器，少用显卡内存。

发明内容

半导体芯片制造过程中，产品成品率是衡量芯片制造工艺的重要指标。为提高成品率，芯片制造过程中，需要用到光学关键尺寸(Optical Critical Dimension，OCD)检测和缺陷检测。

光学关键尺寸(OCD)的测量，通过获取被测区域周期性结构的散射信号以及结构的模型从而估计出结构的形貌参数。OCD测量原理总体上可概括为两个步骤：光谱获取过程--获取样品的散射信号并处理为测量光谱；光谱匹配过程-根据样品的形貌模型寻找特定的形貌参数使其对应的理论光谱与测量光谱实现最佳匹配。

随着集成电路工艺技术进入45纳米之后，技术路线图向32纳米以下技术节点挺进时，图形密度不断增加、关键尺寸不断微缩，工艺控制窗口非常狭小，以前可以被忽略的缺陷现在可能导致器件不能正常工作，成为影响成品率的致命缺陷。检测方法通常有成品检测，中间过程检测等。中间过程检测要求较高，要求快速和无损。光学成像检测能满足这些要求。光学成像检测是用宽带组合光源照射电路。为了增强缺陷信号强度，提高信噪比，需要通过对入射光束有针对性的控制和对散射场作有针对性的选择滤波来实现优化。

光学检测中无论光学关键尺寸(OCD)测量或光学成像缺陷检测都离不开严格精确的电磁场模拟计算。在这一领域，数值仿真计算代表性的方法有：严格波耦合分析理论(Rigorous Coupled-Wave Analysis，简称RCWA)。

如图2，所示，设介质在x，y方向呈周期性变化。z方向通常情况下并非不变，光刻掩模板通常在z方向上均匀，或者z方向上分成几层，每层内均匀不变。晶片上的微细结构通常在z方向变化，但严格波耦合分析方法在z方向将介质划分若干薄片.薄片的厚度如果足够小，则可认为光散射特性方面在z方向介质分布均匀。这样，整个介质的光散射效果可以看成若干个叠加在一起的z方向介质分布均匀的薄片的光散射效果。求解出每个介质薄片上平面和下平面处的电磁场分布就可以得出整个介质的光散射仿真。

这里仅以TE平面波垂直入射二维光栅为例，对RCWA算法作简单介绍。如图2，结构分为三层，I，光栅上层空气层，II，光栅层，III，光栅下衬底层。

在第一和第三层中，

其中， (a1)中，第一项是入射场部分。入射光从(θ，φ)方向入射，其中，θ是入射光与z轴的夹角，φ是入射面与x-z面的夹角，如图2。在二维情况，φ＝0。

光栅层通常在z方向有变化，但严格波耦合分析方法在z方向将介质划分若干薄片.薄片的厚度如果足够小，则可认为光散射特性在z方向介质分布均匀。在薄片中，

将Maxwell方程中的介电常数作Fourier展开，并解一特征值问题，得

其中，是特征值问题的一特征向量，γ_q为对应的特征值。在薄片与薄片间的分界面，光栅与空气层的分界面及光栅与衬底层的分界面上匹配切向电磁场，可获得矩阵方程组。解矩阵方程组可得散射矩阵方程：

[R]＝[S][I](5)

其中[R]是各模式反射系数R_m组成的矢量，[I]是各入射光束的模I_m组成的矢量。[S]是散射矩阵。解上述散射矩阵方程就可获得某一特定入射光束的散射结果。图3是计算OCD一条光谱的流程图。缺陷检测中的理论计算与OCD光谱计算稍有不同，用(5)式计算从各Littrow-Mounting入射方向入射的光产生的各模式反射系数。Littrow-Mounting入射方向(θ_i，φ_i)与几何周期结构的周期有关。设几何结构在x和y方向的周期长度分别为Λ_x和Λ_y，(θ_i，φ_i)必须满足kcosθ_isinφ_i＝2πn_i/Λ_x，和kcosθ_icosφ_i＝2πm_i/Λ_y。其中，k是波数，整数对(n_i，m_i)i＝1，2，…，M又必须满足传播条件k²＞(2πn_i/Λ_x)²+(2πmⁱ/Λ_y)²。

以上的计算中有矩阵的运算如计算矩阵与向量的乘积，矩阵与矩阵相乘，矩阵求逆和求矩阵的特征值与特征向量等。在模拟计算中，计算量很大。为了获得较快的计算速度，满足生产需要，目前一般为光学检测设备配备功能强大的工作站或服务器来完成RCWA计算，这提高了光学检测的成本。因此，提供一种所需硬件较为简单的矩阵运算方法，对大规模矩阵运算，特别是光学检测中的RCWA运算是十分有利的。

本发明的发明构思在于，使用例如CUDA等并行处理架构来进行矩阵运算，将矩阵运算中能够并行执行并且互不相关的运算部分在并行处理架构上分别同时进行，缩短计算时间。并且，具有CUDA架构的显卡的价格相较工作站或服务器来说十分低廉，因此使得实现成本较低。

根据本发明的一个方面，提供了一种使用并行处理架构进行矩阵运算的方法，该并行处理架构包括多个处理器模块，各处理器模块分别与独立的共享存储器关联并执行独立的线程块，该方法包括如下步骤：

i.各处理器模块分别执行并行处理架构调用中的运算指令，其中，各运算指令一一对应于矩阵运算中的各运算部分，各运算部分能够并行执行并且互不相关；

ii.将该运算部分所用的数据分别读入相应处理器模块的共享存储器中；

iii.各处理器模块基于相应的运算指令，读取共享存储器中的相应数据，并行地执行线程块来完成该矩阵运算中的该运算部分。

该方面的优点在于，提供了在并行处理架构上进行矩阵运算的技术方案，使得矩阵运算的运算速度相比单线程来说要快很多。

根据一个进一步的实施方式，所述矩阵运算包括矩阵-列向量的乘法，

所述步骤i中，所述各运算指令包括计算两个向量的内积；

所述步骤ii中，分别将该矩阵的一行向量与该列向量读入处理器模块的共享存储器中；

所述步骤iii中，各处理器模块分别执行该矩阵的该行向量与该列向量的内积，并将各内积依次组成该矩阵-列向量乘法的结果向量。

本实施方式提供了对于矩阵-列向量的乘法的具体的并行计算方案。

根据一个更具体的实施方式，所述处理器模块包括多个处理器，每个处理器执行独立的线程，所述步骤i中，各运算指令包括：

-各个线程分别将一个向量中的一个位置的元素与另一个向量中的相同位置的元素相乘，得到多个积；

-各个线程分别将多个积中的相邻积不重复地相加，得到多个和；

-各个线程分别将多个和中的相邻和不重复地相加，并基于所得到的和重复该步骤，直至得到最终的单个和。

该实施方式提供了对矩阵-列向量的乘法中的向量内积的并行计算方案，其中各个线程并行地计算向量中元素的积，并且并行地计算元素的积的和，提高了矩阵运算的运算速度。

根据一个进一步的实施方式，所述矩阵运算包括两个矩阵的乘积，

所述步骤i中，所述各运算指令包括计算分块矩阵的乘积及矩阵和；

所述步骤ii中，分别将该两个矩阵的各分块矩阵读入处理器模块的共享存储器中；

所述步骤iii中，各处理器模块分别计算该两个矩阵的分块矩阵的乘积及矩阵和，并将结果组成该两个矩阵的乘积矩阵。

本实施方式提供了对于两个矩阵的乘积的具体的并行计算方案。

根据一个更具体的实施方式，所述处理器模块包括多个处理器，每个处理器执行独立的线程，该分块矩阵的维度的平方不大于线程块中的最大线程数，并且，

将两个矩阵的最末列的分块矩阵在列方向进行零填充至该维度，和将两个矩阵的最末行的分块矩阵在行方向进行零填充至该维度。

在该分块矩阵的维度的平方不大于线程块中的最大线程数的情况下，一个处理器模块独自能够完成分块矩阵的相乘，不需要引入其他处理器模块以及引入显卡内存，这样具有较高的效率。对分块矩阵进行零填充使得分块矩阵的大小一致，能够直接计算，不需要再对矩阵的维度进行检测判断。

根据一个进一步的实施方式，所述矩阵运算包括矩阵求逆，

所述步骤i中，各运算指令包括使用第j行来对另一行做初等行变换，消去另一行中的第j列的元素；

所述步骤ii中，将该矩阵的第j行和剩余行中的分别各行读入处理器模块的共享存储器中；

所述步骤iii中，各处理器模块并行地使用该矩阵的第j行来对剩余行中的分别各行做初等行变换，消去另一行中的第j列的元素；

步骤i至iii分别依次对该矩阵的每一行执行，各行的执行可由不同的线程块同时完成。

本实施方式提供了对矩阵求逆的具体的并行计算方案。

根据一个进一步的实施方式，所述矩阵运算包括矩阵的QR分解，

所述步骤i中，各运算指令包括使用一行来消去另一行中的相应列的元素；

所述步骤ii中，将该矩阵的一行和另一行读入处理器模块的共享存储器中；

所述步骤iii中，各处理器模块使用该矩阵的一行来消去另一行中的相应列的元素；

步骤i至iii循环进行，在每一次循环中，并行地分别对于前i-1个元素都已消去的各第i行，步骤i至iii使用该行来消去第j行中第i列的元素，其中，j大于i且前j-1行中第i列元素都已被消去。

本实施方式提供了对矩阵QR分解的具体的并行计算方案。

根据一个进一步的实施方式，所述并行处理架构包括基于具有多核图形处理单元的图形显示卡的统一计算设备架构CUDA。由于基于CUDA的显卡价格较低，因此该实施方式的实现成本较低。

根据本发明的另一个方面，提供了一种使用并行处理架构计算矩阵的特征值的方法，该方法包括如下步骤：

-使用前述的方法对矩阵进行QR分解，获得Q矩阵和R矩阵；

-使用前述的方法计算R矩阵和Q矩阵的乘积，得到新的矩阵；

-判断新矩阵中下三角非对角元的最大模数；

-当该最大模数大于给定值时，使用新的矩阵重复以上步骤；当该最大模数小于给定值或重复次数超过给定次数时，以该新的矩阵的对角元为该矩阵的特征值。

本方面提供了对计算矩阵的特征值的具体的并行计算方案。

根据本发明的又一个方面，提供了一种使用并行处理架构进行严格波耦合分析的方法，该方法包括如下步骤：

-接收几何结构参数和入射光束参数和波长参数；

-使用前述的方法来计算光栅区的本征模；

-匹配边界条件；

-使用前述的方法来解矩阵方程组得到散射矩阵；

-使用前述的方法来计算入射光束产生的所有散射模系数。

根据本发明的又一个方面，提供了一种将严格波耦合分析用于光学关键尺寸测量的方法，该方法包括如下步骤：

-根据接收波长范围和波长步长，确定各波长；

-使用前述的方法来计算在各波长的入射光束产生的所有散射模系数。

根据本发明的又一个方面，提供了一种将严格波耦合分析用于缺陷检测的方法，该方法包括如下步骤：

-接收待测器件的三维几何结构参数和入射光源的光谱波长参数；

-使用前述的方法来计算三维结构的本征模；

-匹配边界条件；

-使用前述的方法来解矩阵方程组得到散射矩阵；

-根据接收的结构参数，决定θ和φ最大张角的阶次M；

-使用前述的方法来计算各入射张角方向入射光束产生的所有散射模系数。

附图说明

通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更加明显：

图1是根据基于多个多核处理器的图形显示卡的CUDA硬件架构的示意图；

图2是一种典型的二维光栅周期结构；

图3是用严格波耦合分析的方法计算一条OCD光谱的流程图；

图4是并行地计算多个数的和的、树状加流程图；

图5是基于本发明所提出的基于并行处理架构来用严格波耦合分析的方法计算一条OCD光谱的流程图。

具体实施方式

本发明提供了一种使用并行处理架构进行矩阵运算的方法，该并行处理架构包括多个处理器模块，各处理器模块分别与独立的共享存储器关联并执行独立的线程块，该方法包括如下步骤：

i.各处理器模块分别执行CUDA调用中的运算指令，其中，各运算指令一一对应于矩阵运算中的各运算部分，各运算部分能够并行执行并且互不相关；

其中，“能够并行执行并且互不相关”的运算部分是指一个运算部分的执行的结果不会对其他运算部分的执行产生影响，例如该运算部分的执行得到的标量、向量或矩阵不会作为其他运算部分的输入，下面将通过多个例子来对这种运算部分进行说明。可以理解，这种运算部分并不限于以下所举的例子，取决于各种矩阵运算的过程，本领域的一般技术人员分析得到其中“能够并行执行并且互不相关”的运算部分，并采用本发明的实施方式来进行并行计算。

以下将分别说明本发明所提供的方法实现各种具体的矩阵运算的技术方案，其中，以CUDA作为并行处理架构为例进行说明，可以理解，其他并行处理架构也可以被用于本发明的实施方式中。

一、矩阵与向量的乘积

矩阵-向量乘法是将一个n×n阶方阵A＝[a_ij]乘以n×1的向量B＝[b₁，b₂，...，b_n]^T得到一个具有n个元素的列向量C＝[c₁，c₂，...，c_n]^T。

在计算前，A矩阵和B向量都先复制到显卡内存中。

矩阵-向量乘法以行划分矩阵，矩阵A的每一行与B向量的相乘在一个处理器模块的同一线程块内完成，能充分利用处理器模块内的共享存储器。将B向量读入每个处理器模块内的共享存储器中，且将矩阵A的各行一一读进相应处理器模块的共享存储器中。另外，再在共享存储器中，声明一与B向量同样长度的向量D，存储A矩阵的行向量与B向量相乘的结果。例如，在第i个线程块内，第j个线程读入的指令做如下操作：

d[j]＝a[i][j]*b[j](6)

并且，在做完乘法后，优选地利用树状加法来计算各个积的和，如图4。第一步，各线程分别将向量D中相邻的元素不重复地相加；第二步，各线程分别将第一步得到的相加的和中相邻的和不重复地相加，...，直到做完。将最终的和存入c[i]中。向量[c[0]，......，c[n]]为所求出的矩阵-向量的乘积。

显卡上的内存是DRAM，因此最有效率的存取方式，是以连续的方式存取。考虑到实际上线程的执行方式。当一个线程在等待内存的数据时，GPU会切换到下一个线程。也就是说，实际上执行的顺序是类似线程0-＞线程1-＞线程2-＞...

因此，在同一个线程中连续存取内存，在实际执行时反而不是连续了。要让实际执行结果是连续的存取，应该要让线程0读取第一个数字，线程1读取第二个数字...依此类推。在一个线程块内，要做矩阵行与向量的乘积。因此，在将矩阵A复制到显卡内存中时，应以行为序，使矩阵的同一行元素连续存放在一起。

具体并行算法框架描述如下：

矩阵-向量乘CUDA并行算法

输入：A_n×n，B_n×1

输出：C_n×1

假设一次乘法运算时间为一个单位时间，不难得出基于CUDA的矩阵-向量乘算法的并行计算时间：若处理器个数和向量维数相当，则其时间复杂度为O(n)。

二、矩阵与矩阵的乘积

一个m×n阶矩阵A＝[a_ij]乘以一个n×k的矩阵B＝[b_ij]就可以得到一个m×k的矩阵C＝[c_ij]，它的元素c_ij为A的第i行向量与B的第j列向量的内积。矩阵相乘的关键是相乘的两个元素的下标要满足一定的要求(即对准)。

A矩阵和B矩阵都先复制到显卡内存中。如简单从显卡内存中读取A矩阵和B矩阵的元素进行相乘，需要访问显卡内存2mnk次。

为减少访问显卡内存的次数，矩阵乘法用分块法，即分块矩阵乘法(Block MatrixMultiplication)，也就是把整个矩阵乘法的动作，切割成很多小矩阵的乘法。例如，

要计算C矩阵的子块C₁₁，可以把块矩阵的运算当成一般矩阵的运算：

C₁₁＝A₁₁B₁₁+A₁₂B₂₁+A₁₃B₃₁(8)

这样一来，我们就可以把两个小矩阵加载到共享内存，则小矩阵本身的乘法就不需要再存取任何外部的内存了。假设小矩阵的大小是p×p，则实际上需要的内存读取次数就会变成约2mnk/p。

由于目前CUDA每个线程块的线程数目最多是512，因此，如p取16，小矩阵块的元素数量为256。一个线程块内就能完成一个矩阵块的相乘。能充分利用共享内存。理论上，这样应该可以让读取效率提高16倍(假设没有遇到别的瓶颈)。

在程序中，因为矩阵的大小不一定会是16的倍数，如果需要使用if判断式检查是否超出矩阵范围，那么会降低运行效率。要把那些if判断式去掉，有一个方法是，在配置内存时，就配置成p＝16的倍数，并在复制矩阵到显卡内存之前，先将它清为0。也就是，矩阵A和矩阵B最右侧的一列的子块，如他们的列维数不足16，补到16，所补部分，都设为零。同样，A和矩阵B最下面的一行的子块，如他们的行维数不足16，补到16，所补部分，也都设为零。

将矩阵A按行划分为u×v块，u＝[m/p]+1，v＝[n/p]+1。这些矩阵块依次记为A_ij(i＝0，1，…，u，j＝0，1，…，v)。将矩阵B按行划分为v×w块，矩阵块依次记为B_ij(i＝0，1，…，v，j＝0，1，…，w)。

C矩阵划分为u×w个子块。每个线程块处理A矩阵的一行子块和B矩阵的一列子块的相乘，并将结果存入矩阵C相应的子块中，如(8)式。

具体CUDA并行算法框架描述如下：

矩阵并行分块乘法算法流程

输入：A_m×n，B_n×k，

输出：C_m×k

三、矩阵求逆

矩阵求逆(Matrix Inversion)是一常用的矩阵运算。对于一个n×n阶的非奇异方阵A＝[a_ij]，其逆矩阵是指满足A^-1A＝AA^-1＝I的n×n阶方阵，其中I为单位方阵。

矩阵求逆的过程中，主程序利用一循环，依次用各主行i(i＝0，1，...，n-1)对其余各行j(j≠i)作初等行变换，消去该行中第i列的元素。由于各行计算之间没有数据相关关系，各线程块可以分别为其余各行独立作初等行变换。

由于各线程块不能作同步，用a_kk消去a_jk(j＝k+1，…n)，可在一次CUDA调用中完成。每一线程块完成一次行初等变换。总共需要n次调用。具体算法框架描述如下：

矩阵求逆的并行算法流程：

输入：矩阵A_n×n，

输出：矩阵A^-1 _n×n

四、矩阵的QR分解

H＝[a_ij]为一个n阶矩阵，对H进行QR分解，就是求一个非奇异方阵Q与上三角方阵R，使得H＝QR。其中方阵Q满足：Q^H＝Q^-1，称为正交矩阵，因此QR分解又称为正交三角分解。

由于QR分解中消去h_ij时，同时要改变第i行及第j行两行的元素，而在LU分解中，仅利用主行i(i＜j)变更第j行的元素。因此QR分解并行计算中对数据的划分与分布与LU分解就不一样。以第i列为例，需要消去h_ij(i＜j)，每消去一个元素，都要改变第i行的元素，并用改变后的第i行去消去下一元素。所以消去同一列中的不同元素，不能并行计算。

但，改变第i行及第j行两行的元素，可以用一个线程块进行。在h_ij消去后，在消去第i行其它元素时，第j行不会再改变。可用i+1行开始消去i+1列的元素。这可用另一线程块完成。

具体步骤如下：

1，用一个线程块消去h_i，i+1

2，用一个线程块消去h_i，i+2

3，用一个线程块消去h_i，i+3，用另一个线程块消去h_i+1，i+2

4，用一个线程块消去h_i，i+4，用另一个线程块消去h_i+1，i+3

5，用第一线程块消去h_i，i+5，用第二线程块消去h_i+1，i+4，用第三线程块消去h_i+2，i+3

……

具体并行算法框架描述如下：

矩阵QR分解并行算法流程：

输入：矩阵H_n×n，单位矩阵Q

输出：矩阵Q_n×n，矩阵R_n×n

五、矩阵特征值求解

对给定的A₀＝A∈C^n×n，QR算法的基本迭代格式如下：

A_m-1＝Q_mR_m，

m＝1，2，...，

A_m＝R_mQ_m，

(9)

其中Q_m为酉矩阵，R_m为上三角阵.为了下面的理论分析方便起见，我们这里暂且要求R_m的对角元都是非负的.由(7)容易推出

即矩阵序列{A_m}中的每个矩阵都与A相似.由(10)可得

其中将A_m＝Q_m+1R_m+1代入上式即有

从而有

即

其中由此即知

反复进行这一过程，即得到矩阵序列：A₁，A₂，...，A_m，A_m+1，...，它们满足如下递推关系：A_i＝Q_iR_i；A_i+1＝R_iQ_i(i＝1，2，...，m，...)其中Q_i均为正交阵，R_i均为上三角方阵。这样得到的矩阵序列{A_i}或者将收敛于一个以A的特征值为对角线元素的上三角矩阵，形如：

或者将收敛于一个特征值容易计算的块上三角矩阵。

并行QR分解求矩阵特征值的思想就是反复运用并行QR分解和并行矩阵相乘算法进行迭代，直到矩阵序列{A_i}收敛于一个上三角矩阵或块上三角矩阵为止。具体的并行算法描述如下：

QR方法求一般矩阵全部特征值的CUDA并行算法流程

输入：矩阵A_n×n，单位矩阵Q，ε

输出：矩阵特征值Eigenvalue

-接收几何结构参数和入射光束参数和波长参数；

-使用前述的方法来计算光栅区的本征模；

-匹配边界条件；

-使用前述的方法来解矩阵方程组得到散射矩阵；

-使用前述的方法来计算入射光束产生的所有散射模系数。

-根据接收波长范围和波长步长，确定各波长；

-使用前述的方法来计算三维结构的本征模；

-匹配边界条件；

-使用前述的方法来解矩阵方程组得到散射矩阵；

-根据接收的结构参数，决定θ和φ最大张角的阶次M；

具体的用严格波耦合分析的方法计算一条OCD光谱的过程如图5所示。其中，计算光栅区的本征模、解矩阵方程组得到散射矩阵以及计算入射光束产生的所有散射模系数的技术是本领域的一般技术人员所熟知的。本发明的改进包括使用前面提到的使用并行处理架构来进行矩阵运算的方法来实现严格波耦合分析中的各个矩阵部分。例如，在计算光栅区的本征模时，使用到前面已经描述的并行的矩阵QR分解方法及其他方法，在解矩阵方程组时，使用前面已经描述的并行的矩阵求逆方法及其他方法，在计算散射模系数时，使用前面已经描述的并行的矩阵与向量积和矩阵与矩阵乘积方法。在实际使用中，在开始计算前，可以对GPU进行设置，在计算结束后，可以关闭GPU的设置。

当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。例如，对于其他矩阵运算，该运算中的互不相关的运算部分能够根据本发明被多个处理器模块并行地执行。

本领域普通技术人员可以理解上述方法中的全部或部分步骤可通过程序来指令相关硬件完成，所述程序可以存储于计算机可读存储介质中，如只读存储器、磁盘或光盘等。可选地，上述实施例的全部或部分步骤也可以使用一个或多个集成电路来实现。相应地，上述实施例中的各模块/单元可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。本发明不限制于任何特定形式的硬件和软件的结合。

Claims

1.一种使用并行处理架构进行矩阵运算并用于严格波耦合分析的方法，该方法包括如下步骤：

-接收几何结构参数和入射光束参数和波长参数；

-采用所述并行处理架构进行矩阵运算来计算光栅区的本征模；

-匹配边界条件，包括匹配切向电磁场，以获得矩阵方程组；

-采用所述并行处理架构进行矩阵运算来解所述矩阵方程组得到散射矩阵；

-采用所述并行处理架构进行矩阵运算来计算入射光束产生的所有散射模系数，

其中，所述并行处理架构包括多个处理器模块，各处理器模块分别与独立的共享存储器关联并执行独立的线程块，其中采用所述并行处理架构进行矩阵运算包括如下步骤：

2.根据权利要求1所述的方法，其特征在于，采用所述并行处理架构进行的所述矩阵运算包括矩阵-列向量的乘法，

所述步骤i中，所述各运算指令包括计算两个向量的内积；

3.根据权利要求2所述的方法，其特征在于，所述处理器模块包括多个处理器，每个处理器执行独立的线程，所述步骤i中，各运算指令包括：

4.根据权利要求1所述的方法，其特征在于，采用所述并行处理架构进行的所述矩阵运算包括两个矩阵的乘积，

5.根据权利要求4所述的方法，其特征在于，所述步骤ii中，所述处理器模块包括多个处理器，每个处理器执行独立的线程，该分块矩阵的维度的平方不大于线程块中的最大线程数，并且，

将两个矩阵的最末列的分块矩阵在列方向进行零填充至该维度，将两个矩阵的最末行的分块矩阵在行方向进行零填充至该维度。

6.根据权利要求1所述的方法，其特征在于，采用所述并行处理架构进行的所述矩阵运算包括矩阵求逆，

所述步骤ii中，将该矩阵的第j行和剩余行中的各行分别读入处理器模块的共享存储器中；

所述步骤iii中，各处理器模块并行地使用该矩阵的第j行来对剩余行中的各行分别做初等行变换，消去另一行中的第j列的元素；

对该矩阵的每一行分别依次执行步骤i至iii，各行的执行可由不同的线程块同时完成。

7.根据权利要求1所述的方法，其特征在于，采用所述并行处理架构进行的所述矩阵运算包括矩阵的QR分解，

8.根据权利要求1所述的方法，其特征在于，所述并行处理架构包括基于具有多核图形处理单元的图形显示卡的统一计算设备架构。

9.一种使用并行处理架构计算矩阵的特征值的方法，该方法包括如下步骤：

-使用根据权利要求7所述的方法对矩阵进行QR分解，获得Q矩阵和R矩阵；

-使用根据权利要求4所述的方法计算R矩阵和Q矩阵的乘积，得到新的矩阵；

-判断新矩阵中下三角非对角元的最大模数；

10.一种将严格波耦合分析用于光学关键尺寸测量的方法，该方法包括如下步骤：

-根据接收波长范围和波长步长，确定各波长；

-使用根据权利要求1的方法来计算在各波长的入射光束产生的所有散射模系数。

11.一种将严格波耦合分析用于缺陷检测的方法，该方法包括如下步骤：

-使用根据权利要求1至9中任一项所述的方法来计算三维结构的本征模；

-匹配边界条件；

-使用根据权利要求1至9中任一项所述的方法来解矩阵方程组得到散射矩阵；

-根据接收的结构参数，决定θ和φ最大张角的阶次M，其中θ是入射光与z轴的夹角，φ是入射面与x-z面的夹角，z轴指示所述待测器件的表面的法线方向，x轴指示所述待测器件的表面中的三维结构的周期性方向；

-使用根据权利要求1至9中任一项所述的方法来计算各入射张角方向入射光束产生的所有散射模系数。