Skip three different ways to do this calculation (including the opposite) on OS X or iOS.
First, let's do what you more or less tried to do for starters: we will do the calculation using LAPACK and BLAS from the Accelerate view. Note that I used the BLAS function cblas_sgemv instead of vDSP_mmul to multiply the matrix vector. I did this for three reasons. Firstly, it is more idiomatic in code that uses LAPACK. Secondly, LAPACK really wants the matrices to be stored in the main column, which supports BLAS, but vDSP does not. Finally, vDSP_mmul is actually just a wrapper around BLAS matrix multiplications, so we can cut out the middleman. This will work on OS X up to 10.2 and iOS up to 4.0 - i.e. for any purpose you can count on today.
#include <Accelerate/Accelerate.h> #include <stdio.h> #include <string.h> void using_blas_and_lapack(void) { printf("Using BLAS and LAPACK:\n"); // Note: you'll want to store A in *column-major* order to use it with // LAPACK (even though it not strictly necessary for this simple example, // if you try to do something more complex you'll need it). float A[3][3] = {{-4,-3,-5}, {6,-7, 9}, {8,-2,-1}}; float x[3] = { -1, 10, -3}; // Compute b = Ax using cblas_sgemv. float b[3]; cblas_sgemv(CblasColMajor, CblasNoTrans, 3, 3, 1.f, &A[0][0], 3, x, 1, 0.f, b, 1); printf("b := A x = [ %g, %g, %g ]\n", b[0], b[1], b[2]); // You probably don't actually want to compute A^-1; instead you simply // want to solve the equation Ay = b for y (like y = A\b in matlab). To // do this with LAPACK, you typically use the routines sgetrf and sgetrs. // sgetrf will overwrite its input matrix with the LU factorization, so // we'll make a copy first in case you need to use A again. float factored[3][3]; memcpy(factored, A, sizeof factored); // Now that we have our copy, go ahead and factor it using sgetrf. __CLPK_integer n = 3; __CLPK_integer info = 0; __CLPK_integer ipiv[3]; sgetrf_(&n, &n, &factored[0][0], &n, ipiv, &info); if (info != 0) { printf("Something went wrong factoring A\n"); return; } // Finally, use the factored matrix to solve y = A\b via sgetrs. Just // like sgetrf overwrites the matrix with its factored form, sgetrs // overwrites the input vector with the solution, so we'll make a copy // before calling the function. float y[3]; memcpy(y, b, sizeof y); __CLPK_integer nrhs = 1; sgetrs_("No transpose", &n, &nrhs, &factored[0][0], &n, ipiv, y, &n, &info); printf("y := A\\b = [ %g, %g, %g ]\n\n", y[0], y[1], y[2]); }
Then we will do the calculation using LinearAlgebra.h, which is a new feature of the accelerator system in iOS 8.0 and OS X 10.10. It abstracts almost all the accounting information needed for these calculations, but it offers only a tiny piece of full functionality, available in LAPACK and BLAS. Note that while there is some scheme needed to convert our raw arrays to and from la_objects, as soon as we have them, the actual calculation is very simple.
#include <Accelerate/Accelerate.h> #include <stdio.h> void using_la(void) { printf("Using LA:\n"); // LA accepts row-major as well as column-major data, but requires a short // two-step dance to get our data in. float Adata[3][3] = {{-4, 6, 8}, {-3,-7,-2}, {-5, 9,-1}}; float xdata[3] = { -1, 10, -3}; la_object_t A = la_matrix_from_float_buffer(&Adata[0][0], 3, 3, 3, LA_NO_HINT, LA_DEFAULT_ATTRIBUTES); la_object_t x = la_vector_from_float_buffer(xdata, 3, 1, LA_DEFAULT_ATTRIBUTES); // Once our data is stored as LA objects, it easy to do the computation: la_object_t b = la_matrix_product(A, x); la_object_t y = la_solve(A, b); // And finally we need to get our data back out: float bdata[3]; if (la_vector_to_float_buffer(bdata, 1, b) != LA_SUCCESS) { printf("Something went wrong computing b.\n"); return; } else printf("b := A x = [ %g, %g, %g ]\n", bdata[0], bdata[1], bdata[2]); float ydata[3]; if (la_vector_to_float_buffer(ydata, 1, y) != LA_SUCCESS) { printf("Something went wrong computing y.\n"); return; } else printf("y := A\\b = [ %g, %g, %g ]\n\n", ydata[0], ydata[1], ydata[2]); }
And finally, another method. If your matrices are really only 3x3, this is the method I would use, because the overhead associated with any of the previous methods will swamp the actual calculation. However, this only works for matrices 4x4 or smaller. Another new title in iOS 8.0 and OS X 10.10 is specifically focused on small matrix and vector maths, which makes this very simple and effective:
#include <simd/simd.h> #include <stdio.h> void using_simd(void) { printf("Using <simd/simd.h>:\n"); matrix_float3x3 A = matrix_from_rows((vector_float3){-4, 6, 8}, (vector_float3){-3, -7, -2}, (vector_float3){-5, 9, -1}); vector_float3 x = { -1, 10, -3 }; vector_float3 b = matrix_multiply(A, x); printf("b := A x = [ %g, %g, %g ]\n", b[0], b[1], b[2]); vector_float3 y = matrix_multiply(matrix_invert(A),b); printf("y := A^-1 b = [ %g, %g, %g ]\n\n", y[0], y[1], y[2]); }
And finally, let’s just double check that they all give the same results (up to small rounding differences):
scanon$ xcrun -sdk macosx clang matrix.m -framework Accelerate && ./a.out Using BLAS and LAPACK: b := A x = [ 40, -61, 98 ] y := A\b = [ -0.999999, 10, -3 ] Using LA: b := A x = [ 40, -61, 98 ] y := A\b = [ -1, 10, -3 ] Using <simd/simd.h>: b := A x = [ 40, -61, 98 ] y := A^-1 b = [ -1, 10, -3 ]