MKLスパース行列ベクトル乗算を使用する

Question

MKLスパース行列ベクトル乗算ルーチンを使用する単純なC++コード例はありますか？複雑なベクトルで複素対称行列（下三角で格納）を乗算するには "mkl_zcsrsymv"を使用する必要がありますが、これに関する単一のデモの例は見つかりませんでした。複雑なケースのために自分のコードをコンパイルできません。

Answer 1

mkl_zcsrsymvのドキュメントをIntel's homepageに読みます。ここで対称は文字通り取られます！ （エルミート語を意味するものではありません）

最大限の便宜のためにicpc -mkl test.cとコンパイルしてください。

#include <stdio.h> 
#include "mkl_spblas.h" 

int main() 
{ 
    /* Matrix in CRS format 
    * 
    * { { 0, 0, i } 
    * { 0, -1, 2 } 
    * { i, 2, 0 } } 
    */ 
    int m = 3; 
    MKL_Complex16 a[] = { {0,1}, {-1,0}, {2,0}, {0,1}, {2,0} }; 
    int ia[] = { 0, 1, 3, 5 }; 
    int ja[] = { 2, 1, 2, 0, 1 }; 

    MKL_Complex16 x[] = { {1,0}, {2,0}, {3,0} }; 
    MKL_Complex16 y[] = { {0,0}, {0,0}, {0,0} }; 

    char uplo = 'L'; 
    // Use MKL to compute 
    // y = A*x 
    mkl_cspblas_zcsrsymv(&uplo, &m, a, ia, ja, x, y); 

    printf("y = { (%g,%g), (%g,%g), (%g,%g) }\n", 
     y[0].real, y[0].imag, 
     y[1].real, y[1].imag, 
     y[2].real, y[2].imag 
    ); 
} 

出力はy = { (0,3), (4,0), (4,1) }です。 WolframAlphaで確認してください。

---

mkl_dcsrmvの例でもあります。

#include <stdio.h> 
#include "mkl_spblas.h" 

int main() 
{ 
    /* Matrix in CRS format 
    * 
    * { { 0, 0, 1 } 
    * { 0, -1, 2 } 
    * { 1, 0, 0 } } 
    */ 
    int m = 3; 
    int k = 3; 
    double val[] = { 1, -1, 2, 1 }; 
    int indx[] = { 2, 1, 2, 0 }; 
    int pntrb[] = { 0, 1, 3 }; 
    int pntre[] = { 1, 3, 4 }; 

    double x[] = { 1, 2, 3 }; 
    double y[] = { 0, 0, 0 }; 

    double alpha = 1; 
    double beta = 0; 
    char transa = 'N'; 
    char matdescra[] = { 
    'G', // type of matrix 
    ' ', // triangular indicator (ignored in multiplication) 
    ' ', // diagonal indicator (ignored in multiplication) 
    'C' // type of indexing 
    }; 

    // Use MKL to compute 
    // y = alpha*A*x + beta*y 
    mkl_dcsrmv(&transa, &m, &k, &alpha, matdescra, val, indx, pntrb, pntre, x, &beta, y); 

    printf("y = { %g, %g, %g }\n", y[0], y[1], y[2]); 
} 

出力がy = { 3, 4, 1 }です。 WolframAlphaで確認してください。

---

これを再生しているうちに、これはArmadilloと直接互換性があることがわかりました。これにより、C++での使用が非常に便利になります。ここではまず、Armadilloでランダム対称行列を生成し、それを疎に変換します。これにはランダムなベクトルが乗算されます。最後に、結果をArmadilloの疎なマトリックスベクター産物と比較します。精度はかなり大幅に異なります。

#include <iostream> 
#include <armadillo> 
#define MKL_Complex16 arma::cx_double 
#include "mkl_spblas.h" 

int main() 
{ 
    /* Matrix in CRS format 
    * 
    * { { 0, 0, i } 
    * { 0, -1, 2 } 
    * { i, 2, 0 } } 
    */ 
    int dim = 1000; 
    arma::sp_cx_mat a(arma::randu<arma::cx_mat>(dim,dim)); 
    a += a.st(); 
    arma::cx_vec x = arma::randu<arma::cx_vec>(dim); 
    arma::cx_vec y(dim); 

    char uplo = 'L'; 
    // Use MKL to compute 
    // y = A*x 
    mkl_cspblas_zcsrsymv(&uplo, &dim, 
         a.values, (int*)a.col_ptrs, (int*)a.row_indices, 
         x.memptr(), y.memptr()); 

    std::cout << std::boolalpha 
      << arma::approx_equal(y, a*x, "absdiff", 1e-10) 
      << '\n'; 
}