%-*- mode: Noweb; noweb-code-mode: c-mode -*-

\index{utilities}
\section{The files}

\subsection{Header}

<<utilities.h>>=
#ifndef __UTILITIES_H__
#define __UTILITIES_H__

#include <sys/stat.h>
#include <stdio.h>
#include <cuda.h>
#include <cuda_runtime.h>
#include <curand.h>
#include <curand_kernel.h>
#include <cufft.h>
#include "cublas_v2.h"
#include "cusparse_v2.h"
#include <cusolverDn.h>
  
#define PI 3.141592653589793
#define RADIAN2ARCSEC (180*3600/PI)
#define MIN(a, b)  (((a) < (b)) ? (a) : (b))
#define MAX(a, b)  (((a) > (b)) ? (a) : (b))
#define CUBLAS_VERBOSE 0
#define RAND_SEED 2013
#define N_THREAD 16
#define N_THREAD2 256
#define SQRT_DBL_MAX   1.3407807929942596e+154
  //#define M_PI       3.14159265358979323846264338328      /* pi */
#define DBL_EPSILON        2.2204460492503131e-16

static void Error( const char *errMsg,
		   const char *file,
		   int line ) {
  fprintf(stderr,"\n\x1B[31m@(CEO)>ERROR: %s in %s at line %d\x1B[0m\n", 
	  errMsg, file, line );
  exit( EXIT_FAILURE );
}
#define ERROR( errMsg ) (Error( errMsg, __FILE__, __LINE__ ))

#ifdef SILENT
#define INFO( infoMsg, args... );
#else
#define INFO( infoMsg , args...)   fprintf(stdout,infoMsg, ##args);
#endif

static void HandleError( cudaError_t err,
                         const char *file,
                         int line ) {
  if (err != cudaSuccess)  {
    fprintf(stderr,"\n\x1B[31m@(CEO)>ERROR: %s in %s at line %d\x1B[0m\n", 
	    cudaGetErrorString( err ), file, line );
    exit( EXIT_FAILURE );
  }
}
#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))

static void HandleErrorCUFFT( cufftResult_t err,
			      const char *msg,
			      const char *file,
			      int line ) {
  if (err != CUFFT_SUCCESS)  {
    fprintf(stderr,"\n\x1B[31m@(CEO)>CUFFT ERROR: %s in %s at line %d\x1B[0m\n", 
	    msg, file, line );
    exit( EXIT_FAILURE );
  }
}
#define HANDLE_ERROR_CUFFT( err, msg ) (HandleErrorCUFFT( err, msg, __FILE__, __LINE__ ))

static void HandleErrorCUSPARSE( cusparseStatus_t err,
			      const char *msg,
			      const char *file,
			      int line ) {
  if (err != CUSPARSE_STATUS_SUCCESS)  {
    fprintf(stderr,"\n\x1B[31m@(CEO)>CUSPARSE ERROR: %s in %s at line %d\x1B[0m\n", 
	    msg, file, line );
    exit( EXIT_FAILURE );
  }
}
#define HANDLE_ERROR_CUSPARSE( err, msg ) (HandleErrorCUSPARSE( err, msg, __FILE__, __LINE__ ))   


static void HandleErrorCURAND( curandStatus_t err,
			      const char *msg,
			      const char *file,
			      int line ) {
  if (err != CURAND_STATUS_SUCCESS)  {
    fprintf(stderr,"\n\x1B[31m@(CEO)>CURAND ERROR: %s in %s at line %d\x1B[0m\n", 
	    msg, file, line );
    exit( EXIT_FAILURE );
  }
}
#define HANDLE_ERROR_CUSOLVER( err ) (HandleErrorCUSOLVER( err, __FILE__, __LINE__ ))   

static void HandleErrorCUSOLVER( cusolverStatus_t err,
                               const char *file,
                               int line ) {
  if (err != CUSOLVER_STATUS_SUCCESS)  {
    if (err==CUSOLVER_STATUS_NOT_INITIALIZED)
      fprintf(stderr,"\n\x1B[31m@(CEO)>CUSOLVER ERROR: %s in %s at line %d:.\x1B[0m\n", 
              "the library was not initialized", file, line );
    else if (err==CUSOLVER_STATUS_INVALID_VALUE)
      fprintf(stderr,"\n\x1B[31m@(CEO)>CUSOLVER ERROR: %s in %s at line %d\x1B[0m\n", 
              "invalid parameters were passed (m,n<0 or lda<max(1,m))", file, line );
    else if (err==CUSOLVER_STATUS_ARCH_MISMATCH)
      fprintf(stderr,"\n\x1B[31m@(CEO)>CUSOLVER ERROR: %s in %s at line %d\x1B[0m\n", 
              "the device only supports compute capability 2.0 and above", file, line );
    else
      fprintf(stderr,"\n\x1B[31m@(CEO)>CURAND ERROR: %s in %s at line %d\x1B[0m\n",
              "an internal operation failed.", file, line );
    exit( EXIT_FAILURE );
  }
}
#define HANDLE_ERROR_CURAND( err, msg ) (HandleErrorCURAND( err, msg, __FILE__, __LINE__ ))   


inline float sign(float x) {
  return (float) ( (x > 0) - (x < 0) );
}

int round_up_to_nhp2(int n);
void set_device(int id);
int get_device_count();

#define ARCSEC(a) (a*PI/180/3600)

<<ray tracing data type>>

<<vector data type>>

<<lenslet2array>>
<<lenslet2array (overlap)>>

<<lenslet2row>>

<<threads2lenslet>>

<<tic-toc>>

void CUBLAS_ERROR(cublasStatus_t status);

<<mask structure>>
<<statistics definition>>

<<lenslet 2 actuator>>
__global__ void fill_ones_char(char *ones, int n_data) ;
__global__ void circular_pupil(char *pupil, int N, float scale);
__global__ void square_pupil(char *pupil, int N, int N_S, int I_S, int J_S);
__global__ void rot_square_pupil(char *pupil, int N, float N_S, 
				 float I_0, float J_0, float theta, float I_S, float J_S);
__global__ void GMT_pupil(char *pupil, int *piston_mask , const int N, 
			  const float S, float Dsh, const float c, const char v);
//__global__ void valid_lenslet(char *mask, int N_pupil, float threshold);
__global__ void fried_geometry(char *lenslet_mask, char *actuator_mask, 
                               int NL, int n, float threshold);
__global__ void fried_geometry_vs_pupil(char *lenslet_mask, char *actuator_mask, 
                                        const int NL, const int n, 
					const float threshold, const char *pupil);
__global__ void wavefront_finite_difference_kernel(float *sx, float *sy, int NL, 
                                                   float *phi, int n, float d);
__global__ void wavefront_finite_difference_kernel_masked(
                                 float *sx, float *sy, int NL, 
                                 float *phi, int n, float d,
                                 char *M);
__global__ void wavefront_finite_difference_sparse_matrix_kernel(
                                                   float *csrValH, 
                                                   int *csrColIndH, 
                                                   int *csrRowPtrH, 
                                                   const int NL, const int n, const float d);
__global__ void set_gmt_piston_kernel(float *phase, float *p, 
				      char *m, int *piston_mask,
				      const int N);
 __global__ void double2float_kern(float *d__single_data, double *d__double_data, const int N);

char fried_geometry_setup(char *lenslet_mask, char *actuator_mask, 
                          int NL, int n, float threshold);
char fried_geometry_setup_vs_pupil(char *lenslet_mask, char *actuator_mask, 
			  int NL, int n, float threshold, char *pupil);

void dev2file(const char* filename, float* d__data, int n_data);
void dev2file(const char* filename, float* d__data, int n_data, 
              int n_data_page, int n_page);
void dev2file(const char* filename, float* d__data, int n_data, 
              int *n_data_page, int n_page);
void dev2file(const char* filename, int* d__data, int n_data);
void dev2file(const char* filename, float2* d__data, int n_data);
void dev2file(const char* filename, char* d__data, int n_data);

double atmosphere_refractive_index(float wavelength, float altitude, float temperature, float humidity);
double atmospheric_dispersion(float wavelength, float delta_wavelength,
			      float zenith, float altitude, 
			      float temperature, float humidity);

void dev2host( float *host_data, float *dev_data, int N);
void dev2host_int( int *host_data, int *dev_data, int N);
void host2dev( float *dev_data, float *host_data, int N);
void host2dev_char( char *dev_data, char *host_data, int N);
void freedev( float **dev_data );

<<template gpu data structure>>

void geqrf(float *d_tau, float *a, int m, int n);
void ormqr(float *b, int m, float *q, float *tau, int n);

#if CUDA_VERSION < 11000
<<sparse matrix>>
#endif

void wavefront_finite_difference(float *sx, float *sy, int NL, 
                                 float *phi, int n, float d, 
				 int N_GS);
void wavefront_finite_difference(float *sx, float *sy, int NL, 
                                 float *phi, int n, float d, 
                                 mask *valid_lenset, int N_GS);

struct cheb_series_struct {
  double * c;   /* coefficients                */
  int order;    /* order of expansion          */
  double a;     /* lower interval point        */
  double b;     /* upper interval point        */
  int order_sp; /* effective single precision order */
};
typedef struct cheb_series_struct cheb_series;

__host__ __device__ void
  cheb_eval_e(const cheb_series * cs,
              const double x,
              double * result);
__host__ __device__ void
temme_gamma(const double nu, double * g_1pnu, double * g_1mnu, double * g1, double * g2);
__host__ __device__ void
K_scaled_temme(const double nu, const double x,
               double * K_nu, double * K_nup1, double * Kp_nu);
__host__ __device__ void
K_scaled_steed_temme_CF2(const double nu, const double x,
                         double * K_nu, double * K_nup1,
                         double * Kp_nu);
__host__ __device__ double _K_nu_(const double nu, const double x);

__host__ __device__ int polywind(double Px, double Py, double *Vx, double *Vy, int NV);
void polywinds(int *W, double *Px, double *Py, int NP, double *Vx, double *Vy, int NV);

#endif  // __UTILITIES_H__
@

\subsection{Source}

<<utilities.cu>>=
#ifndef __UTILITIES_H__
#include "utilities.h"
#endif

int round_up_to_nhp2(int n) {
  int v; // compute the next highest power of 2 of 32-bit v
  v = n;
  v--;
  v |= v >> 1;
  v |= v >> 2;
  v |= v >> 4;
  v |= v >> 8;
  v |= v >> 16;
  v++;
  return v;
}
void set_device(int id) {
  cudaSetDevice(id);
}
int get_device_count() {
  int n;
  cudaGetDeviceCount(&n);
  return n;
}

<<cublas error>>

<<vector functions>>

<<double2float kernel>>

<<fill ones>>
<<fill ones (char)>>
<<fill ones filtered>>
<<fill ones altered>>
<<add mask kernel>>
<<non zeros index>>
<<statistics functions>>
<<mask functions>>
<<circular pupil>>
<<square pupil>>
<<GMT pupil>>
<<GMT segment piston kernel>>
<<polygon winding number kernel>>

<<GPU float methods>>
<<GPU double methods>>
<<QR factorization>>
<<compute Q^T*B>>

<<Fried geometry masks>>
<<Fried geometry>>
<<Fried geometry masks (arbitrary pupil)>>
<<Fried geometry (arbitrary pupil)>>

<<device to file>>

<<atmosphere refractive index>>
<<atmospheric dispersion>>

<<data to host from device>>
<<data to device from host>>
<<free device>>

<<wavefront differentiation>>
<<wavefront differentiation (with mask)>>
<<wavefront differentiation kernel>>
<<wavefront differentiation kernel (with mask)>>
<<wavefront differentiation sparse matrix kernel>>
<<bilinear interpolation sparse matrix kernel>>
<<bilinear interpolation sparse matrix kernel with pupil mask>>

#if CUDA_VERSION < 11000
<<sparse matrix methods>>
#endif

<<Knu temme>>
<<Knu>>
<<polygon winding number>>
@
\section{Functions}
\label{sec:functions}


\subsection{Custom types}
\label{sec:custom-types}

The accuracy of the ray tracing is partly depending upon the data type used to hold data values of rays and surfaces.
GPU are faster in single precision than in double precision but double precision may still be required for accuracy.
So lets defined a ray tracing type:
\index{utilities!vector!=}
\index{utilities!vector!+}
\index{utilities!vector!+=}
\index{utilities!vector!-}
\index{utilities!vector!-=}
\index{utilities!vector!*}
<<ray tracing data type>>=
typedef double rtd;
@ 

\subsection{Vector}
\label{sec:vector}

\index{utilities!vector}

