2 Getting started

2.1 Quasar high-level programming concepts

1. Variables: Quasar variables are (by default) weakly-typed, although some mechanisms exist to enforce strong-typing. Variable names and function names are case sensitive.
2. Data types: some of the built-in data types are listed here:\begin_inset Separator latexpar\end_inset
   • scalar: specifies a floating-point number. The used precision depends on the settings of the computation engine.
   • cscalar: specifies a complex-valued scalar number
   • vec: a dense vector (1D array) of arbitrary length (the size is limited by the system resources)
   • mat: a dense matrix (2D array) of arbitrary size (the size is limited by the system resources)
   • cube: a dense cube (3D array) of arbitrary size (the size is limited by the system resources)
   • cvec: a complex valued dense vector
   • cmat: a complex valued dense matrix
   • ccube: a complex valued dense cube
   • string: a string expression
   • cell: a cell matrix object
   • kernel_function: represents a reference to a kernel or device function (see further).
   • function: a reference to a Quasar (user) function or lambda expression
   • object: a user-defined object (see function object())
   Note that specific functions (see further) need to be used to create variables of a given type. There are also some special built-in datatypes: vecx and cvecx, with x=1,...,32 specify a vector of length x.
3. Scalar numbers: scalar numbers can be entered in decimal notation (5.678) as well as in scientific notation (-1.9e-4). Imaginary numbers are defined by adding the suffix j (or i), hence 1+1j or 1-1j represent complex numbers. Non-decimal numbers are also supported: for example, binary numbers 1011011b (suffix b or B), octal numbers 123456o (suffix o or O) and hexadecimal numbers 1Fh or 0ECD3Fh (suffix h or H).
4. Integer numbers: Quasar supports integer numbers (type int). The bit length of the int type depends on the computation engine, but is typically 32-bit. Also integer types with specified bit length exist: int8, uint8, int16, uint16, uint32.
5. % Comments are cool
   However, note that multi-line comments are currently not (yet) supported.
6. Assignment expressions:
   cool = 1
   quasar = cool
   The separation of lines using “;” is optional, and only mandatory when multiple statements are placed on the same line. One can assign to multiple variables at once, similar to C/C++:a = b = 1also, the result of an assignment is a value (in this case, 1). Multiple variable assignment is also possible (i.e. assigning multiple values to multiple variables at once). For example: [a, b] = [1, 2]will assign 1 to a and 2 to b. It is equivalent to:a=1; b=2The multiple variable assignment is mostly useful for 1) assigning multiple return values from functions, for interchanging values:[a, b]=[b, a]will swap the values of a and b. In some cases, it may be useful to neglect a certain return value. This can be done using the placeholder _:
   [a, _] = [1,2]
   [_, b] = [1,2]
7. Arrays (vectors, matrices, cubes etc.) The data type used may depend on the settings of the computation engine (currently, only 32-bit floating point is allowed, for efficiency). The following program illustrates how to create vectors and perform operations:
   a = [0, 1, 2, 3] + 4
   b = [3, 3, 3, 3]
   print "a = ", a, "a .* b = ", a .* b, "sum(a.*b) = ", sum(a .* b)
   print "a^2 = ", a.^2
   String expressions are defined by double quotes (“”), the function “print” allows to print several comma separated values to the console. The “sum” function computes the sum of all components of the vector. The first line evaluates toa = [4, 5, 6, 7]i.e., 4 is added to every component of the vector. [4, 5, 6, 7] then represents a row vector. A matrix can be defined as follows:a = [[1,2],[2,1]]Statements and expressions can be split across multiple code lines. The following is also valid:
   a = [[1,2],
   [2,1]]
   However, for readability, it is adviced to put an underscore at the line break:
   a = [[1,2], _
   [2,1]]
   Similarly, a 3D matrix can be defined by:a = [[[1,2],[3,4]],[[5,6],[7,8]]]An alternative way of defining matrices is using the function zeros(.) or ones(.), which will initialize the values of the matrix to 0 and 1, respectively:
   a = zeros(5)
   b = zeros(6, 4)
   c = zeros(8, 6, 4)
   d = ones(5)
   e = ones(6, 4)
   f = ones(8, 6, 4)
   Alternatively, a vector with the dimensions can be used:
   dims = [4, 5, 6]
   a = ones(dims)
   The function size(.) returns the size of a vector/matrix/cube:
   dims = size(a)
   dim_y = size(a,0)
   dim_x = size(a,1)
   dim_z = size(a,2)
   [dim_y,dim_x] = size(a,0..1)
   [dim_y,dim_x,dim_z] = size(a)
   Note that the dimensions are zero-based. By convention, y is the first dimension (corresponding to index 0), x is the second dimension and z is the third dimension. Internally, matrices are stored in row-major order. [C]  [C] This is in contrast to MATLAB, which uses column-major order (i.e. FORTRAN order). An n × n identity matrix can be created by using the function eye(.):Q = eye(n)Another example:
   a = [[1,2],[2,1]]
   b = [[3],[4]]
   a[0,0] = 2
   print a * b, ",", eye(3)
   print a[0,0]
   print "size(a)=", size(a), "size(a,1)=", size(a,1)
   The function numel(.) gives the number of elements of a vector/matrix/cube. Practically: numel(X) = prod(size(X)).
