#line 64 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
#ifndef __LMMSE_H__
#include "LMMSE.h"
#endif

#line 641 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
__global__ void zeroPadding(float *d__zp_x, int NP, float *d__x, int N_SIDE_LENSLET, int osf)
{
int i, j, ip, jp, k, kp, iSource, NP2, N_LENSLET;
i = blockIdx.x * blockDim.x + threadIdx.x;
j = blockIdx.y * blockDim.y + threadIdx.y;
iSource = blockIdx.z;
if ( (i<N_SIDE_LENSLET) && (j<N_SIDE_LENSLET) ){
    NP2 = NP*NP;
    N_LENSLET = N_SIDE_LENSLET*N_SIDE_LENSLET; 
    ip = osf*i + osf/2;
    jp = osf*j + osf/2;
    k = i*N_SIDE_LENSLET + j;
    k += N_LENSLET*2*iSource;
    kp = ip*NP + jp;
    kp += NP2*2*iSource;
    d__zp_x[kp] = d__x[k];
    d__zp_x[kp+NP2] = d__x[k+N_LENSLET];
}
}
#line 738 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
__global__ void bilinearSparseOperator(float *csrValH, int *csrColIndH, int *csrRowPtrH,
				  int* MT, int NI, layer *turb, int N_LAYER, float z)
{
  int i, j, k, idx, ndx, NP, k_LAYER, offset;
  float scale, s, t, fs, ft, onemt, onems, delta, delta_p;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (i<NI) && (j<(NI)) ) {
    delta   = 1.0;
    k = i*NI + j;
    idx = N_LAYER*4*k - 1;
    offset = 0;
    for (k_LAYER=0;k_LAYER<N_LAYER;k_LAYER++) {
    	delta_p = 1 - turb[k_LAYER].altitude/z;
	NP = MT[k_LAYER];
        
#line 709 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
scale = delta/(delta_p*(NP-1));
s = (scale*(i + (1-NI)/2) + 0.5)*(NP-1);
t = (scale*(j + (1-NI)/2) + 0.5)*(NP-1);
#line 714 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
fs  = floorf(s);
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*NP );
#line 719 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
if (s==(NP-1)) { s += 1 - fs; ndx -= NP; } else { s -= fs; }    
if (t==(NP-1)) { t += 1 - ft; ndx -= 1; }   else { t -= ft; }
#line 723 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
onems = 1 - s;
onemt = 1 - t;
ndx += offset;
csrValH[++idx]  = onems*onemt;
csrColIndH[idx] = ndx;
csrValH[++idx]  = onems*t;
csrColIndH[idx] = ndx + 1;
csrValH[++idx]  = s*onemt;
csrColIndH[idx] = ndx + NP;
csrValH[++idx]  = s*t;
csrColIndH[idx] = ndx + NP + 1;


#line 754 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
	offset += NP*NP;
    }	
    csrRowPtrH[k+1] = N_LAYER*4*(k+1);
  }
}
#line 762 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
__global__ void bilinearSparseOperatorMask(float *csrValH, int *csrColIndH, int *csrRowPtrH,
				  int* MT, int NI, layer *turb, int N_LAYER, float z, 
                                  const char *mask, int nnz)
{
  int i, j, k, idx, ndx, NP, k_LAYER, offset, l, pos;
  float scale, s, t, fs, ft, onemt, onems, delta, delta_p;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (i<NI) && (j<(NI)) ) {
    delta   = 1.0;
    k = i*NI + j;
    if (mask[k]>0) {
       
#line 798 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
pos = 0;
for (l=0;l<=k;l++)
    pos += (mask[l]) ? 1 : 0;
#line 775 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
       idx =  N_LAYER*4*(pos-1) - 1;
       offset = 0;
       for (k_LAYER=0;k_LAYER<N_LAYER;k_LAYER++) {
         delta_p = 1 - turb[k_LAYER].altitude/z;
         NP = MT[k_LAYER];
         
#line 709 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
scale = delta/(delta_p*(NP-1));
s = (scale*(i + (1-NI)/2) + 0.5)*(NP-1);
t = (scale*(j + (1-NI)/2) + 0.5)*(NP-1);
#line 714 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
fs  = floorf(s);
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*NP );
#line 719 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
if (s==(NP-1)) { s += 1 - fs; ndx -= NP; } else { s -= fs; }    
if (t==(NP-1)) { t += 1 - ft; ndx -= 1; }   else { t -= ft; }
#line 723 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
onems = 1 - s;
onemt = 1 - t;
ndx += offset;
csrValH[++idx]  = onems*onemt;
csrColIndH[idx] = ndx;
csrValH[++idx]  = onems*t;
csrColIndH[idx] = ndx + 1;
csrValH[++idx]  = s*onemt;
csrColIndH[idx] = ndx + NP;
csrValH[++idx]  = s*t;
csrColIndH[idx] = ndx + NP + 1;


#line 781 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
         offset += NP*NP;
       }	
       csrRowPtrH[k] = N_LAYER*4*(pos-1);
       
#line 804 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
l = k-1;
 while ( (l>=0) && (!mask[l]) )
    csrRowPtrH[l--] = csrRowPtrH[k];



#line 785 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
       if ( k==(NI*NI-1) )
          csrRowPtrH[k+1] = N_LAYER*4*nnz;
    } else { 
         if ( k==(NI*NI-1) ) {
           ++k;
           csrRowPtrH[k] = N_LAYER*4*nnz;
           
#line 804 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
l = k-1;
 while ( (l>=0) && (!mask[l]) )
    csrRowPtrH[l--] = csrRowPtrH[k];



#line 792 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
        } 
     }
   }
 }
