Swift BLAS cblas_cgemv 复数

Swift BLAS cblas_cgemv complex numbers

我已经能够使用 cblas_sgemv,其中所有值都是真实的。但是,如果不出现 'EXC_BAD_ACCESS' 错误,我就无法使用 cblas_cgemv。对于这个函数,我是否正确地假设复数部分直接出现在函数作为参数的数组中的实数部分之后?例如,如果我有一个矩阵:

1 + 2i, 3 + 4i
5 + 6i, 7 + 8i

则表示为[1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0]。 此外,对于任何给定的矩阵和向量,N、M、lda、incX 和 incY 应该根据其在数组中的表示形式取什么值?

完整示例: 对于 cblas_sgemv,函数原型为:

func cblas_sgemv(_ __Order: CBLAS_ORDER, _ __TransA: CBLAS_TRANSPOSE, _ __M: Int32, _ __N: Int32, _ __alpha: Float, _ __A: UnsafePointer<Float>, _ __lda: Int32, _ __X: UnsafePointer<Float>, _ __incX: Int32, _ __beta: Float, _ __Y: UnsafeMutablePointer<Float>, _ __incY: Int32)

此代码有效:

let matrix: [Float] = [1.0,2.0,3.0,4.0]
let vector: [Float] = [1.0,2.0]
let matrixConverted: UnsafePointer<Float> = UnsafePointer<Float>(matrix)
let vectorConverted: UnsafePointer<Float> = UnsafePointer<Float>(vector)
let matrixSize = sqrt(Double(matrix.count)
let total: [Float] = [Float](count: Int(matrixSize), repeatedValue: 0.0)
let totalConverted: UnsafeMutablePointer<Float> = UnsafeMutablePointer<Float>(total)
cblas_sgemv(CblasRowMajor, CblasNoTrans, Int32(matrixSize), Int32(matrixSize), 1.0, matrixConverted, Int32(matrixSize), vectorConverted, 1, 0, totalConverted, 1)

答案合计存储。那么如果使用上面的复杂矩阵和带有原型的函数 cblas_cgemv 做同样的事情会是什么样子:

func cblas_cgemv(_ __Order: CBLAS_ORDER, _ __TransA: CBLAS_TRANSPOSE, _ __M: Int32, _ __N: Int32, _ __alpha: UnsafePointer<Void>, _ __A: UnsafePointer<Void>, _ __lda: Int32, _ __X: UnsafePointer<Void>, _ __incX: Int32, _ __beta: UnsafePointer<Void>, _ __Y: UnsafeMutablePointer<Void>, _ __incY: Int32) 

可以在此处找到更多信息和参考:https://developer.apple.com/library/prerelease/ios/documentation/Accelerate/Reference/BLAS_Ref/index.html#//apple_ref/doc/uid/TP30000414-SW55

你现有的代码有点太复杂,有错误:

  • 输入矩阵和向量可以直接传递给cblas_sgemv() 函数,不需要 matrixConvertedvectorConverted
  • 创建指向常量数组total可变指针是 不允许。结果向量必须是可变的(并且 totalConverted 也不需要)。

所以你的代码可以简化为:

let matrix: [Float] = [1.0,2.0,3.0,4.0]
let vector: [Float] = [1.0,2.0]

let matrixSize = sqrt(Double(matrix.count))
var total = [Float](count: Int(matrixSize), repeatedValue: 0.0)

cblas_sgemv(CblasRowMajor, CblasNoTrans, Int32(matrixSize), Int32(matrixSize), 1.0, matrix, Int32(matrixSize), vector, 1, 0, &total, 1)

记录了 BLAS 例程中 复杂 数字的布局 在 <cblas.h>:

 * A note on complex data layouts:
 *
 * In order to allow straightforward interoperation with other libraries and
 * complex types in C and C++, complex data in BLAS is passed through an opaque
 * pointer (void *).  The layout requirements on this complex data are that
 * the real and imaginary parts are stored consecutively in memory, and have
 * the alignment of the corresponding real type (float or double).  The BLAS
 * complex interfaces are compatible with the following types:
 *
 *     - The C complex types, defined in <complex.h>.
 *     - The C++ std::complex types, defined in <complex>.
 *     - The LAPACK complex types, defined in <Accelerate/vecLib/clapack.h>.
 *     - The vDSP types DSPComplex and DSPDoubleComplex, defined in <Accelerate/vecLib/vDSP.h>.
 *     - An array of size two of the corresponding real type.
 *     - A structure containing two elements, each of the corresponding real type.

因此,相乘

| 1 + 2i  3 + 4i |    | 1 + 2i |
|                |  * |        |
| 5 + 6i  7 + 8i |    | 3 + 4i |

你可以将每个复数表示为两个浮点数 连续存储的点数:

let matrix: [Float] = [1.0,2.0, 3.0,4.0, 5.0,6.0, 7.0,8.0]
let vector: [Float] = [1.0,2.0, 3.0,4.0]

let matrixSize = sqrt(Double(matrix.count/2))
var total = [Float](count: vector.count, repeatedValue: 0.0)

let alpha : [Float] = [1.0, 0.0]
let beta : [Float] = [1.0, 0.0]

cblas_cgemv(CblasRowMajor, CblasNoTrans, Int32(matrixSize), Int32(matrixSize), beta, matrix, Int32(matrixSize), vector, 1, alpha, &total, 1)

或者您可以使用 DSPComplex 表示复数, COMPLEX__CLPK_complex 结构(它们都有 相同的布局):

let matrix = [DSPComplex(real: 1.0, imag: 2.0), DSPComplex(real: 3.0, imag: 4.0),
              DSPComplex(real: 5.0, imag: 6.0), DSPComplex(real: 7.0, imag: 8.0)]

let vector = [DSPComplex(real: 1.0, imag: 2.0), DSPComplex(real: 3.0, imag: 4.0)]

let matrixSize = sqrt(Double(matrix.count))
var total  =  [DSPComplex](count: Int(matrixSize), repeatedValue: DSPComplex())

var alpha = [DSPComplex(real: 1.0, imag: 0.0)]
var beta  = [DSPComplex(real: 1.0, imag: 0.0)]

cblas_cgemv(CblasRowMajor, CblasNoTrans, Int32(matrixSize), Int32(matrixSize), alpha, matrix, Int32(matrixSize), vector, 1, beta, &total, 1)

在任何一种情况下,维度 MN 等都是指复数的计数,因此它们与示例中的值 M=N=2 相同 的实数,并且 alphabeta 也是一个数组 代表一个复杂的因素。