A new vector data type is defined 
<<vector data type>>=
struct vector{
  rtd x;
  rtd y;
  rtd z;
  __host__ __device__ vector(rtd x=0.0, rtd y=0.0, rtd z=0.0)
    : x(x), y(y), z(z)
  {
  }
  __host__ __device__ rtd rho2(void);
  __host__ __device__ rtd rho2_shift(rtd x0, rtd y0);
  __host__ __device__ rtd mag(void);
  __host__ __device__ rtd mag(rtd R);
  __host__ __device__ rtd angle(void);
  __host__ __device__ rtd norm(void);
  __host__ __device__ rtd dot(vector *u);
  __host__ __device__ void unit(void);
  __host__ __device__ void left_cross(vector *w, vector *u);
  __host__ __device__ void right_cross(vector *w, vector *u);
  __host__ __device__ vector& operator=(const vector& v){
     x = v.x;
     y = v.y;
     z = v.z;
     return *this;
   }
  __host__ __device__ vector& operator-=(const vector& rhs){
     x -= rhs.x;
     y -= rhs.y;
     z -= rhs.z;
     return *this;
   }
  __host__ __device__ vector& operator+=(const vector& rhs){
     x += rhs.x;
     y += rhs.y;
     z += rhs.z;
     return *this;
   }
  __host__ __device__ vector operator-(const vector& rhs) const
  {
    return vector(x-rhs.x,y-rhs.y,z-rhs.z);
  }
  __host__ __device__ vector operator+(const vector& rhs) const
  {
    return vector(x+rhs.x,y+rhs.y,z+rhs.z);
  }
  __host__ __device__ rtd operator*(const vector& rhs) const
  {
    return x*rhs.x + y*rhs.y + z*rhs.z;
  }
  __host__ __device__ vector operator*(const rtd& rhs) const
  {
    return vector(x*rhs,y*rhs,z*rhs);
  }
};
@ that represents the 3 coordinates in space of the tip of vector.
The square of the magnitude of a vector projected in the x,y--plane i.e. $[[x]]^2+[[y]]^2$ is computed with 
\index{utilities!vector!rho2}
<<vector functions>>=
__host__ __device__ rtd vector::rho2(void) {
  return x*x + y*y;
}
@
The square of the magnitude of a vector projected in the x,y--plane with respect to origin ([[x0]],[[y0]]) i.e. $[[x-x0]]^2+[[y-y0]]^2$ is computed with 
\index{utilities!vector!rho2\_shift}
<<vector functions>>=
__host__ __device__ rtd vector::rho2_shift(rtd x0, rtd y0) {
  rtd xp, yp;
  xp = x - x0;
  yp = y - y0;
  return xp*xp + yp*yp;
}
@ 
The magnitude and angle of the projection of the vector in the x,y--plane is
\index{utilities!vector!mag}
\index{utilities!vector!angle}
<<vector functions>>=
__host__ __device__ rtd vector::mag(void) {
  return hypot(x,y);
}
__host__ __device__ rtd vector::angle(void) {
  return atan2(y,x);
}
@  
\index{utilities!vector!dot}
The dot product of this vector $\mathbf v$ with a vector $\mathbf u$ i.e $\mathbf{v}\cdot\mathbf{u}$ is
<<vector functions>>=
__host__ __device__ rtd vector::dot(vector *u) {
  return x*u->x + y*u->y + z*u->z;
}
@ 
The norm of the vector is
\index{utilities!vector!norm}
<<vector functions>>=
__host__ __device__ rtd vector::norm(void) {
  return sqrt( x*x + y*y + z*z );
}
@ 
The vector is normalized such as $\left\|\mathbf{v}\right\|=1$ with
\index{utilities!vector!unit}
<<vector functions>>=
__host__ __device__ void vector::unit(void) {
  rtd n = norm();
  x /= n;
  y /= n;
  z /= n;
}

@ 
The magnitude can also be normalized by radius [[R]]:
\index{utilities!vector!mag}
<<vector functions>>=
__host__ __device__ rtd vector::mag(rtd R) {
  return hypot(x,y)/R;
}
@ 
The cross product of a vector $\mathbf{u}$ with this vector $\mathbf{v}$ i.e. $\mathbf{w}=\mathbf{u}\wedge\mathbf{v}$ 
\index{utilities!vector!left\_cross}
<<vector functions>>=
__host__ __device__ void vector::left_cross(vector *w, vector *u) {
  w->x = u->y*z - u->z*y;
  w->y = u->z*x - u->x*z;
  w->z = u->x*y - u->y*x;
}
@ 
The cross product of this vector $\mathbf{v}$ with a vector $\mathbf{u}$ i.e. $\mathbf{w}=\mathbf{v}\wedge\mathbf{u}$ 
\index{utilities!vector!right\_cross}
<<vector functions>>=
__host__ __device__ void vector::right_cross(vector *w, vector *u) {
  w->x = u->z*y - u->y*z;
  w->y = u->x*z - u->z*x;
  w->z = u->y*x - u->x*y;
}
@
\subsection{Mapping to lenslet array}
\label{sec:mapp-lensl-array}

When dealing with lenslet arrays, a thread in a given thread block correspond to a pixel coordinate $[i,j]$ in a given lenslet $[i_l,j_L]$.
The following function computes the coordinates $[i,j]$ and $[i_l,j_L]$ based on the thread and block index, the block size and the number of pixel per lenslet [[n_px_lenslet]].
\index{utilities!threads2lenslet}
<<threads2lenslet>>=
__device__ inline char threads2lenslet(dim3 threadIdx, dim3 blockIdx, dim3 blockDim, 
				      int *i_px_lenslet, int *j_px_lenslet, int n_px_lenslet, 
				       int *i_lenslet, int *j_lenslet, int n_lenslet) 
{
  <<from threads to lenslet coordinates>>
  return 1;
}
@ The pixel coordinate in the image is given by
<<from threads to lenslet coordinates>>=
*i_px_lenslet = threadIdx.x + blockIdx.x * blockDim.x;
*j_px_lenslet = threadIdx.y + blockIdx.y * blockDim.y;
@ Then the lenslet coordinate $[i_l,j_L]$ are derived from
<<from threads to lenslet coordinates>>=
*i_lenslet = *i_px_lenslet/n_px_lenslet;
if (*i_lenslet>=n_lenslet)
  return 0;
*j_lenslet = *j_px_lenslet/n_px_lenslet;
if (*j_lenslet>=n_lenslet)
  return 0;
@ and finally the pixel coordinate $[i,j]$ in lenslet $[i_l,j_L]$ is
<<from threads to lenslet coordinates>>=
*i_px_lenslet -=  *i_lenslet*n_px_lenslet;
*j_px_lenslet -=  *j_lenslet*n_px_lenslet;
@ 
Fig.~\ref{fig:1} shows the arrangements for $ 2\times2\times3$ blocks of $16\times16$ threads mapping a $6\times6$ lenslet array with $5\times5$ pixels per lenslet.
The third dimension of the block grid gives the source index.
\begin{figure}
  \centering
  \begin{tikzpicture}
    \foreach \z in {0,4,8} { 
    \begin{scope}[xshift=-\z mm,yshift=-\z mm]
      \fill[orange!50] (-2mm,-2mm) rectangle (34mm,34mm);
      \foreach \x in {0,16} { 
        \foreach \y in {0,16} {
          \begin{scope}[xshift=\x mm,yshift=\y mm]
            \fill[green!50,draw=black,thick] (0,0) rectangle (16mm,16mm);
            \draw[step=1mm,gray,thin] (0,0) grid (16mm,16mm);
          \end{scope}
        }
      }
      \draw[step=5mm,red,very thin] (0,0) grid (30mm,30mm);
    \end{scope}
}
  \end{tikzpicture}
  \caption{$2\times2\times3$ blocks of $16\times16$ threads mapping a $6\times6$ lenslet array with $5\times5$ pixels per lenslet}
  \label{fig:1}
\end{figure}
@
Once the coordinates $[i,j]$ and $[i_l,j_L]$ have been computed, the linear index in the image is computed with
\index{utilities!lenslet2array}
<<lenslet2array>>=
__device__ inline int lenslet2array(int i_px_lenslet, int j_px_lenslet, int n_px_lenslet, 
			 int i_lenslet, int j_lenslet, int n_lenslet, int i_source)
{
  int index;
  index =  i_lenslet*n_px_lenslet + i_px_lenslet;
  index += i_source*n_px_lenslet*n_lenslet;
  index *= n_lenslet*n_px_lenslet;
  index += j_lenslet*n_px_lenslet + j_px_lenslet;
  return index;
}
@ 
\index{utilities!lenslet2array\_overlap}
<<lenslet2array (overlap)>>=
__device__ inline int lenslet2array_overlap(int i_px_lenslet, int j_px_lenslet, int n_px_lenslet, 
			 int i_lenslet, int j_lenslet, int n_lenslet, int i_source)
{
  int index;
  index =  i_lenslet*(n_px_lenslet-1) + i_px_lenslet;
  index += i_source*((n_px_lenslet-1)*n_lenslet+1);
  index *= n_lenslet*(n_px_lenslet-1)+1;
  index += j_lenslet*(n_px_lenslet-1) + j_px_lenslet;
  return index;
}
@ Fig.~\ref{fig:2} shows the global coordinates of pixel $[i_g,j_g]$ in an image made of the horizontal concatenation of 3 lenslet arrays.
$i_g$ and $j_g$ are related to the coordinates $[i,j]$ and $[i_l,j_L]$ through the following relation
\begin{eqnarray}
  \label{eq:1}
  i_g &=& i_SN_LN_P + i_LN_p+i \\
  j_g &=& j_LN_p+i \\
\end{eqnarray}
where $i_S$ is the source index, $N_L$ is the size of one lenslet array and $N_P$ is the number of pixel per lenslet.
The linear global index $k_g$ is given by
\begin{equation}
  \label{eq:2}
  k_g = i_gN_LN_P + j_g
\end{equation}
where $N_s$ is the number of sources.
\begin{figure}
  \centering
  \begin{tikzpicture}[x=1mm,y=1mm]
    \foreach \z in {0,30,60} { 
      \begin{scope}[xshift=\z mm]
        \fill[orange!50,draw=black,thick] (0,0) rectangle (30,30);
      \end{scope}
    }
    \draw[step=1mm,gray,thin] (0,0) grid (90,30);
    \draw[step=5mm,red,very thin] (0,0) grid (90,30);
    \fill[blue!70] (36,16) rectangle (37,17);
    \draw[blue!70] (-2,16.5) -- (36.5,16.5);
    \draw[blue!70] (36.5,-2) -- (36.5,16.5);
    \node[anchor=east] at (-2,16.5) {$i_SN_LN_P + i_LN_p+i$};
    \node[anchor=north] at (36.5,-2) {$j_LN_p+j$};
  \end{tikzpicture}
  \caption{3 $6\times6$ lenslet arrays with $5\times5$ pixels per lenslet}
  \label{fig:2}
\end{figure}

To compute the Fourier transform of the wavefront of each lenslet in a batch, all the lenslets are concatenated into a single column:
\index{utilities!lenslet2row}
<<lenslet2row>>=
__device__ inline int lenslet2row(int i_px_lenslet, int j_px_lenslet, int n_px_lenslet, 
			 int i_lenslet, int j_lenslet, int n_lenslet, 
			 int i_source)
{
  int index;
  index = i_lenslet*n_lenslet + j_lenslet;
  index += i_source*n_lenslet*n_lenslet;
  index *= n_px_lenslet;
  index += i_px_lenslet;
  index *= n_px_lenslet;
  index += j_px_lenslet;
  return index;
}
@ Fig.~\ref{fig:3} shows the global coordinates of pixel $[i_g,j_g]$ in an image made of the vertical concatenation of all the lenslets.
$i_g$ and $j_g$ are related to the coordinates $[i,j]$ and $[i_l,j_L]$ through the following relation
\begin{eqnarray}
  \label{eq:3}
  i_g &=& (i_LN_L+j_L + i_SN_L^2)N_P + i \\
  j_g &=& j \\
\end{eqnarray}
where $i_S$ is the source index, $N_L$ is the size of one lenslet array and $N_P$ is the number of pixel per lenslet.
The linear global index $k_g$ is given by
\begin{equation}
  \label{eq:4}
  k_g = i_gN_P + j_g
\end{equation}
where $N_s$ is the number of sources.
\begin{figure}
  \centering
  \begin{tikzpicture}[x=1mm,y=1mm]
    \fill[orange!50,draw=black,thick] (0,0) rectangle (5,33);
    \draw[step=1mm,gray,thin] (0,0) grid (5,33);
    \draw[step=5mm,red,very thin] (0,0) grid (5,33);
    \fill[blue!70] (3,16) rectangle (4,17);
    \draw[blue!70] (-2,16.5) -- (3.5,16.5);
    \draw[blue!70] (3.5,-2) -- (3.5,16.5);
    \node[anchor=east] at (-2,16.5) {$(i_LN_L+j_L + i_SN_L^2)N_P + i$};
    \node[anchor=north] at (3.5,-2) {$j$};
  \end{tikzpicture}
  \caption{3 $6\times6$ lenslet arrays with $5\times5$ pixels per lenslet}
  \label{fig:3}
\end{figure}

\subsection{Stopwatch}
\label{sec:stopwatch}

In the following, a structure to time a chunk of CUDA code is defined:
\index{utilities!stopwatch}
<<tic-toc>>=
#ifdef SILENT
struct stopwatch{
  float   elapsedTime;
  void tic(void) { elapsedTime=0.0; }
  void toc(void) {}
  void toc(const char *message) {}
  void toc(float *elapsedTime, const char *message) {}
  void toc(float *elapsedTime) {}
};
#else
struct stopwatch{
  float   elapsedTime;
  cudaEvent_t     start, stop;
  void tic(void) {
    HANDLE_ERROR( cudaEventCreate( &start ) );
    HANDLE_ERROR( cudaEventCreate( &stop ) );
    HANDLE_ERROR( cudaEventRecord( start, 0 ) );
  }
  void toc(void) {
    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime,
                                        start, stop ) );
    printf("\x1B[33mElapsed time: %8.2E ms\x1B[0m\n", elapsedTime );
    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );
  }
  void toc(const char *message) {
    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime,
                                        start, stop ) );
    printf("\x1B[33m%s: Elapsed time: %8.2E ms\x1B[0m\n", message, elapsedTime );
    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );
  }
  void toc(float *elapsedTime, const char *message) {
    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    HANDLE_ERROR( cudaEventElapsedTime( elapsedTime,
                                        start, stop ) );
    printf("\x1B[33m%s: Elapsed time: %8.2E ms\x1B[0m\n", message, *elapsedTime );
    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );
  }
  void toc(float *elapsedTime) {
    HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( stop ) );
    HANDLE_ERROR( cudaEventElapsedTime( elapsedTime,
                                        start, stop ) );
    HANDLE_ERROR( cudaEventDestroy( start ) );
    HANDLE_ERROR( cudaEventDestroy( stop ) );
  }
};
#endif
@

\subsection{Statistics definitions}
\label{sec:stat-defin}

