2012-11-05 19 views
4

我需要在C++中實現一個矩陣類,其中一個操作必須通過dgemm進行矩陣乘法。我的教授在C語言課上做過一個例子,但由於某種原因,我無法使用C++來實現它。這是我的頭文件matrix.h:在C++中訪問dgemm

#include <iostream> 
#include <stdlib.h> 

extern "C" { 
#include "blas.h" 
} 

[blah blah blah, matrix class here; overloaded * operator will do the matrix  
multiplication] 

matrix operator* (const matrix &other){ 

matrix AxB(Nrows, other.Ncolumns, "(" + name + "*" + "other.name" + ")"); 

char TRANSA = 'N'; 
char TRANSB = 'N'; 
int M = Nrows; 
int N = other.Ncolumns; 
int K = Ncolumns; 
double alpha = 1.; 
double beta = 0.; 

dgemm_ (&TRANSA, 
     &TRANSB, 
     &M, 
     &N, 
     &K, 
     &alpha, 
     data, 
     &M, 
     other.data, 
     &K, 
     &beta, 
     AxB.data, 
     &M); 

return AxB; 

[blah blah blah, overloaded = operator here; i'm positive this is not the problem 
since it works for matrix addition] 

主要功能:

#include "matrix.h" 

main(){ 

// entries of A and B are randomized between 0 and 1 
matrix A(5,5); 
matrix B(5,5); 

matrix C = A*B; 

} 

現在blas.h頭文件:

我使用這裏列出的子程序:http://www.netlib.org/blas/dgemm.f

基本上,我們在C類實現類中構建矩陣的方式就像一個long a使用(double *)calloc(行*列,sizeof(double))進行rray。

我的C++實行的作用:

double **data; 

data = new double[rows]; 

for(int i=1; i<=rows; ++i){ 
    data[i-1] = new double[columns]; 
} 

,然後我可以使用索引數據[i] [j]。但是由於dgemm子程序應該採用double * A,double * B等,但是我的矩陣是double ** A等等,我該如何解決這個問題?

這是錯誤信息,我從Valgrind的獲得:

==18845== Invalid write of size 8 
==18845== at 0x4165C1C: ATL_dgezero (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4149DA7: ATL_dNCmmJIK (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x415BE8E: ATL_dgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4182F54: ATL_dptgemm_nt (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4183140: ATL_dptgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x40818F6: atl_f77wrap_dgemm_ (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4E4E0004: ??? 
==18845== Address 0x478d680 is 8 bytes after a block of size 16 alloc'd 
==18845== at 0x402B454: operator new[](unsigned int) (in /usr/lib/valgrind 
/vgpreload_memcheck-x86-linux.so) 
==18845== by 0x8049045: dmatrix::dmatrix(int, int, std::string) (dmatrix.hpp:77) 
==18845== by 0x8049532: dmatrix::operator*(dmatrix const&) (dmatrix.hpp:218) 
==18845== by 0x8048DCB: main (main.cpp:17) 
==18845== 
==18845== Invalid write of size 8 
==18845== at 0x4165C1F: ATL_dgezero (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4149DA7: ATL_dNCmmJIK (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x415BE8E: ATL_dgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4182F54: ATL_dptgemm_nt (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4183140: ATL_dptgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x40818F6: atl_f77wrap_dgemm_ (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4E4E0004: ??? 
==18845== Address 0x478d678 is 0 bytes after a block of size 16 alloc'd 
==18845== at 0x402B454: operator new[](unsigned int) (in /usr/lib/valgrind 
/vgpreload_memcheck-x86-linux.so) 
==18845== by 0x8049045: dmatrix::dmatrix(int, int, std::string) (dmatrix.hpp:77) 
==18845== by 0x8049532: dmatrix::operator*(dmatrix const&) (dmatrix.hpp:218) 
==18845== by 0x8048DCB: main (main.cpp:17) 
==18845== 
==18845== Invalid write of size 8 
==18845== at 0x4165C22: ATL_dgezero (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4149DA7: ATL_dNCmmJIK (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x415BE8E: ATL_dgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4182F54: ATL_dptgemm_nt (in /usr/lib/atlas-base/atlas  
/libblas.so.3gf.0) 
==18845== by 0x4183140: ATL_dptgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x40818F6: atl_f77wrap_dgemm_ (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4E4E0004: ??? 
==18845== Address 0x478d690 is not stack'd, malloc'd or (recently) free'd 
==18845== 
==18845== Invalid write of size 8 
==18845== at 0x4165C25: ATL_dgezero (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4149DA7: ATL_dNCmmJIK (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x415BE8E: ATL_dgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4182F54: ATL_dptgemm_nt (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4183140: ATL_dptgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x40818F6: atl_f77wrap_dgemm_ (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4E4E0004: ??? 
==18845== Address 0x478d688 is not stack'd, malloc'd or (recently) free'd 
==18845== 
==18845== Invalid read of size 8 
==18845== at 0x4143D60: ATL_dJIK0x0x0NN0x0x0_aX_bX (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4149A9D: ATL_dNCmmJIK (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x415BE8E: ATL_dgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4182F54: ATL_dptgemm_nt (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4183140: ATL_dptgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x40818F6: atl_f77wrap_dgemm_ (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4E4E0004: ??? 
==18845== Address 0x478cf80 is not stack'd, malloc'd or (recently) free'd 
==18845== 
==18845== Invalid read of size 8 
==18845== at 0x4143D64: ATL_dJIK0x0x0NN0x0x0_aX_bX (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4149A9D: ATL_dNCmmJIK (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x415BE8E: ATL_dgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4182F54: ATL_dptgemm_nt (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4183140: ATL_dptgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x40818F6: atl_f77wrap_dgemm_ (in /usr/lib/atlas-base/atlas  
/libblas.so.3gf.0) 
==18845== by 0x4E4E0004: ??? 
==18845== Address 0x478d130 is 0 bytes after a block of size 16 alloc'd 
==18845== at 0x402B454: operator new[](unsigned int) (in /usr/lib/valgrind 
/vgpreload_memcheck-x86-linux.so) 
==18845== by 0x8049045: dmatrix::dmatrix(int, int, std::string) (dmatrix.hpp:77) 
==18845== by 0x8048D49: main (main.cpp:6) 
==18845== 
==18845== Invalid read of size 8 
==18845== at 0x4143D4C: ATL_dJIK0x0x0NN0x0x0_aX_bX (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4149A9D: ATL_dNCmmJIK (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x415BE8E: ATL_dgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4182F54: ATL_dptgemm_nt (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4183140: ATL_dptgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x40818F6: atl_f77wrap_dgemm_ (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4E4E0004: ??? 
==18845== Address 0x478d678 is 0 bytes after a block of size 16 alloc'd 
==18845== at 0x402B454: operator new[](unsigned int) (in /usr/lib/valgrind 
/vgpreload_memcheck-x86-linux.so) 
==18845== by 0x8049045: dmatrix::dmatrix(int, int, std::string) (dmatrix.hpp:77) 
==18845== by 0x8049532: dmatrix::operator*(dmatrix const&) (dmatrix.hpp:218) 
==18845== by 0x8048DCB: main (main.cpp:17) 
==18845== 
==18845== Invalid write of size 8 
==18845== at 0x4143D93: ATL_dJIK0x0x0NN0x0x0_aX_bX (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4149A9D: ATL_dNCmmJIK (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x415BE8E: ATL_dgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x4182F54: ATL_dptgemm_nt (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4183140: ATL_dptgemm (in /usr/lib/atlas-base/atlas/libblas.so.3gf.0) 
==18845== by 0x40818F6: atl_f77wrap_dgemm_ (in /usr/lib/atlas-base/atlas 
/libblas.so.3gf.0) 
==18845== by 0x4E4E0004: ??? 
==18845== Address 0x478d678 is 0 bytes after a block of size 16 alloc'd 
==18845== at 0x402B454: operator new[](unsigned int) (in /usr/lib/valgrind 
/vgpreload_memcheck-x86-linux.so) 
==18845== by 0x8049045: dmatrix::dmatrix(int, int, std::string) (dmatrix.hpp:77) 
==18845== by 0x8049532: dmatrix::operator*(dmatrix const&) (dmatrix.hpp:218) 
==18845== by 0x8048DCB: main (main.cpp:17) 
==18845== 
==18845== Invalid read of size 8 
==18845== at 0x8048DEA: main (main.cpp:19) 
==18845== Address 0x0 is not stack'd, malloc'd or (recently) free'd 
==18845== 
==18845== 
==18845== Process terminating with default action of signal 11 (SIGSEGV) 
==18845== Access not within mapped region at address 0x0 
==18845== at 0x8048DEA: main (main.cpp:19) 
==18845== If you believe this happened as a result of a stack 
==18845== overflow in your program's main thread (unlikely but 
==18845== possible), you can try to increase the size of the 
==18845== main thread stack using the --main-stacksize= flag. 
==18845== The main thread stack size used in this run was 8388608. 
==18845== 
==18845== HEAP SUMMARY: 
==18845==  in use at exit: 3,581 bytes in 39 blocks 
==18845== total heap usage: 49 allocs, 10 frees, 7,760 bytes allocated 
==18845== 
==18845== LEAK SUMMARY: 
==18845== definitely lost: 128 bytes in 4 blocks 
==18845== indirectly lost: 0 bytes in 0 blocks 
==18845==  possibly lost: 88 bytes in 4 blocks 
==18845== still reachable: 3,365 bytes in 31 blocks 
==18845==   suppressed: 0 bytes in 0 blocks 
==18845== Rerun with --leak-check=full to see details of leaked memory 
==18845== 
==18845== For counts of detected and suppressed errors, rerun with: -v 
==18845== ERROR SUMMARY: 111 errors from 9 contexts (suppressed: 0 from 0) 
Segmentation fault (core dumped) 

,直到我試圖實現DGEMM矩陣乘法,我沒有得到任何錯誤,無泄漏,所以我敢肯定,我所有的煩惱說謊使用dgemm實現

回答

1

C版本的內存佈局與C++版本中的內存佈局不同。由於BLAS預計C版本使用的佈局類型,您的C++版本將不起作用。

因此,您需要在您的C++版本中分配一個大的一維數組;你可以重載operator()以獲得二維數組的索引。或者對於生產代碼,請使用庫,如Eigen

+0

嗯,這就是我想的,謝謝。我想我會回去改變我的索引。 – user1799323

+0

要做到這一點,我應該這樣做: – user1799323

+0

應該做什麼? – janneb

0

對不起,我沒有意識到在響應會發送響應時按下Enter鍵。

我改變了索引:

double *data; 

data = new double[rows*columns]; 

那麼每次我想指數,我只是將數據[(I-1)*列數+(J-1)]爲行主索引。該代碼現在可以正常工作,並編譯0錯誤和0內存泄漏!唯一的情況是我現在得到完全亂碼的結果。

這是否與這樣一個事實有關,即fortran進行列專業排序而我正在進行專業排序?我會如何協調這一點?

+0

將TRANSA和TRANSB設置爲'Y'。然後dgemm會考慮輸入數組的轉置。 – janneb

+0

非常感謝,工作 – user1799323