Least squares, QR and pivoting

[italo1337]

I am trying to solve an overdetermined system in least squares sense with QR decomposition with pivoting.

[cpp]#include 
#include "mkl_lapacke.h"
#include "mkl_cblas.h"
#include "mkl_types.h"
#include "least_squares.h"

void solve_least_squares(double *A, double *b, double *result, int m, int n) {

  int i;
  int *jpvt = calloc(n, sizeof(int));
  double *tau = calloc(n, sizeof(double));

  LAPACKE_dgeqp3(LAPACK_ROW_MAJOR, m, n, A, n, jpvt, tau);
  LAPACKE_dormqr(LAPACK_ROW_MAJOR, 'L', 'T', m, 1, n, A, n, tau, b, 1);
  cblas_dtrsm(CblasRowMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, n, 1, 1.0, A, n, b, 1);

  for ( i = 0 ; i < n ; i++ ) {
    result = b[jpvt-1];
  }

  free(jpvt);
  free(tau);

}
[/cpp]

The result is correct, but sometimes in wrong order. Am I doing anything wrong when mapping the result?

It is related to pivoting. When I add to the beggining of the function:

[cpp]  for ( i = 0 ; i < n ; i++ ) {
    jpvt = i+1;
  }
[/cpp]

This prevents pivoting and yields the correct result in the correct order.

[mecej4]

See this example:

f08bffe.f

In particular, see how the permutation is undone for the solution in the DO 100 loop. More generally, undoing a sequence of permutations requires that the individual permutations be undone in reverse order.