8. Operators: see Table 2.1↓.\begin_inset Separator latexpar\end_inset
   Notes:
   • There are no bit-wise integer operators. Instead, use the functions and (bitwise conjunction), or (bitwise disjunction), xor (exclusive or), not (bitwise negation), shl (bit-wise left shift), shr (bitwise right shift).
   • Within kernel or device functions, the operators +=, -= have a special meaning: they specify atomic operations (i.e. these operations are free of data races). There are currently 13 atomic operators (see table below).
     Table 2.1 Some operators

     =	assignment	!	inversion (of Boolean values)
     +	add	&&	Boolean AND
     -	subtract / negation	||	Boolean OR
     *	matrix multiplication (or multiplication of scalar values)	? :	Conditional expression (similar to C/C++)
     /	division of scalar values	+=	a += b is a shorthand for a = a + b
     .*	point-wise multiplication (vec, mat, cube data types)	-=	a -= b is a shorthand for a = a - b
     ./	point-wise division (vec, mat, cube data types)	*=	a *= b is a shorthand for a = a * b
     ^	exponentiation (scalar values currently)	/=	a /= b is a shorthand for a = a / b
     .^	point-wise exponentiation	^=	a ^= b is a shorthand for a = a ^ b
     <	smaller than	.*=	a .*= b is a shorthand for a = a .* b
     <=	smaller than or equal	./=	a ./= b is a shorthand for a = a ./ b
     >	greater than	.^=	a .^= b is a shorthand for a = a .^ b
     >=	greater than or equal	^^=	Atomic maximum
     ==	equality	__=	Atomic minimum
     !=	inequality	~=	Atomic bitwise exclusive or (XOR)
     ..	Defines a sequence (see further)	|=	Atomic bitwise OR
     		&=	Atomic bitwise AND

9. Sequences: a sequence defines a row vector:
   a=0..9
   b=0..2..6
   The middle argument defines the step size. Generally, the sequence includes the specified lower and upper bounds. [D]  [D] Except when the step size is too large, such as 0..2..3 = [0, 2] or 0..100..10 = [0]. Hence, the above statements are equivalent to:
   a=[0,1,2,3,4,5,6,7,8,9]
   b=[0,2,4,6]
   The sequences can subsequently be used for matrix indexing:
   A=randn(64,64)
   A_sub = A[a,b]
   Example:
   a = 1..2..10
   b = sum(a)
   c = linspace(1, 2, 5)
   print "a = ", a, "c = ", c
   print sum = ", [b, sum(c)]
   The linspace function creates an uniformly spaced row vector of 5 values between 1 and 2, hence c = [1,1.25,1.5,1.75,2]. Implicit sequences (:) can be used to quickly index matrices:print A[:,0], A[0,:]This statement prints the first column of A, followed by the first row of A. Note that for the Matlab keyword “end”, there is no Quasar equivalent. However, it is still possible to use A[0..size(A,0)-1,0].
