1//===- MatmulOptimizer.cpp -----------------------------------------------===//
2//
3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4// See https://llvm.org/LICENSE.txt for license information.
5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6//
7//===----------------------------------------------------------------------===//
8
9#include "polly/MatmulOptimizer.h"
10#include "polly/DependenceInfo.h"
11#include "polly/Options.h"
12#include "polly/ScheduleTreeTransform.h"
13#include "polly/ScopInfo.h"
14#include "polly/ScopPass.h"
15#include "polly/Simplify.h"
16#include "polly/Support/GICHelper.h"
17#include "polly/Support/ISLTools.h"
18#include "llvm/ADT/ArrayRef.h"
19#include "llvm/ADT/DenseSet.h"
20#include "llvm/ADT/Sequence.h"
21#include "llvm/ADT/SetOperations.h"
22#include "llvm/ADT/SmallVector.h"
23#include "llvm/ADT/StringRef.h"
24#include "llvm/ADT/iterator_range.h"
25#include "llvm/Analysis/TargetTransformInfo.h"
26#include "llvm/IR/DataLayout.h"
27#include "llvm/IR/Function.h"
28#include "llvm/IR/Module.h"
29#include "llvm/Support/CommandLine.h"
30#include "llvm/Support/Debug.h"
31#include "llvm/Support/TypeSize.h"
32#include "llvm/Support/raw_ostream.h"
33#include "isl/ctx.h"
34#include "isl/schedule_node.h"
35#include "isl/schedule_type.h"
36#include "isl/union_map.h"
37#include "isl/union_set.h"
38#include <algorithm>
39#include <cassert>
40#include <cmath>
41#include <cstdint>
42#include <string>
43#include <vector>
44
45#include "polly/Support/PollyDebug.h"
46#define DEBUG_TYPE "polly-opt-isl"
47
48using namespace llvm;
49using namespace polly;
50
51namespace llvm {
52class Value;
53}
54
55static cl::opt<int> LatencyVectorFma(
56 "polly-target-latency-vector-fma",
57 cl::desc("The minimal number of cycles between issuing two "
58 "dependent consecutive vector fused multiply-add "
59 "instructions."),
60 cl::Hidden, cl::init(Val: 8), cl::cat(PollyCategory));
61
62static cl::opt<int> ThroughputVectorFma(
63 "polly-target-throughput-vector-fma",
64 cl::desc("A throughput of the processor floating-point arithmetic units "
65 "expressed in the number of vector fused multiply-add "
66 "instructions per clock cycle."),
67 cl::Hidden, cl::init(Val: 1), cl::cat(PollyCategory));
68
69static cl::opt<int> FirstCacheLevelSize(
70 "polly-target-1st-cache-level-size",
71 cl::desc("The size of the first cache level specified in bytes."),
72 cl::Hidden, cl::init(Val: -1), cl::cat(PollyCategory));
73
74static cl::opt<int> FirstCacheLevelDefaultSize(
75 "polly-target-1st-cache-level-default-size",
76 cl::desc("The default size of the first cache level specified in bytes"
77 " (if not enough were provided by the TargetTransformInfo)."),
78 cl::Hidden, cl::init(Val: 32768), cl::cat(PollyCategory));
79
80static cl::opt<int> SecondCacheLevelSize(
81 "polly-target-2nd-cache-level-size",
82 cl::desc("The size of the second level specified in bytes."), cl::Hidden,
83 cl::init(Val: -1), cl::cat(PollyCategory));
84
85static cl::opt<int> SecondCacheLevelDefaultSize(
86 "polly-target-2nd-cache-level-default-size",
87 cl::desc("The default size of the second cache level specified in bytes"
88 " (if not enough were provided by the TargetTransformInfo)."),
89 cl::Hidden, cl::init(Val: 262144), cl::cat(PollyCategory));
90
91// This option, along with --polly-target-2nd-cache-level-associativity,
92// --polly-target-1st-cache-level-size, and --polly-target-2st-cache-level-size
93// represent the parameters of the target cache, which do not have typical
94// values that can be used by default. However, to apply the pattern matching
95// optimizations, we use the values of the parameters of Intel Core i7-3820
96// SandyBridge in case the parameters are not specified or not provided by the
97// TargetTransformInfo.
98static cl::opt<int> FirstCacheLevelAssociativity(
99 "polly-target-1st-cache-level-associativity",
100 cl::desc("The associativity of the first cache level."), cl::Hidden,
101 cl::init(Val: -1), cl::cat(PollyCategory));
102
103static cl::opt<int> FirstCacheLevelDefaultAssociativity(
104 "polly-target-1st-cache-level-default-associativity",
105 cl::desc("The default associativity of the first cache level"
106 " (if not enough were provided by the TargetTransformInfo)."),
107 cl::Hidden, cl::init(Val: 8), cl::cat(PollyCategory));
108
109static cl::opt<int> SecondCacheLevelAssociativity(
110 "polly-target-2nd-cache-level-associativity",
111 cl::desc("The associativity of the second cache level."), cl::Hidden,
112 cl::init(Val: -1), cl::cat(PollyCategory));
113
114static cl::opt<int> SecondCacheLevelDefaultAssociativity(
115 "polly-target-2nd-cache-level-default-associativity",
116 cl::desc("The default associativity of the second cache level"
117 " (if not enough were provided by the TargetTransformInfo)."),
118 cl::Hidden, cl::init(Val: 8), cl::cat(PollyCategory));
119
120static cl::opt<int> VectorRegisterBitwidth(
121 "polly-target-vector-register-bitwidth",
122 cl::desc("The size in bits of a vector register (if not set, this "
123 "information is taken from LLVM's target information."),
124 cl::Hidden, cl::init(Val: -1), cl::cat(PollyCategory));
125
126static cl::opt<int> PollyPatternMatchingNcQuotient(
127 "polly-pattern-matching-nc-quotient",
128 cl::desc("Quotient that is obtained by dividing Nc, the parameter of the"
129 "macro-kernel, by Nr, the parameter of the micro-kernel"),
130 cl::Hidden, cl::init(Val: 256), cl::cat(PollyCategory));
131
132static cl::opt<bool>
133 PMBasedTCOpts("polly-tc-opt",
134 cl::desc("Perform optimizations of tensor contractions based "
135 "on pattern matching"),
136 cl::init(Val: false), cl::ZeroOrMore, cl::cat(PollyCategory));
137
138static cl::opt<bool>
139 PMBasedMMMOpts("polly-matmul-opt",
140 cl::desc("Perform optimizations of matrix multiplications "
141 "based on pattern matching"),
142 cl::init(Val: true), cl::ZeroOrMore, cl::cat(PollyCategory));
143
144static cl::opt<int> OptComputeOut(
145 "polly-tc-dependences-computeout",
146 cl::desc("Bound the dependence analysis by a maximal amount of "
147 "computational steps (0 means no bound)"),
148 cl::Hidden, cl::init(Val: 500000), cl::ZeroOrMore, cl::cat(PollyCategory));
149
150namespace {
151/// Parameters of the micro kernel.
152///
153/// Parameters, which determine sizes of rank-1 (i.e., outer product) update
154/// used in the optimized matrix multiplication.
155struct MicroKernelParamsTy {
156 int Mr;
157 int Nr;
158};
159
160/// Parameters of the macro kernel.
161///
162/// Parameters, which determine sizes of blocks of partitioned matrices
163/// used in the optimized matrix multiplication.
164struct MacroKernelParamsTy {
165 int Mc;
166 int Nc;
167 int Kc;
168};
169
170/// Parameters of the matrix multiplication operands.
171///
172/// Parameters, which describe access relations that represent operands of the
173/// matrix multiplication.
174struct MatMulInfoTy {
175 MemoryAccess *A = nullptr;
176 MemoryAccess *B = nullptr;
177 MemoryAccess *ReadFromC = nullptr;
178 MemoryAccess *WriteToC = nullptr;
179 int i = -1;
180 int j = -1;
181 int k = -1;
182};
183
184/// Parameters of the tensor contraction operands.
185///
186/// A general d-dimensional tensor T ∈ R ^ Nu0 x ... x Nud−1 can be defined
187/// as the set of scalar elements indexed by the set of indices u0 ... ud,
188///
189/// T ≡ {Anu0...nud−1 ∈ R | (u0,...,ud−1) ∈ Nu0 x ... x Nud−1}.
190///
191/// Let A, B, and C be dA, dB, and dC-dimensional tensors, respectively.
192/// Let the free and the contracted indices of the tensor A be grouped into
193/// two bundles I = i0...ir−1 and P = p0...pt−1, respectively. Similarly,
194/// the free and the contracted indices of B are grouped into bundles
195/// J = j0..js−1 and P and the free indices of C are grouped into
196/// bundles I and J.
197///
198/// Tensor contraction (TC) of tensors A, B into tensor C can be represented as
199/// C(shuffle(I,J))=∑α·A(shuffle(I,P))·B(shuffle(P,J))+β·C(shuffle(I,J)),
200/// where ∑ is a summation over all contracted indices of P,
201/// α, β ∈ R, Npi is the length of the tensor dimension that corresponds
202/// to the index pi, A(shuffle(I, P)), B(shuffle(P, J)), C(shuffle(I, J)) are
203/// accesses to tensors A, B, C, respectively,
204/// shuffle(I, J), shuffle(I, P), and shuffle(P, J) are permutations of
205/// the enclosed indices.
206///
207/// Multiplication of C(shuffle(I,J)) by β can be moved into a different SCoP
208/// statement by loop distribution, which is done by the isl scheduler.
209// If β is not equal to one, the optimization of TC of Polly requires
210/// such a transformation.
211///
212/// TCInfoTy contains parameters, which describe access relations that represent
213/// operands of the tensor contraction.
214struct TCInfoTy {
215 /// @{
216 /// Memory accesses that represent reading from tensors, which are operands of
217 /// the tensor contraction.
218 MemoryAccess *A = nullptr;
219 MemoryAccess *B = nullptr;
220 /// @}
221
222 /// @{
223 /// Memory accesses that represent reading from and writing into the tensor,
224 /// which contains the result of the tensor contraction.
225 MemoryAccess *ReadFromC = nullptr;
226 MemoryAccess *WriteToC = nullptr;
227 /// @}
228
229 /// @{
230 /// Input dimensions of the schedule space, which represent free
231 /// indices of tensors.
232 SmallDenseSet<int> I;
233 SmallDenseSet<int> J;
234 /// @}
235
236 /// Input dimension of the schedule space, which represents contracted
237 /// indices of tensors.
238 SmallDenseSet<int> P;
239
240 /// @{
241 /// Sizes of tensor dimensions for corresponding input dimensions of
242 /// the schedule space. The size of the tensor dimension can be larger than
243 /// the size of the corresponding input dimension of the schedule space.
244 /// This does not correspond to a tensor contraction. However, such a pattern
245 /// will be optimized by the transformation.
246 SmallVector<int> DimensionSizes;
247 SmallVector<int> ADimensions;
248 SmallVector<int> BDimensions;
249 SmallVector<int> CDimensions;
250 /// @}
251
252 /// @{
253 /// Permutations of indices of I, J, and P, which describe operands of
254 /// the tensor contraction and its result.
255 SmallVector<int> OrderedI;
256 SmallVector<int> OrderedJ;
257 SmallVector<int> OrderedP;
258 /// @}
259};
260
261/// Create an isl::union_set, which describes the option of the form
262/// [isolate[] -> unroll[x]].
263///
264/// @param Ctx An isl::ctx, which is used to create the isl::union_set.
265static isl::union_set getUnrollIsolatedSetOptions(isl::ctx Ctx) {
266 isl::space Space = isl::space(Ctx, 0, 0, 1);
267 isl::map UnrollIsolatedSetOption = isl::map::universe(space: Space);
268 isl::id DimInId = isl::id::alloc(ctx: Ctx, name: "isolate", user: nullptr);
269 isl::id DimOutId = isl::id::alloc(ctx: Ctx, name: "unroll", user: nullptr);
270 UnrollIsolatedSetOption =
271 UnrollIsolatedSetOption.set_tuple_id(type: isl::dim::in, id: DimInId);
272 UnrollIsolatedSetOption =
273 UnrollIsolatedSetOption.set_tuple_id(type: isl::dim::out, id: DimOutId);
274 return UnrollIsolatedSetOption.wrap();
275}
276
277/// Permute the two dimensions of the isl map.
278///
279/// Permute @p DstPos and @p SrcPos dimensions of the isl map @p Map that
280/// have type @p DimType.
281///
282/// @param Map The isl map to be modified.
283/// @param DimType The type of the dimensions.
284/// @param DstPos The first dimension.
285/// @param SrcPos The second dimension.
286/// @return The modified map.
287static isl::map permuteDimensions(isl::map Map, isl::dim DimType,
288 unsigned DstPos, unsigned SrcPos) {
289 assert(DstPos < unsignedFromIslSize(Map.dim(DimType)) &&
290 SrcPos < unsignedFromIslSize(Map.dim(DimType)));
291 if (DstPos == SrcPos)
292 return Map;
293 isl::id DimId;
294 if (Map.has_tuple_id(type: DimType))
295 DimId = Map.get_tuple_id(type: DimType);
296 auto FreeDim = DimType == isl::dim::in ? isl::dim::out : isl::dim::in;
297 isl::id FreeDimId;
298 if (Map.has_tuple_id(type: FreeDim))
299 FreeDimId = Map.get_tuple_id(type: FreeDim);
300 auto MaxDim = std::max(a: DstPos, b: SrcPos);
301 auto MinDim = std::min(a: DstPos, b: SrcPos);
302 Map = Map.move_dims(dst_type: FreeDim, dst_pos: 0, src_type: DimType, src_pos: MaxDim, n: 1);
303 Map = Map.move_dims(dst_type: FreeDim, dst_pos: 0, src_type: DimType, src_pos: MinDim, n: 1);
304 Map = Map.move_dims(dst_type: DimType, dst_pos: MinDim, src_type: FreeDim, src_pos: 1, n: 1);
305 Map = Map.move_dims(dst_type: DimType, dst_pos: MaxDim, src_type: FreeDim, src_pos: 0, n: 1);
306 if (!DimId.is_null())
307 Map = Map.set_tuple_id(type: DimType, id: DimId);
308 if (!FreeDimId.is_null())
309 Map = Map.set_tuple_id(type: FreeDim, id: FreeDimId);
310 return Map;
311}
312
313/// Check the form of the access relation.
314///
315/// Check that the access relation @p AccMap has the form M[i][j], where i
316/// is a @p FirstPos and j is a @p SecondPos.
317///
318/// @param AccMap The access relation to be checked.
319/// @param FirstPos The index of the input dimension that is mapped to
320/// the first output dimension.
321/// @param SecondPos The index of the input dimension that is mapped to the
322/// second output dimension.
323/// @return True in case @p AccMap has the expected form and false,
324/// otherwise.
325static bool isMatMulOperandAcc(isl::set Domain, isl::map AccMap, int &FirstPos,
326 int &SecondPos) {
327 isl::space Space = AccMap.get_space();
328 isl::map Universe = isl::map::universe(space: Space);
329
330 if (unsignedFromIslSize(Size: Space.dim(type: isl::dim::out)) != 2)
331 return false;
332
333 // MatMul has the form:
334 // for (i = 0; i < N; i++)
335 // for (j = 0; j < M; j++)
336 // for (k = 0; k < P; k++)
337 // C[i, j] += A[i, k] * B[k, j]
338 //
339 // Permutation of three outer loops: 3! = 6 possibilities.
340 int FirstDims[] = {0, 0, 1, 1, 2, 2};
341 int SecondDims[] = {1, 2, 2, 0, 0, 1};
342 for (int i = 0; i < 6; i += 1) {
343 auto PossibleMatMul =
344 Universe.equate(type1: isl::dim::in, pos1: FirstDims[i], type2: isl::dim::out, pos2: 0)
345 .equate(type1: isl::dim::in, pos1: SecondDims[i], type2: isl::dim::out, pos2: 1);
346
347 AccMap = AccMap.intersect_domain(set: Domain);
348 PossibleMatMul = PossibleMatMul.intersect_domain(set: Domain);
349
350 // If AccMap spans entire domain (Non-partial write),
351 // compute FirstPos and SecondPos.
352 // If AccMap != PossibleMatMul here (the two maps have been gisted at
353 // this point), it means that the writes are not complete, or in other
354 // words, it is a Partial write and Partial writes must be rejected.
355 if (AccMap.is_equal(map2: PossibleMatMul)) {
356 if (FirstPos != -1 && FirstPos != FirstDims[i])
357 continue;
358 FirstPos = FirstDims[i];
359 if (SecondPos != -1 && SecondPos != SecondDims[i])
360 continue;
361 SecondPos = SecondDims[i];
362 return true;
363 }
364 }
365
366 return false;
367}
368
369/// Does the memory access represent a non-scalar operand of the matrix
370/// multiplication.
371///
372/// Check that the memory access @p MemAccess is the read access to a non-scalar
373/// operand of the matrix multiplication or its result.
374///
375/// @param MemAccess The memory access to be checked.
376/// @param MMI Parameters of the matrix multiplication operands.
377/// @return True in case the memory access represents the read access
378/// to a non-scalar operand of the matrix multiplication and
379/// false, otherwise.
380static bool isMatMulNonScalarReadAccess(MemoryAccess *MemAccess,
381 MatMulInfoTy &MMI) {
382 if (!MemAccess->isLatestArrayKind() || !MemAccess->isRead())
383 return false;
384 auto AccMap = MemAccess->getLatestAccessRelation();
385 isl::set StmtDomain = MemAccess->getStatement()->getDomain();
386 if (isMatMulOperandAcc(Domain: StmtDomain, AccMap, FirstPos&: MMI.i, SecondPos&: MMI.j) && !MMI.ReadFromC) {
387 MMI.ReadFromC = MemAccess;
388 return true;
389 }
390 if (isMatMulOperandAcc(Domain: StmtDomain, AccMap, FirstPos&: MMI.i, SecondPos&: MMI.k) && !MMI.A) {
391 MMI.A = MemAccess;
392 return true;
393 }
394 if (isMatMulOperandAcc(Domain: StmtDomain, AccMap, FirstPos&: MMI.k, SecondPos&: MMI.j) && !MMI.B) {
395 MMI.B = MemAccess;
396 return true;
397 }
398 return false;
399}
400
401/// Check accesses to operands of the matrix multiplication.
402///
403/// Check that accesses of the SCoP statement, which corresponds to
404/// the partial schedule @p PartialSchedule, are scalar in terms of loops
405/// containing the matrix multiplication, in case they do not represent
406/// accesses to the non-scalar operands of the matrix multiplication or
407/// its result.
408///
409/// @param PartialSchedule The partial schedule of the SCoP statement.
410/// @param MMI Parameters of the matrix multiplication operands.
411/// @return True in case the corresponding SCoP statement
412/// represents matrix multiplication and false,
413/// otherwise.
414static bool containsOnlyMatrMultAcc(isl::map PartialSchedule,
415 MatMulInfoTy &MMI) {
416 auto InputDimId = PartialSchedule.get_tuple_id(type: isl::dim::in);
417 auto *Stmt = static_cast<ScopStmt *>(InputDimId.get_user());
418 unsigned OutDimNum = unsignedFromIslSize(Size: PartialSchedule.range_tuple_dim());
419 assert(OutDimNum > 2 && "In case of the matrix multiplication the loop nest "
420 "and, consequently, the corresponding scheduling "
421 "functions have at least three dimensions.");
422 auto MapI =
423 permuteDimensions(Map: PartialSchedule, DimType: isl::dim::out, DstPos: MMI.i, SrcPos: OutDimNum - 1);
424 auto MapJ =
425 permuteDimensions(Map: PartialSchedule, DimType: isl::dim::out, DstPos: MMI.j, SrcPos: OutDimNum - 1);
426 auto MapK =
427 permuteDimensions(Map: PartialSchedule, DimType: isl::dim::out, DstPos: MMI.k, SrcPos: OutDimNum - 1);
428
429 auto Accesses = getAccessesInOrder(Stmt&: *Stmt);
430 for (auto *MemA = Accesses.begin(); MemA != Accesses.end() - 1; MemA++) {
431 auto *MemAccessPtr = *MemA;
432 if (MemAccessPtr->isLatestArrayKind() && MemAccessPtr != MMI.WriteToC &&
433 !isMatMulNonScalarReadAccess(MemAccess: MemAccessPtr, MMI) &&
434 !(MemAccessPtr->isStrideZero(Schedule: MapI) &&
435 MemAccessPtr->isStrideZero(Schedule: MapJ) && MemAccessPtr->isStrideZero(Schedule: MapK)))
436 return false;
437 }
438 return true;
439}
440
441/// Check for dependencies corresponding to the matrix multiplication.
442///
443/// Check that there is only true dependence of the form
444/// S(..., k, ...) -> S(..., k + 1, …), where S is the SCoP statement
445/// represented by @p Schedule and k is @p Pos. Such a dependence corresponds
446/// to the dependency produced by the matrix multiplication.
447///
448/// @param Schedule The schedule of the SCoP statement.
449/// @param D The SCoP dependencies.
450/// @param Pos The parameter to describe an acceptable true dependence.
451/// In case it has a negative value, try to determine its
452/// acceptable value.
453/// @return True in case dependencies correspond to the matrix multiplication
454/// and false, otherwise.
455static bool containsOnlyMatMulDep(isl::map Schedule, const Dependences *D,
456 int &Pos) {
457 isl::union_map Dep = D->getDependences(Kinds: Dependences::TYPE_RAW);
458 isl::union_map Red = D->getDependences(Kinds: Dependences::TYPE_RED);
459 if (!Red.is_null())
460 Dep = Dep.unite(umap2: Red);
461 auto DomainSpace = Schedule.get_space().domain();
462 auto Space = DomainSpace.map_from_domain_and_range(range: DomainSpace);
463 auto Deltas = Dep.extract_map(space: Space).deltas();
464 int DeltasDimNum = unsignedFromIslSize(Size: Deltas.dim(type: isl::dim::set));
465 for (int i = 0; i < DeltasDimNum; i++) {
466 auto Val = Deltas.plain_get_val_if_fixed(type: isl::dim::set, pos: i);
467 Pos = Pos < 0 && Val.is_one() ? i : Pos;
468 if (Val.is_nan() || !(Val.is_zero() || (i == Pos && Val.is_one())))
469 return false;
470 }
471 if (DeltasDimNum == 0 || Pos < 0)
472 return false;
473 return true;
474}
475
476/// Check if the SCoP statement could probably be optimized with analytical
477/// modeling.
478///
479/// containsMatrMult tries to determine whether the following conditions
480/// are true:
481/// 1. The last memory access modeling an array, MA1, represents writing to
482/// memory and has the form S(..., i1, ..., i2, ...) -> M(i1, i2) or
483/// S(..., i2, ..., i1, ...) -> M(i1, i2), where S is the SCoP statement
484/// under consideration.
485/// 2. There is only one loop-carried true dependency, and it has the
486/// form S(..., i3, ...) -> S(..., i3 + 1, ...), and there are no
487/// loop-carried or anti dependencies.
488/// 3. SCoP contains three access relations, MA2, MA3, and MA4 that represent
489/// reading from memory and have the form S(..., i3, ...) -> M(i1, i3),
490/// S(..., i3, ...) -> M(i3, i2), S(...) -> M(i1, i2), respectively,
491/// and all memory accesses of the SCoP that are different from MA1, MA2,
492/// MA3, and MA4 have stride 0, if the innermost loop is exchanged with any
493/// of loops i1, i2 and i3.
494///
495/// @param PartialSchedule The PartialSchedule that contains a SCoP statement
496/// to check.
497/// @D The SCoP dependencies.
498/// @MMI Parameters of the matrix multiplication operands.
499static bool containsMatrMult(isl::map PartialSchedule, const Dependences *D,
500 MatMulInfoTy &MMI) {
501 auto InputDimsId = PartialSchedule.get_tuple_id(type: isl::dim::in);
502 auto *Stmt = static_cast<ScopStmt *>(InputDimsId.get_user());
503 if (Stmt->size() <= 1)
504 return false;
505
506 auto Accesses = getAccessesInOrder(Stmt&: *Stmt);
507 for (auto *MemA = Accesses.end() - 1; MemA != Accesses.begin(); MemA--) {
508 auto *MemAccessPtr = *MemA;
509 if (!MemAccessPtr->isLatestArrayKind())
510 continue;
511 if (!MemAccessPtr->isWrite())
512 return false;
513 auto AccMap = MemAccessPtr->getLatestAccessRelation();
514 if (!isMatMulOperandAcc(Domain: Stmt->getDomain(), AccMap, FirstPos&: MMI.i, SecondPos&: MMI.j))
515 return false;
516 MMI.WriteToC = MemAccessPtr;
517 break;
518 }
519
520 if (!containsOnlyMatMulDep(Schedule: PartialSchedule, D, Pos&: MMI.k))
521 return false;
522
523 if (!MMI.WriteToC || !containsOnlyMatrMultAcc(PartialSchedule, MMI))
524 return false;
525
526 if (!MMI.A || !MMI.B || !MMI.ReadFromC)
527 return false;
528 return true;
529}
530
531/// Permute two dimensions of the band node.
532///
533/// Permute FirstDim and SecondDim dimensions of the Node.
534///
535/// @param Node The band node to be modified.
536/// @param FirstDim The first dimension to be permuted.
537/// @param SecondDim The second dimension to be permuted.
538static isl::schedule_node permuteBandNodeDimensions(isl::schedule_node Node,
539 unsigned FirstDim,
540 unsigned SecondDim) {
541 assert(isl_schedule_node_get_type(Node.get()) == isl_schedule_node_band &&
542 (unsigned)isl_schedule_node_band_n_member(Node.get()) >
543 std::max(FirstDim, SecondDim));
544 auto PartialSchedule =
545 isl::manage(ptr: isl_schedule_node_band_get_partial_schedule(node: Node.get()));
546 auto PartialScheduleFirstDim = PartialSchedule.at(pos: FirstDim);
547 auto PartialScheduleSecondDim = PartialSchedule.at(pos: SecondDim);
548 PartialSchedule =
549 PartialSchedule.set_union_pw_aff(pos: SecondDim, el: PartialScheduleFirstDim);
550 PartialSchedule =
551 PartialSchedule.set_union_pw_aff(pos: FirstDim, el: PartialScheduleSecondDim);
552 Node = isl::manage(ptr: isl_schedule_node_delete(node: Node.release()));
553 return Node.insert_partial_schedule(schedule: PartialSchedule);
554}
555
556static isl::schedule_node
557createMicroKernel(isl::schedule_node Node,
558 MicroKernelParamsTy MicroKernelParams) {
559 Node = applyRegisterTiling(Node, TileSizes: {MicroKernelParams.Mr, MicroKernelParams.Nr},
560 DefaultTileSize: 1);
561 Node = Node.parent().parent();
562 return permuteBandNodeDimensions(Node, FirstDim: 0, SecondDim: 1).child(pos: 0).child(pos: 0);
563}
564
565/// Create the BLIS macro-kernel.
566///
567/// We create the BLIS macro-kernel by applying a combination of tiling
568/// of dimensions of the band node and interchanging of two innermost
569/// modified dimensions. The values of MacroKernelParams's fields are used
570/// as tile sizes.
571///
572/// @param Node The schedule node to be modified.
573/// @param MacroKernelParams Parameters of the macro kernel
574/// to be used as tile sizes.
575static isl::schedule_node
576createMacroKernel(isl::schedule_node Node,
577 MacroKernelParamsTy MacroKernelParams) {
578 assert(isl_schedule_node_get_type(Node.get()) == isl_schedule_node_band);
579 if (MacroKernelParams.Mc == 1 && MacroKernelParams.Nc == 1 &&
580 MacroKernelParams.Kc == 1)
581 return Node;
582 int DimOutNum = isl_schedule_node_band_n_member(node: Node.get());
583 std::vector<int> TileSizes(DimOutNum, 1);
584 TileSizes[DimOutNum - 3] = MacroKernelParams.Mc;
585 TileSizes[DimOutNum - 2] = MacroKernelParams.Nc;
586 TileSizes[DimOutNum - 1] = MacroKernelParams.Kc;
587 Node = tileNode(Node, Identifier: "1st level tiling", TileSizes, DefaultTileSize: 1);
588 Node = Node.parent().parent();
589 Node = permuteBandNodeDimensions(Node, FirstDim: DimOutNum - 2, SecondDim: DimOutNum - 1);
590 Node = permuteBandNodeDimensions(Node, FirstDim: DimOutNum - 3, SecondDim: DimOutNum - 1);
591
592 return Node.child(pos: 0).child(pos: 0);
593}
594
595/// Get the size of the widest type of the matrix multiplication operands
596/// in bytes, including alignment padding.
597///
598/// @param MMI Parameters of the matrix multiplication operands.
599/// @return The size of the widest type of the matrix multiplication operands
600/// in bytes, including alignment padding.
601static uint64_t getMatMulAlignTypeSize(const MatMulInfoTy &MMI) {
602 auto *S = MMI.A->getStatement()->getParent();
603 auto &DL = S->getFunction().getParent()->getDataLayout();
604 auto ElementSizeA = DL.getTypeAllocSize(Ty: MMI.A->getElementType());
605 auto ElementSizeB = DL.getTypeAllocSize(Ty: MMI.B->getElementType());
606 auto ElementSizeC = DL.getTypeAllocSize(Ty: MMI.WriteToC->getElementType());
607 return std::max(l: {ElementSizeA, ElementSizeB, ElementSizeC});
608}
609
610/// Get the size of the widest type of the matrix multiplication operands
611/// in bits.
612///
613/// @param MMI Parameters of the matrix multiplication operands.
614/// @return The size of the widest type of the matrix multiplication operands
615/// in bits.
616static uint64_t getMatMulTypeSize(const MatMulInfoTy &MMI) {
617 auto *S = MMI.A->getStatement()->getParent();
618 auto &DL = S->getFunction().getParent()->getDataLayout();
619 auto ElementSizeA = DL.getTypeSizeInBits(Ty: MMI.A->getElementType());
620 auto ElementSizeB = DL.getTypeSizeInBits(Ty: MMI.B->getElementType());
621 auto ElementSizeC = DL.getTypeSizeInBits(Ty: MMI.WriteToC->getElementType());
622 return std::max(l: {ElementSizeA, ElementSizeB, ElementSizeC});
623}
624
625/// Get parameters of the BLIS micro kernel.
626///
627/// We choose the Mr and Nr parameters of the micro kernel to be large enough
628/// such that no stalls caused by the combination of latencies and dependencies
629/// are introduced during the updates of the resulting matrix of the matrix
630/// multiplication. However, they should also be as small as possible to
631/// release more registers for entries of multiplied matrices.
632///
633/// @param TTI Target Transform Info.
634/// @param MMI Parameters of the matrix multiplication operands.
635/// @return The structure of type MicroKernelParamsTy.
636/// @see MicroKernelParamsTy
637static MicroKernelParamsTy getMicroKernelParams(const TargetTransformInfo *TTI,
638 const MatMulInfoTy &MMI) {
639 assert(TTI && "The target transform info should be provided.");
640
641 // Nvec - Number of double-precision floating-point numbers that can be hold
642 // by a vector register. Use 2 by default.
643 long RegisterBitwidth = VectorRegisterBitwidth;
644
645 if (RegisterBitwidth == -1)
646 RegisterBitwidth =
647 TTI->getRegisterBitWidth(K: TargetTransformInfo::RGK_FixedWidthVector);
648 auto ElementSize = getMatMulTypeSize(MMI);
649 assert(ElementSize > 0 && "The element size of the matrix multiplication "
650 "operands should be greater than zero.");
651 auto Nvec = RegisterBitwidth / ElementSize;
652 if (Nvec == 0)
653 Nvec = 2;
654 int Nr = ceil(x: sqrt(x: (double)(Nvec * LatencyVectorFma * ThroughputVectorFma)) /
655 Nvec) *
656 Nvec;
657 int Mr = ceil(x: (double)(Nvec * LatencyVectorFma * ThroughputVectorFma / Nr));
658 return {.Mr: Mr, .Nr: Nr};
659}
660
661/// Determine parameters of the target cache.
662///
663/// @param TTI Target Transform Info.
664static void getTargetCacheParameters(const llvm::TargetTransformInfo *TTI) {
665 auto L1DCache = llvm::TargetTransformInfo::CacheLevel::L1D;
666 auto L2DCache = llvm::TargetTransformInfo::CacheLevel::L2D;
667 if (FirstCacheLevelSize == -1) {
668 if (TTI->getCacheSize(Level: L1DCache))
669 FirstCacheLevelSize = TTI->getCacheSize(Level: L1DCache).value();
670 else
671 FirstCacheLevelSize = static_cast<int>(FirstCacheLevelDefaultSize);
672 }
673 if (SecondCacheLevelSize == -1) {
674 if (TTI->getCacheSize(Level: L2DCache))
675 SecondCacheLevelSize = TTI->getCacheSize(Level: L2DCache).value();
676 else
677 SecondCacheLevelSize = static_cast<int>(SecondCacheLevelDefaultSize);
678 }
679 if (FirstCacheLevelAssociativity == -1) {
680 if (TTI->getCacheAssociativity(Level: L1DCache))
681 FirstCacheLevelAssociativity =
682 TTI->getCacheAssociativity(Level: L1DCache).value();
683 else
684 FirstCacheLevelAssociativity =
685 static_cast<int>(FirstCacheLevelDefaultAssociativity);
686 }
687 if (SecondCacheLevelAssociativity == -1) {
688 if (TTI->getCacheAssociativity(Level: L2DCache))
689 SecondCacheLevelAssociativity =
690 TTI->getCacheAssociativity(Level: L2DCache).value();
691 else
692 SecondCacheLevelAssociativity =
693 static_cast<int>(SecondCacheLevelDefaultAssociativity);
694 }
695}
696
697/// Get parameters of the BLIS macro kernel.
698///
699/// During the computation of matrix multiplication, blocks of partitioned
700/// matrices are mapped to different layers of the memory hierarchy.
701/// To optimize data reuse, blocks should be ideally kept in cache between
702/// iterations. Since parameters of the macro kernel determine sizes of these
703/// blocks, there are upper and lower bounds on these parameters.
704///
705/// @param TTI Target Transform Info.
706/// @param MicroKernelParams Parameters of the micro-kernel
707/// to be taken into account.
708/// @param MMI Parameters of the matrix multiplication operands.
709/// @return The structure of type MacroKernelParamsTy.
710/// @see MacroKernelParamsTy
711/// @see MicroKernelParamsTy
712static MacroKernelParamsTy
713getMacroKernelParams(const llvm::TargetTransformInfo *TTI,
714 const MicroKernelParamsTy &MicroKernelParams,
715 const MatMulInfoTy &MMI) {
716 getTargetCacheParameters(TTI);
717 // According to www.cs.utexas.edu/users/flame/pubs/TOMS-BLIS-Analytical.pdf,
718 // it requires information about the first two levels of a cache to determine
719 // all the parameters of a macro-kernel. It also checks that an associativity
720 // degree of a cache level is greater than two. Otherwise, another algorithm
721 // for determination of the parameters should be used.
722 if (!(MicroKernelParams.Mr > 0 && MicroKernelParams.Nr > 0 &&
723 FirstCacheLevelSize > 0 && SecondCacheLevelSize > 0 &&
724 FirstCacheLevelAssociativity > 2 && SecondCacheLevelAssociativity > 2))
725 return {.Mc: 1, .Nc: 1, .Kc: 1};
726 // The quotient should be greater than zero.
727 if (PollyPatternMatchingNcQuotient <= 0)
728 return {.Mc: 1, .Nc: 1, .Kc: 1};
729 int Car = floor(
730 x: (FirstCacheLevelAssociativity - 1) /
731 (1 + static_cast<double>(MicroKernelParams.Nr) / MicroKernelParams.Mr));
732
733 // Car can be computed to be zero since it is floor to int.
734 // On Mac OS, division by 0 does not raise a signal. This causes negative
735 // tile sizes to be computed. Prevent division by Cac==0 by early returning
736 // if this happens.
737 if (Car == 0)
738 return {.Mc: 1, .Nc: 1, .Kc: 1};
739
740 auto ElementSize = getMatMulAlignTypeSize(MMI);
741 assert(ElementSize > 0 && "The element size of the matrix multiplication "
742 "operands should be greater than zero.");
743 int Kc = (Car * FirstCacheLevelSize) /
744 (MicroKernelParams.Mr * FirstCacheLevelAssociativity * ElementSize);
745 double Cac =
746 static_cast<double>(Kc * ElementSize * SecondCacheLevelAssociativity) /
747 SecondCacheLevelSize;
748 int Mc = floor(x: (SecondCacheLevelAssociativity - 2) / Cac);
749 int Nc = PollyPatternMatchingNcQuotient * MicroKernelParams.Nr;
750
751 assert(Mc > 0 && Nc > 0 && Kc > 0 &&
752 "Matrix block sizes should be greater than zero");
753 return {.Mc: Mc, .Nc: Nc, .Kc: Kc};
754}
755
756/// Create an access relation that is specific to
757/// the matrix multiplication pattern.
758///
759/// Create an access relation of the following form:
760/// [O0, O1, O2, O3, O4, O5, O6, O7, O8] -> [OI, O5, OJ]
761/// where I is @p FirstDim, J is @p SecondDim.
762///
763/// It can be used, for example, to create relations that helps to consequently
764/// access elements of operands of a matrix multiplication after creation of
765/// the BLIS micro and macro kernels.
766///
767/// @see ScheduleTreeOptimizer::createMicroKernel
768/// @see ScheduleTreeOptimizer::createMacroKernel
769///
770/// Subsequently, the described access relation is applied to the range of
771/// @p MapOldIndVar, that is used to map original induction variables to
772/// the ones, which are produced by schedule transformations. It helps to
773/// define relations using a new space and, at the same time, keep them
774/// in the original one.
775///
776/// @param MapOldIndVar The relation, which maps original induction variables
777/// to the ones, which are produced by schedule
778/// transformations.
779/// @param FirstDim, SecondDim The input dimensions that are used to define
780/// the specified access relation.
781/// @return The specified access relation.
782static isl::map getMatMulAccRel(isl::map MapOldIndVar, unsigned FirstDim,
783 unsigned SecondDim) {
784 auto AccessRelSpace = isl::space(MapOldIndVar.ctx(), 0, 9, 3);
785 auto AccessRel = isl::map::universe(space: AccessRelSpace);
786 AccessRel = AccessRel.equate(type1: isl::dim::in, pos1: FirstDim, type2: isl::dim::out, pos2: 0);
787 AccessRel = AccessRel.equate(type1: isl::dim::in, pos1: 5, type2: isl::dim::out, pos2: 1);
788 AccessRel = AccessRel.equate(type1: isl::dim::in, pos1: SecondDim, type2: isl::dim::out, pos2: 2);
789 return MapOldIndVar.apply_range(map2: AccessRel);
790}
791
792static isl::schedule_node createExtensionNode(isl::schedule_node Node,
793 isl::map ExtensionMap) {
794 auto Extension = isl::union_map(ExtensionMap);
795 auto NewNode = isl::schedule_node::from_extension(extension: Extension);
796 return Node.graft_before(graft: NewNode);
797}
798
799static isl::schedule_node optimizePackedB(isl::schedule_node Node,
800 ScopStmt *Stmt, isl::map MapOldIndVar,
801 MicroKernelParamsTy MicroParams,
802 MacroKernelParamsTy MacroParams,
803 MatMulInfoTy &MMI) {
804 Scop *S = Stmt->getParent();
805 isl::set Domain = Stmt->getDomain();
806
807 // Create packed array.
808 unsigned FirstDimSize = MacroParams.Nc / MicroParams.Nr;
809 unsigned SecondDimSize = MacroParams.Kc;
810 unsigned ThirdDimSize = MicroParams.Nr;
811 ScopArrayInfo *PackedB =
812 S->createScopArrayInfo(ElementType: MMI.B->getElementType(), BaseName: "Packed_B",
813 Sizes: {FirstDimSize, SecondDimSize, ThirdDimSize});
814
815 // Compute the access relation for copying from B to PackedB.
816 isl::map AccRelB = MMI.B->getLatestAccessRelation();
817 isl::map AccRelPackedB = getMatMulAccRel(MapOldIndVar, FirstDim: 3, SecondDim: 7);
818 AccRelPackedB =
819 AccRelPackedB.set_tuple_id(type: isl::dim::out, id: PackedB->getBasePtrId());
820
821 // Create the copy statement and redirect access.
822 ScopStmt *CopyStmt = S->addScopStmt(SourceRel: AccRelB, TargetRel: AccRelPackedB, Domain);
823 MMI.B->setNewAccessRelation(AccRelPackedB);
824
825 unsigned Dim = unsignedFromIslSize(Size: MapOldIndVar.range_tuple_dim());
826 assert(Dim >= 2);
827 // Insert into the schedule tree.
828 isl::map ExtMap = MapOldIndVar.project_out(type: isl::dim::out, first: 2, n: Dim - 2);
829 ExtMap = ExtMap.reverse();
830 ExtMap = ExtMap.fix_si(type: isl::dim::out, pos: MMI.i, value: 0);
831 ExtMap = ExtMap.intersect_range(set: Domain);
832 ExtMap = ExtMap.set_tuple_id(type: isl::dim::out, id: CopyStmt->getDomainId());
833 return createExtensionNode(Node, ExtensionMap: ExtMap);
834}
835
836static isl::schedule_node optimizePackedA(isl::schedule_node Node, ScopStmt *,
837 isl::map MapOldIndVar,
838 MicroKernelParamsTy MicroParams,
839 MacroKernelParamsTy MacroParams,
840 MatMulInfoTy &MMI) {
841 isl::id InputDimsId = MapOldIndVar.get_tuple_id(type: isl::dim::in);
842 ScopStmt *Stmt = static_cast<ScopStmt *>(InputDimsId.get_user());
843 isl::set Domain = Stmt->getDomain();
844 isl::id DomainId = Domain.get_tuple_id();
845
846 // Create the packed array.
847 unsigned FirstDimSize = MacroParams.Mc / MicroParams.Mr;
848 unsigned SecondDimSize = MacroParams.Kc;
849 unsigned ThirdDimSize = MicroParams.Mr;
850 ScopArrayInfo *PackedA = Stmt->getParent()->createScopArrayInfo(
851 ElementType: MMI.A->getElementType(), BaseName: "Packed_A",
852 Sizes: {FirstDimSize, SecondDimSize, ThirdDimSize});
853
854 // Compute the access relation for copying from A to PackedA.
855 isl::map AccRelA = MMI.A->getLatestAccessRelation();
856 isl::map AccRelPackedA = getMatMulAccRel(MapOldIndVar, FirstDim: 4, SecondDim: 6);
857 AccRelPackedA =
858 AccRelPackedA.set_tuple_id(type: isl::dim::out, id: PackedA->getBasePtrId());
859 // { MemrefA[] -> PackedA[] }
860 isl::map PackedATranslator = AccRelPackedA.apply_domain(map2: AccRelA);
861
862 // Compute the domain for the copy statement.
863 // Construct the copy statement domain out of the 3 outermost scatter
864 // dimensions (to match the 3 band nodes surrounding the extension node) and
865 // the array elements to copy (one statement instance per array element).
866 // { Scatter[] }
867 isl::set ScatterDomain = MapOldIndVar.intersect_domain(set: Domain).range();
868 // { Scatter[] -> OutermostScatter[] }
869 isl::map OuterDomainMap =
870 makeIdentityMap(Set: ScatterDomain, RestrictDomain: true).project_out(type: isl::dim::out, first: 3, n: 6);
871 // { Scatter[] -> MemrefA[] }
872 isl::map CopyFrom = MapOldIndVar.reverse().apply_range(map2: AccRelA);
873 // { Scatter[] -> CopyStmt[] }
874 isl::map DomainTranslator = OuterDomainMap.range_product(map2: CopyFrom);
875 // { CopyStmt[] }
876 isl::set CopyDomain = DomainTranslator.range();
877
878 // Translate the access relations to the new domain.
879 // { CopyStmt[] -> MemrefA[] }
880 CopyFrom = CopyFrom.apply_domain(map2: DomainTranslator);
881 // { CopyStmt[] -> PackedA[] }
882 isl::map CopyTo = CopyFrom.apply_range(map2: PackedATranslator);
883
884 // Create the copy statement and redirect access.
885 ScopStmt *CopyStmt =
886 Stmt->getParent()->addScopStmt(SourceRel: CopyFrom, TargetRel: CopyTo, Domain: CopyDomain);
887 MMI.A->setNewAccessRelation(AccRelPackedA);
888
889 // Insert into the schedule tree.
890 // { Scatter[] -> CopyStmt[] }
891 isl::map ExtScatterCopy = makeIdentityMap(Set: CopyStmt->getDomain(), RestrictDomain: true);
892 ExtScatterCopy = ExtScatterCopy.project_out(type: isl::dim::in, first: 3, n: 2);
893 return createExtensionNode(Node, ExtensionMap: ExtScatterCopy);
894}
895
896/// Apply the packing transformation.
897///
898/// The packing transformation can be described as a data-layout
899/// transformation that requires to introduce a new array, copy data
900/// to the array, and change memory access locations to reference the array.
901/// It can be used to ensure that elements of the new array are read in-stride
902/// access, aligned to cache lines boundaries, and preloaded into certain cache
903/// levels.
904///
905/// As an example let us consider the packing of the array A that would help
906/// to read its elements with in-stride access. An access to the array A
907/// is represented by an access relation that has the form
908/// S[i, j, k] -> A[i, k]. The scheduling function of the SCoP statement S has
909/// the form S[i,j, k] -> [floor((j mod Nc) / Nr), floor((i mod Mc) / Mr),
910/// k mod Kc, j mod Nr, i mod Mr].
911///
912/// To ensure that elements of the array A are read in-stride access, we add
913/// a new array Packed_A[Mc/Mr][Kc][Mr] to the SCoP, using
914/// Scop::createScopArrayInfo, change the access relation
915/// S[i, j, k] -> A[i, k] to
916/// S[i, j, k] -> Packed_A[floor((i mod Mc) / Mr), k mod Kc, i mod Mr], using
917/// MemoryAccess::setNewAccessRelation, and copy the data to the array, using
918/// the copy statement created by Scop::addScopStmt.
919///
920/// @param Node The schedule node to be optimized.
921/// @param MapOldIndVar The relation, which maps original induction variables
922/// to the ones, which are produced by schedule
923/// transformations.
924/// @param MicroParams, MacroParams Parameters of the BLIS kernel
925/// to be taken into account.
926/// @param MMI Parameters of the matrix multiplication operands.
927/// @return The optimized schedule node.
928static isl::schedule_node
929optimizeDataLayoutMatrMulPattern(isl::schedule_node Node, isl::map MapOldIndVar,
930 MicroKernelParamsTy MicroParams,
931 MacroKernelParamsTy MacroParams,
932 MatMulInfoTy &MMI) {
933 isl::id InputDimsId = MapOldIndVar.get_tuple_id(type: isl::dim::in);
934 ScopStmt *Stmt = static_cast<ScopStmt *>(InputDimsId.get_user());
935
936 Node = Node.parent().parent().parent().parent().parent().parent();
937 Node = isl::manage(ptr: isl_schedule_node_band_split(node: Node.release(), pos: 2));
938
939 Node = Node.child(pos: 0);
940 Node =
941 optimizePackedB(Node, Stmt, MapOldIndVar, MicroParams, MacroParams, MMI);
942
943 Node = Node.child(pos: 0);
944 Node =
945 optimizePackedA(Node, Stmt, MapOldIndVar, MicroParams, MacroParams, MMI);
946
947 return Node.child(pos: 0).child(pos: 0).child(pos: 0).child(pos: 0).child(pos: 0);
948}
949
950/// Get a relation mapping induction variables produced by schedule
951/// transformations to the original ones.
952///
953/// @param Node The schedule node produced as the result of creation
954/// of the BLIS kernels.
955/// @param MicroKernelParams, MacroKernelParams Parameters of the BLIS kernel
956/// to be taken into account.
957/// @return The relation mapping original induction variables to the ones
958/// produced by schedule transformation.
959/// @see ScheduleTreeOptimizer::createMicroKernel
960/// @see ScheduleTreeOptimizer::createMacroKernel
961/// @see getMacroKernelParams
962static isl::map
963getInductionVariablesSubstitution(isl::schedule_node Node,
964 MicroKernelParamsTy MicroKernelParams,
965 MacroKernelParamsTy MacroKernelParams) {
966 auto Child = Node.child(pos: 0);
967 auto UnMapOldIndVar = Child.get_prefix_schedule_union_map();
968 auto MapOldIndVar = isl::map::from_union_map(umap: UnMapOldIndVar);
969 unsigned Dim = unsignedFromIslSize(Size: MapOldIndVar.range_tuple_dim());
970 if (Dim > 9u)
971 return MapOldIndVar.project_out(type: isl::dim::out, first: 0, n: Dim - 9);
972 return MapOldIndVar;
973}
974
975/// Isolate a set of partial tile prefixes and unroll the isolated part.
976///
977/// The set should ensure that it contains only partial tile prefixes that have
978/// exactly Mr x Nr iterations of the two innermost loops produced by
979/// the optimization of the matrix multiplication. Mr and Nr are parameters of
980/// the micro-kernel.
981///
982/// In case of parametric bounds, this helps to auto-vectorize the unrolled
983/// innermost loops, using the SLP vectorizer.
984///
985/// @param Node The schedule node to be modified.
986/// @param MicroKernelParams Parameters of the micro-kernel
987/// to be taken into account.
988/// @return The modified isl_schedule_node.
989static isl::schedule_node
990isolateAndUnrollMatMulInnerLoops(isl::schedule_node Node,
991 MicroKernelParamsTy MicroKernelParams) {
992 isl::schedule_node Child = Node.child(pos: 0);
993 isl::union_map UnMapOldIndVar = Child.get_prefix_schedule_relation();
994 isl::set Prefix = isl::map::from_union_map(umap: UnMapOldIndVar).range();
995 unsigned Dims = unsignedFromIslSize(Size: Prefix.tuple_dim());
996 assert(Dims >= 1);
997 Prefix = Prefix.project_out(type: isl::dim::set, first: Dims - 1, n: 1);
998 Prefix = getPartialTilePrefixes(ScheduleRange: Prefix, VectorWidth: MicroKernelParams.Nr);
999 Prefix = getPartialTilePrefixes(ScheduleRange: Prefix, VectorWidth: MicroKernelParams.Mr);
1000
1001 isl::union_set IsolateOption =
1002 getIsolateOptions(IsolateDomain: Prefix.add_dims(type: isl::dim::set, n: 3), OutDimsNum: 3);
1003 isl::ctx Ctx = Node.ctx();
1004 auto Options = IsolateOption.unite(uset2: getDimOptions(Ctx, Option: "unroll"));
1005 Options = Options.unite(uset2: getUnrollIsolatedSetOptions(Ctx));
1006 Node = Node.as<isl::schedule_node_band>().set_ast_build_options(Options);
1007 Node = Node.parent().parent().parent();
1008 IsolateOption = getIsolateOptions(IsolateDomain: Prefix, OutDimsNum: 3);
1009 Options = IsolateOption.unite(uset2: getDimOptions(Ctx, Option: "separate"));
1010 Node = Node.as<isl::schedule_node_band>().set_ast_build_options(Options);
1011 Node = Node.child(pos: 0).child(pos: 0).child(pos: 0);
1012 return Node;
1013}
1014
1015/// Insert "Loop Vectorizer Disabled" mark node.
1016///
1017/// @param Node The child of the mark node to be inserted.
1018/// @return The modified isl_schedule_node.
1019static isl::schedule_node markLoopVectorizerDisabled(isl::schedule_node Node) {
1020 auto Id = isl::id::alloc(ctx: Node.ctx(), name: "Loop Vectorizer Disabled", user: nullptr);
1021 return Node.insert_mark(mark: Id).child(pos: 0);
1022}
1023
1024/// Restore the initial ordering of dimensions of the band node
1025///
1026/// In case the band node represents all the dimensions of the iteration
1027/// domain, recreate the band node to restore the initial ordering of the
1028/// dimensions.
1029///
1030/// @param Node The band node to be modified.
1031/// @return The modified schedule node.
1032static isl::schedule_node
1033getBandNodeWithOriginDimOrder(isl::schedule_node Node) {
1034 assert(isl_schedule_node_get_type(Node.get()) == isl_schedule_node_band);
1035 if (isl_schedule_node_get_type(node: Node.child(pos: 0).get()) != isl_schedule_node_leaf)
1036 return Node;
1037 auto Domain = Node.get_universe_domain();
1038 assert(isl_union_set_n_set(Domain.get()) == 1);
1039 if (Node.get_schedule_depth().release() != 0 ||
1040 (unsignedFromIslSize(Size: isl::set(Domain).tuple_dim()) !=
1041 unsignedFromIslSize(Size: Node.as<isl::schedule_node_band>().n_member())))
1042 return Node;
1043 Node = isl::manage(ptr: isl_schedule_node_delete(node: Node.copy()));
1044 auto PartialSchedulePwAff = Domain.identity_union_pw_multi_aff();
1045 auto PartialScheduleMultiPwAff =
1046 isl::multi_union_pw_aff(PartialSchedulePwAff);
1047 PartialScheduleMultiPwAff =
1048 PartialScheduleMultiPwAff.reset_tuple_id(type: isl::dim::set);
1049 return Node.insert_partial_schedule(schedule: PartialScheduleMultiPwAff);
1050}
1051
1052static isl::schedule_node optimizeMatMulPattern(isl::schedule_node Node,
1053 const TargetTransformInfo *TTI,
1054 MatMulInfoTy &MMI) {
1055 assert(TTI && "The target transform info should be provided.");
1056 int DimOutNum = isl_schedule_node_band_n_member(node: Node.get());
1057 assert(DimOutNum > 2 && "In case of the matrix multiplication the loop nest "
1058 "and, consequently, the corresponding scheduling "
1059 "functions have at least three dimensions.");
1060 Node = getBandNodeWithOriginDimOrder(Node);
1061 Node = permuteBandNodeDimensions(Node, FirstDim: MMI.i, SecondDim: DimOutNum - 3);
1062 int NewJ = MMI.j == DimOutNum - 3 ? MMI.i : MMI.j;
1063 int NewK = MMI.k == DimOutNum - 3 ? MMI.i : MMI.k;
1064 Node = permuteBandNodeDimensions(Node, FirstDim: NewJ, SecondDim: DimOutNum - 2);
1065 NewK = NewK == DimOutNum - 2 ? NewJ : NewK;
1066 Node = permuteBandNodeDimensions(Node, FirstDim: NewK, SecondDim: DimOutNum - 1);
1067 auto MicroKernelParams = getMicroKernelParams(TTI, MMI);
1068 auto MacroKernelParams = getMacroKernelParams(TTI, MicroKernelParams, MMI);
1069 Node = createMacroKernel(Node, MacroKernelParams);
1070 Node = createMicroKernel(Node, MicroKernelParams);
1071 if (MacroKernelParams.Mc == 1 || MacroKernelParams.Nc == 1 ||
1072 MacroKernelParams.Kc == 1)
1073 return Node;
1074 auto MapOldIndVar = getInductionVariablesSubstitution(Node, MicroKernelParams,
1075 MacroKernelParams);
1076 if (MapOldIndVar.is_null())
1077 return Node;
1078 Node = markLoopVectorizerDisabled(Node: Node.parent()).child(pos: 0);
1079 Node = isolateAndUnrollMatMulInnerLoops(Node, MicroKernelParams);
1080 return optimizeDataLayoutMatrMulPattern(Node, MapOldIndVar, MicroParams: MicroKernelParams,
1081 MacroParams: MacroKernelParams, MMI);
1082}
1083
1084/// Check if this node contains a partial schedule that could
1085/// probably be optimized with analytical modeling.
1086///
1087/// isMatrMultPattern tries to determine whether the following conditions
1088/// are true:
1089/// 1. the partial schedule contains only one statement.
1090/// 2. there are exactly three input dimensions.
1091/// 3. all memory accesses of the statement will have stride 0 or 1, if we
1092/// interchange loops (switch the variable used in the inner loop to
1093/// the outer loop).
1094/// 4. all memory accesses of the statement except from the last one, are
1095/// read memory access and the last one is write memory access.
1096/// 5. all subscripts of the last memory access of the statement don't
1097/// contain the variable used in the inner loop.
1098/// If this is the case, we could try to use an approach that is similar to
1099/// the one used to get close-to-peak performance of matrix multiplications.
1100///
1101/// @param Node The node to check.
1102/// @param D The SCoP dependencies.
1103/// @param MMI Parameters of the matrix multiplication operands.
1104static bool isMatrMultPattern(isl::schedule_node Node, const Dependences *D,
1105 MatMulInfoTy &MMI) {
1106 auto PartialSchedule = isl::manage(
1107 ptr: isl_schedule_node_band_get_partial_schedule_union_map(node: Node.get()));
1108 if (isl_schedule_node_band_n_member(node: Node.get()) < 3 ||
1109 Node.get_schedule_depth().release() != 0 ||
1110 isl_union_map_n_map(umap: PartialSchedule.get()) != 1)
1111 return false;
1112 auto NewPartialSchedule = isl::map::from_union_map(umap: PartialSchedule);
1113 if (containsMatrMult(PartialSchedule: NewPartialSchedule, D, MMI))
1114 return true;
1115 return false;
1116}
1117
1118/// Get the dimension size.
1119///
1120/// Return the size of the dimension @p Pos, which is obtained from @p SAI.
1121/// Return -1 in the case of the first dimension of a multi-dimensional array,
1122/// since the ScopArrayInfo class does not carry size information.
1123///
1124/// @param SAI The information about the array.
1125/// @param Pos The position of the dimension.
1126/// @return The size of the dimension.
1127static int getDimSize(const ScopArrayInfo *SAI, unsigned Pos) {
1128 if (Pos == 0)
1129 return -1;
1130 const llvm::SCEV *SCEVDimSize = SAI->getDimensionSize(Dim: Pos);
1131 assert(SCEVDimSize);
1132 auto *ConstantDimSize = dyn_cast<const SCEVConstant>(Val: SCEVDimSize);
1133 assert(ConstantDimSize);
1134 auto *IntDimSize = dyn_cast<ConstantInt>(Val: ConstantDimSize->getValue());
1135 assert(IntDimSize);
1136 return IntDimSize->getSExtValue();
1137}
1138
1139/// Check whether the access relation has the specified form.
1140///
1141/// Check that the access relation @p AccMap has the form T[I0, …, In], where
1142/// indexes I0, …, In are specified by @p Dimensions.
1143///
1144/// @param Domain The domain of the access relation.
1145/// @param AccMap The access relation to be checked.
1146/// @param Dimensions The permutation of the subset of the input dimensions.
1147/// @return True if @p AccMap has the expected form and false,
1148/// otherwise.
1149static bool isCorrectAccessMap(isl::set Domain, isl::map AccMap,
1150 ArrayRef<int> Dimensions) {
1151 isl::space Space = AccMap.get_space();
1152 if (unsignedFromIslSize(Size: Space.dim(type: isl::dim::out)) != Dimensions.size())
1153 return false;
1154
1155 // Create an access relation of the following form:
1156 // [I0, …, Im] -> [Il, …, In], where indexes
1157 // Il, …, In are specified by @p Dimensions.
1158 isl::map PossibleTensor = isl::map::universe(space: Space);
1159 unsigned DimInSize = unsignedFromIslSize(Size: Space.dim(type: isl::dim::in));
1160 for (unsigned i = 0; i < Dimensions.size(); i++) {
1161 const int InPos = Dimensions[i];
1162 if ((InPos >= static_cast<int>(DimInSize)) || (InPos < 0))
1163 return false;
1164 PossibleTensor =
1165 PossibleTensor.equate(type1: isl::dim::in, pos1: InPos, type2: isl::dim::out, pos2: i);
1166 }
1167
1168 AccMap = AccMap.intersect_domain(set: Domain);
1169 PossibleTensor = PossibleTensor.intersect_domain(set: Domain);
1170
1171 // If AccMap != PossibleTensor here (the two maps have been gisted at
1172 // this point), it means that the writes are not complete, or in other
1173 // words, it is a Partial write and Partial writes must be rejected.
1174 return AccMap.is_equal(map2: PossibleTensor);
1175}
1176
1177/// Check whether the access represents the tensor contraction operand.
1178///
1179/// Check that the access relation @p AccMap has the form T[i1, …, in].
1180/// Obtained indexes i1, …, in, their sizes and their permutation are stored
1181/// into @p IndexSet, @p DimensionSizes, and @p Dimensions, respectively.
1182///
1183/// @param Domain The domain of the access relation.
1184/// @param AccMap The access relation to be checked.
1185/// @param IndexSet The subset of the input dimensions.
1186/// @param DimensionSizes Sizes of the input dimensions of @p Dimensions.
1187/// @param Dimensions The permutation of the subset of the input dimensions.
1188/// @return True if @p AccMap has the expected form and false,
1189/// otherwise.
1190static bool isTCOperandAcc(isl::set Domain, isl::map AccMap,
1191 SmallDenseSet<int> &IndexSet,
1192 SmallVectorImpl<int> &DimensionSizes,
1193 SmallVectorImpl<int> &Dimensions) {
1194 isl::id Id = AccMap.get_tuple_id(type: isl::dim::out);
1195 const ScopArrayInfo *SAI = ScopArrayInfo::getFromId(Id);
1196 assert(SAI && "AccMap should represent memory access");
1197
1198 // Fix values of output dimensions with respect to their positions.
1199 // In the case of the tensor contraction, values of output dimensions are
1200 // fixed and form a permutation of a subset of values of input dimensions.
1201 //
1202 // For example, in the case of Stmt[i][j][k] -> A[k][i], which represents
1203 // the operand of the tensor contraction, we get the following map by fixing
1204 // the output dimensions Stmt[1][j][0] -> A[0][1].
1205 //
1206 // We store the permutation of the subset of the input dimensions {2, 0} into
1207 // @p Dimensions.
1208 //
1209 // The obtained permutation and the isCorrectAccessMap function are used to
1210 // check whether the access relation @p AccMap represents the tensor
1211 // contraction operand. For example, in the case of
1212 // Stmt[i][j][k] -> A[i-1][j+1], we get Stmt[1][0][k] -> A[0][1] and,
1213 // consequently, {1, 0}, which is rejected by isCorrectAccessMap,
1214 // since it corresponds to Stmt[i][j][k] -> A[j][i].
1215 isl::map CheckMap = isl::manage(ptr: AccMap.copy());
1216 unsigned OutDimNum = unsignedFromIslSize(Size: CheckMap.dim(type: isl::dim::out));
1217 for (unsigned i = 0; i < OutDimNum; i++)
1218 CheckMap = CheckMap.fix_si(type: isl::dim::out, pos: i, value: i);
1219
1220 // Try to obtain the permutation and sizes of corresponding input dimensions.
1221 Dimensions.assign(NumElts: OutDimNum, Elt: -1);
1222 for (unsigned i : rangeIslSize(Begin: 0, End: CheckMap.dim(type: isl::dim::in))) {
1223 isl::val Val = getConstant(Map: CheckMap, Dim: isl::dim::in, Pos: i);
1224 if (!Val.is_int())
1225 continue;
1226 int OutPos = -1;
1227 llvm::APInt ValAPInt = APIntFromVal(V: Val);
1228 if (ValAPInt.isSignedIntN(N: 32))
1229 OutPos = ValAPInt.getSExtValue();
1230 if ((OutPos < 0) || (OutPos >= static_cast<int>(OutDimNum)) ||
1231 IndexSet.count(V: i))
1232 return false;
1233 IndexSet.insert(V: i);
1234 Dimensions[OutPos] = i;
1235 if (DimensionSizes[i] <= 0)
1236 DimensionSizes[i] = getDimSize(SAI, Pos: OutPos);
1237 }
1238
1239 return isCorrectAccessMap(Domain, AccMap, Dimensions);
1240}
1241
1242/// Find the intersection of two sets.
1243///
1244/// Find the intersection of the set @p A and the set @p B.
1245///
1246/// @param A, B Sets to intersect.
1247/// @return The set intersection.
1248static SmallDenseSet<int> intersect(const SmallDenseSet<int> &A,
1249 const SmallDenseSet<int> &B) {
1250 SmallDenseSet<int> Intersection = A;
1251 set_intersect(S1&: Intersection, S2: B);
1252 return Intersection;
1253}
1254
1255/// Check whether the set is a superset.
1256///
1257/// Check that the set @p A is a superset of @p B.
1258///
1259/// @param A, B Sets to be checked.
1260/// @return True if the set A is a superset of B.
1261static bool isSuperset(const SmallDenseSet<int> &A,
1262 const SmallDenseSet<int> &B) {
1263 return intersect(A, B).size() == B.size();
1264}
1265
1266/// Find the union of two sets.
1267///
1268/// Find the union of the set @p A and the set @p B.
1269///
1270/// @param A, B Sets to unite.
1271/// @return The set union.
1272static SmallDenseSet<int> unite(const SmallDenseSet<int> &A,
1273 const SmallDenseSet<int> &B) {
1274 SmallDenseSet<int> Union = A;
1275 set_union(S1&: Union, S2: B);
1276 return Union;
1277}
1278
1279/// Determine the access that writes to the tensor, which contains
1280/// the result of the tensor contraction.
1281///
1282/// @param Domain The domain of the statement.
1283/// @param Stmt The statement, which writes to memory.
1284/// @param TCI The information about the tensor contraction.
1285/// @param IandJIndexSet The set, which contains free indexes of tensors.
1286/// @return The determined MemoryAccess, or nullptr if there is no necessary
1287/// access within the SCoP.
1288static MemoryAccess *getWriteAccess(isl::set Domain, ScopStmt *Stmt,
1289 TCInfoTy &TCI,
1290 SmallDenseSet<int> &IandJIndexSet) {
1291 TCI.WriteToC = nullptr;
1292 SmallVector<MemoryAccess *, 32> Accesses = getAccessesInOrder(Stmt&: *Stmt);
1293 for (MemoryAccess *MemA : reverse(C&: Accesses)) {
1294 // A TC-like does not contain write scalar memory accesses
1295 if (!MemA->isLatestArrayKind())
1296 return nullptr;
1297 // The last memory access should be a write memory access.
1298 if (!MemA->isWrite())
1299 return nullptr;
1300
1301 isl::map AccMap = MemA->getLatestAccessRelation();
1302 if (!isTCOperandAcc(Domain, AccMap, IndexSet&: IandJIndexSet, DimensionSizes&: TCI.DimensionSizes,
1303 Dimensions&: TCI.CDimensions))
1304 return nullptr;
1305
1306 return MemA;
1307 }
1308 return nullptr;
1309}
1310
1311/// Determine an access, which reads elements of an operand of the tensor
1312/// contraction
1313///
1314/// @param MemAccessPtr The access, which reads elements of the tensor.
1315/// @param IndexSet The set, which contains indexes of the tensors.
1316/// @param IandJIndexSet The set, which contains free indexes of tensors.
1317/// @param Dimensions The permutation of the subset of the input dimensions.
1318/// @param TCI The information about the tensor contraction.
1319/// @return True if the memory access @p MemAccessPtr corresponds
1320/// to the tensor contraction.
1321static bool setReadAccess(MemoryAccess *MemAccessPtr,
1322 const SmallDenseSet<int> &IndexSet,
1323 const SmallDenseSet<int> &IandJIndexSet,
1324 ArrayRef<int> Dimensions, TCInfoTy &TCI) {
1325 if (!TCI.A) {
1326 // Probably IndexSet is a union of I and P sets.
1327 if (!isSuperset(A: IndexSet, B: TCI.P))
1328 return false;
1329
1330 // Obtain the set I.
1331 TCI.I = set_difference(S1: IndexSet, S2: TCI.P);
1332 if (!isSuperset(A: IandJIndexSet, B: TCI.I))
1333 return false;
1334
1335 // Obtain the set J.
1336 TCI.J = set_difference(S1: IandJIndexSet, S2: TCI.I);
1337
1338 // Set the first operand of the tensor contraction.
1339 TCI.A = MemAccessPtr;
1340 llvm::replace(Cont&: TCI.ADimensions, ContIt: TCI.ADimensions.begin(),
1341 ContEnd: TCI.ADimensions.end(), ValIt: Dimensions.begin(), ValEnd: Dimensions.end());
1342 return true;
1343 }
1344
1345 if (!TCI.B) {
1346 // IndexSet should be a union of J and P sets.
1347 if (unite(A: TCI.P, B: TCI.J) != IndexSet)
1348 return false;
1349
1350 // Set the second operand of the tensor contraction.
1351 TCI.B = MemAccessPtr;
1352 llvm::replace(Cont&: TCI.BDimensions, ContIt: TCI.BDimensions.begin(),
1353 ContEnd: TCI.BDimensions.end(), ValIt: Dimensions.begin(), ValEnd: Dimensions.end());
1354 return true;
1355 }
1356
1357 return false;
1358}
1359
1360/// Check that all memory accesses of the statement, except from the last
1361/// one, are read memory accesses, which read elements of operands of the tensor
1362/// contraction and its result.
1363///
1364/// @param Domain The domain of the statement.
1365/// @param Stmt The statement, which writes to memory.
1366/// @param TCI The information about the tensor contraction.
1367/// @param IandJIndexSet The set, which contains free indexes of tensors.
1368/// @return True if all read memory accesses of the statement @p Stmt correspond
1369/// to the tensor contraction.
1370static bool setReadAccesses(isl::set Domain, ScopStmt *Stmt, TCInfoTy &TCI,
1371 SmallDenseSet<int> &IandJIndexSet) {
1372 TCI.A = nullptr;
1373 TCI.B = nullptr;
1374 TCI.ReadFromC = nullptr;
1375 SmallVector<MemoryAccess *, 32> Accesses = getAccessesInOrder(Stmt&: *Stmt);
1376 for (auto *MemA = Accesses.begin(); *MemA != TCI.WriteToC; MemA++) {
1377 MemoryAccess *MemAccessPtr = *MemA;
1378
1379 // All memory accesses, except from the last one, should be read memory
1380 // accesses.
1381 if (MemAccessPtr->isWrite())
1382 return false;
1383
1384 isl::map AccMap = MemAccessPtr->getLatestAccessRelation();
1385
1386 if (!MemAccessPtr->isLatestArrayKind()) {
1387 // Check whether the scalar read memory access is not partial.
1388 if (!Domain.is_subset(set2: AccMap.domain()))
1389 return false;
1390 continue;
1391 return false;
1392 }
1393
1394 // There is only one memory access, which reads elements of the result of
1395 // the tensor contraction.
1396 if (AccMap.is_equal(map2: TCI.WriteToC->getLatestAccessRelation())) {
1397 if (TCI.ReadFromC)
1398 return false;
1399 TCI.ReadFromC = MemAccessPtr;
1400 continue;
1401 }
1402
1403 SmallVector<int> Dimensions;
1404 SmallDenseSet<int> IndexSet;
1405 if (!isTCOperandAcc(Domain, AccMap, IndexSet, DimensionSizes&: TCI.DimensionSizes,
1406 Dimensions))
1407 return false;
1408
1409 if (!setReadAccess(MemAccessPtr, IndexSet, IandJIndexSet, Dimensions, TCI))
1410 return false;
1411 }
1412
1413 // Check that there are read memory accesses, which read elements of operands
1414 // of the tensor contraction and its result.
1415 return TCI.ReadFromC && TCI.A && TCI.B;
1416}
1417
1418/// Check accesses to operands of the tensor contraction.
1419///
1420/// Check that accesses of the SCoP statement, which corresponds to
1421/// the partial schedule @p PartialSchedule, represent accesses
1422/// to the non-scalar operands of the tensor contraction.
1423///
1424/// @param Domain The domain of the SCoP statement.
1425/// @param PartialSchedule The partial schedule of the SCoP statement.
1426/// @param TCI Parameters of the tensor contraction operands.
1427/// @return True if the corresponding SCoP statement
1428/// represents tensor contraction and false,
1429/// otherwise.
1430static bool containsOnlyTCAcc(isl::set Domain, isl::map PartialSchedule,
1431 TCInfoTy &TCI) {
1432 isl::id InputDimsId = PartialSchedule.get_tuple_id(type: isl::dim::in);
1433 ScopStmt *Stmt = static_cast<ScopStmt *>(InputDimsId.get_user());
1434
1435 // In region statements, the order of memory accesses execution is not
1436 // predictable at compile-time.
1437 if ((Stmt->size() <= 1) || Stmt->isRegionStmt())
1438 return false;
1439
1440 unsigned DimNum = unsignedFromIslSize(Size: PartialSchedule.dim(type: isl::dim::in));
1441 TCI.DimensionSizes.resize(N: DimNum);
1442 SmallDenseSet<int> IandJIndexSet;
1443
1444 TCI.WriteToC = getWriteAccess(Domain, Stmt, TCI, IandJIndexSet);
1445 if (!TCI.WriteToC)
1446 return false;
1447
1448 if (intersect(A: IandJIndexSet, B: TCI.P).size() != 0)
1449 return false;
1450
1451 if (!setReadAccesses(Domain, Stmt, TCI, IandJIndexSet))
1452 return false;
1453
1454 return true;
1455}
1456
1457/// Check that dependency corresponds to the tensor contraction carried over
1458/// loop dimension @p Dim.
1459///
1460/// Check that the dependency has the form
1461/// S(..., ki, max(k(i + 1)), ..., max(kn), ...) ->
1462/// S(..., ki + 1, min(k(i + 1)), ..., min(kn), ...), where S is the SCoP
1463/// statement. For this purpose, we analyze the set @p DepDelta, which
1464/// represents the differences between image elements and domain elements of
1465/// the corresponding map.
1466///
1467/// @param DepDelta The set contains the differences between image elements
1468/// and corresponding domain elements of the map, which
1469/// represents the dependency.
1470/// @param Dim The position of the index ki.
1471/// @param BoundDeltas In the case of indexes of ki, the difference between
1472/// image elements and corresponding domain elements
1473/// corresponds to the difference between lexicographic
1474/// minimum and lexicographic maximum of the corresponding
1475/// dimension of the domain of the statement.
1476/// @param IndexSet Obtained indexes ki, which describe the dependency.
1477/// @return True if dependencies correspond to the tensor contraction
1478/// and false, otherwise.
1479static bool isReductionCarriedOverDim(isl::set DepDelta, unsigned Dim,
1480 isl::pw_multi_aff BoundDeltas,
1481 const SmallDenseSet<int> &IndexSet) {
1482 isl::space Space = DepDelta.get_space();
1483 isl::set Superset = isl::set::universe(space: Space);
1484 for (unsigned i = 0; i < Dim; i += 1)
1485 Superset = Superset.fix_si(type: isl::dim::set, pos: i, value: 0);
1486 Superset = Superset.fix_si(type: isl::dim::set, pos: Dim, value: 1);
1487
1488 // Check that the difference between the image element and the domain element
1489 // is equal to one in the case of the index ki. Image elements and
1490 // corresponding domain elements should be equal in the case of positions,
1491 // which are lower than the specified position.
1492 if (!DepDelta.is_subset(set2: Superset))
1493 return false;
1494
1495 // Compute a set, which is used to analyze how values of
1496 // the domain are related to the map that describes the dependency.
1497 isl_pw_multi_aff *DepDeltaPW = isl_pw_multi_aff_from_set(set: DepDelta.copy());
1498 BoundDeltas = BoundDeltas.add(pma2: isl::manage(ptr: DepDeltaPW));
1499 isl_set *ComplementRawSet = isl_set_from_pw_multi_aff(pma: BoundDeltas.release());
1500 isl::set Complement = isl::manage(ptr: ComplementRawSet);
1501
1502 for (unsigned i : rangeIslSize(Begin: Dim + 1, End: DepDelta.dim(type: isl::dim::set))) {
1503 if (!IndexSet.count(V: i)) {
1504 // Check the difference between the image element and the domain element
1505 // in the case of indexes, which do not describe the dependency.
1506 if (DepDelta.plain_get_val_if_fixed(type: isl::dim::set, pos: i).is_zero())
1507 continue;
1508 return false;
1509 }
1510
1511 // In the case of other indexes, which describe the dependency,
1512 // the difference between the image element and the domain element
1513 // should be equal to the difference between lexicographic minimum and
1514 // lexicographic maximum of the domain of the statement.
1515 if (!Complement.plain_get_val_if_fixed(type: isl::dim::set, pos: i).is_zero())
1516 return false;
1517 }
1518
1519 return true;
1520}
1521
1522/// Check whether dependencies are over the complete domain.
1523///
1524/// In the case of the tensor contraction RAW, WAW, WAR dependencies
1525/// have the form
1526/// S(..., ki, max(k(i + 1)), ..., max(kn), ...) ->
1527/// S(..., ki + 1, min(k(i + 1)), ..., min(kn), ...), where S is the SCoP
1528/// statement. Consequently, the domain of the dependencies
1529/// can be described as
1530/// Domain / Domain ∩ S(…, max(kn),…) ∩ S(…, max(k(i + 1)),…),
1531/// where Domain is the domain of the statement S.
1532///
1533/// For example, in the case of the following tensor contraction,
1534/// corresponding domains will have the following form.
1535///
1536/// An example of the tensor contraction:
1537/// for (i = 0; i < 1024; i++)
1538/// for (j = 0; j < 1024; j++)
1539/// for (l = 0; l < 64; ++l)
1540/// for (w = 0; w < 64; ++w)
1541/// C[i][j] += A[i][l][w] * B[w][j][l];
1542///
1543/// The domain of the statement:
1544/// { S[i0, i1, i2, i3] : i0 >= 0 and i0 <= 1023 and
1545/// i1 >= 0 and i1 <= 1023 and
1546/// i2 >= 0 and i2 <= 63 and
1547/// i3 >= 0 and i3 <= 63 }
1548///
1549/// The domain of the dependencies:
1550/// { S[i0, i1, i2, i3] : (i0 >= 0 and i0 <= 1023 and
1551/// i1 >= 0 and i1 <= 1023 and
1552/// i2 >= 0 and i2 <= 63 and
1553/// i3 >= 0 and i3 <= 62) or
1554/// (i3 = 63 and i0 >= 0 and i0 <= 1023 and
1555/// i1 >= 0 and i1 <= 1023 and
1556/// i2 >= 0 and i2 <= 62) }
1557///
1558/// @param Domain The domain of the statement.
1559/// @param DepsForStmt RAW and RED dependencies for the statement.
1560/// @param UpperBound The lexicographic maximum of the elements in
1561/// the @p Domain.
1562/// @param IndexSet Obtained indexes ki, which describe the dependencies.
1563/// @return True if dependencies are over the complete domain
1564/// and false, otherwise.
1565static bool areDepsOverCompleteDomain(isl::set Domain, isl::map DepsForStmt,
1566 isl::pw_multi_aff UpperBound,
1567 SmallDenseSet<int> &IndexSet) {
1568 isl_set *UpperBoundRawSet = isl_set_from_pw_multi_aff(pma: UpperBound.copy());
1569 isl::set UpperBoundSet = isl::manage(ptr: UpperBoundRawSet);
1570
1571 isl::set DomainRed = isl::manage(ptr: Domain.copy());
1572 for (const auto It : IndexSet) {
1573 isl::val FixedVal = UpperBoundSet.plain_get_val_if_fixed(type: isl::dim::set, pos: It);
1574 if (FixedVal.is_nan())
1575 return false;
1576 DomainRed = isl::manage(
1577 ptr: isl_set_fix_val(set: DomainRed.copy(), type: isl_dim_set, pos: It, v: FixedVal.release()));
1578 }
1579 return DepsForStmt.domain().intersect(set2: Domain).is_equal(
1580 set2: Domain.subtract(set2: DomainRed));
1581}
1582
1583/// Check that dependencies correspond to the tensor contraction.
1584///
1585/// Check that there are only true dependencies of the form
1586/// S(..., ki, max(k(i + 1)), ..., max(kn), ...) ->
1587/// S(..., ki + 1, min(k(i + 1)), ..., min(kn), ...), where S is the SCoP
1588/// statement represented by @p Schedule. Such dependencies are produced by
1589/// the tensor contraction. Obtained indexes ki are stored into @p IndexSet.
1590///
1591/// The form of anti and output dependencies is specified implicitly by
1592/// the form the SCoP statement, which is checked by subsequent analysis.
1593///
1594/// @param Schedule The schedule of the SCoP statement.
1595/// @param D The SCoP dependencies.
1596/// @param Domain The domain of the statement.
1597/// @param IndexSet Obtained indexes ki, which describe the dependencies.
1598/// @return True if dependencies correspond to the tensor contraction
1599/// and false, otherwise.
1600static bool containsOnlyTcDeps(isl::map Schedule, const Dependences *D,
1601 SmallDenseSet<int> &IndexSet, isl::set Domain) {
1602 IslMaxOperationsGuard MaxOpGuard(Schedule.ctx().get(), OptComputeOut);
1603
1604 isl::union_map Dep =
1605 D->getDependences(Kinds: Dependences::TYPE_RAW | Dependences::TYPE_RED);
1606
1607 isl::space DomainSpace = Schedule.get_space().domain();
1608 isl::space Space = DomainSpace.map_from_domain_and_range(range: DomainSpace);
1609 isl::map DepsForStmt = Dep.extract_map(space: Space);
1610 isl::set DepDeltas = DepsForStmt.deltas();
1611 isl::size DeltasDimNum = DepDeltas.dim(type: isl::dim::set);
1612 isl::pw_multi_aff LowerBound = Domain.lexmin_pw_multi_aff();
1613 isl::pw_multi_aff UpperBound = Domain.lexmax_pw_multi_aff();
1614 isl::pw_multi_aff BoundDeltas = UpperBound.sub(pma2: LowerBound);
1615
1616 for (int i : reverse(C: rangeIslSize(Begin: 0, End: DeltasDimNum))) {
1617 // In the case of the tensor contraction, the difference between image
1618 // elements and domain elements lies on a hyperplane where a dimension
1619 // has the fixed value one.
1620 isl::set Intersection = DepDeltas.fix_si(type: isl::dim::set, pos: i, value: 1);
1621 if (Intersection.is_empty())
1622 continue;
1623
1624 if (!isReductionCarriedOverDim(DepDelta: Intersection, Dim: i, BoundDeltas, IndexSet))
1625 return false;
1626
1627 IndexSet.insert(V: i);
1628 DepDeltas = DepDeltas.subtract(set2: Intersection);
1629 }
1630
1631 // In the case of the tensor contraction, all dependencies should have
1632 // the previously described form.
1633 if ((unsignedFromIslSize(Size: DeltasDimNum) == 0) || !DepDeltas.is_empty())
1634 return false;
1635
1636 return areDepsOverCompleteDomain(Domain, DepsForStmt, UpperBound, IndexSet);
1637}
1638
1639/// Check if the SCoP statement could probably be optimized with analytical
1640/// modeling.
1641///
1642/// containsTCInfoTy tries to determine whether the following conditions
1643/// are true:
1644///
1645/// 1. The last memory access modeling an array, MA1, represents writing to
1646/// memory and has the form S(..., I, ..., J, ...) -> M(shuffle(I, J)),
1647/// where S is the SCoP statement under consideration and shuffle(I, J)
1648/// is a permutation of indexes of sets I and J.
1649/// 2. There are only true dependencies of the form
1650/// S(..., ki, max(k(i + 1)), ..., max(kn), ...) ->
1651/// S(..., ki + 1, min(k(i + 1)), ..., min(kn), ...), where S is the SCoP
1652/// statement represented by @p Schedule and ki are indexes of the set P.
1653/// 3. SCoP contains an arbitrary number of reads from constants and only three
1654/// access relations, MA2, MA3, and MA4 that represent reading from memory
1655/// and have the form
1656/// S(..., I, ..., P, ...) -> M(shuffle(I, P)),
1657/// S(..., P, ..., J, ...) -> M(shuffle(J, P)),
1658/// S(...) -> M(shuffle(I, J)), respectively.
1659///
1660/// @param PartialSchedule The PartialSchedule that contains a SCoP statement
1661/// to check.
1662/// @param D The SCoP dependencies.
1663/// @param TCI Parameters of the tensor contraction operands.
1664/// @param Domain The domain of the statement.
1665/// @return True if dependencies and memory accesses correspond to the tensor
1666/// contraction and false, otherwise.
1667static bool containsTCInfoTy(isl::map PartialSchedule, const Dependences *D,
1668 TCInfoTy &TCI, isl::set Domain) {
1669 if (!containsOnlyTcDeps(Schedule: PartialSchedule, D, IndexSet&: TCI.P, Domain))
1670 return false;
1671
1672 // TODO: handle cases of scalar multiplication if needed.
1673 if (TCI.P.size() == 0)
1674 return false;
1675
1676 if (!containsOnlyTCAcc(Domain, PartialSchedule, TCI))
1677 return false;
1678
1679 // TODO: handle cases of GEMV if needed.
1680 if ((TCI.I.size() == 0) || (TCI.J.size() == 0))
1681 return false;
1682
1683 return true;
1684}
1685
1686/// Check if this node contains a partial schedule that could
1687/// probably be optimized with analytical modeling.
1688///
1689/// isTCPattern is used to determine whether the SCoP represents a TC-like
1690/// kernel [1], which is a perfectly nested set of loops, with a data usage
1691/// pattern that is similar to that produced by the tensor contraction.
1692///
1693/// A TC-like kernel can be defined as follows:
1694///
1695/// 1. It satisfies the requirements of the polyhedral model.
1696/// 2. Without loss of generality, it contains three nonempty bundles of
1697/// one-dimensional for-loops with induction variables that are grouped into
1698/// bundles I = i0...i(r-1), J = j0..j(s-1), and P = p0...p(t-1), and they
1699/// are incremented by one.
1700/// 3. The innermost loop body can be represented as a statement of the form
1701/// C(shuffle(I, J)) = E(A(shuffle(I, P)), B(shuffle(P, J)),
1702/// C(shuffle(I, J))), where A(shuffle(I, P)), B(shuffle(P, J)),
1703/// C(shuffle(I, J)) are accesses to tensors A, B, C, respectively,
1704/// shuffle(I, J), shuffle(I, P), and shuffle(P, J) are permutations of the
1705/// enclosed indices, and E is an expression that contains reads from
1706/// the tensors A, B, C, and an arbitrary number of reads from constants
1707/// with respect to bundles I, J, and P.
1708///
1709/// TC can be considered as a particular case of a TC-like kernel.
1710///
1711/// The order of loops with indexes from P should be preserved. Otherwise,
1712/// isTCPattern should check if a commutative operation is used.
1713///
1714/// isTCPattern performs the following steps to check whether the SCoP
1715/// corresponds to a definition of a TC-like kernel:
1716///
1717/// 1. Checks that the node is the innermost band node.
1718/// 2. Checks that the partial schedule contains only one statement.
1719/// 3. Check that all ancestors of the node contain all band nodes for
1720/// the statement and only mark nodes interleave such band nodes. This
1721/// corresponds to a straightforward implementation of TC.
1722/// 4. Analyses the dependencies to determine contraction dimensions.
1723/// 5. Check that the last memory access modeling an array, represents writing
1724/// to the result of the TC-like kernel.
1725/// 6. Check that SCoP contains only three access relations that represent
1726/// reading of the operands of the TC-like kernel and an arbitrary number of
1727/// reads from constants.
1728///
1729/// [1] - Gareev R., Grosser T., Kruse M. High-Performance Generalized Tensor
1730/// Operations: A Compiler-Oriented Approach // ACM Transactions
1731/// Architecture and Code Optimization (TACO). 2018.
1732/// Vol. 15, no. 3. P. 34:1–34:27. DOI: 10.1145/3235029.
1733///
1734/// If this is the case, we could logically represent tensors as matrices and
1735/// apply algorithms, which are used to get close-to-peak performance of
1736/// matrix multiplications in manually tuned BLAS libraries (e.g., BLIS).
1737///
1738/// @param Node The node to check.
1739/// @param D The SCoP dependencies.
1740/// @param TCI Parameters of the tensor contraction operands.
1741static bool isTCPattern(isl::schedule_node Node, const Dependences *D,
1742 TCInfoTy &TCI) {
1743 Node = Node.child(pos: 0);
1744 isl::union_map PartialSchedule = Node.get_prefix_schedule_union_map();
1745 isl::union_set Domain = Node.domain();
1746 Node = Node.parent();
1747
1748 // The partial schedule should contain only one statement.
1749 // TODO: This constraint should not be intrinsic to the algorithm.
1750 if (isl_union_set_n_set(uset: Domain.get()) != 1)
1751 return false;
1752
1753 isl_schedule_node_type NodeType = isl_schedule_node_get_type(node: Node.get());
1754
1755 // Check that all ancestors of the node contain all band nodes for
1756 // the statement, which represents the TC-like kernel, and only mark nodes
1757 // interleave such band nodes. This corresponds to a straightforward
1758 // implementation of TC with/without DeLICM applied.
1759 //
1760 // For example, this covers the matrix multiplication pattern after a full
1761 // run of -polly-optree and -polly-delicm, where the write access is not
1762 // through the original memory access, but trough a PHI node that was
1763 // delicmed. Subsequently, such band nodes will be replaced by a single band
1764 // node.
1765 //
1766 // The corresponding schedule can be the following, where Stmt_for_body8
1767 // contains the matrix multiplication:
1768 //
1769 // domain: "{ Stmt_for_body8[i0, i1, i2] : 0 <= i0 <= 1599 and
1770 // 0 <= i1 <= 1799 and
1771 // 0 <= i2 <= 2199;
1772 // Stmt_for_body3[i0, i1] : 0 <= i0 <= 1599 and
1773 // 0 <= i1 <= 1799;
1774 // Stmt_for_body3_last[i0, i1] : 0 <= i0 <= 1599 and
1775 // 0 <= i1 <= 1799 }"
1776 // child:
1777 // sequence:
1778 // - filter: "{ Stmt_for_body3[i0, i1] }"
1779 // child:
1780 // schedule: "[{ Stmt_for_body3[i0, i1] -> [(i0)] },
1781 // { Stmt_for_body3[i0, i1] -> [(i1)] }]"
1782 // permutable: 1
1783 // coincident: [ 1, 1 ]
1784 // - filter: "{ Stmt_for_body3_last[i0, i1] }"
1785 // child:
1786 // schedule: "[{ Stmt_for_body3_last[i0, i1] -> [(i0)] },
1787 // { Stmt_for_body3_last[i0, i1] -> [(i1)] }]"
1788 // permutable: 1
1789 // coincident: [ 1, 1 ]
1790 // - filter: "{ Stmt_for_body8[i0, i1, i2] }"
1791 // child:
1792 // schedule: "[{ Stmt_for_body8[i0, i1, i2] -> [(i0)] },
1793 // { Stmt_for_body8[i0, i1, i2] -> [(i1)] },
1794 // { Stmt_for_body8[i0, i1, i2] -> [(i2)] }]"
1795 // permutable: 1
1796 // coincident: [ 1, 1, 0 ]
1797 //
1798 while (NodeType != isl_schedule_node_domain) {
1799 if (NodeType == isl_schedule_node_filter) {
1800 if (!Node.parent().isa<isl::schedule_node_sequence>() ||
1801 !Node.parent().parent().isa<isl::schedule_node_domain>())
1802 return false;
1803 break;
1804 }
1805
1806 if ((NodeType != isl_schedule_node_band) &&
1807 (NodeType != isl_schedule_node_mark))
1808 return false;
1809
1810 Node = Node.parent();
1811 NodeType = isl_schedule_node_get_type(node: Node.get());
1812 }
1813
1814 isl::map PartialScheduleMap = isl::map::from_union_map(umap: PartialSchedule);
1815 if (containsTCInfoTy(PartialSchedule: PartialScheduleMap, D, TCI, Domain: isl::set(Domain)))
1816 return true;
1817
1818 return false;
1819}
1820
1821} // namespace
1822
1823isl::schedule_node
1824polly::tryOptimizeMatMulPattern(isl::schedule_node Node,
1825 const llvm::TargetTransformInfo *TTI,
1826 const Dependences *D) {
1827 TCInfoTy TCI;
1828 if (PMBasedTCOpts && isTCPattern(Node, D, TCI))
1829 POLLY_DEBUG(dbgs() << "The tensor contraction pattern was detected\n");
1830 MatMulInfoTy MMI;
1831 if (PMBasedMMMOpts && isMatrMultPattern(Node, D, MMI)) {
1832 POLLY_DEBUG(dbgs() << "The matrix multiplication pattern was detected\n");
1833 return optimizeMatMulPattern(Node, TTI, MMI);
1834 }
1835 return {};
1836}
1837

source code of polly/lib/Transform/MatmulOptimizer.cpp