Proyecto: motor matricial · C++

Un motor matricial convierte una idea matemática en algo que la máquina ejecuta: una estructura que guarda números ordenados y operaciones —crear, sumar, multiplicar, leer— con reglas explícitas. En Python esa pieza es NumPy y la das por hecha. Aquí la escribes, porque es la forma más rápida de entender qué ocurre dentro de np.dot.

El objetivo no es competir con Eigen ni con BLAS: es tener un componente propio, con errores que se entienden, pruebas que corren en un segundo y números que puedes defender. Al final tendrás una biblioteca que el servidor de inferencia de la lección siguiente usará sin cambios.

¿Qué es un motor matricial?

Formalmente, una matriz es una tabla de números con dos índices. En memoria no existen las tablas: existe una tira de bytes. El motor es el conjunto de decisiones que traduce (i, j) a una posición de esa tira y que define qué significa sumar o multiplicar dos de esas tiras. Todo lo demás —rendimiento, seguridad, legibilidad— se deriva de esas decisiones.

La analogía útil es un estacionamiento con filas y columnas numeradas. Si guardas los autos por filas, recorrer una fila es caminar hacia adelante; recorrer una columna es saltar de un extremo a otro. El mismo conjunto de autos, dos formas de caminarlo, dos velocidades distintas.

¿Para qué sirve?

Capas de una red: cada capa densa es y = W * x + b. Sin matrices no hay red.
Transformaciones geométricas y normalización: rotar un conjunto de puntos, estandarizar columnas, aplicar una matriz de covarianza.
Distancias y Kernels: el producto X * Xᵀ aparece en k-NN, SVM y en la atención de un Transformer.
Sistemas de ecuaciones: resolver A x = b para calibrar un modelo pequeño o un ajuste por mínimos cuadrados.
Prototipos embebidos: un dispositivo sin Python necesita una implementación mínima y predecible.

En la lección anterior entrenaste un modelo con bucles y vectores. Cuando el modelo crece a varias features, esos bucles se vuelven matrices. Este proyecto es el paso intermedio entre «un bucle por feature» y «una biblioteca de álgebra lineal»: suficiente para un modelo real, pequeño como para leerlo entero.

Cómo funciona por dentro: almacenamiento contiguo

La decisión central es el orden de filas (row-major), el mismo que usa C y NumPy por defecto: la fila 0 completa, luego la fila 1, y así. El elemento (i, j) vive en la posición i * columnas + j. Una matriz de 3x4 ocupa doce double consecutivos, 96 bytes, y ninguna fila está partida.

Esa fórmula parece un detalle y gobierna el rendimiento. Un double mide 8 bytes y una línea de caché 64: caben ocho. Cuando sumas dos matrices recorriendo por filas, cada línea que traes se usa ocho veces. Cuando las recorres por columnas, cada acceso cae en una línea distinta y pagas un fallo de caché por elemento. El resultado es el mismo; el tiempo puede diferir en un factor de cinco o más.

El producto tiene su propia jerarquía. Multiplicar dos matrices de n x n cuesta n³ multiplicaciones y n³ sumas, es decir 2n³ operaciones en coma flotante. Para n = 512 son 268 millones de operaciones: si tu implementación sostiene 12 GFLOP/s, tarda unos 22 ms. Ese número —y no la elegancia del código— es lo que separa un juguete de una herramienta. Un BLAS multihilo bien ajustado hace lo mismo en bastantes menos milisegundos, y saberlo evita promesas exageradas.

El orden de los bucles decide cuántas veces lees cada dato. Con i-k-j —fila de A, columna común, fila de B— cada a(i, k) se lee una vez y se usa para recorrer entera la fila k de B y acumular en la fila i del resultado. Los tres accesos del bucle interno son contiguos y el compilador puede vectorizarlos. Con i-j-k leerías la columna j de B saltando de fila en fila en cada iteración.

Estructura del proyecto

SHELL árbol de archivos

12345678910

