WebSVN – nexmon – Blame – Rev 1 – /buildtools/isl-0.10/isl_transitive

1

office

1

/*

2

3

*

4

* Use of this software is governed by the GNU LGPLv2.1 license

5

*

6

* Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,

7

* Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,

8

* 91893 Orsay, France

9

*/

10

11

#include <isl_ctx_private.h>

12

#include <isl_map_private.h>

13

#include <isl/map.h>

14

#include <isl/seq.h>

15

#include <isl_space_private.h>

16

#include <isl/lp.h>

17

#include <isl/union_map.h>

18

#include <isl_mat_private.h>

19

#include <isl_options_private.h>

20

21

int isl_map_is_transitively_closed(__isl_keep isl_map *map)

{

isl_map *map2;

int closed;

map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map));

27

closed = isl_map_is_subset(map2, map);

28

isl_map_free(map2);

return closed;

}

int isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap)

34

{

35

isl_union_map *umap2;

36

int closed;

37

38

umap2 = isl_union_map_apply_range(isl_union_map_copy(umap),

39

isl_union_map_copy(umap));

40

closed = isl_union_map_is_subset(umap2, umap);

41

isl_union_map_free(umap2);

return closed;

}

/* Given a map that represents a path with the length of the path

47

* encoded as the difference between the last output coordindate

48

* and the last input coordinate, set this length to either

49

* exactly "length" (if "exactly" is set) or at least "length"

50

* (if "exactly" is not set).

51

*/

52

static __isl_give isl_map *set_path_length(__isl_take isl_map *map,

53

int exactly, int length)

54

{

55

isl_space *dim;

56

struct isl_basic_map *bmap;

unsigned d;

unsigned nparam;

int k;

isl_int *c;

if (!map)

return NULL;

dim = isl_map_get_space(map);

66

d = isl_space_dim(dim, isl_dim_in);

67

nparam = isl_space_dim(dim, isl_dim_param);

68

bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);

69

if (exactly) {

70

k = isl_basic_map_alloc_equality(bmap);

71

c = bmap->eq[k];

72

} else {

73

k = isl_basic_map_alloc_inequality(bmap);

c = bmap->ineq[k];

}

if (k < 0)

goto error;

isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));

79

isl_int_set_si(c[0], -length);

80

isl_int_set_si(c[1 + nparam + d - 1], -1);

81

isl_int_set_si(c[1 + nparam + d + d - 1], 1);

82

83

bmap = isl_basic_map_finalize(bmap);

84

map = isl_map_intersect(map, isl_map_from_basic_map(bmap));

return map;

error:

isl_basic_map_free(bmap);

isl_map_free(map);

return NULL;

}

/* Check whether the overapproximation of the power of "map" is exactly

94

* the power of "map". Let R be "map" and A_k the overapproximation.

95

* The approximation is exact if

96

*

97

* A_1 = R

98

* A_k = A_{k-1} \circ R k >= 2

99

*

100

* Since A_k is known to be an overapproximation, we only need to check

101

*

102

* A_1 \subset R

103

* A_k \subset A_{k-1} \circ R k >= 2

104

*

105

* In practice, "app" has an extra input and output coordinate

106

* to encode the length of the path. So, we first need to add

107

* this coordinate to "map" and set the length of the path to

108

* one.

109

*/

110

static int check_power_exactness(__isl_take isl_map *map,

111

__isl_take isl_map *app)

{

int exact;

isl_map *app_1;

isl_map *app_2;

map = isl_map_add_dims(map, isl_dim_in, 1);

118

map = isl_map_add_dims(map, isl_dim_out, 1);

119

map = set_path_length(map, 1, 1);

120

121

app_1 = set_path_length(isl_map_copy(app), 1, 1);

122

123

exact = isl_map_is_subset(app_1, map);

124

isl_map_free(app_1);

125

126

if (!exact || exact < 0) {

isl_map_free(app);

isl_map_free(map);

return exact;

}

app_1 = set_path_length(isl_map_copy(app), 0, 1);

133

app_2 = set_path_length(app, 0, 2);

134

app_1 = isl_map_apply_range(map, app_1);

135

136

exact = isl_map_is_subset(app_2, app_1);

137

138

isl_map_free(app_1);

139

isl_map_free(app_2);

return exact;

}

/* Check whether the overapproximation of the power of "map" is exactly

145

* the power of "map", possibly after projecting out the power (if "project"

146

* is set).

147

*

148

* If "project" is set and if "steps" can only result in acyclic paths,

149

* then we check

150

*

151

* A = R \cup (A \circ R)

152

*

153

* where A is the overapproximation with the power projected out, i.e.,

154

* an overapproximation of the transitive closure.

155

* More specifically, since A is known to be an overapproximation, we check

156

*

157

* A \subset R \cup (A \circ R)

158

*

159

* Otherwise, we check if the power is exact.

160

*

161

* Note that "app" has an extra input and output coordinate to encode

162

* the length of the part. If we are only interested in the transitive

163

* closure, then we can simply project out these coordinates first.

164

*/

165

static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,

int project)

{

isl_map *test;

int exact;

unsigned d;

if (!project)

return check_power_exactness(map, app);

174

175

d = isl_map_dim(map, isl_dim_in);

176

app = set_path_length(app, 0, 1);

177

app = isl_map_project_out(app, isl_dim_in, d, 1);

178

app = isl_map_project_out(app, isl_dim_out, d, 1);

179

180

app = isl_map_reset_space(app, isl_map_get_space(map));

181

182

test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));

183

test = isl_map_union(test, isl_map_copy(map));

184

185

exact = isl_map_is_subset(app, test);

186

187

isl_map_free(app);

188

isl_map_free(test);

isl_map_free(map);

return exact;

}

/*

* The transitive closure implementation is based on the paper

197

* "Computing the Transitive Closure of a Union of Affine Integer

198

* Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and

* Albert Cohen.

*/

/* Given a set of n offsets v_i (the rows of "steps"), construct a relation

203

* of the given dimension specification (Z^{n+1} -> Z^{n+1})

204

* that maps an element x to any element that can be reached

205

* by taking a non-negative number of steps along any of

206

* the extended offsets v'_i = [v_i 1].

207

* That is, construct

208

*

209

* { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }

210

*

211

* For any element in this relation, the number of steps taken

212

* is equal to the difference in the final coordinates.

213

*/

214

static __isl_give isl_map *path_along_steps(__isl_take isl_space *dim,

215

__isl_keep isl_mat *steps)

216

{

217

int i, j, k;

218

struct isl_basic_map *path = NULL;

unsigned d;

unsigned n;

unsigned nparam;

if (!dim || !steps)

224

goto error;

225

226

d = isl_space_dim(dim, isl_dim_in);

227

n = steps->n_row;

228

nparam = isl_space_dim(dim, isl_dim_param);

229

230

path = isl_basic_map_alloc_space(isl_space_copy(dim), n, d, n);

231

232

for (i = 0; i < n; ++i) {

233

k = isl_basic_map_alloc_div(path);

234

if (k < 0)

235

goto error;

236

isl_assert(steps->ctx, i == k, goto error);

237

isl_int_set_si(path->div[k][0], 0);

238

}

239

240

for (i = 0; i < d; ++i) {

241

k = isl_basic_map_alloc_equality(path);

242

if (k < 0)

243

goto error;

244

isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));

245

isl_int_set_si(path->eq[k][1 + nparam + i], 1);

246

isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);

247

if (i == d - 1)

248

for (j = 0; j < n; ++j)

249

isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);

250

else

251

for (j = 0; j < n; ++j)

252

isl_int_set(path->eq[k][1 + nparam + 2 * d + j],

steps->row[j][i]);

}

for (i = 0; i < n; ++i) {

257

k = isl_basic_map_alloc_inequality(path);

258

if (k < 0)

259

goto error;

260

isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));

261

isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);

262

}

263

264

isl_space_free(dim);

265

266

path = isl_basic_map_simplify(path);

267

path = isl_basic_map_finalize(path);

268

return isl_map_from_basic_map(path);

269

error:

270

isl_space_free(dim);

271

isl_basic_map_free(path);

return NULL;

}

#define IMPURE 0

#define PURE_PARAM 1

#define PURE_VAR 2

#define MIXED 3

/* Check whether the parametric constant term of constraint c is never

281

* positive in "bset".

282

*/

283

static int parametric_constant_never_positive(__isl_keep isl_basic_set *bset,

284

isl_int *c, int *div_purity)

{

unsigned d;

unsigned n_div;

unsigned nparam;

int i;

int k;

int empty;

n_div = isl_basic_set_dim(bset, isl_dim_div);

294

d = isl_basic_set_dim(bset, isl_dim_set);

295

nparam = isl_basic_set_dim(bset, isl_dim_param);

296

297

bset = isl_basic_set_copy(bset);

298

bset = isl_basic_set_cow(bset);

299

bset = isl_basic_set_extend_constraints(bset, 0, 1);

300

k = isl_basic_set_alloc_inequality(bset);

301

if (k < 0)

302

goto error;

303

isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset));

304

isl_seq_cpy(bset->ineq[k], c, 1 + nparam);

305

for (i = 0; i < n_div; ++i) {

306

if (div_purity[i] != PURE_PARAM)

307

continue;

308

isl_int_set(bset->ineq[k][1 + nparam + d + i],

309

c[1 + nparam + d + i]);

310

}

311

isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);

312

empty = isl_basic_set_is_empty(bset);

313

isl_basic_set_free(bset);

return empty;

error:

isl_basic_set_free(bset);

return -1;

}

/* Return PURE_PARAM if only the coefficients of the parameters are non-zero.

322

* Return PURE_VAR if only the coefficients of the set variables are non-zero.

323

* Return MIXED if only the coefficients of the parameters and the set

324

* variables are non-zero and if moreover the parametric constant

325

* can never attain positive values.

326

* Return IMPURE otherwise.

327

*

328

* If div_purity is NULL then we are dealing with a non-parametric set

329

* and so the constraint is obviously PURE_VAR.

330

*/

331

static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity,

int eq)

{

unsigned d;

unsigned n_div;

unsigned nparam;

int empty;

int i;

int p = 0, v = 0;

if (!div_purity)

return PURE_VAR;

n_div = isl_basic_set_dim(bset, isl_dim_div);

345

d = isl_basic_set_dim(bset, isl_dim_set);

346

nparam = isl_basic_set_dim(bset, isl_dim_param);

347

348

for (i = 0; i < n_div; ++i) {

349

if (isl_int_is_zero(c[1 + nparam + d + i]))

350

continue;

351

switch (div_purity[i]) {

352

case PURE_PARAM: p = 1; break;

353

case PURE_VAR: v = 1; break;

354

default: return IMPURE;

355

}

356

}

357

if (!p && isl_seq_first_non_zero(c + 1, nparam) == -1)

358

return PURE_VAR;

359

if (!v && isl_seq_first_non_zero(c + 1 + nparam, d) == -1)

360

return PURE_PARAM;

361

362

empty = parametric_constant_never_positive(bset, c, div_purity);

363

if (eq && empty >= 0 && !empty) {

364

isl_seq_neg(c, c, 1 + nparam + d + n_div);

365

empty = parametric_constant_never_positive(bset, c, div_purity);

366

}

367

368

return empty < 0 ? -1 : empty ? MIXED : IMPURE;

369

}

370

371

/* Return an array of integers indicating the type of each div in bset.

372

* If the div is (recursively) defined in terms of only the parameters,

373

* then the type is PURE_PARAM.

374

* If the div is (recursively) defined in terms of only the set variables,

375

* then the type is PURE_VAR.

376

* Otherwise, the type is IMPURE.

377

*/

378

static __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset)

