I am a physicist with limited numerical methods knowledge and I am trying to speed up the inversion of a very ill-conditioned problem ($rcond>10^{30}$). The same sparse square matrix is used repeatedly, so our current setup relies on a one time LU factorisation and a triangular solver.
I have at my disposal a very powerful GPU and would like to use it (using CUDA with Cusparse, Magma...). My current understanding is that triangular solver are very sequential by nature and do not work well on parallel SIMD architecture, but iterative solver on the other hand do work well on GPU, sadly my problem is far to ill-condtionned.
Since the matrix is used a large number of time, taking some time to compute a preconditioner seems like the way to go. I am therefore looking for ways to design a preconditioner suitable for parallel solver.
If you guys want to take look at it, i have included a link to the kind of matrix I am working with.
Big Matrix
size=114708*114708
One-based_numbering
34Mb
Small Matrix
size=12890*12890
One-based_numbering
3Mb
Also a script to load the matrix with 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);
Edit: Despite the very high condition number, the direct solvers I have used (matlab, mumps,umfpack and superlu) have no problem solving it. In matlab,
cond(s)
x=rand(n,1);
y=s\x;
norm(s*y-x)
Returns,
1.4996e+34
2.7318e-10
Edit2: Here is the kind of structure my matrix might have (Very big value are always on the diagonal):
\begin{pmatrix} 0.1 & 2.1 & 3.2 \\ 4.1 & -1.1 & 21\\ 8 &-0.2 & 10^{30} \end{pmatrix}
And the the RHS: \begin{pmatrix} 3.1 \\ 4 \\ 0.3 \end{pmatrix} And the solution is: \begin{pmatrix} 1.3544 \\ 1.4117 \\ 0 \end{pmatrix}
It is very Ill-conditionned but every direct solver is completely okay with it. It is just a easy way to impose boundary conditions (not very graceful, i have to admit)
Edit3: Following Brian Borchers and Tobias advice, I have removed from my matrix and my RHS all the rows and columns where my very large values appear. It works great in the sense that the solution to a given RHS is the same and the condition number is a more reasonable $10^{5}$. Here is the script used to do that:
%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)=[];
With that out of the way, my main interrogation remains, the condition number is still too high for a iterative method to converge.