Squashed 'third_party/blasfeo/' content from commit 2a828ca
Change-Id: If1c3caa4799b2d4eb287ef83fa17043587ef07a3
git-subtree-dir: third_party/blasfeo
git-subtree-split: 2a828ca5442108c4c58e4b42b061a0469043f6ea
diff --git a/test_problems/test_blas_d.c b/test_problems/test_blas_d.c
new file mode 100644
index 0000000..1e71494
--- /dev/null
+++ b/test_problems/test_blas_d.c
@@ -0,0 +1,480 @@
+/**************************************************************************************************
+* *
+* This file is part of BLASFEO. *
+* *
+* BLASFEO -- BLAS For Embedded Optimization. *
+* Copyright (C) 2016-2017 by Gianluca Frison. *
+* Developed at IMTEK (University of Freiburg) under the supervision of Moritz Diehl. *
+* All rights reserved. *
+* *
+* HPMPC is free software; you can redistribute it and/or *
+* modify it under the terms of the GNU Lesser General Public *
+* License as published by the Free Software Foundation; either *
+* version 2.1 of the License, or (at your option) any later version. *
+* *
+* HPMPC is distributed in the hope that it will be useful, *
+* but WITHOUT ANY WARRANTY; without even the implied warranty of *
+* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
+* See the GNU Lesser General Public License for more details. *
+* *
+* You should have received a copy of the GNU Lesser General Public *
+* License along with HPMPC; if not, write to the Free Software *
+* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA *
+* *
+* Author: Gianluca Frison, giaf (at) dtu.dk *
+* gianluca.frison (at) imtek.uni-freiburg.de *
+* *
+**************************************************************************************************/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <sys/time.h>
+
+//#if defined(TARGET_X64_INTEL_SANDY_BRIDGE)
+//#include <xmmintrin.h> // needed to flush to zero sub-normals with _MM_SET_FLUSH_ZERO_MODE (_MM_FLUSH_ZERO_ON); in the main()
+//#endif
+
+#include "../include/blasfeo_common.h"
+#include "../include/blasfeo_d_aux_ext_dep.h"
+#include "../include/blasfeo_d_aux.h"
+#include "../include/blasfeo_i_aux_ext_dep.h"
+#include "../include/blasfeo_v_aux_ext_dep.h"
+#include "../include/blasfeo_d_kernel.h"
+#include "../include/blasfeo_d_blas.h"
+
+#ifndef D_PS
+#define D_PS 1
+#endif
+#ifndef D_NC
+#define D_NC 1
+#endif
+
+
+
+#if defined(REF_BLAS_OPENBLAS)
+void openblas_set_num_threads(int num_threads);
+#endif
+#if defined(REF_BLAS_BLIS)
+void omp_set_num_threads(int num_threads);
+#endif
+#if defined(REF_BLAS_MKL)
+#include "mkl.h"
+#endif
+
+
+
+#include "cpu_freq.h"
+
+
+
+int main()
+ {
+
+#if defined(REF_BLAS_OPENBLAS)
+ openblas_set_num_threads(1);
+#endif
+#if defined(REF_BLAS_BLIS)
+ omp_set_num_threads(1);
+#endif
+#if defined(REF_BLAS_MKL)
+ mkl_set_num_threads(1);
+#endif
+
+//#if defined(TARGET_X64_INTEL_SANDY_BRIDGE)
+// _MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON); // flush to zero subnormals !!! works only with one thread !!!
+//#endif
+
+ printf("\n");
+ printf("\n");
+ printf("\n");
+
+ printf("BLAS performance test - double precision\n");
+ printf("\n");
+
+ // maximum frequency of the processor
+ const float GHz_max = GHZ_MAX;
+ printf("Frequency used to compute theoretical peak: %5.1f GHz (edit test_param.h to modify this value).\n", GHz_max);
+ printf("\n");
+
+ // maximum flops per cycle, double precision
+#if defined(TARGET_X64_INTEL_HASWELL)
+ const float flops_max = 16;
+ printf("Testing BLAS version for AVX2 and FMA instruction sets, 64 bit (optimized for Intel Haswell): theoretical peak %5.1f Gflops\n", flops_max*GHz_max);
+#elif defined(TARGET_X64_INTEL_SANDY_BRIDGE)
+ const float flops_max = 8;
+ printf("Testing BLAS version for AVX instruction set, 64 bit (optimized for Intel Sandy Bridge): theoretical peak %5.1f Gflops\n", flops_max*GHz_max);
+#elif defined(TARGET_X64_INTEL_CORE)
+ const float flops_max = 4;
+ printf("Testing BLAS version for SSE3 instruction set, 64 bit (optimized for Intel Core): theoretical peak %5.1f Gflops\n", flops_max*GHz_max);
+#elif defined(TARGET_X64_AMD_BULLDOZER)
+ const float flops_max = 8;
+ printf("Testing BLAS version for SSE3 and FMA instruction set, 64 bit (optimized for AMD Bulldozer): theoretical peak %5.1f Gflops\n", flops_max*GHz_max);
+#elif defined(TARGET_ARMV8A_ARM_CORTEX_A57)
+ const float flops_max = 4;
+ printf("Testing BLAS version for NEONv2 instruction set, 64 bit (optimized for ARM Cortex A57): theoretical peak %5.1f Gflops\n", flops_max*GHz_max);
+#elif defined(TARGET_ARMV7A_ARM_CORTEX_A15)
+ const float flops_max = 2;
+ printf("Testing BLAS version for VFPv4 instruction set, 32 bit (optimized for ARM Cortex A15): theoretical peak %5.1f Gflops\n", flops_max*GHz_max);
+#elif defined(TARGET_GENERIC)
+ const float flops_max = 2;
+ printf("Testing BLAS version for generic scalar instruction set: theoretical peak %5.1f Gflops ???\n", flops_max*GHz_max);
+#endif
+
+// FILE *f;
+// f = fopen("./test_problems/results/test_blas.m", "w"); // a
+
+#if defined(TARGET_X64_INTEL_HASWELL)
+// fprintf(f, "C = 'd_x64_intel_haswell';\n");
+// fprintf(f, "\n");
+#elif defined(TARGET_X64_INTEL_SANDY_BRIDGE)
+// fprintf(f, "C = 'd_x64_intel_sandybridge';\n");
+// fprintf(f, "\n");
+#elif defined(TARGET_X64_INTEL_CORE)
+// fprintf(f, "C = 'd_x64_intel_core';\n");
+// fprintf(f, "\n");
+#elif defined(TARGET_X64_AMD_BULLDOZER)
+// fprintf(f, "C = 'd_x64_amd_bulldozer';\n");
+// fprintf(f, "\n");
+#elif defined(TARGET_ARMV8A_ARM_CORTEX_A57)
+// fprintf(f, "C = 'd_armv8a_arm_cortex_a57';\n");
+// fprintf(f, "\n");
+#elif defined(TARGET_ARMV7A_ARM_CORTEX_A15)
+// fprintf(f, "C = 'd_armv7a_arm_cortex_a15';\n");
+// fprintf(f, "\n");
+#elif defined(TARGET_GENERIC)
+// fprintf(f, "C = 'd_generic';\n");
+// fprintf(f, "\n");
+#endif
+
+// fprintf(f, "A = [%f %f];\n", GHz_max, flops_max);
+// fprintf(f, "\n");
+
+// fprintf(f, "B = [\n");
+
+
+
+ int i, j, rep, ll;
+
+ const int bsd = D_PS;
+ const int ncd = D_NC;
+
+/* int info = 0;*/
+
+ printf("\nn\t dgemm_blasfeo\t dgemm_blas\n");
+ printf("\nn\t Gflops\t %%\t Gflops\n\n");
+
+#if 1
+ int nn[] = {4, 8, 12, 16, 20, 24, 28, 32, 36, 40, 44, 48, 52, 56, 60, 64, 68, 72, 76, 80, 84, 88, 92, 96, 100, 104, 108, 112, 116, 120, 124, 128, 132, 136, 140, 144, 148, 152, 156, 160, 164, 168, 172, 176, 180, 184, 188, 192, 196, 200, 204, 208, 212, 216, 220, 224, 228, 232, 236, 240, 244, 248, 252, 256, 260, 264, 268, 272, 276, 280, 284, 288, 292, 296, 300, 304, 308, 312, 316, 320, 324, 328, 332, 336, 340, 344, 348, 352, 356, 360, 364, 368, 372, 376, 380, 384, 388, 392, 396, 400, 404, 408, 412, 416, 420, 424, 428, 432, 436, 440, 444, 448, 452, 456, 460, 500, 550, 600, 650, 700};
+ int nnrep[] = {10000, 10000, 10000, 10000, 10000, 10000, 10000, 10000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 400, 400, 400, 400, 400, 200, 200, 200, 200, 200, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 20, 20, 20, 20, 20, 20, 20, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 4, 4, 4, 4, 4};
+
+// for(ll=0; ll<24; ll++)
+ for(ll=0; ll<75; ll++)
+// for(ll=0; ll<115; ll++)
+// for(ll=0; ll<120; ll++)
+
+ {
+
+ int n = nn[ll];
+ int nrep = nnrep[ll];
+// int n = ll+1;
+// int nrep = nnrep[0];
+// n = n<12 ? 12 : n;
+// n = n<8 ? 8 : n;
+
+#else
+ int nn[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24};
+
+ for(ll=0; ll<24; ll++)
+
+ {
+
+ int n = nn[ll];
+ int nrep = 40000; //nnrep[ll];
+#endif
+
+ double *A; d_zeros(&A, n, n);
+ double *B; d_zeros(&B, n, n);
+ double *C; d_zeros(&C, n, n);
+ double *M; d_zeros(&M, n, n);
+
+ char c_n = 'n';
+ char c_l = 'l';
+ char c_r = 'r';
+ char c_t = 't';
+ char c_u = 'u';
+ int i_1 = 1;
+ int i_t;
+ double d_1 = 1;
+ double d_0 = 0;
+
+ for(i=0; i<n*n; i++)
+ A[i] = i;
+
+ for(i=0; i<n; i++)
+ B[i*(n+1)] = 1;
+
+ for(i=0; i<n*n; i++)
+ M[i] = 1;
+
+ int n2 = n*n;
+ double *B2; d_zeros(&B2, n, n);
+ for(i=0; i<n*n; i++)
+ B2[i] = 1e-15;
+ for(i=0; i<n; i++)
+ B2[i*(n+1)] = 1;
+
+ int pnd = ((n+bsd-1)/bsd)*bsd;
+ int cnd = ((n+ncd-1)/ncd)*ncd;
+ int cnd2 = 2*((n+ncd-1)/ncd)*ncd;
+
+ double *x; d_zeros_align(&x, pnd, 1);
+ double *y; d_zeros_align(&y, pnd, 1);
+ double *x2; d_zeros_align(&x2, pnd, 1);
+ double *y2; d_zeros_align(&y2, pnd, 1);
+ double *diag; d_zeros_align(&diag, pnd, 1);
+ int *ipiv; int_zeros(&ipiv, n, 1);
+
+ for(i=0; i<pnd; i++) x[i] = 1;
+ for(i=0; i<pnd; i++) x2[i] = 1;
+
+ // matrix struct
+#if 0
+ struct d_strmat sA; d_allocate_strmat(n+4, n+4, &sA);
+ struct d_strmat sB; d_allocate_strmat(n+4, n+4, &sB);
+ struct d_strmat sC; d_allocate_strmat(n+4, n+4, &sC);
+ struct d_strmat sD; d_allocate_strmat(n+4, n+4, &sD);
+ struct d_strmat sE; d_allocate_strmat(n+4, n+4, &sE);
+#else
+ struct d_strmat sA; d_allocate_strmat(n, n, &sA);
+ struct d_strmat sB; d_allocate_strmat(n, n, &sB);
+ struct d_strmat sB2; d_allocate_strmat(n, n, &sB2);
+ struct d_strmat sB3; d_allocate_strmat(n, n, &sB3);
+ struct d_strmat sC; d_allocate_strmat(n, n, &sC);
+ struct d_strmat sD; d_allocate_strmat(n, n, &sD);
+ struct d_strmat sE; d_allocate_strmat(n, n, &sE);
+#endif
+ struct d_strvec sx; d_allocate_strvec(n, &sx);
+ struct d_strvec sy; d_allocate_strvec(n, &sy);
+ struct d_strvec sz; d_allocate_strvec(n, &sz);
+
+ d_cvt_mat2strmat(n, n, A, n, &sA, 0, 0);
+ d_cvt_mat2strmat(n, n, B, n, &sB, 0, 0);
+ d_cvt_mat2strmat(n, n, B2, n, &sB2, 0, 0);
+ d_cvt_vec2strvec(n, x, &sx, 0);
+ int ii;
+ for(ii=0; ii<n; ii++)
+ {
+ DMATEL_LIBSTR(&sB3, ii, ii) = 1.0;
+// DMATEL_LIBSTR(&sB3, n-1, ii) = 1.0;
+ DMATEL_LIBSTR(&sB3, ii, n-1) = 1.0;
+ DVECEL_LIBSTR(&sx, ii) = 1.0;
+ }
+// d_print_strmat(n, n, &sB3, 0, 0);
+// if(n==20) return;
+
+ int qr_work_size = 0;//dgeqrf_work_size_libstr(n, n);
+ void *qr_work;
+ v_zeros_align(&qr_work, qr_work_size);
+
+ int lq_work_size = 0;//dgelqf_work_size_libstr(n, n);
+ void *lq_work;
+ v_zeros_align(&lq_work, lq_work_size);
+
+ // create matrix to pivot all the time
+// dgemm_nt_libstr(n, n, n, 1.0, &sA, 0, 0, &sA, 0, 0, 1.0, &sB, 0, 0, &sD, 0, 0);
+
+ double *dummy;
+
+ int info;
+
+ /* timing */
+ struct timeval tvm1, tv0, tv1, tv2, tv3, tv4, tv5, tv6, tv7, tv8, tv9, tv10, tv11, tv12, tv13, tv14, tv15, tv16;
+
+ /* warm up */
+ for(rep=0; rep<nrep; rep++)
+ {
+ dgemm_nt_libstr(n, n, n, 1.0, &sA, 0, 0, &sA, 0, 0, 1.0, &sB, 0, 0, &sC, 0, 0);
+ }
+
+ double alpha = 1.0;
+ double beta = 0.0;
+
+ gettimeofday(&tv0, NULL); // stop
+
+ for(rep=0; rep<nrep; rep++)
+ {
+
+// dgemm_nt_lib(n, n, n, 1.0, pA, cnd, pB, cnd, 0.0, pC, cnd, pC, cnd);
+// dgemm_nn_lib(n, n, n, 1.0, pA, cnd, pB, cnd, 0.0, pC, cnd, pC, cnd);
+// dsyrk_nt_l_lib(n, n, n, 1.0, pA, cnd, pB, cnd, 1.0, pC, cnd, pD, cnd);
+// dtrmm_nt_ru_lib(n, n, pA, cnd, pB, cnd, 0, pC, cnd, pD, cnd);
+// dpotrf_nt_l_lib(n, n, pB, cnd, pD, cnd, diag);
+// dsyrk_dpotrf_nt_l_lib(n, n, n, pA, cnd, pA, cnd, 1, pB, cnd, pD, cnd, diag);
+// dsyrk_nt_l_lib(n, n, n, pA, cnd, pA, cnd, 1, pB, cnd, pD, cnd);
+// dpotrf_nt_l_lib(n, n, pD, cnd, pD, cnd, diag);
+// dgetrf_nn_nopivot_lib(n, n, pB, cnd, pB, cnd, diag);
+// dgetrf_nn_lib(n, n, pB, cnd, pB, cnd, diag, ipiv);
+// dtrsm_nn_ll_one_lib(n, n, pD, cnd, pB, cnd, pB, cnd);
+// dtrsm_nn_lu_inv_lib(n, n, pD, cnd, diag, pB, cnd, pB, cnd);
+ }
+
+ gettimeofday(&tv1, NULL); // stop
+
+ for(rep=0; rep<nrep; rep++)
+ {
+// kernel_dgemm_nt_12x4_lib4(n, &alpha, sA.pA, sA.cn, sB.pA, &beta, sD.pA, sD.cn, sD.pA, sD.cn);
+// kernel_dgemm_nt_8x8_lib4(n, &alpha, sA.pA, sA.cn, sB.pA, sB.cn, &beta, sD.pA, sD.cn, sD.pA, sD.cn);
+// kernel_dsyrk_nt_l_8x8_lib4(n, &alpha, sA.pA, sA.cn, sB.pA, sB.cn, &beta, sD.pA, sD.cn, sD.pA, sD.cn);
+// kernel_dgemm_nt_8x4_lib4(n, &alpha, sA.pA, sA.cn, sB.pA, &beta, sD.pA, sD.cn, sD.pA, sD.cn);
+// kernel_dgemm_nt_4x8_lib4(n, &alpha, sA.pA, sB.pA, sB.cn, &beta, sD.pA, sD.pA);
+// kernel_dgemm_nt_4x4_lib4(n, &alpha, sA.pA, sB.pA, &beta, sD.pA, sD.pA);
+// kernel_dger4_12_sub_lib4(n, sA.pA, sA.cn, sB.pA, sD.pA, sD.cn);
+// kernel_dger4_sub_12r_lib4(n, sA.pA, sA.cn, sB.pA, sD.pA, sD.cn);
+// kernel_dger4_sub_8r_lib4(n, sA.pA, sA.cn, sB.pA, sD.pA, sD.cn);
+// kernel_dger12_add_4r_lib4(n, sA.pA, sB.pA, sB.cn, sD.pA);
+// kernel_dger8_add_4r_lib4(n, sA.pA, sB.pA, sB.cn, sD.pA);
+// kernel_dger4_sub_4r_lib4(n, sA.pA, sB.pA, sD.pA);
+// kernel_dger2_sub_4r_lib4(n, sA.pA, sB.pA, sD.pA);
+// kernel_dger4_sub_8c_lib4(n, sA.pA, sA.cn, sB.pA, sD.pA, sD.cn);
+// kernel_dger4_sub_4c_lib4(n, sA.pA, sA.cn, sB.pA, sD.pA, sD.cn);
+// kernel_dgemm_nn_4x12_lib4(n, &alpha, sA.pA, 0, sB.pA, sB.cn, &beta, sD.pA, sD.pA);
+// kernel_dgemm_nn_4x8_lib4(n, &alpha, sA.pA, 0, sB.pA, sB.cn, &beta, sD.pA, sD.pA);
+// kernel_dgemm_nn_4x4_lib4(n, &alpha, sA.pA, 0, sB.pA, sB.cn, &beta, sD.pA, sD.pA);
+
+ dgemm_nt_libstr(n, n, n, 1.0, &sA, 0, 0, &sB, 0, 0, 0.0, &sC, 0, 0, &sD, 0, 0);
+// dgemm_nn_libstr(n, n, n, 1.0, &sA, 0, 0, &sB, 0, 0, 0.0, &sC, 0, 0, &sD, 0, 0);
+// dsyrk_ln_libstr(n, n, 1.0, &sA, 0, 0, &sA, 0, 0, 0.0, &sC, 0, 0, &sD, 0, 0);
+// dsyrk_ln_mn_libstr(n, n, n, 1.0, &sA, 0, 0, &sA, 0, 0, 0.0, &sC, 0, 0, &sD, 0, 0);
+// dpotrf_l_mn_libstr(n, n, &sB, 0, 0, &sB, 0, 0);
+// dpotrf_l_libstr(n, &sB, 0, 0, &sB, 0, 0);
+// dgetrf_nopivot_libstr(n, n, &sB, 0, 0, &sB, 0, 0);
+// dgetrf_libstr(n, n, &sB, 0, 0, &sB, 0, 0, ipiv);
+// dgeqrf_libstr(n, n, &sC, 0, 0, &sD, 0, 0, qr_work);
+// dcolin_libstr(n, &sx, 0, &sB3, 0, n-1);
+// dgelqf_libstr(n, n, &sB3, 0, 0, &sB3, 0, 0, lq_work);
+// dtrmm_rlnn_libstr(n, n, 1.0, &sA, 0, 0, &sD, 0, 0, &sD, 0, 0); //
+// dtrmm_rutn_libstr(n, n, 1.0, &sA, 0, 0, &sB, 0, 0, &sD, 0, 0);
+// dtrsm_llnu_libstr(n, n, 1.0, &sD, 0, 0, &sB, 0, 0, &sB, 0, 0);
+// dtrsm_lunn_libstr(n, n, 1.0, &sD, 0, 0, &sB, 0, 0, &sB, 0, 0);
+// dtrsm_rltn_libstr(n, n, 1.0, &sB2, 0, 0, &sD, 0, 0, &sD, 0, 0); //
+// dtrsm_rltu_libstr(n, n, 1.0, &sD, 0, 0, &sB, 0, 0, &sB, 0, 0);
+// dtrsm_rutn_libstr(n, n, 1.0, &sD, 0, 0, &sB, 0, 0, &sB, 0, 0);
+// dgemv_n_libstr(n, n, 1.0, &sA, 0, 0, &sx, 0, 0.0, &sy, 0, &sz, 0);
+// dgemv_t_libstr(n, n, 1.0, &sA, 0, 0, &sx, 0, 0.0, &sy, 0, &sz, 0);
+// dsymv_l_libstr(n, n, 1.0, &sA, 0, 0, &sx, 0, 0.0, &sy, 0, &sz, 0);
+// dgemv_nt_libstr(n, n, 1.0, 1.0, &sA, 0, 0, &sx, 0, &sx, 0, 0.0, 0.0, &sy, 0, &sy, 0, &sz, 0, &sz, 0);
+ }
+
+// d_print_strmat(n, n, &sD, 0, 0);
+
+ gettimeofday(&tv2, NULL); // stop
+
+ for(rep=0; rep<nrep; rep++)
+ {
+#if defined(REF_BLAS_OPENBLAS) || defined(REF_BLAS_NETLIB) || defined(REF_BLAS_MKL)
+// dgemm_(&c_n, &c_t, &n, &n, &n, &d_1, A, &n, M, &n, &d_0, C, &n);
+// dpotrf_(&c_l, &n, B2, &n, &info);
+// dgemm_(&c_n, &c_n, &n, &n, &n, &d_1, A, &n, M, &n, &d_0, C, &n);
+// dsyrk_(&c_l, &c_n, &n, &n, &d_1, A, &n, &d_0, C, &n);
+// dtrmm_(&c_r, &c_u, &c_t, &c_n, &n, &n, &d_1, A, &n, C, &n);
+// dgetrf_(&n, &n, B2, &n, ipiv, &info);
+// dtrsm_(&c_l, &c_l, &c_n, &c_u, &n, &n, &d_1, B2, &n, B, &n);
+// dtrsm_(&c_l, &c_u, &c_n, &c_n, &n, &n, &d_1, B2, &n, B, &n);
+// dtrtri_(&c_l, &c_n, &n, B2, &n, &info);
+// dlauum_(&c_l, &n, B, &n, &info);
+// dgemv_(&c_n, &n, &n, &d_1, A, &n, x, &i_1, &d_0, y, &i_1);
+// dgemv_(&c_t, &n, &n, &d_1, A, &n, x2, &i_1, &d_0, y2, &i_1);
+// dtrmv_(&c_l, &c_n, &c_n, &n, B, &n, x, &i_1);
+// dtrsv_(&c_l, &c_n, &c_n, &n, B, &n, x, &i_1);
+// dsymv_(&c_l, &n, &d_1, A, &n, x, &i_1, &d_0, y, &i_1);
+
+// for(i=0; i<n; i++)
+// {
+// i_t = n-i;
+// dcopy_(&i_t, &B[i*(n+1)], &i_1, &C[i*(n+1)], &i_1);
+// }
+// dsyrk_(&c_l, &c_n, &n, &n, &d_1, A, &n, &d_1, C, &n);
+// dpotrf_(&c_l, &n, C, &n, &info);
+
+#endif
+
+#if defined(REF_BLAS_BLIS)
+// dgemm_(&c_n, &c_t, &n77, &n77, &n77, &d_1, A, &n77, B, &n77, &d_0, C, &n77);
+// dgemm_(&c_n, &c_n, &n77, &n77, &n77, &d_1, A, &n77, B, &n77, &d_0, C, &n77);
+// dsyrk_(&c_l, &c_n, &n77, &n77, &d_1, A, &n77, &d_0, C, &n77);
+// dtrmm_(&c_r, &c_u, &c_t, &c_n, &n77, &n77, &d_1, A, &n77, C, &n77);
+// dpotrf_(&c_l, &n77, B, &n77, &info);
+// dtrtri_(&c_l, &c_n, &n77, B, &n77, &info);
+// dlauum_(&c_l, &n77, B, &n77, &info);
+#endif
+ }
+
+ gettimeofday(&tv3, NULL); // stop
+
+ float Gflops_max = flops_max * GHz_max;
+
+// float flop_operation = 4*16.0*2*n; // kernel 12x4
+// float flop_operation = 3*16.0*2*n; // kernel 12x4
+// float flop_operation = 2*16.0*2*n; // kernel 8x4
+// float flop_operation = 1*16.0*2*n; // kernel 4x4
+// float flop_operation = 0.5*16.0*2*n; // kernel 2x4
+
+ float flop_operation = 2.0*n*n*n; // dgemm
+// float flop_operation = 1.0*n*n*n; // dsyrk dtrmm dtrsm
+// float flop_operation = 1.0/3.0*n*n*n; // dpotrf dtrtri
+// float flop_operation = 2.0/3.0*n*n*n; // dgetrf
+// float flop_operation = 4.0/3.0*n*n*n; // dgeqrf
+// float flop_operation = 2.0*n*n; // dgemv dsymv
+// float flop_operation = 1.0*n*n; // dtrmv dtrsv
+// float flop_operation = 4.0*n*n; // dgemv_nt
+
+// float flop_operation = 4.0/3.0*n*n*n; // dsyrk+dpotrf
+
+ float time_hpmpc = (float) (tv1.tv_sec-tv0.tv_sec)/(nrep+0.0)+(tv1.tv_usec-tv0.tv_usec)/(nrep*1e6);
+ float time_blasfeo = (float) (tv2.tv_sec-tv1.tv_sec)/(nrep+0.0)+(tv2.tv_usec-tv1.tv_usec)/(nrep*1e6);
+ float time_blas = (float) (tv3.tv_sec-tv2.tv_sec)/(nrep+0.0)+(tv3.tv_usec-tv2.tv_usec)/(nrep*1e6);
+
+ float Gflops_hpmpc = 1e-9*flop_operation/time_hpmpc;
+ float Gflops_blasfeo = 1e-9*flop_operation/time_blasfeo;
+ float Gflops_blas = 1e-9*flop_operation/time_blas;
+
+
+// printf("%d\t%7.2f\t%7.2f\t%7.2f\t%7.2f\t%7.2f\t%7.2f\n", n, Gflops_hpmpc, 100.0*Gflops_hpmpc/Gflops_max, Gflops_blasfeo, 100.0*Gflops_blasfeo/Gflops_max, Gflops_blas, 100.0*Gflops_blas/Gflops_max);
+// fprintf(f, "%d\t%7.2f\t%7.2f\t%7.2f\t%7.2f\t%7.2f\t%7.2f\n", n, Gflops_hpmpc, 100.0*Gflops_hpmpc/Gflops_max, Gflops_blasfeo, 100.0*Gflops_blasfeo/Gflops_max, Gflops_blas, 100.0*Gflops_blas/Gflops_max);
+ printf("%d\t%7.2f\t%7.2f\t%7.2f\t%7.2f\n", n, Gflops_blasfeo, 100.0*Gflops_blasfeo/Gflops_max, Gflops_blas, 100.0*Gflops_blas/Gflops_max);
+// fprintf(f, "%d\t%7.2f\t%7.2f\t%7.2f\t%7.2f\n", n, Gflops_blasfeo, 100.0*Gflops_blasfeo/Gflops_max, Gflops_blas, 100.0*Gflops_blas/Gflops_max);
+
+
+ d_free(A);
+ d_free(B);
+ d_free(B2);
+ d_free(M);
+ d_free_align(x);
+ d_free_align(y);
+ d_free_align(x2);
+ d_free_align(y2);
+ int_free(ipiv);
+ free(qr_work);
+ free(lq_work);
+
+ d_free_strmat(&sA);
+ d_free_strmat(&sB);
+ d_free_strmat(&sB2);
+ d_free_strmat(&sB3);
+ d_free_strmat(&sC);
+ d_free_strmat(&sD);
+ d_free_strmat(&sE);
+ d_free_strvec(&sx);
+ d_free_strvec(&sy);
+ d_free_strvec(&sz);
+
+ }
+
+ printf("\n");
+
+// fprintf(f, "];\n");
+// fclose(f);
+
+ return 0;
+
+ }