{

int i, j;

int *div_purity;

unsigned d;

unsigned n_div;

unsigned nparam;

if (!bset)

return NULL;

n_div = isl_basic_set_dim(bset, isl_dim_div);

390

d = isl_basic_set_dim(bset, isl_dim_set);

391

nparam = isl_basic_set_dim(bset, isl_dim_param);

392

393

div_purity = isl_alloc_array(bset->ctx, int, n_div);

if (!div_purity)

return NULL;

for (i = 0; i < bset->n_div; ++i) {

398

int p = 0, v = 0;

399

if (isl_int_is_zero(bset->div[i][0])) {

400

div_purity[i] = IMPURE;

401

continue;

402

}

403

if (isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -1)

404

p = 1;

405

if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1)

406

v = 1;

407

for (j = 0; j < i; ++j) {

408

if (isl_int_is_zero(bset->div[i][2 + nparam + d + j]))

409

continue;

410

switch (div_purity[j]) {

411

case PURE_PARAM: p = 1; break;

412

case PURE_VAR: v = 1; break;

413

default: p = v = 1; break;

414

}

415

}

416

div_purity[i] = v ? p ? IMPURE : PURE_VAR : PURE_PARAM;

}

return div_purity;

}

/* Given a path with the as yet unconstrained length at position "pos",

423

* check if setting the length to zero results in only the identity

424

* mapping.

425

*/

426

static int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)

427

{

428

isl_basic_map *test = NULL;

429

isl_basic_map *id = NULL;

int k;

int is_id;

test = isl_basic_map_copy(path);

434

test = isl_basic_map_extend_constraints(test, 1, 0);

435

k = isl_basic_map_alloc_equality(test);

436

if (k < 0)

437

goto error;

438

isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test));

439

isl_int_set_si(test->eq[k][pos], 1);

440

id = isl_basic_map_identity(isl_basic_map_get_space(path));

441

is_id = isl_basic_map_is_equal(test, id);

442

isl_basic_map_free(test);

443

isl_basic_map_free(id);

444

return is_id;

445

error:

446

isl_basic_map_free(test);

return -1;

}

/* If any of the constraints is found to be impure then this function

451

* sets *impurity to 1.

452

*/

453

static __isl_give isl_basic_map *add_delta_constraints(

454

__isl_take isl_basic_map *path,

455

__isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,

456

unsigned d, int *div_purity, int eq, int *impurity)

457

{

458

int i, k;

459

int n = eq ? delta->n_eq : delta->n_ineq;

460

isl_int **delta_c = eq ? delta->eq : delta->ineq;

461

unsigned n_div;

462

463

n_div = isl_basic_set_dim(delta, isl_dim_div);

464

465

for (i = 0; i < n; ++i) {

466

isl_int *path_c;

467

int p = purity(delta, delta_c[i], div_purity, eq);

468

if (p < 0)

469

goto error;

470

if (p != PURE_VAR && p != PURE_PARAM && !*impurity)

*impurity = 1;

if (p == IMPURE)

continue;

if (eq && p != MIXED) {

475

k = isl_basic_map_alloc_equality(path);

476

path_c = path->eq[k];

477

} else {

478

k = isl_basic_map_alloc_inequality(path);

479

path_c = path->ineq[k];

}

if (k < 0)

goto error;

isl_seq_clr(path_c, 1 + isl_basic_map_total_dim(path));

484

if (p == PURE_VAR) {

485

isl_seq_cpy(path_c + off,

486

delta_c[i] + 1 + nparam, d);

487

isl_int_set(path_c[off + d], delta_c[i][0]);

488

} else if (p == PURE_PARAM) {

489

isl_seq_cpy(path_c, delta_c[i], 1 + nparam);

490

} else {

491

isl_seq_cpy(path_c + off,

492

delta_c[i] + 1 + nparam, d);

493

isl_seq_cpy(path_c, delta_c[i], 1 + nparam);

494

}

495

isl_seq_cpy(path_c + off - n_div,

496

delta_c[i] + 1 + nparam + d, n_div);

}

return path;

error:

isl_basic_map_free(path);

return NULL;

}

/* Given a set of offsets "delta", construct a relation of the

506

* given dimension specification (Z^{n+1} -> Z^{n+1}) that

507

* is an overapproximation of the relations that

508

* maps an element x to any element that can be reached

509

* by taking a non-negative number of steps along any of

510

* the elements in "delta".

511

* That is, construct an approximation of

512

*

513

* { [x] -> [y] : exists f \in \delta, k \in Z :

514

* y = x + k [f, 1] and k >= 0 }

515

*

516

* For any element in this relation, the number of steps taken

517

* is equal to the difference in the final coordinates.

518

*

519

* In particular, let delta be defined as

520

*

521

* \delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and

522

* C x + C'p + c >= 0 and

523

* D x + D'p + d >= 0 }

524

*

525

* where the constraints C x + C'p + c >= 0 are such that the parametric

526

* constant term of each constraint j, "C_j x + C'_j p + c_j",

527

* can never attain positive values, then the relation is constructed as

528

*

529

* { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and

530

* A f + k a >= 0 and B p + b >= 0 and

531

* C f + C'p + c >= 0 and k >= 1 }

532

* union { [x] -> [x] }

533

*

534

* If the zero-length paths happen to correspond exactly to the identity

535

* mapping, then we return

536

*

537

* { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and

538

* A f + k a >= 0 and B p + b >= 0 and

539

* C f + C'p + c >= 0 and k >= 0 }

*

* instead.

*

* Existentially quantified variables in \delta are handled by

544

* classifying them as independent of the parameters, purely

545

* parameter dependent and others. Constraints containing

546

* any of the other existentially quantified variables are removed.

547

* This is safe, but leads to an additional overapproximation.

548

*

549

* If there are any impure constraints, then we also eliminate

550

* the parameters from \delta, resulting in a set

551

*

552

* \delta' = { [x] : E x + e >= 0 }

553

*

554

* and add the constraints

*

* E f + k e >= 0

*

* to the constructed relation.

559

*/

560

static __isl_give isl_map *path_along_delta(__isl_take isl_space *dim,

561

__isl_take isl_basic_set *delta)

562

{

563

isl_basic_map *path = NULL;

unsigned d;

unsigned n_div;

unsigned nparam;

unsigned off;

int i, k;

int is_id;

int *div_purity = NULL;

int impurity = 0;

if (!delta)

goto error;

n_div = isl_basic_set_dim(delta, isl_dim_div);

576

d = isl_basic_set_dim(delta, isl_dim_set);

577

nparam = isl_basic_set_dim(delta, isl_dim_param);

578

path = isl_basic_map_alloc_space(isl_space_copy(dim), n_div + d + 1,

579

d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1);

580

off = 1 + nparam + 2 * (d + 1) + n_div;

581

582

for (i = 0; i < n_div + d + 1; ++i) {

583

k = isl_basic_map_alloc_div(path);

584

if (k < 0)

585

goto error;

586

isl_int_set_si(path->div[k][0], 0);

587

}

588

589

for (i = 0; i < d + 1; ++i) {

590

k = isl_basic_map_alloc_equality(path);

591

if (k < 0)

592

goto error;

593

isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));

594

isl_int_set_si(path->eq[k][1 + nparam + i], 1);

595

isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);

596

isl_int_set_si(path->eq[k][off + i], 1);

597

}

598

599

div_purity = get_div_purity(delta);

if (!div_purity)

goto error;

path = add_delta_constraints(path, delta, off, nparam, d,

604

div_purity, 1, &impurity);

605

path = add_delta_constraints(path, delta, off, nparam, d,

606

div_purity, 0, &impurity);

607

if (impurity) {

608

isl_space *dim = isl_basic_set_get_space(delta);

609

delta = isl_basic_set_project_out(delta,

610

isl_dim_param, 0, nparam);

611

delta = isl_basic_set_add(delta, isl_dim_param, nparam);

612

delta = isl_basic_set_reset_space(delta, dim);

613

if (!delta)

614

goto error;

615

path = isl_basic_map_extend_constraints(path, delta->n_eq,

616

delta->n_ineq + 1);

617

path = add_delta_constraints(path, delta, off, nparam, d,

618

NULL, 1, &impurity);

619

path = add_delta_constraints(path, delta, off, nparam, d,

620

NULL, 0, &impurity);

621

path = isl_basic_map_gauss(path, NULL);

622

}

623

624

is_id = empty_path_is_identity(path, off + d);

if (is_id < 0)

goto error;

k = isl_basic_map_alloc_inequality(path);

629

if (k < 0)

630

goto error;

631

isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));

632

if (!is_id)

633

isl_int_set_si(path->ineq[k][0], -1);

634

isl_int_set_si(path->ineq[k][off + d], 1);

635

636

free(div_purity);

637

isl_basic_set_free(delta);

638

path = isl_basic_map_finalize(path);

639

if (is_id) {

640

isl_space_free(dim);

641

return isl_map_from_basic_map(path);

642

}

643

return isl_basic_map_union(path, isl_basic_map_identity(dim));

644

error:

645

free(div_purity);

646

isl_space_free(dim);

647

isl_basic_set_free(delta);

648

isl_basic_map_free(path);

return NULL;

}

/* Given a dimension specification Z^{n+1} -> Z^{n+1} and a parameter "param",

653

* construct a map that equates the parameter to the difference

654

* in the final coordinates and imposes that this difference is positive.

655

* That is, construct

656

*

657

* { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }

658

*/

659

static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_space *dim,

660

unsigned param)

661

{

662

struct isl_basic_map *bmap;

unsigned d;

unsigned nparam;

int k;

d = isl_space_dim(dim, isl_dim_in);

668

nparam = isl_space_dim(dim, isl_dim_param);

669

bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);

670

k = isl_basic_map_alloc_equality(bmap);

671

if (k < 0)

672

goto error;

673

isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));

674

isl_int_set_si(bmap->eq[k][1 + param], -1);

675

isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);

676

isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);

677

678

k = isl_basic_map_alloc_inequality(bmap);

679

if (k < 0)

680

goto error;

681

isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));

682

isl_int_set_si(bmap->ineq[k][1 + param], 1);

683

isl_int_set_si(bmap->ineq[k][0], -1);

684

685

bmap = isl_basic_map_finalize(bmap);

686

return isl_map_from_basic_map(bmap);

687

error:

688

isl_basic_map_free(bmap);

return NULL;

}

/* Check whether "path" is acyclic, where the last coordinates of domain

693

* and range of path encode the number of steps taken.

694

* That is, check whether

695

*

696

* { d | d = y - x and (x,y) in path }

697

*

698

* does not contain any element with positive last coordinate (positive length)

699

* and zero remaining coordinates (cycle).

700

*/

701

static int is_acyclic(__isl_take isl_map *path)

{

int i;

int acyclic;

unsigned dim;

struct isl_set *delta;

707

708

delta = isl_map_deltas(path);

709

dim = isl_set_dim(delta, isl_dim_set);

710

for (i = 0; i < dim; ++i) {

711

if (i == dim -1)

712

delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);

713

else

714

delta = isl_set_fix_si(delta, isl_dim_set, i, 0);

715

}

716

717

acyclic = isl_set_is_empty(delta);

718

isl_set_free(delta);

return acyclic;

}

/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D

724

* and a dimension specification (Z^{n+1} -> Z^{n+1}),

725

* construct a map that is an overapproximation of the map

726

* that takes an element from the space D \times Z to another

727

* element from the same space, such that the first n coordinates of the

728

* difference between them is a sum of differences between images

729

* and pre-images in one of the R_i and such that the last coordinate

730

* is equal to the number of steps taken.

731

* That is, let

732

*

733

* \Delta_i = { y - x | (x, y) in R_i }

734

*

735

* then the constructed map is an overapproximation of

736

*

737

* { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :

738

* d = (\sum_i k_i \delta_i, \sum_i k_i) }

739

*

740

* The elements of the singleton \Delta_i's are collected as the

741

* rows of the steps matrix. For all these \Delta_i's together,

742

* a single path is constructed.

743

* For each of the other \Delta_i's, we compute an overapproximation

744

* of the paths along elements of \Delta_i.

745

* Since each of these paths performs an addition, composition is

746

* symmetric and we can simply compose all resulting paths in any order.

747

*/

748

static __isl_give isl_map *construct_extended_path(__isl_take isl_space *dim,

749

__isl_keep isl_map *map, int *project)