Statistical moment are computed with the structure [[stats]]
\index{utilities!stats}
<<statistics definition>>=
struct stats {
  cublasHandle_t handle;
  cublasStatus_t status;
  void setup(void);
  void cleanup(void);
  float mean(const float *data, int n_data);
  float mean(const float *data, mask *M, int n_data);
  float var(const float *data, int n_data);
  float std(const float *data, int n_data);
  float diff_var(const float *data_1, const float *data_2, int n_data);
  float diff_std(const float *data_1, const float *data_2, int n_data);
  float var(const float *data1, mask *M, int n_data);
  float std(const float *data1, mask *M, int n_data);
  float diff_var(const float *data1, const float *data2, mask *M, int n_data);
  float diff_std(const float *data1, const float *data2, mask *M, int n_data);
};
@ 
\index{utilities!stats!setup}
<<statistics functions>>=
void stats::setup(void) {
  cublasCreate(&handle);
}
@ 
\index{utilities!stats!cleanup}
<<statistics functions>>=
void stats::cleanup(void) {
  cublasDestroy(handle);
}
<<statistics functions>>=
<<mean>> 
<<mean filtered data>>
<<variance>> 
<<differential variance>> 
<<standart deviation>>
<<differential standart deviation>>
<<variance filtered data>>
<<standart deviation filtered data>>
<<differential variance filtered data>>
<<differential standard deviation filtered data>>
@  
Computing the mean of a vector [[data]]$=x$, $\bar x = \langle x \rangle$:
\index{utilities!stats!mean}
<<mean>>=
float stats::mean(const float *data, int n_data) {
  float results, *ones;
  HANDLE_ERROR( cudaMalloc((void**)&ones, sizeof(float)*n_data ) );
  dim3 blockDim(256);
  dim3 gridDim(n_data/256+1);
  fill_ones LLL gridDim,blockDim RRR (ones, n_data);
  CUBLAS_ERROR( cublasSdot(handle, n_data, data, 1, ones, 1, &results) );
  HANDLE_ERROR( cudaFree( ones ) );
  return results/n_data;
}
@  with the kernel
<<fill ones>>=
__global__ void fill_ones(float *ones, int n_data) 
{
int i;
i = blockIdx.x * blockDim.x + threadIdx.x;
 if (i<n_data)
       ones[i] = 1.;
}
<<fill ones (char)>>=
__global__ void fill_ones_char(char *ones, int n_data) 
{
int i;
i = blockIdx.x * blockDim.x + threadIdx.x;
 if (i<n_data)
       ones[i] = 1;
}
@ %def fill_ones fill_ones_char  
Computing the mean of a vector [[data]]$=x\times f$, $\bar xf = \langle xf \rangle$ where $f$ is a filtering vector the same size than $x$ and made of 0 and 1:
\index{utilities!stats!mean}
<<mean filtered data>>=
float stats::mean(const float *data, mask *M, int n_data) {
  float results;
  CUBLAS_ERROR( cublasSdot(handle, n_data, data, 1, M->f, 1, &results) );
  return results/M->nnz;
}

@  
Computing the variance of a vector [[data]]$=x$, $\sigma^2_x=\langle x^2 \rangle - \langle x \rangle^2$:
\index{utilities!stats!var}
<<variance>>=
float stats::var(const float *data, int n_data) {
  <<variance algorithm>>
  return results;
}
@  with 
<<variance algorithm>>=
float results, mean_data, *ones;
HANDLE_ERROR( cudaMalloc((void**)&ones, sizeof(float)*n_data ) );
dim3 blockDim(256);
dim3 gridDim(n_data/256+1);
fill_ones LLL gridDim,blockDim RRR (ones, n_data);
CUBLAS_ERROR( cublasSdot(handle, n_data, data, 1, ones, 1, &mean_data) );
HANDLE_ERROR( cudaFree( ones ) );
mean_data /= n_data;
CUBLAS_ERROR( cublasSdot(handle, n_data, data, 1, data, 1, &results) );
results /= n_data;
results -= (mean_data*mean_data);
@  
Computing the standart deviation $\sigma_x$ of a vector [[data]]$=x$:
\index{utilities!stats!std}
<<standart deviation>>=
float stats::std(const float *data, int n_data) {
  <<variance algorithm>>
  return sqrtf(results);
}
@  
Computing the variance of a vector [[data]]$=x\times f$, $\sigma^2_{xf}=\langle (xf)^2 \rangle - \langle xf \rangle^2$ where $f$ is a filtering vector the same size than $x$ and made of 0 and 1:
\index{utilities!stats!var}
<<variance filtered data>>=
float stats::var(const float *data1, mask *M, int n_data) {
  float results, mean_data, *data;
  HANDLE_ERROR( cudaMalloc((void**)&data, sizeof(float)*n_data ) );
  CUBLAS_ERROR( cublasScopy(handle, n_data, data1, 1, data, 1) );
  <<variance filtered data steps>>
  return results;
}
@ 
\index{utilities!stats!diff\_var}
<<differential variance filtered data>>=
float stats::diff_var(const float *data1, const float *data2, mask *M, int n_data) {
  float results, mean_data, *data;
  HANDLE_ERROR( cudaMalloc((void**)&data, sizeof(float)*n_data ) );
  float alpha = -1;
  CUBLAS_ERROR( cublasScopy(handle, n_data, data1, 1, data, 1) );
  CUBLAS_ERROR( cublasSaxpy(handle, n_data, &alpha, data2, 1, data, 1) );
  <<variance filtered data steps>>
  return results;
}
@
\index{utilities!stats!diff\_std}
<<differential standard deviation filtered data>>=
float stats::diff_std(const float *data1, const float *data2, mask *M, int n_data) {
  float results, mean_data, *data;
  HANDLE_ERROR( cudaMalloc((void**)&data, sizeof(float)*n_data ) );
  float alpha = -1;
  CUBLAS_ERROR( cublasScopy(handle, n_data, data1, 1, data, 1) );
  CUBLAS_ERROR( cublasSaxpy(handle, n_data, &alpha, data2, 1, data, 1) );
  <<variance filtered data steps>>
  return sqrtf(results);
}
@ the variance algorithm can now be applied on the filtered data
<<variance filtered data steps>>=
CUBLAS_ERROR( cublasSdot(handle, n_data, data, 1, M->f, 1, &mean_data) );
mean_data /= M->nnz;
CUBLAS_ERROR( cublasSdot(handle, n_data, data, 1, data, 1, &results) );
results /= M->nnz;
results -= (mean_data*mean_data);
HANDLE_ERROR( cudaFree( data ) );
@  
Computing the standard deviation of a vector [[data]]$=x\times f$, $\sigma_{xf}=\left(\langle (xf)^2 \rangle - \langle xf \rangle^2\right)^{1/2}$ where $f$ is a filtering vector the same size than $x$ and made of 0 and 1:
\index{utilities!stats!std}
<<standart deviation filtered data>>=
float stats::std(const float *data1, mask *M, int n_data) {
  float results, mean_data, *data;
  HANDLE_ERROR( cudaMalloc((void**)&data, sizeof(float)*n_data ) );
  CUBLAS_ERROR( cublasScopy(handle, n_data, data1, 1, data, 1) );
  <<variance filtered data steps>>
  return sqrtf(results);
}

@  
Computing the differential variance of 2 vectors [[data1]]$=x_1$ and [[data2]]$=x_2$, $\sigma^2_{x_1-x_2}=\langle (x_1-x_2)^2 \rangle - \langle (x_1-x_2) \rangle^2$:
<<differential variance>>=
float stats::diff_var(const float *data1, const float *data2, int n_data) {
  <<differential variance algorithm>>
  return results;
}
@ with
<<differential variance algorithm>>=
float *data;
HANDLE_ERROR( cudaMalloc((void**)&data, sizeof(float)*n_data) );
float alpha = -1;
CUBLAS_ERROR( cublasScopy(handle, n_data, data1, 1, data, 1) );
CUBLAS_ERROR( cublasSaxpy(handle, n_data, &alpha, data2, 1, data, 1) );
<<variance algorithm>>
HANDLE_ERROR( cudaFree( data ) );
@  
Computing the differential standart deviation $\sigma_{x_1-x_2}$ of vector [[data1]]$=x_1$ and [[data2]]$=x_2$:
\index{utilities!stats!diff\_std}
<<differential standart deviation>>=
float stats::diff_std(const float *data1, const float *data2, int n_data) {
  <<differential variance algorithm>>
  return sqrtf(results);
}
@
The error function;
<<cublas error>>=
void CUBLAS_ERROR(cublasStatus_t status) {
  if (CUBLAS_VERBOSE) {
      printf("\nCUBLAS: ");
      if (status==CUBLAS_STATUS_SUCCESS)
	printf("the operation completed successfully\n");
      else if (status==CUBLAS_STATUS_NOT_INITIALIZED)
	printf("the library was not initialized\n");
      else if (status==CUBLAS_STATUS_ALLOC_FAILED)
	printf("the reduction buffer could not be allocated\n");
      else if (status==CUBLAS_STATUS_ARCH_MISMATCH)
	printf("the device does not support double-precision\n");
      else if (status==CUBLAS_STATUS_EXECUTION_FAILED)
	printf("the function failed to launch on the GPU\n");
    }
}

@ 
\subsection{Pupil definition}
\label{sec:pup-def}

\subsubsection{Circular pupil}
\label{sec:circular-pupil}


In the following, a circular pupil in a [[N]]$\times$[[N]] array is computed:
<<circular pupil>>=
__global__ void circular_pupil(char *pupil, int N, float scale)
{
  int i, j, k, l2;
  float h,x,y,r;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  k = i*N + j;
  if ( (i<N) && (j<N) ) {
    h = (float) (N-1)/2;
    x = i - h;
    y = j - h;
    r = x*x + y*y;
    l2 = N*N/4;
    l2 *= scale*scale;
    pupil[k] = (r>l2) ? 0 : 1;
  }
}
@
\subsubsection{Square pupil}
\label{sec:square-pupil}

A square pupil $P(i,j)$ of size [[N_S]] pixel and center of the pixels [[I_S]] and [[J_S]] inside a [[N]]$\times$[[N]] array is defined with
\begin{eqnarray}
  \label{eq:9}
  P(i,j) &=& 1,\quad \left| i-[[I_S]]\over [[N_S]] \right| \leq {1\over 2} \& \left| j-[[J_S]]\over [[N_S]] \right| \leq {1\over 2} \\
         &=& 0,\quad \text{elsewhere}
\end{eqnarray}
and computed with the kernel:
<<square pupil>>=
__global__ void square_pupil(char *pupil, int N, int N_S, int I_S, int J_S)
{
  int i, j, k, h,x,y;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  k = i*N + j;
  if ( (i<N) && (j<N) ) {
    h = N_S/2;
    x = i - I_S + h;
    y = j - J_S + h;
    pupil[k] = ( ( (x>=0) && (x<N_S) ) && ( (y>=0) && (y<N_S) ) ) ? 1 : 0;
  }
}
@ 
If the square pupil is rotated of angle $\theta$ degree, the kernel becomes
<<square pupil>>=
__global__ void rot_square_pupil(char *pupil, int N, float N_S, float I_0, float J_0, float theta, float I_S, float J_S)
{
  int i, j, k;
  float c, s, h, x, y, xr, yr;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  k = i*N + j;
  if ( (i<N) && (j<N) ) {
    sincosf(theta,&s,&c);
    h = (1.0*N_S)/N;
    x = (2.0*(i-I_0))/(N-1) - 1.0;
    y = (2.0*(j-J_0))/(N-1) - 1.0;
    xr = c*x + y*s;
    yr = -s*x + y*c;
    xr -= (2.0*I_S)/(N-1) - 1.0;
    yr -= (2.0*J_S)/(N-1) - 1.0;
    pupil[k] = ( ( (xr>=-h) && (xr<=h) ) && ( (yr>=-h) && (yr<=h) ) ) ? 1 : 0;
  }
}
@ 
\subsubsection{Valid lenslets}
\label{sec:valid-lenslets}

We compute the mask of lenslets mapping the pupil:
<<valid lenslet>>=
__global__ void valid_lenslet(char *mask, int N_pupil, float threshold)
{
  int i, j, n, iL, jL, sum, ai, bi, aj, bj, kL;
  float h,x,y,r,R;
  // pixel per lenslet
  n = N_pupil/N_SIDE_LENSLET;
  // lenslet index
  iL = blockIdx.x * blockDim.x + threadIdx.x;
  jL = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (iL<N_SIDE_LENSLET) & (jL<N_SIDE_LENSLET)) {

    sum = 0;
    ai = iL*n;
    bi = (iL+1)*n;
    aj = jL*n;
    bj = (jL+1)*n;
    h = (float) (N_pupil-1)/2;
    R = N_pupil*N_pupil/4;

    for (i=ai;i<bi;i++) {
      for (j=aj;j<bj;j++) {
	x = i - h;
	y = j - h;
	r = x*x + y*y;
    	sum += (r>R) ? 0 : 1;
      }
    }
 
    kL = iL*N_SIDE_LENSLET + jL;
    threshold *= n*n;
    mask[kL] = (sum>threshold) ? 1 : 0;

  }
}
@  The following routine computes the masks of valid lenslet [[lenslet_mask]] and actuators [[actuator_mask]] for a circular pupil corresponding to a square [[NL]]$\times$[[NL]] lenslet array with [[n]] pixel per lenslet.
\index{utilities!fried\_geometry\_setup}
<<Fried geometry masks>>=
char fried_geometry_setup(char *lenslet_mask, char *actuator_mask, 
			  int NL, int n, float threshold)
{
HANDLE_ERROR( cudaMemset(lenslet_mask,  0, sizeof(char)*NL*NL         ) );
HANDLE_ERROR( cudaMemset(actuator_mask, 0, sizeof(char)*(NL+1)*(NL+1) ) );
dim3 blockDim(16,16);
dim3 gridDim(NL/16+1,NL/16+1);
fried_geometry LLL gridDim,blockDim RRR (lenslet_mask, actuator_mask, NL, n, threshold);
return 1;
}
@
<<Fried geometry>>=
  __global__ void fried_geometry(char *lenslet_mask, char *actuator_mask, 
				 int NL, int n, float threshold)
{
  int i, j, iL, jL, sum, ai, bi, aj, bj, kL, kA, N_pupil;
  float h,x,y,r,R;
  // total pixel
  N_pupil = n*NL;
  // lenslet index
  iL = blockIdx.x * blockDim.x + threadIdx.x;
  jL = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (iL<NL) & (jL<NL)) {

    sum = 0;
    ai = iL*n;
    bi = (iL+1)*n;
    aj = jL*n;
    bj = (jL+1)*n;
    h = (float) (N_pupil-1)/2;
    R = N_pupil*N_pupil/4;

    for (i=ai;i<bi;i++) {
      for (j=aj;j<bj;j++) {
	x = i - h;
	y = j - h;
	r = x*x + y*y;
    	sum += (r>R) ? 0 : 1;
      }
    }
 
    threshold *= n*n;

    if (sum>threshold) {

      kL = iL*NL + jL;
      lenslet_mask[kL] = 1;

      kA = lenset2actuator(iL,jL,NL, 1, 1);     
      actuator_mask[kA] = 1;
      kA = lenset2actuator(iL,jL,NL,-1, 1);     
      actuator_mask[kA] = 1;
      kA = lenset2actuator(iL,jL,NL,-1,-1);     
      actuator_mask[kA] = 1;
      kA = lenset2actuator(iL,jL,NL, 1,-1);     
      actuator_mask[kA] = 1;

    }

  }
}
@ 
The following routine computes the masks of valid lenslet [[lenslet_mask]] and actuators [[actuator_mask]] for an arbitrary pupil corresponding to a square [[NL]]$\times$[[NL]] lenslet array with [[n]] pixel per lenslet.
\index{utilities!fried\_geometry\_setup\_vs\_pupil}
<<Fried geometry masks (arbitrary pupil)>>=
char fried_geometry_setup_vs_pupil(char *lenslet_mask, char *actuator_mask, 
			  int NL, int n, float threshold, char *pupil)
{
HANDLE_ERROR( cudaMemset(lenslet_mask,  0, sizeof(char)*NL*NL         ) );
HANDLE_ERROR( cudaMemset(actuator_mask, 0, sizeof(char)*(NL+1)*(NL+1) ) );
dim3 blockDim(16,16);
dim3 gridDim(NL/16+1,NL/16+1);
fried_geometry_vs_pupil LLL gridDim,blockDim RRR (lenslet_mask, actuator_mask, NL, n, threshold, pupil);
return 1;
}
@
<<Fried geometry (arbitrary pupil)>>=
__global__ void fried_geometry_vs_pupil(char *lenslet_mask, char *actuator_mask, 
				 const int NL, const int n, const float threshold, const char *pupil)
{
  int i, j, k, iL, jL, sum, ai, bi, aj, bj, kL, kA, N_pupil;
  float threshold_px;
  // total pixel
  N_pupil = n*NL;
  // lenslet index
  iL = blockIdx.x * blockDim.x + threadIdx.x;
  jL = blockIdx.y * blockDim.y + threadIdx.y;
  if ( (iL<NL) & (jL<NL)) {

    sum = 0;
    ai = iL*n;
    bi = (iL+1)*n;
    aj = jL*n;
    bj = (jL+1)*n;

    for (i=ai;i<bi;i++) {
      for (j=aj;j<bj;j++) {
        k = i*N_pupil + j;
    	sum += pupil[k];
      }
    }
 
    threshold_px = threshold*n*n;

    if (sum>threshold_px) {

      kL = iL*NL + jL;
      lenslet_mask[kL] = 1;

      kA = lenset2actuator(iL,jL,NL, 1, 1);     
      actuator_mask[kA] = 1;
      kA = lenset2actuator(iL,jL,NL,-1, 1);     
      actuator_mask[kA] = 1;
      kA = lenset2actuator(iL,jL,NL,-1,-1);     
      actuator_mask[kA] = 1;
      kA = lenset2actuator(iL,jL,NL, 1,-1);     
      actuator_mask[kA] = 1;

    }

  }
}
@
<<lenslet 2 actuator>>=
__device__ inline int lenset2actuator(int iL, int jL, int NL, int io, int jo)
{
  int iL_p, jL_p, iL_pp, jL_pp, iA, jA;

  iL_p = 2*iL + 1;
  jL_p = 2*jL + 1;
  
  iL_pp = iL_p + io;
  jL_pp = jL_p + jo;
  
  iA = iL_pp/2;
  jA = jL_pp/2;
  
  return iA*(NL+1) + jA;

}
@

