Information about complex matrix to be solved

Open discussion regarding features, bugs, issues, vendors, etc.

Information about complex matrix to be solved

Postby thanhvanchi » Wed Dec 28, 2016 5:12 am

Hello all member!

I am trying to solve Ax = b system where A is a square, highly ill-conditioned, unsymmetric, partly dense, partly sparse, complex matrix. The linear system arises out of a physical problem where the exact solution to the problem is known and thus the accuracy of the solution x can be verified. I have been able to solve this system using ZGELSS with a fixed RCOND = -1.0D0 so that LAPACK uses machine precision as the SVD threshold.
As I cannot predict what the condition number of A is going to be in advance, choosing RCOND = -1.0D0 seemed a good idea and the results with this setting in ZGELSS are satisfactory. The only problem with ZGELSS is that it is very slow.
I was aware that ZGELSY will be faster as it uses QR decomposition but do not know how to choose RCOND for ZGELSY. I did try ZGELSY for my linear system and found that the results depend greatly on the input value of RCOND. So the questions are:

a/How do I choose RCOND if I do not know the condition number of A in advance for ZGELSY or
b/ Alternatively, is there a faster solver than ZGELSS that relies on SVD which could be used instead? I want to stick to SVD as the condition number of A is likely to remain between 10^16 to 10^22 all the time. Some search on the internet indicates that Jacobi SVD (work of Drmac et al) will be faster and more accurate (?) but I could not find the solver based on Jacobi SVD in LAPACK.

Thank's a lot!
Posts: 1
Joined: Fri Dec 23, 2016 5:01 am

Re: Information about complex matrix to be solved

Postby Julien Langou » Mon Jan 02, 2017 10:18 am

a) If you use xGELSS with RCOND=-1, xGELSS essentially sets RCOND to EPS for you and that is that.

So a good starting point to convert a code that uses GELSS with RCOND=-1 to a code that uses GELSY is to set RCOND to 1.1E-016 (in double precision). Unfortunately the meaning of RCOND in GELSS and GELSY is not exactly the same. One RCOND is based on the SVD, the other is based on the R factor of the QR factorization. Slightly different. Yet starting by setting RCOND to 1.1E-016 would be a good start.

(b) xGELSS relies on xBDSQR to compute the singular value decomposition. xBDSQR uses the SVD implicit zero-shift QR algorithm and is pretty slow indeed.

A faster SVD-based least squares solver should be xGELSD. It uses the Divide-and-Conquer algorithm to compute the SVD. You might want to try.

(c) Jacobi SVD is xGESVJ (pure one-sided Jacobi) or xGEJSV (preconditioned Jacobi SVD), but we do not have the linear least squares solver associated with them.

Julien Langou
Posts: 835
Joined: Thu Dec 09, 2004 12:32 pm
Location: Denver, CO, USA

Return to User Discussion

Who is online

Users browsing this forum: No registered users and 4 guests