750

{

751

struct isl_mat *steps = NULL;

752

struct isl_map *path = NULL;

unsigned d;

int i, j, n;

d = isl_map_dim(map, isl_dim_in);

757

758

path = isl_map_identity(isl_space_copy(dim));

759

760

steps = isl_mat_alloc(map->ctx, map->n, d);

if (!steps)

goto error;

n = 0;

for (i = 0; i < map->n; ++i) {

766

struct isl_basic_set *delta;

767

768

delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));

769

770

for (j = 0; j < d; ++j) {

771

int fixed;

772

773

fixed = isl_basic_set_plain_dim_is_fixed(delta, j,

774

&steps->row[n][j]);

775

if (fixed < 0) {

776

isl_basic_set_free(delta);

goto error;

}

if (!fixed)

break;

}

if (j < d) {

path = isl_map_apply_range(path,

786

path_along_delta(isl_space_copy(dim), delta));

787

path = isl_map_coalesce(path);

788

} else {

789

isl_basic_set_free(delta);

++n;

}

}

if (n > 0) {

steps->n_row = n;

path = isl_map_apply_range(path,

797

path_along_steps(isl_space_copy(dim), steps));

798

}

799

800

if (project && *project) {

801

*project = is_acyclic(isl_map_copy(path));

if (*project < 0)

goto error;

}

isl_space_free(dim);

807

isl_mat_free(steps);

808

return path;

809

error:

810

isl_space_free(dim);

811

isl_mat_free(steps);

812

isl_map_free(path);

return NULL;

}

static int isl_set_overlaps(__isl_keep isl_set *set1, __isl_keep isl_set *set2)

{

isl_set *i;

int no_overlap;

if (!isl_space_tuple_match(set1->dim, isl_dim_set, set2->dim, isl_dim_set))

822

return 0;

823

824

i = isl_set_intersect(isl_set_copy(set1), isl_set_copy(set2));

825

no_overlap = isl_set_is_empty(i);

826

isl_set_free(i);

827

828

return no_overlap < 0 ? -1 : !no_overlap;

829

}

830

831

/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D

832

* and a dimension specification (Z^{n+1} -> Z^{n+1}),

833

* construct a map that is an overapproximation of the map

834

* that takes an element from the dom R \times Z to an

835

* element from ran R \times Z, such that the first n coordinates of the

836

* difference between them is a sum of differences between images

837

* and pre-images in one of the R_i and such that the last coordinate

838

* is equal to the number of steps taken.

839

* That is, let

840

*

841

* \Delta_i = { y - x | (x, y) in R_i }

842

*

843

* then the constructed map is an overapproximation of

844

*

845

* { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :

846

* d = (\sum_i k_i \delta_i, \sum_i k_i) and

847

* x in dom R and x + d in ran R and

848

* \sum_i k_i >= 1 }

849

*/

850

static __isl_give isl_map *construct_component(__isl_take isl_space *dim,

851

__isl_keep isl_map *map, int *exact, int project)

852

{

853

struct isl_set *domain = NULL;

854

struct isl_set *range = NULL;

855

struct isl_map *app = NULL;

856

struct isl_map *path = NULL;

857

858

domain = isl_map_domain(isl_map_copy(map));

859

domain = isl_set_coalesce(domain);

860

range = isl_map_range(isl_map_copy(map));

861

range = isl_set_coalesce(range);

862

if (!isl_set_overlaps(domain, range)) {

863

isl_set_free(domain);

864

isl_set_free(range);

865

isl_space_free(dim);

866

867

map = isl_map_copy(map);

868

map = isl_map_add_dims(map, isl_dim_in, 1);

869

map = isl_map_add_dims(map, isl_dim_out, 1);

870

map = set_path_length(map, 1, 1);

871

return map;

872

}

873

app = isl_map_from_domain_and_range(domain, range);

874

app = isl_map_add_dims(app, isl_dim_in, 1);

875

app = isl_map_add_dims(app, isl_dim_out, 1);

876

877

path = construct_extended_path(isl_space_copy(dim), map,

878

exact && *exact ? &project : NULL);

879

app = isl_map_intersect(app, path);

880

881

if (exact && *exact &&

882

(*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),

project)) < 0)

goto error;

isl_space_free(dim);

887

app = set_path_length(app, 0, 1);

888

return app;

889

error:

890

isl_space_free(dim);

isl_map_free(app);

return NULL;

}

/* Call construct_component and, if "project" is set, project out

896

* the final coordinates.

897

*/

898

static __isl_give isl_map *construct_projected_component(

899

__isl_take isl_space *dim,

900

__isl_keep isl_map *map, int *exact, int project)

{

isl_map *app;

unsigned d;

if (!dim)

return NULL;

d = isl_space_dim(dim, isl_dim_in);

908

909

app = construct_component(dim, map, exact, project);

910

if (project) {

911

app = isl_map_project_out(app, isl_dim_in, d - 1, 1);

912

app = isl_map_project_out(app, isl_dim_out, d - 1, 1);

}

return app;

}

/* Compute an extended version, i.e., with path lengths, of

918

* an overapproximation of the transitive closure of "bmap"

919

* with path lengths greater than or equal to zero and with

920

* domain and range equal to "dom".

921

*/

922

static __isl_give isl_map *q_closure(__isl_take isl_space *dim,

923

__isl_take isl_set *dom, __isl_keep isl_basic_map *bmap, int *exact)

{

int project = 1;

isl_map *path;

isl_map *map;

isl_map *app;

dom = isl_set_add_dims(dom, isl_dim_set, 1);

931

app = isl_map_from_domain_and_range(dom, isl_set_copy(dom));

932

map = isl_map_from_basic_map(isl_basic_map_copy(bmap));

933

path = construct_extended_path(dim, map, &project);

934

app = isl_map_intersect(app, path);

935

936

if ((*exact = check_exactness(map, isl_map_copy(app), project)) < 0)

goto error;

return app;

error:

isl_map_free(app);

return NULL;

}

/* Check whether qc has any elements of length at least one

946

* with domain and/or range outside of dom and ran.

947

*/

948

static int has_spurious_elements(__isl_keep isl_map *qc,

949

__isl_keep isl_set *dom, __isl_keep isl_set *ran)

{

isl_set *s;

int subset;

unsigned d;

if (!qc || !dom || !ran)

956

return -1;

957

958

d = isl_map_dim(qc, isl_dim_in);

959

960

qc = isl_map_copy(qc);

961

qc = set_path_length(qc, 0, 1);

962

qc = isl_map_project_out(qc, isl_dim_in, d - 1, 1);

963

qc = isl_map_project_out(qc, isl_dim_out, d - 1, 1);

964

965

s = isl_map_domain(isl_map_copy(qc));

966

subset = isl_set_is_subset(s, dom);

isl_set_free(s);

if (subset < 0)

goto error;

if (!subset) {

isl_map_free(qc);

return 1;

}

s = isl_map_range(qc);

976

subset = isl_set_is_subset(s, ran);

977

isl_set_free(s);

978

979

return subset < 0 ? -1 : !subset;

error:

isl_map_free(qc);

return -1;

}

#define LEFT 2

#define RIGHT 1

/* For each basic map in "map", except i, check whether it combines

989

* with the transitive closure that is reflexive on C combines

990

* to the left and to the right.

991

*

992

* In particular, if

993

*

994

* dom map_j \subseteq C

995

*

996

* then right[j] is set to 1. Otherwise, if

997

*

998

* ran map_i \cap dom map_j = \emptyset

999

*

1000

* then right[j] is set to 0. Otherwise, composing to the right

1001

* is impossible.

1002

*

1003

* Similar, for composing to the left, we have if

1004

*

1005

* ran map_j \subseteq C

1006

*

1007

* then left[j] is set to 1. Otherwise, if

1008

*

1009

* dom map_i \cap ran map_j = \emptyset

1010

*

1011

* then left[j] is set to 0. Otherwise, composing to the left

1012

* is impossible.

1013

*

1014

* The return value is or'd with LEFT if composing to the left

1015

* is possible and with RIGHT if composing to the right is possible.

1016

*/

1017

static int composability(__isl_keep isl_set *C, int i,

1018

isl_set **dom, isl_set **ran, int *left, int *right,

1019

__isl_keep isl_map *map)

{

int j;

int ok;

ok = LEFT | RIGHT;

for (j = 0; j < map->n && ok; ++j) {

1026

int overlaps, subset;

if (j == i)

continue;

if (ok & RIGHT) {

if (!dom[j])

dom[j] = isl_set_from_basic_set(

1033

isl_basic_map_domain(

1034

isl_basic_map_copy(map->p[j])));

1035

if (!dom[j])

1036

return -1;

1037

overlaps = isl_set_overlaps(ran[i], dom[j]);

if (overlaps < 0)

return -1;

if (!overlaps)

right[j] = 0;

else {

subset = isl_set_is_subset(dom[j], C);

if (subset < 0)

return -1;

if (subset)

right[j] = 1;

else

ok &= ~RIGHT;

}

}

if (ok & LEFT) {

if (!ran[j])

ran[j] = isl_set_from_basic_set(

1056

isl_basic_map_range(

1057

isl_basic_map_copy(map->p[j])));

1058

if (!ran[j])

1059

return -1;

1060

overlaps = isl_set_overlaps(dom[i], ran[j]);

if (overlaps < 0)

return -1;

if (!overlaps)

left[j] = 0;

else {

subset = isl_set_is_subset(ran[j], C);

if (subset < 0)

return -1;

if (subset)

left[j] = 1;

else

ok &= ~LEFT;

}

}

}

return ok;

}

static __isl_give isl_map *anonymize(__isl_take isl_map *map)

1081

{

1082

map = isl_map_reset(map, isl_dim_in);

1083

map = isl_map_reset(map, isl_dim_out);

return map;

}

/* Return a map that is a union of the basic maps in "map", except i,

1088

* composed to left and right with qc based on the entries of "left"

1089

* and "right".

1090

*/

1091

static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,

1092

__isl_take isl_map *qc, int *left, int *right)

{

int j;

isl_map *comp;

comp = isl_map_empty(isl_map_get_space(map));

1098

for (j = 0; j < map->n; ++j) {

isl_map *map_j;

if (j == i)

continue;

map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j]));

1105

map_j = anonymize(map_j);

1106

if (left && left[j])

1107

map_j = isl_map_apply_range(map_j, isl_map_copy(qc));

1108

if (right && right[j])

1109

map_j = isl_map_apply_range(isl_map_copy(qc), map_j);

1110

comp = isl_map_union(comp, map_j);

1111

}

1112

1113

comp = isl_map_compute_divs(comp);

1114

comp = isl_map_coalesce(comp);

isl_map_free(qc);

return comp;

}

/* Compute the transitive closure of "map" incrementally by

1122

* computing

1123

*

1124

* map_i^+ \cup qc^+

*

* or

*

* map_i^+ \cup ((id \cup map_i^) \circ qc^+)

*

* or

*

* map_i^+ \cup (qc^+ \circ (id \cup map_i^))

1133

*

1134

* depending on whether left or right are NULL.

1135

*/

1136

static __isl_give isl_map *compute_incremental(

1137

__isl_take isl_space *dim, __isl_keep isl_map *map,

1138

int i, __isl_take isl_map *qc, int *left, int *right, int *exact)

{

isl_map *map_i;

isl_map *tc;

isl_map *rtc = NULL;

if (!map)

goto error;

isl_assert(map->ctx, left || right, goto error);

1147

1148

map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));

1149

tc = construct_projected_component(isl_space_copy(dim), map_i,

1150

exact, 1);

1151

isl_map_free(map_i);

1152

1153

if (*exact)

1154

qc = isl_map_transitive_closure(qc, exact);

1155

1156

if (!*exact) {

1157

isl_space_free(dim);

1158

isl_map_free(tc);

1159

isl_map_free(qc);

1160

return isl_map_universe(isl_map_get_space(map));

1161

}

1162

1163

if (!left || !right)

1164

