c++ - 使用 cblas_dgemv() 的矩阵 vector 积的结果不正确

Question

我正在尝试使用 cblas.h 中声明的例程“cblas_dgemv()”来计算矩阵 vector 积。代码如下:

#include <cblas.h>
#include <iostream>

using namespace std;

/* compile with: g++ test.cpp -o test -lblas */

int main (int argc, char** argv)
{
    /* We store the following 3x4 array in column major form:
        A = [1.83292  0.267964 0.382422 0.520162
             0.428562 0.720323 0.839606 1.30816
             0.388731 0.619452 1.01375  0.229333 ];
     */

     const double A[12] = {  1.83292,
                            0.428562,
                            0.388731,
                            0.267964,
                            0.720323,
                            0.619452,
                            0.382422,
                            0.839606,
                            1.01375,
                            0.520162,
                            1.30816,
                            0.229333};

    /* The vector X is:
        X = [ 0.570695, 0.670179, 0.927146, 0.297343 ]';
        where ' is transpose.
     */
    const double X[4] = { 0.570695, 0.670179, 0.927146, 0.297343};
    double Y[4] = {0, 0, 0, 0};

    /* Calculate Y = A*X (alpha = 1 and beta = 0)
     */
    cblas_dgemv(CblasColMajor, CblasNoTrans, 3, 4, 1, A, 3, X, 1, 0, Y, 1);

    /* If I was correct, I should have got:
            Y =[1.69366 1.20999 1.72082 1.38618] = A*X;
        but I get:
            Y = [1.73485 1.89473 1.64507 0] = A'*X;
     */
     cout<<"Y = [";
     for ( unsigned int i = 0; i < 4; i++ )
        cout<<Y[i]<<" ";
    cout<<"]";
}

但是，我始终得到 Y = A'*X，而不是得到 Y = A*X，其中 ' 代表转置。我不确定我是否在某个地方犯了一些经典的愚蠢错误，但经过数小时的尝试，我无法找出问题所在。请帮忙!!

如果需要，我使用的是 Linux 版本 3.2.0-4-amd64 (Debian 4.6.3-14) ) #1 SMP Debian 3.2.57-3+deb7u2 和 g++ (Debian 4.4.7-2) 4.4.7.提前致谢。

Answer 1

你的数学有问题。 3x4 矩阵和 4 分量 vector 的乘积是 3 分量 vector 。在您的例子中，A*X 的乘积是 [1.73485 1.89473 1.64507]。

如果您乘以 A(4x3)的转置，您需要一个 3 分量 vector 与之相乘，乘积是一个 4 分量 vector 。我们称 X' 为 X 的前三个分量 -> X' = [0.570695, 0.670179, 0.927146]。那么 A'*X' = [1.693662 1.209994 1.720827 1.386180]。

#include <stdio.h>
void dgemv(const double *A, const double *u, double *v, const int n, const int m) {
    for(int i=0; i<n; i++) {
        double sum = 0;
        for(int j=0; j<m; j++) {
            sum += A[m*i+j]*u[j];
        }
        v[i] = sum;
    }
}

int main() {                          
    const double A[12] = {1.83292 , 0.267964, 0.382422, 0.520162,
                          0.428562, 0.720323, 0.839606, 1.30816 ,
                          0.388731, 0.619452, 1.01375 , 0.229333};            

    const double AT[12] = {1.83292 , 0.428562, 0.388731,
                           0.267964, 0.720323, 0.619452,
                           0.382422, 0.839606, 1.01375 ,
                           0.520162, 1.30816 , 0.229333};

    const double X[4]  = { 0.570695, 0.670179, 0.927146, 0.297343};
    const double X2[4] = { 0.570695, 0.670179, 0.927146};
    double Y[3], Y2[4];

    dgemv(A, X, Y,  3,4);
    for(int i=0; i<3; i++) printf("%f ", Y[i]); printf("\n");
    dgemv(AT, X2, Y2, 4,3);
    for(int i=0; i<4; i++) printf("%f ", Y2[i]); printf("\n");

}