Logo Questions Linux Laravel Mysql Ubuntu Git Menu
HTML CSS JAVASCRIPT SQL PYTHON PHP BOOTSTRAP JAVA JQUERY R React Kotlin
 

Accessing dgemm in C++

Tags:

c++

fortran

I need to implement a matrix class in C++ and one of the operations must be matrix multiplication via dgemm. My professor had done an example in class in C, but for some reason I can't get it to work in C++. This is my header file matrix.h:

#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]

Main function:

#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;

}

Now the blas.h header file:

void dgemm_ (char *TRANSA,
             char *TRANSB,
             int *M,
             int *N,
             int *K,
             double *ALPHA,
             double **A,
             int *LDA,
             double **B,
             int *LDB,
             double *BETA,
             double **C,
             int *LDC);

I'm to use the subroutine outlined here: http://www.netlib.org/blas/dgemm.f

Basically, the way we did it in class for the C implementation built a matrix as one long array using (double*) calloc(rows*columns, sizeof(double)).

My C++ implemenation does:

double **data;

data = new double[rows];

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

and then I can index using data[i][j]. But since the dgemm subroutine is supposed to take in double *A, double *B, etc. but my matrices are double **A and so on, how do I fix this?

This is the error message I get from 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)

up until I try to implement the dgemm matrix multiplication, I get no errors and no leaks, so I'm pretty sure all my troubles lie with the dgemm implementation

like image 905
user1799323 Avatar asked Nov 12 '22 18:11

user1799323

1 Answers

The memory layout of the C version is different from the memory layout in your C++ version. As BLAS expects the kind of layout used by the C version, your C++ version won't work.

Thus, you need to allocate a big 1-D array in your C++ version as well; you could overload operator() in order to get indexing for 2D arrays. Or for production code, use a library such as Eigen.

like image 62
janneb Avatar answered Nov 15 '22 13:11

janneb
Sign in to Comment

Related questions

Improving efficiency of std::copy() with back_inserter() or inserter()
Finding the most popular likes in my friend network
Passing in function pointer as argument to another function in std::bind
Win32 - determine if path has jumbo icon available
Resuming an a-star search in BGL
recursive permutation generator for characters [duplicate]
utc seconds since midnight to datetime
Portable way of accessing an array of bool (C++11) from Nvidia PTX
OpenCV not initializing USB camera
C++ error distinguishing specialization when passing a const array
How to implement a back-off with Microsoft PPL lightweight task scheduler?
boost::unique_lock and boost::shared_lock for reader writer locks
CancelEvent equivalent in Qt signals
Implementing the DllRegisterServer method
How to use size of in __align__ for cache alignment, while using GCC?
Should std::atomic<int*>::load be doing a compare-and-swap loop?
C++ Library for modifying a ZIP file in place
Union – useless anachronism or useful old school trick?
What is the difference between false and FALSE?

Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!