rtc = isl_map_union(isl_map_copy(tc),

1165

isl_map_identity(isl_map_get_space(tc)));

1166

if (!right)

1167

qc = isl_map_apply_range(rtc, qc);

1168

if (!left)

1169

qc = isl_map_apply_range(qc, rtc);

1170

qc = isl_map_union(tc, qc);

1171

1172

isl_space_free(dim);

return qc;

error:

isl_space_free(dim);

isl_map_free(qc);

return NULL;

}

/* Given a map "map", try to find a basic map such that

1182

* map^+ can be computed as

1183

*

1184

* map^+ = map_i^+ \cup

1185

* \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+

1186

*

1187

* with C the simple hull of the domain and range of the input map.

1188

* map_i^ \cup Id_C is computed by allowing the path lengths to be zero

1189

* and by intersecting domain and range with C.

1190

* Of course, we need to check that this is actually equal to map_i^ \cup Id_C.

1191

* Also, we only use the incremental computation if all the transitive

1192

* closures are exact and if the number of basic maps in the union,

1193

* after computing the integer divisions, is smaller than the number

1194

* of basic maps in the input map.

1195

*/

1196

static int incemental_on_entire_domain(__isl_keep isl_space *dim,

1197

__isl_keep isl_map *map,

1198

isl_set **dom, isl_set **ran, int *left, int *right,

1199

__isl_give isl_map **res)

{

int i;

isl_set *C;

unsigned d;

*res = NULL;

C = isl_set_union(isl_map_domain(isl_map_copy(map)),

1208

isl_map_range(isl_map_copy(map)));

1209

C = isl_set_from_basic_set(isl_set_simple_hull(C));

if (!C)

return -1;

if (C->n != 1) {

isl_set_free(C);

return 0;

}

d = isl_map_dim(map, isl_dim_in);

1218

1219

for (i = 0; i < map->n; ++i) {

1220

isl_map *qc;

1221

int exact_i, spurious;

1222

int j;

1223

dom[i] = isl_set_from_basic_set(isl_basic_map_domain(

1224

isl_basic_map_copy(map->p[i])));

1225

ran[i] = isl_set_from_basic_set(isl_basic_map_range(

1226

isl_basic_map_copy(map->p[i])));

1227

qc = q_closure(isl_space_copy(dim), isl_set_copy(C),

1228

map->p[i], &exact_i);

if (!qc)

goto error;

if (!exact_i) {

isl_map_free(qc);

continue;

}

spurious = has_spurious_elements(qc, dom[i], ran[i]);

if (spurious) {

isl_map_free(qc);

if (spurious < 0)

goto error;

continue;

}

qc = isl_map_project_out(qc, isl_dim_in, d, 1);

1243

qc = isl_map_project_out(qc, isl_dim_out, d, 1);

1244

qc = isl_map_compute_divs(qc);

1245

for (j = 0; j < map->n; ++j)

1246

left[j] = right[j] = 1;

1247

qc = compose(map, i, qc, left, right);

1248

if (!qc)

1249

goto error;

1250

if (qc->n >= map->n) {

isl_map_free(qc);

continue;

}

*res = compute_incremental(isl_space_copy(dim), map, i, qc,

1255

left, right, &exact_i);

if (!*res)

goto error;

if (exact_i)

break;

isl_map_free(*res);

*res = NULL;

}

isl_set_free(C);

return *res != NULL;

error:

isl_set_free(C);

return -1;

}

/* Try and compute the transitive closure of "map" as

1273

*

1274

* map^+ = map_i^+ \cup

1275

* \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+

1276

*

1277

* with C either the simple hull of the domain and range of the entire

1278

* map or the simple hull of domain and range of map_i.

1279

*/

1280

static __isl_give isl_map *incremental_closure(__isl_take isl_space *dim,

1281

__isl_keep isl_map *map, int *exact, int project)

1282

{

1283

int i;

1284

isl_set **dom = NULL;

1285

isl_set **ran = NULL;

int *left = NULL;

int *right = NULL;

isl_set *C;

unsigned d;

isl_map *res = NULL;

1291

1292

if (!project)

1293

return construct_projected_component(dim, map, exact, project);

if (!map)

goto error;

if (map->n <= 1)

return construct_projected_component(dim, map, exact, project);

1299

1300

d = isl_map_dim(map, isl_dim_in);

1301

1302

dom = isl_calloc_array(map->ctx, isl_set *, map->n);

1303

ran = isl_calloc_array(map->ctx, isl_set *, map->n);

1304

left = isl_calloc_array(map->ctx, int, map->n);

1305

right = isl_calloc_array(map->ctx, int, map->n);

1306

if (!ran || !dom || !left || !right)

1307

goto error;

1308

1309

if (incemental_on_entire_domain(dim, map, dom, ran, left, right, &res) < 0)

1310

goto error;

1311

1312

for (i = 0; !res && i < map->n; ++i) {

1313

isl_map *qc;

1314

int exact_i, spurious, comp;

1315

if (!dom[i])

1316

dom[i] = isl_set_from_basic_set(

1317

isl_basic_map_domain(

1318

isl_basic_map_copy(map->p[i])));

if (!dom[i])

goto error;

if (!ran[i])

ran[i] = isl_set_from_basic_set(

1323

isl_basic_map_range(

1324

isl_basic_map_copy(map->p[i])));

1325

if (!ran[i])

1326

goto error;

1327

C = isl_set_union(isl_set_copy(dom[i]),

1328

isl_set_copy(ran[i]));

1329

C = isl_set_from_basic_set(isl_set_simple_hull(C));

if (!C)

goto error;

if (C->n != 1) {

isl_set_free(C);

continue;

}

comp = composability(C, i, dom, ran, left, right, map);

1337

if (!comp || comp < 0) {

isl_set_free(C);

if (comp < 0)

goto error;

continue;

}

qc = q_closure(isl_space_copy(dim), C, map->p[i], &exact_i);

if (!qc)

goto error;

if (!exact_i) {

isl_map_free(qc);

continue;

}

spurious = has_spurious_elements(qc, dom[i], ran[i]);

if (spurious) {

isl_map_free(qc);

if (spurious < 0)

goto error;

continue;

}

qc = isl_map_project_out(qc, isl_dim_in, d, 1);

1358

qc = isl_map_project_out(qc, isl_dim_out, d, 1);

1359

qc = isl_map_compute_divs(qc);

1360

qc = compose(map, i, qc, (comp & LEFT) ? left : NULL,

1361

(comp & RIGHT) ? right : NULL);

1362

if (!qc)

1363

goto error;

1364

if (qc->n >= map->n) {

isl_map_free(qc);

continue;

}

res = compute_incremental(isl_space_copy(dim), map, i, qc,

1369

(comp & LEFT) ? left : NULL,

1370

(comp & RIGHT) ? right : NULL, &exact_i);

if (!res)

goto error;

if (exact_i)

break;

isl_map_free(res);

res = NULL;

}

for (i = 0; i < map->n; ++i) {

1380

isl_set_free(dom[i]);

1381

isl_set_free(ran[i]);

}

free(dom);

free(ran);

free(left);

free(right);

if (res) {

isl_space_free(dim);

return res;

}

return construct_projected_component(dim, map, exact, project);

1394

error:

1395

if (dom)

1396

for (i = 0; i < map->n; ++i)

1397

isl_set_free(dom[i]);

1398

free(dom);

1399

if (ran)

1400

for (i = 0; i < map->n; ++i)

1401

isl_set_free(ran[i]);

free(ran);

free(left);

free(right);

isl_space_free(dim);

return NULL;

}

/* Given an array of sets "set", add "dom" at position "pos"

1410

* and search for elements at earlier positions that overlap with "dom".

1411

* If any can be found, then merge all of them, together with "dom", into

1412

* a single set and assign the union to the first in the array,

1413

* which becomes the new group leader for all groups involved in the merge.

1414

* During the search, we only consider group leaders, i.e., those with

1415

* group[i] = i, as the other sets have already been combined

1416

* with one of the group leaders.

1417

*/

1418

static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos)

{

int i;

group[pos] = pos;

set[pos] = isl_set_copy(dom);

1424

1425

for (i = pos - 1; i >= 0; --i) {

int o;

if (group[i] != i)

continue;

o = isl_set_overlaps(set[i], dom);

if (o < 0)

goto error;

if (!o)

continue;

set[i] = isl_set_union(set[i], set[group[pos]]);

1438

set[group[pos]] = NULL;

1439

if (!set[i])

1440

goto error;

1441

group[group[pos]] = i;

group[pos] = i;

}

isl_set_free(dom);

return 0;

error:

isl_set_free(dom);

return -1;

}

/* Replace each entry in the n by n grid of maps by the cross product

1453

* with the relation { [i] -> [i + 1] }.

1454

*/

1455

static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)

{

int i, j, k;

isl_space *dim;

isl_basic_map *bstep;

isl_map *step;

unsigned nparam;

if (!map)

return -1;

dim = isl_map_get_space(map);

1467

nparam = isl_space_dim(dim, isl_dim_param);

1468

dim = isl_space_drop_dims(dim, isl_dim_in, 0, isl_space_dim(dim, isl_dim_in));

1469

dim = isl_space_drop_dims(dim, isl_dim_out, 0, isl_space_dim(dim, isl_dim_out));

1470

dim = isl_space_add_dims(dim, isl_dim_in, 1);

1471

dim = isl_space_add_dims(dim, isl_dim_out, 1);

1472

bstep = isl_basic_map_alloc_space(dim, 0, 1, 0);

1473

k = isl_basic_map_alloc_equality(bstep);

1474

if (k < 0) {

1475

isl_basic_map_free(bstep);

1476

return -1;

1477

}

1478

isl_seq_clr(bstep->eq[k], 1 + isl_basic_map_total_dim(bstep));

1479

isl_int_set_si(bstep->eq[k][0], 1);

1480

isl_int_set_si(bstep->eq[k][1 + nparam], 1);

1481

isl_int_set_si(bstep->eq[k][1 + nparam + 1], -1);

1482

bstep = isl_basic_map_finalize(bstep);

1483

step = isl_map_from_basic_map(bstep);

1484

1485

for (i = 0; i < n; ++i)

1486

for (j = 0; j < n; ++j)

1487

grid[i][j] = isl_map_product(grid[i][j],

1488

isl_map_copy(step));

1489

1490

isl_map_free(step);

return 0;

}

/* The core of the Floyd-Warshall algorithm.

1496

* Updates the given n x x matrix of relations in place.

1497

*

1498

* The algorithm iterates over all vertices. In each step, the whole

1499

* matrix is updated to include all paths that go to the current vertex,

1500

* possibly stay there a while (including passing through earlier vertices)

1501

* and then come back. At the start of each iteration, the diagonal

1502

* element corresponding to the current vertex is replaced by its

1503

* transitive closure to account for all indirect paths that stay

1504

* in the current vertex.

1505

*/

1506

static void floyd_warshall_iterate(isl_map ***grid, int n, int *exact)

{

int r, p, q;

for (r = 0; r < n; ++r) {

1511

int r_exact;

1512

grid[r][r] = isl_map_transitive_closure(grid[r][r],

1513

(exact && *exact) ? &r_exact : NULL);

1514

if (exact && *exact && !r_exact)

1515

*exact = 0;

1516

1517

for (p = 0; p < n; ++p)

1518

for (q = 0; q < n; ++q) {

1519

isl_map *loop;

1520

if (p == r && q == r)

1521

continue;

1522

loop = isl_map_apply_range(

1523

isl_map_copy(grid[p][r]),

1524

isl_map_copy(grid[r][q]));

1525

grid[p][q] = isl_map_union(grid[p][q], loop);

1526

loop = isl_map_apply_range(

1527

isl_map_copy(grid[p][r]),

1528

isl_map_apply_range(

1529

isl_map_copy(grid[r][r]),

1530

isl_map_copy(grid[r][q])));

1531

grid[p][q] = isl_map_union(grid[p][q], loop);

1532

grid[p][q] = isl_map_coalesce(grid[p][q]);

}

}

}

