```python
import torch
import torch.nn as nn
from torch.utils.cpp_extension import load_inline

_cpp_source = r"""
#include <torch/extension.h>
torch::Tensor triangular_matmul_cuda(torch::Tensor A, torch::Tensor B);
"""

_cuda_source = r"""
#include <torch/extension.h>
#include <ATen/cuda/CUDAContext.h>
#include <c10/cuda/CUDAGuard.h>
#include <c10/cuda/CUDAException.h>
#include <cuda_runtime.h>

constexpr int BM = 64;
constexpr int BN = 128;
constexpr int BK = 16;

__global__ __launch_bounds__(256)
void triangular_matmul_kernel(
    const float* __restrict__ A,
    const float* __restrict__ B,
    float* __restrict__ C,
    int N)
{
    const int tid = threadIdx.x;
    const int ty = tid / 16;
    const int tx = tid % 16;
    const int row_base = blockIdx.y * BM;
    const int col_base = blockIdx.x * BN;

    if (row_base + BM <= col_base) {
        for (int index = tid; index < BM * BN; index += 256) {
            const int row = row_base + index / BN;
            const int col = col_base + index % BN;
            if (row < N && col < N)
                C[(int64_t)row * N + col] = 0.0f;
        }
        return;
    }

    __shared__ float As[BM][BK + 1];
    __shared__ float Bs[BK][BN];

    float acc[4][8] = {};
    const int k_end = min(row_base + BM, N);

    // For lower triangular inputs, only col <= k <= row contributes.
    for (int k_base = col_base; k_base < k_end; k_base += BK) {
        #pragma unroll
        for (int index = tid; index < BM * BK; index += 256) {
            const int r = index / BK;
            const int k = index % BK;
            const int global_r = row_base + r;
            const int global_k = k_base + k;
            As[r][k] = (global_r < N && global_k < N)
                ? A[(int64_t)global_r * N + global_k] : 0.0f;
        }

        #pragma unroll
        for (int index = tid; index < BK * BN; index += 256) {
            const int k = index / BN;
            const int c = index % BN;
            const int global_k = k_base + k;
            const int global_c = col_base + c;
            Bs[k][c] = (global_k < N && global_c < N)
                ? B[(int64_t)global_k * N + global_c] : 0.0f;
        }

        __syncthreads();

        #pragma unroll
        for (int k = 0; k < BK; ++k) {
            float a[4];
            float b[8];

            #pragma unroll
            for (int i = 0; i < 4; ++i)
                a[i] = As[ty * 4 + i][k];

            #pragma unroll
            for (int j = 0; j < 8; ++j)
                b[j] = Bs[k][tx + j * 16];

            #pragma unroll
            for (int i = 0; i < 4; ++i) {
                #pragma unroll
                for (int j = 0; j < 8; ++j)
                    acc[i][j] = fmaf(a[i], b[j], acc[i][j]);
            }
        }

        __syncthreads();
    }

    #pragma unroll
    for (int i = 0; i < 4; ++i) {
        const int row = row_base + ty * 4 + i;
        #pragma unroll
        for (int j = 0; j < 8; ++j) {
            const int col = col_base + tx + j * 16;
            if (row < N && col < N) {
                C[(int64_t)row * N + col] =
                    row >= col ? acc[i][j] : 0.0f;
            }
        }
    }
}

torch::Tensor triangular_matmul_cuda(torch::Tensor A, torch::Tensor B)
{
    TORCH_CHECK(A.is_cuda() && B.is_cuda(), "Inputs must be CUDA tensors");
    TORCH_CHECK(A.device() == B.device(), "Inputs must share a device");
    TORCH_CHECK(A.scalar_type() == at::kFloat &&
                B.scalar_type() == at::kFloat, "Inputs must be float32");
    TORCH_CHECK(A.dim() == 2 && B.dim() == 2 &&
                A.size(0) == A.size(1) && A.sizes() == B.sizes(),
                "Inputs must be square matrices of equal shape");
    TORCH_CHECK(A.is_contiguous() && B.is_contiguous(),
                "Inputs must be contiguous");

    const c10::cuda::CUDAGuard guard(A.device());
    auto C = torch::empty_like(A);
    const int N = static_cast<int>(A.size(0));

    if (N > 0) {
        const dim3 grid((N + BN - 1) / BN, (N + BM - 1) / BM);
        triangular_matmul_kernel<<<
            grid, 256, 0, at::cuda::getCurrentCUDAStream()
        >>>(
            A.data_ptr<float>(),
            B.data_ptr<float>(),
            C.data_ptr<float>(),
            N
        );
        C10_CUDA_KERNEL_LAUNCH_CHECK();
    }
    return C;
}
"""

_triangular_matmul = load_inline(
    name="triangular_matmul_tiled_64x128",
    cpp_sources=_cpp_source,
    cuda_sources=_cuda_source,
    functions=["triangular_matmul_cuda"],
    extra_cflags=["-O3"],
    extra_cuda_cflags=["-O3"],
    verbose=False,
)


class ModelNew(nn.Module):
    def __init__(self):
        super().__init__()

    def forward(self, A, B):
        if (
            not A.is_cuda
            or not B.is_cuda
            or A.dtype != torch.float32
            or B.dtype != torch.float32
        ):
            return torch.tril(torch.matmul(A, B))

        return _triangular_matmul.triangular_matmul_cuda(
            A.contiguous(), B.contiguous()
        )
```