10. Control structures: Quasar supports several control structures:
    for a=0..2..4
    break
    continue
    endfor
    if a==2
    endif
    if a==2
    ...
    elseif a==3
    ...
    else
    ...
    endif
    while expr
    break
    continue
    endwhile
    repeat
    break
    continue
    until expr
    Note that “if” is ended with “endif”. Also “if”, “endif” statements must be spread along several lines of code. This is to improve readability of the code. The following is NOT allowed:if a==2; do_something(); endif % Not allowed! An example of a for-loop:
    for i=1..2..100
    j=i+1
    print i, " ", j
    if i==1
    print "i is one"
    elseif i==3
    print "i is three"
    endif
    endfor
    Non-uniform ranges can be specified as follows:
    for powerOfTwo=[1,2,4,8,16,32,64,128]
    print powerOfTwo
    endfor
    or more conveniently as:
    for powerOfTwo=2.^(1..7)
    print powerOfTwo
    endfor
11. Switches are also possible, the syntax is a little different, for example:
    match a with
    | 1 ->
    print "a=1"
    | 2 ->
    print "a=2"
    | (3, 4) ->
    print "a=3 or a=4"
    | "String" ->
    print "a=String"
    | _ -> print
    "a is something else"
    endmatch
    Note that different data types (i.e. strings and scalar numbers) can be mixed. Multiple case values can be specified (grouped by parentheses).
12. Ternary operators: an inline if is also available, just like in C/C++:
    y = condition ? true_value : false_value
    y = (x > T) ? x - T : 0
    When the condition is true, only the value for true is evaluated. Conversely, when the condition is false, only the false-part is executed.
13. Lambda expressions: simply speaking, lambda expressions define inline functions, for example:
    v = (x,y) -> 2*x+y
    u = x -> 2*x
    w = x -> y -> x + y
    z = w(10)
    print v(1,2), " ", z(5)
    a = [[1,2],[2,1]]
    print v(a,a)
    print w(4)(5)
    Note that here, w is a lambda expression that returns another lambda expression (y -> x + y) when evaluated. As such, partial evaluation is possible, e.g., z=w(10) (see further in Section 4.10↓). Lambda expressions can contain several sub-expressions and can be spread over several lines, as follows:
    print_sum = (a, b) -> (sum=a+b;
    print(sum); sum)
    The different expressions are separated using semicolons (’;’). The return value of the lambda expression is always the last expression (in the above example, sum). Using ternary operators, it is fairly simple to define recursive lambda expressions:factorial = x -> x > 0 ? x * factorial(x - 1) : 1
14. Functions: the syntax for functions is different from the syntax for lambda expressions:function [outarg1, ..., outargM] = name (inarg1, ..., inargN)Here there are M output arguments (outarg1, ..., outargM) and N input arguments (inarg1, ..., inargN). “name” is the name of the function. Note that all output arguments must be assigned, otherwise the function call fails. An example:
    function y = do_something(x)
    y = x * 2
    endfunction
    a = [[1,2],[2,1]]
    b = do_something(a)
    print b
    Calling a function with multiple output arguments requires multiple variable assignment:
    function [x, y] = compute(a, b)
    x = a + b
    y = a * b
    endfunction
    [u, v] = compute(2,3)
    print u, " ", v
    Functions can contain inner functions (up to arbitrary nest depths). The inner functions (direct childs, not siblings) can then only be accessed from the outer function. For example:
    function y = colortransform (x : vec3, cname)
    function a = hsv2rgb (c)
    h = floor(c / 60)
    f = frac(c / 60)
    v = 255 * c
    p = v * (1 - c)
    q = v * (1 - f * c)
    t = v * (1 - (1 - f) * c)
    match h with
    | 0 -> a = [v, t, p]
    | 1 -> a = [q, v, p]
    | 2 -> a = [p, v, t]
    | 3 -> a = [p, q, v]
    | 4 -> a = [t, p, v]
    | _ -> a = [v, p, q]
    endmatch
    endfunction
    if cname=="hsv2rgb"
    y = hsv2rgb(x)
    else
    error "the specified color transform ",cname, " is not supported!"
    endif
    endfunction
    Argument types can optionally be specified, as shown above in x : vec3. Quasar will check at compile-time (and run-time) if the arguments are of the correct type, otherwise an error will be raised. The presence or absence of argument types has no further influence on the execution and end result of the program (except when types do not match and an error is generated). However, specifying argument types can help the Quasar optimizer to generate more efficient code.\begin_inset Separator latexpar\end_inset
    Function handles can also be used, for example:

    my_func = colortransform
    print my_func([0.2, 0.2, 0.3])
