1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
#include <mpi.h>
#include <iostream>
#include <algorithm>
#include <functional>
#include <string>
#include <vector>
#include "CombBLAS/CombBLAS.h"
using namespace combblas;
template <class IT>
struct KTipsSR
{
static IT id() { return static_cast<IT>(0); }
static bool returnedSAID() { return false; }
static MPI_Op mpi_op() { return MPI_LOR; }
static IT add(const IT& arg1, const IT& arg2) { return (arg1 || arg2); }
static IT multiply(const IT& arg1, const IT& arg2) { return (arg1 && arg2); }
static void axpy(IT a, const IT& x, IT& y) { y = add(y, multiply(a, x)); }
};
template <class IT, class NT, class DER>
FullyDistVec<IT,IT> LastNzRowIdxPerCol(const SpParMat<IT,NT,DER>& A)
{
std::shared_ptr<CommGrid> grid = A.getcommgrid();
int myrank = grid->GetRank();
int myproccol = grid->GetRankInProcRow();
int myprocrow = grid->GetRankInProcCol();
MPI_Comm ColWorld = grid->GetColWorld();
IT total_rows = A.getnrow();
IT total_cols = A.getncol();
int procrows = grid->GetGridRows();
int proccols = grid->GetGridCols();
IT rows_perproc = total_rows / procrows;
IT cols_perproc = total_cols / proccols;
IT row_offset = myprocrow * rows_perproc;
IT col_offset = myproccol * cols_perproc;
DER *spSeq = A.seqptr();
IT localcols = spSeq->getncol();
std::vector<IT> local_colidx(localcols, static_cast<IT>(-1));
for (auto colit = spSeq->begcol(); colit != spSeq->endcol(); ++colit)
{
auto nzit = spSeq->begnz(colit);
if (nzit != spSeq->endnz(colit))
local_colidx[colit.colid()] = nzit.rowid() + row_offset;
}
MPI_Allreduce(MPI_IN_PLACE, local_colidx.data(), static_cast<int>(localcols), MPIType<IT>(), MPI_MAX, ColWorld);
std::vector<IT> fillarr;
if (!myprocrow)
for (auto itr = local_colidx.begin(); itr != local_colidx.end(); ++itr)
fillarr.push_back(*itr);
return FullyDistVec<IT,IT>(fillarr, grid);
}
template <class IT, class NT, class DER>
SpParMat<IT,NT,DER> FrontierMat(const SpParMat<IT,NT,DER>& A, const FullyDistSpVec<IT,IT>& sources, const NT& initval)
{
FullyDistVec<IT,IT> ri = sources.FindInds([](int arg1) { return true; });
FullyDistVec<IT,IT> ci(A.getcommgrid());
ci.iota(sources.getnnz(), static_cast<IT>(0));
return SpParMat<IT,NT,DER>(A.getnrow(), sources.getnnz(), ri, ci, initval, false);
}
int main(int argc, char *argv[])
{
int myrank, nprocs;
MPI_Init(&argc, &argv);
MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
if (argc != 3)
{
if (!myrank)
std::cerr << "Usage: ./KTipsTest <Matrix> <l>" << std::endl;
MPI_Finalize();
return -1;
}
{
int l = atoi(argv[2]);
std::shared_ptr<CommGrid> fullWorld;
fullWorld.reset(new CommGrid(MPI_COMM_WORLD, 0, 0));
SpParMat<int,int, SpDCCols<int,int>> A(fullWorld);
A.ParallelReadMM(std::string(argv[1]), false, maximum<int>());
FullyDistVec<int,int> D = A.Reduce(Column, std::plus<int>(), static_cast<int>(0));
FullyDistSpVec<int,int> R = D.Find(std::bind2nd(std::equal_to<int>(), static_cast<int>(1)));
SpParMat<int,int, SpDCCols<int,int>> F0 = FrontierMat(A, R, static_cast<int>(1));
SpParMat<int,int, SpDCCols<int,int>> F1 = PSpGEMM<KTipsSR<int>>(A, F0);
SpParMat<int,int, SpDCCols<int,int>> V = F0;
V += F1;
FullyDistVec<int,int> TipSources(A.getcommgrid(), F0.getncol(), static_cast<int>(-1));
FullyDistVec<int,int> TipDests(A.getcommgrid(), F0.getncol(), static_cast<int>(-1));
for (int k = 1; k <= l; ++k)
{
SpParMat<int,int, SpDCCols<int,int>> F2 = PSpGEMM<KTipsSR<int>>(A, F1);
F2.SetDifference(V);
V += F2;
FullyDistVec<int,int> Ns = F2.Reduce(Column, std::plus<int>(), static_cast<int>(0));
FullyDistSpVec<int,int> Tc = Ns.Find(std::bind2nd(std::greater_equal<int>(), static_cast<int>(2)));
FullyDistSpVec<int,int> Td = Ns.Find(std::bind2nd(std::not_equal_to<int>(), static_cast<int>(1)));
FullyDistVec<int,int> C0 = LastNzRowIdxPerCol(F0);
FullyDistVec<int,int> C1 = LastNzRowIdxPerCol(F1);
FullyDistSpVec<int,int> kSources = C0.GGet(Tc, [](const int arg1, const int arg2) { return arg2; }, static_cast<int>(-1));
FullyDistSpVec<int,int> kDests = C1.GGet(Tc, [](const int arg1, const int arg2) { return arg2; }, static_cast<int>(-1));
TipSources.Set(kSources);
TipDests.Set(kDests);
F1.PruneColumnByIndex(Td);
F2.PruneColumnByIndex(Td);
F0 = F1;
F1 = F2;
}
TipSources.DebugPrint();
TipDests.DebugPrint();
}
MPI_Finalize();
return 0;
}
|
