计算科学 - 为一个非常病态的矩阵设计一个预条件子 - 吾爱随笔录

我是一位数值方法知识有限的物理学家，我正试图加快一个非常病态问题的反演（ $rcond>10^{30}$ ）。重复使用相同的稀疏方阵，因此我们当前的设置依赖于一次性 LU 分解和三角求解器。

我有一个非常强大的 GPU 可供我使用，并且想使用它（将 CUDA 与 Cusparse、Magma 一起使用……）。我目前的理解是，三角形求解器本质上是非常连续的，在并行 SIMD 架构上不能很好地工作，但另一方面，迭代求解器在 GPU 上工作得很好，遗憾的是我的问题远非条件恶劣。

由于矩阵使用了大量时间，因此花费一些时间来计算预处理器似乎是要走的路。因此，我正在寻找设计适合并行求解器的预处理器的方法。

如果你们想看看它，我已经包含了一个指向我正在使用的矩阵的链接。

Big Matrix
size=114708*114708
one-based_numbering
34Mb

小矩阵
大小=12890*12890
one-based_numbering
3Mb

还有一个用 Matlab 加载矩阵的脚本：

n=12890;
%n=114708;
fileID = fopen('SmallMatrix.txt','r');
%fileID = fopen('BigMatrix.txt','r');
dataArray = textscan(fileID,'%f%f%f%[^\n\r]', 'Delimiter',' ', 'MultipleDelimsAsOne', true);
fclose(fileID);
U1 = dataArray{:, 1};
U2 = dataArray{:, 2};
V = dataArray{:, 3};
s=sparse(U1,U2,V,n,n);
figure;
spy(s);

编辑：尽管条件数非常高，但我使用的直接求解器（matlab、mumps、umfpack 和 superlu）解决它没有问题。在matlab中，

cond(s)
x=rand(n,1);
y=s\x;
norm(s*y-x)

退货，

1.4996e+34
2.7318e-10

Edit2：这是我的矩阵可能具有的结构（非常大的值总是在对角线上）：

(\begin{matrix} 0.1 & 2.1 & 3.2 \\ 4.1 & - 1.1 & 21 \\ 8 & - 0.2 & 10^{30} \end{matrix})

$\begin{pmatrix} 0.1 & 2.1 & 3.2 \\ 4.1 & -1.1 & 21\\ 8 &-0.2 & 10^{30} \end{pmatrix}$

和 RHS：

(\begin{matrix} 3.1 \\ 4 \\ 0.3 \end{matrix})

$\begin{pmatrix} 3.1 \\ 4 \\ 0.3 \end{pmatrix}$ 解决方案是：

(\begin{matrix} 1.3544 \\ 1.4117 \\ 0 \end{matrix})

$\begin{pmatrix} 1.3544 \\ 1.4117 \\ 0 \end{pmatrix}$

这是非常病态的，但每个直接求解器都完全可以接受。这只是施加边界条件的一种简单方法（我不得不承认，不是很优雅）

Edit3：按照 Brian Borchers 和 Tobias 的建议，我已经从我的矩阵和 RHS 中删除了我的非常大的值出现的所有行和列。从某种意义上说，它的效果很好，给定 RHS 的解决方案是相同的，并且条件数更合理 $10^{5}$ . 这是用于执行此操作的脚本：

%Creating the new matrix
sCopy=s;
nRemoved=nnz(abs(sCopy)>10^29);
for index=1:n-nRemoved
    while abs(sCopy(index,index))>10^29
        sCopy(index,:)=[];   
        sCopy(:,index)=[];  
    end
end
nnz(abs(sCopy)>10^29)

%Updating the  RHS
[i,j,]=find(abs(s)>10^29);
NewRHS=RHS;
NewRHS(i)=[];

有了这个，我的主要询问仍然存在，条件数仍然太高，迭代方法无法收敛。