/* Given a partition of the domains and ranges of the basic maps in "map",

1538

* apply the Floyd-Warshall algorithm with the elements in the partition

1539

* as vertices.

1540

*

1541

* In particular, there are "n" elements in the partition and "group" is

1542

* an array of length 2 * map->n with entries in [0,n-1].

1543

*

1544

* We first construct a matrix of relations based on the partition information,

1545

* apply Floyd-Warshall on this matrix of relations and then take the

1546

* union of all entries in the matrix as the final result.

1547

*

1548

* If we are actually computing the power instead of the transitive closure,

1549

* i.e., when "project" is not set, then the result should have the

1550

* path lengths encoded as the difference between an extra pair of

1551

* coordinates. We therefore apply the nested transitive closures

1552

* to relations that include these lengths. In particular, we replace

1553

* the input relation by the cross product with the unit length relation

1554

* { [i] -> [i + 1] }.

1555

*/

1556

static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_space *dim,

1557

__isl_keep isl_map *map, int *exact, int project, int *group, int n)

1558

{

1559

int i, j, k;

1560

isl_map ***grid = NULL;

isl_map *app;

if (!map)

goto error;

if (n == 1) {

free(group);

return incremental_closure(dim, map, exact, project);

1569

}

1570

1571

grid = isl_calloc_array(map->ctx, isl_map **, n);

1572

if (!grid)

1573

goto error;

1574

for (i = 0; i < n; ++i) {

1575

grid[i] = isl_calloc_array(map->ctx, isl_map *, n);

1576

if (!grid[i])

1577

goto error;

1578

for (j = 0; j < n; ++j)

1579

grid[i][j] = isl_map_empty(isl_map_get_space(map));

1580

}

1581

1582

for (k = 0; k < map->n; ++k) {

1583

i = group[2 * k];

1584

j = group[2 * k + 1];

1585

grid[i][j] = isl_map_union(grid[i][j],

1586

isl_map_from_basic_map(

1587

isl_basic_map_copy(map->p[k])));

1588

}

1589

1590

if (!project && add_length(map, grid, n) < 0)

1591

goto error;

1592

1593

floyd_warshall_iterate(grid, n, exact);

1594

1595

app = isl_map_empty(isl_map_get_space(map));

1596

1597

for (i = 0; i < n; ++i) {

1598

for (j = 0; j < n; ++j)

1599

app = isl_map_union(app, grid[i][j]);

free(grid[i]);

}

free(grid);

free(group);

isl_space_free(dim);

return app;

error:

if (grid)

for (i = 0; i < n; ++i) {

1611

if (!grid[i])

1612

continue;

1613

for (j = 0; j < n; ++j)

1614

isl_map_free(grid[i][j]);

free(grid[i]);

}

free(grid);

free(group);

isl_space_free(dim);

return NULL;

}

/* Partition the domains and ranges of the n basic relations in list

1624

* into disjoint cells.

1625

*

1626

* To find the partition, we simply consider all of the domains

1627

* and ranges in turn and combine those that overlap.

1628

* "set" contains the partition elements and "group" indicates

1629

* to which partition element a given domain or range belongs.

1630

* The domain of basic map i corresponds to element 2 * i in these arrays,

1631

* while the domain corresponds to element 2 * i + 1.

1632

* During the construction group[k] is either equal to k,

1633

* in which case set[k] contains the union of all the domains and

1634

* ranges in the corresponding group, or is equal to some l < k,

1635

* with l another domain or range in the same group.

1636

*/

1637

static int *setup_groups(isl_ctx *ctx, __isl_keep isl_basic_map **list, int n,

1638

isl_set ***set, int *n_group)

{

int i;

int *group = NULL;

int g;

*set = isl_calloc_array(ctx, isl_set *, 2 * n);

1645

group = isl_alloc_array(ctx, int, 2 * n);

1646

1647

if (!*set || !group)

1648

goto error;

1649

1650

for (i = 0; i < n; ++i) {

1651

isl_set *dom;

1652

dom = isl_set_from_basic_set(isl_basic_map_domain(

1653

isl_basic_map_copy(list[i])));

1654

if (merge(*set, group, dom, 2 * i) < 0)

1655

goto error;

1656

dom = isl_set_from_basic_set(isl_basic_map_range(

1657

isl_basic_map_copy(list[i])));

1658

if (merge(*set, group, dom, 2 * i + 1) < 0)

goto error;

}

g = 0;

for (i = 0; i < 2 * n; ++i)

1664

if (group[i] == i) {

1665

if (g != i) {

1666

(*set)[g] = (*set)[i];

(*set)[i] = NULL;

}

group[i] = g++;

} else

group[i] = group[group[i]];

*n_group = g;

return group;

error:

if (*set) {

for (i = 0; i < 2 * n; ++i)

1679

isl_set_free((*set)[i]);

free(*set);

*set = NULL;

}

free(group);

return NULL;

}

/* Check if the domains and ranges of the basic maps in "map" can

1688

* be partitioned, and if so, apply Floyd-Warshall on the elements

1689

* of the partition. Note that we also apply this algorithm

1690

* if we want to compute the power, i.e., when "project" is not set.

1691

* However, the results are unlikely to be exact since the recursive

1692

* calls inside the Floyd-Warshall algorithm typically result in

1693

* non-linear path lengths quite quickly.

1694

*/

1695

static __isl_give isl_map *floyd_warshall(__isl_take isl_space *dim,

1696

__isl_keep isl_map *map, int *exact, int project)

1697

{

1698

int i;

1699

isl_set **set = NULL;

int *group = NULL;

int n;

if (!map)

goto error;

if (map->n <= 1)

return incremental_closure(dim, map, exact, project);

1707

1708

group = setup_groups(map->ctx, map->p, map->n, &set, &n);

if (!group)

goto error;

for (i = 0; i < 2 * map->n; ++i)

1713

isl_set_free(set[i]);

free(set);

return floyd_warshall_with_groups(dim, map, exact, project, group, n);

1718

error:

1719

isl_space_free(dim);

return NULL;

}

/* Structure for representing the nodes in the graph being traversed

1724

* using Tarjan's algorithm.

1725

* index represents the order in which nodes are visited.

1726

* min_index is the index of the root of a (sub)component.

1727

* on_stack indicates whether the node is currently on the stack.

1728

*/

1729

struct basic_map_sort_node {

int index;

int min_index;

int on_stack;

};

/* Structure for representing the graph being traversed

1735

* using Tarjan's algorithm.

1736

* len is the number of nodes

1737

* node is an array of nodes

1738

* stack contains the nodes on the path from the root to the current node

1739

* sp is the stack pointer

1740

* index is the index of the last node visited

1741

* order contains the elements of the components separated by -1

1742

* op represents the current position in order

1743

*

1744

* check_closed is set if we may have used the fact that

1745

* a pair of basic maps can be interchanged

1746

*/

1747

struct basic_map_sort {

1748

int len;

1749

struct basic_map_sort_node *node;

int *stack;

int sp;

int index;

int *order;

int op;

int check_closed;

};

static void basic_map_sort_free(struct basic_map_sort *s)

{

if (!s)

return;

free(s->node);

free(s->stack);

free(s->order);

free(s);

}

static struct basic_map_sort *basic_map_sort_alloc(struct isl_ctx *ctx, int len)

1769

{

1770

struct basic_map_sort *s;

1771

int i;

1772

1773

s = isl_calloc_type(ctx, struct basic_map_sort);

if (!s)

return NULL;

s->len = len;

s->node = isl_alloc_array(ctx, struct basic_map_sort_node, len);

1778

if (!s->node)

1779

goto error;

1780

for (i = 0; i < len; ++i)

1781

s->node[i].index = -1;

1782

s->stack = isl_alloc_array(ctx, int, len);

1783

if (!s->stack)

1784

goto error;

1785

s->order = isl_alloc_array(ctx, int, 2 * len);

if (!s->order)

goto error;

s->sp = 0;

s->index = 0;

s->op = 0;

s->check_closed = 0;

return s;

error:

basic_map_sort_free(s);

return NULL;

}

/* Check whether in the computation of the transitive closure

1802

* "bmap1" (R_1) should follow (or be part of the same component as)

1803

* "bmap2" (R_2).

1804

*

1805

* That is check whether

*

* R_1 \circ R_2

*

* is a subset of

*

* R_2 \circ R_1

*

* If so, then there is no reason for R_1 to immediately follow R_2

1814

* in any path.

1815

*

1816

* *check_closed is set if the subset relation holds while

1817

* R_1 \circ R_2 is not empty.

1818

*/

1819

static int basic_map_follows(__isl_keep isl_basic_map *bmap1,

1820

__isl_keep isl_basic_map *bmap2, int *check_closed)

1821

{

1822

struct isl_map *map12 = NULL;

1823

struct isl_map *map21 = NULL;

1824

int subset;

1825

1826

if (!isl_space_tuple_match(bmap1->dim, isl_dim_in, bmap2->dim, isl_dim_out))

1827

return 0;

1828

1829

map21 = isl_map_from_basic_map(

1830

isl_basic_map_apply_range(

1831

isl_basic_map_copy(bmap2),

1832

isl_basic_map_copy(bmap1)));

1833

subset = isl_map_is_empty(map21);

if (subset < 0)

goto error;

if (subset) {

isl_map_free(map21);

return 0;

}

if (!isl_space_tuple_match(bmap1->dim, isl_dim_in, bmap1->dim, isl_dim_out) ||

1842

!isl_space_tuple_match(bmap2->dim, isl_dim_in, bmap2->dim, isl_dim_out)) {

1843

isl_map_free(map21);

return 1;

}

map12 = isl_map_from_basic_map(

1848

isl_basic_map_apply_range(

1849

isl_basic_map_copy(bmap1),

1850

isl_basic_map_copy(bmap2)));

1851

1852

subset = isl_map_is_subset(map21, map12);

1853

1854

isl_map_free(map12);

1855

isl_map_free(map21);

if (subset)

*check_closed = 1;

return subset < 0 ? -1 : !subset;

1861

error:

1862

isl_map_free(map21);

return -1;

}

/* Perform Tarjan's algorithm for computing the strongly connected components

1867

* in the graph with the disjuncts of "map" as vertices and with an

1868

* edge between any pair of disjuncts such that the first has

1869

* to be applied after the second.

1870

*/

1871

static int power_components_tarjan(struct basic_map_sort *s,

1872

__isl_keep isl_basic_map **list, int i)

{

int j;

s->node[i].index = s->index;

1877

s->node[i].min_index = s->index;

1878

s->node[i].on_stack = 1;

1879

s->index++;

1880

s->stack[s->sp++] = i;

1881

1882

for (j = s->len - 1; j >= 0; --j) {

int f;

if (j == i)

continue;

if (s->node[j].index >= 0 &&

1888

(!s->node[j].on_stack ||

1889

s->node[j].index > s->node[i].min_index))

1890

continue;

1891

1892

f = basic_map_follows(list[i], list[j], &s->check_closed);

if (f < 0)

return -1;

if (!f)

continue;

if (s->node[j].index < 0) {

1899

power_components_tarjan(s, list, j);

1900

if (s->node[j].min_index < s->node[i].min_index)

1901

s->node[i].min_index = s->node[j].min_index;

1902

} else if (s->node[j].index < s->node[i].min_index)

1903

s->node[i].min_index = s->node[j].index;

1904

}

1905

1906

if (s->node[i].index != s->node[i].min_index)

return 0;

do {

j = s->stack[--s->sp];

1911

s->node[j].on_stack = 0;

1912

s->order[s->op++] = j;

1913

} while (j != i);

1914

s->order[s->op++] = -1;

return 0;

}

/* Decompose the "len" basic relations in "list" into strongly connected

1920

* components.

1921

*/

1922

static struct basic_map_sort *basic_map_sort_init(isl_ctx *ctx, int len,

1923

__isl_keep isl_basic_map **list)

1924