#line 953 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
__global__ void bilinearInterpSparseOperator(float *csrValH, int *csrColIndH, 
                                             int *csrRowPtrH, int NP, int NI)
{
  int i, j, k, idx, ndx, offset;
  float scale, s, t, fs, ft, onemt, onems, delta, delta_p;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (i<NI) && (j<(NI)) ) {
    delta   = 1.0;
    k = i*NI + j;
    idx = 4*k - 1;
    offset = 0;

    delta_p = delta;
    
#line 1053 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
scale = delta/(delta_p*(NP-1));
s = (scale*(i + (1-NI)/2) + 0.5)*(NP-1);
t = (scale*(j + (1-NI)/2) + 0.5)*(NP-1);
#line 1058 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
fs  = floorf(s);
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*NP );
#line 1063 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
if (s==(NP-1)) { s += 1 - fs; ndx -= NP; } else { s -= fs; }    
if (t==(NP-1)) { t += 1 - ft; ndx -= 1; }   else { t -= ft; }
#line 1067 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
onems = 1 - s;
onemt = 1 - t;
ndx += offset;
csrValH[++idx]  = onems*onemt;
csrColIndH[idx] = ndx;
csrValH[++idx]  = onems*t;
csrColIndH[idx] = ndx + 1;
csrValH[++idx]  = s*onemt;
csrColIndH[idx] = ndx + NP;
csrValH[++idx]  = s*t;
csrColIndH[idx] = ndx + NP + 1;

#line 969 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
    csrRowPtrH[k+1] = 4*(k+1);
  }
}
#line 975 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
__global__ void bilinearInterpSparseOperatorMask(float *csrValH, int *csrColIndH, 
                                                  int *csrRowPtrH, int NP, int NI,
                                                  const char *mask, int nnz,
                                                  const float i0, const float j0)
{
  int ii, jj, k, idx, ndx, offset, l, pos;
  float scale, s, t, fs, ft, onemt, onems, delta, delta_p, i, j;
  ii = blockIdx.x * blockDim.x + threadIdx.x;
  jj = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (ii<NI) && (jj<(NI)) ) {
    delta   = 1.0;
    k = ii*NI + jj;
    if (mask[k]>0) {
       i = ii - i0;
       j = jj - j0;
       
#line 798 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
pos = 0;
for (l=0;l<=k;l++)
    pos += (mask[l]) ? 1 : 0;
#line 991 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
       idx =  4*(pos-1) - 1;
       offset = 0;

         delta_p = delta;

         
#line 1053 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
scale = delta/(delta_p*(NP-1));
s = (scale*(i + (1-NI)/2) + 0.5)*(NP-1);
t = (scale*(j + (1-NI)/2) + 0.5)*(NP-1);
#line 1058 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
fs  = floorf(s);
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*NP );
#line 1063 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
if (s==(NP-1)) { s += 1 - fs; ndx -= NP; } else { s -= fs; }    
if (t==(NP-1)) { t += 1 - ft; ndx -= 1; }   else { t -= ft; }
#line 1067 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
onems = 1 - s;
onemt = 1 - t;
ndx += offset;
csrValH[++idx]  = onems*onemt;
csrColIndH[idx] = ndx;
csrValH[++idx]  = onems*t;
csrColIndH[idx] = ndx + 1;
csrValH[++idx]  = s*onemt;
csrColIndH[idx] = ndx + NP;
csrValH[++idx]  = s*t;
csrColIndH[idx] = ndx + NP + 1;
#line 997 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
         offset += NP*NP;

       csrRowPtrH[k] = 4*(pos-1);
       
#line 804 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
l = k-1;
 while ( (l>=0) && (!mask[l]) )
    csrRowPtrH[l--] = csrRowPtrH[k];



#line 1001 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
       if ( k==(NI*NI-1) )
          csrRowPtrH[k+1] = 4*nnz;
    } else { 
         if ( k==(NI*NI-1) ) {
           ++k;
           csrRowPtrH[k] = 4*nnz;
           
#line 804 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
l = k-1;
 while ( (l>=0) && (!mask[l]) )
    csrRowPtrH[l--] = csrRowPtrH[k];



#line 1008 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
        } 
     }
   }
 }

#line 920 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void bilinearInterpolation::setup(int _NI_, int NP)
{
NI = _NI_;
nnz = 4*NI*NI;
#line 838 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
alpha = 1;
beta  = 0;
HANDLE_ERROR_CUSPARSE( cusparseCreate(&handle), 
		       "CUSPARSE Library initialization failed!");
HANDLE_ERROR_CUSPARSE( cusparseCreateMatDescr(&descr), 
		       "Matrix descriptor initialization failed!");

cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr,CUSPARSE_INDEX_BASE_ZERO);

HANDLE_ERROR( cudaMalloc((void**)&csrValH,    sizeof(float)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrColIndH, sizeof(int)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrRowPtrH, sizeof(int)*(NI*NI+1) ) );
HANDLE_ERROR( cudaMemset(csrRowPtrH, 0, sizeof(int)*(NI*NI+1) ) );

HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_i, sizeof(float)*NI*NI ) );

#line 940 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearInterpSparseOperator <<< gridDim,blockDim >>> 
  (csrValH, csrColIndH, csrRowPtrH, NP, NI);
#line 926 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}
#line 930 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void bilinearInterpolation::setup(int _NI_, int NP, mask *pupil,
                                  float i0, float j0)
{
NI = _NI_;
nnz = 4*pupil->nnz;
#line 838 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
alpha = 1;
beta  = 0;
HANDLE_ERROR_CUSPARSE( cusparseCreate(&handle), 
		       "CUSPARSE Library initialization failed!");
HANDLE_ERROR_CUSPARSE( cusparseCreateMatDescr(&descr), 
		       "Matrix descriptor initialization failed!");

cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr,CUSPARSE_INDEX_BASE_ZERO);