\subsection{The Giant Magellan Telescope pupil}
\label{sec:gmt-pupil}

The GMT pupil is made of 7 circular segments of diameter  $D_s=8.417$m.
The clear aperture diameter is $D_{sca}=8.365$m.
There is one segment centered on the optical axis with a $D_h=1.78$m diameter hole in the middle.
The six other segments are evenly located around the center one and tilted of $\theta_s = 13.522^{\circ}$.
The distance from the optical axis to the center of an off--axis segment is $L=8.774$m.
The ASM and the FSM off--axis segments have a central hole of diameter $D_{sh}=0.559m$ and $D_{sh}=0.07$m, respectively.
<<GMT pupil>>=
__global__ void GMT_pupil(char *pupil, int *piston_mask, const int N, 
			  const float S, float _Dsh, const float c, const char v)
{
  int i, j, k, seg;
  float h,x,y,xp,yp, co_seg, so_seg, cOs, Dh, Dsca, Ds, Dsh, L;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  k = i*N + j;
  seg = blockIdx.z;
  if ( (i<N) && (j<N) ) {
    Dh   = c*2.300;
    Dsh  = _Dsh;
    Dsca = c*8.365;
    Ds   = c*8.417;
    L    = c*8.774;
    <<GMT pupil coordinates>>
    if (seg==6) {
      <<GMT pupil definition 0>>
    } else {
      <<GMT pupil definition>>
    }
  }
}
@
Lets call $x$ and $y$ the coordinates in the pupil plane.
<<GMT pupil coordinates>>=
h = S/(N-1);
x = (j - 0.5*(N-1))*h;
y = (i - 0.5*(N-1))*h;
Dsca *= Dsca;
@ 
The center segment is defined with
\begin{equation*}
   D_h^2 \le  4\left( x^2 + y^2 \right)  \le D_{sca}^2.
\end{equation*}
<<GMT pupil definition 0>>=
h = (x*x + y*y)*4.0;
Dh   *= Dh;
if ( (h>=Dh) && (h<=Dsca) ) {
  piston_mask[k] = seg+1;
  pupil[k] = v;
}
@ For the off--axis segments, the coordinates need to be projected on the tilted-planes of the segments.
The segments are tilted with respect to a line tangent to the segment at the point where the line joining the center of the segment to the optical axis crossed the segment edge.
The coordinate transformations are the following
\begin{enumerate}
\item the X--axis is aligned with the line joining the center of the segment to the optical axis:
  \begin{eqnarray*}
    x^{\prime} &=& x\cos(\theta_k) - y\sin(\theta_k) \\
    y^{\prime} &=& x\sin(\theta_k) + y\cos(\theta_k) 
  \end{eqnarray*}
where $$\theta_k=\left(2k-3\right){\pi\over 6},\quad\forall k=1,\dots,6$$ is the angle between the X--axis and the line joining the center of the $k$th segment,
<<GMT pupil definition>>=
sincospif( (2*seg-3)/6.0, &so_seg, &co_seg);
xp = x*co_seg - y*so_seg;
yp = x*so_seg + y*co_seg;
@ \item the origin of coordinates is moved to the rotation point:
  \begin{equation*}
    x^{\prime} \leftarrow x^{\prime} - \left( L - { D_s \over 2 } \right)
  \end{equation*}
<<GMT pupil definition>>=
xp -= L - 0.5*Ds;
@ \item the coordinates are projected on the tilted planes:
  \begin{eqnarray*}
    x^{\prime} &\leftarrow & {x^{\prime} \over \cos(\theta_s) } \\
    y^{\prime} &\leftarrow &  {y^{\prime} \over \cos(\theta_s) }
  \end{eqnarray*}
<<GMT pupil definition>>=
cOs = cospif( 13.522/180.0 );
xp /= cOs;
yp /= cOs;
@ \item the projected coordinates are centered on the center of the segment
  \begin{equation*}
    x^{\prime} \leftarrow x^{\prime} - \left( { D_s \over 2 } \right)
  \end{equation*}
<<GMT pupil definition>>=
xp -= Ds*0.5;
@ \end{enumerate}
and the off--axis segment is given with
\begin{equation*}
  4\left( {x^{\prime}}^2 + {y^{\prime}}^2 \right) \le D_{sca}^2.
\end{equation*}
<<GMT pupil definition>>=
h = (xp*xp + yp*yp)*4.0;
Dsh *= Dsh;
if ( (h>=Dsh) && (h<=Dsca) ) {
  piston_mask[k] = seg+1;
  pupil[k] = v;
 }
@

\subsection{Mask structure}
\label{sec:mask-structure}

\index{utilities!mask}
Lets define a mask structure with the following parameters:
\begin{itemize}
\item the mask [[m]] has a type of char and contains 0 and 1,
<<mask parameters>>=
char *m;  
@ \item the filter [[f]] has a type of float and contains 0 and 1,
<<mask parameters>>=
float *f;  
@ \item the index of the non zeros elements [[idx]],
<<mask parameters>>=
int *idx;
@ \item the mask x and y size in pixel [[size_px]],
<<mask parameters>>=
int size_px[2];
@ \item the total number of element $[[nel]]=[[size_px]][0] \times [[size_px]][1]$,
<<mask parameters>>=
int nel;  
@ \item the number of non zeros elements [[nnz]],
<<mask parameters>>=
float nnz;  
@ \item the mask x and y length in meter [[size_m]],
<<mask parameters>>=
float size_m[2];
@ \item the surface area of the mask in $m^2$ [[mask_area]],
<<mask parameters>>=
float area;  
@ \item the x and y pixel scale [[delta]],
<<mask parameters>>=
float delta[2];  
@ \item the CUBLAS handle
<<mask parameters>>=
cublasHandle_t handle;
@ \item the GMT segment pupil mask
<<mask parameters>>=
int *d__piston_mask;
@ \end{itemize}

The mask structure definition is:
<<mask structure>>=
struct mask {
  <<mask parameters>>
  void setup(int n);
  void setup(int n, float L);
  void setup(int n, float L, int i_s, int j_s, int n_out);
  void setup(float n, float L, float i_0, float j_0, 
	     float theta, float i_s, float j_s, int n_out );
  void setup_circular(int n);
  void setup_circular(int n, float D);
  void setup_circular(int n, float D, float scale);
  void setup_GMT(int n, float S);
  void set_filter(void);
  void set_filter_quiet(void);
  void set_index(void);
  void set_gmt_piston(float *phase, float *d__p);
  void alter(float *filter);
  void add(mask *other);
  void add(char *other, int other_nel);
  void add(mask *other, int offset);
  void add(char *other, int other_nel, int offset);
  void reset(void);
  void cleanup(void);
};
@ 
The mask functions are:
\begin{itemize}
\item [[setup]]
  \begin{itemize}
  \item with the total number of element in the mask [[n]]:
\index{utilities!mask!setup}
<<mask functions>>=
void mask::setup(int n) {
  nel = n;
  area = 1.0;
  <<setup common>>
  dim3 blockDim(256);
  dim3 gridDim(nel/256+1);
  fill_ones_char LLL gridDim,blockDim RRR (m,nel);
  set_filter();
}
@ with 
<<setup common>>=
HANDLE_ERROR( cudaMalloc((void**)&m, sizeof(char)*nel ) );
HANDLE_ERROR( cudaMalloc((void**)&f, sizeof(float)*nel ) );
HANDLE_ERROR( cudaMemset( m, 0, sizeof(char)*nel) );
HANDLE_ERROR( cudaMemset( f, 0, sizeof(float)*nel) );
HANDLE_ERROR( cudaMalloc((void**)&idx, sizeof(int)*nel ) );
d__piston_mask = NULL;
cublasCreate(&handle);
@ \item with a square mask of [[n]] pixel across and [[L]] meter width:
\index{utilities!mask!setup}
<<mask functions>>=
void mask::setup(int n, float L) {
  size_px[0] = size_px[1] = n;
  size_m[0]  = size_m[1]  = L;
  delta[0]   = delta[1]   = L/(n-1);
  nel = n*n;
  area = L*L;
  <<setup common>>
  dim3 blockDim(256);
  dim3 gridDim(nel/256+1);
  fill_ones_char LLL gridDim,blockDim RRR (m,nel);
  set_filter();
}
@ \item with a square mask of [[n]] pixel across and [[L]] meter width inside a [[n_out]]$\times$[[n_out]] array and centered on pixel coordinates [[i_s]] and [[j_s]]:
\index{utilities!mask!setup}
<<mask functions>>=
void mask::setup(int n, float L, int i_s, int j_s, int n_out) {
  size_px[0] = size_px[1] = n_out;
  size_m[0]  = size_m[1]  = L;
  delta[0]   = delta[1]   = L/(n-1);
  nel = n_out*n_out;
  nnz = n*n;
  area = L*L;
  <<setup common>>
  dim3 blockDim(16,16);
  dim3 gridDim(n_out/16+1,n_out/16+1);
  square_pupil LLL gridDim,blockDim RRR (m,n_out,n,i_s,j_s);
  set_filter();
}
@ \item with a square mask of [[n]] pixel across and [[L]] meter width inside a [[n_out]]$\times$[[n_out]] array and rotated around [[i_0]] and [[j_0]] of an angle $\theta$ and centered on pixel coordinates [[i_s]] and [[j_s]] :
\index{utilities!mask!setup}
<<mask functions>>=
void mask::setup(float n, float L, float i_0, float j_0, float theta, 
                 float i_s, float j_s, int n_out ) {
  size_px[0] = size_px[1] = n_out;
  size_m[0]  = size_m[1]  = L;
  delta[0]   = delta[1]   = L/(n-1);
  nel = n_out*n_out;
  nnz = n*n;
  area = L*L;
  <<setup common>>
  dim3 blockDim(16,16);
  dim3 gridDim(n_out/16+1,n_out/16+1);
  rot_square_pupil LLL gridDim,blockDim RRR (m,n_out,n,i_0,j_0,theta,i_s,j_s);
  set_filter();
}
@ \item with a circular mask inside an [[n]]$\times$[[n]] array
\index{utilities!mask!setup\_circular}
<<mask functions>>=
void mask::setup_circular(int n) {
  nel = n*n;
  area = 1.0;
  <<setup common>>
  dim3 blockDim(16,16);
  dim3 gridDim(n/16+1,n/16+1);
  circular_pupil LLL gridDim,blockDim RRR ( m , n , 1.0 );
  set_filter();
}
@ \item with a circular mask of diameter [[D]] (in meter) inside an [[n]]$\times$[[n]] array
\index{utilities!mask!setup\_circular}
<<mask functions>>=
void mask::setup_circular(int n, float D) {
  size_px[0] = size_px[1] = n;
  size_m[0]  = size_m[1]  = D;
  delta[0]   = delta[1]   = D/(n-1);
  nel = n*n;
  area = 0.25*PI*D*D;
  <<setup common>>
  dim3 blockDim(16,16);
  dim3 gridDim(n/16+1,n/16+1);
  circular_pupil LLL gridDim,blockDim RRR ( m , n , 1.0 );
  set_filter();
}
@ \item with a circular mask of diameter [[D]] (in meter) inside an [[n]]$\times$[[n]] array, [[D]] is sampled with [[n]]$\times$[[scale]] pixels
\index{utilities!mask!setup\_circular}
<<mask functions>>=
void mask::setup_circular(int n, float D, float scale) {
  size_px[0] = size_px[1] = n;
  size_m[0]  = size_m[1]  = D;
  delta[0]   = delta[1]   = D/(n-1);
  nel = n*n;
  area = 0.25*PI*D*D;
  <<setup common>>
  dim3 blockDim(16,16);
  dim3 gridDim(n/16+1,n/16+1);
  circular_pupil LLL gridDim,blockDim RRR ( m , n , scale );
  set_filter();
}
@ \item with the Giant Magellan Telescope mask inside an [[n]]$\times$[[n]] array of length [[S]] in meter
\index{utilities!mask!setup\_GMT}
<<mask functions>>=
void mask::setup_GMT(int n, float S) {
  size_px[0] = size_px[1] = n;
  size_m[0]  = size_m[1]  = S;
  delta[0]   = delta[1]   = S/(n-1);
  nel = n*n;
  area = 368.0;
  <<setup common>>
  HANDLE_ERROR( cudaMalloc((void**)&d__piston_mask, sizeof(int)*nel ) );
  HANDLE_ERROR( cudaMemset( d__piston_mask, 0, sizeof(int)*nel) );
  dim3 blockDim(16,16);
  dim3 gridDim(n/16+1,n/16+1,7);
  // PRIMARY + ASM
  GMT_pupil LLL gridDim,blockDim RRR ( m, d__piston_mask , n, S, 0.559 , 1.0        , 1);
  // ASM SHADOW
  GMT_pupil LLL gridDim,blockDim RRR ( m, d__piston_mask , n, S, 0.000 , 1.042/8.365, 0);
  /*
  // PRIMARY + FSM
  GMT_pupil LLL gridDim,blockDim RRR ( m , n, S, 0.070 , 1.0        , 1);
  // FSM SHADOW
  GMT_pupil LLL gridDim,blockDim RRR ( m , n, S, 0.000 , 1.049/8.365, 0);
  */
  set_filter();
}
@ \item
A wavefront with a different piston for each GMT segment is set with
<<mask functions>>=
void mask::set_gmt_piston(float *phase, float *d__p) {
  int n = size_px[0];
  dim3 blockDim(16,16);
  dim3 gridDim(n/16+1,n/16+1,7);
  set_gmt_piston_kernel LLL gridDim,blockDim RRR (phase, d__p, m, d__piston_mask, n);
}
@  with
<<GMT segment piston kernel>>=
__global__ void set_gmt_piston_kernel(float *phase, float *p, 
				      char *m, int *piston_mask,
				      const int N)
{
  int i, j, k, seg;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  j = blockIdx.y * blockDim.y + threadIdx.y;
  k = i*N + j;
  seg = blockIdx.z;
    if ( ( (i<N) && (j<N) ) && (m[k] && (piston_mask[k]==(seg+1)) ) ) {
    phase[k] += p[seg];
  }
}
@
\end{itemize}
\item set the filter and compute the number of non--zeros elements:
\index{utilities!mask!set\_filter}
<<mask functions>>=
void mask::set_filter(void) {
  <<set filter common>>
  INFO("@(CEO)>mask: nnz = %.0f (density: %4.2f) \n",nnz,nnz/nel);
}
@ 
\index{utilities!mask!set\_filter\_quiet}
<<mask functions>>=
void mask::set_filter_quiet(void) {
  <<set filter common>>
}
@ %def set_filter set_filter_quiet
@ with
<<set filter common>>=
dim3 blockDim(256);
dim3 gridDim(nel/256+1);
fill_ones_filtered LLL gridDim,blockDim RRR (f, m, nel);
CUBLAS_ERROR( cublasSasum(handle, nel, f, 1, &nnz) );
@
\item computes the array index of the non zeros values
\index{utilities!mask!set\_index}
<<mask functions>>=
void mask::set_index(void) {
  nnz_idx LLL 1,1 RRR (idx, m, nel); //time comsuming!!!
}
@
\item [[cleanup]]
\index{utilities!mask!cleanup}
<<mask functions>>=
void mask::cleanup(void) {
  INFO("@(CEO)>mask: freeing memory!\n");
  cublasDestroy(handle);
  HANDLE_ERROR( cudaFree( m ) );
  HANDLE_ERROR( cudaFree( f ) );
  HANDLE_ERROR( cudaFree( idx ) );
  if (d__piston_mask!=NULL)
    HANDLE_ERROR( cudaFree( d__piston_mask ) );
}
@  with the kernels
<<fill ones filtered>>=
__global__ void fill_ones_filtered(float *ones, const char* mask, int n_data) 
{
  int i;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i<n_data)
    ones[i] = (mask[i]) ? 1.0 : 0.0;
}
@ and
<<non zeros index>>=
__global__ void nnz_idx(int *idx, const char* mask, int n_data)
{
  int i, k;
  k = 0;
  for (i=0;i<n_data;i++)
    if (mask[i])
      idx[k++] = i;
}
@ %def fill_ones_filtered nnz_idx
\item the mask is wiped clean to 0 with [[reset]]
\index{utilities!reset}
<<mask functions>>=
void mask::reset(void)
{
  HANDLE_ERROR( cudaMemset( m, 0, sizeof(char)*nel) );
  HANDLE_ERROR( cudaMemset( f, 0, sizeof(float)*nel) );
}
@
\item the mask can be altered by passing an array of 1 and 0 that is multiplied by the existing mask
\index{utilities!mask!alter}
<<mask functions>>=
void mask::alter(float *filter)
{
  dim3 blockDim(256);
  dim3 gridDim(nel/256+1);
  fill_ones_altered LLL gridDim,blockDim RRR (f, m, nel, filter);
  CUBLAS_ERROR( cublasSasum(handle, nel, f, 1, &nnz) );
  INFO("@(CEO)>mask: nnz = %.0f (density: %4.2f) \n",nnz,nnz/nel);
}
@  with the kernel
<<fill ones altered>>=
__global__ void fill_ones_altered(float *ones, char* mask, int n_data, float *filter) 
{
  int i;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i<n_data) {
    ones[i] *= filter[i];
    mask[i] = (ones[i]>0) ? 1 : 0;
  }
}
@
\item masks can be added together with 
\index{utilities!mask!add}
<<mask functions>>=
void mask::add(mask *other)
{
  if (other->nel>nel)
    ERROR("Addition of masks failed because of size mismatch!");
  dim3 blockDim(256);
  dim3 gridDim(other->nel/256+1);
  add_kernel LLL gridDim,blockDim RRR (m, other->m, other->nel);  
  set_filter();
}
void mask::add(char *other, int other_nel)
{
  if (other_nel>nel)
    ERROR("Addition of masks failed because of size mismatch!");
  dim3 blockDim(256);
  dim3 gridDim(other_nel/256+1);
  add_kernel LLL gridDim,blockDim RRR (m, other, other_nel);  
  set_filter();
}
@ with the kernel
<<add mask kernel>>=
__global__ void add_kernel(char *mask, char *other_mask, int n_data)
{
  int i;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  if ( (i<n_data) && (other_mask[i]) ) {
    mask[i] = other_mask[i];
  }
}
@
\item the mask addition can be done with an offset on the mask
\index{utilities!mask!add}
<<mask functions>>=
void mask::add(mask *other, int offset)
{
  if (other->nel>nel-offset)
    ERROR("Addition of masks failed because of size mismatch!");
  dim3 blockDim(256);
  dim3 gridDim(other->nel/256+1);
  add_kernel LLL gridDim,blockDim RRR (m+offset, other->m, other->nel);  
  set_filter();
}
void mask::add(char *other, int other_nel, int offset)
{
  if (other_nel>nel-offset)
    ERROR("Addition of masks failed because of size mismatch!");
  dim3 blockDim(256);
  dim3 gridDim(other_nel/256+1);
  add_kernel LLL gridDim,blockDim RRR (m+offset, other, other_nel);  
  set_filter();
}
@
\end{itemize}

\subsection{Input/Output}
\label{sec:inputoutput}

In the following, the routines to write data on the device into a file are described.
Each file contains first an integer with the number of pages [[n_page]] and for each page the number of elements [[n_data]] and the data [[data]].
\index{utilities!dev2file}
<<device to file>>=
void dev2file(const char* filename, float* d__data, int n_data) {
  int n_page=1;
   <<device to file common start>>
  float *data;
  data = (float*)malloc(sizeof(float)*n_data);
  HANDLE_ERROR( cudaMemcpy( data, d__data,
			  sizeof(float)*n_data,
			  cudaMemcpyDeviceToHost ) );
  fwrite(&n_data,sizeof(int),1,fid);
  fwrite(data,sizeof(float),n_data,fid);
  <<device to file common end>>
}
void dev2file(const char* filename, float* d__data, int n_data, 
	      int n_data_page, int n_page) {
   <<device to file common start>>
  float *data;
  data = (float*)malloc(sizeof(float)*n_data);
  HANDLE_ERROR( cudaMemcpy( data, d__data,
			  sizeof(float)*n_data,
			  cudaMemcpyDeviceToHost ) );
  int pos=0;
  for (int k_page=0;k_page<n_page;k_page++) {
    fwrite(&n_data_page,sizeof(int),1,fid);
    fwrite(data+pos,sizeof(float),n_data_page,fid);
    pos += n_data_page;
  }
  <<device to file common end>>
}
void dev2file(const char* filename, float* d__data, int n_data, 
	      int *n_data_page, int n_page) {
   <<device to file common start>>
  float *data;
  data = (float*)malloc(sizeof(float)*n_data);
  HANDLE_ERROR( cudaMemcpy( data, d__data,
			  sizeof(float)*n_data,
			  cudaMemcpyDeviceToHost ) );
  int pos=0;
  for (int k_page=0;k_page<n_page;k_page++) {
    fwrite(n_data_page+k_page,sizeof(int),1,fid);
    fwrite(data+pos,sizeof(float),n_data_page[k_page],fid);
    pos += n_data_page[k_page];
  }
  <<device to file common end>>
}
void dev2file(const char* filename, float2* d__data, int n_data) {
  int n_page=1;
   <<device to file common start>>
  float2 *data;
  data = (float2*)malloc(sizeof(float2)*n_data);
  HANDLE_ERROR( cudaMemcpy( data, d__data,
			  sizeof(float2)*n_data,
			  cudaMemcpyDeviceToHost ) );
  fwrite(&n_data,sizeof(int),1,fid);
  fwrite(data,sizeof(float2),n_data,fid);
  <<device to file common end>>
}
void dev2file(const char* filename, int* d__data, int n_data) {
  int n_page=1;
   <<device to file common start>>
  int *data;
  data = (int*)malloc(sizeof(int)*n_data);
  HANDLE_ERROR( cudaMemcpy( data, d__data,
			  sizeof(int)*n_data,
			  cudaMemcpyDeviceToHost ) );
  fwrite(&n_data,sizeof(int),1,fid);
  fwrite(data,sizeof(int),n_data,fid);
  <<device to file common end>>
}
void dev2file(const char* filename, char* d__data, int n_data) {
  int n_page=1;
   <<device to file common start>>
  char *data;
  data = (char*)malloc(sizeof(char)*n_data);
  HANDLE_ERROR( cudaMemcpy( data, d__data,
			  sizeof(char)*n_data,
			  cudaMemcpyDeviceToHost ) );
  fwrite(&n_data,sizeof(int),1,fid);
  fwrite(data,sizeof(char),n_data,fid);
  <<device to file common end>>
}

<<device to file common start>>=
struct stat sb;
char path[50];

if (stat("bins",&sb)!=0)
  if (mkdir("bins", S_IRWXU)==0)
    fprintf(stdout,"@(CEO)>utilities: bins directory created!\n");
  else {
    fprintf(stdout,"\n\x1B[31m@(CEO)>utilities: bins mkdir failed!\x1B[0m\n");
    return;
  }
strcpy(path,"bins/");
strcat(path,filename);
FILE *fid;
fid = fopen(path,"wb");
if (fid==NULL) {
  fprintf(stdout,"\n\x1B[31m@(CEO)>utilities: Unable to open the file: %s!\x1B[0m\n",path);
  return;
 }
fwrite(&n_page,sizeof(int),1,fid);
<<device to file common end>>=
fclose(fid);
free(data);
fprintf(stdout,"@(CEO)>utilities: Data saved to %s!\n",filename);
@

\subsection{Device to host}
\label{sec:dev-to-host}

\index{utilities!dev2host}
<<data to host from device>>=
void dev2host( float *host_data, float *dev_data, int N) {
     HANDLE_ERROR( cudaMemcpy( host_data, dev_data, 
			       N*sizeof(float), cudaMemcpyDeviceToHost ) ) ;
}
void dev2host_int( int *host_data, int *dev_data, int N) {
  HANDLE_ERROR( cudaMemcpy( host_data, dev_data, 
                            N*sizeof(int), cudaMemcpyDeviceToHost ) ) ;
}
@
\index{utilities!host2dev}
<<data to device from host>>=
void host2dev( float *dev_data, float *host_data, int N) {
    HANDLE_ERROR( cudaMemcpy( dev_data, host_data, 
			       N*sizeof(float), cudaMemcpyHostToDevice ) ) ;
}
void host2dev_char( char *dev_data, char *host_data, int N) {
  HANDLE_ERROR( cudaMemcpy( dev_data, host_data, 
                            N*sizeof(char), cudaMemcpyHostToDevice ) ) ;
}
@ 
\index{utilities!freedev}
<<free device>>=
void freedev( float **dev_data ) {
    fprintf(stdout,"__ FREEDEV __\n");
    HANDLE_ERROR( cudaFree( *dev_data) );
}
@ 
\index{utilities!gpu}
<<template gpu data structure>>=
  //template<typename T> 