15. Optional function arguments: functions can have optional arguments. In case an argument is missing, the default value is used. For example:
    function [y, k] = my_func(b : mat, a : scalar = 4)
    print a + b
    y = k = 0
    endfunction
    my_func(2)
    Since my_func is called with one argument, the default value for the second argument will be used (4 in this case).
    Hence, functions can have multiple outputs and optional arguments, whereas lambda expressions can not. Note that the optional function arguments can - on their turn - be expressions and even function calls:
    function [y, k] = my_func(b : mat, a : mat = eye(4))
    function [y, k] = my_func(b : mat, a : mat = A .* B)
    function [y, k] = my_func(b : mat, a : mat = 2 * b)
    Note that by default, variable references (if the name does not correspond to another input argument) refer to the outer context in which the function is defined. They capture the value at the time the function is defined. The variables are defined in the order that they are put as argument. The following would lead to an error:function [y, k] = my_func(a : mat = 2 * b, b : mat) Here, b is not defined at the time a = 2 * b is evaluated.
16. Cell matrices: vectors, matrices and cubes can be grouped in a cell-structure. Cell matrices are either created using the function cell or using the special designated quotes ‘’. copy(.) performs a deep copy of a cell matrix (i.e. the function recursively applies copy(.) to all its elements). Some examples are given below:
    A = cell(2,2)
    G = cell(3)
    A[0,0] = eye(4)
    A[1,1] = 3
    A[0,1] = cell(1,4)
    A[0,1][1] = ones(3,3)
    A[0,1][1][0,0] = 2
    print A[0,1][1]*3
    print size(A)*2
    B = copy(A)
    C = B-A
    print C[0,0]
    D = {A,B,C}
    D_names = {"A","B","C"}
    print D_names[1]
    print size(’’)
    One important special feature is that operations on cell matrices are supported when the different operands have the same structure. It is possible to compute the sum of two cell matrices using:C = B+AAlthernatively, we can multiply all elements of a cell matrix by a constant:C = B*4Or, we can use cell matrices in function calls (note that this is only allowed with built-in functions).C = max(B,4)
    Cell matrices are convenient structures especially for rapid prototyping. However, because cell matrices can store any data, type inference that is required for efficient parallelization may fail on cell matrices. For this purpose, it is useful to use fully typed cell matrices (see Section 3.5↓).
17. Dynamic evaluation: string expressions can be parsed and evaluated at runtime using the eval(.) function:val = eval("(x) -> 3*eye(x)")(8)Here, the eval function parses the string expression ’’(x) -> 3*eye(x)´´ and returns a corresponding lambda expression. This lambda expression can then be evaluated at the same speed as “regular” lambda expressions. This can be useful for simulations (e.g. passing functions through the command line).
18. Reading an input image:img_in = imread("lena_big.tif")Grayscale images return a two-dimensional matrix, color images return a three-dimensional cube, in which the length of the third dimension is either 3 (RGB) or 4 (RGBA - RGB with an alpha channel).
19. The spread operator: the spread operator “...´´ allows to unpack vectors to arbitrary indices or function parameters. Using the spread operator, the following lines of code can be simplified:
    pos = [0,1,2]
    y = im[pos[0],pos[1],pos[2],0] % Before
    y = im[...pos, 0] % After
    luminance = (R,G,B) -> 0.2126 * R + 0.7152 * G + 0.0722 * B
    c = [128, 42, 96]
    lum = luminance(c[0],c[1],c[2]) % Before
    lum = luminance(...c) % After
    The spread operator is in particularly useful in combination with variadic functions (see Section 4.8↓).

2.2  A brief introduction of the type system

isreal      = x -> type(x, "scalar")  || type(x,"vec")  || type(x,"mat")  
                || type(x, "cube")
iscomplex   = x -> type(x, "cscalar") || type(x,"cvec") || type(x,"cmat") 
                || type(x, "ccube")
isscalar    = x -> type(x, "scalar")  || type(x, "cscalar")
isvector    = x -> type(x, "vec")     || type(x, "cvec")
ismatrix    = x -> type(x, "mat")     || type(x, "cmat")  || isvector(x)
iscube      = x -> type(x, "cube")    || type(x, "ccube") || ismatrix(x)