HANDLE_ERROR( cudaMalloc((void**)&csrValH,    sizeof(float)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrColIndH, sizeof(int)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrRowPtrH, sizeof(int)*(NI*NI+1) ) );
HANDLE_ERROR( cudaMemset(csrRowPtrH, 0, sizeof(int)*(NI*NI+1) ) );

HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_i, sizeof(float)*NI*NI ) );

#line 946 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearInterpSparseOperatorMask <<< gridDim,blockDim >>> 
  (csrValH, csrColIndH, csrRowPtrH, NP, NI, 
  pupil->m, pupil->nnz, i0, j0);
#line 937 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}
#line 1081 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void bilinearInterpolation::cleanup(void)
{
#line 891 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
HANDLE_ERROR_CUSPARSE( cusparseDestroyMatDescr(descr),
		       "Matrix descriptor destruction failed!");
HANDLE_ERROR_CUSPARSE( cusparseDestroy(handle),
		       "CUSPARSE Library release of resources failed!");
HANDLE_ERROR( cudaFree( csrValH ) );
HANDLE_ERROR( cudaFree( csrColIndH ) );
HANDLE_ERROR( cudaFree( csrRowPtrH ) );
HANDLE_ERROR( cudaFree( d__phase_est_i ) );
#line 1084 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}

#line 351 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::setup(atmosphere *atm, source *guide_star, source *mmse_star,
		  float sampling, int _N_SIDE_LENSLET_, char *solver_id)
{
  N_SIDE_LENSLET_ = _N_SIDE_LENSLET_;
  osf = 1; // The wavefront over-sampling factor (1 means the same sampling that the WFS)
  
#line 451 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  N_guide_star = guide_star->N_SRC;
  N_mmse_star  = mmse_star->N_SRC;

  
#line 528 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aa.setup(N_SIDE_LENSLET_,atm,sampling,guide_star, N_guide_star);
#line 531 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aaCov.setup(2*N_guide_star,2*N_guide_star,N_SIDE_LENSLET_,aa.d__cov);

#line 456 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  NP = osf*N_SIDE_LENSLET_+1; // The number of wavefront sample
  NS = (osf>1) ? NP : N_SIDE_LENSLET_;

  pa.setup(NP,NS,osf,atm,sampling,
           mmse_star, N_mmse_star, guide_star, N_guide_star);

#line 464 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  if (mmse_star[0].height==guide_star[0].height) {

     paCov.setup(N_mmse_star,2*N_guide_star,pa.M_LAYER,N_SIDE_LENSLET_,pa.d__cov);

     PS_E_N_PX = paCov.MT2_TOTAL;
     nnz = 0;

  } else {
 
    paCov.setup(N_mmse_star*atm->N_LAYER,2*N_guide_star,pa.M_LAYER,N_SIDE_LENSLET_,pa.d__cov);
  
     PS_E_N_PX = paCov.MT2_TOTAL;

     NI = NP;
     nnz = 4*NI*NI*atm->N_LAYER;
     
#line 838 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
alpha = 1;
beta  = 0;
HANDLE_ERROR_CUSPARSE( cusparseCreate(&handle), 
		       "CUSPARSE Library initialization failed!");
HANDLE_ERROR_CUSPARSE( cusparseCreateMatDescr(&descr), 
		       "Matrix descriptor initialization failed!");

cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr,CUSPARSE_INDEX_BASE_ZERO);

HANDLE_ERROR( cudaMalloc((void**)&csrValH,    sizeof(float)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrColIndH, sizeof(int)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrRowPtrH, sizeof(int)*(NI*NI+1) ) );
HANDLE_ERROR( cudaMemset(csrRowPtrH, 0, sizeof(int)*(NI*NI+1) ) );

HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_i, sizeof(float)*NI*NI ) );

#line 857 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearSparseOperator <<< gridDim,blockDim >>> 
  (csrValH, csrColIndH, csrRowPtrH, paCov.d__MT, NI, 
   atm->d__layers, atm->N_LAYER, guide_star[0].height);
#line 872 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
/*
// For debugging purposes:

dev2file("csrValH.bin",csrValH,nnz);
dev2file("csrColIndH.bin",csrColIndH,nnz);
dev2file("csrRowPtrH.bin",csrRowPtrH,NI*NI+1);

float *H;
HANDLE_ERROR( cudaMalloc( (void**)&H, sizeof(float)*NI*NI*paCov.MT2_TOTAL ) );
HANDLE_ERROR_CUSPARSE(cusparseScsr2dense(handle, NI*NI, paCov.MT2_TOTAL, descr, 
					 csrValH, csrRowPtrH, csrColIndH,
					 H, NI*NI),
		      "Sparse to dense conversion failed!");
dev2file("H.bin",H,NI*NI*paCov.MT2_TOTAL);
HANDLE_ERROR( cudaFree( H ) );
*/

#line 480 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  }
#line 503 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  int N = N_SIDE_LENSLET_*N_SIDE_LENSLET_*2;
  if (strcmp(solver_id,"CG")==0)
    iSolve.cg_setup(N*N_guide_star);
  if (strcmp(solver_id,"MINRES")==0)
    iSolve.minres_setup(N*N_guide_star);

  iSolve.RTOL = 5E-2;

