GrB_Index* pos = new GrB_Index[2];

pos[0] = 0;

pos[1] = size;

taco_tensor_t ind;

ind.indices = new uint8_t**[1];

ind.indices[0] = new uint8_t*[2];

ind.indices[0][0] = (uint8_t*)pos;

ind.indices[0][1] = (uint8_t*)indices;

return ind;

taco_tensor_t vec;

vec.dimensions = new int32_t[1];

vec.dimensions[0] = N;

vec.indices = new uint8_t**[1];

vec.indices[0] = new uint8_t*[1];

vec.indices[0][0] = nullptr;

vec.vals = nullptr;

return vec;

taco_tensor_t mat;

mat.dimensions = new int32_t[2];

mat.dimensions[0] = M;

mat.dimensions[1] = N;

mat.indices = new uint8_t**[2];

mat.indices[1] = new uint8_t*[2];

mat.indices[1][0] = nullptr;

mat.indices[1][1] = nullptr;

mat.vals = nullptr;

return mat;

free(mat.indices[1][0]);

free(mat.indices[1][1]);

free(mat.vals);

int sparsity;

GxB_Vector_Option_get(*vec, GxB_SPARSITY_STATUS, &sparsity);

return (sparsity == GxB_BITMAP);

int sparsity;

GxB_Vector_Option_get(*vec, GxB_SPARSITY_STATUS, &sparsity);

return (sparsity == GxB_FULL);

int sparsity;

GxB_Format_Value fmt;

GxB_Matrix_Option_get(*mat, GxB_SPARSITY_STATUS, &sparsity);

GxB_Matrix_Option_get(*mat, GxB_FORMAT, &fmt);

return (sparsity == GxB_SPARSE && fmt == GxB_BY_ROW);

GrB_Index m;

int8_t* valid = nullptr;

void* vals = nullptr;

BitmapArrays vec;

vec.m = vector.dimensions[0];

vec.valid = (int8_t*)vector.indices[0][0];

vec.vals = vector.vals;

return vec;

if (a.m != b.m) {

return std::numeric_limits<double>::infinity();

}

double* avals = (double*)a.vals;

double* bvals = (double*)b.vals;

double ret = 0.0;

for (int i = 0 ; i < a.m; ++i) {

bool avalid = a.valid[i] && !std::isinf(avals[i]);

bool bvalid = b.valid[i] && !std::isinf(bvals[i]);

if (avalid != bvalid) {

return std::numeric_limits<double>::infinity();

}

if (avalid) {

if (avals[i] != 0.0) {

const double diff = std::abs(bvals[i] / avals[i] - 1.0);

if (diff > ret) {

ret = diff;

}

} else if (bvals[i] != 0.0) {

return std::numeric_limits<double>::infinity();

}

}

}

return ret;

if (a.m != b.m) {

return 1.0;

}

bool* avals = (bool*)a.vals;

bool* bvals = (bool*)b.vals;

for (int i = 0 ; i < a.m; ++i) {

bool avalid = a.valid[i] && avals[i];

bool bvalid = b.valid[i] && bvals[i];

if (avalid != bvalid) {

return 1.0;

}

}

return 0.0;

GrB_Index m, n;

GrB_Index* pos = nullptr;

GrB_Index* crd = nullptr;

void* vals = nullptr;

CSRArrays csr;

csr.m = matrix.dimensions[0];

csr.n = matrix.dimensions[1];

csr.pos = (GrB_Index*)matrix.indices[1][0];

csr.crd = (GrB_Index*)matrix.indices[1][1];

csr.vals = matrix.vals;

return csr;

//std::cout << a.m << " " << b.m << " " << a.n << " " << b.n << std::endl;

if (a.m != b.m || a.n != b.n) {

return std::numeric_limits<double>::infinity();

}

double* avals = (double*)a.vals;

double* bvals = (double*)b.vals;

double ret = 0.0;

for (int i = 0; i < a.m; ++i) {

int pA = a.pos[i];

int pB = b.pos[i];

while (pA < a.pos[i + 1] && pB < b.pos[i + 1]) {

while (pA < a.pos[i + 1] && std::isinf(avals[pA])) pA++;

while (pB < b.pos[i + 1] && std::isinf(bvals[pB])) pB++;

if (pA < a.pos[i + 1] && pB < b.pos[i + 1]) {

//std::cout << a.crd[pA] << " " << b.crd[pB] << " " << avals[pA] << " " << bvals[pB] << std::endl;

if (a.crd[pA] != b.crd[pB]) {

return std::numeric_limits<double>::infinity();

} else if (avals[pA] != 0.0) {

const double diff = std::abs(bvals[pB] / avals[pA] - 1.0);

if (diff > ret) {

ret = diff;

//std::cout << i << " " << a.crd[pA] << " " << b.crd[pB] << " " << avals[pA] << " " << bvals[pB] << std::endl;

}

} else if (bvals[pB] != 0.0) {

return std::numeric_limits<double>::infinity();

}

pA++;

pB++;

}

}

while (pA < a.pos[i + 1] && std::isinf(avals[pA])) pA++;

while (pB < b.pos[i + 1] && std::isinf(bvals[pB])) pB++;

if (pA != a.pos[i + 1] || pB != b.pos[i + 1]) {

return std::numeric_limits<double>::infinity();

}

}