motor/
├── CMakeLists.txt
├── include/
│   └── matrix.hpp      # interfaz publica: lo unico que incluyen los demas
├── src/
│   ├── matrix.cpp      # implementacion
│   ├── main.cpp        # demo ejecutable
│   └── bench.cpp       # medicion de rendimiento
└── tests/
    └── test_matrix.cpp # pruebas sin framework externo

La separación no es decorativa: include/ contiene el contrato —tipos, firmas, invariantes— y src/ cómo se cumple. Quien use la biblioteca incluye matrix.hpp y no necesita saber que dentro hay un std::vector. Si cambias la implementación, nadie más recompila su lógica.

Sintaxis: la interfaz de Matrix

CPP include/matrix.hpp

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657

#pragma once

#include <cstddef>
#include <initializer_list>
#include <iosfwd>
#include <vector>

// Matriz densa de doubles con almacenamiento contiguo en orden de filas.
class Matrix {
public:
    Matrix() = default;

    // Crea filas x columnas celdas inicializadas en `valor`.
    Matrix(std::size_t filas, std::size_t columnas, double valor = 0.0);

    // Crea desde una lista de filas: {{1, 2}, {3, 4}}.
    Matrix(std::initializer_list<std::initializer_list<double>> filas);

    static Matrix identidad(std::size_t n);

    std::size_t filas() const noexcept { return filas_; }
    std::size_t columnas() const noexcept { return columnas_; }
    std::size_t tamano() const noexcept { return datos_.size(); }
    bool vacia() const noexcept { return datos_.empty(); }

    // Acceso al buffer contiguo, sin comprobaciones.
    double* data() noexcept { return datos_.data(); }
    const double* data() const noexcept { return datos_.data(); }

    // Acceso con comprobacion de rango.
    double& operator()(std::size_t i, std::size_t j);
    double operator()(std::size_t i, std::size_t j) const;

    Matrix& operator+=(const Matrix& otra);
    Matrix& operator-=(const Matrix& otra);
    Matrix& operator*=(double escalar);

    bool casi_igual(const Matrix& otra, double tolerancia = 1e-9) const noexcept;

    friend Matrix operator+(Matrix a, const Matrix& b);
    friend Matrix operator-(Matrix a, const Matrix& b);
    friend Matrix operator*(Matrix a, double escalar);
    friend Matrix operator*(double escalar, Matrix a);
    friend Matrix operator*(const Matrix& a, const Matrix& b);
    friend std::ostream& operator<<(std::ostream& os, const Matrix& m);

private:
    std::size_t indice(std::size_t i, std::size_t j) const noexcept {
        return i * columnas_ + j;
    }
    void comprobar_indice(std::size_t i, std::size_t j) const;
    void comprobar_forma(const Matrix& otra, const char* operacion) const;

    std::size_t filas_ = 0;
    std::size_t columnas_ = 0;
    std::vector<double> datos_;
};

Dos constructores: uno por dimensiones con valor inicial, otro desde listas anidadas. El segundo valida que todas las filas midan lo mismo.
operator() en lugar de operator[]: una matriz necesita dos índices, y m(i, j) además acepta la comprobación de rango en un solo sitio.
data() expone el buffer para los bucles críticos, donde comprobar índices por elemento sería absurdo.
casi_igual en lugar de operator==: comparar double con == es una trampa; una tolerancia explícita obliga a decidir qué significa «igual».
friend para los operadores binarios: así a + b y 2.0 * a funcionan con la misma naturalidad que los tipos built-in.

Los constructores son el único lugar donde se valida el tamaño. A partir de ahí, todo objeto Matrix que exista tiene filas y columnas mayores que cero y un buffer del tamaño exacto. Ese invariante es lo que permite que el resto del código no repita comprobaciones: si sostienes la garantía en la frontera, el interior se simplifica.

Ejemplo básico: crear, sumar e imprimir

CPP src/main.cpp

1234567891011121314151617181920212223242526272829

#include <exception>
#include <iostream>

#include "matrix.hpp"

