#include <cblas.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char **argv)
{
    const double a[4] = {1, 2, 3, 4};
    const double b[4] = {5, 6, 7, 8};
    const double expected[4] = {19, 22, 43, 50};
    double result[4] = {0};
    int threads = openblas_get_num_threads64_();
    printf("%s\n", openblas_get_config64_());
    printf("Threading mode: %d; configured threads: %d; BLAS integer bytes: %zu\n",
           openblas_get_parallel64_(), threads, sizeof(blasint));
    if (argc == 2 && threads != atoi(argv[1])) return EXIT_FAILURE;
    cblas_dgemm64_(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                2, 2, 2, 1.0, a, 2, b, 2, 0.0, result, 2);
    for (int i = 0; i < 4; ++i)
        if (!(fabs(result[i] - expected[i]) < 1e-12)) return EXIT_FAILURE;
    printf("Product: [[%.0f, %.0f], [%.0f, %.0f]]\n",
           result[0], result[1], result[2], result[3]);

    const blasint n = 256;
    size_t count = (size_t)n * (size_t)n;
    double *left = malloc(count * sizeof(double));
    double *right = malloc(count * sizeof(double));
    double *product = calloc(count, sizeof(double));
    if (!left || !right || !product) return EXIT_FAILURE;
    for (size_t i = 0; i < count; ++i) left[i] = right[i] = 1.0;
    cblas_dgemm64_(CblasRowMajor, CblasNoTrans, CblasNoTrans,
                n, n, n, 1.0, left, n, right, n, 0.0, product, n);
    for (size_t i = 0; i < count; ++i)
        if (!(fabs(product[i] - (double)n) < 1e-10)) return EXIT_FAILURE;
    free(product);
    free(right);
    free(left);
    printf("PASS: matrix products with %d configured threads\n", threads);
    return EXIT_SUCCESS;
}