return ret;

//std::cout << a.m << " " << b.m << " " << a.n << " " << b.n << std::endl;

if (a.m != b.m || a.n != b.n) {

return 1.0;

}

bool* avals = (bool*)a.vals;

bool* bvals = (bool*)b.vals;

for (int i = 0; i < a.m; ++i) {

int pA = a.pos[i];

int pB = b.pos[i];

while (pA < a.pos[i + 1] && pB < b.pos[i + 1]) {

while (pA < a.pos[i + 1] && !avals[pA]) pA++;

while (pB < b.pos[i + 1] && !bvals[pB]) pB++;

if (pA < a.pos[i + 1] && pB < b.pos[i + 1]) {

//std::cout << a.crd[pA] << " " << b.crd[pB] << " " << avals[pA] << " " << bvals[pB] << std::endl;

if (a.crd[pA] != b.crd[pB]) {

return 1.0;

}

pA++;

pB++;

}

}

while (pA < a.pos[i + 1] && !avals[pA]) pA++;

while (pB < b.pos[i + 1] && !bvals[pB]) pB++;

if (pA != a.pos[i + 1] || pB != b.pos[i + 1]) {

return 1.0;

}

}

return 0.0;

GraphBLASFixture() {

//const auto path = "/data/scratch/s3chou/formats-bench/data/pwtk.mtx";

const auto path = "/data/scratch/s3chou/formats-bench/data/webbase_1M.mtx";

//const auto path = "/data/scratch/s3chou/formats-bench/data/coPapersDBLP/coPapersDBLP.mtx";

//const auto path = "/data/scratch/changwan/florida_all/soc-LiveJournal1/soc-LiveJournal1.mtx";

//const auto path = "/data/scratch/changwan/florida_all/com-LiveJournal/com-LiveJournal.mtx";

//const auto path = "/data/scratch/changwan/florida_all/indochina-2004/indochina-2004.mtx";

// TODO: Only run for square matrices

//double bsw[GxB_NBITMAP_SWITCH] = {0};

//GrB_init(GrB_BLOCKING);

GrB_init(GrB_NONBLOCKING);

//GxB_Global_Option_set(GxB_HYPER_SWITCH, 1.0);

//GxB_Global_Option_set(GxB_HYPER_SWITCH, GxB_NEVER_HYPER);

GxB_Global_Option_set(GxB_FORMAT, GxB_BY_ROW);

GxB_Global_Option_set(GxB_NTHREADS, nthreads);

//GxB_Global_Option_set(GxB_BITMAP_SWITCH, bsw);

//GxB_Global_Option_set(GxB_BURBLE, 1);

read_matrix(path);

if (is_bool) {

GrB_Vector_new(&x_gb, GrB_BOOL, N);

} else {

GrB_Vector_new(&x_gb, GrB_FP64, N);

}

GxB_Vector_Option_set(x_gb, GxB_SPARSITY_CONTROL, GxB_BITMAP);

//GrB_Vector_assign_FP64(x_gb, NULL, NULL, 1.0, GrB_ALL, N, NULL);

for (GrB_Index i = 0; i < M; i += 4) {

if (is_bool) {

GrB_Vector_setElement_BOOL(x_gb, 1, i);

} else {

GrB_Vector_setElement_FP64(x_gb, (double)i, i);

}

}

GrB_Vector_wait(&x_gb);

taco_uassert(is_bitmap_vector(&x_gb)) << "x is not bitmap";

GrB_Vector_new(&m_gb, GrB_BOOL, M);

GxB_Vector_Option_set(m_gb, GxB_SPARSITY_CONTROL, GxB_FULL);

GrB_Vector_assign_BOOL(m_gb, NULL, NULL, true, GrB_ALL, N, NULL);

for (GrB_Index i = 0; i < M; i += 4) {

GrB_Vector_setElement_BOOL(m_gb, false, i);

}

GrB_Vector_wait(&m_gb);

taco_uassert(is_dense_vector(&m_gb)) << "m is not dense";

//GrB_Index stride = (GrB_Index)std::sqrt(M);

//for (GrB_Index i = 0; i < M; i += stride) {

// indices.push_back(i);

//}

//indices_taco = indices_to_taco_tensor(indices.data(), indices.size());

}

