added new routines [S,D,C,Z]GELST and test code for the least squares

[scr2016]

Description
added new routines [S,D,C,Z]GELST and test code for the least squares
squares algorithm that solves overdetermined or underdetermined systems for GE matrices using QR or LQ factorization with compact WY representation of Q inside the routine.

*GELST is similar to *GELS, but unlike *GELS does not compute triangular blocks twice, i.e. runs faster.

modified: SRC/Makefile
new file: SRC/cgelst.f
new file: SRC/dgelst.f
new file: SRC/sgelst.f
new file: SRC/zgelst.f
modified: TESTING/LIN/alahd.f
modified: TESTING/LIN/cdrvls.f
modified: TESTING/LIN/cerrls.f
modified: TESTING/LIN/ddrvls.f
modified: TESTING/LIN/derrls.f
modified: TESTING/LIN/sdrvls.f
modified: TESTING/LIN/serrls.f
modified: TESTING/LIN/zdrvls.f

[angsch]

@scr2016 Is there a particular reason why your changes are a separate routine and not an upgrade of GELS? Your changes appear fully compatible with GELS and could be a great modernization that could be widely picked up by, say, octave or Julia, which likely call GELS as the de facto standard LS solver.

Edit: I mixed up geqrt and geqrt3. Sorry. Have you considered interleaving the QR factorization and the transform of B, which reduces the workspace to a single triangular factor? Sketched for QR, $m \geq n$:

do i = 1, nb
  ib = min(nb,n-i+1) ! No special treatment of the last block
  ! Compute QR factorization of panel A(i:end, i:i+ib-1)
  call geqrt()
  ! Trailing matrix update H**T * A(i:m,i+ib:n)
  call larfb()
  ! Transform H**T * B(i:m,1:nrhs)
  call larfb()

Then your change is really fully compatible with GELS

[dbielich]

I see edits to all TESTING/LIN/Xerrls.f files except TESTING/LIN/zerrls.f.

So there is no way to check any errors from the double complex code currently.