struct gpu_float {
  float *dev_data, *host_data, *d_tau;
  <<gpu template common>>
  void dev2dev(gpu_float *other);
  void axpy(const gpu_float *x, const float alpha);
  void scale(const float alpha);
  void mv(gpu_float *y, const gpu_float *x);
  void qr_solve(gpu_float *x, gpu_float *b);
};
struct gpu_double {
  double *dev_data, *host_data;
  <<gpu template common>>
  void mv(gpu_float *y, const gpu_float *x);
};
@
with
<<gpu template common>>=
int N, nb;
stats S;
cublasStatus_t stat;
cublasHandle_t handle;
cusolverDnHandle_t cusolverH;
void setup(void);
void setup(int N_T);
void dev_malloc(void);
void free_dev(void);
void dev2host(void);
void host2dev(void);
void reset(void);
void qr(int m);
@
<<GPU float methods>>=
void gpu_float::setup(void) {
  dev_data = NULL;
  host_data = NULL;
  d_tau = NULL;
  nb = sizeof(float);
}
void gpu_float::setup(int N_T) {
  dev_data = NULL;
  host_data = NULL;
  d_tau = NULL;
  nb = sizeof(float);
  N = N_T;
}
void gpu_float::dev_malloc(void) {
  HANDLE_ERROR( cudaMalloc( (void**)&dev_data, N*nb) );
}
void gpu_float::free_dev(void) {
  if (dev_data) HANDLE_ERROR( cudaFree( dev_data ) );
  if (d_tau) HANDLE_ERROR( cudaFree( d_tau ) );
}
void gpu_float::dev2host(void) {
  HANDLE_ERROR( cudaMemcpy( host_data, dev_data, 
                            N*nb, cudaMemcpyDeviceToHost ) ) ;
}
void gpu_float::dev2dev(gpu_float *other) {
  HANDLE_ERROR( cudaMemcpy( other->dev_data, dev_data, 
                            N*nb, cudaMemcpyDeviceToDevice ) ) ;
}
void gpu_float::host2dev(void) {
  HANDLE_ERROR( cudaMemcpy( dev_data, host_data, 
                            N*nb, cudaMemcpyHostToDevice ) ) ;
}
void gpu_float::reset(void) {
  HANDLE_ERROR( cudaMemset( dev_data, 0, N*nb) );
}
void gpu_float::mv(gpu_float *y, const gpu_float *x) {
  cublasCreate(&handle);
  float alpha=1, beta=0;
  CUBLAS_ERROR(cublasSgemv(handle,CUBLAS_OP_N,
                           y->N,x->N,
                           &alpha,
                           dev_data,y->N,
                           x->dev_data,1,
                           &beta,
                           y->dev_data,1));
  cublasDestroy(handle);
}
void gpu_float::axpy(const gpu_float *x, const float alpha) {
  cublasCreate(&handle);
  cublasSaxpy(handle, N,
              &alpha,
              x->dev_data, 1,
              dev_data, 1);
  cublasDestroy(handle);
}
void gpu_float::scale(const float alpha) {
  cublasCreate(&handle);
  cublasSscal(handle, N,
              &alpha,
              dev_data, 1);
  cublasDestroy(handle);
}
void gpu_float::qr(int m) {

  int *devInfo = NULL;
  float *d_work = NULL;
  int  lwork_geqrf = 0;
  const int lda = m;

  cusolverDnCreate(&cusolverH);
  HANDLE_ERROR( cudaMalloc ((void**)&d_tau, sizeof(float) * m) );
  HANDLE_ERROR( cudaMalloc ((void**)&devInfo, sizeof(int)) );

  HANDLE_ERROR_CUSOLVER( cusolverDnSgeqrf_bufferSize(
                                                cusolverH,
                                                m,
                                                N/m,
                                                dev_data,
                                                lda,
                                                &lwork_geqrf) );

  HANDLE_ERROR( cudaMalloc((void**)&d_work, sizeof(double)*lwork_geqrf) );

  /* compute QR factorization */
  HANDLE_ERROR_CUSOLVER( cusolverDnSgeqrf(
                                     cusolverH, 
                                     m, 
                                     N/m, 
                                     dev_data, 
                                     lda, 
                                     d_tau, 
                                     d_work, 
                                     lwork_geqrf, 
                                     devInfo) );
  HANDLE_ERROR( cudaDeviceSynchronize());

  int info_gpu = 0;
  HANDLE_ERROR( cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost));

  INFO("after geqrf: info_gpu = %d\n", info_gpu);

  if (d_work) HANDLE_ERROR( cudaFree( d_work ) );
  if (devInfo) HANDLE_ERROR( cudaFree( devInfo ) );
  cusolverDnDestroy(cusolverH);

}
void gpu_float::qr_solve(gpu_float *x, gpu_float *b) {
  int *devInfo = NULL;
  float *d_work = NULL;
  int  lwork_ormqr = 0;
  int info_gpu = 0;
  const int m = b->N;
  const int lda = m;
  const int ldb = m;
  const int nrhs = 1;

  cusolverDnCreate(&cusolverH);
  cublasCreate(&handle);

  HANDLE_ERROR( cudaMalloc ((void**)&devInfo, sizeof(int)) );

  HANDLE_ERROR_CUSOLVER( cusolverDnSormqr_bufferSize(
      cusolverH,
      CUBLAS_SIDE_LEFT,
      CUBLAS_OP_T,
      m,
      nrhs,
      N/m,
      dev_data,
      lda,
      d_tau,
      b->dev_data,
      ldb,
      &lwork_ormqr));

  HANDLE_ERROR( cudaMalloc((void**)&d_work, sizeof(float)*lwork_ormqr) );

  /* compute Q^T*B */
  HANDLE_ERROR_CUSOLVER(cusolverDnSormqr(
                                    cusolverH,
                                    CUBLAS_SIDE_LEFT,
                                    CUBLAS_OP_T,
                                    m,
                                    nrhs,
                                    N/m,
                                    dev_data,
                                    lda,
                                    d_tau,
                                    b->dev_data,
                                    ldb,
                                    d_work,
                                    lwork_ormqr,
                                    devInfo));
  HANDLE_ERROR( cudaDeviceSynchronize());

  /* check if QR is good or not */
  HANDLE_ERROR( cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost));

   INFO("after ormqr: info_gpu = %d\n", info_gpu);

  /* compute x = R \ Q^T*B */
  float one = 1;
  //assert(cudaSuccess == cudaStat1);
  ( cublasStrsm(
                              handle,
                              CUBLAS_SIDE_LEFT,
                              CUBLAS_FILL_MODE_UPPER,
                              CUBLAS_OP_N, 
                              CUBLAS_DIAG_NON_UNIT,
                              N/m,
                              nrhs,
                              &one,
                              dev_data,
                              lda,
                              b->dev_data,
                              N/m));
  HANDLE_ERROR( cudaDeviceSynchronize());

  HANDLE_ERROR(cudaMemcpy(x->dev_data, b->dev_data, sizeof(float)*N/m, cudaMemcpyDeviceToDevice));


  if (d_work)  HANDLE_ERROR( cudaFree( d_work ) );
  if (devInfo) HANDLE_ERROR( cudaFree( devInfo ) );
  cusolverDnDestroy(cusolverH);
  cublasDestroy(handle);
}
@
<<GPU double methods>>=
void gpu_double::setup(void) {
  nb = sizeof(double);
}
void gpu_double::setup(int N_T) {
  nb = sizeof(double);
  N = N_T;
}
void gpu_double::dev_malloc(void) {
  HANDLE_ERROR( cudaMalloc( (void**)&dev_data, N*nb) );
}
void gpu_double::free_dev(void) {
  HANDLE_ERROR( cudaFree( dev_data ) );
}
void gpu_double::dev2host(void) {
  HANDLE_ERROR( cudaMemcpy( host_data, dev_data, 
                            N*nb, cudaMemcpyDeviceToHost ) ) ;
}
void gpu_double::host2dev(void) {
  HANDLE_ERROR( cudaMemcpy( dev_data, host_data, 
                            N*nb, cudaMemcpyHostToDevice ) ) ;
}
void gpu_double::reset(void) {
  HANDLE_ERROR( cudaMemset( dev_data, 0, N*nb) );
}
@
<<GPU methods>>=
  /*
  void double2float(gpu_float *cuFloatArray ) {
    dim3 blockDim(256);
    dim3 gridDim(N/256+1);
    double2float_kern LLL gridDim,blockDim RRR (cuFloatArray->dev_data, dev_data, N);
    }*/

@ with
<<double2float kernel>>=
__global__ void double2float_kern(float *d__single_data, double *d__double_data, const int N)
{
  int i;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i<N)
    d__single_data[i] =  __double2float_rn ( d__double_data[i] ) ;
}
@

\subsubsection{QR factorization}
\label{sec:qr_factor}

<<QR factorization>>=
  void geqrf(float *d_tau, float *a, int m, int n) {

  int *devInfo = NULL;
  float *d_work = NULL;
  int  lwork_geqrf = 0;
  const int lda = m;
  cusolverDnHandle_t cusolverH;
  cusolverDnCreate(&cusolverH);
  //HANDLE_ERROR( cudaMalloc ((void**)&d_tau, sizeof(float) * m) );
  HANDLE_ERROR( cudaMalloc ((void**)&devInfo, sizeof(int)) );

  HANDLE_ERROR_CUSOLVER( cusolverDnSgeqrf_bufferSize(
                                                     cusolverH,
                                                     m,
                                                     n,
                                                     a,
                                                     lda,
                                                     &lwork_geqrf));

  HANDLE_ERROR( cudaMalloc((void**)&d_work, sizeof(double)*lwork_geqrf) );

  /* compute QR factorization */
  HANDLE_ERROR_CUSOLVER( cusolverDnSgeqrf(
                                          cusolverH, 
                                          m, 
                                          n, 
                                          a, 
                                          lda, 
                                          d_tau, 
                                          d_work, 
                                          lwork_geqrf, 
                                          devInfo));
  HANDLE_ERROR( cudaDeviceSynchronize());

  int info_gpu = 0;
  HANDLE_ERROR( cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost));

  INFO("after geqrf: info_gpu = %d\n", info_gpu);
  //assert(0 == info_gpu);

  if (d_work) HANDLE_ERROR( cudaFree( d_work ) );
  if (devInfo) HANDLE_ERROR( cudaFree( devInfo ) );
  cusolverDnDestroy(cusolverH);
}
<<compute Q^T*B>>=
  void ormqr(float *b, int m, float *q, float *tau, int n)
{
  int *devInfo = NULL;
  float *d_work = NULL;
  int  lwork_ormqr = 0;
  int info_gpu = 0;
  const int lda = m;
  const int ldb = m;
  const int nrhs = 1;
  cusolverDnHandle_t cusolverH;
  cusolverDnCreate(&cusolverH);

  HANDLE_ERROR( cudaMalloc ((void**)&devInfo, sizeof(int)) );

  HANDLE_ERROR_CUSOLVER( cusolverDnSormqr_bufferSize(
                                                     cusolverH,
                                                     CUBLAS_SIDE_LEFT,
                                                     CUBLAS_OP_T,
                                                     m,
                                                     1,
                                                     n,
                                                     q,
                                                     lda,
                                                     tau,
                                                     b,
                                                     m,
                                                     &lwork_ormqr));


  //lwork = (lwork_geqrf > lwork_ormqr)? lwork_geqrf : lwork_ormqr;
  HANDLE_ERROR( cudaMalloc((void**)&d_work, sizeof(float)*lwork_ormqr) );

  /* compute Q^T*B */
  HANDLE_ERROR_CUSOLVER(cusolverDnSormqr(
                                         cusolverH,
                                         CUBLAS_SIDE_LEFT,
                                         CUBLAS_OP_T,
                                         m,
                                         1,
                                         n,
                                         q,
                                         lda,
                                         tau,
                                         b,
                                         m,
                                         d_work,
                                         lwork_ormqr,
                                         devInfo));
  HANDLE_ERROR( cudaDeviceSynchronize());

  /* check if QR is good or not */
  HANDLE_ERROR(cudaMemcpy(&info_gpu, devInfo, sizeof(int), cudaMemcpyDeviceToHost));

   INFO("after ormqr: info_gpu = %d\n", info_gpu);

   if (d_work) HANDLE_ERROR( cudaFree( d_work ) );
   if (devInfo) HANDLE_ERROR( cudaFree( devInfo ) );
   cusolverDnDestroy(cusolverH);
}
@

\subsection{Atmospheric dispersion}
\label{sec:atmosph-disp}

