3
$\begingroup$

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.

Improve this question
$\endgroup$
16
  • 3
    $\begingroup$ With a condition number of $10^{30}$ compared to machine epsilon of about $10^{-16}$ for double precision, the results you get are basically going to be amplified numerical noise regardless of how you solve it. $\endgroup$
    – Nick Alger
    Nov 22, 2014 at 23:57
  • 3
    $\begingroup$ Any method of solution (using double precision floating point arithmetic) is hopeless. You should be looking at strategies for regularizing the problem. $\endgroup$
    – Brian Borchers
    Nov 23, 2014 at 0:11
  • 3
    $\begingroup$ If that's the only reason for the ill-conditioning, then just set all of those variables equal to 0 and substitute them into the other equations and look at the condition number of the resulting system of equations. It may be fine. $\endgroup$
    – Brian Borchers
    Nov 23, 2014 at 3:41
  • 2
    $\begingroup$ Note that having a small residual norm is not necessarily equivalent to having an extremely accurate solution- that's the nature of ill-conditioned systems of equations. $\endgroup$
    – Brian Borchers
    Nov 23, 2014 at 3:42
  • 2
    $\begingroup$ Note, that in your case you can easily avoid the ill-condition by balancing the matrix (already row-scaling is suffient in your case). As Brian already noted this just leads to zeroing the unknowns with the huge diagonal elements. (@BrianBorchers) Brian should write that as an answer to your question and you should accept it. $\endgroup$
    – Tobias
    Nov 23, 2014 at 10:13

1 Answer 1

Reset to default
0
$\begingroup$

I found a solver that seems to be suitable for my problem. My very ill-conditioned matrix is actually the FEM discretisation of a Poisson equation on a bidimentionnal unstructured mesh. First results with a geometrical multi-grid approach (multiple V-cycle) show very good weak-scaling on CPU. And since most of the time is spent multiplying big sparse matrices with vectors (on the fine mesh), I expect the GPU performance to be okay.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.