• [int -> int]: a function that takes "int" as input and gives “int” as output
• [(mat, scalar) -> (mat, mat)]: a function that takes two input arguments (of type mat and scalar) and that has two output arguments (both of type mat)
• [int -> int -> int]: a function that projects an integer input onto a lambda expression of type int -> int.
• [...scalar -> scalar]: a variadic function that takes an arbitrary number of scalar values as input and that returns a scalar number.

2.2.1  Floating point representation

• When numerical accuracy is an issue: remark that results obtained using double-precision arithmetic may differ from the same operations obtained using single-precision arithmetic, due to the greater precision and due to rounding errors. Therefore, it is important to compare and express the results within a certain tolerance rather than expecting them to be exact. Moreover, GPU devices typically flush numbers smaller than the minimum representable value (in absolute sense) to zero. Correspondingly, by using double precision FP numbers it may be possible to reduce some of the error introduced by underflow, as the minimal representable value is of the order 10^− 308 for double, while 10^− 38 for single precision (see Table 2.5↑).
• For comparing the results of the algorithms to MATLAB/C++ implementations using double precision FP numbers.

2.2.2  Mixed precision floating point computations

• scalar’single: represents a 32-bit IEEE single precision floating point type (FP32)
• scalar’double: represents a 64-bit IEEE double precision floating point type (FP64)
• scalar’half: represents a 16-bit IEEE half precision floating point type (FP16)

2.2.3  Integer types

• int8: a signed 8-bit integer (with range -128..127)
• int16: a signed 16-bit integer (with range -32768..32767)
• int32: a signed 32-bit integer (with range  − 2³¹..2³¹ − 1)
• int64: (not fully implemented yet)
• uint8: an unsigned 8-bit integer (with range 0..255)
• uint16: an unsigned 16-bit integer (with range 0..65535)
• uint32: an unsigned 32-bit integer (with range 0..2³² − 1)
• uint64: an unsigned 64-bit integer (with range 0..2⁶⁴ − 1)

• int’checked (default): generates an error when the integer can not be represented using the current type (note: not implemented yet)
• int’sat: in case of overflow, the integer is saturated (clipped) to the highest (or lowest) possible value that can be represented.
• int’unchecked: performs no integer overflow checking. This may often be the fastest.

Important note

2.2.4  Fixed sized datatypes

2.2.5  Higher dimensional matrices

2.2.6  User-defined types, type definitions and pointers

• Multiple indirections (^^point) are not allowed.
• All pointer values should be initialized, either used a constructor of the class, or using a null pointer (nullptr). For example, the above class can be initialized using:r = pt_rectangle(p1:=nullptr, p2:=point(1,2))
• Pointers are only allowed to be used for UDTs, not for scalars (scalar, cscalar, ...) or matrices (vec, mat, cube, ...).
• All pointer values are typed. It is for example not allowed to declare a pointer to an unknown type (^??).
• Pointer arithmetic is also not allowed.

2.3  Automatic parallelization

2.4  Writing parallel code using kernel functions

• basic usage (does not require any prior knowledge on GPU programming) - see Subsection 2.4.2↓.
• advanced usage (for “experienced” GPU users) - see Subsection 2.4.4↓.

2.4.1  Basic usage: kernel functions

• vec(X) (or shortcut [F]  [F] replace X by the exact value, for example vec2, vec3, ... vecX): corresponds to a vector of length X.
• ivec(X) (or shortcut ivecX): corresponds to an integer vector of length X.
• cvec(X) (or shortcut ivecX): corresponds to a complex-valued vector of length X.
• mat(X,Y) (or shortcut matXxY): a matrix of size X × Y.
• cube(X,Y,Z)(or shortcut cubeXxYxZ): a cube of size X × Y × Z.
• int: integer data type

