This paper introduces a novel framework for Gaussian process (GP) inference called Blackbox Matrix-Matrix (BBMM) inference, which leverages GPU acceleration to significantly improve the efficiency of GP inference. BBMM uses a modified batched conjugate gradients (mBCG) algorithm to perform all necessary computations for training and inference in a single call. This approach reduces the asymptotic complexity of exact GP inference from O(n³) to O(n²), making it much faster and more efficient than traditional methods like Cholesky decomposition. BBMM also incorporates a specialized preconditioner based on the pivoted Cholesky decomposition to further accelerate convergence. The framework is implemented in GPyTorch, a software platform built on PyTorch that enables scalable GP inference via BBMM. Experiments show that BBMM achieves substantial speedups on both exact GP inference and scalable approximations, with exact GPs being up to 20 times faster than Cholesky-based approaches on large datasets. Additionally, BBMM is shown to be effective for complex models like multi-output GPs and scalable GP approximations, requiring only efficient matrix-matrix multiplication routines. The paper also discusses the benefits of preconditioning, demonstrating that it significantly improves the efficiency of the BBMM approach. Overall, BBMM provides a more efficient and scalable solution for GP inference, particularly on large datasets and with GPU acceleration.This paper introduces a novel framework for Gaussian process (GP) inference called Blackbox Matrix-Matrix (BBMM) inference, which leverages GPU acceleration to significantly improve the efficiency of GP inference. BBMM uses a modified batched conjugate gradients (mBCG) algorithm to perform all necessary computations for training and inference in a single call. This approach reduces the asymptotic complexity of exact GP inference from O(n³) to O(n²), making it much faster and more efficient than traditional methods like Cholesky decomposition. BBMM also incorporates a specialized preconditioner based on the pivoted Cholesky decomposition to further accelerate convergence. The framework is implemented in GPyTorch, a software platform built on PyTorch that enables scalable GP inference via BBMM. Experiments show that BBMM achieves substantial speedups on both exact GP inference and scalable approximations, with exact GPs being up to 20 times faster than Cholesky-based approaches on large datasets. Additionally, BBMM is shown to be effective for complex models like multi-output GPs and scalable GP approximations, requiring only efficient matrix-matrix multiplication routines. The paper also discusses the benefits of preconditioning, demonstrating that it significantly improves the efficiency of the BBMM approach. Overall, BBMM provides a more efficient and scalable solution for GP inference, particularly on large datasets and with GPU acceleration.