int main() {
    const Matrix a{{1.0, 2.0}, {3.0, 4.0}};
    const Matrix b{{5.0, 6.0}, {7.0, 8.0}};

    std::cout << "A =\n" << a;
    std::cout << "A + B =\n" << a + b;
    std::cout << "A * B =\n" << a * b;
    std::cout << "2 * A =\n" << 2.0 * a;
    std::cout << "A * I =\n" << a * Matrix::identidad(2);

    try {
        const Matrix c(3, 2, 1.0);
        std::cout << a * c;
    } catch (const std::exception& error) {
        std::cout << "error: " << error.what() << '\n';
    }

    try {
        std::cout << a(5, 0);
    } catch (const std::exception& error) {
        std::cout << "error: " << error.what() << '\n';
    }
    return 0;
}

SHELL compilar y ejecutar

12345678910111213141516171819202122

g++ -O2 -std=c++20 -Wall -Wextra -Iinclude \
    src/matrix.cpp src/main.cpp -o demo

./demo

# A =
# 1 2
# 3 4
# A + B =
# 6 8
# 10 12
# A * B =
# 19 22
# 43 50
# 2 * A =
# 2 4
# 6 8
# A * I =
# 1 2
# 3 4
# error: Matrix: A*B requiere columnas(A) == filas(B)
# error: Matrix: indice (5, 0) fuera de rango

Fíjate en las dos últimas líneas: el programa no se cae ni imprime basura, informa qué operación era incompatible y qué índice se salió. Esa diferencia entre un error diagnosticable y un segmentation fault separa una biblioteca usable de un problema ajeno. A * I se construye con Matrix::identidad(2), una función estática que devuelve una matriz nueva.

Ejemplo intermedio: multiplicar sin desperdiciar memoria

CPP src/matrix.cpp

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145

#include "matrix.hpp"

#include <cmath>
#include <ostream>
#include <stdexcept>
#include <string>

namespace {

std::size_t celdas_validas(std::size_t filas, std::size_t columnas) {
    if (filas == 0 || columnas == 0) {
        throw std::invalid_argument("Matrix: las dimensiones deben ser mayores que cero");
    }
    return filas * columnas;
}

}  // namespace

Matrix::Matrix(std::size_t filas, std::size_t columnas, double valor)
    : filas_(filas),
      columnas_(columnas),
      datos_(celdas_validas(filas, columnas), valor) {}

Matrix::Matrix(std::initializer_list<std::initializer_list<double>> filas)
    : filas_(filas.size()), columnas_(filas.size() ? filas.begin()->size() : 0) {
    if (filas_ == 0 || columnas_ == 0) {
        throw std::invalid_argument("Matrix: la lista de filas no puede estar vacia");
    }
    datos_.reserve(filas_ * columnas_);
    for (const auto& fila : filas) {
        if (fila.size() != columnas_) {
            throw std::invalid_argument("Matrix: todas las filas deben tener la misma longitud");
        }
        datos_.insert(datos_.end(), fila.begin(), fila.end());
    }
}

Matrix Matrix::identidad(std::size_t n) {
    Matrix resultado(n, n, 0.0);
    for (std::size_t i = 0; i < n; ++i) resultado(i, i) = 1.0;
    return resultado;
}

void Matrix::comprobar_indice(std::size_t i, std::size_t j) const {
    if (i >= filas_ || j >= columnas_) {
        throw std::out_of_range("Matrix: indice (" + std::to_string(i) + ", " +
                                std::to_string(j) + ") fuera de rango");
    }
}

double& Matrix::operator()(std::size_t i, std::size_t j) {
    comprobar_indice(i, j);
    return datos_[indice(i, j)];
}

double Matrix::operator()(std::size_t i, std::size_t j) const {
    comprobar_indice(i, j);
    return datos_[indice(i, j)];
}

void Matrix::comprobar_forma(const Matrix& otra, const char* operacion) const {
    if (filas_ != otra.filas_ || columnas_ != otra.columnas_) {
        throw std::invalid_argument(std::string("Matrix: dimensiones incompatibles en ") + operacion);
    }
}

Matrix& Matrix::operator+=(const Matrix& otra) {
    comprobar_forma(otra, "A + B");
    for (std::size_t k = 0; k < datos_.size(); ++k) datos_[k] += otra.datos_[k];
    return *this;
}

