simple matrix multiplication result error using cuBLAS lib

seehakoo · December 19, 2009, 11:07pm

hi all,

I got some problem about using cublas library.  Now, I have two matrices to multiply, but result is error.

I can’t find out problem in my code. Any somebody can help me, please.

Matrix P is : (NxN)

[0] [I]

[0] [-J]

[J] [0]

[I] [0]

here [0] is zero matrix, [I] is identity,

[J] = 0 0 0 0 1

[0], [I] and [J] both is (N/4 x N/4) matrices.

here is my code

[codebox]

#include <math.h>

#include <cuda.h>

#include <cuda_runtime.h>

#include

#include <cutil.h>

#include <cutil_inline.h>

#include “cublas.h”

using namespace std;

#define IDX2C(i,j,ld) (((j)*(ld))+(i))

/* Allocate host memory for the matrices */\

#define HOSTmalloc(H_ptr,N, id);{\

H_ptr = (float*)malloc(N * sizeof(H_ptr[0]));\

if (H_ptr == 0) {\

	fprintf (stderr, "!!!! (%dth)host memory allocation error.\n", id);\

	return EXIT_FAILURE;\

}\

memset(H_ptr, 0, N * sizeof(H_ptr[0]));\

}

int first_algo1();

int LOG_N;

long N;

float* h_S;

float* h_P;

unsigned int n2;

unsigned int np;

int main(int argc, char* argv)

{

LOG_N = 6;

N = 1;

N = N << LOG_N;

/****************  start ************/

n2 = N/2 * N/2;

np = N * N/2;

/* Allocate host memory for the matrices */

HOSTmalloc(h_S, n2, 1);

HOSTmalloc(h_P, np, 2);

long n4 = N/2 * N/2;

float* I;

float* J;

HOSTmalloc(I, n4, 7);

HOSTmalloc(J, n4, 8);

/* Fill the matrices with test data */

/*!!!!!!!   column major !!!!!!! */

for(long i=0; i < N/2; i++){

	for(long j=0; j< N/2; j++){

		long idx = IDX2C(i, j, N/2);

		h_S[idx] = 1;

	}

}

/*

 *      |  [0] [I]   |

 *   P= |  [0] [-J]  |

 *      |  [J] [0]   |

 *      |  [-I] [0]   | N by N  matrix

 *

 *

 *  I = | 1 0 0 |

 *      | 0 : 0 |

 *      | 0 0 1 |

 *

 *  J = | 0 0 1 |

 *      | 0 : 0 |

 *      | 1 0 0 |

 * */

for(long i=0; i < N/4; i++){

	for(long j=0; j< N/4; j++){

		long idx =IDX2C(i, j, N/4);

		I[idx] = 0.0f;

		J[idx] = 0.0f;

		if(i == j)

			I[idx] = 1.0;

		if((i+j) == N/4-1)

			J[idx] = 1.0;

	}

}

for(long i=0; i < N; i++){

	for(long j=0; j< N/2; j++){

		long idx = IDX2C(i, j, N);

		h_P[idx] = 0.0f;

		if(j>= N/4){

			if(i < N/4){

				//h_P[idx] = I[i][j-N/4];

				long Midx = IDX2C(i, j-N/4, N/4);

				h_P[idx] = I[Midx];

			}else if(i >= N/4 && i < N/2){

				//h_P[idx] = -1*J[i-N/4][j-N/4];

				long Midx = IDX2C(i-N/4, j-N/4, N/4);

				if(J[Midx] != 0)

					h_P[idx] = -1*J[Midx];

			}

		}else if(j < N/4){

			if(i >= (N/2) && i < 3*N/4){

				//h_P[idx] = J[i-N/2][j];

				long Midx = IDX2C(i-N/2, j, N/4);

				h_P[idx] = J[Midx];

			}else if(i >= 3*N/4){

				//h_P[idx] = I[i-3*N/4][j];

				long Midx = IDX2C(i-3*N/4, j, N/4);

				h_P[idx] = -I[Midx];

			}

		}

	}

}

first_algo1();

free(h_S);

free(h_P);

return 0;

}

#define CUBLASmalloc(N, D_ptr, ID);{\

/* Allocate device memory for the matrices */\

status = cublasAlloc(N, sizeof(D_ptr[0]), (void**)&D_ptr);\

if (status != CUBLAS_STATUS_SUCCESS) {\

	fprintf (stderr, "!!!!  (%dth)device memory allocation error\n",ID);\

	return EXIT_FAILURE;\

}\

}

#define CUBLASset(N, H_ptr, D_ptr, ID);{\

/* Initialize the device matrices with the host matrices */\

status = cublasSetVector(N, sizeof(H_ptr[0]), H_ptr, 1, D_ptr, 1);\

if (status != CUBLAS_STATUS_SUCCESS) {\

	fprintf (stderr, "!!!! device access error (write :%dth)\n", ID);\

	return EXIT_FAILURE;\

}\

}

