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

int dpotrf_(char *uplo, int *n, double *a, int *lda, int *info);

int dpotri_(char *uplo, int *n, double *a, int *lda, int *info);

void dsyrk_(char *uplo, char *trans, int *N, int *K, double *alpha,
            const double *A, int *lda, double *beta, double *C,
            int *ldc);

int invert_pd_matrix (double *A, int n, int lda)
{
    int info = 0;
    char uplo = 'L';

    dpotrf_(&uplo, &n, A, &lda, &info);

    if (info == 0) {
        dpotri_(&uplo, &n, A, &n, &info);
    }

    return info;
}

void make_pd_matrix (double *A, double *M, int n)
{
    char uplo = 'L';
    char trans = 'T';
    double alpha = 1.0;
    double beta = 0;
    int N = n*n;
    int i;

    for (i=0; i<N; i++) {
        A[i] = rand() / (double) RAND_MAX;
    }
    dsyrk_(&uplo, &trans, &n, &n, &alpha, A, &n, &beta, M, &n);
}

#define NDIM 6

void print_invpd_times (int n[], double T[][NDIM])
{
    int i;

    puts("\nTimes in seconds plus ratio default/single:");
    puts("\n  n   default    single     d/s");
    for (i=0; i<NDIM; i++) {
        printf("%3d %9.4f %9.4f %7.2f\n", n[i], T[0][i], T[1][i],
               T[0][i] / T[1][i]);
    }
    putchar('\n');
}

int record_invpd_times (const double *M, double *A,
                        int lda, int K, int n[],
                        double T[][NDIM], int Tcol)
{
    double t0, t1;
    int N = lda*lda;
    int i, j, err = 0;
    
    for (i=0; i<NDIM && !err; i++) {
        t0 = omp_get_wtime();
        for (j=0; j<K; j++) {
            memcpy(A, M, N * sizeof *A);
            err = invert_pd_matrix(A, n[i], lda);
        }
        t1 = omp_get_wtime();
        T[Tcol][i] = t1 - t0;
    }

    return err;
}

int main (void)
{
    double *M, *A;
    int K = 50000;
    int lda = 64;
    int N = lda*lda;
    int n[NDIM] = {4, 8, 16, 17, 32, 64};
    double T[2][NDIM];
    int err = 0;

    A = malloc(N * sizeof *A);
    M = malloc(N * sizeof *M);
    make_pd_matrix(A, M, lda);

    puts("Using default threading...");
    err = record_invpd_times(M, A, lda, K, n, T, 0);
    if (err) {
        fprintf(stderr, "error %d in inversion (1)\n", err);
    }

    omp_set_num_threads(1);

    puts("Forcing single-threaded behavior...");
    err = record_invpd_times(M, A, lda, K, n, T, 1);
    if (err) {
        fprintf(stderr, "error %d in inversion (2)\n", err);
    }

    print_invpd_times(n, T);

    free(M);
    free(A);

    return 0;
}