Matrix& Matrix::operator-=(const Matrix& otra) {
    comprobar_forma(otra, "A - B");
    for (std::size_t k = 0; k < datos_.size(); ++k) datos_[k] -= otra.datos_[k];
    return *this;
}

Matrix& Matrix::operator*=(double escalar) {
    for (double& valor : datos_) valor *= escalar;
    return *this;
}

Matrix operator+(Matrix a, const Matrix& b) {
    a += b;
    return a;
}

Matrix operator-(Matrix a, const Matrix& b) {
    a -= b;
    return a;
}

Matrix operator*(Matrix a, double escalar) {
    a *= escalar;
    return a;
}

Matrix operator*(double escalar, Matrix a) {
    a *= escalar;
    return a;
}

Matrix operator*(const Matrix& a, const Matrix& b) {
    if (a.columnas_ != b.filas_) {
        throw std::invalid_argument("Matrix: A*B requiere columnas(A) == filas(B)");
    }
    Matrix resultado(a.filas_, b.columnas_, 0.0);
    const std::size_t comun = a.columnas_;
    const std::size_t columnas_b = b.columnas_;
    double* out = resultado.data();
    const double* pa = a.data();
    const double* pb = b.data();

    for (std::size_t i = 0; i < a.filas_; ++i) {
        double* fila_resultado = out + i * columnas_b;
        for (std::size_t k = 0; k < comun; ++k) {
            const double aik = pa[i * comun + k];
            const double* fila_b = pb + k * columnas_b;
            for (std::size_t j = 0; j < columnas_b; ++j) {
                fila_resultado[j] += aik * fila_b[j];
            }
        }
    }
    return resultado;
}

bool Matrix::casi_igual(const Matrix& otra, double tolerancia) const noexcept {
    if (filas_ != otra.filas_ || columnas_ != otra.columnas_) return false;
    for (std::size_t k = 0; k < datos_.size(); ++k) {
        if (std::abs(datos_[k] - otra.datos_[k]) > tolerancia) return false;
    }
    return true;
}

std::ostream& operator<<(std::ostream& os, const Matrix& m) {
    for (std::size_t i = 0; i < m.filas_; ++i) {
        for (std::size_t j = 0; j < m.columnas_; ++j) {
            if (j != 0) os << ' ';
            os << m(i, j);
        }
        os << '\n';
    }
    return os;
}

El archivo completo implementa el encabezado. En el producto, tres decisiones merecen atención. Primera: la comprobación de dimensiones ocurre antes de reservar, así que una multiplicación imposible no gasta memoria ni tiempo. Segunda: comun y los punteros se extraen una sola vez, de modo que el triple bucle no comprueba índices por elemento. Tercera: el resultado se inicializa a cero porque el producto acumula; olvidarlo produce números absurdos sin ningún error de compilación.

La suma sigue el mismo patrón: recibe Matrix a por valor —una copia que se va a modificar— y aplica +=. Así a + b + c no crea copias de más y la semántica de valor evita que x = y + z altere y. El compilador puede elidir la copia del parámetro; incluso sin esa optimización el diseño es correcto.

Visualización: el producto celda a celda

Cada celda del resultado es el producto escalar de una fila de A por una columna de B. Recorre los pasos y verifica las cuentas: C(0,0) = 1·5 + 2·7 = 19 y C(1,1) = 3·6 + 4·8 = 50. Si un resultado no coincide con esta tabla, el error está en los índices, no en la aritmética.

Mapa de datos ·C = A · B con A = [[1, 2], [3, 4]] y B = [[5, 6], [7, 8]]

2 × 2 rango [0, 19] C(0,0) = 1·5 + 2·7 = 19. Se acumulan dos productos en la misma celda.

col 0

col 1

fila 0

fila 1

Recorre los pasos para ver cómo cambian las celdas.

Visualización: dónde vive cada número

Una Matrix no guarda sus números dentro de sí misma: guarda tres datos en el stack y un buffer en el heap. Avanza por los pasos y comprueba cómo i * columnas + j convierte un par de índices en una dirección concreta. Esa multiplicación es la que hace que las filas sean contiguas y las columnas no.