void read_matrix(const std::string& matrix_path) {

std::fstream stream;

stream.open(matrix_path, std::fstream::in);

if (!stream) {

stream.close();

return;

}

std::string line;

std::getline(stream, line);

// Read Header

std::stringstream lineStream(line);

std::string head, type, formats, field, symmetry;

lineStream >> head >> type >> formats >> field >> symmetry;

assert(head=="%%MatrixMarket");

// type = [matrix tensor]

// formats = [coordinate array]

assert((type == "matrix") || (type == "tensor"));

// field = [real integer complex pattern]

bool isreal = false;

bool isint = false;

if (field == "complex") {

stream.close();

return;

} else if (field == "real") {

isreal = true;

} else if (field == "integer") {

isint = true;

}

// symmetry = [general symmetric skew-symmetric Hermitian]

if ((symmetry != "general") && (symmetry != "symmetric") &&

(symmetry != "skew-symmetric")) {

stream.close();

return;

}

const bool symm = ((symmetry == "symmetric") ||

(symmetry == "skew-symmetric"));

const bool skew = (symmetry == "skew-symmetric");

std::getline(stream, line);

// Skip comments at the top of the file

std::string token;

do {

std::stringstream lineStream(line);

lineStream >> token;

if (token[0] != '%') {

break;

}

} while (std::getline(stream, line));

// The first non-comment line is the header with dimensions

std::vector<GrB_Index> dimensions;

char* linePtr = (char*)line.data();

while (auto dimension = std::strtoull(linePtr, &linePtr, 10)) {

dimensions.push_back(dimension);

}

assert(dimensions.size() == 3);

nnz = dimensions[2];

GrB_Index* rows = (GrB_Index*)malloc(sizeof(GrB_Index) * nnz * (1 + symm));

GrB_Index* cols = (GrB_Index*)malloc(sizeof(GrB_Index) * nnz * (1 + symm));

double* fvals = !is_bool ? (double*)malloc(sizeof(double) * nnz * (1 + symm)) : nullptr;

bool* bvals = is_bool ? (bool*)malloc(sizeof(bool) * nnz * (1 + symm)) : nullptr;

for (nnz = 0; std::getline(stream, line); nnz++) {

//if (nnz % 10000000 == 0) std::cout << nnz << std::endl;

linePtr = (char*)line.data();

const GrB_Index i = strtoull(linePtr, &linePtr, 10) - 1;

const GrB_Index j = strtoull(linePtr, &linePtr, 10) - 1;

double fval = 1.0;

bool bval = true;;

if (isreal) {

fval = strtod(linePtr, &linePtr);

bval = (fval != 0.0);

} else if (isint) {

fval = strtoll(linePtr, &linePtr, 10);

bval = (fval != 0.0);

}

rows[nnz] = i;

cols[nnz] = j;

if (fvals) {

fvals[nnz] = fval;

}

if (bvals) {

bvals[nnz] = bval;

}

if (symm && i != j) {

nnz++;

if (skew) {

fval = -1.0 * fval;

}

rows[nnz] = j;

cols[nnz] = i;

if (fvals) {

fvals[nnz] = fval;

}

if (bvals) {

bvals[nnz] = bval;

}

}

}

stream.close();

GrB_Matrix_new(&A_gb, is_bool ? GrB_BOOL : GrB_FP64, dimensions[0], dimensions[1]);

GxB_Matrix_Option_set(A_gb, GxB_SPARSITY_CONTROL, GxB_SPARSE);

if (is_bool) {

GrB_Matrix_build_BOOL(A_gb, rows, cols, bvals, nnz, GrB_LOR);

} else {

GrB_Matrix_build_FP64(A_gb, rows, cols, fvals, nnz, GrB_PLUS_FP64);

}

GrB_Matrix_wait(&A_gb);

taco_uassert(is_csr_matrix(&A_gb)) << "A is not CSR";

free(rows);

free(cols);

free(fvals);

free(bvals);

M = dimensions[0];

N = dimensions[1];

}

bool is_bool = true;

bool validate = false;

GrB_Index M, N, nnz;

GrB_Matrix A_gb = nullptr;

GrB_Matrix C_gb = nullptr;

GrB_Vector x_gb = nullptr;

GrB_Vector m_gb = nullptr;

GrB_Vector y_gb = nullptr;

taco_tensor_t* A_taco_t = nullptr;

taco_tensor_t* x_taco_t = nullptr;

taco_tensor_t* m_taco_t = nullptr;

//std::vector<GrB_Index> indices;

//taco_tensor_t indices_taco;

GrB_Descriptor desc;

GrB_Descriptor_new(&desc);

