dgemmtest.c

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

void dtrmm_(char *side, char *uplo, char *transa, char *diag, int *m, int *n, double *alpha, double *A, int *LDA, double *B, int *LDB);

double *mult(double *A, double *B, int n)
{
    double *X = (double *)calloc(n*n, sizeof(double));
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
            {
                X[i+j*n] += A[i+k*n] * B[k+j*n];
            }
    return X;
}

double *sub(double *A, double *B, int n)
{
    double *X = (double *)calloc(n*n, sizeof(double));
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            X[i+j*n] = A[i+j*n] - B[i+j*n];
    return X;
}

void print(double* A, int n)
{
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
            printf("%f ", A[i+j*n]);
        printf("\n");
    }
}

int main()
{
    int n = 10;
    double *A = (double *)malloc(sizeof(double)*n*n);
    double *B = (double *)malloc(sizeof(double)*n*n);
    double *Bcopy = (double *)malloc(sizeof(double)*n*n);
    double *U = (double *)calloc(n*n, sizeof(double));

    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
        {
            A[i+j*n] = 1.0;
            B[i+j*n] = 1.0;
            Bcopy[i+j*n] = B[i+j*n];
            if (i <= j)
              U[i+j*n] = A[i+j*n];
        }

    print(B, n);
    print(A, n);

    double zero = 0;
    double one = 1;
    dtrmm_("Right", "U", "No transpose", "Non-unit", &n, &n, &one, A, &n, B, &n);

    print(A,n);

    double *R = sub(B, mult(Bcopy, U, n), n);
    printf("\nDTRMM residual\n\n");
    print(R, n);
    return 0;
}