#define CUBLASget(N, H_ptr, D_ptr);{\

/* Read the result back */\

status = cublasGetVector(N, sizeof(H_ptr[0]), D_ptr, 1, H_ptr, 1);\

if (status != CUBLAS_STATUS_SUCCESS) {\

    fprintf (stderr, "!!!! device access error (read C)\n");\

    return EXIT_FAILURE;\

}\

}

#define CUBLASfree(D_ptr, id);{\

status = cublasFree(D_ptr);\

if (status != CUBLAS_STATUS_SUCCESS) {\

    fprintf (stderr, "!!!! memory free error (%dth)\n", id);\

    return EXIT_FAILURE;\

}\

}

int first_algo1()

{

cublasStatus status;

float* d_S = 0;

float* d_P = 0;

float alpha = 1.0f;

float beta = 0.0f;

/* Initialize CUBLAS */

status = cublasInit();

if (status != CUBLAS_STATUS_SUCCESS) {

    fprintf (stderr, "!!!! CUBLAS initialization error\n");

    return EXIT_FAILURE;

}

status = cublasGetError();

if (status != CUBLAS_STATUS_SUCCESS) {

    fprintf (stderr, "!!!! kernel execution error.\n");

    return EXIT_FAILURE;

}

/***************************/

// x = P * S * L * DCT' * D * X

CUBLASmalloc(np, d_P, 5);

CUBLASset(np, h_P, d_P, 5);

CUBLASmalloc(n2, d_S, 3);

CUBLASset(n2, h_S, d_S, 3);

// P = P * S (n by N/2)

cublasSgemm('n', 'n', N, N/2, N/2, alpha, d_P, N, d_S, N/2, beta, d_P, N);    // [A]_m*k  X  [B]_k*n  -->  ('n', 'n', m, n, k, alpha, A, m, B, k, beta, C, m)   // is column major

status = cublasGetError();

if (status != CUBLAS_STATUS_SUCCESS) {

    fprintf (stderr, "!!!! kernel execution error.\n");

    return EXIT_FAILURE;

}

float *tmp;

HOSTmalloc(tmp, np, 3);

int i, j;

cout << "Get Value:" << endl;

CUBLASget(np, tmp, d_P);

for(i=0; i< N; i++){

	for(j=0; j< N/2; j++){

		int idx = IDX2C(i, j, N);

		cout << " " << tmp[idx];

	}

	cout << endl;

}

free(tmp);

CUBLASfree(d_S, 1);

CUBLASfree(d_P, 2);

/* Shutdown */

status = cublasShutdown();

if (status != CUBLAS_STATUS_SUCCESS) {

    fprintf (stderr, "!!!! shutdown error (A)\n");

    return EXIT_FAILURE;

}

return EXIT_SUCCESS;

}

[/codebox]

And here is my result ( when my LOG_N = 6) ( when LOG_N < 6 is no error)

[codebox]

Get Value:

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16 16

-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -16 -16 -16 -16 -16 -16 -16 -16 -16 -16 -16 -16 -16 -16 -16 -16

[/codebox]

here is strange about right and down side of the matrix.

when LOG_N = 7;

[codebox]

Get Value:

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33 33

-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33 -33

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32 32

-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32 -32

[/codebox]

the right side of the matrix is error.

here is my device intomation

[codebox]

Device 0: “GeForce 9300M G”

CUDA Driver Version: 2.30

CUDA Runtime Version: 2.30

CUDA Capability Major revision number: 1

CUDA Capability Minor revision number: 1

Total amount of global memory: 267714560 bytes

Number of multiprocessors: 2

Number of cores: 16

Total amount of constant memory: 65536 bytes

Total amount of shared memory per block: 16384 bytes

Total number of registers available per block: 8192

Warp size: 32

Maximum number of threads per block: 512

Maximum sizes of each dimension of a block: 512 x 512 x 64

Maximum sizes of each dimension of a grid: 65535 x 65535 x 1

Maximum memory pitch: 262144 bytes

Texture alignment: 256 bytes

Clock rate: 0.80 GHz

Concurrent copy and execution: Yes

Run time limit on kernels: Yes

Integrated: No

Support host page-locked memory mapping: No

Compute mode: Default (multiple host threads can use this device simultaneously)

[/codebox]

Thanks everyone!

avidday · December 20, 2009, 10:12am

Are these matrices stored in column or row major order? CUBLAS requires data stored in column order (ie. like Fortran) rather than row ordered (like C).

seehakoo · December 20, 2009, 12:39pm

I use column major order, and refer to “CUBLAS_Library_2.3.pdf”. So I use macro define “IDX2C(i,j,ld) (((j)*(ld))+(i))” in my code.

More information, my OS is ubuntu 9.04.

And what something I miss or lost?

Thinks your help.