1
$\begingroup$

I am trying to use Least Squares Minimization to solve a the matrix problem: b = A*x for x. The system is overdetermined, and A is a dense matrix.

In the LAPACK library, I believe the routine DGELSY to be the best suited for this problem (or whatever is closest to Matlab's LSQMINNORM function). However, I am a relative amateur at coding, and Fortran in particular, and am having issues with the inputs for this function.

Thy syntax is as follows:

call dgelsy(m, n, nrhs, a, lda, b, ldb, jpvt, rcond, rank, work, lwork, info)

(Here is the site for the documentation: https://software.intel.com/en-us/mkl-developer-reference-fortran-gelsy).

For example, let A be a 40,000 x 3,000 matrix, b is a 40,000 element vector, and x is a 3,000 element vector.

Question: What should the inputs for DGELSY should be in this case?

[Aside: If there is a LAPACK function that is better suited to solving this problem, feel free to let me know!]

Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ The documentation looks pretty clear to me at first glance. What specifically don't you understand? $\endgroup$
    – LedHead
    Jun 7, 2019 at 3:28
  • $\begingroup$ The inputs that I don't understand are: work (and correspondingly, lwork), rcond, and rank. My understanding is that 'a' corresponds to the matrix A, 'b' corresponds to the matrix b, and 'jpvt' corresponds to the output matrix x. 'lda' and 'ldb', for the case of a typically defined matrix (starting with index 1), are just the number of rows in 'a' and 'b'. Please correct me if I am wrong, though. $\endgroup$
    – jecht300
    Jun 7, 2019 at 3:36
  • $\begingroup$ @jecht300 work is a workspace array (i.e. a block of allocated memory) you pass to the routine for it to do calculation in (see your own link for size hints) rcond is a tolerance (i.e. small number) on how far to push the QR factorisation (necessary since computer arithmetic isn't exact). Rank is an output integer telling you how many rows your specified rcond actually meant the routine generated. $\endgroup$
    – origimbo
    Jun 7, 2019 at 12:31

1 Answer 1

Reset to default
1
$\begingroup$

You're most of the way there.

  • On input, a is the matrix A and b is the right-hand side vector or matrix b. When the subroutine finishes, a has been completely changed (it now holds the QR factorization) and b now contains the result x.
  • jpvt does not contain the output matrix x, as you thought, it contains pivoting information which you probably don't need.
  • I'd just set rcond to some small number like 1.e-8
  • As @origimbo commented, work is a workspace array needed by the DGELSY subroutine. LAPACK is supposed to be really fast, and memory allocation is really slow, but this subroutine needs extra space for its calculations. Instead of allocating space itself it requires you to allocate the space. If you were going to call this subroutine lots of times it makes sense that you allocated memory just once, instead of doing it inside DGELSY every single time.
  • Setting lwork (and correspondingly the size of work) can be confusing. You can call DGELSY first in "workspace query" mode just to get the suggested value for lwork, or set it to a value that's going to be generous enough. For your case lwork=30000 should be plenty. Look at the documentation for more info.
  • lda and ldb are usually the same as a and b. They're there in case you're operating on a subblock of a larger matrix. Since the data must be in column-major format, this tells LAPACK how many values it must "skip over" between reading/writing columns of data.
Improve this answer
$\endgroup$
0

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.