{

1925

int i;

1926

struct basic_map_sort *s = NULL;

1927

1928

s = basic_map_sort_alloc(ctx, len);

1929

if (!s)

1930

return NULL;

1931

for (i = len - 1; i >= 0; --i) {

1932

if (s->node[i].index >= 0)

1933

continue;

1934

if (power_components_tarjan(s, list, i) < 0)

goto error;

}

return s;

error:

basic_map_sort_free(s);

return NULL;

}

/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D

1945

* and a dimension specification (Z^{n+1} -> Z^{n+1}),

1946

* construct a map that is an overapproximation of the map

1947

* that takes an element from the dom R \times Z to an

1948

* element from ran R \times Z, such that the first n coordinates of the

1949

* difference between them is a sum of differences between images

1950

* and pre-images in one of the R_i and such that the last coordinate

1951

* is equal to the number of steps taken.

1952

* If "project" is set, then these final coordinates are not included,

1953

* i.e., a relation of type Z^n -> Z^n is returned.

1954

* That is, let

1955

*

1956

* \Delta_i = { y - x | (x, y) in R_i }

1957

*

1958

* then the constructed map is an overapproximation of

1959

*

1960

* { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :

1961

* d = (\sum_i k_i \delta_i, \sum_i k_i) and

1962

* x in dom R and x + d in ran R }

*

* or

*

* { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :

1967

* d = (\sum_i k_i \delta_i) and

1968

* x in dom R and x + d in ran R }

1969

*

1970

* if "project" is set.

1971

*

1972

* We first split the map into strongly connected components, perform

1973

* the above on each component and then join the results in the correct

1974

* order, at each join also taking in the union of both arguments

1975

* to allow for paths that do not go through one of the two arguments.

1976

*/

1977

static __isl_give isl_map *construct_power_components(__isl_take isl_space *dim,

1978

__isl_keep isl_map *map, int *exact, int project)

1979

{

1980

int i, n, c;

1981

struct isl_map *path = NULL;

1982

struct basic_map_sort *s = NULL;

int *orig_exact;

int local_exact;

if (!map)

goto error;

if (map->n <= 1)

return floyd_warshall(dim, map, exact, project);

1990

1991

s = basic_map_sort_init(map->ctx, map->n, map->p);

if (!s)

goto error;

orig_exact = exact;

1996

if (s->check_closed && !exact)

1997

exact = &local_exact;

c = 0;

i = 0;

n = map->n;

if (project)

path = isl_map_empty(isl_map_get_space(map));

2004

else

2005

path = isl_map_empty(isl_space_copy(dim));

2006

path = anonymize(path);

2007

while (n) {

2008

struct isl_map *comp;

2009

isl_map *path_comp, *path_comb;

2010

comp = isl_map_alloc_space(isl_map_get_space(map), n, 0);

2011

while (s->order[i] != -1) {

2012

comp = isl_map_add_basic_map(comp,

2013

isl_basic_map_copy(map->p[s->order[i]]));

--n;

++i;

}

path_comp = floyd_warshall(isl_space_copy(dim),

2018

comp, exact, project);

2019

path_comp = anonymize(path_comp);

2020

path_comb = isl_map_apply_range(isl_map_copy(path),

2021

isl_map_copy(path_comp));

2022

path = isl_map_union(path, path_comp);

2023

path = isl_map_union(path, path_comb);

2024

isl_map_free(comp);

++i;

++c;

}

if (c > 1 && s->check_closed && !*exact) {

2030

int closed;

2031

2032

closed = isl_map_is_transitively_closed(path);

if (closed < 0)

goto error;

if (!closed) {

basic_map_sort_free(s);

2037

isl_map_free(path);

2038

return floyd_warshall(dim, map, orig_exact, project);

}

}

basic_map_sort_free(s);

2043

isl_space_free(dim);

return path;

error:

basic_map_sort_free(s);

2048

isl_space_free(dim);

2049

isl_map_free(path);

return NULL;

}

/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,

2054

* construct a map that is an overapproximation of the map

2055

* that takes an element from the space D to another

2056

* element from the same space, such that the difference between

2057

* them is a strictly positive sum of differences between images

2058

* and pre-images in one of the R_i.

2059

* The number of differences in the sum is equated to parameter "param".

2060

* That is, let

2061

*

2062

* \Delta_i = { y - x | (x, y) in R_i }

2063

*

2064

* then the constructed map is an overapproximation of

2065

*

2066

* { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :

2067

* d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }

2068

* or

2069

*

2070

* { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :

2071

* d = \sum_i k_i \delta_i and \sum_i k_i > 0 }

2072

*

2073

* if "project" is set.

2074

*

2075

* If "project" is not set, then

2076

* we construct an extended mapping with an extra coordinate

2077

* that indicates the number of steps taken. In particular,

2078

* the difference in the last coordinate is equal to the number

2079

* of steps taken to move from a domain element to the corresponding

2080

* image element(s).

2081

*/

2082

static __isl_give isl_map *construct_power(__isl_keep isl_map *map,

2083

int *exact, int project)

2084

{

2085

struct isl_map *app = NULL;

2086

isl_space *dim = NULL;

unsigned d;

if (!map)

return NULL;

dim = isl_map_get_space(map);

2093

2094

d = isl_space_dim(dim, isl_dim_in);

2095

dim = isl_space_add_dims(dim, isl_dim_in, 1);

2096

dim = isl_space_add_dims(dim, isl_dim_out, 1);

2097

2098

app = construct_power_components(isl_space_copy(dim), map,

2099

exact, project);

2100

2101

isl_space_free(dim);

return app;

}

/* Compute the positive powers of "map", or an overapproximation.

2107

* If the result is exact, then *exact is set to 1.

2108

*

2109

* If project is set, then we are actually interested in the transitive

2110

* closure, so we can use a more relaxed exactness check.

2111

* The lengths of the paths are also projected out instead of being

2112

* encoded as the difference between an extra pair of final coordinates.

2113

*/

2114

static __isl_give isl_map *map_power(__isl_take isl_map *map,

2115

int *exact, int project)

2116

{

2117

struct isl_map *app = NULL;

if (exact)

*exact = 1;

if (!map)

return NULL;

isl_assert(map->ctx,

2126

isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),

2127

goto error);

2128

2129

app = construct_power(map, exact, project);

isl_map_free(map);

return app;

error:

isl_map_free(map);

isl_map_free(app);

return NULL;

}

/* Compute the positive powers of "map", or an overapproximation.

2140

* The result maps the exponent to a nested copy of the corresponding power.

2141

* If the result is exact, then *exact is set to 1.

2142

* map_power constructs an extended relation with the path lengths

2143

* encoded as the difference between the final coordinates.

2144

* In the final step, this difference is equated to an extra parameter

2145

* and made positive. The extra coordinates are subsequently projected out

2146

* and the parameter is turned into the domain of the result.

2147

*/

2148

__isl_give isl_map *isl_map_power(__isl_take isl_map *map, int *exact)

2149

{

2150

isl_space *target_dim;

isl_space *dim;

isl_map *diff;

unsigned d;

unsigned param;

if (!map)

return NULL;

d = isl_map_dim(map, isl_dim_in);

2160

param = isl_map_dim(map, isl_dim_param);

2161

2162

map = isl_map_compute_divs(map);

2163

map = isl_map_coalesce(map);

2164

2165

if (isl_map_plain_is_empty(map)) {

2166

map = isl_map_from_range(isl_map_wrap(map));

2167

map = isl_map_add_dims(map, isl_dim_in, 1);

2168

map = isl_map_set_dim_name(map, isl_dim_in, 0, "k");

return map;

}

target_dim = isl_map_get_space(map);

2173

target_dim = isl_space_from_range(isl_space_wrap(target_dim));

2174

target_dim = isl_space_add_dims(target_dim, isl_dim_in, 1);

2175

target_dim = isl_space_set_dim_name(target_dim, isl_dim_in, 0, "k");

2176

2177

map = map_power(map, exact, 0);

2178

2179

map = isl_map_add_dims(map, isl_dim_param, 1);

2180

dim = isl_map_get_space(map);

2181

diff = equate_parameter_to_length(dim, param);

2182

map = isl_map_intersect(map, diff);

2183

map = isl_map_project_out(map, isl_dim_in, d, 1);

2184

map = isl_map_project_out(map, isl_dim_out, d, 1);

2185

map = isl_map_from_range(isl_map_wrap(map));

2186

map = isl_map_move_dims(map, isl_dim_in, 0, isl_dim_param, param, 1);

2187

2188

map = isl_map_reset_space(map, target_dim);

return map;

}

/* Compute a relation that maps each element in the range of the input

2194

* relation to the lengths of all paths composed of edges in the input

2195

* relation that end up in the given range element.

2196

* The result may be an overapproximation, in which case *exact is set to 0.

2197

* The resulting relation is very similar to the power relation.

2198

* The difference are that the domain has been projected out, the

2199

* range has become the domain and the exponent is the range instead

2200

* of a parameter.

2201

*/

2202

__isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map,

int *exact)

{

isl_space *dim;

isl_map *diff;

unsigned d;

unsigned param;

if (!map)

return NULL;

d = isl_map_dim(map, isl_dim_in);

2214

param = isl_map_dim(map, isl_dim_param);

2215

2216

map = isl_map_compute_divs(map);

2217

map = isl_map_coalesce(map);

2218

2219

if (isl_map_plain_is_empty(map)) {

2220

if (exact)

2221

*exact = 1;

2222

map = isl_map_project_out(map, isl_dim_out, 0, d);

2223

map = isl_map_add_dims(map, isl_dim_out, 1);

return map;

}

map = map_power(map, exact, 0);

2228

2229

map = isl_map_add_dims(map, isl_dim_param, 1);

2230

dim = isl_map_get_space(map);

2231

diff = equate_parameter_to_length(dim, param);

2232

map = isl_map_intersect(map, diff);

2233

map = isl_map_project_out(map, isl_dim_in, 0, d + 1);

2234

map = isl_map_project_out(map, isl_dim_out, d, 1);

2235

map = isl_map_reverse(map);

2236

map = isl_map_move_dims(map, isl_dim_out, 0, isl_dim_param, param, 1);

return map;

}

/* Check whether equality i of bset is a pure stride constraint

2242

* on a single dimensions, i.e., of the form

*

* v = k e

*

* with k a constant and e an existentially quantified variable.

2247

*/

2248

static int is_eq_stride(__isl_keep isl_basic_set *bset, int i)

{

unsigned nparam;

unsigned d;

unsigned n_div;

int pos1;

int pos2;

if (!bset)

return -1;

if (!isl_int_is_zero(bset->eq[i][0]))

2260

return 0;

2261

2262

nparam = isl_basic_set_dim(bset, isl_dim_param);

2263

d = isl_basic_set_dim(bset, isl_dim_set);

2264

n_div = isl_basic_set_dim(bset, isl_dim_div);

2265

2266

if (isl_seq_first_non_zero(bset->eq[i] + 1, nparam) != -1)

2267

return 0;

2268

pos1 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam, d);

2269

if (pos1 == -1)

2270

return 0;

2271

if (isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + pos1 + 1,

2272

d - pos1 - 1) != -1)

2273

return 0;

2274

2275

pos2 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d, n_div);

2276

if (pos2 == -1)

2277

return 0;

2278

if (isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d + pos2 + 1,

2279

n_div - pos2 - 1) != -1)

2280

return 0;

2281

if (!isl_int_is_one(bset->eq[i][1 + nparam + pos1]) &&

2282

!isl_int_is_negone(bset->eq[i][1 + nparam + pos1]))

return 0;

return 1;

}

/* Given a map, compute the smallest superset of this map that is of the form

2289

*

2290

* { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }

2291

*

2292

* (where p ranges over the (non-parametric) dimensions),

2293

* compute the transitive closure of this map, i.e.,

2294

*

2295

* { i -> j : exists k > 0:

2296

* k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }

2297

*

2298

* and intersect domain and range of this transitive closure with

2299

* the given domain and range.

2300

*

2301

* If with_id is set, then try to include as much of the identity mapping

2302

* as possible, by computing

2303

*

2304

* { i -> j : exists k >= 0:

2305

* k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }

2306

*

2307

* instead (i.e., allow k = 0).

2308

*

2309

* In practice, we compute the difference set

2310

*

2311

* delta = { j - i | i -> j in map },

2312

*

2313

* look for stride constraint on the individual dimensions and compute

2314

* (constant) lower and upper bounds for each individual dimension,

2315

* adding a constraint for each bound not equal to infinity.

2316

*/

