任何人都可以向我展示openblas中函数dspr的示例

Question

我正在尝试将Openblas中的函数dspr与Rcpp一起使用。DSPR 的目的是：

A := alpha*x*x**T + A

就我而言，首先我将A定义为一个矩阵，所有元素都0, alpha=1, x=(1,3)，所以，最终的矩阵 A 应该是{(1,3),(3,9)}的，但我从来没有得到正确的结果，我设置参数如下：

cblas_dspr(CblasColMajor,CblasUpper,2, 1, &(x[0]),1, &(A[0]));

谁能告诉我如何设置正确的dspr参数？谢谢。

Answer 1

我的机器上/usr/include/cblas,h的文件显示了 BLAS 的 C 接口的以下签名：

void cblas_dspr(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,
                const int N, const double alpha, const double *X,
                const int incX, double *Ap);

试试看。 您可以通过x.begin()或&(x[0])获得 Rcpp 载体的开头，就像您所做的那样。

不过，这里与 Rcpp 无关。

Answer 2

从您自己的帖子中重复：BLAS 二元产品执行A := alpha*x*x' + A

因此，A 需要用零值初始化。此外，不要忘记A是上三角矩阵。

为了进一步阅读，我推荐以下链接：

https://developer.apple.com/library/ios/documentation/Accelerate/Reference/BLAS_Ref/Reference/reference.html

https://github.com/foundintranslation/Kaldi/blob/master/tools/ATLAS/interfaces/blas/C/src/cblas_dspr.c

但是，您想要一个例子。这里是：

/** dspr_demo.cpp */
#include <cblas.h>
#include <iostream>
#include <cstdlib>
using namespace std;
int main(int argc, char** argv)
{
  int n=2;
  double* x = (double*)malloc(n*sizeof(double));
  double* upperTriangleResult = (double*)malloc(n*(n+1)*sizeof(double)/2);
  for (int j=0;j<n*(n+1)/2;j++) upperTriangleResult[j] = 0;
  x[0] = 1; x[1] = 3;
  cblas_dspr(CblasRowMajor,CblasUpper,n,1,x,1,upperTriangleResult);
  double*& A = upperTriangleResult;
  cout << A[0] << "t" << A[1] << endl << "*t" << A[2] << endl;
  free(upperTriangleResult); free(x);
  return EXIT_SUCCESS;
}