GrB_Descriptor_set(desc, GrB_MASK, GrB_COMP);

GrB_Descriptor_set(desc, GrB_OUTP, GrB_REPLACE);

for (auto _ : state) {

state.PauseTiming();

GrB_Vector_free(&fixture.y_gb);

GrB_Vector_new(&fixture.y_gb, fixture.is_bool ? GrB_BOOL : GrB_FP64, fixture.M);

GxB_Vector_Option_set(fixture.y_gb, GxB_SPARSITY_CONTROL, GxB_BITMAP);

state.ResumeTiming();

if (fixture.is_bool) {

GrB_mxv(fixture.y_gb, fixture.m_gb, NULL, GrB_LOR_LAND_SEMIRING_BOOL, fixture.A_gb, fixture.x_gb, desc);

} else {

GrB_mxv(fixture.y_gb, fixture.m_gb, NULL, GrB_MIN_PLUS_SEMIRING_FP64, fixture.A_gb, fixture.x_gb, desc);

}

}

taco_uassert(is_bitmap_vector(&fixture.y_gb)) << "y is not bitmap";

if (!fixture.validate) {

GrB_Vector_free(&fixture.y_gb);

}

GrB_Descriptor desc;

GrB_Descriptor_new(&desc);

GrB_Descriptor_set(desc, GrB_OUTP, GrB_REPLACE);

GrB_Descriptor_set(desc, GxB_AxB_METHOD, GxB_AxB_GUSTAVSON);

for (auto _ : state) {

state.PauseTiming();

GrB_Matrix_free(&fixture.C_gb);

GrB_Matrix_new(&fixture.C_gb, fixture.is_bool ? GrB_BOOL : GrB_FP64, fixture.M, fixture.N);

GxB_Matrix_Option_set(fixture.C_gb, GxB_SPARSITY_CONTROL, GxB_SPARSE);

state.ResumeTiming();

if (fixture.is_bool) {

GrB_mxm(fixture.C_gb, NULL, NULL, GrB_LOR_LAND_SEMIRING_BOOL, fixture.A_gb, fixture.A_gb, desc);

} else {

GrB_mxm(fixture.C_gb, NULL, NULL, GrB_MIN_PLUS_SEMIRING_FP64, fixture.A_gb, fixture.A_gb, desc);

}

//GrB_Matrix_wait(&fixture.C_gb);

}

taco_uassert(is_csr_matrix(&fixture.C_gb)) << "C is not CSR";

if (!fixture.validate) {

GrB_Matrix_free(&fixture.C_gb);

}

//GrB_init(GrB_BLOCKING);

//GxB_Global_Option_set(GxB_HYPER_SWITCH, GxB_NEVER_HYPER);

//GxB_Global_Option_set(GxB_FORMAT, GxB_BY_ROW);

//GxB_Global_Option_set(GxB_NTHREADS, nthreads);

GrB_Descriptor desc;

GrB_Descriptor_new(&desc);

GrB_Descriptor_set(desc, GrB_OUTP, GrB_REPLACE);

GrB_Index* indices = fixture.indices.data();

GrB_Index size = fixture.indices.size();

GrB_Matrix C = NULL;

for (auto _ : state) {

state.PauseTiming();

GrB_Matrix_free(&C);

state.ResumeTiming();

//GrB_Matrix_new(&C, GrB_FP64, fixture.A_taco.getDimension(0), fixture.A_taco.getDimension(1));

GrB_Matrix_new(&C, GrB_FP64, fixture.indices.size(), fixture.indices.size());

GrB_Matrix_extract(C, NULL, NULL, fixture.A_gb, indices, size, indices, size, desc);

}

//GrB_Index nnz;

//GrB_Matrix_nvals(&nnz, C);

//std::cout << "nnz: " << nnz << std::endl;

GrB_Matrix_free(&C);

#if 0

GrB_Index y1_dimension = (GrB_Index)(y->dimensions[0]);

double* restrict y_vals = (double*)(y->vals);

GrB_Index A1_dimension = (GrB_Index)(A->dimensions[0]);

GrB_Index* restrict A2_pos = (GrB_Index*)(A->indices[1][0]);

GrB_Index* restrict A2_crd = (GrB_Index*)(A->indices[1][1]);

double* restrict A_vals = (double*)(A->vals);

GrB_Index x1_dimension = (GrB_Index)(x->dimensions[0]);

double* restrict x_vals = (double*)(x->vals);

GrB_Index m1_dimension = (GrB_Index)(m->dimensions[0]);

bool* restrict m_vals = (bool*)(m->vals);

//y_vals = (double*)calloc(y1_dimension, sizeof(double));

y_vals = (double*)malloc(sizeof(double) * y1_dimension);

