Comparing bgq and ref implementation

bgq/Makefile.in

@@ -30,7 +30,7 @@ LDADD =
 COMPILE = ${CC} $(DEFS) ${INCLUDES} ${CPPFLAGS} ${CFLAGS}
 
 LIBRARIES = libbgq
-libbgq_TARGETS = bgq_gaugefield bgq_spinorfield bgq_comm bgq_utils bgq_field bgq_HoppingMatrix bgq_dispatch mypapi
+libbgq_TARGETS = bgq_qpx bgq_gaugefield bgq_spinorfield bgq_comm bgq_utils bgq_field bgq_HoppingMatrix bgq_dispatch mypapi
 libbgq_STARGETS = 
 
 libbgq_OBJECTS = $(addsuffix .o, ${libbgq_TARGETS})

bgq/bgq_HoppingMatrix.h

@@ -111,7 +111,7 @@ EXTERN_INLINE void bgq_HoppingMatrix_loadWeyllayout_raw(bgq_su3_spinor_params(*t
 	// Z- //////////////////////////////////////////////////////////////////////////
 	{
 		bgq_su3_weyl_decl(weyl_zdown);
-		bgq_su3_weyl_load_double(weyl_zdown, &spinorsite->d[ZDOWN]);
+		bgq_su3_weyl_load_double(weyl_zdown, &spinorsite->d[ZDOWN]);//TODO: Valgrind says invalid read
 		bgq_su3_accum_weyl_zdown(result, weyl_zdown);
 	}
 

bgq/bgq_field.c

@@ -983,6 +983,7 @@ void bgq_spinorfield_setup(bgq_weylfield_controlblock *field, bool isOdd, bool r
 			field->sec_recv[d] = (bgq_weyl_vec*) (weylbase + bgq_weyl_section_offset(bgq_direction2section(d,false)));
 		}
 		field->sec_surface = (bgq_weylsite*) (weylbase + bgq_weyl_section_offset(sec_surface));
+		field->sec_collapsed = field->sec_surface;
 		field->sec_body = (bgq_weylsite*) (weylbase + bgq_weyl_section_offset(sec_body));
 		field->sec_end = weylbase + bgq_weyl_section_offset(sec_end);
 

bgq/bgq_field.h

@@ -508,6 +508,7 @@ typedef struct {
 	uint8_t *sec_weyl; //TODO: can be made bgq_weyl_vec // corresponds to offset 0 for the following fields
 	bgq_weyl_vec *sec_send[PHYSICAL_LD];
 	bgq_weyl_vec *sec_recv[PHYSICAL_LD];
+	bgq_weylsite *sec_collapsed;
 	bgq_weylsite *sec_surface;
 	bgq_weylsite *sec_body;
 	uint8_t *sec_end;
@@ -536,6 +537,7 @@ EXTERN_FIELD size_t *g_bgq_index_body2halfvolume[PHYSICAL_LP];
 EXTERN_FIELD size_t *g_bgq_index_halfvolume2surface[PHYSICAL_LP]; // -1 if not surface
 EXTERN_FIELD size_t *g_bgq_index_halfvolume2body[PHYSICAL_LP]; // -1 if not body
 EXTERN_FIELD size_t *g_bgq_index_halfvolume2surfacebody[PHYSICAL_LP];
+//EXTERN_FIELD size_t *g_bgq_index_ordered2halfvolume[PHYSICAL_LP];
 
 // Mapping of dst weyls to memory offsets
 EXTERN_FIELD bgq_weyl_offsets_t *g_bgq_ih_dst2offset[PHYSICAL_LP];
@@ -582,6 +584,7 @@ EXTERN_INLINE size_t bgq_local2global_z(size_t z) {
 }
 
 
+
 EXTERN_INLINE bool bgq_local2isOdd(size_t t, size_t x, size_t y, size_t z) {
 	assert(0 <= t && t < LOCAL_LT);
 	assert(0 <= x && x < LOCAL_LX);
@@ -879,6 +882,72 @@ EXTERN_INLINE size_t bgq_local2volume(size_t t, size_t x, size_t y, size_t z) {
 	return bgq_halfvolume2volume(isOdd, ih, k);
 }
 
+EXTERN_INLINE size_t bgq_halfvolume2collapsed(bool isOdd, size_t ih) {
+	bool isSurface = bgq_halfvolume2isSurface(isOdd, ih);
+	size_t ic;
+	if (isSurface) {
+		size_t is = bgq_halfvolume2surface(isOdd, ih);
+		ic = is;
+	} else {
+		size_t ib = bgq_halfvolume2body(isOdd, ih);
+		ic = PHYSICAL_SURFACE + ib;
+	}
+	assert(0 <= ic && ic < PHYSICAL_VOLUME);
+	return ic;
+}
+
+EXTERN_INLINE bool bgq_collapsed2isSurface(size_t ic) {
+	assert(0 <= ic && ic < PHYSICAL_VOLUME);
+	return ic < PHYSICAL_SURFACE;
+}
+
+EXTERN_INLINE size_t bgq_collapsed2halfvolume(bool isOdd, size_t ic) {
+	assert(0 <= ic && ic < PHYSICAL_VOLUME);
+	size_t ih;
+	if (bgq_collapsed2isSurface(ic)) {
+		size_t is = ic;
+		ih = bgq_surface2halfvolume(isOdd, is);
+	} else {
+		size_t ib  = ic - PHYSICAL_SURFACE;
+		ih = bgq_body2halfvolume(isOdd, ib);
+	}
+	assert(0 <= ih && ih < PHYSICAL_VOLUME);
+	return ih;
+}
+
+EXTERN_INLINE size_t bgq_surface2collapsed(size_t is) {
+	assert(0 <= is && is < PHYSICAL_SURFACE);
+	return is;
+}
+EXTERN_INLINE size_t bgq_collapsed2surface(size_t ic) {
+	assert(0 <= ic && ic < PHYSICAL_VOLUME);
+	assert(bgq_collapsed2isSurface(ic));
+	return ic;
+}
+
+EXTERN_INLINE size_t bgq_body2collapsed(size_t ib) {
+	assert(0 <= ib && ib < PHYSICAL_BODY);
+	size_t ic = PHYSICAL_SURFACE + ib;
+	assert(0 <= ic && ic < PHYSICAL_VOLUME);
+	return ic;
+}
+EXTERN_INLINE size_t bgq_collapsed2body(size_t ic) {
+	assert(0 <= ic && ic < PHYSICAL_VOLUME);
+	assert(!bgq_collapsed2isSurface(ic));
+	return ic - PHYSICAL_SURFACE;
+}
+
+EXTERN_INLINE size_t bgq_local2collapsed(size_t t, size_t x, size_t y, size_t z) {
+	assert(0 <= t && t < LOCAL_LT);
+	assert(0 <= x && x < LOCAL_LX);
+	assert(0 <= y && y < LOCAL_LY);
+	assert(0 <= z && z < LOCAL_LZ);
+	bool isOdd = bgq_local2isOdd(t,x,y,z);
+	size_t ih = bgq_local2halfvolume(t,x,y,z);
+	return bgq_halfvolume2collapsed(isOdd, ih);
+}
+
+
 EXTERN_INLINE size_t bgq_local2halfvolume_neighbor(size_t t, size_t x, size_t y, size_t z, bgq_direction d) {
 	switch (d) {
 	case TUP:

bgq/bgq_qpx.c

@@ -0,0 +1,9 @@
+/*
+ * bgq_qpx.c
+ *
+ *  Created on: Oct 26, 2012
+ *      Author: meinersbur
+ */
+
+#define BGQ_QPX_C_
+#include "bgq_qpx.h"

bgq/bgq_qpx.h

@@ -12,6 +12,23 @@
 //#include <mpi.h>
 //#include "complex_c99.h"
 #include <string.h>
+#include <complex.h>
+
+#ifndef BGQ_QPX_C_
+#define EXTERN_INLINE EXTERN_INLINE_DECLARATION
+#define EXTERN_FIELD extern
+#define EXTERN_INIT(val)
+#else
+#define EXTERN_INLINE EXTERN_INLINE_DEFINITION
+#define EXTERN_FIELD
+#define EXTERN_INIT(val) = (val)
+#endif
+
+
+
+
+
+
 
 #ifndef BGQ_QPX
 #define BGQ_QPX 0
@@ -1252,6 +1269,75 @@ do {\
 
 
 
+EXTERN_INLINE complexdouble bgq_su3_spinor_getcomplex_double(bgq_su3_spinor_params(spinor), size_t v, size_t c, size_t k) {
+	assert(0 <= v && v < 4);
+	assert(0 <= c && c < 4);
+	assert(0 <= k && c < 2);
+
+	bgq_vector4double_decl(result);
+	switch (v) {
+	case 0:
+		switch (c) {
+		case 0:
+			bgq_mov(result,spinor_v0_c0);
+			break;
+		case 1:
+			bgq_mov(result,spinor_v0_c1);
+			break;
+		case 2:
+			bgq_mov(result,spinor_v0_c2);
+			break;
+		}
+		break;
+	case 1:
+		switch (c) {
+		case 0:
+			bgq_mov(result,spinor_v1_c0);
+			break;
+		case 1:
+			bgq_mov(result,spinor_v1_c1);
+			break;
+		case 2:
+			bgq_mov(result,spinor_v1_c2);
+			break;
+		}
+		break;
+	case 2:
+		switch (c) {
+		case 0:
+			bgq_mov(result,spinor_v2_c0);
+			break;
+		case 1:
+			bgq_mov(result,spinor_v2_c1);
+			break;
+		case 2:
+			bgq_mov(result,spinor_v2_c2);
+			break;
+		}
+		break;
+	case 3:
+		switch (c) {
+		case 0:
+			bgq_mov(result,spinor_v3_c0);
+			break;
+		case 1:
+			bgq_mov(result,spinor_v3_c1);
+			break;
+		case 2:
+			bgq_mov(result,spinor_v3_c2);
+			break;
+		}
+		break;
+	}
+
+	if (k==0) {
+		return bgq_elem0(result) + bgq_elem1(result) * _Complex_I;
+	} else {
+		return bgq_elem2(result) + bgq_elem3(result) * _Complex_I;
+	}
+}
+
+
 #if !BGQ_QPX
 
 // No semantic effects
@@ -1580,4 +1666,8 @@ static inline void mbar() {
 #endif
 
 
+#undef EXTERN_INLINE
+#undef EXTERN_FIELD
+#undef EXTERN_INIT
+
 #endif /* BGQ_QPX_H_ */

bgq/bgq_spinorfield.c

@@ -7,3 +7,89 @@
 
 #define BGQ_SPINORFIELD_C_
 #include "bgq_spinorfield.h"
+
+#include "bgq_HoppingMatrix.h"
+#include "bgq_qpx.h"
+
+#include "../geometry_eo.h"
+
+#include <mpi.h>
+
+
+typedef struct {
+	double _Complex c[3];
+} su3_vector_array64;
+
+typedef struct {
+	su3_vector_array64 v[4];
+} spinor_array64;
+
+double bgq_spinorfield_compare(bool isOdd, bgq_weylfield_controlblock *bgqfield, spinor *reffield, bool silent) {
+	assert(bgqfield);
+	assert(reffield);
+
+	bool readFulllayout = bgqfield->hasFullspinorData;
+	bgq_spinorfield_setup(bgqfield, isOdd, readFulllayout, false, !readFulllayout, false);
+
+	double diff_max = 0;
+	size_t count = 0;
+	for (size_t z = 0; z < LOCAL_LZ ; z += 1) {
+		for (size_t y = 0; y < LOCAL_LY ; y += 1) {
+			for (size_t x = 0; x < LOCAL_LX ; x += 1) {
+				for (size_t t = 0; t < LOCAL_LT ; t += 1) {
+					if (bgq_local2isOdd(t, x, y, z) != isOdd)
+						continue;
+
+					const int ix = g_ipt[t][x][y][z]; /* lexic coordinate */
+					assert(ix == Index(t,x,y,z));
+					int iy = g_lexic2eo[ix]; /* even/odd coordinate (even and odd sites in two different fields of size VOLUME/2, first even field followed by odd) */
+					assert(0 <= iy && iy < (VOLUME+RAND));
+					int icx = g_lexic2eosub[ix]; /*  even/odd coordinate relative to field base */
+					assert(0 <= icx && icx < VOLUME/2);
+					assert(icx == iy - (isOdd ? (VOLUME+RAND)/2 : 0));
+					spinor_array64 *sp = (spinor_array64*) &reffield[icx];
+
+					size_t ic = bgq_local2collapsed(t, x, y, z);
+					size_t k = bgq_local2k(t, x, y, z);
+
+					bgq_su3_spinor_decl(bgqspinor);
+					if (readFulllayout) {
+						bgq_HoppingMatrix_loadFulllayout_raw(bgq_su3_spinor_vars(&bgqspinor), &bgqfield->sec_fullspinor[ic]);
+					} else {
+						bgq_HoppingMatrix_loadWeyllayout_raw(bgq_su3_spinor_vars(&bgqspinor), &bgqfield->sec_collapsed[ic]);
+					}
+
+					for (size_t v = 0; v < 4; v += 1) {
+						for (size_t c = 0; c < 3; c += 1) {
+							complexdouble bgqvalue = bgq_su3_spinor_getcomplex_double(bgq_su3_spinor_vars(bgqspinor), v, c, k);
+							complexdouble refvalue = sp->v[v].c[c];
+
+							double diff = cabs(bgqvalue - refvalue);
+
+							if (diff > 0.01) {
+								if (!silent)
+									master_print("Coordinate (%zu,%zu,%zu,%zu)(%zu,%zu): ref=(%f + %fi) != bgb=(%f + %fi) off by %f\n", t, x, y, z, v, c, creal(refvalue), cimag(refvalue), creal(bgqvalue), cimag(bgqvalue), diff);
+								count += 1;
+							}
+							if (diff > diff_max) {
+								diff_max = diff;
+							}
+						}
+					}
+				}
+			}
+		}
+	}
+
+	double global_diff_max;
+	MPI_Allreduce(&diff_max, &global_diff_max, 2, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+	return diff_max;
+
+	if (count > 0) {
+		if (!silent)
+			master_print("%zu sites of %d wrong\n", count, VOLUME/2);
+	}
+
+	return global_diff_max;
+}
+

bgq/bgq_spinorfield.h

@@ -8,8 +8,11 @@
 #ifndef BGQ_SPINORFIELD_H_
 #define BGQ_SPINORFIELD_H_
 
+#include "bgq_field.h"
 #include "bgq_utils.h"
 
+#include <stdbool.h>
+
 #ifndef BGQ_SPINORFIELD_C_
 #define EXTERN_INLINE EXTERN_INLINE_DECLARATION
 #define EXTERN_FIELD extern
@@ -20,6 +23,7 @@
 #define EXTERN_INIT(val) = (val)
 #endif
 
+double bgq_spinorfield_compare(bool isOdd, bgq_weylfield_controlblock *bgqfield, spinor *reffield, bool silent);
 
 #undef EXTERN_INLINE
 #undef EXTERN_FIELD