       // This case corresponds to a precision of ~20000 bits in GMP
     }
     gettimeofday(&t1, NULL);
-    cudaMemcpy(tmp, d_res, nr_rowsA * nr_rowsA * sizeof(long long), cudaMemcpyDeviceToHost);
+    cudaMemcpy(tmp, d_res, nr_rowsA * nr_colsB * sizeof(long long), cudaMemcpyDeviceToHost);
     gettimeofday(&t2, NULL);
     etBackAlloc += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
     
     gettimeofday(&t1, NULL);
-    addToGMPMatrixSymm(c, tmp, nr_rowsA, 2 * expA - (i + 2) * ownLimbSize);
+    addToGMPMatrix(c, tmp, nr_rowsA, nr_colsB, expA + expB - (i + 2) * ownLimbSize);
     gettimeofday(&t2, NULL);
     etaddBack += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
     #pragma omp parallel for schedule(dynamic)
-    for(int l = 0; l < nr_rowsA; ++l){
-      for(int k = 0; k < nr_rowsA; ++k){
+    for(int k = 0; k < nr_rowsA; ++k){
+      for(int l = 0; l < nr_colsB; ++l){
 	tmp[l *  nr_rowsA + k] = 0;
       }
     }
   }
- 
   printf("Transfer GPU = %fms\n", etTransfer);
   printf("Multiplication GPU = %fms\n", etMult);
   printf("Addition GPU = %fms\n", etAdd);

void toom3(const cublasHandle_t handle, mpf_class *c, mpf_class *a, mpf_class *b, 
	   double *d_aS, double *d_bS, double *tmpTransferLongToGMP, 
	   long long *tmp, int *sizeMatrix, int *realExpMatrix, int *signMatrix,  
	   double *d_prodRes, long long *d_res, const int nr_rowsA, const int nr_colsA, const int nr_colsB) {
}

// Implements Karatsuba algorithm for matrix multiplication c = a.a^T

void toom2Symm(const cublasHandle_t handle, mpf_class *c, mpf_class *a, 
	       double *d_aS, double *tmpTransferLongToGMP, 
	       long long *tmp, int *sizeMatrix, int *realExpMatrix, int *signMatrix,  
	       double *d_prodRes, long long *d_res, const int nr_rowsA, const int nr_colsA, const int prec) {
  mpf_class *a1, *a2;
  int sizeA1 = (size_aS / 2 + size_aS % 2), sizeA2 = size_aS / 2; 
 
  for(int i = 0; i < nr_colsA; ++i){
    for(int k = 0; k < nr_rowsA; ++k){
      if (sizeA1 + sizeA2 != abs(a[i *  nr_rowsA + k].get_mpf_t()->_mp_size))
	std::cout << "Error?" << std::endl;
      int sign =  a[i *  nr_rowsA + k].get_mpf_t()->_mp_size /
	abs(a[i *  nr_rowsA + k].get_mpf_t()->_mp_size);
      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_size = sign * sizeA1;
      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_size = sign * sizeA2;
      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_exp = a[i *  nr_rowsA + k].get_mpf_t()->_mp_exp; 
      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_exp = (a[i *  nr_rowsA + k].get_mpf_t()->_mp_exp) - sizeA1;
      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_d = &(a[i *  nr_rowsA + k].get_mpf_t()->_mp_d[sizeA2]); 
      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_d = &(a[i *  nr_rowsA + k].get_mpf_t()->_mp_d[0]);
    }
  }
  matrixMultSymm_cuBlas(handle, tmpC1, a1, d_aS, tmpTransferLongToGMP,
			tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
			nr_rowsA, nr_colsA, prec);
  matrixMultSymm_cuBlas(handle, tmpC2, a2, d_aS, tmpTransferLongToGMP,
                        tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
                        nr_rowsA, nr_colsA, prec);
  for(int i = 0; i < nr_colsA; ++i){
    for(int k = 0; k < nr_rowsA; ++k){
      a1[i *  nr_rowsA + k] += a2[i *  nr_rowsA + k]; 
    }
  }
  
  matrixMultSymm_cuBlas(handle, C, a1, d_aS, d_bS,  tmpTransferLongToGMP,
			tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
			nr_rowsA, nr_colsA, prec);
  for (int i = 0; i < nr_rowsA; ++i) {
    for (int j = 0; j < nr_rowsA; ++j) {
      C[i *  nr_rowsA + j] -= tmpC1[i *  nr_rowsA + j];
      C[i *  nr_rowsA + j] -= tmpC2[i *  nr_rowsA + j];
      C[i *  nr_rowsA + j].get_mpf_t()->_mp_exp -= sizeA1;
      tmpC2[i *  nr_rowsA + j].get_mpf_t()->_mp_exp -= size;
      C[i * nr_rowsA + j] += (tmpC1[i * nr_rowsA + j] + tmpC2[i * nr_rowsA + j]);
    }
  }  
}

// Implements Toom-3 algorithm for matrix multiplication c = a.a^T
void toom3Symm(const cublasHandle_t handle, mpf_class *c, mpf_class *a, 
	       double *d_aS, double *tmpTransferLongToGMP, 
	       long long *tmp, int *sizeMatrix, int *realExpMatrix, int *signMatrix,  
	       double *d_prodRes, long long *d_res, const int nr_rowsA, const int nr_colsA, const int prec) {
  mpf_class *a1, *a2, *a3;
  int sizeA1 = (size_aS / 3 + (size_aS % 3)), sizeA2 = (size_aS / 3), sizeA3 = size_aS / 3; 
 
  for(int i = 0; i < nr_colsA; ++i){
    for(int k = 0; k < nr_rowsA; ++k){
      if (sizeA1 + sizeA2 + sizeA3 != abs(a[i *  nr_rowsA + k].get_mpf_t()->_mp_size))
	std::cout << "Error?" << std::endl;
      int sign =  a[i *  nr_rowsA + k].get_mpf_t()->_mp_size /
	abs(a[i *  nr_rowsA + k].get_mpf_t()->_mp_size);
      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_size = sign * sizeA1;
      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_size = sign * sizeA2;
      a3[i *  nr_rowsA + k].get_mpf_t()->_mp_size = sign * sizeA3;
      
      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_exp = a[i *  nr_rowsA + k].get_mpf_t()->_mp_exp; 
      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_exp = (a[i *  nr_rowsA + k].get_mpf_t()->_mp_exp) - sizeA1;
      a3[i *  nr_rowsA + k].get_mpf_t()->_mp_exp = (a[i *  nr_rowsA + k].get_mpf_t()->_mp_exp) - sizeA1 - sizeA2;
      
      
      a1[i *  nr_rowsA + k].get_mpf_t()->_mp_d = &(a[i *  nr_rowsA + k].get_mpf_t()->_mp_d[sizeA2 + sizeA3]); 
      a2[i *  nr_rowsA + k].get_mpf_t()->_mp_d = &(a[i *  nr_rowsA + k].get_mpf_t()->_mp_d[sizeA3]);
      a3[i *  nr_rowsA + k].get_mpf_t()->_mp_d = &(a[i *  nr_rowsA + k].get_mpf_t()->_mp_d[0]);
    }
  }
  
  matrixMultSymm_cuBlas(handle, tmpC1, a1, d_aS, tmpTransferLongToGMP,
			tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
			nr_rowsA, nr_colsA, prec);
  matrixMultSymm_cuBlas(handle, tmpC5, a3, d_aS, tmpTransferLongToGMP,
                        tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
                        nr_rowsA, nr_colsA, prec);
  
  for(int i = 0; i < nr_colsA; ++i){
    for(int k = 0; k < nr_rowsA; ++k){
      a1[i *  nr_rowsA + k] += (a2[i *  nr_rowsA + k] + a3[i *  nr_rowsA + k]); 
    }
  }
  matrixMultSymm_cuBlas(handle, tmpC2, a1, d_aS, d_bS,  tmpTransferLongToGMP,
			tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
			nr_rowsA, nr_colsA, prec);
 
  for(int i = 0; i < nr_colsA; ++i){
    for(int k = 0; k < nr_rowsA; ++k){
      a1[i *  nr_rowsA + k] -= (2 * a2[i *  nr_rowsA + k]); 
    }
  }
  matrixMultSymm_cuBlas(handle, tmpC3, a1, d_aS, d_bS,  tmpTransferLongToGMP,
			tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
			nr_rowsA, nr_colsA, prec);
  
  for(int i = 0; i < nr_colsA; ++i){
    for(int k = 0; k < nr_rowsA; ++k){
      a1[i *  nr_rowsA + k] += 3 * (a2[i *  nr_rowsA + k] + a3[i * nr_rowsA + k]); 
    }
  }
  matrixMultSymm_cuBlas(handle, tmpC4, a1, d_aS, d_bS,  tmpTransferLongToGMP,
			tmp, sizeMatrix, realExpMatrix, signMatrix, d_prodRes, d_res,
			nr_rowsA, nr_colsA, prec);
  
  for (int i = 0; i < nr_rowsA; ++i) {
    for (int j = 0; j < nr_rowsA; ++j) {
      tmpC2[i * nr_rowsA + j] = (- tmpC1[i * nr_rowsA + j] / 2 + tmpC2[i * nr_rowsA + j] 
				 - tmpC3[i * nr_rowsA + j] /3 - tmpC4[i * nr_rowsA + j] / 6 
				 + 2 * tmpC5[i * nr_rowsA + j]); 
      tmpC3[i * nr_rowsA + j] = (-3/4 * tmpC1[i * nr_rowsA + j] + tmpC2[i * nr_rowsA + j] / 2 
				 + 2/3 * tmpC3[i * nr_rowsA + j] + 1/12 * tmpC4[i * nr_rowsA + j] 
				 - 2 * tmpC5[i * nr_rowsA + j]);
      tmpC4[i * nr_rowsA + j] = (- tmpC1[i * nr_rowsA + j] / 8 - tmpC2[i * nr_rowsA + j] / 2 
				 - tmpC3[i * nr_rowsA + j] / 4 + tmpC4[i * nr_rowsA + j] / 8 
				 - 2 * tmpC5[i * nr_rowsA + j]);
  
      tmpC2[i *  nr_rowsA + j].get_mpf_t()->_mp_exp -= sizeA;
      tmpC3[i *  nr_rowsA + j].get_mpf_t()->_mp_exp -= 2 * sizeA;
      tmpC4[i *  nr_rowsA + j].get_mpf_t()->_mp_exp -= 3 * sizeA;
      tmpC5[i *  nr_rowsA + j].get_mpf_t()->_mp_exp -= 4 * sizeA;
      
      C[i *  nr_rowsA + j] = (tmpC1[i *  nr_rowsA + j] + tmpC2[i *  nr_rowsA + j] 
			      + tmpC3[i *  nr_rowsA + j] + tmpC4[i *  nr_rowsA + j] 
			      + tmpC5[i *  nr_rowsA + j]);
    }
  }
}


// Returns matrix product c = a.a^T where each entry is of the type mpf_class
//
// All arrays need to already have allocated memory:
// handle : Handle for cuBlas computations that needs to be previously allocated.                                                                   
// *c : array of mpf_class with size nr_rowsA * nr_rowsA                  
// *a : array of mpf_class with size nr_rowsA * nr_colsA
// *d_aS: array allocated in GPU with size sizeOfArray * nr_rowsA * nr_colsA. Note that we flatten this               
//        array of matrices in order to speed up the access in the GPU.                                              
//        Thus, to access the k-th limb from the (i, j) matrix entry one calls d_aS[k * nr_rowsA * nr_colsA + j * nr_rowsA + i]
// *tmpTransferLongToGMP : temporary array of doubles used for transfers that needs nr_rowsA * nr_colsA entries                    
// *tmp: temporary array of 64-bit vars used for transfers. Needs to be allocated nr_rowsA * nr_rowsA entries
// *d_prodRes: temporary matrix allocated in GPU memory which handles individual multiplications of matrices
// *d_res : temporary matrix allocated in GPU memory in which we sum up all individual multiplications of matrices
// *d_rem : in case we encounter overflow in the matrix d_res we use the entries in the matrix d_rem to store the remainders. 
//          Note that this matrix is only used when using the function vecAdd__wRem. When using vecAdd__wSign we assume that 
//          overflow never occurs. This is the case if the number of own generated limbs is greater than 1024. This correpons 
//          to a GMP precision ~20,000. For current bootstrap applications such precision is not needed.
// nr_rows_A: number of rows for matrix A as well as for matrix C                                                                                
// nr_cols_A: number of columns for matrix A as well as number of rows for matrix B
// nr_cols_B: number of columns for matrix B as well as for matrix C                                                                            
// maxExp:  maximum power of 2 for the leading most limb among all the entries of the array *a 
void matrixMultSymmArbOrder_cuBlas(const cublasHandle_t handle, mpf_class *c, mpf_class *a, 
				   double *d_aS, double *tmpTransferLongToGMP, 
				   long long *tmp, int *sizeMatrix, int *realExpMatrix, int *signMatrix,  
				   double *d_prodRes, long long *d_res, const int nr_rowsA, const int nr_colsA, 
				   const int size_aS, const int expA, const int ownLimbSize, const int prec) {
  timeval t1, t2;
  
  for(int i = 0; i < nr_rowsA; ++i){
    for(int j = 0; j < nr_rowsA ++j){
      c[j * nr_rowsA +  i] = mpf_class("0.0");
      tmp[j * nr_rowsA +  i] = 0;
    }
  }

  // Number of threads in each thread block                                                                          
  int blockSize = 256; // Make sure this number is good? I don't know how
  // Number of thread blocks in grid                                                                                 
  int gridSize = (int)ceil((float)(nr_rowsA * nr_rowsA)/blockSize);
  
  double etMult = 0;
  double etAdd = 0;
  double etAlloc = 0;
  double etBackAlloc = 0;
  double etaddBack = 0;
  
  long sizeToCompute =  (whatOrder + expA) / ownLimbSize;
  
  for (int i = 0; i < min(size, sizeToCompute); i++) {
    gettimeofday(&t1, NULL);
    cudaMemcpy(d_res, tmp, nr_rowsA * nr_rowsA * sizeof(long long), cudaMemcpyHostToDevice);
    gettimeofday(&t2, NULL);
    etAlloc += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
    //cudaMemcpy(d_rem, tmp, nr_rowsA * nr_colsB * sizeof(long long), cudaMemcpyHostToDevice);
    for (int j = 0; j < i + 1; j++) {
      gettimeofday(&t1, NULL);
      if (i - j != j)
	gpu_blas_mulWithTranspAndSum(handle, &d_aS[j * nr_rowsA * nr_rowsA], &d_aS[(i - j) * nr_colsA * nr_colsB], d_prodRes, nr_rowsA, nr_colsA); 
      else
	gpu_blas_mulWithTransp(handle, &d_aS[j * nr_rowsA * nr_rowsA], d_prodRes, nr_rowsA, nr_colsA); 
      cudaThreadSynchronize();
      gettimeofday(&t2, NULL);
      etMult += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
      
      gettimeofday(&t1, NULL);
      vecAdd__wSign<<<gridSize, blockSize>>>(d_prodRes, d_res,  nr_rowsA * nr_rowsA);
      cudaThreaMdSynchronize();
      gettimeofday(&t2, NULL);
      etAdd += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
      // This is safe with overflow until there are 1024 of our own limbs that need to be summed up
      // This case corresponds to a precision of ~20000 bits in GMP
    }
    gettimeofday(&t1, NULL);
    cudaMemcpy(tmp, d_res, nr_rowsA * nr_rowsA * sizeof(long long), cudaMemcpyDeviceToHost);
    gettimeofday(&t2, NULL);
    etBackAlloc += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
    
    gettimeofday(&t1, NULL);
    addToGMPMatrixSymm(c, tmp, nr_rowsA, 2 * expA - (i + 2) * ownLimbSize);
    gettimeofday(&t2, NULL);
    etaddBack += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
    
    #pragma omp parallel for schedule(dynamic)
    for(int l = 0; l < nr_rowsA; ++l){
      for(int k = 0; k < nr_rowsA; ++k){
	tmp[l *  nr_rowsA + k] = 0;
      }
    }
  }
  
   
  for (int i = size; i < min(2 * size - 1, sizeToCompute); i++) {
    gettimeofday(&t1, NULL);
    cudaMemcpy(d_res, tmp, nr_rowsA * nr_rowsA * sizeof(long long), cudaMemcpyHostToDevice);
    gettimeofday(&t2, NULL);
    etAlloc += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
    //cudaMemcpy(d_rem, tmp, nr_rowsA * nr_colsB * sizeof(long long), cudaMemcpyHostToDevice);
    for (int j = i - size + 1; j < size; j++) {
      gettimeofday(&t1, NULL);
      if (i - j != j)
	gpu_blas_mulWithTranspAndSum(handle, &d_aS[j * nr_rowsA * nr_rowsA], &d_aS[(i - j) * nr_colsA * nr_colsB], d_prodRes, nr_rowsA, nr_colsA); 
      else
	gpu_blas_mulWithTransp(handle, &d_aS[j * nr_rowsA * nr_rowsA], d_prodRes, nr_rowsA, nr_colsA); 
      cudaThreadSynchronize();
      gettimeofday(&t2, NULL);
      etMult += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
      
      gettimeofday(&t1, NULL);
      vecAdd__wSign<<<gridSize, blockSize>>>(d_prodRes, d_res,  nr_rowsA * nr_rowsA);
      cudaThreaMdSynchronize();
      gettimeofday(&t2, NULL);
      etAdd += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
      // This is safe with overflow until there are 1024 of our own limbs that need to be summed up
      // This case corresponds to a precision of ~20000 bits in GMP
    }
    gettimeofday(&t1, NULL);
    cudaMemcpy(tmp, d_res, nr_rowsA * nr_rowsA * sizeof(long long), cudaMemcpyDeviceToHost);
    gettimeofday(&t2, NULL);
    etBackAlloc += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
    
    gettimeofday(&t1, NULL);
    addToGMPMatrixSymm(c, tmp, nr_rowsA, 2 * expA - (i + 2) * ownLimbSize);
    gettimeofday(&t2, NULL);
    etaddBack += (((t2.tv_sec*uS_PER_SEC)+t2.tv_usec) - ((t1.tv_sec*uS_PER_SEC)+t1.tv_usec))/(float)uS_PER_mS;
    #pragma omp parallel for schedule(dynamic)
    for(int l = 0; l < nr_rowsA; ++l){
      for(int k = 0; k < nr_rowsA; ++k){
	tmp[l *  nr_rowsA + k] = 0;
      }
    }
  }
 
  printf("Transfer GPU = %fms\n", etTransfer);
  printf("Multiplication GPU = %fms\n", etMult);
  printf("Addition GPU = %fms\n", etAdd);
  printf("Alloc of zero to GPU = %fms\n", etAlloc);
  printf("Transfer from GPU to host = %fms\n", etBackAlloc);
  printf("Addition on CPU = %fms\n", etaddBack);
}


// Returns matrix product c = a.a^T where each entry is of the type mpf_class
//