Mapa de memoria ·Matrix(2, 3): cabecera en el stack y celdas en el heap

STACK memoria automática

m.filas_ 0x7ffd1000

size_t

m.columnas_ 0x7ffd1008

size_t

m.datos_ 0x7ffd1010

0x6000

vector<double>

↕

HEAP memoria dinámica

m(0,0) 0x6000

0.0

double

m(0,1) 0x6008

0.0

double

m(0,2) 0x6010

0.0

double

m(1,0) 0x6018

0.0

double

m(1,1) 0x6020

0.0

double

m(1,2) 0x6028

0.0

double

m.datos_ ──▶ 0x6000 primer elemento

El constructor valida las dimensiones y reserva 2 · 3 = 6 doubles contiguos, inicializados a 0.0. Son 48 bytes en el heap.

última acción: El constructor valida las dimensiones y reserva 2 · 3 = 6 doubles contiguos, inicializados a 0.0. Son 48 bytes en el heap.

Ver código de esta simulación

Matrix m(2, 3);
m(0, 0) = 1.5;
m(0, 1) = 2.5;
m(1, 2) = 9.0;

Cada paso muestra qué vive en el stack, qué se reserva en el heap y cómo apunta cada variable.

Caso práctico: probar, medir y empaquetar

Una biblioteca sin pruebas es una hipótesis. Estas comprueban propiedades, no valores aislados: A * I debe devolver A sea cual sea A, y lo imposible debe lanzar. No hace falta framework; comprobar y comprobar_lanza bastan, y el programa devuelve un código distinto de cero si algo falla, que es lo único que necesita la integración continua.

CPP tests/test_matrix.cpp

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051

#include <exception>
#include <iostream>
#include <string>

#include "matrix.hpp"

namespace {

int fallos = 0;

void comprobar(bool condicion, const std::string& nombre) {
    if (condicion) {
        std::cout << "ok    " << nombre << '\n';
    } else {
        std::cout << "FALLO " << nombre << '\n';
        ++fallos;
    }
}

template <typename Accion>
void comprobar_lanza(Accion accion, const std::string& nombre) {
    try {
        accion();
        std::cout << "FALLO " << nombre << " (no lanzo excepcion)\n";
        ++fallos;
    } catch (const std::exception&) {
        std::cout << "ok    " << nombre << '\n';
    }
}

}  // namespace

int main() {
    const Matrix a{{1.0, 2.0}, {3.0, 4.0}};

    comprobar(a.filas() == 2 && a.columnas() == 2, "dimensiones 2x2");
    comprobar((a * Matrix::identidad(2)).casi_igual(a), "A * I == A");
    comprobar((a + a).casi_igual(Matrix{{2.0, 4.0}, {6.0, 8.0}}), "A + A");
    comprobar((a * a).casi_igual(Matrix{{7.0, 10.0}, {15.0, 22.0}}), "A * A");
    comprobar((2.0 * a).casi_igual(Matrix{{2.0, 4.0}, {6.0, 8.0}}), "2 * A");

    comprobar_lanza([&a] { [[maybe_unused]] const Matrix r = a * Matrix(3, 1); },
                    "A*B con dimensiones incompatibles");
    comprobar_lanza([&a] { [[maybe_unused]] const double v = a(2, 0); },
                    "indice fuera de rango");
    comprobar_lanza([] { [[maybe_unused]] const Matrix m(0, 3); },
                    "dimensiones cero");

    std::cout << (fallos == 0 ? "todas las pruebas pasaron\n" : "hubo fallos\n");
    return fallos == 0 ? 0 : 1;
}

SHELL compilar y ejecutar las pruebas

1234567891011121314

g++ -O2 -std=c++20 -Wall -Wextra -Iinclude \
    src/matrix.cpp tests/test_matrix.cpp -o run_tests

./run_tests

# ok    dimensiones 2x2
# ok    A * I == A
# ok    A + A
# ok    A * A
# ok    2 * A
# ok    A*B con dimensiones incompatibles
# ok    indice fuera de rango
# ok    dimensiones cero
# todas las pruebas pasaron