2317

static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map,

2318

__isl_take isl_set *dom, __isl_take isl_set *ran, int with_id)

{

int i;

int k;

unsigned d;

unsigned nparam;

unsigned total;

isl_space *dim;

isl_set *delta;

isl_map *app = NULL;

2328

isl_basic_set *aff = NULL;

2329

isl_basic_map *bmap = NULL;

2330

isl_vec *obj = NULL;

isl_int opt;

isl_int_init(opt);

delta = isl_map_deltas(isl_map_copy(map));

2336

2337

aff = isl_set_affine_hull(isl_set_copy(delta));

2338

if (!aff)

2339

goto error;

2340

dim = isl_map_get_space(map);

2341

d = isl_space_dim(dim, isl_dim_in);

2342

nparam = isl_space_dim(dim, isl_dim_param);

2343

total = isl_space_dim(dim, isl_dim_all);

2344

bmap = isl_basic_map_alloc_space(dim,

2345

aff->n_div + 1, aff->n_div, 2 * d + 1);

2346

for (i = 0; i < aff->n_div + 1; ++i) {

2347

k = isl_basic_map_alloc_div(bmap);

2348

if (k < 0)

2349

goto error;

2350

isl_int_set_si(bmap->div[k][0], 0);

2351

}

2352

for (i = 0; i < aff->n_eq; ++i) {

2353

if (!is_eq_stride(aff, i))

2354

continue;

2355

k = isl_basic_map_alloc_equality(bmap);

2356

if (k < 0)

2357

goto error;

2358

isl_seq_clr(bmap->eq[k], 1 + nparam);

2359

isl_seq_cpy(bmap->eq[k] + 1 + nparam + d,

2360

aff->eq[i] + 1 + nparam, d);

2361

isl_seq_neg(bmap->eq[k] + 1 + nparam,

2362

aff->eq[i] + 1 + nparam, d);

2363

isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d,

2364

aff->eq[i] + 1 + nparam + d, aff->n_div);

2365

isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0);

2366

}

2367

obj = isl_vec_alloc(map->ctx, 1 + nparam + d);

2368

if (!obj)

2369

goto error;

2370

isl_seq_clr(obj->el, 1 + nparam + d);

2371

for (i = 0; i < d; ++ i) {

2372

enum isl_lp_result res;

2373

2374

isl_int_set_si(obj->el[1 + nparam + i], 1);

2375

2376

res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt,

2377

NULL, NULL);

2378

if (res == isl_lp_error)

2379

goto error;

2380

if (res == isl_lp_ok) {

2381

k = isl_basic_map_alloc_inequality(bmap);

2382

if (k < 0)

2383

goto error;

2384

isl_seq_clr(bmap->ineq[k],

2385

1 + nparam + 2 * d + bmap->n_div);

2386

isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1);

2387

isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1);

2388

isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);

2389

}

2390

2391

res = isl_set_solve_lp(delta, 1, obj->el, map->ctx->one, &opt,

2392

NULL, NULL);

2393

if (res == isl_lp_error)

2394

goto error;

2395

if (res == isl_lp_ok) {

2396

k = isl_basic_map_alloc_inequality(bmap);

2397

if (k < 0)

2398

goto error;

2399

isl_seq_clr(bmap->ineq[k],

2400

1 + nparam + 2 * d + bmap->n_div);

2401

isl_int_set_si(bmap->ineq[k][1 + nparam + i], 1);

2402

isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1);

2403

isl_int_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);

2404

}

2405

2406

isl_int_set_si(obj->el[1 + nparam + i], 0);

2407

}

2408

k = isl_basic_map_alloc_inequality(bmap);

2409

if (k < 0)

2410

goto error;

2411

isl_seq_clr(bmap->ineq[k],

2412

1 + nparam + 2 * d + bmap->n_div);

2413

if (!with_id)

2414

isl_int_set_si(bmap->ineq[k][0], -1);

2415

isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1);

2416

2417

app = isl_map_from_domain_and_range(dom, ran);

2418

2419

isl_vec_free(obj);

2420

isl_basic_set_free(aff);

2421

isl_map_free(map);

2422

bmap = isl_basic_map_finalize(bmap);

2423

isl_set_free(delta);

2424

isl_int_clear(opt);

2425

2426

map = isl_map_from_basic_map(bmap);

2427

map = isl_map_intersect(map, app);

return map;

error:

isl_vec_free(obj);

isl_basic_map_free(bmap);

2433

isl_basic_set_free(aff);

isl_set_free(dom);

isl_set_free(ran);

isl_map_free(map);

isl_set_free(delta);

2438

isl_int_clear(opt);

return NULL;

}

/* Given a map, compute the smallest superset of this map that is of the form

2443

*

2444

* { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }

2445

*

2446

* (where p ranges over the (non-parametric) dimensions),

2447

* compute the transitive closure of this map, i.e.,

2448

*

2449

* { i -> j : exists k > 0:

2450

* k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }

2451

*

2452

* and intersect domain and range of this transitive closure with

2453

* domain and range of the original map.

2454

*/

2455

static __isl_give isl_map *box_closure(__isl_take isl_map *map)

{

isl_set *domain;

isl_set *range;

domain = isl_map_domain(isl_map_copy(map));

2461

domain = isl_set_coalesce(domain);

2462

range = isl_map_range(isl_map_copy(map));

2463

range = isl_set_coalesce(range);

2464

2465

return box_closure_on_domain(map, domain, range, 0);

2466

}

2467

2468

/* Given a map, compute the smallest superset of this map that is of the form

2469

*

2470

* { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }

2471

*

2472

* (where p ranges over the (non-parametric) dimensions),

2473

* compute the transitive and partially reflexive closure of this map, i.e.,

2474

*

2475

* { i -> j : exists k >= 0:

2476

* k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }

2477

*

2478

* and intersect domain and range of this transitive closure with

2479

* the given domain.

2480

*/

2481

static __isl_give isl_map *box_closure_with_identity(__isl_take isl_map *map,

2482

__isl_take isl_set *dom)

2483

{

2484

return box_closure_on_domain(map, dom, isl_set_copy(dom), 1);

2485

}

2486

2487

/* Check whether app is the transitive closure of map.

2488

* In particular, check that app is acyclic and, if so,

2489

* check that

2490

*

2491

* app \subset (map \cup (map \circ app))

2492

*/

2493

static int check_exactness_omega(__isl_keep isl_map *map,

2494

__isl_keep isl_map *app)

{

isl_set *delta;

int i;

int is_empty, is_exact;

unsigned d;

isl_map *test;

delta = isl_map_deltas(isl_map_copy(app));

2503

d = isl_set_dim(delta, isl_dim_set);

2504

for (i = 0; i < d; ++i)

2505

delta = isl_set_fix_si(delta, isl_dim_set, i, 0);

2506

is_empty = isl_set_is_empty(delta);

2507

isl_set_free(delta);

if (is_empty < 0)

return -1;

if (!is_empty)

return 0;

test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map));

2514

test = isl_map_union(test, isl_map_copy(map));

2515

is_exact = isl_map_is_subset(app, test);

2516

isl_map_free(test);

return is_exact;

}

/* Check if basic map M_i can be combined with all the other

2522

* basic maps such that

*

* (\cup_j M_j)^+

*

* can be computed as

2527

*

2528

* M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+

2529

*

2530

* In particular, check if we can compute a compact representation

2531

* of

2532

*

2533

* M_i^* \circ M_j \circ M_i^*

2534

*

2535

* for each j != i.

2536

* Let M_i^? be an extension of M_i^+ that allows paths

2537

* of length zero, i.e., the result of box_closure(., 1).

2538

* The criterion, as proposed by Kelly et al., is that

2539

* id = M_i^? - M_i^+ can be represented as a basic map

2540

* and that

2541

*

2542

* id \circ M_j \circ id = M_j

*

* for each j != i.

*

* If this function returns 1, then tc and qc are set to

2547

* M_i^+ and M_i^?, respectively.

2548

*/

2549

static int can_be_split_off(__isl_keep isl_map *map, int i,

2550

__isl_give isl_map **tc, __isl_give isl_map **qc)

2551

{

2552

isl_map *map_i, *id = NULL;

int j = -1;

isl_set *C;

*tc = NULL;

*qc = NULL;

C = isl_set_union(isl_map_domain(isl_map_copy(map)),

2560

isl_map_range(isl_map_copy(map)));

2561

C = isl_set_from_basic_set(isl_set_simple_hull(C));

if (!C)

goto error;

map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));

2566

*tc = box_closure(isl_map_copy(map_i));

2567

*qc = box_closure_with_identity(map_i, C);

2568

id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc));

if (!id || !*qc)

goto error;

if (id->n != 1 || (*qc)->n != 1)

2573

goto done;

2574

2575

for (j = 0; j < map->n; ++j) {

2576

isl_map *map_j, *test;

int is_ok;

if (i == j)

continue;

map_j = isl_map_from_basic_map(

2582

isl_basic_map_copy(map->p[j]));

2583

test = isl_map_apply_range(isl_map_copy(id),

2584

isl_map_copy(map_j));

2585

test = isl_map_apply_range(test, isl_map_copy(id));

2586

is_ok = isl_map_is_equal(test, map_j);

2587

isl_map_free(map_j);

2588

isl_map_free(test);

if (is_ok < 0)

goto error;

if (!is_ok)

break;

}

done:

isl_map_free(id);

if (j == map->n)

return 1;

isl_map_free(*qc);

isl_map_free(*tc);

*qc = NULL;

*tc = NULL;

return 0;

error:

isl_map_free(id);

isl_map_free(*qc);

isl_map_free(*tc);

*qc = NULL;

*tc = NULL;

return -1;

}

static __isl_give isl_map *box_closure_with_check(__isl_take isl_map *map,

int *exact)

{

isl_map *app;

app = box_closure(isl_map_copy(map));

2621

if (exact)

2622

*exact = check_exactness_omega(map, app);

isl_map_free(map);

return app;

}

/* Compute an overapproximation of the transitive closure of "map"

2629

* using a variation of the algorithm from

2630

* "Transitive Closure of Infinite Graphs and its Applications"

2631

* by Kelly et al.

2632

*

2633

* We first check whether we can can split of any basic map M_i and

* compute

*

* (\cup_j M_j)^+

*

* as

*

* M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+

2641

*

2642

* using a recursive call on the remaining map.

2643

*

2644

* If not, we simply call box_closure on the whole map.

2645

*/

2646

static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map,

int *exact)

{

int i, j;

int exact_i;

isl_map *app;

if (!map)

return NULL;

if (map->n == 1)

return box_closure_with_check(map, exact);

2657

2658

for (i = 0; i < map->n; ++i) {

2659

int ok;

2660

isl_map *qc, *tc;

2661

ok = can_be_split_off(map, i, &tc, &qc);

if (ok < 0)

goto error;

if (!ok)

continue;

app = isl_map_alloc_space(isl_map_get_space(map), map->n - 1, 0);

2668

2669

for (j = 0; j < map->n; ++j) {

2670

if (j == i)

2671

continue;

2672

app = isl_map_add_basic_map(app,

2673

isl_basic_map_copy(map->p[j]));

2674

}

2675

2676

app = isl_map_apply_range(isl_map_copy(qc), app);

2677

app = isl_map_apply_range(app, qc);

2678

2679

app = isl_map_union(tc, transitive_closure_omega(app, NULL));

2680

exact_i = check_exactness_omega(map, app);

2681

if (exact_i == 1) {

if (exact)

*exact = exact_i;

isl_map_free(map);

return app;

}

isl_map_free(app);

if (exact_i < 0)

goto error;

}

return box_closure_with_check(map, exact);

error:

isl_map_free(map);

return NULL;

}

