#include <stdio.h>
#include <stdlib.h>
/*
#include <stdio.h>
#include <pthread.h>
#include <cblas.h>
#include <cublas_v2.h>
*/

#define BLOCK_SIZE 16

void dgemm_(char*, char*, int*, int*, int*, double*, double*, int*, double*, int*,
		double*, double*, int*);



void* compute_r2_space_chol_gpu(const int nO, const int nV, const int cholesky_mo_num,
                double* t1,
                double* tau,
                double* cc_space_v_vo_chol,
                double* cc_space_v_vv_chol,
                double* r2)
{
    int m,n,k, lda, ldb, ldc;
    double alpha, beta;
    double* A;
    double* B;
    double* C;

    double* tmp_cc = malloc(cholesky_mo_num*nV*nV*sizeof(double));


    m=cholesky_mo_num*nV;  n=nV;  k=nO;
    alpha=1.0; beta=0.0;
    lda=m ; ldb=k ; ldc=m;
    A=cc_space_v_vo_chol; B=t1; C=tmp_cc;
    dgemm_("N","N", &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);


    #pragma omp parallel
    {

        double* tmp_cc2 = malloc(cholesky_mo_num*nV*sizeof(double));
        double* B1 = malloc(nV*nV*BLOCK_SIZE*sizeof(double));
        double* tmpB1 = malloc(nV*BLOCK_SIZE*nV*sizeof(double));

        #pragma omp for
        for (size_t gam=0 ; gam<nV ; ++gam)
        {
           double* cc_space_v_vv_chol_ = &(cc_space_v_vv_chol[gam*nV*cholesky_mo_num]);
           double* tmp_cc_ = &(tmp_cc[gam*nV*cholesky_mo_num]);

           for (size_t l=0 ; l<nV*cholesky_mo_num ; ++l) {
             tmp_cc2[l] = cc_space_v_vv_chol_[l] - tmp_cc_[l];
           }

           for (size_t iblock=0 ; iblock<nV ; iblock += BLOCK_SIZE)
           {
                const size_t mbs = BLOCK_SIZE < nV-iblock ? BLOCK_SIZE : nV-iblock;


                alpha=-1.0; beta=0.0;
                m=nV*mbs;  n=nV;  k=cholesky_mo_num;
                A=&(tmp_cc[iblock*cholesky_mo_num*nV]); lda=cholesky_mo_num;
                B=cc_space_v_vv_chol_; ldb=cholesky_mo_num;
                C=tmpB1 ; ldc=nV*BLOCK_SIZE;
                dgemm_("T","N", &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);



                alpha=1.0; beta=1.0;
                m=nV*mbs;  n=nV;  k=cholesky_mo_num;
                A=&(cc_space_v_vv_chol[iblock*cholesky_mo_num*nV]); lda=cholesky_mo_num;
                B=tmp_cc2; ldb=cholesky_mo_num;
                C=tmpB1 ; ldc=nV*BLOCK_SIZE;
                dgemm_("T","N", &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);



                for (size_t beta=iblock ; beta<(nV < iblock+BLOCK_SIZE ? nV : iblock+BLOCK_SIZE) ; ++beta) {
                  for (size_t b=0 ; b<nV ; ++b) {
                    for (size_t a=0 ; a<nV ; ++a) {
                       B1[a + nV*(b + nV*(beta-iblock))] =
                           tmpB1[a + nV*(beta-iblock + BLOCK_SIZE*b)];
                    }
                  }
                }

                alpha=1.0; beta=1.0;
                m=nO*nO;  n=mbs;  k=nV*nV;
                A=tau; lda=nO*nO;
                B=B1 ; ldb=nV*nV;
                C=&(r2[nO*nO*(iblock + nV*gam)]); ldc=nO*nO;
                dgemm_("N","N", &m, &n, &k, &alpha, A, &lda, B, &ldb, &beta, C, &ldc);

           }
        }

        free(tmp_cc2);
        free(B1);
        free(tmpB1);
    }
    free(tmp_cc);
}