Hi,

So you have a 5-by-6 matrix A. The matrix A has full rank 5.

You have a right-hand side b of length 5. You wish to solve Ax =b.

In your case, we have

- Code: Select all
`a =[ -1/2 1/2 3/2 5/2 7/2 9/2 ]; A = [ a.^0; a.^1; a.^2; a.^3; a.^4 ];`

b = [ 0; 1; 0; 0; 0];

In this case, the set of solution is of the form ( x1 + Null( A ) ). So you are looking for

(1) one particular solution. In your case, you chose x1 = [-11/12; 17/24; 3/8; -5/24; 1/24; 0]. Fine.

(2) and a nonzero vector in the null space. In your case you chose u = [-1; 5; -10; 10; -5; 1];

so that all solution writes x1 + k * u where k is any scalar. All good.

So now how LAPACK can help you solve your problem.

First to find a particular solution you can use LAPACK DGELS. Now LAPACK DGELS will give the particular solution with minimum norm. In general it has no zero in it. In your case, the answer returned by LAPACK is:

- Code: Select all
`x2 = `

-0.922949735449755

0.739748677248692

0.312169312169325

-0.14550264550266

0.0102513227513299

0.0062830687830675

Note that Ax2 = b (as Ax1 =b) but || x2 || < || x1 ||. This is a consequence of the fact that the x2 returned by LAPACK not only is a solution, it is the solution with minimum 2-norm.

A point of reference to understand the lapack answer as opposed to say the matlab A\b answer (the matlab's answer will have some zeros in it) is my post on the forum:

viewtopic.php?f=5&t=160Second to get the null space, several ways would be OK. You can do an SVD of A but that would be expensive. The best way is to do a FULL QR factorization of A^T and get the Null space from the last column of Q. (In your case, we know that the null space is of dimension 1, so we know we are interested in the last column of Q.) Then Null(A) = Span( w ).

So now you can both (1) and (2) at once. By "at once", I mean with only one QR factorization of A.

This is pretty cool. The problem is that you do not have an LAPACK code to do this. I can provide a Matlab pseudocode below:

- Code: Select all
`[Q,R]=qr(A');`

w=Q(:,6);

x2 = Q(1:6,1:5)*((R(1:5,1:5)')\b);

Then solutions are of the form x2 + k*w where k is any scalar. x2 is a particular solution of Ax =b, to be precise, x2 is the solution of Ax=b with minimum 2-norm. w is a nonzero vector in Null(A).

Julien.