d__idx = 0;
d__ce  = 0;
int N_LENSLET = N_SIDE_LENSLET_*N_SIDE_LENSLET_;
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_c, sizeof(float)*PS_E_N_PX ) );
HANDLE_ERROR( cudaMalloc((void**)&d__x, sizeof(float)* N_LENSLET*2*N_guide_star ) );
HANDLE_ERROR( cudaMemset(d__x, 0, sizeof(float)* N_LENSLET*2*N_guide_star ) );
if (osf>1) {
   int NS2 = NS*NS;
   HANDLE_ERROR( cudaMalloc((void**)&d__zp_x, sizeof(float)* NS2*2*N_guide_star ) );
   HANDLE_ERROR( cudaMemset(d__zp_x, 0, sizeof(float)* NS2*2*N_guide_star ) );
}
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
#line 357 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}
#line 361 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::setup(atmosphere *atm, source *guide_star, source *mmse_star,
		       float sampling, int _N_SIDE_LENSLET_,
                       mask *pupil, char *solver_id)
{
  N_SIDE_LENSLET_ = _N_SIDE_LENSLET_;
osf = 1; // The wavefront over-sampling factor (1 means the same sampling that the WFS)
#line 451 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  N_guide_star = guide_star->N_SRC;
  N_mmse_star  = mmse_star->N_SRC;

  
#line 528 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aa.setup(N_SIDE_LENSLET_,atm,sampling,guide_star, N_guide_star);
#line 531 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aaCov.setup(2*N_guide_star,2*N_guide_star,N_SIDE_LENSLET_,aa.d__cov);

#line 456 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  NP = osf*N_SIDE_LENSLET_+1; // The number of wavefront sample
  NS = (osf>1) ? NP : N_SIDE_LENSLET_;

  pa.setup(NP,NS,osf,atm,sampling,
           mmse_star, N_mmse_star, guide_star, N_guide_star);

#line 483 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  if (mmse_star[0].height==guide_star[0].height) {

     paCov.setup(N_mmse_star,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
     paCov.mask = pupil->m;

     PS_E_N_PX = paCov.MT2_TOTAL;
     nnz = 0;

  } else {
 
    paCov.setup(N_mmse_star*atm->N_LAYER,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
  
     PS_E_N_PX = paCov.MT2_TOTAL;

     NI = NP;
     nnz = 4*atm->N_LAYER*pupil->nnz;
     
#line 838 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
alpha = 1;
beta  = 0;
HANDLE_ERROR_CUSPARSE( cusparseCreate(&handle), 
		       "CUSPARSE Library initialization failed!");
HANDLE_ERROR_CUSPARSE( cusparseCreateMatDescr(&descr), 
		       "Matrix descriptor initialization failed!");

cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr,CUSPARSE_INDEX_BASE_ZERO);

HANDLE_ERROR( cudaMalloc((void**)&csrValH,    sizeof(float)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrColIndH, sizeof(int)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrRowPtrH, sizeof(int)*(NI*NI+1) ) );
HANDLE_ERROR( cudaMemset(csrRowPtrH, 0, sizeof(int)*(NI*NI+1) ) );

HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_i, sizeof(float)*NI*NI ) );

#line 864 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearSparseOperatorMask <<< gridDim,blockDim >>> 
  (csrValH, csrColIndH, csrRowPtrH, paCov.d__MT, NI, 
   atm->d__layers, atm->N_LAYER, guide_star[0].height, pupil->m, pupil->nnz);

#line 872 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
/*
// For debugging purposes:

dev2file("csrValH.bin",csrValH,nnz);
dev2file("csrColIndH.bin",csrColIndH,nnz);
dev2file("csrRowPtrH.bin",csrRowPtrH,NI*NI+1);

float *H;
HANDLE_ERROR( cudaMalloc( (void**)&H, sizeof(float)*NI*NI*paCov.MT2_TOTAL ) );
HANDLE_ERROR_CUSPARSE(cusparseScsr2dense(handle, NI*NI, paCov.MT2_TOTAL, descr, 
					 csrValH, csrRowPtrH, csrColIndH,
					 H, NI*NI),
		      "Sparse to dense conversion failed!");
dev2file("H.bin",H,NI*NI*paCov.MT2_TOTAL);
HANDLE_ERROR( cudaFree( H ) );
*/


#line 500 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  }
#line 503 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  int N = N_SIDE_LENSLET_*N_SIDE_LENSLET_*2;
  if (strcmp(solver_id,"CG")==0)
    iSolve.cg_setup(N*N_guide_star);
  if (strcmp(solver_id,"MINRES")==0)
    iSolve.minres_setup(N*N_guide_star);

  iSolve.RTOL = 5E-2;

