Como punto de partida, analicemos una implementación secuencial tradicional:

void multiplicarSecuencial(float *M1, float *M2, float *Resultado, int n) {
    for (int fila = 0; fila < n; ++fila) {
        for (int col = 0; col < n; ++col) {
            float acumulado = 0.0f;
            for (int pos = 0; pos < n; ++pos) {
                acumulado += M1[fila * n + pos] * M2[pos * n + col];
            }
            Resultado[fila * n + col] = acumulado;
        }
    }
}

Para transformar este algoritmo en su contraparte paralela con SYCL, el primer paso consiste en incluir los encabezados necesarios y establecer un dispositivo de ejecución:

#include <CL/sycl.hpp>
#include <cmath>

namespace sycl = cl::sycl;

La selección del dispositivo puede adaptarse según disponibilidad. En este caso, se empleará un selector que priorice aceleradores GPU cuando estén disponibles:

sycl::queue colaEjecucion(sycl::gpu_selector_v);

La gestión de memoria en SYCL utiliza buffers que abstraen la transferencia entre host y dispositivo. A continuación se ilustra la preparación de datos:

const int dimension = 512;

sycl::buffer<float, 2> bufM1(sycl::range<2>(dimension, dimension));
sycl::buffer<float, 2> bufM2(sycl::range<2>(dimension, dimension));
sycl::buffer<float, 2> bufSalida(sycl::range<2>(dimension, dimension));

colaEjecucion.submit([&](sycl::handler& manejador) {
    auto accM1 = bufM1.get_access<sycl::access::mode::write>(manejador);
    auto accM2 = bufM2.get_access<sycl::access::mode::write>(manejador);
    
    manejador.parallel_for(sycl::range<2>(dimension, dimension), 
        [=](sycl::id<2> idx) {
            accM1[idx] = static_cast<float>(idx[0] % 16) / 16.0f;
            accM2[idx] = static_cast<float>(idx[1] % 16) / 16.0f;
        });
});

El núcleo computacional paralelo se define mediante parallel_for, donde cada ítem de trabajo calcula un elemento de la matriz resultado. La clave reside en que el bucle interno de acumulación permenece dentro del kernel, mientras que las iteraciones externas se distribuyen espacialmente:

colaEjecucion.submit([&](sycl::handler& manejador) {
    auto accM1 = bufM1.get_access<sycl::access::mode::read>(manejador);
    auto accM2 = bufM2.get_access<sycl::access::mode::read>(manejador);
    auto accRes = bufSalida.get_access<sycl::access::mode::write>(manejador);
    
    manejador.parallel_for<class nucleoMatMul>(
        sycl::range<2>(dimension, dimension),
        [=](sycl::item<2> elemento) {
            int i = elemento[0];
            int j = elemento[1];
            float sumaParcial = 0.0f;
            
            for (int k = 0; k < dimension; ++k) {
                sumaParcial += accM1[i][k] * accM2[k][j];
            }
            
            accRes[i][j] = sumaParcial;
        }
    );
});

Para verificar la corrección del resultado, se puede comparar contra la implementación de referencia. El siguiente fragmento completo integra inicialización, ejecución paralela y validación:

#include <CL/sycl.hpp>
#include <iostream>
#include <vector>
#include <cmath>

void multiplicarSecuencial(const std::vector<float>& M1, 
                           const std::vector<float>& M2,
                           std::vector<float>& Resultado, int n) {
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < n; ++j) {
            float valor = 0.0f;
            for (int k = 0; k < n; ++k) {
                valor += M1[i * n + k] * M2[k * n + j];
            }
            Resultado[i * n + j] = valor;
        }
    }
}

int main() {
    const int N = 256;
    std::vector<float> datosA(N * N), datosB(N * N);
    
    for (int i = 0; i < N * N; ++i) {
        datosA[i] = static_cast<float>(i % 7) / 7.0f;
        datosB[i] = static_cast<float>(i % 5) / 5.0f;
    }
    
    sycl::queue q(sycl::default_selector_v);
    
    sycl::buffer<float, 2> bA(datosA.data(), sycl::range<2>(N, N));
    sycl::buffer<float, 2> bB(datosB.data(), sycl::range<2>(N, N));
    sycl::buffer<float, 2> bC(sycl::range<2>(N, N));
    
    q.submit([&](sycl::handler& h) {
        auto a = bA.get_access<sycl::access::mode::read>(h);
        auto b = bB.get_access<sycl::access::mode::read>(h);
        auto c = bC.get_access<sycl::access::mode::write>(h);
        
        h.parallel_for<class mm_paralelo>(sycl::range<2>(N, N), 
            [=](sycl::item<2> it) {
                int fila = it[0], col = it[1];
                float acum = 0.0f;
                for (int p = 0; p < N; ++p) {
                    acum += a[fila][p] * b[p][col];
                }
                c[fila][col] = acum;
            });
    });
    
    std::vector<float> resultadoParalelo(N * N);
    {
        auto accFinal = bC.get_access<sycl::access::mode::read>();
        for (int i = 0; i < N; ++i)
            for (int j = 0; j < N; ++j)
                resultadoParalelo[i * N + j] = accFinal[i][j];
    }
    
    std::vector<float> resultadoReferencia(N * N);
    multiplicarSecuencial(datosA, datosB, resultadoReferencia, N);
    
    bool coincide = true;
    for (size_t idx = 0; idx < resultadoParalelo.size(); ++idx) {
        if (std::abs(resultadoParalelo[idx] - resultadoReferencia[idx]) > 1e-4f) {
            coincide = false;
            break;
        }
    }
    
    std::cout << (coincide ? "Verificación exitosa" : "Discrepancia detectada") << std::endl;
    
    return 0;
}

El modelo de ejecución de SYCL permite que el runtime de oneAPI optimice automáticamente la distribución de trabajo entre unidades de cómputo, aprovechando la jerarquía de memoria del dispositivo seleccionado. La especificación del rango bidimensional en parallel_for preserva la localidad espacial de los accesos a memoria, factor crítico para el rendimiento en arquitecturas vectoriales.