lapack

#include <stdlib.h>
#include <stdio.h>
#include "cpuidh.h"

void dgesv_(int *n, int *nrhs, double *a, int *lda, int *ipiv, double *b, int *ldb, int *info );
void sgesv_(int *n, int *nrhs, float *a, int *lda, int *ipiv, float *b, int *ldb, int *info );
int dgemm_(char *transa, char *transb, int *m, int *n,
    int *k, double *alpha, double *a, int *lda,
    double *b, int *ldb, double *beta, double *c,
    int *ldc);
int sgemm_(char *transa, char *transb, int *m, int *n,
    int *k, float *alpha, float *a, int *lda,
    float *b, int *ldb, float *beta, float *c,
    int *ldc);

#ifdef DP
#define REAL double
#define xgesv dgesv_
#define xgemm dgemm_
#else
#define REAL float
#define xgesv sgesv_
#define xgemm sgemm_
#endif

#define N 200
#define M (N + N/2)
#define K (N - N/2)
#define NN 2000
#define LDAA 4000
#define MOPS ((2.0*N*N*N/3.0 + 2.0*N*N)/1000000)
#define LDA 200
#define TIMES 1000
#define PASS 5

void matgen (REAL a[], int lda, int n, REAL b[], REAL *norma)
{
        int init, i, j;

        init = 1325;
        *norma = 0.0;
        for (j = 0; j < n; j++) {
                for (i = 0; i < n; i++) {
                        init = 3125*init % 65536;
                        a[lda*j+i] = (init - 32768.0)/16384.0;
                        *norma = (a[lda*j+i] > *norma) ? a[lda*j+i] : *norma;

                        /* alternative for some compilers
                        if (fabs(a[lda*j+i]) > *norma) *norma = fabs(a[lda*j+i]);
                        */
                }
        }
        for (i = 0; i < n; i++) {
          b[i] = 0.0;
        }
        for (j = 0; j < n; j++) {
                for (i = 0; i < n; i++) {
                        b[i] = b[i] + a[lda*j+i];
                }
        }
        return;
}

void rmatgen (REAL a[], int m, int n, REAL *norma)
{
        int init, i, j;

        init = 1325;
        *norma = 0.0;
        for (j = 0; j < m; j++) {
                for (i = 0; i < n; i++) {
                        init = 3125*init % 65536;
                        a[n*j+i] = (init - 32768.0)/16384.0;
                        *norma = (a[n*j+i] > *norma) ? a[n*j+i] : *norma;

                        /* alternative for some compilers
                        if (fabs(a[lda*j+i]) > *norma) *norma = fabs(a[lda*j+i]);
                        */
                }
        }
        return;
}

int main (int argc, char * argv[])
{
    static REAL a[LDA * N], b[LDA], norma;
    static int ipiv[N];
    static REAL aa[LDAA * NN], bb[LDAA];
    static int ipivx[NN];
    static REAL mma[M * K], mmb[K * N], mmc[M * N];
    int i, j, m, n, k, nrhs, ntimes, info, lda, ldb;
    double gentime;
    REAL alpha, beta;

    // TEST

    n = N;
    lda = LDA;
    ldb = LDA;
    nrhs = 1;
    matgen(a, LDA, N, b, &norma);
    start_time();
    xgesv(&n, &nrhs, a, &lda, ipiv, b, &ldb, &info);
    end_time();
#ifdef DP
    printf ("DGESV returned %d (time %lf, %lf MFLOPS)\n", info, secs, MOPS/secs);
#else
    printf ("SGESV returned %d (time %lf, %lf MFLOPS)\n", info, secs, MOPS/secs);
#endif
    // matgen
    printf ("matgen overhead calc\n");

    start_time();
    for (i=0; i < TIMES * 10; i++)
        matgen(a, LDA, N, b, &norma);
    end_time();
    gentime = secs/10;

    printf ("matgen time is %lf (%lf)\n", gentime/TIMES, gentime);

    // PASS passes

    for (i=0; i<PASS; i++)
    {
        start_time();
        for (j=0; j< TIMES; j++){
            matgen(a, LDA, N, b, &norma);
            xgesv(&n, &nrhs, a, &lda, ipiv, b, &ldb, &info);
        }
        end_time();
        printf ("n = %d, lda = %d, ldb = %d, nrhs = %d, info = %d\n", n, lda, ldb, nrhs, info);
        printf ("%d. %d xGESV's took %lf sec, %lf MFLOPS\n",
            i, TIMES, (secs-gentime), TIMES * MOPS/(secs - gentime));
    }

    // FAT
    matgen(aa, LDAA, NN, bb, &norma);
    n = NN;
    lda = ldb = LDAA;
    start_time();
    xgesv(&n, &nrhs, aa, &lda, ipivx, bb, &ldb, &info);
    end_time();
    printf ("n = %d, lda = %d, ldb = %d, nrhs = %d, info = %d\n", n, lda, ldb, nrhs, info);
    printf ("This xGESV's took %lf sec, %lf MFLOPS\n", secs, TIMES * MOPS/secs);

    // xGEMM

    printf ("Timing rmatgen\n");
    start_time();
    for (i=0;i<TIMES*10;i++){
        rmatgen(mma, M, K, &norma);
        rmatgen(mmb, K, N, &norma);
    }
    end_time();
    gentime = secs/10;
    printf ("two TIMES rmatgen takes %lf sec\n", gentime);
    alpha = 1.0;
    beta = 1.0;
    m = M;
    n = N;
    k = K;
    for (i=0;i<PASS; i++){
        start_time();
        for (j=0;j<TIMES;j++) {
            rmatgen(mma, M, K, &norma);
            rmatgen(mmb, K, N, &norma);
            xgemm("n", "n", &m, &n, &k, &alpha, mma, &m, mmb, &k, &beta, mmc, &m);
        }
        end_time();
        printf ("%d. %d xGEMM's took %lf sec, %lf MFLOPS\n",
            i, TIMES, (secs-gentime), TIMES * (M*N*K*2.0/1000000)/(secs - gentime));
    }
    return 0;
}