Ahora la medición. El benchmark genera dos matrices de 512x512 con una semilla fija —reproducible—, cronometra solo el producto y calcula GFLOP/s a partir de 2n³. Medir con -O2 y sin instrumentación es obligatorio: un binario de depuración puede ser diez veces más lento y llevarte a optimizar lo que no lo necesita.

CPP src/bench.cpp

123456789101112131415161718192021222324252627282930313233343536

#include <chrono>
#include <cstddef>
#include <iostream>
#include <random>

#include "matrix.hpp"

namespace {

Matrix aleatoria(std::size_t n, std::mt19937& rng) {
    std::uniform_real_distribution<double> dist(-1.0, 1.0);
    Matrix m(n, n);
    for (std::size_t k = 0; k < m.tamano(); ++k) m.data()[k] = dist(rng);
    return m;
}

}  // namespace

int main() {
    const std::size_t n = 512;
    std::mt19937 rng(42);
    const Matrix a = aleatoria(n, rng);
    const Matrix b = aleatoria(n, rng);

    const auto inicio = std::chrono::steady_clock::now();
    const Matrix c = a * b;
    const auto fin = std::chrono::steady_clock::now();

    const double ms = std::chrono::duration<double, std::milli>(fin - inicio).count();
    const double flops = 2.0 * static_cast<double>(n) * n * n;
    std::cout << "n = " << n << '\n';
    std::cout << "tiempo = " << ms << " ms\n";
    std::cout << "GFLOP/s = " << flops / (ms * 1e6) << '\n';
    std::cout << "c(0, 0) = " << c(0, 0) << '\n';
    return 0;
}

SHELL medición

123456789

g++ -O2 -std=c++20 -Wall -Wextra -Iinclude \
    src/matrix.cpp src/bench.cpp -o bench

./bench

# n = 512
# tiempo = 21.4 ms
# GFLOP/s = 12.5
# c(0, 0) = -1.80495

Con CMake el proyecto deja de depender de la línea que recordabas: add_library produce la biblioteca, target_include_directories propaga include/ y enable_testing conecta las pruebas con ctest.

SHELL CMakeLists.txt

12345678910111213141516171819202122232425

cmake_minimum_required(VERSION 3.20)
project(motor_matricial LANGUAGES CXX)

set(CMAKE_CXX_STANDARD 20)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)

if(NOT CMAKE_BUILD_TYPE)
  set(CMAKE_BUILD_TYPE Release)
endif()

add_library(matriz src/matrix.cpp)
target_include_directories(matriz PUBLIC include)
target_compile_options(matriz PRIVATE -Wall -Wextra)

add_executable(demo src/main.cpp)
target_link_libraries(demo PRIVATE matriz)

add_executable(bench src/bench.cpp)
target_link_libraries(bench PRIVATE matriz)

enable_testing()
add_executable(tests tests/test_matrix.cpp)
target_link_libraries(tests PRIVATE matriz)
add_test(NAME matriz COMMAND tests)

SHELL configurar, compilar y probar

12345678910111213

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release
cmake --build build -j
ctest --test-dir build --output-on-failure

# -- Configuring done
# -- Generating done
# -- Build files have been written to: .../motor/build
# [100%] Built target demo
#     Start 1: matriz
# 1/1 Test #1: matriz ...........   Passed
# 100% tests passed, 0 tests failed out of 1

./build/bench

Operación	Coste	Acceso a memoria	Riesgo típico
Suma `A + B`	O(n²), 1 suma por celda	Contiguo y predecible	Dimensiones distintas
Escala `k * A`	O(n²)	Contiguo	Confundir `k * A` con `A * B`
Producto `A * B`	O(n³), 2n³ flops	Sensible al orden de bucles	Acumulador sin inicializar
Traspuesta	O(n²)	Un recorrido salta de fila	Escribir `t(i, j) = m(i, j)`
Acceso `m(i, j)`	O(1)	Aritmética de índices	Intercambiar fila y columna

El mismo producto en los dos lenguajes

Python

PYTHON producto.py

12345678