d__idx = 0;
d__ce  = 0;
int N_LENSLET = N_SIDE_LENSLET_*N_SIDE_LENSLET_;
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_c, sizeof(float)*PS_E_N_PX ) );
HANDLE_ERROR( cudaMalloc((void**)&d__x, sizeof(float)* N_LENSLET*2*N_guide_star ) );
HANDLE_ERROR( cudaMemset(d__x, 0, sizeof(float)* N_LENSLET*2*N_guide_star ) );
if (osf>1) {
   int NS2 = NS*NS;
   HANDLE_ERROR( cudaMalloc((void**)&d__zp_x, sizeof(float)* NS2*2*N_guide_star ) );
   HANDLE_ERROR( cudaMemset(d__zp_x, 0, sizeof(float)* NS2*2*N_guide_star ) );
}
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
#line 368 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}
#line 372 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::setup(atmosphere *atm, source *guide_star, source *mmse_star,
                  float sampling, int _N_SIDE_LENSLET_,
                  mask *pupil, char *solver_id, int wavefront_osf)
{
  N_SIDE_LENSLET_ = _N_SIDE_LENSLET_;
osf = wavefront_osf; // The wavefront over-sampling factor (1 means the same sampling that the WFS)
#line 451 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  N_guide_star = guide_star->N_SRC;
  N_mmse_star  = mmse_star->N_SRC;

  
#line 528 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aa.setup(N_SIDE_LENSLET_,atm,sampling,guide_star, N_guide_star);
#line 531 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aaCov.setup(2*N_guide_star,2*N_guide_star,N_SIDE_LENSLET_,aa.d__cov);

#line 456 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  NP = osf*N_SIDE_LENSLET_+1; // The number of wavefront sample
  NS = (osf>1) ? NP : N_SIDE_LENSLET_;

  pa.setup(NP,NS,osf,atm,sampling,
           mmse_star, N_mmse_star, guide_star, N_guide_star);

#line 483 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  if (mmse_star[0].height==guide_star[0].height) {

     paCov.setup(N_mmse_star,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
     paCov.mask = pupil->m;

     PS_E_N_PX = paCov.MT2_TOTAL;
     nnz = 0;

  } else {
 
    paCov.setup(N_mmse_star*atm->N_LAYER,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
  
     PS_E_N_PX = paCov.MT2_TOTAL;

     NI = NP;
     nnz = 4*atm->N_LAYER*pupil->nnz;
     
#line 838 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
alpha = 1;
beta  = 0;
HANDLE_ERROR_CUSPARSE( cusparseCreate(&handle), 
		       "CUSPARSE Library initialization failed!");
HANDLE_ERROR_CUSPARSE( cusparseCreateMatDescr(&descr), 
		       "Matrix descriptor initialization failed!");

cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr,CUSPARSE_INDEX_BASE_ZERO);

HANDLE_ERROR( cudaMalloc((void**)&csrValH,    sizeof(float)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrColIndH, sizeof(int)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrRowPtrH, sizeof(int)*(NI*NI+1) ) );
HANDLE_ERROR( cudaMemset(csrRowPtrH, 0, sizeof(int)*(NI*NI+1) ) );

HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_i, sizeof(float)*NI*NI ) );

#line 864 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearSparseOperatorMask <<< gridDim,blockDim >>> 
  (csrValH, csrColIndH, csrRowPtrH, paCov.d__MT, NI, 
   atm->d__layers, atm->N_LAYER, guide_star[0].height, pupil->m, pupil->nnz);

#line 872 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
/*
// For debugging purposes:

dev2file("csrValH.bin",csrValH,nnz);
dev2file("csrColIndH.bin",csrColIndH,nnz);
dev2file("csrRowPtrH.bin",csrRowPtrH,NI*NI+1);

float *H;
HANDLE_ERROR( cudaMalloc( (void**)&H, sizeof(float)*NI*NI*paCov.MT2_TOTAL ) );
HANDLE_ERROR_CUSPARSE(cusparseScsr2dense(handle, NI*NI, paCov.MT2_TOTAL, descr, 
					 csrValH, csrRowPtrH, csrColIndH,
					 H, NI*NI),
		      "Sparse to dense conversion failed!");
dev2file("H.bin",H,NI*NI*paCov.MT2_TOTAL);
HANDLE_ERROR( cudaFree( H ) );
*/


#line 500 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  }
#line 503 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  int N = N_SIDE_LENSLET_*N_SIDE_LENSLET_*2;
  if (strcmp(solver_id,"CG")==0)
    iSolve.cg_setup(N*N_guide_star);
  if (strcmp(solver_id,"MINRES")==0)
    iSolve.minres_setup(N*N_guide_star);

  iSolve.RTOL = 5E-2;