The optical refractive index of air $n-1$ is given as a function of the wavelength [$\mu$m], the altitude [m], the temperature [$^{\circ}$C] and the humidity [\%].
It is defined with the model of Owens\cite{Owens67} (\url{http://www.astrosurf.com/prostjp/Dispersion_en.html}) given by
\begin{eqnarray*}
  \label{eq:5}
  \lefteqn{ \left(n - 1\right) \times 10^8 = } \\
  && \left[ 2371.34 + {683939.7 \over 130 - \sigma^2} + {4547.3 \over 38.9 - \sigma^2} \right]D_s + \\
  && \left[ 6487.31 + 58.05\sigma^2 - 0.71150\sigma^4 + 0.08851\sigma^6 \right]D_w,
\end{eqnarray*}
in which $\sigma=1/\lambda$, $\lambda$ is the wavelength in micron.
\index{utilities!atmosphere\_refractive\_index}
<<atmosphere refractive index>>=
double atmosphere_refractive_index(float wavelength, float altitude, 
				   float temperature, float humidity)
  {
    double sigma2, n_minus_one, Ds, Dw, T, Ps, Pw, H, P0;

    sigma2 = 1.0/wavelength/wavelength;
    T = temperature + 273;
    H = humidity/100.0;
    P0 = 1013.25;

    <<dry air partial pressure>>
    <<density of dry air>>
    <<water vapor partial pressure>>
    <<density of water vapor>>

    n_minus_one = ( 2371.34 + 683939.7 / ( 130.0 - sigma2) + 
		    4547.3 / ( 38.9 - sigma2 ) ) * Ds +
    ( 6487.31 + sigma2 * ( 58.05 - sigma2 * ( 0.71150 - 0.08851*sigma2 ) ) ) * Dw;

    return n_minus_one*1e-8;
  }					
@ 
$D_s$ is the density of dry air:
\begin{equation*}
  \label{eq:6}
  D_s = { P_s \over T} \left[ 1 + P_s\left( 57.90 \times 10^{-8} - { 9.3250 \times 10^{-4} \over T } + {0.25844 \over T^2} \right) \right],
\end{equation*}
<<density of dry air>>=
Ds = (Ps/T)*( 1 + Ps*( 57.90e-8 - ( 9.3250e-4  - 0.25844/T )/T ) );
@ 
$D_w$ is the density of water vapor:
\begin{eqnarray*}
  \lefteqn{ D_w =  {P_w\over T} \Bigg\{ 1 + P_w \left[ 1 + \left( 3.7\times 10^{-4} \right)P_w \right] \times }\\
  && \left[ -2.37321 \times 10^{-3} + {2.23366 \over T } - {710.792 \over T^2 } + {7.75141\times 10^4 \over T^3} \right] \Bigg\}
\end{eqnarray*}
<<density of water vapor>>=
Dw = (Pw/T)*( 1 + Pw*( 1 + 3.7e-4*Pw ) )*
    ( -2.37321e-3 + ( 2.23366  - ( 710.792 - 7.75141e4/T )/T )/T );
@ 
where $P_s$ is the dry air partial pressure
\begin{equation*}
  \label{eq:7}
  P_s = P_0\left( 1 -  (2.25577 \times 10^{-5} ) h \right)^{5.25588},
\end{equation*}
with $P_0=1013.25$mb and $h$ the altitude in meter
<<dry air partial pressure>>=
Ps = P0*pow( 1 - 2.25577e-5*altitude , 5.25588 );
@ 
and $P_w$ is the water vapor partial pressure
\begin{eqnarray*}
  P_w &=& P_{ws} { H \over 100 } \\
  P_{ws} &=& T^{-8.2}\exp\left( 77.3450 + 0.0057 T - {7235 \over T} \right)
\end{eqnarray*}
with $H$ the relative humidity in percent and $T$ is the temperature in degree Kelvin.
<<water vapor partial pressure>>=
Pw = H*exp( 77.3450 + 0.0057*T - 7235.0/T )/pow(T,8.2);
@ 
The refraction angle is given by
\begin{equation*}
  \varepsilon(\lambda) = \left[ n(\lambda)-1 \right] \tan(z)
\end{equation*}
where $z$ is the zenith angle.
The atmospheric dispersion $\Xi$ corresponding to a spectral filter of width $\Delta\lambda$ with central wavelength $\lambda_c$ is written
\begin{equation*}
  \Xi = \left[ n\left( \lambda_c - {\Delta\lambda\over2} \right) - n\left( \lambda_c + {\Delta\lambda\over2} \right) \right] \tan(z)
\end{equation*}
\index{utilities!atmosphere\_dispersion}
<<atmospheric dispersion>>=
double atmospheric_dispersion(float wavelength, float delta_wavelength,
                                   float zenith, float altitude, 
				   float temperature, float humidity)
{
double n1, n2, xi;
float lambda1, lambda2;
lambda1 = wavelength - 0.5*delta_wavelength;
lambda2 = wavelength + 0.5*delta_wavelength;
n1 =  atmosphere_refractive_index(lambda1, altitude, 
				   temperature, humidity);
n2 =  atmosphere_refractive_index(lambda2, altitude, 
				   temperature, humidity);
xi = (n1-n2)*tan(zenith);
return xi;
}
@


\subsection{Wavefront differentiation}
\label{sec:wavefr-diff}

The following routine is computing the finite difference of a given wavefront on a [[N_L]]$\times$[[N_L]] lenslet array of pitch [[d]].
The wavefront is sample with [[n+1]]$\times$[[n+1]] points per lenslet ([[n]]$\geq 2$) meaning that the full wavefront is of size $[[n]][[N_L]]+1$.

The x and y slopes corresponding to the average gradient of the wavefront $\varphi$ on the lenslet $(i_L,j_L)$ with $0\leq i_l,j_l <$[[N_L]] are given by:
\begin{eqnarray}
  \label{eq:8}
  s_x(i_L,j_L) &=& {1\over (n+1)d} \sum_{k=0}^{n} \varphi\left[ i_Ln+k,j_Ln \right] -  \varphi\left[ i_Ln+k,(j_L+1)n\right] \\
  s_x(i_L,j_L) &=& {1\over (n+1)d} \sum_{k=0}^{n} \varphi\left[ i_Ln,j_Ln+k \right] -  \varphi\left[ (i_L+1)n,j_Ln+k\right] 
\end{eqnarray}

\index{utilities!wavefront\_finite\_difference}
<<wavefront differentiation>>=
void wavefront_finite_difference(float *sx, float *sy, int NL, 
                                 float *phi, int n, float d, int N_GS)
  {
    dim3 blockDim(16,16);
    dim3 gridDim(NL/16+1,NL/16+1,N_GS);
    wavefront_finite_difference_kernel LLL gridDim,blockDim RRR (sx, sy, NL, phi, n, d);
  }
@ %def wavefront_finite_difference
@ with the kernel
<<wavefront differentiation kernel>>=
__global__ void wavefront_finite_difference_kernel(float *sx, float *sy, int NL, 
                                                   float *phi, int n, float d)
  {
    int iL, jL, kL, k, u, v, w, w0, NP, iSource;
    iL = blockIdx.x * blockDim.x + threadIdx.x;
    jL = blockIdx.y * blockDim.y + threadIdx.y;
    iSource = blockIdx.z;
    if ( (iL<NL) && (jL<NL) ) {
      <<wavefront differentiation core>>
    }
  }
@  and
<<wavefront differentiation core>>=
 kL = iL*NL+jL;
 NP = n*NL + 1;
 kL += 2*iSource*NL*NL;
 u = iL*n;
 v = jL*n;
 w0 = iSource*NP*NP;
 sx[kL] = sy[kL] = 0.0;
 for (k=0;k<=n;k++) {
   // x slope
   w = w0 + u*NP + v + k;
   sx[kL] += phi[w];
   w += n*NP;
   sx[kL] -= phi[w];
   // y slope
   w = w0 + (u+k)*NP + v;
   sy[kL] += phi[w];
   w += n;
   sy[kL] -= phi[w];
 }
 sx[kL] /= (n+1)*d;
 sy[kL] /= (n+1)*d;
@ 
The wavefront averaged-gradient can also be derived with a sparse gradient matrix operator,
<<wavefront differentiation sparse matrix kernel>>=
__global__ void wavefront_finite_difference_sparse_matrix_kernel(
                                                   float *csrValH, 
                                                   int *csrColIndH, 
                                                   int *csrRowPtrH, 
                                                   const int NL, 
                                                   const int n, 
                                                   const float d)
  {
    int iL, jL, kL, k, u, v, w, w0, NP, 
        NL2, iSource, idx, nnz_per_axis;
    float norm;
    iL = blockIdx.x * blockDim.x + threadIdx.x;
    jL = blockIdx.y * blockDim.y + threadIdx.y;
    iSource = blockIdx.z;
    if ( (iL<NL) && (jL<NL) ) {
      <<wavefront differentiation sparse matrix core>>
    }
  }
@  where
<<wavefront differentiation sparse matrix core>>=
 NL2 = NL*NL;
 norm = 1.0/((n+1)*d);
 kL = iL*NL+jL;
 NP = n*NL + 1;
 kL += iSource*NL2;
 u = iL*n;
 v = jL*n;
 w0 = iSource*NP*NP;
 nnz_per_axis = NL2*(n+1)*2;
 idx = kL*2*(n+1) + nnz_per_axis*2*iSource;
   // x slope
 for (k=0;k<=n;k++) {
   w = w0 + u*NP + v + k;
   csrValH[idx]      = norm;
   csrColIndH[idx] = w;
   w += n*NP;
   csrValH[n+1+idx]      = -norm;
   csrColIndH[n+1+idx++] = w;
 }
csrRowPtrH[kL+1] = 2*(n+1)*(kL+1);
 idx = (NL2+kL)*2*(n+1) + nnz_per_axis*2*iSource;
  // y slope
 for (k=0;k<=n;k++) {
   w = w0 + (u+k)*NP + v;
   csrValH[idx]      = norm;
   csrColIndH[idx] = w;
   w += n;
   csrValH[n+1+idx]      = -norm;
   csrColIndH[n+1+idx++] = w;
 }
csrRowPtrH[kL+1+NL2] = 2*(n+1)*(kL+1+NL2);
   
@
\subsubsection{Masked lenslet array}
\label{sec:masked-lenslet-array}

When the lenslet array is partially illuminated, the valid lenslet mask [[valid_lenslet]] can be passed to the routine.
The slopes outside the mask are set to 0.
\index{utilities!wavefront\_finite\_difference}
<<wavefront differentiation (with mask)>>=
void wavefront_finite_difference(float *sx, float *sy, int NL, 
                                 float *phi, int n, float d,
                                 mask *valid_lenslet, int N_GS)
  {
    int sxy_nbyte = sizeof(float)*NL*NL*N_GS;
    HANDLE_ERROR( cudaMemset(sx,  0, sxy_nbyte ) );
    HANDLE_ERROR( cudaMemset(sy,  0, sxy_nbyte ) );
    dim3 blockDim(16,16);
    dim3 gridDim(NL/16+1,NL/16+1,N_GS);
    wavefront_finite_difference_kernel_masked LLL gridDim,blockDim RRR 
                                (sx, sy, NL, phi, n, d, valid_lenslet->m);
  }
@ with the kernel
<<wavefront differentiation kernel (with mask)>>=
__global__ void wavefront_finite_difference_kernel_masked(
                                 float *sx, float *sy, int NL, 
                                 float *phi, int n, float d,
                                 char *M)
  {
    int iL, jL, kL, k, u, v, w0, w, NP, iSource;
    iL = blockIdx.x * blockDim.x + threadIdx.x;
    jL = blockIdx.y * blockDim.y + threadIdx.y;
    iSource = blockIdx.z;
    kL = iL*NL+jL;
    kL += iSource*NL*NL;
    if ( ( (iL<NL) && (jL<NL) ) && (M[kL]>0) ) {
      <<wavefront differentiation core>>
    }
  }
@ 
\subsection{Sparse matrix}
\label{sec:sparse-matrix}

\index{utilities!sparseMatrix}
The sparse matrix structure holds the data to build a matrix in the compressed row storage.
<<sparse matrix>>=
struct sparseMatrix {
<<sparse matrix parameters>>
void setup(const int _n_row_, const int _n_col_, const int _nnz_);
void gradient( int NL,  int N_LENSLET_PX,  float d);
void interpolation(int NI, int NP);
void interpolation(int NI, int NP, mask *pupil,float i0, float j0);
void cleanup(void);
void MVM(float *y, float *x);
void add(sparseMatrix *C, sparseMatrix *B);
};
@ 
The sparse matrix parameters are
\begin{itemize}
\item the number of row and columns of the sparse matrix
<<sparse matrix parameters>>=
int n_row, n_col;
@
\item the number of non--zeros values
<<sparse matrix parameters>>=
int nnz;
@
\item the non zeros values
<<sparse matrix parameters>>=
float *csrValH;
@
\item the column index of the non--zeros values in each row 
<<sparse matrix parameters>>=
int *csrColIndH;
@
\item the index where a row start in [[csrColIndH]]
<<sparse matrix parameters>>=
int  *csrRowPtrH;
@
\item variables for linear algebra calculation
<<sparse matrix parameters>>=
float alpha, beta;
cudaError_t cudaStat;
cusparseStatus_t status;
cusparseHandle_t handle;
cusparseMatDescr_t descr;
@
\end{itemize}
@ 
The sparse matrix methods are
\begin{itemize}
\item setup
\index{utilities!sparseMatrix!setup}
<<sparse matrix methods>>=
void sparseMatrix::setup(int _n_row_, int _n_col_, int _nnz_) 
{
n_row = _n_row_;
n_col = _n_col_;
nnz   = _nnz_;
<<sparse matrix allocation>>
}
@ 
with
<<sparse matrix allocation>>=
alpha = 1.0;
beta = 0.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)*(n_row+1) ) );
HANDLE_ERROR( cudaMemset(csrRowPtrH, 0, sizeof(int)*(n_row+1) ) );
@
\item sparse gradient matrix operator to compute the average slopes of the wavefront on an array of $[[NL]]\times[[NL]]$ lenslets with [[N_LENSET_PX]] pixel per lenslet and a lenslet pitch [[d]]:,
\index{utilities!sparseMatrix!gradient}
<<sparse matrix methods>>=
void sparseMatrix::gradient( int NL,  int N_LENSLET_PX,  float d) 
{
  n_row = 2*NL*NL;
  n_col = NL*N_LENSLET_PX + 1;
  nnz   = 2*(N_LENSLET_PX+1)*n_row;
  <<sparse matrix allocation>>
  dim3 blockDim(16,16);
  dim3 gridDim(NL/16+1,NL/16+1);
  wavefront_finite_difference_sparse_matrix_kernel LLL gridDim,blockDim RRR
  (csrValH, csrColIndH, csrRowPtrH, NL, N_LENSLET_PX, d);
}
@
\item bi--linear 2D interpolation
<<sparse matrix methods>>=
void sparseMatrix::interpolation(int NI, int NP)
{
n_row = NI;
n_col = NP;
nnz   = 4*NI*NI;
<<sparse matrix allocation>>
<<sparse building>>
}
void sparseMatrix::interpolation(int NI, int NP, mask *pupil,
                                 float i0, float j0)
{
n_row = NI*NI;
n_col = NP*NP;
nnz = 4*pupil->nnz;
<<sparse matrix allocation>>
<<sparse building with pupil mask>>
}
@ creation of the sparse matrix:
<<sparse building>>=
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearInterpSparseOperator2 LLL gridDim,blockDim RRR 
  (csrValH, csrColIndH, csrRowPtrH, NP, NI);
@ Creation of the sparse matrix with the pupil mask:
<<sparse building with pupil mask>>=
dim3 blockDim(16,16);
dim3 gridDim(NI/16+1,NI/16+1);
bilinearInterpSparseOperatorMask2 LLL gridDim,blockDim RRR 
  (csrValH, csrColIndH, csrRowPtrH, NP, NI, 
  pupil->m, pupil->nnz, i0, j0);
@ The sparse matrix is stored in the compressed sparse row format and it is computed with the kernel
<<bilinear interpolation sparse matrix kernel>>=
__global__ void bilinearInterpSparseOperator2(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;
    <<bilinear interpolation modified>>

    csrRowPtrH[k+1] = 4*(k+1);
  }
}
@ The sparse matrix is modified to take into account the pupil mask.
Where the pupil mask is zero, the corresponding row in the sparse matrix must be filled with zeros.
<<bilinear interpolation sparse matrix kernel with pupil mask>>=
__global__ void bilinearInterpSparseOperatorMask2(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;
       <<nnz position>>
       idx =  4*(pos-1) - 1;
       offset = 0;

         delta_p = delta;

         <<bilinear interpolation modified>>
         offset += NP*NP;

       csrRowPtrH[k] = 4*(pos-1);
       <<csrRowPtrH for empty rows>>
       if ( k==(NI*NI-1) )
          csrRowPtrH[k+1] = 4*nnz;
    } else { 
         if ( k==(NI*NI-1) ) {
           ++k;
           csrRowPtrH[k] = 4*nnz;
           <<csrRowPtrH for empty rows>>
        } 
     }
   }
 }