//#pragma omp parallel for schedule(static) num_threads(nthreads)

#pragma omp parallel for schedule(dynamic, 256) num_threads(nthreads)

for (GrB_Index i = 0; i < x1_dimension; i++) {

if (!(m_vals[i] != 0)) {

double tj_val = INFINITY;

//double tj_val = 0.0;

for (GrB_Index jA = A2_pos[i]; jA < A2_pos[(i + 1)]; jA++) {

GrB_Index j = A2_crd[jA];

tj_val = fmin(tj_val,A_vals[jA] + x_vals[j]);

}

y_vals[i] = tj_val;

}

else {

y_vals[i] = INFINITY;

}

}

y->vals = (uint8_t*)y_vals;

return 0;

#else

GrB_Index y1_dimension = (GrB_Index)(y->dimensions[0]);

int8_t* restrict y1_valid = (int8_t*)(y->indices[0][0]);

double* restrict y_vals = (double*)(y->vals);

GrB_Index A1_dimension = (GrB_Index)(A->dimensions[0]);

GrB_Index* restrict A2_pos = (GrB_Index*)(A->indices[1][0]);

GrB_Index* restrict A2_crd = (GrB_Index*)(A->indices[1][1]);

double* restrict A_vals = (double*)(A->vals);

GrB_Index x1_dimension = (GrB_Index)(x->dimensions[0]);

int8_t* restrict x1_valid = (int8_t*)(x->indices[0][0]);

double* restrict x_vals = (double*)(x->vals);

GrB_Index m1_dimension = (GrB_Index)(m->dimensions[0]);

bool* restrict m_vals = (bool*)(m->vals);

y1_valid = (int8_t*)calloc(1, sizeof(int8_t) * y1_dimension);

int32_t y_capacity = y1_dimension;

y_vals = (double*)malloc(sizeof(double) * y_capacity);

#pragma omp parallel for schedule(dynamic, 256) num_threads(nthreads)

for (GrB_Index i = 0; i < m1_dimension; i++) {

if (!(m_vals[i] != 0)) {

//double tj_val = 0.0;

double tj_val = INFINITY;

for (GrB_Index jA = A2_pos[i]; jA < A2_pos[(i + 1)]; jA++) {

GrB_Index j = A2_crd[jA];

if (x1_valid[j] == 1) {

tj_val = fmin(tj_val,A_vals[jA] + x_vals[j]);

}

}

y_vals[i] = tj_val;

y1_valid[i] = 1;

}

}

y->indices[0][0] = (uint8_t*)(y1_valid);

y->vals = (uint8_t*)y_vals;

return 0;

#endif

GrB_Index y1_dimension = (GrB_Index)(y->dimensions[0]);

int8_t* restrict y1_valid = (int8_t*)(y->indices[0][0]);

bool* restrict y_vals = (bool*)(y->vals);

GrB_Index A1_dimension = (GrB_Index)(A->dimensions[0]);

GrB_Index* restrict A2_pos = (GrB_Index*)(A->indices[1][0]);

GrB_Index* restrict A2_crd = (GrB_Index*)(A->indices[1][1]);

bool* restrict A_vals = (bool*)(A->vals);

GrB_Index x1_dimension = (GrB_Index)(x->dimensions[0]);

int8_t* restrict x1_valid = (int8_t*)(x->indices[0][0]);

bool* restrict x_vals = (bool*)(x->vals);

GrB_Index m1_dimension = (GrB_Index)(m->dimensions[0]);

bool* restrict m_vals = (bool*)(m->vals);

y1_valid = (int8_t*)calloc(1, sizeof(int8_t) * y1_dimension);

int32_t y_capacity = y1_dimension;

y_vals = (bool*)malloc(sizeof(bool) * y_capacity);

#pragma omp parallel for schedule(dynamic, 256) num_threads(nthreads)

for (GrB_Index i = 0; i < m1_dimension; i++) {

if (!(m_vals[i] != 0)) {

bool tj_val = 0;

for (GrB_Index jA = A2_pos[i]; jA < A2_pos[(i + 1)]; jA++) {

GrB_Index j = A2_crd[jA];

if (x1_valid[j] == 1) {

tj_val = tj_val || A_vals[jA];

if (tj_val == 1) {

break;

}

}

}

y_vals[i] = tj_val;

y1_valid[i] = 1;

}

}

y->indices[0][0] = (uint8_t*)(y1_valid);

y->vals = (uint8_t*)y_vals;

return 0;

GrB_Index A1_dimension = (GrB_Index)(A->dimensions[0]);

GrB_Index* restrict A2_pos = (GrB_Index*)(A->indices[1][0]);

GrB_Index* restrict A2_crd = (GrB_Index*)(A->indices[1][1]);