d__idx = 0;
d__ce  = 0;
int N_LENSLET = N_SIDE_LENSLET_*N_SIDE_LENSLET_;
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_c, sizeof(float)*PS_E_N_PX ) );
HANDLE_ERROR( cudaMalloc((void**)&d__x, sizeof(float)* N_LENSLET*2*N_guide_star ) );
HANDLE_ERROR( cudaMemset(d__x, 0, sizeof(float)* N_LENSLET*2*N_guide_star ) );
if (osf>1) {
   int NS2 = NS*NS;
   HANDLE_ERROR( cudaMalloc((void**)&d__zp_x, sizeof(float)* NS2*2*N_guide_star ) );
   HANDLE_ERROR( cudaMemset(d__zp_x, 0, sizeof(float)* NS2*2*N_guide_star ) );
}
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
#line 379 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}
#line 383 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::setup(atmosphere *atm, source *guide_star,
                  float sampling, int _N_SIDE_LENSLET_,
                  mask *pupil, char *solver_id, int wavefront_osf,
                  float z_radius)
{
  N_SIDE_LENSLET_ = _N_SIDE_LENSLET_;
osf = wavefront_osf; // The wavefront over-sampling factor (1 means the same sampling that the WFS)

  N_guide_star = guide_star->N_SRC;
  N_mmse_star  = 1;

  
#line 528 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aa.setup(N_SIDE_LENSLET_,atm,sampling,guide_star, N_guide_star);
#line 531 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aaCov.setup(2*N_guide_star,2*N_guide_star,N_SIDE_LENSLET_,aa.d__cov);

#line 396 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  NP = osf*N_SIDE_LENSLET_+1; // The number of wavefront sample
  NS = (osf>1) ? NP : N_SIDE_LENSLET_;

  pa.setup(NP,NS,osf,atm,sampling,
           guide_star, N_guide_star,
           z_radius);

  paCov.setup(N_mmse_star,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
  paCov.mask = pupil->m;

  PS_E_N_PX = paCov.MT2_TOTAL;
  nnz = 0;

#line 503 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  int N = N_SIDE_LENSLET_*N_SIDE_LENSLET_*2;
  if (strcmp(solver_id,"CG")==0)
    iSolve.cg_setup(N*N_guide_star);
  if (strcmp(solver_id,"MINRES")==0)
    iSolve.minres_setup(N*N_guide_star);

  iSolve.RTOL = 5E-2;

d__idx = 0;
d__ce  = 0;
int N_LENSLET = N_SIDE_LENSLET_*N_SIDE_LENSLET_;
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_c, sizeof(float)*PS_E_N_PX ) );
HANDLE_ERROR( cudaMalloc((void**)&d__x, sizeof(float)* N_LENSLET*2*N_guide_star ) );
HANDLE_ERROR( cudaMemset(d__x, 0, sizeof(float)* N_LENSLET*2*N_guide_star ) );
if (osf>1) {
   int NS2 = NS*NS;
   HANDLE_ERROR( cudaMalloc((void**)&d__zp_x, sizeof(float)* NS2*2*N_guide_star ) );
   HANDLE_ERROR( cudaMemset(d__zp_x, 0, sizeof(float)* NS2*2*N_guide_star ) );
}
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
#line 410 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}
#line 414 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::setup(atmosphere *atm, source *guide_star, source *mmse_star,
		       shackHartmann *wfs, char *solver_id)
{
float sampling      = wfs->lenslet_pitch;
N_SIDE_LENSLET_ = wfs->camera.N_SIDE_LENSLET;
mask *pupil;
pupil = &(wfs->valid_actuator);
osf = 1; // The wavefront over-sampling factor (1 means the same sampling that the WFS)
#line 451 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  N_guide_star = guide_star->N_SRC;
  N_mmse_star  = mmse_star->N_SRC;

  
#line 528 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aa.setup(N_SIDE_LENSLET_,atm,sampling,guide_star, N_guide_star);
#line 531 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aaCov.setup(2*N_guide_star,2*N_guide_star,N_SIDE_LENSLET_,aa.d__cov);

#line 456 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  NP = osf*N_SIDE_LENSLET_+1; // The number of wavefront sample
  NS = (osf>1) ? NP : N_SIDE_LENSLET_;

  pa.setup(NP,NS,osf,atm,sampling,
           mmse_star, N_mmse_star, guide_star, N_guide_star);

#line 483 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  if (mmse_star[0].height==guide_star[0].height) {

     paCov.setup(N_mmse_star,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
     paCov.mask = pupil->m;

     PS_E_N_PX = paCov.MT2_TOTAL;
     nnz = 0;

  } else {
 
    paCov.setup(N_mmse_star*atm->N_LAYER,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
  
     PS_E_N_PX = paCov.MT2_TOTAL;

     NI = NP;
     nnz = 4*atm->N_LAYER*pupil->nnz;
     
#line 838 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
alpha = 1;
beta  = 0;
HANDLE_ERROR_CUSPARSE( cusparseCreate(&handle), 
		       "CUSPARSE Library initialization failed!");
HANDLE_ERROR_CUSPARSE( cusparseCreateMatDescr(&descr), 
		       "Matrix descriptor initialization failed!");

cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr,CUSPARSE_INDEX_BASE_ZERO);

HANDLE_ERROR( cudaMalloc((void**)&csrValH,    sizeof(float)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrColIndH, sizeof(int)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrRowPtrH, sizeof(int)*(NI*NI+1) ) );
HANDLE_ERROR( cudaMemset(csrRowPtrH, 0, sizeof(int)*(NI*NI+1) ) );

HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_i, sizeof(float)*NI*NI ) );

#line 864 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearSparseOperatorMask <<< gridDim,blockDim >>> 
  (csrValH, csrColIndH, csrRowPtrH, paCov.d__MT, NI, 
   atm->d__layers, atm->N_LAYER, guide_star[0].height, pupil->m, pupil->nnz);

#line 872 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
/*
// For debugging purposes:

dev2file("csrValH.bin",csrValH,nnz);
dev2file("csrColIndH.bin",csrColIndH,nnz);
dev2file("csrRowPtrH.bin",csrRowPtrH,NI*NI+1);

float *H;
HANDLE_ERROR( cudaMalloc( (void**)&H, sizeof(float)*NI*NI*paCov.MT2_TOTAL ) );
HANDLE_ERROR_CUSPARSE(cusparseScsr2dense(handle, NI*NI, paCov.MT2_TOTAL, descr, 
					 csrValH, csrRowPtrH, csrColIndH,
					 H, NI*NI),
		      "Sparse to dense conversion failed!");
dev2file("H.bin",H,NI*NI*paCov.MT2_TOTAL);
HANDLE_ERROR( cudaFree( H ) );
*/


#line 500 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  }
#line 503 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  int N = N_SIDE_LENSLET_*N_SIDE_LENSLET_*2;
  if (strcmp(solver_id,"CG")==0)
    iSolve.cg_setup(N*N_guide_star);
  if (strcmp(solver_id,"MINRES")==0)
    iSolve.minres_setup(N*N_guide_star);

  iSolve.RTOL = 5E-2;