@ 

\paragraph{Sparse bi--linear interpolation}
\label{sec:sparse-bi-linear-1}


A wavefront $\hat\varphi$ is sampled on a $N_p\times N_p$ square grid with sampling step $\delta_p$.
The wavefront $\varphi_{bli}$, sampled on a square grid $N\times N$ with sampling step $\delta$ is derived from a bi--linear interpolation of $\hat\varphi$.
The coordinates of $\hat\varphi$ samples are given by
\begin{equation}
  \label{eq:9}
  \left(
    \begin{array}{c}
      x_{i_p} \\
      y_{j_p}
    \end{array}
  \right) = \delta_p \left[ \left(
    \begin{array}{c}
      i_p \\
      j_p
    \end{array}
    \right) + {1-N_p \over 2} \right], \forall i,j \in [0,\dots,N_p-1]
\end{equation}
and the samples of $\varphi_{bli}$ are located at
\begin{equation}
  \label{eq:10}
  \left(
    \begin{array}{c}
      x_i \\
      y_j
    \end{array}
  \right) = \delta \left[ \left(
    \begin{array}{c}
      i \\
      j
    \end{array}
    \right) + {1-N \over 2} \right], \forall i,j \in [0,\dots,N-1]  
\end{equation}

@ The interpolation starts with normalizing the coordinates $(x_i,y_j)$ with respect to the coordinates $(x_{i_p},y_{j_p})$
<<bilinear interpolation modified>>=
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);
@  and taking the floor integer part to locate the bottom--left coordinate of the pixels surrounding [[xi]] and [[yi]]
<<bilinear interpolation modified>>=
fs  = floorf(s);
ft  = floorf(t);
ndx = __float2int_rd( ft + fs*NP );
@ Next, one checks if the edges of the array have been reached
<<bilinear interpolation modified>>=
if (s==(NP-1)) { s += 1 - fs; ndx -= NP; } else { s -= fs; }    
if (t==(NP-1)) { t += 1 - ft; ndx -= 1; }   else { t -= ft; }
@ and finally the interpolation is computed as
<<bilinear interpolation modified>>=	 
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;
@  For a given location in the pupil, the position of the non--zero values in the non-zero value vector is given but the sum of the mask up to that point:
<<nnz position>>=
pos = 0;
for (l=0;l<=k;l++)
    pos += (mask[l]) ? 1 : 0;
@ This position is saved in [[csrRowPtrH]].
The consecutive rows above the current row corresponding to zero values in the pupil mask must have their [[csrRowPtrH]] value set to the same as the current row.
<<csrRowPtrH for empty rows>>=
l = k-1;
 while ( (l>=0) && (!mask[l]) )
    csrRowPtrH[l--] = csrRowPtrH[k];
@
\item cleanup
\index{utilities!sparseMatrix!cleanup}
<<sparse matrix methods>>=
void sparseMatrix::cleanup(void)
{
  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 ) );
}
@
\item matrix vector multiplication $y=Sx$
\index{utilities!sparseMatrix!MVM}
<<sparse matrix methods>>=
void sparseMatrix::MVM(float *y, float *x)
{
HANDLE_ERROR_CUSPARSE(cusparseScsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
                       n_row, n_col, nnz, &alpha, descr,
                       csrValH, csrRowPtrH, csrColIndH,
                       x, &beta, y), "Sparse matrix to vector multiplication failed!");
}
@
\item matrix addition $C = S + B$
\index{utilities!sparseMatrix!add}
<<sparse matrix methods>>=
void sparseMatrix::add(sparseMatrix *C, sparseMatrix *B)
{
 HANDLE_ERROR_CUSPARSE( cusparseScsrgeam( handle, n_row, n_col, &alpha,
                        descr, nnz, csrValH, csrRowPtrH, csrColIndH,
                        &beta, 
                        B->descr, B->nnz, B->csrValH, B->csrRowPtrH, B->csrColIndH,
                        C->descr, C->csrValH, C->csrRowPtrH, C->csrColIndH) ,
                        "Sparse matrix addition failed!");
}
@
\end{itemize}

\subsection{Special functions}
\label{sec:spec_fun}

\subsubsection{Irregular Modified Bessel Functions}
\label{sec:irreg_modif_bessel}

\index{utilities!Knu}
<<Knu temme>>=
__device__ void
  cheb_eval_e(const cheb_series * cs,
              const double x,
              double * result)
{
 int j;
 double d  = 0.0;
 double dd = 0.0;

 double y  = (2.0*x - cs->a - cs->b) / (cs->b - cs->a);
 double y2 = 2.0 * y;

 double e = 0.0;

 for(j = cs->order; j>=1; j--) {
   double temp = d;
   d = y2*d - dd + cs->c[j];
   e += fabs(y2*temp) + fabs(dd) + fabs(cs->c[j]);
   dd = temp;
 }

 { 
   double temp = d;
   d = y*d - dd + 0.5 * cs->c[0];
   e += fabs(y*temp) + fabs(dd) + 0.5 * fabs(cs->c[0]);
 }

 *result = d;
 // result->err = GSL_DBL_EPSILON * e + fabs(cs->c[cs->order]);

 // return GSL_SUCCESS;
}
__host__ __device__ void
temme_gamma(const double nu, double * g_1pnu, double * g_1mnu, double * g1, double * g2)
{
/* nu = (x+1)/4, -1<x<1, 1/(2nu)(1/Gamma[1-nu]-1/Gamma[1+nu]) */
 double g1_dat[14] = {
  -1.14516408366268311786898152867,
   0.00636085311347084238122955495,
   0.00186245193007206848934643657,
   0.000152833085873453507081227824,
   0.000017017464011802038795324732,
  -6.4597502923347254354668326451e-07,
  -5.1819848432519380894104312968e-08,
   4.5189092894858183051123180797e-10,
   3.2433227371020873043666259180e-11,
   6.8309434024947522875432400828e-13,
   2.8353502755172101513119628130e-14,
  -7.9883905769323592875638087541e-16,
  -3.3726677300771949833341213457e-17,
  -3.6586334809210520744054437104e-20
};

/* nu = (x+1)/4, -1<x<1,  1/2 (1/Gamma[1-nu]+1/Gamma[1+nu]) */
 double g2_dat[15] = 
{
  1.882645524949671835019616975350,
 -0.077490658396167518329547945212,  
 -0.018256714847324929419579340950,
  0.0006338030209074895795923971731,
  0.0000762290543508729021194461175,
 -9.5501647561720443519853993526e-07,
 -8.8927268107886351912431512955e-08,
 -1.9521334772319613740511880132e-09,
 -9.4003052735885162111769579771e-11,
  4.6875133849532393179290879101e-12,
  2.2658535746925759582447545145e-13,
 -1.1725509698488015111878735251e-15,
 -7.0441338200245222530843155877e-17,
 -2.4377878310107693650659740228e-18,
 -7.5225243218253901727164675011e-20
};
  cheb_series g1_cs = {
                       g1_dat,
                       13,
                       -1, 1,
                       7
  };
  cheb_series g2_cs = {
                       g2_dat,
                       14,
                       -1, 1,
                       8
  };

 const double anu = fabs(nu);    /* functions are even */
 const double x = 4.0*anu - 1.0;
 double r_g1;
 double r_g2;
 cheb_eval_e(&g1_cs, x, &r_g1);
 cheb_eval_e(&g2_cs, x, &r_g2);
 *g1 = r_g1;
 *g2 = r_g2;
 *g_1mnu = 1.0/(r_g2 + nu * r_g1);
 *g_1pnu = 1.0/(r_g2 - nu * r_g1);
 //return GSL_SUCCESS;
}
__host__ __device__ void
K_scaled_temme(const double nu, const double x,
               double * K_nu, double * K_nup1, double * Kp_nu)
{
  const int max_iter = 15000;

  const double half_x    = 0.5 * x;
  const double ln_half_x = log(half_x);
  const double half_x_nu = exp(nu*ln_half_x);
  const double pi_nu   = M_PI * nu;
  const double sigma   = -nu * ln_half_x;
  const double sinrat  = (fabs(pi_nu) < DBL_EPSILON ? 1.0 : pi_nu/sin(pi_nu));
  const double sinhrat = (fabs(sigma) < DBL_EPSILON ? 1.0 : sinh(sigma)/sigma);
  const double ex = exp(x);

  double sum0, sum1;
  double fk, pk, qk, hk, ck;
  int k = 0;
  //  int stat_iter;

  double g_1pnu, g_1mnu, g1, g2;
  temme_gamma(nu, &g_1pnu, &g_1mnu, &g1, &g2);

  fk = sinrat * (cosh(sigma)*g1 - sinhrat*ln_half_x*g2);
  pk = 0.5/half_x_nu * g_1pnu;
  qk = 0.5*half_x_nu * g_1mnu;
  hk = pk;
  ck = 1.0;
  sum0 = fk;
  sum1 = hk;
  while(k < max_iter) {
    double del0;
    double del1;
    k++;
    fk  = (k*fk + pk + qk)/(k*k-nu*nu);
    ck *= half_x*half_x/k;
    pk /= (k - nu);
    qk /= (k + nu);
    hk  = -k*fk + pk;
    del0 = ck * fk;
    del1 = ck * hk;
    sum0 += del0;
    sum1 += del1;
    if(fabs(del0) < 0.5*fabs(sum0)*DBL_EPSILON) break;
  }
  
  *K_nu   = sum0 * ex;
  *K_nup1 = sum1 * 2.0/x * ex;
  *Kp_nu  = - *K_nup1 + nu/x * *K_nu;

  //stat_iter = ( k == max_iter ? GSL_EMAXITER : GSL_SUCCESS );
  //return GSL_ERROR_SELECT_2(stat_iter, stat_g);
}
__host__ __device__ void
K_scaled_steed_temme_CF2(const double nu, const double x,
                         double * K_nu, double * K_nup1,
                         double * Kp_nu)
{
  const int maxiter = 10000;

  int i = 1;
  double bi = 2.0*(1.0 + x);
  double di = 1.0/bi;
  double delhi = di;
  double hi    = di;

  double qi   = 0.0;
  double qip1 = 1.0;

  double ai = -(0.25 - nu*nu);
  double a1 = ai;
  double ci = -ai;
  double Qi = -ai;

  double s = 1.0 + Qi*delhi;

  for(i=2; i<=maxiter; i++) {
    double dels;
    double tmp;
    ai -= 2.0*(i-1);
    ci  = -ai*ci/i;
    tmp  = (qi - bi*qip1)/ai;
    qi   = qip1;
    qip1 = tmp;
    Qi += ci*qip1;
    bi += 2.0;
    di  = 1.0/(bi + ai*di);
    delhi = (bi*di - 1.0) * delhi;
    hi += delhi;
    dels = Qi*delhi;
    s += dels;
    if(fabs(dels/s) < DBL_EPSILON) break;
  }
  
  hi *= -a1;
  
  *K_nu   = sqrt(M_PI/(2.0*x)) / s;
  *K_nup1 = *K_nu * (nu + x + 0.5 - hi)/x;
  *Kp_nu  = - *K_nup1 + nu/x * *K_nu;
  /*
  if(i == maxiter)
    GSL_ERROR ("error", GSL_EMAXITER);
  else
    return GSL_SUCCESS;
  */
}

<<Knu>>=
__host__ __device__ double _K_nu_(const double nu, const double x)
{
 double result = 0.0;
 int N = (int)(nu + 0.5);
 double mu = nu - N;      /* -1/2 <= mu <= 1/2 */
 double K_mu, K_mup1, Kp_mu;
 double K_nu, K_nup1, K_num1;
 int n, e10 = 0;

 if(x < 2.0) {
   K_scaled_temme(mu, x, &K_mu, &K_mup1, &Kp_mu);
 }
 else {
   K_scaled_steed_temme_CF2(mu, x, &K_mu, &K_mup1, &Kp_mu);
 }

 /* recurse forward to obtain K_num1, K_nu */
 K_nu   = K_mu;
 K_nup1 = K_mup1;

 for(n=0; n<N; n++) {
   K_num1 = K_nu;
   K_nu   = K_nup1;
   if (fabs(K_nu) > SQRT_DBL_MAX) {
     double p = floor(log(fabs(K_nu))/M_LN10);
     double factor = pow(10.0, p);
     K_num1 /= factor;
     K_nu /= factor;
     e10 += p;
   };
   K_nup1 = 2.0*(mu+n+1)/x * K_nu + K_num1;
 }
 result = K_nu*exp(-x);
 return result;
}
@

\subsection{Geometry}
\label{sec:geometry}

The following routine computes the winding number of a polygon specified by the
[[NV]] vertex $([[Vx]],[[Vy]])$ coordinates around a given point $([[Px]],[[Py]])$:
<<polygon winding number>>=
__host__ __device__ int polywind(double Px, double Py, double *Vx, double *Vy, int NV)
{
  int i, np, wind=0;
  double d0,d1,p0,p1,pt0,pt1;
  np = NV;
  pt0 = Px;
  pt1 = Py;
  p0 = Vx[np-1];
  p1 = Vy[np-1];
  for (i=0;i<np; i++) {
    d0 = Vx[i];
    d1 = Vy[i];
    if (p1 <= pt1) {
      if (d1 > pt1 && (p0-pt0)*(d1-pt1)-(p1-pt1)*(d0-pt0) > 0 ) wind ++;
    }
    else {
      if (d1 <= pt1 && (p0-pt0)*(d1-pt1)-(p1-pt1)*(d0-pt0) < 0 ) wind --;
    }
    p0=d0;
    p1=d1;
  }
  return wind;
}
@
The version of polywind for CUDA is given below:
<<polygon winding number>>=
void polywinds(int *W, double *Px, double *Py, int NP, double *Vx, double *Vy, int NV)
{
  dim3 blockDim(256);
  dim3 gridDim(NP/256+1);
  polywind_kern LLL gridDim,blockDim RRR (W,Px,Py,NP,Vx,Vy,NV);
}
@
with the device kernel
<<polygon winding number kernel>>=
__global__ void polywind_kern(int *W, double *Px, double *Py, int NP, double *Vx, double *Vy, int NV)
{
  int i;
  i = blockIdx.x * blockDim.x + threadIdx.x;
  if (i<NP) {
    W[i] = polywind(Px[i],Py[i],Vx,Vy,NV);
  }
}
@