def producto(a, b):
    n, m, p = len(a), len(b), len(b[0])
    c = [[0.0] * p for _ in range(n)]
    for i in range(n):
        for k in range(m):
            for j in range(p):
                c[i][j] += a[i][k] * b[k][j]
    return c

C++

CPP producto.cpp

123456789

Matrix resultado(a.filas(), b.columnas(), 0.0);
for (std::size_t i = 0; i < a.filas(); ++i) {
    for (std::size_t k = 0; k < a.columnas(); ++k) {
        const double aik = a.data()[i * a.columnas() + k];
        for (std::size_t j = 0; j < b.columnas(); ++j) {
            resultado.data()[i * b.columnas() + j] += aik * b.data()[k * b.columnas() + j];
        }
    }
}

La versión de Python crea una lista por fila y guarda punteros a objetos `float`: cada suma implica desreferenciar y crear un objeto. La de C++ opera sobre `double` contiguos, y el compilador puede emitir instrucciones SIMD. Además, con NumPy el mismo producto se resuelve en una llamada a BLAS: ahí C++ deja de tener ventaja frente a Python.

Errores comunes

Comprobar dimensiones al principio, no dentro del bucle

Incorrecto

PYTHON Incorrecto: sin comprobación previa

Matrix operator*(const Matrix& a, const Matrix& b) {
    Matrix r(a.filas(), b.columnas(), 0.0);
    for (std::size_t k = 0; k < a.columnas(); ++k)
        r(0, 0) += a(0, k) * b(k, 0);   // puede salirse de b
    return r;
}

Correcto

PYTHON Correcto: fallar antes de calcular

Matrix operator*(const Matrix& a, const Matrix& b) {
    if (a.columnas() != b.filas())
        throw std::invalid_argument("Matrix: A*B requiere columnas(A) == filas(B)");
    Matrix r(a.filas(), b.columnas(), 0.0);
    // ... bucles seguros
    return r;
}

Sin la comprobación, el error aparece como `indice (k, j) fuera de rango` en mitad del cálculo, o peor, como lectura de memoria válida pero ajena si usas `data()`. Comprobar al entrar convierte un fallo difuso en un mensaje que dice exactamente qué operación y qué regla se rompió.

Intercambiar fila y columna: m(j, i) devuelve un número, no un error, y arruina el resultado sin avisar.
No inicializar el acumulador: Matrix r(n, m) ya inicializa a 0.0, pero si reutilizas una matriz existente con += estarás sumando al valor anterior.
Usar == con double: 0.1 + 0.2 != 0.3. Por eso la clase expone casi_igual con tolerancia explícita.
Desbordar el cálculo del tamaño: int n = 100000; n * n * n desborda. Usa std::size_t y sé consciente de que n * columnas también puede desbordar en matrices absurdamente grandes.
Devolver referencias a objetos locales: const Matrix& f() { Matrix m(2, 2); return m; } compila con un aviso y produce basura. Devuelve por valor: la elisión de copia es la norma.
Compilar sin -O2: el benchmark medirá el código de depuración y concluirás que C++ es lento.

Buenas prácticas

Mantén un único invariante —dimensiones mayores que cero y buffer del tamaño exacto— y valídalo solo en los constructores.
Usa operator() para código legible y data() en los bucles internos; documenta que data() no comprueba nada.
Devuelve por valor y apóyate en la elisión de copia; evita new y delete dentro de la clase.
Escribe pruebas de propiedades (A * I == A) además de casos concretos: detectan errores que un ejemplo no cubre.
Compila siempre con -Wall -Wextra y trata los avisos como errores en el código nuevo.
Mide con Release y con datos de tamaño realista; 2x2 no dice nada sobre 512x512.
Cuando necesites velocidad de verdad, cambia a Eigen o a BLAS antes que reescribir tu propio producto bloqueado a mano.

Ejercicios

Ejercicio 1

Traspuesta y una identidad algebraica

Añade una función Matrix traspuesta(const Matrix& m) que devuelva una matriz de dimensiones intercambiadas. Comprueba con casi_igual que (A * B)ᵀ = Bᵀ * Aᵀ para una A de 2x3 y una B de 3x2. Verifica también que trasponer dos veces devuelve la original.

