1 #include "cado.h"
2 #include "bblas.hpp"
3 #include <cinttypes>
4 #include "test_bblas_level4.hpp"
5 #include "time_bblas_common.hpp"
6
7 #ifdef HAVE_M4RI
8 /* m4ri types are defined in m4ri/misc.h
9 * rci_t: rows and column indices
10 * wi_t: "indices for array of words that make up a row"
11 * word: typical data structure to represent packed bits.
12 *
13 * rci_t and wi_t are currently int, and will probably remain
14 * undistinguishable anyway.
15 *
16 * word is currently uint64_t, and we are really assuming that this stays
17 * so
18 */
19 static_assert(std::is_same<word,uint64_t>::value, "m4ri word MUST be uint64_t");
20
my_mzd_randomize(mzd_t * A,gmp_randstate_t rstate)21 void my_mzd_randomize(mzd_t * A, gmp_randstate_t rstate)/*{{{*/
22 {
23 for(rci_t i = 0 ; i < A->nrows ; i++) {
24 uint64_t * ptr = (uint64_t *) A->rows[i];
25 memfill_random(ptr, A->width * sizeof(uint64_t), rstate);
26 }
27 }/*}}}*/
28
mzd_set_mem(mzd_t * M,const uint64_t * s,unsigned int n)29 void test_bblas_base::mzd_set_mem(mzd_t * M, const uint64_t * s, unsigned int n)/*{{{*/
30 {
31 ASSERT_ALWAYS(M->width == 1);
32 ASSERT_ALWAYS(M->nrows == (rci_t) n);
33 for(rci_t i = 0 ; i < M->nrows ; i++) {
34 uint64_t * ptr = (uint64_t *) M->rows[i];
35 /* back in the days (before version 20111203 at least), m4ri had
36 * limbs in the wrong order */
37 ptr[0] = s[i];
38 }
39 }/*}}}*/
40
mzd_set_memT(mzd_t * M,const uint64_t * s,unsigned int n)41 void test_bblas_base::mzd_set_memT(mzd_t * M, const uint64_t * s, unsigned int n)/*{{{*/
42 {
43 mzd_set_mem(M, s, n);
44 mzd_transpose(M, M);
45 }/*}}}*/
46
mzd_check_mem(mzd_t * M,uint64_t * s,unsigned int n)47 void test_bblas_base::mzd_check_mem(mzd_t * M, uint64_t * s, unsigned int n)/*{{{*/
48 {
49 ASSERT_ALWAYS(M->width == 1);
50 ASSERT_ALWAYS(M->nrows >= (rci_t) n);
51 for(unsigned int i = 0 ; i < n ; i++) {
52 uint64_t * ptr = (uint64_t *) M->rows[i];
53 /* back in the days (before version 20111203 at least), m4ri had
54 * limbs in the wrong order */
55 if (ptr[0] != s[i]) {
56 fprintf(stderr, "Rows %d differ: %016" PRIx64 " != %016" PRIx64 "\n",
57 (int) i, ptr[0], s[i]);
58 abort();
59 }
60 }
61 }/*}}}*/
62
63 #if 0
64 void mzd_check_memT(mzd_t * M, uint64_t * s, unsigned int n)/*{{{*/
65 {
66 mzd_t * tmp = mzd_transpose(NULL, M);
67 mzd_check_mem(M, s, n);
68 mzd_free(tmp);
69 }/*}}}*/
70 #endif
71
72 #endif /* HAVE_M4RI */
73
74
75
76
77
78 /*************/
79
80
81 #ifdef HAVE_M4RI
mzd_mypluq(mzd_t * LU,mzd_t * M,mzp_t * P,mzp_t * Q,int c)82 static inline void mzd_mypluq(mzd_t * LU, mzd_t * M, mzp_t * P, mzp_t * Q, int c)
83 {
84 mzd_copy(LU, M);
85 mzd_pluq(LU,P,Q,c);
86 }
mzd_myechelonize_m4ri(mzd_t * E,mzd_t * M,int full,int k)87 static inline void mzd_myechelonize_m4ri(mzd_t * E, mzd_t * M, int full, int k)
88 {
89 mzd_copy(E, M);
90 mzd_echelonize_m4ri(E,full,k);
91 }
mzd_myechelonize_pluq(mzd_t * E,mzd_t * M,int full)92 static inline void mzd_myechelonize_pluq(mzd_t * E, mzd_t * M, int full)
93 {
94 mzd_copy(E, M);
95 mzd_echelonize_pluq(E,full);
96 }
97 #endif
98
99
100 #ifdef HAVE_M4RI
m4ri_plu_tests(int n)101 void test_bblas_level4::m4ri_plu_tests(int n __attribute__((unused)))
102 {
103 mzd_t * M;
104 mzd_t * LU;
105 mzp_t * P, *Q;
106 M = mzd_init(n, n);
107 #if 0
108 mzd_set_mem(M, m, n);
109 uint64_t * m = new uint64_t[n*n/64];
110 memfill_random(m, (64) * sizeof(uint64_t), rstate);
111 delete[] m;
112 #else
113 my_mzd_randomize(M, rstate);
114 #endif
115 LU = mzd_init(n, n);
116 P = mzp_init(n);
117 Q = mzp_init(n);
118 TIME1N(2, mzd_mypluq, LU, M, P, Q, 0);
119 TIME1N(2, mzd_myechelonize_m4ri, LU, M, 0, 0);
120 TIME1N(2, mzd_myechelonize_pluq, LU, M, 0);
121 mzd_free(M);
122 mzd_free(LU);
123 mzp_free(P);
124 mzp_free(Q);
125 }
126 #endif
127
128
129
130