• pos (of type int, (i)vecX): the current position of the work item being processed. Note that a “work item” can be either an individual pixel, or a “group of pixels”, depending on how you specify the “dimensions” argument.
• blkpos (of type int, (i)vecX): the current position within the block (for advanced users, see Subsection 2.4.4↓)
• blkidx (of type int, (i)vecX): the block index (for advanced users, see Subsection 2.4.4↓)
• blkdim (of type int, (i)vecX): the current dimensions a block (for advanced users, see Subsection 2.4.4↓). Internally, the data is processed on a block-by-block basis. The dimensions of a block depend on the computation engine in use. For example, for the CUDA computation engine (with CUDA compute capability 2.0), the block dimensions can be as large as 16 × 32 or 32 × 16. For the CPU computation engine, the block dimensions will rather be 1 × #num_processors.
• blkcnt (of type int, (i)vecx): the number of blocks in each dimension (for advanced users, see Subsection 2.4.4↓)
• warpsize (of type int): the warp size of the device (for advanced users, see Subsection 2.4.4↓)

• ’safe: disregards writes outside the data boundaries, reads outside the data boundaries evaluate to zero.
• ’circular: performs circular boundary handling
• ’mirror: mirrors when accessing outside the data boundaries.
• ’clamped: clamps (saturates) to the data boundaries (y[0] = y[-1] = y[-2] = ... and y[N-1] = y[N] = y[N+1] = ...)
• ’unchecked (warning: dangerous usage - your program may crash if not used properly): specifies no bounds checking on the input/output data. In case of access outside the data boundaries, a runtime error may or may not be generated. Specify this modifier in case you are sure your kernel/device function is 100% correct, and when you want to enjoy a modest extra code speedup.
• ’checked: the opposite of ’unchecked: generates an error when indices are out of the data boundaries. Quasar will give information on which matrices are the prime suspect.
• ’const: indicates that the matrix variable is constant and will not be changed inside the kernel function.

• Kernel functions can not have optional arguments.
• A kernel function can not call a “host” function.
• A kernel function can call a “device” function (see Subsection 2.4.2↓)
• A kernel function can call a lambda expression declared with the __device__ attribute (see Subsection 2.4.2↓).
• A kernel function can call a lambda expression declared without the __device__ attribute, on the premise that the lambda expression can be inlined.
• A kernel function can call other kernel functions, through parallel_do (see further in Section 4.4↓). A kernel function cannot directly call another kernel function using a standard function call.

• periodize(x, N): periodizes the input coordinate, i.e. k + a⋅N, with 0 ≤ k < N becomes k. This function is used automatically when the modifier ’circular is specified.
• mirror_ext(x, N): mirrors the input coordinate between [0, N − 1]. This function is used automatically when the modifier ’mirror is specified.
• clamp(x, N): clamps the input coordinates to [0, N]. This function is used automatically when the modifier ’clamped is specified.
• int(x): converts the input argument to integer (using type casting)
• float(x): converts the input argument to floating-point (using type casting)
• shared(dims), shared_zeros(dims), shared[T](dims): the function has a special meaning - allocation of shared memory (see Subsection 2.4.4↓).

2.4.2  Device functions

• A device function cannot call a host function. This is simply because “default” functions are, by default, not natively compiled. However, in many cases, it is possible to convert the host function to a device function, by adding the __device__ modifier.
• A device function cannot call a kernel function directly, nor can a kernel function call another kernel function (unless parallel_do is used, see further in Section 4.4↓). This is because kernel functions have special facilities for parallelization (e.g. they can use OpenMP etc).
• However, a host function can call a device function. This is useful for declaring functions that can be used both from host code as from kernel code. An example is the sinc function:sinc = __device__ (x:scalar) -> x == 0 ? 1.0 : sin(x)/x print sinc(0) % call the device function

2.4.3  Memory usage inside kernel or device functions

1. local memory: this is memory that is local to the function, and each parallel run of the kernel function (called ’thread’) contains a private copy of this memory. Below are a few examples of the creation of local memory:
   A = [0, 1, 2, 3] % Generates a variable of type ’ivec4’
   B = [0., 1., 2., 3.] % Generates a variable of type ’vec4’
   C = [1 + 1j, 2 - 2j] % Generates a complex-valued variable of type ’cvec2’
   D = ones(6) % Generate a vector of length 6, filled with 1.
   E = zeros(8) % Generates a vector of length 8, filled with 0.
   F = complex(zeros(4)) % Generates a complex-valued vector of length 4
   For the GPU computation engine, there are however a few limitations: first, local memory is internally stored in device registers, is hence very fast, but also scarse. When the maximum number of device registers is exceeded, global memory is used instead (which has a much larger memory access latency).