double* restrict A_vals = (double*)(A->vals);

GrB_Index B1_dimension = (GrB_Index)(B->dimensions[0]);

GrB_Index* restrict B2_pos = (GrB_Index*)(B->indices[1][0]);

GrB_Index* restrict B2_crd = (GrB_Index*)(B->indices[1][1]);

double* restrict B_vals = (double*)(B->vals);

GrB_Index C1_dimension = (GrB_Index)(C->dimensions[0]);

GrB_Index C2_dimension = (GrB_Index)(C->dimensions[1]);

GrB_Index* restrict C2_pos = (GrB_Index*)(C->indices[1][0]);

GrB_Index* restrict C2_crd = (GrB_Index*)(C->indices[1][1]);

double* restrict C_vals = (double*)(C->vals);

#if 0

GrB_Index* restrict A2_nnz = 0;

A2_nnz = (GrB_Index*)malloc(sizeof(GrB_Index) * B1_dimension);

GrB_Index* restrict qw_index_list_all = 0;

qw_index_list_all = (GrB_Index*)malloc(sizeof(GrB_Index) * C2_dimension * omp_get_max_threads());

bool* restrict qw_already_set_all = (bool*)calloc(C2_dimension * omp_get_max_threads(), sizeof(bool));

#pragma omp parallel for schedule(dynamic, 128) num_threads(nthreads)

for (GrB_Index qi = 0; qi < B1_dimension; qi++) {

GrB_Index qw_index_list_size = 0;

GrB_Index* qw_index_list = qw_index_list_all + (C2_dimension * omp_get_thread_num());

bool* qw_already_set = qw_already_set_all + (C2_dimension * omp_get_thread_num());

for (GrB_Index qkB = B2_pos[qi]; qkB < B2_pos[(qi + 1)]; qkB++) {

GrB_Index qk = B2_crd[qkB];

for (GrB_Index qjC = C2_pos[qk]; qjC < C2_pos[(qk + 1)]; qjC++) {

GrB_Index qj = C2_crd[qjC];

if (!qw_already_set[qj]) {

qw_index_list[qw_index_list_size] = qj;

qw_already_set[qj] = 1;

qw_index_list_size++;

}

}

}

GrB_Index tqjA2_nnz_val = 0;

for (GrB_Index qw_index_locator = 0; qw_index_locator < qw_index_list_size; qw_index_locator++) {

GrB_Index qj = qw_index_list[qw_index_locator];

tqjA2_nnz_val += (GrB_Index)1;

qw_already_set[qj] = 0;

}

A2_nnz[qi] = tqjA2_nnz_val;

}

free(qw_index_list_all);

free(qw_already_set_all);

A2_pos = (GrB_Index*)malloc(sizeof(GrB_Index) * (A1_dimension + 1));

A2_pos[0] = 0;

for (GrB_Index i = 0; i < A1_dimension; i++) {

A2_pos[i + 1] = A2_pos[i] + A2_nnz[i];

}

std::cout << A2_pos[A1_dimension] << std::endl;

A2_crd = (GrB_Index*)malloc(sizeof(GrB_Index) * A2_pos[A1_dimension]);

A_vals = (double*)malloc(sizeof(double) * A2_pos[A1_dimension]);

double* restrict w_all = 0;

GrB_Index* restrict w_index_list_all = 0;

w_index_list_all = (GrB_Index*)malloc(sizeof(GrB_Index) * C2_dimension * omp_get_max_threads());

bool* restrict w_already_set_all = (bool*)calloc(C2_dimension * omp_get_max_threads(), sizeof(bool));

w_all = (double*)malloc(sizeof(double) * C2_dimension * omp_get_max_threads());

#pragma omp parallel for schedule(dynamic, 128) num_threads(nthreads)

for (GrB_Index i = 0; i < B1_dimension; i++) {

GrB_Index w_index_list_size = 0;

GrB_Index* w_index_list = w_index_list_all + (C2_dimension * omp_get_thread_num());

bool* w_already_set = w_already_set_all + (C2_dimension * omp_get_thread_num());

double* w = w_all + (C2_dimension * omp_get_thread_num());

for (GrB_Index kB = B2_pos[i]; kB < B2_pos[(i + 1)]; kB++) {

GrB_Index k = B2_crd[kB];

for (GrB_Index jC = C2_pos[k]; jC < C2_pos[(k + 1)]; jC++) {

GrB_Index j = C2_crd[jC];

if (!w_already_set[j]) {

w[j] = B_vals[kB] + C_vals[jC];

w_index_list[w_index_list_size] = j;

w_already_set[j] = 1;

w_index_list_size++;

}

else {

w[j] = fmin(w[j], B_vals[kB] + C_vals[jC]);

//w[j] = w[j] + B_vals[kB] * C_vals[jC]);

}

}

}

//qsort(w_index_list, w_index_list_size, sizeof(GrB_Index), cmp);

for (GrB_Index w_index_locator = 0; w_index_locator < w_index_list_size; w_index_locator++) {

GrB_Index j = w_index_list[w_index_locator];

GrB_Index pA2 = A2_pos[i];

A2_pos[i] = A2_pos[i] + 1;

A2_crd[pA2] = j;

A_vals[pA2] = w[j];

w_already_set[j] = 0;

}

}

