Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Using Accelerate (CLAPACK) to solve an augmented matrix?

Does anyone know what function/method to use in Accelerate (CLAPACK) to solve an augmented matrix such as the one below? Looking for any sample code, links to samples, hints on how to solve the matrix. I've been looking through the documentation but most everything has to do with more complex graphical systems and there are hundreds of seemingly similar methods.

 ____            ____
|                    |
|  4   3  -1   |  2  |
| -2   3   8   |  0  |
|  0   2   6   | -1  |
|____            ____|
like image 967
John Avatar asked Dec 29 '22 04:12

John


1 Answers

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

int main(int argc, char *argv[]) {

    /* Dimension of the matrix */
    __CLPK_integer n = 3;

    /* Number of right-hand side vectors to solve for */
    __CLPK_integer nrhs = 1;

    /* Note the ordering of the entries in A.  LAPACK uses "column major"
       ordering as follows:

        0 3 6
        1 4 7
        2 5 8

       This is a Fortran-ism that persists in CLAPACK.                    */
    double A[9] = {4.0, -2.0, 0.0, 3.0, 3.0, 2.0, -1.0, 8.0, 6.0 };

    /* "Leading dimension" of the matrix; most of the time, this is just the
       matrix dimension (but not always; you'll learn about this by the time
       you need to use it. */
    __CLPK_integer lda = 3;

    /* Integer array to hold information about the matrix factorization */
    __CLPK_integer ipiv[3];

    /* Right hand side to solve for */
    double x[3] = { 2.0, 0.0, -1.0 };

    /* Leading dimension of the right hand side vector */
    __CLPK_integer ldb = 3;

    /* Status variable */
    __CLPK_integer info;

    /* Solve the augmented system with a call to dgesv_.  Note that this
       routine will overwrite the contents of the array A with a factored
       form of the matrix.  If you need the original matrix, you need to
       copy it before calling dgesv_.  Note that all the scalar arguments
       are passed as pointers; this too is a Fortran-ism.                 */
    dgesv_(&n, &nrhs, A, &lda, ipiv, x, &ldb, &info);

    /* Handle error conditions */
    if (info)
        printf("Could not solve system; dgesv exited with error %d\n", info);

    /* If no error, print the result */
    else
        printf("Solution is [%f, %f, %f]\n", x[0], x[1], x[2]);

    return 0;
}

Compile and run:

 scanon$ gcc solve.c -framework Accelerate
 scanon$ ./a.out 
 Solution is [-0.479167, 1.125000, -0.541667]
like image 76
Stephen Canon Avatar answered Jan 05 '23 01:01

Stephen Canon