2. shared memory: this type of memory is shared across threads, and allocated using the functions shared and shared_zeros. Its usage is discussed in Subsection 2.4.4↓.
3. global memory: this type of memory is used for storing vectors and matrices with either large dimensions or dimensions that cannot be determined at compile-time. Global memory is also used for dynamically allocated objects (see Section 8.3↓). For example:
   function [] = __kernel__ my_kernel (X : mat, pos : ivec2)
   % X[pos] is stored in global memory
   endfunction
   X = ones(4096,4096)
   parallel_do(size(X), my_kernel)
   Here, X is allocated outside a kernel function. The values of X, in total 4 × 4096 × 4096 bytes (in case of 32-bit floating point), are stored automatically in global memory in a linear way. The following formula is used for translating the 3D index to a linear index:
   index(dim₁,  dim₂,  dim₃) = ( dim₁ Ndims₂ +  dim₂) Ndims₃ +  dim₃
   The ind2pos(size, index) function performs exactly this calculation. Global memory can reside either in CPU memory, GPU memory or both. When calling a kernel function using parallel_do in the GPU computation engine, the global memory will automatically be transferred to the GPU. Because the maximum amount of local memory and shared memory that can be used is limited by the hardware (e.g., not more than 48K), global memory is the only way to pass large amounts of data to a kernel function. The only premise is: a kernel/device function cannot allocate global memory, the memory should best be allocated in advance and passed to the function.
   In some cases, a Quasar program may run out of global GPU memory. In that case, Quasar will automatically transfer a non-frequently used memory buffer back to the CPU. This memory buffer can be later transferred back to the GPU. By this technique, Quasar programs can use all the available memory in the system (both CPU and GPU).
4. texture memory: texture memory is read-only global memory that is internally optimized for spatial access patterns (whereas the global memory is more optimal for linear accesses). In particular, the data layout is optimized for texture sampling using nearest neighbor interpolation or linear interpolation. CUDA uses space filling curves [G]  [G] http://en.wikipedia.org/wiki/Space-filling_curve for optimizing the data layout. See Section 9.3↓ for more information.

2.4.4  Advanced usage: shared memory and synchronization

1. For the CUDA computation engine, the amount of shared memory per block is typically 32 KB. On CUDA architectures, shared memory is on-chip and much faster than other off-chip memory. Consequently the amount of shared memory is limited. Taking into account that a (single precision) floating point value takes 4 bytes, the maximum dimensions of a square block of shared memory that you can allocate are 64 × 64. The more shared memory a kernel uses, the less blocks that can be executed in parallel on the GPU.
2. To obtain maximal performance benefits when using shared memory, it is important to make sure that the compiler can determine statically the amount of memory that will be used by the kernel function. If not, the compiler will assume that the kernel function will take all of the available shared memory on the GPU, which prevents the hardware from processing multiple blocks in parallel. For example, if you request: x = shared(20,3,6)the compiler will reserve 20 × 3 × 6 × 4 bytes = 1440 bytes for the kernel function. However, often the arguments of the function shared are non-constant. In this case you can use assertions (see further in Chapter 5↓):
   assert(M<8 && N<20 && K<4)
   x = shared(M,N,K)
   With this assertion the compiler is able to infer the amount of required shared memory. In this case: 8 × 20 × 4 × 4 bytes = 2560 bytes. The compiler then gives the following message:Information: Calculated an upper bound for the amount of shared memory: 2560 bytes

• For the CUDA computation engine, the maximum block size is limited by the maximum number of threads per block. For CUDA compute capability 2.0, we should have that the maximum number of elements  ≤ 1024. In practice, this number is often even lower, due to the resources used by the kernel (registers, shared memory, ...).
• For optimal performance, the number of elements in each block should be a multiple of the warp size (typically 32). To optimize the memory access pattern it might be even desirable to set the block width as a multiple of 32.
• For the CPU computation engine, the maximum block size is unlimited, unless synchronization (syncthreads) is used. In this case, the block size is limited depending on the number of multi-processors in the system.
• Performance is not necessarily proportional to the number of threads per block. In some cases, the optimal launch configuration has a block size that is as small as 64 or 128.