d__idx = 0;
d__ce  = 0;
int N_LENSLET = N_SIDE_LENSLET_*N_SIDE_LENSLET_;
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_c, sizeof(float)*PS_E_N_PX ) );
HANDLE_ERROR( cudaMalloc((void**)&d__x, sizeof(float)* N_LENSLET*2*N_guide_star ) );
HANDLE_ERROR( cudaMemset(d__x, 0, sizeof(float)* N_LENSLET*2*N_guide_star ) );
if (osf>1) {
   int NS2 = NS*NS;
   HANDLE_ERROR( cudaMalloc((void**)&d__zp_x, sizeof(float)* NS2*2*N_guide_star ) );
   HANDLE_ERROR( cudaMemset(d__zp_x, 0, sizeof(float)* NS2*2*N_guide_star ) );
}
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
#line 423 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}
#line 427 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::setup(atmosphere *atm, source *guide_star, source *mmse_star,
		       shackHartmann *wfs, char *solver_id, 
                       int osf_, mask *pupil_)
{
float sampling      = wfs->lenslet_pitch;
N_SIDE_LENSLET_ = wfs->camera.N_SIDE_LENSLET;
mask *pupil;
pupil = pupil_;
osf = osf_; // The wavefront over-sampling factor (1 means the same sampling that the WFS)
#line 451 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  N_guide_star = guide_star->N_SRC;
  N_mmse_star  = mmse_star->N_SRC;

  
#line 528 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aa.setup(N_SIDE_LENSLET_,atm,sampling,guide_star, N_guide_star);
#line 531 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
aaCov.setup(2*N_guide_star,2*N_guide_star,N_SIDE_LENSLET_,aa.d__cov);

#line 456 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  NP = osf*N_SIDE_LENSLET_+1; // The number of wavefront sample
  NS = (osf>1) ? NP : N_SIDE_LENSLET_;

  pa.setup(NP,NS,osf,atm,sampling,
           mmse_star, N_mmse_star, guide_star, N_guide_star);

#line 483 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  if (mmse_star[0].height==guide_star[0].height) {

     paCov.setup(N_mmse_star,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
     paCov.mask = pupil->m;

     PS_E_N_PX = paCov.MT2_TOTAL;
     nnz = 0;

  } else {
 
    paCov.setup(N_mmse_star*atm->N_LAYER,2*N_guide_star,pa.M_LAYER,NS,pa.d__cov);
  
     PS_E_N_PX = paCov.MT2_TOTAL;

     NI = NP;
     nnz = 4*atm->N_LAYER*pupil->nnz;
     
#line 838 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
alpha = 1;
beta  = 0;
HANDLE_ERROR_CUSPARSE( cusparseCreate(&handle), 
		       "CUSPARSE Library initialization failed!");
HANDLE_ERROR_CUSPARSE( cusparseCreateMatDescr(&descr), 
		       "Matrix descriptor initialization failed!");

cusparseSetMatType(descr,CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr,CUSPARSE_INDEX_BASE_ZERO);

HANDLE_ERROR( cudaMalloc((void**)&csrValH,    sizeof(float)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrColIndH, sizeof(int)*nnz ) );
HANDLE_ERROR( cudaMalloc((void**)&csrRowPtrH, sizeof(int)*(NI*NI+1) ) );
HANDLE_ERROR( cudaMemset(csrRowPtrH, 0, sizeof(int)*(NI*NI+1) ) );

HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_i, sizeof(float)*NI*NI ) );

#line 864 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearSparseOperatorMask <<< gridDim,blockDim >>> 
  (csrValH, csrColIndH, csrRowPtrH, paCov.d__MT, NI, 
   atm->d__layers, atm->N_LAYER, guide_star[0].height, pupil->m, pupil->nnz);

#line 872 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
/*
// For debugging purposes:

dev2file("csrValH.bin",csrValH,nnz);
dev2file("csrColIndH.bin",csrColIndH,nnz);
dev2file("csrRowPtrH.bin",csrRowPtrH,NI*NI+1);

float *H;
HANDLE_ERROR( cudaMalloc( (void**)&H, sizeof(float)*NI*NI*paCov.MT2_TOTAL ) );
HANDLE_ERROR_CUSPARSE(cusparseScsr2dense(handle, NI*NI, paCov.MT2_TOTAL, descr, 
					 csrValH, csrRowPtrH, csrColIndH,
					 H, NI*NI),
		      "Sparse to dense conversion failed!");
dev2file("H.bin",H,NI*NI*paCov.MT2_TOTAL);
HANDLE_ERROR( cudaFree( H ) );
*/


#line 500 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  }
#line 503 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  int N = N_SIDE_LENSLET_*N_SIDE_LENSLET_*2;
  if (strcmp(solver_id,"CG")==0)
    iSolve.cg_setup(N*N_guide_star);
  if (strcmp(solver_id,"MINRES")==0)
    iSolve.minres_setup(N*N_guide_star);

  iSolve.RTOL = 5E-2;

d__idx = 0;
d__ce  = 0;
int N_LENSLET = N_SIDE_LENSLET_*N_SIDE_LENSLET_;
HANDLE_ERROR( cudaMalloc( (void**)&d__phase_est_c, sizeof(float)*PS_E_N_PX ) );
HANDLE_ERROR( cudaMalloc((void**)&d__x, sizeof(float)* N_LENSLET*2*N_guide_star ) );
HANDLE_ERROR( cudaMemset(d__x, 0, sizeof(float)* N_LENSLET*2*N_guide_star ) );
if (osf>1) {
   int NS2 = NS*NS;
   HANDLE_ERROR( cudaMalloc((void**)&d__zp_x, sizeof(float)* NS2*2*N_guide_star ) );
   HANDLE_ERROR( cudaMemset(d__zp_x, 0, sizeof(float)* NS2*2*N_guide_star ) );
}
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
#line 437 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}

