template <unsigned int blockSize>
__global__ void scalar_7(float *Mat, int ind_x, int ind_y, int Len, int Size, float *resMat ){
    size_t s = threadIdx.x + blockIdx.x * blockDim.x;
    int sID = threadIdx.x;
    int i;
    int stride;
    float sum=0;
    int indx;
    int ad_x=ind_x*Len;
    int ad_y=ind_y*Len;
    int loc_x;
    int loc_y;

    stride = Len / blockDim.x +1;


    // We use a single block to process single pair of vectors

    volatile __shared__ float pom[BlockSize];

    pom[sID] = 0;
    // Najpierw każdy wątek czyta 'stride' wartości z wektorów x i y
    // i wykonuje fragment obliczeń do iloczynu skalarnego
    // Potem wszystkie wątki robią redukcję.
    // Pamiętamy, że każdy wektor jest w rzeczywistości fragmentem pamięci w jednowymiarowej
    // tablicy Mat. Dlatego musimy policzyć odpowiedni adres każdego elementu z tej tablicy
    for (i=0;i<stride;i++) {
        indx = s*stride+i; // Wyliczamy dla każdego wątku numer elementu w wektorach x i y
        loc_x = ad_x + indx; // Wyliczamy numery elementów w tablicy Mat odpowiadające
        loc_y = ad_y + indx; // elementom o numerze indx w wektorach x i y
        if (indx<Len) {
            sum+=Mat[loc_x]*Mat[loc_y];
        }
    }
    pom[sID] = sum;
    __syncthreads(); // synchronizujemy wątki - wszystkie wątki muszą skończyć swoje wczytywanie
    // zanim przejdziemy dalej

    if (blockSize >= 1024) {
        if (sID < 512) { pom[sID] += pom[sID + 512];} __syncthreads();
    }

    if (blockSize >= 512) {
        if (sID < 256) { pom[sID] += pom[sID + 256]; } __syncthreads();
    }
    if (blockSize >= 256) {
        if (sID < 128) { pom[sID] += pom[sID + 128]; } __syncthreads();
    }
    if (blockSize >= 128) {
        if (sID < 64) { pom[sID] += pom[sID + 64]; } __syncthreads();
    }

    if(sID < 32) {
        pom[sID] += pom[sID + 32];
        pom[sID] += pom[sID + 16];
        pom[sID] += pom[sID + 8];
        pom[sID] += pom[sID + 4];
        pom[sID] += pom[sID + 2];
        pom[sID] += pom[sID + 1];
    }

    if (sID==0) {
        // wpisujemy odpowiedni wynik do tablicy wynikowej
        resMat[ind_x*Size+ind_y]=pom[0];
    }

}