Умножение одной прямоугольной матрицы на другую - C (СИ)

#include <stdio.h>
#include <stdlib.h>
 
typedef struct MATRIX {
    double ** data;
    size_t rows;
    size_t columns;
} matrix_t;
 
matrix_t * matrix_new(size_t rows, size_t columns) {
    size_t i;
    matrix_t * pMatrix;
        
    if ( ! ( pMatrix = malloc(sizeof(matrix_t)) ) )
        return NULL;
 
    if ( ! ( pMatrix->data = malloc(sizeof(double*) * rows) ) ) {
        free(pMatrix);
        return NULL;
    }
    for ( i = 0; i < rows; ++i ) {
        if ( ! ( pMatrix->data[i] = malloc(sizeof(double) * columns) ) ) {
            while ( i )
                free(pMatrix->data[--i]);
            free(pMatrix->data);
            free(pMatrix);
            return NULL;
        }
    }
    
    pMatrix->rows = rows;
    pMatrix->columns = columns;
    
    return pMatrix;
}
 
matrix_t * matrix_new_from_array(size_t rows, size_t columns, const double * pArray) {
    size_t i, j;
    matrix_t * pMatrix;
    
    if ( ! ( pMatrix = matrix_new(rows, columns) ) )
        return NULL;
    
    for ( i = 0; i < rows; ++i )
        for ( j = 0; j < columns; ++j )
            pMatrix->data[i][j] = *(pArray + i * columns + j);
    
    return pMatrix;
}
 
void matrix_free(matrix_t * pMatrix) {
    size_t i;
    
    for ( i = 0; i < pMatrix->rows; ++i )
        free(pMatrix->data[i]);
    free(pMatrix->data);
    free(pMatrix);
}
 
void matrix_dump(const matrix_t * pMatrix, size_t width, size_t precision) {
    size_t i, j;
    
    for ( i = 0; i < pMatrix->rows; ++i ) {
        for ( j = 0; j < pMatrix->columns; ++j )
            printf("%*.*f", width, precision, pMatrix->data[i][j]);
        printf("\n");
    }
}
 
matrix_t * matrix_mul(const matrix_t * pMatrixA, const matrix_t * pMatrixB) {
    matrix_t * pMatrixC;
    size_t i, j, k;
    
    if ( pMatrixA->columns != pMatrixB->rows )
        return NULL;
    
    if ( ! ( pMatrixC = matrix_new(pMatrixA->rows, pMatrixB->columns) ) )
        return NULL;
    
    for ( i = 0; i < pMatrixC->rows; ++i ) {
        for ( j = 0; j < pMatrixC->columns; ++j ) {
            pMatrixC->data[i][j] = 0.0;
            
            for ( k = 0; k < pMatrixA->columns; ++k )
                pMatrixC->data[i][j] += pMatrixA->data[i][k] * pMatrixB->data[k][j];
        }
    }
    
    return pMatrixC;
}
 
int main(void) {
    double arrA[] = { 1., 2., 3., 4., 5., 6. };
    double arrB[] = { 7., 8., 9., 10., 11., 12. };
    matrix_t * mA = matrix_new_from_array(2, 3, arrA);
    matrix_t * mB = matrix_new_from_array(3, 2, arrB);
    matrix_t * mC = matrix_mul(mA, mB);
    
    printf("First matrix:\n");
    matrix_dump(mA, 4, 0);
    printf("\nSecond matrix:\n");
    matrix_dump(mB, 4, 0);
    
    printf("\nMultipled matrix:\n");
    matrix_dump(mC, 4, 0);
    
    matrix_free(mA);
    matrix_free(mB);
    matrix_free(mC);
    
    return 0;
}