Pista

La traspuesta de una matriz de f filas y c columnas tiene c filas y f columnas. Recorre la original por filas y escribe en t(j, i).

Ejercicio 2

Producto matriz-vector con validación

Implementa std::vector<double> producto_matriz_vector(const Matrix& m, const std::vector<double>& v). Debe lanzar std::invalid_argument si m.columnas() != v.size() y devolver un vector de tamaño m.filas(). Pruébalo con A = [[1, 2, 3], [4, 5, 6]] y v = [1, 0, -1], cuyo resultado es [-2, -2].

Pista

Reserva el vector de salida con std::vector<double> resultado(m.filas(), 0.0) y acumula por filas; así recorres la matriz de forma contigua.

Completar código Completa el indexado y el acumulador

// m(i, j) vive en esta posicion del buffer contiguo
std::size_t k = i *  + j;

// producto: la celda se acumula, no se asigna
resultado(i, j)  a(i, k) * b(k, j);

El primer hueco es el número de columnas de la matriz. El segundo operador suma al valor existente.

Soluciones explicadas

La traspuesta enseña algo que los libros repiten y el código demuestra: transponer no es «dar la vuelta» a los números, es cambiar el mapeo de índices. Si tu almacenamiento es por filas, la traspuesta obliga a escribir con un salto de filas elementos entre destinos consecutivos; por eso las bibliotecas serias no trasponen físicamente cuando pueden evitarse y en su lugar llevan un stride en la estructura.

El producto matriz-vector ilustra la diferencia entre una operación que se lee de memoria una vez y otra que se relee. Aquí cada elemento de m se usa exactamente una vez y cada v[j] se reutiliza en todas las filas: el patrón de acceso es casi óptimo. Con esos números, la implementación tarda microsegundos incluso para matrices grandes, y envolverla en una llamada a BLAS solo añadiría complejidad. La regla práctica: mide la operación concreta, no la categoría.

Resumen

Un motor matricial es una tira contigua de double más una regla que traduce (i, j) a una posición.
En orden de filas, k = i * columnas + j; las filas son contiguas y las columnas no.
El producto cuesta 2n³ operaciones; el orden i-k-j mantiene los tres accesos en memoria contigua.
Valida dimensiones en el constructor y en las operaciones binarias, antes de reservar o calcular.
Devuelve por valor y usa casi_igual para comparar double; nunca operator==.
Prueba propiedades algebraicas, no solo valores: A * I == A cubre infinitos casos.
Compila con -O2 -Wall -Wextra, ejecuta ctest y mide con datos del tamaño real.
Un motor propio de decenas de GFLOP/s es suficiente para servir un modelo pequeño; para lo demás, Eigen o BLAS.

Evaluación

Evaluación 1

En una Matrix de 3 filas y 4 columnas almacenada por filas, ¿en qué posición vive m(2, 1)?

Evaluación 2

¿Por qué la multiplicación recorre los bucles en orden i-k-j y no i-j-k?

Evaluación 3

Matrix a(2, 2); Matrix b(3, 3); auto c = a * b; ¿Qué ocurre?

Referencia

cppreference · Sobrecarga de operadores reglas y convenciones Eigen · Documentación la alternativa madura para producción CMake · Tutorial bibliotecas, tests y directorios BLAS · Referencia qué hace un producto optimizado por dentro

¿Qué es un motor matricial?

¿Para qué sirve?

Cómo funciona por dentro: almacenamiento contiguo

Estructura del proyecto

Sintaxis: la interfaz de Matrix

Ejemplo básico: crear, sumar e imprimir

Ejemplo intermedio: multiplicar sin desperdiciar memoria

Visualización: el producto celda a celda

Visualización: dónde vive cada número

Caso práctico: probar, medir y empaquetar

Errores comunes

Buenas prácticas

Ejercicios

Traspuesta y una identidad algebraica

Producto matriz-vector con validación

Soluciones explicadas

Resumen

Evaluación

¿Ya lo entendiste?

Más en «Proyectos reales»