free(w_index_list_all);

free(w_already_set_all);

free(w_all);

for (GrB_Index p = 0; p < A1_dimension; p++) {

A2_pos[A1_dimension - p] = A2_pos[((A1_dimension - p) - 1)];

}

A2_pos[0] = 0;

free(A2_nnz);

#else

GrB_Index* restrict A2_nnz = 0;

A2_nnz = (GrB_Index*)malloc(sizeof(GrB_Index) * B1_dimension);

GrB_Index* restrict qw_index_list_all = 0;

qw_index_list_all = (GrB_Index*)malloc(sizeof(GrB_Index) * (C2_dimension * omp_get_max_threads()));

bool* restrict qw_already_set_all = calloc((C2_dimension * omp_get_max_threads()), sizeof(bool));

//#pragma omp parallel for schedule(runtime)

#pragma omp parallel for schedule(dynamic, 128) num_threads(nthreads)

for (GrB_Index qi = 0; qi < B1_dimension; qi++) {

GrB_Index qw_index_list_all_size = 0;

GrB_Index* restrict qw_index_list = qw_index_list_all + C2_dimension * omp_get_thread_num();

//GrB_Index* restrict qw_index_list = qw_index_list_all + qw_index_list_all_size * omp_get_thread_num();

bool* restrict qw_already_set = qw_already_set_all + C2_dimension * omp_get_thread_num();

for (GrB_Index qkB = B2_pos[qi]; qkB < B2_pos[(qi + 1)]; qkB++) {

GrB_Index qk = B2_crd[qkB];

for (GrB_Index qjC = C2_pos[qk]; qjC < C2_pos[(qk + 1)]; qjC++) {

GrB_Index qj = C2_crd[qjC];

if (!qw_already_set[qj]) {

qw_index_list[qw_index_list_all_size] = qj;

qw_already_set[qj] = 1;

qw_index_list_all_size++;

}

}

}

GrB_Index tqjA2_nnz_val = 0;

for (GrB_Index qw_index_locator = 0; qw_index_locator < qw_index_list_all_size; qw_index_locator++) {

GrB_Index qj = qw_index_list[qw_index_locator];

tqjA2_nnz_val += (GrB_Index)1;

qw_already_set[qj] = 0;

}

A2_nnz[qi] = tqjA2_nnz_val;

}

free(qw_index_list_all);

free(qw_already_set_all);

A2_pos = (GrB_Index*)malloc(sizeof(GrB_Index) * (A1_dimension + 1));

A2_pos[0] = 0;

for (GrB_Index i = 0; i < A1_dimension; i++) {

A2_pos[i + 1] = A2_pos[i] + A2_nnz[i];

}

//std::cout << A2_pos[A1_dimension] << std::endl;

A2_crd = (GrB_Index*)malloc(sizeof(GrB_Index) * A2_pos[A1_dimension]);

A_vals = (double*)malloc(sizeof(double) * A2_pos[A1_dimension]);

double* restrict w_all = 0;

GrB_Index* restrict w_index_list_all = 0;

w_index_list_all = (GrB_Index*)malloc(sizeof(GrB_Index) * (C2_dimension * omp_get_max_threads()));

bool* restrict w_already_set_all = calloc((C2_dimension * omp_get_max_threads()), sizeof(bool));

w_all = (double*)malloc(sizeof(double) * (C2_dimension * omp_get_max_threads()));

//#pragma omp parallel for schedule(runtime)

#pragma omp parallel for schedule(dynamic, 128) num_threads(nthreads)