#line 536 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::cleanup(void)
{
  INFO("@(CEO)>LMMSE: freeing memory!\n");

  INFO(" |-");
  aa.cleanup();
  INFO(" |-");
  pa.cleanup();
  INFO(" |-");
  aaCov.cleanup();
  INFO(" |-");
  paCov.cleanup();
  INFO(" |-");
  iSolve.cleanup();

  HANDLE_ERROR( cudaEventDestroy( start ) );
  HANDLE_ERROR( cudaEventDestroy( stop ) );

  if (d__idx) HANDLE_ERROR( cudaFree( d__idx ) );
  if (d__ce)  HANDLE_ERROR( cudaFree( d__ce ) );
  HANDLE_ERROR( cudaFree( d__phase_est_c ) );
  HANDLE_ERROR( cudaFree(d__x) );
  if (osf>1) {
     HANDLE_ERROR( cudaFree(d__zp_x) );
  }
  if (nnz) {
     
#line 891 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
HANDLE_ERROR_CUSPARSE( cusparseDestroyMatDescr(descr),
		       "Matrix descriptor destruction failed!");
HANDLE_ERROR_CUSPARSE( cusparseDestroy(handle),
		       "CUSPARSE Library release of resources failed!");
HANDLE_ERROR( cudaFree( csrValH ) );
HANDLE_ERROR( cudaFree( csrColIndH ) );
HANDLE_ERROR( cudaFree( csrRowPtrH ) );
HANDLE_ERROR( cudaFree( d__phase_est_i ) );
#line 563 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
  }

}

#line 571 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::estimation(const centroiding *cog)
{
  
#line 586 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
HANDLE_ERROR( cudaEventRecord( start, 0 ) );

#line 615 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
if (cog->MASK_SET) aaCov.mask = cog->lenslet_mask;

//tid.tic();

iSolve.minres_vorst(d__x, &aaCov, cog->d__c, d__x);

//tid.toc("WAVEFRONT ESTIMATION (MINRES)");

//tid.tic();

if (osf>1) {
  dim3 blockDim(16,16);
  dim3 gridDim(N_SIDE_LENSLET_/16+1,N_SIDE_LENSLET_/16+1,N_guide_star);
  zeroPadding <<< gridDim,blockDim >>>(d__zp_x, NP, d__x, N_SIDE_LENSLET_, osf);
  paCov.MVM(d__phase_est_c,d__zp_x);
} else {
  paCov.MVM(d__phase_est_c,d__x);
}

d__phase_est = d__phase_est_c;

//tid.toc("WAVEFRONT ESTIMATION (PA MVM)");
//tid.toc("WAVEFRONT ESTIMATION");

#line 590 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
#if CUDA_VERSION < 11000
if (nnz) {
//tid.tic();

HANDLE_ERROR_CUSPARSE( cusparseScsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
				      NI*NI, paCov.MT2_TOTAL, nnz, &alpha,
				      descr, csrValH, csrRowPtrH, csrColIndH,
				      d__phase_est_c, &beta, d__phase_est_i),
		       "Sparse to vector product failed!");

d__phase_est = d__phase_est_i;

//tid.toc("WAVEFRONT ESTIMATION (SPARSE MVM)");
}
#endif

//printf("Solver residue norm           : %.2E\n",iSolve.rnorm);
//printf("\nSolver mean time per iteration: %.2E\n",iSolve.mean_time_per_iteration);

HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
HANDLE_ERROR( cudaEventSynchronize( stop ) );
HANDLE_ERROR( cudaEventElapsedTime( &elapsed_time, start, stop ) );


#line 574 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}
#line 578 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::estimation(shackHartmann *wfs)
{
  centroiding *cog;
  cog = &(wfs->data_proc);
  
#line 586 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
HANDLE_ERROR( cudaEventRecord( start, 0 ) );

#line 615 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
if (cog->MASK_SET) aaCov.mask = cog->lenslet_mask;

//tid.tic();

iSolve.minres_vorst(d__x, &aaCov, cog->d__c, d__x);

//tid.toc("WAVEFRONT ESTIMATION (MINRES)");

//tid.tic();

if (osf>1) {
  dim3 blockDim(16,16);
  dim3 gridDim(N_SIDE_LENSLET_/16+1,N_SIDE_LENSLET_/16+1,N_guide_star);
  zeroPadding <<< gridDim,blockDim >>>(d__zp_x, NP, d__x, N_SIDE_LENSLET_, osf);
  paCov.MVM(d__phase_est_c,d__zp_x);
} else {
  paCov.MVM(d__phase_est_c,d__x);
}

d__phase_est = d__phase_est_c;

//tid.toc("WAVEFRONT ESTIMATION (PA MVM)");
//tid.toc("WAVEFRONT ESTIMATION");

#line 590 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
#if CUDA_VERSION < 11000
if (nnz) {
//tid.tic();

HANDLE_ERROR_CUSPARSE( cusparseScsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
				      NI*NI, paCov.MT2_TOTAL, nnz, &alpha,
				      descr, csrValH, csrRowPtrH, csrColIndH,
				      d__phase_est_c, &beta, d__phase_est_i),
		       "Sparse to vector product failed!");

d__phase_est = d__phase_est_i;

//tid.toc("WAVEFRONT ESTIMATION (SPARSE MVM)");
}
#endif

//printf("Solver residue norm           : %.2E\n",iSolve.rnorm);
//printf("\nSolver mean time per iteration: %.2E\n",iSolve.mean_time_per_iteration);

HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
HANDLE_ERROR( cudaEventSynchronize( stop ) );
HANDLE_ERROR( cudaEventElapsedTime( &elapsed_time, start, stop ) );


#line 583 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
}

#line 665 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::reset(void)
{
  HANDLE_ERROR( cudaMemset(d__x, 0, 
                           sizeof(float)*N_SIDE_LENSLET_*N_SIDE_LENSLET_*2*N_guide_star ) );
}

#line 1092 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::set_phase_est_ptr(float *data_ptr)
{
if (nnz)
  d__phase_est_i = data_ptr;
else
  d__phase_est_c = data_ptr;
}

#line 1103 "/home/ubuntu/projects/crseo/sys/CEO/LMMSE/LMMSE.nw"
void LMMSE::toFile(const char *filename){
  dev2file(filename,d__phase_est,NP*NP);
}