让我们走过三种不同的方式来计算这个问题（包括反向）在OS X或iOS上。首先，让我们做你开始时所尝试的：使用加速框架中的LAPACK和BLAS进行计算。请注意，我使用了BLAS函数 cblas_sgemv 而不是 vDSP_mmul 来执行矩阵向量乘法。我这样做有三个原因。首先，在使用LAPACK的代码中，这更符合惯例。其次，LAPACK真正希望矩阵以
列主顺序存储，而BLAS支持，而vDSP则不支持。最后， vDSP_mmul 实际上只是BLAS矩阵乘法例程的包装器，因此我们可以切掉中间人。这将在OS X回到10.2和iOS回到4.0上工作-即任何您今天可以合理预期遇到的目标。
#include <Accelerate/Accelerate.h>
#include <stdio.h>
#include <string.h>
void using_blas_and_lapack(void) {
printf("Using BLAS and LAPACK:\n");
float A[3][3] = {{-4,-3,-5}, {6,-7, 9}, {8,-2,-1}};
float x[3] = { -1, 10, -3};
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]);
float factored[3][3];
memcpy(factored, A, sizeof factored);
__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; }
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]);
}
接下来，让我们使用 LinearAlgebra.h 进行计算，这是 iOS 8.0 和 OS X 10.10 中加速框架的一个新功能。它在几乎所有这些计算所需的簿记方面进行了抽象处理，但它仅提供了 LAPACK 和 BLAS 中可用完整功能的一小部分。请注意，虽然需要一些样板文件将我们的原始数组转换为和从 la_objects 转换回来，但一旦拥有了它们，实际计算是非常简单的。
#include <Accelerate/Accelerate.h>
#include <stdio.h>
void using_la(void) {
printf("Using LA:\n");
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);
la_object_t b = la_matrix_product(A, x);
la_object_t y = la_solve(A, b);
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]);
}
最后，还有一种方法。如果你的矩阵只是3x3的话，我会使用这种方法，因为之前的两种方法所涉及到的开销将会淹没实际计算。然而，这种方法仅适用于大小为4x4或更小的矩阵。在iOS 8.0和OS X 10.10中有一个专门针对小矩阵和向量数学的新头文件，使这个过程变得非常简单和高效：
#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]);
}
最后，让我们仔细检查一下这些结果是否完全相同（除了小的四舍五入差异）：
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 ]