for (GrB_Index i = 0; i < B1_dimension; i++) {

GrB_Index w_index_list_all_size = 0;

double* restrict w = w_all + C2_dimension * omp_get_thread_num();

GrB_Index* restrict w_index_list = w_index_list_all + C2_dimension * omp_get_thread_num();

//GrB_Index* restrict w_index_list = w_index_list_all + w_index_list_all_size * omp_get_thread_num();

bool* restrict w_already_set = w_already_set_all + C2_dimension * omp_get_thread_num();

for (GrB_Index kB = B2_pos[i]; kB < B2_pos[(i + 1)]; kB++) {

GrB_Index k = B2_crd[kB];

for (GrB_Index jC = C2_pos[k]; jC < C2_pos[(k + 1)]; jC++) {

GrB_Index j = C2_crd[jC];

if (!w_already_set[j]) {

w[j] = B_vals[kB] + C_vals[jC];

//w[j] = B_vals[kB] * C_vals[jC];

w_index_list[w_index_list_all_size] = j;

w_already_set[j] = 1;

w_index_list_all_size++;

}

else {

w[j] = fmin(w[j], B_vals[kB] + C_vals[jC]);

//w[j] = w[j] + B_vals[kB] * C_vals[jC];

}

}

}

//qsort(w_index_list, w_index_list_all_size, sizeof(GrB_Index), cmp);

for (GrB_Index w_index_locator = 0; w_index_locator < w_index_list_all_size; w_index_locator++) {

GrB_Index j = w_index_list[w_index_locator];

GrB_Index pA2 = A2_pos[i];

A2_pos[i] = A2_pos[i] + 1;

A2_crd[pA2] = j;

A_vals[pA2] = w[j];

w_already_set[j] = 0;

}

}

free(w_index_list_all);

free(w_already_set_all);

free(w_all);

for (GrB_Index p = 0; p < A1_dimension; p++) {

A2_pos[A1_dimension - p] = A2_pos[((A1_dimension - p) - 1)];

}

A2_pos[0] = 0;

free(A2_nnz);

#endif

A->indices[1][0] = (uint8_t*)(A2_pos);

A->indices[1][1] = (uint8_t*)(A2_crd);

A->vals = (uint8_t*)A_vals;

return 0;

GrB_Index A1_dimension = (GrB_Index)(A->dimensions[0]);

GrB_Index* restrict A2_pos = (GrB_Index*)(A->indices[1][0]);

GrB_Index* restrict A2_crd = (GrB_Index*)(A->indices[1][1]);

bool* restrict A_vals = (bool*)(A->vals);

GrB_Index B1_dimension = (GrB_Index)(B->dimensions[0]);

GrB_Index* restrict B2_pos = (GrB_Index*)(B->indices[1][0]);

GrB_Index* restrict B2_crd = (GrB_Index*)(B->indices[1][1]);

bool* restrict B_vals = (bool*)(B->vals);

GrB_Index C1_dimension = (GrB_Index)(C->dimensions[0]);

GrB_Index C2_dimension = (GrB_Index)(C->dimensions[1]);

GrB_Index* restrict C2_pos = (GrB_Index*)(C->indices[1][0]);

GrB_Index* restrict C2_crd = (GrB_Index*)(C->indices[1][1]);

bool* restrict C_vals = (bool*)(C->vals);

GrB_Index* restrict A2_nnz = 0;

A2_nnz = (GrB_Index*)malloc(sizeof(GrB_Index) * B1_dimension);

GrB_Index* restrict qw_index_list_all = 0;

qw_index_list_all = (GrB_Index*)malloc(sizeof(GrB_Index) * (C2_dimension * omp_get_max_threads()));

bool* restrict qw_already_set_all = calloc((C2_dimension * omp_get_max_threads()), sizeof(bool));

//#pragma omp parallel for schedule(runtime)

#pragma omp parallel for schedule(dynamic, 128) num_threads(nthreads)

for (GrB_Index qi = 0; qi < B1_dimension; qi++) {

GrB_Index qw_index_list_all_size = 0;

GrB_Index* restrict qw_index_list = qw_index_list_all + C2_dimension * omp_get_thread_num();

//GrB_Index* restrict qw_index_list = qw_index_list_all + qw_index_list_all_size * omp_get_thread_num();

bool* restrict qw_already_set = qw_already_set_all + C2_dimension * omp_get_thread_num();

for (GrB_Index qkB = B2_pos[qi]; qkB < B2_pos[(qi + 1)]; qkB++) {

GrB_Index qk = B2_crd[qkB];

for (GrB_Index qjC = C2_pos[qk]; qjC < C2_pos[(qk + 1)]; qjC++) {

GrB_Index qj = C2_crd[qjC];

if (!qw_already_set[qj]) {

qw_index_list[qw_index_list_all_size] = qj;

qw_already_set[qj] = 1;

qw_index_list_all_size++;

}

}

}

GrB_Index tqjA2_nnz_val = 0;

for (GrB_Index qw_index_locator = 0; qw_index_locator < qw_index_list_all_size; qw_index_locator++) {

GrB_Index qj = qw_index_list[qw_index_locator];

tqjA2_nnz_val += (GrB_Index)1;

qw_already_set[qj] = 0;

}

A2_nnz[qi] = tqjA2_nnz_val;

}

free(qw_index_list_all);

free(qw_already_set_all);