     etaddBack += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
     //cudaMemcpy(tmp, d_rem, nr_rowsA * nr_colsB * sizeof(long long),cudaMemcpyDeviceToHost);
     //addToGMPMatrix(c, tmp, nr_rowsA, nr_colsB, 2 * exp - (i + 2) * ownLimbSize + (INT64L - 1));
+    #pragma omp parallel for schedule(dynamic)
     for(int k = 0; k < nr_rowsA; ++k){
       for(int l = 0; l < nr_colsB; ++l){
 	tmp[l *  nr_rowsA + k] = 0;

 // show memory usage of GPU

     size_t free_byte ;

     size_t total_byte ;

     cudaError_t cuda_status = cudaMemGetInfo( &free_byte, &total_byte ) ;

     if ( cudaSuccess != cuda_status ){


// Multiply the array A with its transpse. 
// B(n, n) =  A(n, k) A^T(k, n)
void gpu_blas_mulWithTransp(const cublasHandle_t handle, const double *A, double *B, const int n, const int k) {
  int lda = n;
  const double alf = 1; 
  const double bet = 0;
  const double *alpha = &alf;
  const double *beta = &bet;
  
     
  // Do the actual multiplication
  cublasDsyrk(handle, CUBLAS_FILL_MODE_LOWER, CUBLAS_OP_N, n, k, alpha, A, lda, beta, B, lda);
}


// Multiply the array A with its transpse. 
// B(n, n) =  A(n, k) A^T(k, n)
void gpu_blas_mulWithTranspAndSum(const cublasHandle_t handle, const double *A, const double *B, double *C, const int m, const int k) {
  int lda = n, ldb = k, ldc = k;
  const double alf = 1; 
  const double bet = 0;
  const double *alpha = &alf;
  const double *beta = &bet;

     
  // Do the actual multiplication
  cublasDsyrk(handle, CUBLAS_FILL_MODE_LOWER, CUBLAS_OP_N, n, k, alpha, A, lda, B, ldb, beta, C, ldc);
}