/* Compute the transitive closure of "map", or an overapproximation.

2699

* If the result is exact, then *exact is set to 1.

2700

* Simply use map_power to compute the powers of map, but tell

2701

* it to project out the lengths of the paths instead of equating

2702

* the length to a parameter.

2703

*/

2704

__isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,

2705

int *exact)

2706

{

2707

isl_space *target_dim;

int closed;

if (!map)

goto error;

if (map->ctx->opt->closure == ISL_CLOSURE_BOX)

2714

return transitive_closure_omega(map, exact);

2715

2716

map = isl_map_compute_divs(map);

2717

map = isl_map_coalesce(map);

2718

closed = isl_map_is_transitively_closed(map);

if (closed < 0)

goto error;

if (closed) {

if (exact)

*exact = 1;

return map;

}

target_dim = isl_map_get_space(map);

2728

map = map_power(map, exact, 1);

2729

map = isl_map_reset_space(map, target_dim);

return map;

error:

isl_map_free(map);

return NULL;

}

static int inc_count(__isl_take isl_map *map, void *user)

{

int *n = user;

*n += map->n;

isl_map_free(map);

return 0;

}

static int collect_basic_map(__isl_take isl_map *map, void *user)

2749

{

2750

int i;

2751

isl_basic_map ***next = user;

2752

2753

for (i = 0; i < map->n; ++i) {

2754

**next = isl_basic_map_copy(map->p[i]);

if (!**next)

goto error;

(*next)++;

}

isl_map_free(map);

return 0;

error:

isl_map_free(map);

return -1;

}

/* Perform Floyd-Warshall on the given list of basic relations.

2768

* The basic relations may live in different dimensions,

2769

* but basic relations that get assigned to the diagonal of the

2770

* grid have domains and ranges of the same dimension and so

2771

* the standard algorithm can be used because the nested transitive

2772

* closures are only applied to diagonal elements and because all

2773

* compositions are peformed on relations with compatible domains and ranges.

2774

*/

2775

static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx,

2776

__isl_keep isl_basic_map **list, int n, int *exact)

{

int i, j, k;

int n_group;

int *group = NULL;

isl_set **set = NULL;

2782

isl_map ***grid = NULL;

2783

isl_union_map *app;

2784

2785

group = setup_groups(ctx, list, n, &set, &n_group);

if (!group)

goto error;

grid = isl_calloc_array(ctx, isl_map **, n_group);

2790

if (!grid)

2791

goto error;

2792

for (i = 0; i < n_group; ++i) {

2793

grid[i] = isl_calloc_array(ctx, isl_map *, n_group);

2794

if (!grid[i])

2795

goto error;

2796

for (j = 0; j < n_group; ++j) {

2797

isl_space *dim1, *dim2, *dim;

2798

dim1 = isl_space_reverse(isl_set_get_space(set[i]));

2799

dim2 = isl_set_get_space(set[j]);

2800

dim = isl_space_join(dim1, dim2);

2801

grid[i][j] = isl_map_empty(dim);

}

}

for (k = 0; k < n; ++k) {

2806

i = group[2 * k];

2807

j = group[2 * k + 1];

2808

grid[i][j] = isl_map_union(grid[i][j],

2809

isl_map_from_basic_map(

2810

isl_basic_map_copy(list[k])));

2811

}

2812

2813

floyd_warshall_iterate(grid, n_group, exact);

2814

2815

app = isl_union_map_empty(isl_map_get_space(grid[0][0]));

2816

2817

for (i = 0; i < n_group; ++i) {

2818

for (j = 0; j < n_group; ++j)

2819

app = isl_union_map_add_map(app, grid[i][j]);

free(grid[i]);

}

free(grid);

for (i = 0; i < 2 * n; ++i)

2825

isl_set_free(set[i]);

free(set);

free(group);

return app;

error:

if (grid)

for (i = 0; i < n_group; ++i) {

2833

if (!grid[i])

2834

continue;

2835

for (j = 0; j < n_group; ++j)

2836

isl_map_free(grid[i][j]);

free(grid[i]);

}

free(grid);

if (set) {

for (i = 0; i < 2 * n; ++i)

2842

isl_set_free(set[i]);

free(set);

}

free(group);

return NULL;

}

/* Perform Floyd-Warshall on the given union relation.

2850

* The implementation is very similar to that for non-unions.

2851

* The main difference is that it is applied unconditionally.

2852

* We first extract a list of basic maps from the union map

2853

* and then perform the algorithm on this list.

2854

*/

2855

static __isl_give isl_union_map *union_floyd_warshall(

2856

__isl_take isl_union_map *umap, int *exact)

{

int i, n;

isl_ctx *ctx;

isl_basic_map **list = NULL;

2861

isl_basic_map **next;

2862

isl_union_map *res;

2863

2864

n = 0;

2865

if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)

2866

goto error;

2867

2868

ctx = isl_union_map_get_ctx(umap);

2869

list = isl_calloc_array(ctx, isl_basic_map *, n);

if (!list)

goto error;

next = list;

if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)

2875

goto error;

2876

2877

res = union_floyd_warshall_on_list(ctx, list, n, exact);

2878

2879

if (list) {

2880

for (i = 0; i < n; ++i)

2881

isl_basic_map_free(list[i]);

free(list);

}

isl_union_map_free(umap);

return res;

error:

if (list) {

for (i = 0; i < n; ++i)

2890

isl_basic_map_free(list[i]);

2891

free(list);

2892

}

2893

isl_union_map_free(umap);

return NULL;

}

/* Decompose the give union relation into strongly connected components.

2898

* The implementation is essentially the same as that of

2899

* construct_power_components with the major difference that all

2900

* operations are performed on union maps.

2901

*/

2902

static __isl_give isl_union_map *union_components(

2903

__isl_take isl_union_map *umap, int *exact)

{

int i;

int n;

isl_ctx *ctx;

isl_basic_map **list;

2909

isl_basic_map **next;

2910

isl_union_map *path = NULL;

2911

struct basic_map_sort *s = NULL;

int c, l;

int recheck = 0;

n = 0;

if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)

goto error;

if (n <= 1)

return union_floyd_warshall(umap, exact);

2921

2922

ctx = isl_union_map_get_ctx(umap);

2923

list = isl_calloc_array(ctx, isl_basic_map *, n);

if (!list)

goto error;

next = list;

if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)

2929

goto error;

2930

2931

s = basic_map_sort_init(ctx, n, list);

if (!s)

goto error;

c = 0;

i = 0;

l = n;

path = isl_union_map_empty(isl_union_map_get_space(umap));

2939

while (l) {

2940

isl_union_map *comp;

2941

isl_union_map *path_comp, *path_comb;

2942

comp = isl_union_map_empty(isl_union_map_get_space(umap));

2943

while (s->order[i] != -1) {

2944

comp = isl_union_map_add_map(comp,

2945

isl_map_from_basic_map(

2946

isl_basic_map_copy(list[s->order[i]])));

--l;

++i;

}

path_comp = union_floyd_warshall(comp, exact);

2951

path_comb = isl_union_map_apply_range(isl_union_map_copy(path),

2952

isl_union_map_copy(path_comp));

2953

path = isl_union_map_union(path, path_comp);

2954

path = isl_union_map_union(path, path_comb);

++i;

++c;

}

if (c > 1 && s->check_closed && !*exact) {

2960

int closed;

2961

2962

closed = isl_union_map_is_transitively_closed(path);

if (closed < 0)

goto error;

recheck = !closed;

}

basic_map_sort_free(s);

2969

2970

for (i = 0; i < n; ++i)

2971

isl_basic_map_free(list[i]);

free(list);

if (recheck) {

isl_union_map_free(path);

2976

return union_floyd_warshall(umap, exact);

2977

}

2978

2979

isl_union_map_free(umap);

return path;

error:

basic_map_sort_free(s);

2984

if (list) {

2985

for (i = 0; i < n; ++i)

2986

isl_basic_map_free(list[i]);

2987

free(list);

2988

}

2989

isl_union_map_free(umap);

2990

isl_union_map_free(path);

return NULL;

}

/* Compute the transitive closure of "umap", or an overapproximation.

2995

* If the result is exact, then *exact is set to 1.

2996

*/

2997

__isl_give isl_union_map *isl_union_map_transitive_closure(

2998

__isl_take isl_union_map *umap, int *exact)

{

int closed;

if (!umap)

return NULL;

if (exact)

*exact = 1;

umap = isl_union_map_compute_divs(umap);

3009

umap = isl_union_map_coalesce(umap);

3010

closed = isl_union_map_is_transitively_closed(umap);

if (closed < 0)

goto error;

if (closed)

return umap;

umap = union_components(umap, exact);

3016

return umap;

3017

error:

3018

isl_union_map_free(umap);

return NULL;

}

struct isl_union_power {

3023

isl_union_map *pow;

int *exact;

};

static int power(__isl_take isl_map *map, void *user)

3028

{

3029

struct isl_union_power *up = user;

3030

3031

map = isl_map_power(map, up->exact);

3032

up->pow = isl_union_map_from_map(map);

return -1;

}

/* Construct a map [x] -> [x+1], with parameters prescribed by "dim".

3038

*/

3039

static __isl_give isl_union_map *increment(__isl_take isl_space *dim)

3040

{

3041

int k;

3042

isl_basic_map *bmap;

3043

3044

dim = isl_space_add_dims(dim, isl_dim_in, 1);

3045

dim = isl_space_add_dims(dim, isl_dim_out, 1);

3046

bmap = isl_basic_map_alloc_space(dim, 0, 1, 0);

3047

k = isl_basic_map_alloc_equality(bmap);

3048

if (k < 0)

3049

goto error;

3050

isl_seq_clr(bmap->eq[k], isl_basic_map_total_dim(bmap));

3051

isl_int_set_si(bmap->eq[k][0], 1);

3052

isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1);

3053

isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1);

3054

return isl_union_map_from_map(isl_map_from_basic_map(bmap));

3055

error:

3056

isl_basic_map_free(bmap);

return NULL;

}

/* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "dim".

3061

*/

3062

static __isl_give isl_union_map *deltas_map(__isl_take isl_space *dim)

3063

{

3064

isl_basic_map *bmap;

3065

3066

dim = isl_space_add_dims(dim, isl_dim_in, 1);

3067

dim = isl_space_add_dims(dim, isl_dim_out, 1);

3068

bmap = isl_basic_map_universe(dim);

3069

bmap = isl_basic_map_deltas_map(bmap);

3070

3071

return isl_union_map_from_map(isl_map_from_basic_map(bmap));

3072

}

3073

3074

/* Compute the positive powers of "map", or an overapproximation.

3075

* The result maps the exponent to a nested copy of the corresponding power.

3076

* If the result is exact, then *exact is set to 1.

3077

*/

3078

__isl_give isl_union_map *isl_union_map_power(__isl_take isl_union_map *umap,

int *exact)

{

int n;

isl_union_map *inc;

isl_union_map *dm;

if (!umap)

return NULL;

n = isl_union_map_n_map(umap);

if (n == 0)

return umap;

if (n == 1) {

struct isl_union_power up = { NULL, exact };

3092

isl_union_map_foreach_map(umap, &power, &up);

3093

isl_union_map_free(umap);

3094

return up.pow;

3095

}

3096

inc = increment(isl_union_map_get_space(umap));

3097

umap = isl_union_map_product(inc, umap);

3098

umap = isl_union_map_transitive_closure(umap, exact);

3099

umap = isl_union_map_zip(umap);

3100

dm = deltas_map(isl_union_map_get_space(umap));

3101

umap = isl_union_map_apply_domain(umap, dm);

return umap;

}

#undef TYPE

#define TYPE isl_map

3108

#include "isl_power_templ.c"

3109

3110

#undef TYPE

3111

#define TYPE isl_union_map

3112

#include "isl_power_templ.c"

nexmon – Blame information for rev 1