LAPACK sgelss usage

Question

I would like to use the sgelss function from the LAPACK framework to solve an overdetermined system of linear equations.

minimize || A*x-b||

For this I chose to use the sgelss function from LAPACK, the header file can be found here: SGELSS

My problem is that I seam to use it wrong, since the results don't fit.

An example case, I use these matricies:

double a[] =
    {
       25,5,1,
       16,4,1,
        9,3,1,
        4,2,1,
        1,1,1,
    };
double b[] =
     {
       31,21,13,7,3
     }

And then I would call:

sgelss_(&(5), &(3), &(1), a, (5), b, &(5), 3, -1, &(3), malloc(3000), &(3000), &output);

[[ Since sgelss expects most input variables not to be normal variables but pointers I use &(...) here in order to skip instantiating variables and then referencing their memory location for sake of readability here! ]]

I would expect that afterwards b would be [1,1,1,XXXXX] since my inputs a and b were set up this way. Unfortunately this is not the case.

I have also tried rotating a (switching rows and columns), without success.

Answer 1

There are a lot of problems in the code that you posted; I assume that it’s abstracted from what you’re really doing, as it clearly could never even compile as is. There are two major issues, one or both of which I suspect to be the cause of the trouble in your actual code:

First, you’re attempting to use sgelss_ with double data. sgelss_ operates on float (single-precision). If you have double-precision data, you need to use dgelss_ instead.

Second, your parameters describe a 5x3 matrix, but remember that LAPACK uses column-major ordering for matrix elements. This means that the matrix as described by your code is:

25  1  2 
 5  9  1
 1  3  1
16  1  1
 4  4  1

Somehow I doubt that’s the matrix you really want. More likely you want the matrix:

25  5  1
16  4  1
 9  3  1
 4  2  1
 1  1  1

I went ahead and put together a working version that assumes this is the matrix you’re actually trying to use:

#include <Accelerate/Accelerate.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[]) {
    double a[] = /* column major storage order! */
    {
        25, 16,  9,  4,  1,
         5,  4,  3,  2,  1,
         1,  1,  1,  1,  1,
    };
    double b[] =
    {
        31,21,13,7,3
    };

    // Setup parameters
    __CLPK_integer m = 5;
    __CLPK_integer n = 3;
    __CLPK_integer nrhs = 1;
    __CLPK_integer lda = 5;
    __CLPK_integer ldb = 5;
    double *s = malloc(3 * sizeof*s);
    double rcond = -1.0f; // use machine precision
    __CLPK_integer rank;
    __CLPK_integer info;
    // Query correct worksize
    double worksize;
    __CLPK_integer lwork = -1;
    dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, &worksize, &lwork, &info);
    // Allocate workspace
    lwork = worksize;
    double *work = malloc(lwork * sizeof *work);
    // Do computation
    dgelss_(&m, &n, &nrhs, a, &lda, b, &ldb, s, &rcond, &rank, work, &lwork, &info);
    // Free workspace
    free(work);
    // Print result vector
    for (int i=0; i<3; ++i)
        printf("%g\t", b[i]);
    printf("\n");
    return 0;
}