Coverage Report

Created: 2018-08-14 02:14

/Users/buildslave/jenkins/workspace/clang-stage2-coverage-R/llvm/tools/polly/lib/External/isl/isl_mat.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 * Copyright 2010      INRIA Saclay
4
 * Copyright 2014      Ecole Normale Superieure
5
 * Copyright 2017      Sven Verdoolaege
6
 *
7
 * Use of this software is governed by the MIT license
8
 *
9
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
10
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
11
 * and INRIA Saclay - Ile-de-France, Parc Club Orsay Universite,
12
 * ZAC des vignes, 4 rue Jacques Monod, 91893 Orsay, France
13
 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
14
 */
15
16
#include <isl_ctx_private.h>
17
#include <isl_map_private.h>
18
#include <isl/space.h>
19
#include <isl_seq.h>
20
#include <isl_mat_private.h>
21
#include <isl_vec_private.h>
22
#include <isl_space_private.h>
23
#include <isl_val_private.h>
24
25
isl_ctx *isl_mat_get_ctx(__isl_keep isl_mat *mat)
26
3.54M
{
27
3.54M
  return mat ? mat->ctx : NULL;
28
3.54M
}
29
30
/* Return a hash value that digests "mat".
31
 */
32
uint32_t isl_mat_get_hash(__isl_keep isl_mat *mat)
33
0
{
34
0
  int i;
35
0
  uint32_t hash;
36
0
37
0
  if (!mat)
38
0
    return 0;
39
0
40
0
  hash = isl_hash_init();
41
0
  isl_hash_byte(hash, mat->n_row & 0xFF);
42
0
  isl_hash_byte(hash, mat->n_col & 0xFF);
43
0
  for (i = 0; i < mat->n_row; ++i) {
44
0
    uint32_t row_hash;
45
0
46
0
    row_hash = isl_seq_get_hash(mat->row[i], mat->n_col);
47
0
    isl_hash_hash(hash, row_hash);
48
0
  }
49
0
50
0
  return hash;
51
0
}
52
53
struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
54
  unsigned n_row, unsigned n_col)
55
5.87M
{
56
5.87M
  int i;
57
5.87M
  struct isl_mat *mat;
58
5.87M
59
5.87M
  mat = isl_alloc_type(ctx, struct isl_mat);
60
5.87M
  if (!mat)
61
0
    return NULL;
62
5.87M
63
5.87M
  mat->row = NULL;
64
5.87M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
65
5.87M
  if (isl_blk_is_error(mat->block))
66
0
    goto error;
67
5.87M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
68
5.87M
  if (n_row && 
!mat->row4.36M
)
69
0
    goto error;
70
5.87M
71
32.1M
  
for (i = 0; 5.87M
i < n_row;
++i26.3M
)
72
26.3M
    mat->row[i] = mat->block.data + i * n_col;
73
5.87M
74
5.87M
  mat->ctx = ctx;
75
5.87M
  isl_ctx_ref(ctx);
76
5.87M
  mat->ref = 1;
77
5.87M
  mat->n_row = n_row;
78
5.87M
  mat->n_col = n_col;
79
5.87M
  mat->max_col = n_col;
80
5.87M
  mat->flags = 0;
81
5.87M
82
5.87M
  return mat;
83
0
error:
84
0
  isl_blk_free(ctx, mat->block);
85
0
  free(mat);
86
0
  return NULL;
87
5.87M
}
88
89
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
90
  unsigned n_row, unsigned n_col)
91
102k
{
92
102k
  int i;
93
102k
  isl_int *old;
94
102k
  isl_int **row;
95
102k
96
102k
  if (!mat)
97
0
    return NULL;
98
102k
99
102k
  if (mat->max_col >= n_col && 
mat->n_row >= n_row99.0k
) {
100
7.80k
    if (mat->n_col < n_col)
101
42
      mat->n_col = n_col;
102
7.80k
    return mat;
103
7.80k
  }
104
95.1k
105
95.1k
  if (mat->max_col < n_col) {
106
3.84k
    struct isl_mat *new_mat;
107
3.84k
108
3.84k
    if (n_row < mat->n_row)
109
9
      n_row = mat->n_row;
110
3.84k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
111
3.84k
    if (!new_mat)
112
0
      goto error;
113
29.9k
    
for (i = 0; 3.84k
i < mat->n_row;
++i26.0k
)
114
26.0k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
115
3.84k
    isl_mat_free(mat);
116
3.84k
    return new_mat;
117
3.84k
  }
118
91.2k
119
91.2k
  mat = isl_mat_cow(mat);
120
91.2k
  if (!mat)
121
0
    goto error;
122
91.2k
123
91.2k
  old = mat->block.data;
124
91.2k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
125
91.2k
  if (isl_blk_is_error(mat->block))
126
0
    goto error;
127
91.2k
  row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
128
91.2k
  if (n_row && !row)
129
0
    goto error;
130
91.2k
  mat->row = row;
131
91.2k
132
936k
  for (i = 0; i < mat->n_row; 
++i845k
)
133
845k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
134
318k
  for (i = mat->n_row; i < n_row; 
++i227k
)
135
227k
    mat->row[i] = mat->block.data + i * mat->max_col;
136
91.2k
  mat->n_row = n_row;
137
91.2k
  if (mat->n_col < n_col)
138
0
    mat->n_col = n_col;
139
91.2k
140
91.2k
  return mat;
141
0
error:
142
0
  isl_mat_free(mat);
143
0
  return NULL;
144
91.2k
}
145
146
__isl_give isl_mat *isl_mat_sub_alloc6(isl_ctx *ctx, isl_int **row,
147
  unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
148
2.79M
{
149
2.79M
  int i;
150
2.79M
  struct isl_mat *mat;
151
2.79M
152
2.79M
  mat = isl_alloc_type(ctx, struct isl_mat);
153
2.79M
  if (!mat)
154
0
    return NULL;
155
2.79M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
156
2.79M
  if (n_row && 
!mat->row1.57M
)
157
0
    goto error;
158
9.58M
  
for (i = 0; 2.79M
i < n_row;
++i6.78M
)
159
6.78M
    mat->row[i] = row[first_row+i] + first_col;
160
2.79M
  mat->ctx = ctx;
161
2.79M
  isl_ctx_ref(ctx);
162
2.79M
  mat->ref = 1;
163
2.79M
  mat->n_row = n_row;
164
2.79M
  mat->n_col = n_col;
165
2.79M
  mat->block = isl_blk_empty();
166
2.79M
  mat->flags = ISL_MAT_BORROWED;
167
2.79M
  return mat;
168
0
error:
169
0
  free(mat);
170
0
  return NULL;
171
2.79M
}
172
173
__isl_give isl_mat *isl_mat_sub_alloc(__isl_keep isl_mat *mat,
174
  unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
175
651k
{
176
651k
  if (!mat)
177
0
    return NULL;
178
651k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
179
651k
          first_col, n_col);
180
651k
}
181
182
void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
183
  unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
184
409k
{
185
409k
  int i;
186
409k
187
1.34M
  for (i = 0; i < n_row; 
++i932k
)
188
932k
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
189
409k
}
190
191
void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
192
  unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
193
147k
{
194
147k
  int i;
195
147k
196
454k
  for (i = 0; i < n_row; 
++i306k
)
197
306k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
198
147k
}
199
200
__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
201
2.09M
{
202
2.09M
  if (!mat)
203
0
    return NULL;
204
2.09M
205
2.09M
  mat->ref++;
206
2.09M
  return mat;
207
2.09M
}
208
209
__isl_give isl_mat *isl_mat_dup(__isl_keep isl_mat *mat)
210
366k
{
211
366k
  int i;
212
366k
  struct isl_mat *mat2;
213
366k
214
366k
  if (!mat)
215
2.67k
    return NULL;
216
363k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
217
363k
  if (!mat2)
218
0
    return NULL;
219
1.17M
  
for (i = 0; 363k
i < mat->n_row;
++i811k
)
220
811k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
221
363k
  return mat2;
222
363k
}
223
224
__isl_give isl_mat *isl_mat_cow(__isl_take isl_mat *mat)
225
3.95M
{
226
3.95M
  struct isl_mat *mat2;
227
3.95M
  if (!mat)
228
0
    return NULL;
229
3.95M
230
3.95M
  if (mat->ref == 1 && 
!3.91M
ISL_F_ISSET3.91M
(mat, ISL_MAT_BORROWED))
231
3.95M
    
return mat3.59M
;
232
360k
233
360k
  mat2 = isl_mat_dup(mat);
234
360k
  isl_mat_free(mat);
235
360k
  return mat2;
236
360k
}
237
238
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
239
12.1M
{
240
12.1M
  if (!mat)
241
1.36M
    return NULL;
242
10.7M
243
10.7M
  if (--mat->ref > 0)
244
2.09M
    return NULL;
245
8.66M
246
8.66M
  if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
247
8.66M
    
isl_blk_free(mat->ctx, mat->block)5.86M
;
248
8.66M
  isl_ctx_deref(mat->ctx);
249
8.66M
  free(mat->row);
250
8.66M
  free(mat);
251
8.66M
252
8.66M
  return NULL;
253
8.66M
}
254
255
int isl_mat_rows(__isl_keep isl_mat *mat)
256
493k
{
257
493k
  return mat ? mat->n_row : 
-10
;
258
493k
}
259
260
int isl_mat_cols(__isl_keep isl_mat *mat)
261
55.6k
{
262
55.6k
  return mat ? mat->n_col : 
-10
;
263
55.6k
}
264
265
/* Check that "col" is a valid column position for "mat".
266
 */
267
static isl_stat check_col(__isl_keep isl_mat *mat, int col)
268
4.96k
{
269
4.96k
  if (!mat)
270
0
    return isl_stat_error;
271
4.96k
  if (col < 0 || col >= mat->n_col)
272
4.96k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
273
4.96k
      "column out of range", return isl_stat_error);
274
4.96k
  return isl_stat_ok;
275
4.96k
}
276
277
/* Check that "row" is a valid row position for "mat".
278
 */
279
static isl_stat check_row(__isl_keep isl_mat *mat, int row)
280
4.97k
{
281
4.97k
  if (!mat)
282
0
    return isl_stat_error;
283
4.97k
  if (row < 0 || row >= mat->n_row)
284
4.97k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
285
4.97k
      "row out of range", return isl_stat_error);
286
4.97k
  return isl_stat_ok;
287
4.97k
}
288
289
/* Check that there are "n" columns starting at position "first" in "mat".
290
 */
291
static isl_stat check_col_range(__isl_keep isl_mat *mat, unsigned first,
292
  unsigned n)
293
1.57M
{
294
1.57M
  if (!mat)
295
0
    return isl_stat_error;
296
1.57M
  if (first + n > mat->n_col || first + n < first)
297
1.57M
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
298
1.57M
      "column position or range out of bounds",
299
1.57M
      return isl_stat_error);
300
1.57M
  return isl_stat_ok;
301
1.57M
}
302
303
/* Check that there are "n" rows starting at position "first" in "mat".
304
 */
305
static isl_stat check_row_range(__isl_keep isl_mat *mat, unsigned first,
306
  unsigned n)
307
4.11M
{
308
4.11M
  if (!mat)
309
0
    return isl_stat_error;
310
4.11M
  if (first + n > mat->n_row || first + n < first)
311
4.11M
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
312
4.11M
      "row position or range out of bounds",
313
4.11M
      return isl_stat_error);
314
4.11M
  return isl_stat_ok;
315
4.11M
}
316
317
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
318
3.39k
{
319
3.39k
  if (check_row(mat, row) < 0)
320
0
    return -1;
321
3.39k
  if (check_col(mat, col) < 0)
322
0
    return -1;
323
3.39k
  isl_int_set(*v, mat->row[row][col]);
324
3.39k
  return 0;
325
3.39k
}
326
327
/* Extract the element at row "row", oolumn "col" of "mat".
328
 */
329
__isl_give isl_val *isl_mat_get_element_val(__isl_keep isl_mat *mat,
330
  int row, int col)
331
24
{
332
24
  isl_ctx *ctx;
333
24
334
24
  if (check_row(mat, row) < 0)
335
0
    return NULL;
336
24
  if (check_col(mat, col) < 0)
337
0
    return NULL;
338
24
  ctx = isl_mat_get_ctx(mat);
339
24
  return isl_val_int_from_isl_int(ctx, mat->row[row][col]);
340
24
}
341
342
__isl_give isl_mat *isl_mat_set_element(__isl_take isl_mat *mat,
343
  int row, int col, isl_int v)
344
1.48k
{
345
1.48k
  mat = isl_mat_cow(mat);
346
1.48k
  if (check_row(mat, row) < 0)
347
0
    return isl_mat_free(mat);
348
1.48k
  if (check_col(mat, col) < 0)
349
0
    return isl_mat_free(mat);
350
1.48k
  isl_int_set(mat->row[row][col], v);
351
1.48k
  return mat;
352
1.48k
}
353
354
__isl_give isl_mat *isl_mat_set_element_si(__isl_take isl_mat *mat,
355
  int row, int col, int v)
356
43
{
357
43
  mat = isl_mat_cow(mat);
358
43
  if (check_row(mat, row) < 0)
359
0
    return isl_mat_free(mat);
360
43
  if (check_col(mat, col) < 0)
361
0
    return isl_mat_free(mat);
362
43
  isl_int_set_si(mat->row[row][col], v);
363
43
  return mat;
364
43
}
365
366
/* Replace the element at row "row", column "col" of "mat" by "v".
367
 */
368
__isl_give isl_mat *isl_mat_set_element_val(__isl_take isl_mat *mat,
369
  int row, int col, __isl_take isl_val *v)
370
0
{
371
0
  if (!v)
372
0
    return isl_mat_free(mat);
373
0
  if (!isl_val_is_int(v))
374
0
    isl_die(isl_val_get_ctx(v), isl_error_invalid,
375
0
      "expecting integer value", goto error);
376
0
  mat = isl_mat_set_element(mat, row, col, v->n);
377
0
  isl_val_free(v);
378
0
  return mat;
379
0
error:
380
0
  isl_val_free(v);
381
0
  return isl_mat_free(mat);
382
0
}
383
384
__isl_give isl_mat *isl_mat_diag(isl_ctx *ctx, unsigned n_row, isl_int d)
385
986k
{
386
986k
  int i;
387
986k
  struct isl_mat *mat;
388
986k
389
986k
  mat = isl_mat_alloc(ctx, n_row, n_row);
390
986k
  if (!mat)
391
0
    return NULL;
392
5.41M
  
for (i = 0; 986k
i < n_row;
++i4.43M
) {
393
4.43M
    isl_seq_clr(mat->row[i], i);
394
4.43M
    isl_int_set(mat->row[i][i], d);
395
4.43M
    isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
396
4.43M
  }
397
986k
398
986k
  return mat;
399
986k
}
400
401
/* Create an "n_row" by "n_col" matrix with zero elements.
402
 */
403
__isl_give isl_mat *isl_mat_zero(isl_ctx *ctx, unsigned n_row, unsigned n_col)
404
149
{
405
149
  int i;
406
149
  isl_mat *mat;
407
149
408
149
  mat = isl_mat_alloc(ctx, n_row, n_col);
409
149
  if (!mat)
410
0
    return NULL;
411
190
  
for (i = 0; 149
i < n_row;
++i41
)
412
41
    isl_seq_clr(mat->row[i], n_col);
413
149
414
149
  return mat;
415
149
}
416
417
__isl_give isl_mat *isl_mat_identity(isl_ctx *ctx, unsigned n_row)
418
986k
{
419
986k
  if (!ctx)
420
0
    return NULL;
421
986k
  return isl_mat_diag(ctx, n_row, ctx->one);
422
986k
}
423
424
/* Is "mat" a (possibly scaled) identity matrix?
425
 */
426
int isl_mat_is_scaled_identity(__isl_keep isl_mat *mat)
427
0
{
428
0
  int i;
429
0
430
0
  if (!mat)
431
0
    return -1;
432
0
  if (mat->n_row != mat->n_col)
433
0
    return 0;
434
0
435
0
  for (i = 0; i < mat->n_row; ++i) {
436
0
    if (isl_seq_first_non_zero(mat->row[i], i) != -1)
437
0
      return 0;
438
0
    if (isl_int_ne(mat->row[0][0], mat->row[i][i]))
439
0
      return 0;
440
0
    if (isl_seq_first_non_zero(mat->row[i] + i + 1,
441
0
              mat->n_col - (i + 1)) != -1)
442
0
      return 0;
443
0
  }
444
0
445
0
  return 1;
446
0
}
447
448
__isl_give isl_vec *isl_mat_vec_product(__isl_take isl_mat *mat,
449
  __isl_take isl_vec *vec)
450
328k
{
451
328k
  int i;
452
328k
  struct isl_vec *prod;
453
328k
454
328k
  if (!mat || !vec)
455
0
    goto error;
456
328k
457
328k
  isl_assert(mat->ctx, mat->n_col == vec->size, goto error);
458
328k
459
328k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
460
328k
  if (!prod)
461
0
    goto error;
462
328k
463
2.55M
  
for (i = 0; 328k
i < prod->size;
++i2.22M
)
464
2.22M
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
465
2.22M
          &prod->block.data[i]);
466
328k
  isl_mat_free(mat);
467
328k
  isl_vec_free(vec);
468
328k
  return prod;
469
0
error:
470
0
  isl_mat_free(mat);
471
0
  isl_vec_free(vec);
472
0
  return NULL;
473
328k
}
474
475
__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
476
  __isl_take isl_vec *vec)
477
188
{
478
188
  struct isl_mat *vec_mat;
479
188
  int i;
480
188
481
188
  if (!mat || !vec)
482
0
    goto error;
483
188
  vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
484
188
  if (!vec_mat)
485
0
    goto error;
486
1.51k
  
for (i = 0; 188
i < vec->size;
++i1.32k
)
487
1.32k
    isl_int_set(vec_mat->row[i][0], vec->el[i]);
488
188
  vec_mat = isl_mat_inverse_product(mat, vec_mat);
489
188
  isl_vec_free(vec);
490
188
  if (!vec_mat)
491
0
    return NULL;
492
188
  vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
493
188
  if (vec)
494
1.51k
    
for (i = 0; 188
i < vec->size;
++i1.32k
)
495
1.32k
      isl_int_set(vec->el[i], vec_mat->row[i][0]);
496
188
  isl_mat_free(vec_mat);
497
188
  return vec;
498
0
error:
499
0
  isl_mat_free(mat);
500
0
  isl_vec_free(vec);
501
0
  return NULL;
502
188
}
503
504
__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,
505
  __isl_take isl_mat *mat)
506
21.8k
{
507
21.8k
  int i, j;
508
21.8k
  struct isl_vec *prod;
509
21.8k
510
21.8k
  if (!mat || !vec)
511
0
    goto error;
512
21.8k
513
21.8k
  isl_assert(mat->ctx, mat->n_row == vec->size, goto error);
514
21.8k
515
21.8k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
516
21.8k
  if (!prod)
517
0
    goto error;
518
21.8k
519
122k
  
for (i = 0; 21.8k
i < prod->size;
++i100k
) {
520
100k
    isl_int_set_si(prod->el[i], 0);
521
772k
    for (j = 0; j < vec->size; 
++j671k
)
522
671k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
523
100k
  }
524
21.8k
  isl_mat_free(mat);
525
21.8k
  isl_vec_free(vec);
526
21.8k
  return prod;
527
0
error:
528
0
  isl_mat_free(mat);
529
0
  isl_vec_free(vec);
530
0
  return NULL;
531
21.8k
}
532
533
__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
534
  __isl_take isl_mat *right)
535
147k
{
536
147k
  int i;
537
147k
  struct isl_mat *sum;
538
147k
539
147k
  if (!left || !right)
540
0
    goto error;
541
147k
542
147k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
543
147k
  isl_assert(left->ctx, left->n_row >= 1, goto error);
544
147k
  isl_assert(left->ctx, left->n_col >= 1, goto error);
545
147k
  isl_assert(left->ctx, right->n_col >= 1, goto error);
546
147k
  isl_assert(left->ctx,
547
147k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
548
147k
      goto error);
549
147k
  isl_assert(left->ctx,
550
147k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
551
147k
      goto error);
552
147k
553
147k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
554
147k
  if (!sum)
555
0
    goto error;
556
147k
  isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
557
147k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
558
147k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
559
147k
560
147k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
561
1.12M
  for (i = 1; i < sum->n_row; 
++i976k
) {
562
976k
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
563
976k
    isl_int_addmul(sum->row[i][0],
564
976k
        right->row[0][0], right->row[i][0]);
565
976k
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
566
976k
        left->n_col-1);
567
976k
    isl_seq_scale(sum->row[i]+left->n_col,
568
976k
        right->row[i]+1, right->row[0][0],
569
976k
        right->n_col-1);
570
976k
  }
571
147k
572
147k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
573
147k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
574
147k
  isl_mat_free(left);
575
147k
  isl_mat_free(right);
576
147k
  return sum;
577
0
error:
578
0
  isl_mat_free(left);
579
0
  isl_mat_free(right);
580
0
  return NULL;
581
147k
}
582
583
static void exchange(struct isl_mat *M, struct isl_mat **U,
584
  struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
585
890k
{
586
890k
  int r;
587
6.72M
  for (r = row; r < M->n_row; 
++r5.83M
)
588
5.83M
    isl_int_swap(M->row[r][i], M->row[r][j]);
589
890k
  if (U) {
590
8.70M
    for (r = 0; r < (*U)->n_row; 
++r7.85M
)
591
7.85M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
592
851k
  }
593
890k
  if (Q)
594
327k
    isl_mat_swap_rows(*Q, i, j);
595
890k
}
596
597
static void subtract(struct isl_mat *M, struct isl_mat **U,
598
  struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
599
717k
{
600
717k
  int r;
601
2.70M
  for (r = row; r < M->n_row; 
++r1.98M
)
602
1.98M
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
603
717k
  if (U) {
604
4.43M
    for (r = 0; r < (*U)->n_row; 
++r3.77M
)
605
3.77M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
606
660k
  }
607
717k
  if (Q) {
608
1.22M
    for (r = 0; r < (*Q)->n_col; 
++r1.04M
)
609
1.04M
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
610
186k
  }
611
717k
}
612
613
static void oppose(struct isl_mat *M, struct isl_mat **U,
614
  struct isl_mat **Q, unsigned row, unsigned col)
615
278k
{
616
278k
  int r;
617
1.18M
  for (r = row; r < M->n_row; 
++r905k
)
618
905k
    isl_int_neg(M->row[r][col], M->row[r][col]);
619
278k
  if (U) {
620
1.51M
    for (r = 0; r < (*U)->n_row; 
++r1.28M
)
621
1.28M
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
622
233k
  }
623
278k
  if (Q)
624
134k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
625
278k
}
626
627
/* Given matrix M, compute
628
 *
629
 *    M U = H
630
 *    M   = H Q
631
 *
632
 * with U and Q unimodular matrices and H a matrix in column echelon form
633
 * such that on each echelon row the entries in the non-echelon column
634
 * are non-negative (if neg == 0) or non-positive (if neg == 1)
635
 * and strictly smaller (in absolute value) than the entries in the echelon
636
 * column.
637
 * If U or Q are NULL, then these matrices are not computed.
638
 */
639
__isl_give isl_mat *isl_mat_left_hermite(__isl_take isl_mat *M, int neg,
640
  __isl_give isl_mat **U, __isl_give isl_mat **Q)
641
479k
{
642
479k
  isl_int c;
643
479k
  int row, col;
644
479k
645
479k
  if (U)
646
433k
    *U = NULL;
647
479k
  if (Q)
648
231k
    *Q = NULL;
649
479k
  if (!M)
650
0
    goto error;
651
479k
  M = isl_mat_cow(M);
652
479k
  if (!M)
653
0
    goto error;
654
479k
  if (U) {
655
433k
    *U = isl_mat_identity(M->ctx, M->n_col);
656
433k
    if (!*U)
657
0
      goto error;
658
479k
  }
659
479k
  if (Q) {
660
231k
    *Q = isl_mat_identity(M->ctx, M->n_col);
661
231k
    if (!*Q)
662
0
      goto error;
663
479k
  }
664
479k
665
479k
  col = 0;
666
479k
  isl_int_init(c);
667
2.16M
  for (row = 0; row < M->n_row; 
++row1.68M
) {
668
1.68M
    int first, i, off;
669
1.68M
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
670
1.68M
    if (first == -1)
671
460k
      continue;
672
1.22M
    first += col;
673
1.22M
    if (first != col)
674
639k
      exchange(M, U, Q, row, first, col);
675
1.22M
    if (isl_int_is_neg(M->row[row][col]))
676
1.22M
      
oppose(M, U, Q, row, col)232k
;
677
1.22M
    first = col+1;
678
1.58M
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
679
1.58M
                   M->n_col-first)) != -1) {
680
360k
      first += off;
681
360k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
682
360k
      subtract(M, U, Q, row, col, first, c);
683
360k
      if (!isl_int_is_zero(M->row[row][first]))
684
360k
        
exchange(M, U, Q, row, first, col)22.6k
;
685
338k
      else
686
338k
        ++first;
687
360k
    }
688
5.08M
    for (i = 0; i < col; 
++i3.86M
) {
689
3.86M
      if (isl_int_is_zero(M->row[row][i]))
690
3.86M
        
continue3.73M
;
691
135k
      if (neg)
692
135k
        
isl_int_cdiv_q0
(c, M->row[row][i], M->row[row][col]);
693
135k
      else
694
135k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
695
135k
      if (isl_int_is_zero(c))
696
135k
        
continue42.4k
;
697
92.9k
      subtract(M, U, Q, row, col, i, c);
698
92.9k
    }
699
1.22M
    ++col;
700
1.22M
  }
701
479k
  isl_int_clear(c);
702
479k
703
479k
  return M;
704
0
error:
705
0
  if (Q) {
706
0
    isl_mat_free(*Q);
707
0
    *Q = NULL;
708
0
  }
709
0
  if (U) {
710
0
    isl_mat_free(*U);
711
0
    *U = NULL;
712
0
  }
713
0
  isl_mat_free(M);
714
0
  return NULL;
715
479k
}
716
717
/* Use row "row" of "mat" to eliminate column "col" from all other rows.
718
 */
719
static __isl_give isl_mat *eliminate(__isl_take isl_mat *mat, int row, int col)
720
1.18k
{
721
1.18k
  int k, nr, nc;
722
1.18k
  isl_ctx *ctx;
723
1.18k
724
1.18k
  if (!mat)
725
0
    return NULL;
726
1.18k
727
1.18k
  ctx = isl_mat_get_ctx(mat);
728
1.18k
  nr = isl_mat_rows(mat);
729
1.18k
  nc = isl_mat_cols(mat);
730
1.18k
731
3.28k
  for (k = 0; k < nr; 
++k2.09k
) {
732
2.09k
    if (k == row)
733
1.18k
      continue;
734
908
    if (isl_int_is_zero(mat->row[k][col]))
735
908
      continue;
736
0
    mat = isl_mat_cow(mat);
737
0
    if (!mat)
738
0
      return NULL;
739
0
    isl_seq_elim(mat->row[k], mat->row[row], col, nc, NULL);
740
0
    isl_seq_normalize(ctx, mat->row[k], nc);
741
0
  }
742
1.18k
743
1.18k
  return mat;
744
1.18k
}
745
746
/* Perform Gaussian elimination on the rows of "mat", but start
747
 * from the final row and the final column.
748
 * Any zero rows that result from the elimination are removed.
749
 *
750
 * In particular, for each column from last to first,
751
 * look for the last row with a non-zero coefficient in that column,
752
 * move it last (but before other rows moved last in previous steps) and
753
 * use it to eliminate the column from the other rows.
754
 */
755
__isl_give isl_mat *isl_mat_reverse_gauss(__isl_take isl_mat *mat)
756
1.26k
{
757
1.26k
  int k, row, last, nr, nc;
758
1.26k
759
1.26k
  if (!mat)
760
0
    return NULL;
761
1.26k
762
1.26k
  nr = isl_mat_rows(mat);
763
1.26k
  nc = isl_mat_cols(mat);
764
1.26k
765
1.26k
  last = nc - 1;
766
2.45k
  for (row = nr - 1; row >= 0; 
--row1.18k
) {
767
1.57k
    for (; last >= 0; 
--last384
) {
768
2.03k
      for (k = row; k >= 0; 
--k462
)
769
1.65k
        if (!isl_int_is_zero(mat->row[k][last]))
770
1.65k
          
break1.18k
;
771
1.57k
      if (k >= 0)
772
1.18k
        break;
773
1.57k
    }
774
1.18k
    if (last < 0)
775
0
      break;
776
1.18k
    if (k != row)
777
0
      mat = isl_mat_swap_rows(mat, k, row);
778
1.18k
    if (!mat)
779
0
      return NULL;
780
1.18k
    if (isl_int_is_neg(mat->row[row][last]))
781
1.18k
      
mat = isl_mat_row_neg(mat, row)4
;
782
1.18k
    mat = eliminate(mat, row, last);
783
1.18k
    if (!mat)
784
0
      return NULL;
785
1.18k
  }
786
1.26k
  mat = isl_mat_drop_rows(mat, 0, row + 1);
787
1.26k
788
1.26k
  return mat;
789
1.26k
}
790
791
/* Negate the lexicographically negative rows of "mat" such that
792
 * all rows in the result are lexicographically non-negative.
793
 */
794
__isl_give isl_mat *isl_mat_lexnonneg_rows(__isl_take isl_mat *mat)
795
1.26k
{
796
1.26k
  int i, nr, nc;
797
1.26k
798
1.26k
  if (!mat)
799
0
    return NULL;
800
1.26k
801
1.26k
  nr = isl_mat_rows(mat);
802
1.26k
  nc = isl_mat_cols(mat);
803
1.26k
804
2.45k
  for (i = 0; i < nr; 
++i1.18k
) {
805
1.18k
    int pos;
806
1.18k
807
1.18k
    pos = isl_seq_first_non_zero(mat->row[i], nc);
808
1.18k
    if (pos < 0)
809
0
      continue;
810
1.18k
    if (isl_int_is_nonneg(mat->row[i][pos]))
811
1.18k
      
continue1.17k
;
812
12
    mat = isl_mat_row_neg(mat, i);
813
12
    if (!mat)
814
0
      return NULL;
815
12
  }
816
1.26k
817
1.26k
  return mat;
818
1.26k
}
819
820
/* Given a matrix "H" is column echelon form, what is the first
821
 * zero column?  That is how many initial columns are non-zero?
822
 * Start looking at column "first_col" and only consider
823
 * the columns to be of size "n_row".
824
 * "H" is assumed to be non-NULL.
825
 *
826
 * Since "H" is in column echelon form, the first non-zero entry
827
 * in a column is always in a later position compared to the previous column.
828
 */
829
static int hermite_first_zero_col(__isl_keep isl_mat *H, int first_col,
830
  int n_row)
831
68
{
832
68
  int row, col;
833
68
834
215
  for (col = first_col, row = 0; col < H->n_col; 
++col147
) {
835
316
    for (; row < n_row; 
++row127
)
836
274
      if (!isl_int_is_zero(H->row[row][col]))
837
274
        
break147
;
838
189
    if (row == n_row)
839
42
      return col;
840
189
  }
841
68
842
68
  
return H->n_col26
;
843
68
}
844
845
/* Return the rank of "mat", or -1 in case of error.
846
 */
847
int isl_mat_rank(__isl_keep isl_mat *mat)
848
68
{
849
68
  int rank;
850
68
  isl_mat *H;
851
68
852
68
  H = isl_mat_left_hermite(isl_mat_copy(mat), 0, NULL, NULL);
853
68
  if (!H)
854
0
    return -1;
855
68
856
68
  rank = hermite_first_zero_col(H, 0, H->n_row);
857
68
  isl_mat_free(H);
858
68
859
68
  return rank;
860
68
}
861
862
__isl_give isl_mat *isl_mat_right_kernel(__isl_take isl_mat *mat)
863
0
{
864
0
  int rank;
865
0
  struct isl_mat *U = NULL;
866
0
  struct isl_mat *K;
867
0
868
0
  mat = isl_mat_left_hermite(mat, 0, &U, NULL);
869
0
  if (!mat || !U)
870
0
    goto error;
871
0
872
0
  rank = hermite_first_zero_col(mat, 0, mat->n_row);
873
0
  K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);
874
0
  if (!K)
875
0
    goto error;
876
0
  isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
877
0
  isl_mat_free(mat);
878
0
  isl_mat_free(U);
879
0
  return K;
880
0
error:
881
0
  isl_mat_free(mat);
882
0
  isl_mat_free(U);
883
0
  return NULL;
884
0
}
885
886
__isl_give isl_mat *isl_mat_lin_to_aff(__isl_take isl_mat *mat)
887
697k
{
888
697k
  int i;
889
697k
  struct isl_mat *mat2;
890
697k
891
697k
  if (!mat)
892
0
    return NULL;
893
697k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
894
697k
  if (!mat2)
895
0
    goto error;
896
697k
  isl_int_set_si(mat2->row[0][0], 1);
897
697k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
898
4.27M
  for (i = 0; i < mat->n_row; 
++i3.57M
) {
899
3.57M
    isl_int_set_si(mat2->row[1+i][0], 0);
900
3.57M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
901
3.57M
  }
902
697k
  isl_mat_free(mat);
903
697k
  return mat2;
904
0
error:
905
0
  isl_mat_free(mat);
906
0
  return NULL;
907
697k
}
908
909
/* Given two matrices M1 and M2, return the block matrix
910
 *
911
 *  [ M1  0  ]
912
 *  [ 0   M2 ]
913
 */
914
__isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1,
915
  __isl_take isl_mat *mat2)
916
132k
{
917
132k
  int i;
918
132k
  isl_mat *mat;
919
132k
920
132k
  if (!mat1 || !mat2)
921
0
    goto error;
922
132k
923
132k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
924
132k
               mat1->n_col + mat2->n_col);
925
132k
  if (!mat)
926
0
    goto error;
927
264k
  
for (i = 0; 132k
i < mat1->n_row;
++i132k
) {
928
132k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
929
132k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
930
132k
  }
931
845k
  for (i = 0; i < mat2->n_row; 
++i713k
) {
932
713k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
933
713k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
934
713k
                mat2->row[i], mat2->n_col);
935
713k
  }
936
132k
  isl_mat_free(mat1);
937
132k
  isl_mat_free(mat2);
938
132k
  return mat;
939
0
error:
940
0
  isl_mat_free(mat1);
941
0
  isl_mat_free(mat2);
942
0
  return NULL;
943
132k
}
944
945
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
946
1.00M
{
947
1.00M
  int i;
948
1.00M
949
1.71M
  for (i = 0; i < n_row; 
++i703k
)
950
1.20M
    if (!isl_int_is_zero(row[i][col]))
951
1.20M
      
return i505k
;
952
1.00M
  
return -1501k
;
953
1.00M
}
954
955
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
956
500k
{
957
500k
  int i, min = row_first_non_zero(row, n_row, col);
958
500k
  if (min < 0)
959
0
    return -1;
960
1.20M
  
for (i = min + 1; 500k
i < n_row;
++i704k
) {
961
704k
    if (isl_int_is_zero(row[i][col]))
962
704k
      
continue700k
;
963
4.00k
    if (isl_int_abs_lt(row[i][col], row[min][col]))
964
4.00k
      
min = i63
;
965
4.00k
  }
966
500k
  return min;
967
500k
}
968
969
static isl_stat inv_exchange(__isl_keep isl_mat **left,
970
  __isl_keep isl_mat **right, unsigned i, unsigned j)
971
440
{
972
440
  *left = isl_mat_swap_rows(*left, i, j);
973
440
  *right = isl_mat_swap_rows(*right, i, j);
974
440
975
440
  if (!*left || !*right)
976
0
    return isl_stat_error;
977
440
  return isl_stat_ok;
978
440
}
979
980
static void inv_oppose(
981
  struct isl_mat *left, struct isl_mat *right, unsigned row)
982
133
{
983
133
  isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
984
133
  isl_seq_neg(right->row[row], right->row[row], right->n_col);
985
133
}
986
987
static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
988
  unsigned row, unsigned i, isl_int m)
989
4.01k
{
990
4.01k
  isl_int_neg(m, m);
991
4.01k
  isl_seq_combine(left->row[i]+row,
992
4.01k
      left->ctx->one, left->row[i]+row,
993
4.01k
      m, left->row[row]+row,
994
4.01k
      left->n_col-row);
995
4.01k
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
996
4.01k
      m, right->row[row], right->n_col);
997
4.01k
}
998
999
/* Compute inv(left)*right
1000
 */
1001
__isl_give isl_mat *isl_mat_inverse_product(__isl_take isl_mat *left,
1002
  __isl_take isl_mat *right)
1003
167k
{
1004
167k
  int row;
1005
167k
  isl_int a, b;
1006
167k
1007
167k
  if (!left || !right)
1008
0
    goto error;
1009
167k
1010
167k
  isl_assert(left->ctx, left->n_row == left->n_col, goto error);
1011
167k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
1012
167k
1013
167k
  if (left->n_row == 0) {
1014
0
    isl_mat_free(left);
1015
0
    return right;
1016
0
  }
1017
167k
1018
167k
  left = isl_mat_cow(left);
1019
167k
  right = isl_mat_cow(right);
1020
167k
  if (!left || !right)
1021
0
    goto error;
1022
167k
1023
167k
  isl_int_init(a);
1024
167k
  isl_int_init(b);
1025
668k
  for (row = 0; row < left->n_row; 
++row500k
) {
1026
500k
    int pivot, first, i, off;
1027
500k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
1028
500k
    if (pivot < 0) {
1029
0
      isl_int_clear(a);
1030
0
      isl_int_clear(b);
1031
0
      isl_assert(left->ctx, pivot >= 0, goto error);
1032
0
    }
1033
500k
    pivot += row;
1034
500k
    if (pivot != row)
1035
438
      if (inv_exchange(&left, &right, pivot, row) < 0)
1036
0
        goto error;
1037
500k
    if (isl_int_is_neg(left->row[row][row]))
1038
500k
      
inv_oppose(left, right, row)133
;
1039
500k
    first = row+1;
1040
504k
    while ((off = row_first_non_zero(left->row+first,
1041
504k
          left->n_row-first, row)) != -1) {
1042
4.01k
      first += off;
1043
4.01k
      isl_int_fdiv_q(a, left->row[first][row],
1044
4.01k
          left->row[row][row]);
1045
4.01k
      inv_subtract(left, right, row, first, a);
1046
4.01k
      if (!isl_int_is_zero(left->row[first][row])) {
1047
2
        if (inv_exchange(&left, &right, row, first) < 0)
1048
0
          goto error;
1049
4.00k
      } else {
1050
4.00k
        ++first;
1051
4.00k
      }
1052
4.01k
    }
1053
1.20M
    
for (i = 0; 500k
i < row;
++i705k
) {
1054
705k
      if (isl_int_is_zero(left->row[i][row]))
1055
705k
        
continue705k
;
1056
497
      isl_int_gcd(a, left->row[row][row], left->row[i][row]);
1057
497
      isl_int_divexact(b, left->row[i][row], a);
1058
497
      isl_int_divexact(a, left->row[row][row], a);
1059
497
      isl_int_neg(b, b);
1060
497
      isl_seq_combine(left->row[i] + i,
1061
497
          a, left->row[i] + i,
1062
497
          b, left->row[row] + i,
1063
497
          left->n_col - i);
1064
497
      isl_seq_combine(right->row[i], a, right->row[i],
1065
497
          b, right->row[row], right->n_col);
1066
497
    }
1067
500k
  }
1068
167k
  isl_int_clear(b);
1069
167k
1070
167k
  isl_int_set(a, left->row[0][0]);
1071
500k
  for (row = 1; row < left->n_row; 
++row333k
)
1072
333k
    isl_int_lcm(a, a, left->row[row][row]);
1073
167k
  if (isl_int_is_zero(a)){
1074
0
    isl_int_clear(a);
1075
0
    isl_assert(left->ctx, 0, goto error);
1076
0
  }
1077
668k
  
for (row = 0; 167k
row < left->n_row;
++row500k
) {
1078
500k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
1079
500k
    if (isl_int_is_one(left->row[row][row]))
1080
500k
      
continue492k
;
1081
8.29k
    isl_seq_scale(right->row[row], right->row[row],
1082
8.29k
        left->row[row][row], right->n_col);
1083
8.29k
  }
1084
167k
  isl_int_clear(a);
1085
167k
1086
167k
  isl_mat_free(left);
1087
167k
  return right;
1088
0
error:
1089
0
  isl_mat_free(left);
1090
0
  isl_mat_free(right);
1091
0
  return NULL;
1092
167k
}
1093
1094
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
1095
59.7k
{
1096
59.7k
  int i;
1097
59.7k
1098
296k
  for (i = 0; i < mat->n_row; 
++i236k
)
1099
236k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
1100
59.7k
}
1101
1102
void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
1103
  isl_int m1, unsigned src1, isl_int m2, unsigned src2)
1104
168k
{
1105
168k
  int i;
1106
168k
  isl_int tmp;
1107
168k
1108
168k
  isl_int_init(tmp);
1109
782k
  for (i = 0; i < mat->n_row; 
++i613k
) {
1110
613k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
1111
613k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
1112
613k
    isl_int_set(mat->row[i][dst], tmp);
1113
613k
  }
1114
168k
  isl_int_clear(tmp);
1115
168k
}
1116
1117
__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
1118
65.2k
{
1119
65.2k
  struct isl_mat *inv;
1120
65.2k
  int row;
1121
65.2k
  isl_int a, b;
1122
65.2k
1123
65.2k
  mat = isl_mat_cow(mat);
1124
65.2k
  if (!mat)
1125
0
    return NULL;
1126
65.2k
1127
65.2k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
1128
65.2k
  inv = isl_mat_cow(inv);
1129
65.2k
  if (!inv)
1130
0
    goto error;
1131
65.2k
1132
65.2k
  isl_int_init(a);
1133
65.2k
  isl_int_init(b);
1134
277k
  for (row = 0; row < mat->n_row; 
++row211k
) {
1135
211k
    int pivot, first, i, off;
1136
211k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
1137
211k
    if (pivot < 0) {
1138
0
      isl_int_clear(a);
1139
0
      isl_int_clear(b);
1140
0
      isl_assert(mat->ctx, pivot >= 0, goto error);
1141
0
    }
1142
211k
    pivot += row;
1143
211k
    if (pivot != row)
1144
32.3k
      exchange(mat, &inv, NULL, row, pivot, row);
1145
211k
    if (isl_int_is_neg(mat->row[row][row]))
1146
211k
      
oppose(mat, &inv, NULL, row, row)45.6k
;
1147
211k
    first = row+1;
1148
475k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
1149
475k
                mat->n_col-first)) != -1) {
1150
263k
      first += off;
1151
263k
      isl_int_fdiv_q(a, mat->row[row][first],
1152
263k
                mat->row[row][row]);
1153
263k
      subtract(mat, &inv, NULL, row, row, first, a);
1154
263k
      if (!isl_int_is_zero(mat->row[row][first]))
1155
263k
        
exchange(mat, &inv, NULL, row, row, first)196k
;
1156
67.4k
      else
1157
67.4k
        ++first;
1158
263k
    }
1159
556k
    for (i = 0; i < row; 
++i345k
) {
1160
345k
      if (isl_int_is_zero(mat->row[row][i]))
1161
345k
        
continue260k
;
1162
84.4k
      isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
1163
84.4k
      isl_int_divexact(b, mat->row[row][i], a);
1164
84.4k
      isl_int_divexact(a, mat->row[row][row], a);
1165
84.4k
      isl_int_neg(a, a);
1166
84.4k
      isl_mat_col_combine(mat, i, a, i, b, row);
1167
84.4k
      isl_mat_col_combine(inv, i, a, i, b, row);
1168
84.4k
    }
1169
211k
  }
1170
65.2k
  isl_int_clear(b);
1171
65.2k
1172
65.2k
  isl_int_set(a, mat->row[0][0]);
1173
211k
  for (row = 1; row < mat->n_row; 
++row146k
)
1174
146k
    isl_int_lcm(a, a, mat->row[row][row]);
1175
65.2k
  if (isl_int_is_zero(a)){
1176
0
    isl_int_clear(a);
1177
0
    goto error;
1178
0
  }
1179
277k
  
for (row = 0; 65.2k
row < mat->n_row;
++row211k
) {
1180
211k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1181
211k
    if (isl_int_is_one(mat->row[row][row]))
1182
211k
      
continue152k
;
1183
59.7k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1184
59.7k
  }
1185
65.2k
  isl_int_clear(a);
1186
65.2k
1187
65.2k
  isl_mat_free(mat);
1188
65.2k
1189
65.2k
  return inv;
1190
0
error:
1191
0
  isl_mat_free(mat);
1192
0
  isl_mat_free(inv);
1193
0
  return NULL;
1194
65.2k
}
1195
1196
__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1197
4.53k
{
1198
4.53k
  struct isl_mat *transpose = NULL;
1199
4.53k
  int i, j;
1200
4.53k
1201
4.53k
  if (!mat)
1202
0
    return NULL;
1203
4.53k
1204
4.53k
  if (mat->n_col == mat->n_row) {
1205
4.53k
    mat = isl_mat_cow(mat);
1206
4.53k
    if (!mat)
1207
0
      return NULL;
1208
26.0k
    
for (i = 0; 4.53k
i < mat->n_row;
++i21.5k
)
1209
80.9k
      
for (j = i + 1; 21.5k
j < mat->n_col;
++j59.4k
)
1210
59.4k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1211
4.53k
    return mat;
1212
4.53k
  }
1213
0
  transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
1214
0
  if (!transpose)
1215
0
    goto error;
1216
0
  for (i = 0; i < mat->n_row; ++i)
1217
0
    for (j = 0; j < mat->n_col; ++j)
1218
0
      isl_int_set(transpose->row[j][i], mat->row[i][j]);
1219
0
  isl_mat_free(mat);
1220
0
  return transpose;
1221
0
error:
1222
0
  isl_mat_free(mat);
1223
0
  return NULL;
1224
0
}
1225
1226
__isl_give isl_mat *isl_mat_swap_cols(__isl_take isl_mat *mat,
1227
  unsigned i, unsigned j)
1228
703k
{
1229
703k
  int r;
1230
703k
1231
703k
  mat = isl_mat_cow(mat);
1232
703k
  if (check_col_range(mat, i, 1) < 0 ||
1233
703k
      check_col_range(mat, j, 1) < 0)
1234
0
    return isl_mat_free(mat);
1235
703k
1236
20.1M
  
for (r = 0; 703k
r < mat->n_row;
++r19.4M
)
1237
19.4M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1238
703k
  return mat;
1239
703k
}
1240
1241
__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1242
  unsigned i, unsigned j)
1243
1.99M
{
1244
1.99M
  isl_int *t;
1245
1.99M
1246
1.99M
  if (!mat)
1247
0
    return NULL;
1248
1.99M
  mat = isl_mat_cow(mat);
1249
1.99M
  if (check_row_range(mat, i, 1) < 0 ||
1250
1.99M
      check_row_range(mat, j, 1) < 0)
1251
0
    return isl_mat_free(mat);
1252
1.99M
1253
1.99M
  t = mat->row[i];
1254
1.99M
  mat->row[i] = mat->row[j];
1255
1.99M
  mat->row[j] = t;
1256
1.99M
  return mat;
1257
1.99M
}
1258
1259
/* Calculate the product of two matrices.
1260
 *
1261
 * This function is optimized for operand matrices that contain many zeros and
1262
 * skips multiplications where we know one of the operands is zero.
1263
 */
1264
__isl_give isl_mat *isl_mat_product(__isl_take isl_mat *left,
1265
  __isl_take isl_mat *right)
1266
1.98M
{
1267
1.98M
  int i, j, k;
1268
1.98M
  struct isl_mat *prod;
1269
1.98M
1270
1.98M
  if (!left || !right)
1271
0
    goto error;
1272
1.98M
  isl_assert(left->ctx, left->n_col == right->n_row, goto error);
1273
1.98M
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1274
1.98M
  if (!prod)
1275
0
    goto error;
1276
1.98M
  if (left->n_col == 0) {
1277
51
    for (i = 0; i < prod->n_row; 
++i12
)
1278
12
      isl_seq_clr(prod->row[i], prod->n_col);
1279
39
    isl_mat_free(left);
1280
39
    isl_mat_free(right);
1281
39
    return prod;
1282
39
  }
1283
6.54M
  
for (i = 0; 1.98M
i < prod->n_row;
++i4.56M
) {
1284
38.0M
    for (j = 0; j < prod->n_col; 
++j33.4M
)
1285
33.4M
      isl_int_mul(prod->row[i][j],
1286
4.56M
            left->row[i][0], right->row[0][j]);
1287
43.8M
    for (k = 1; k < left->n_col; 
++k39.3M
) {
1288
39.3M
      if (isl_int_is_zero(left->row[i][k]))
1289
39.3M
        
continue34.4M
;
1290
45.7M
      
for (j = 0; 4.84M
j < prod->n_col;
++j40.9M
)
1291
40.9M
        isl_int_addmul(prod->row[i][j],
1292
4.84M
              left->row[i][k], right->row[k][j]);
1293
4.84M
    }
1294
4.56M
  }
1295
1.98M
  isl_mat_free(left);
1296
1.98M
  isl_mat_free(right);
1297
1.98M
  return prod;
1298
0
error:
1299
0
  isl_mat_free(left);
1300
0
  isl_mat_free(right);
1301
0
  return NULL;
1302
1.98M
}
1303
1304
/* Replace the variables x in the rows q by x' given by x = M x',
1305
 * with M the matrix mat.
1306
 *
1307
 * If the number of new variables is greater than the original
1308
 * number of variables, then the rows q have already been
1309
 * preextended.  If the new number is smaller, then the coefficients
1310
 * of the divs, which are not changed, need to be shifted down.
1311
 * The row q may be the equalities, the inequalities or the
1312
 * div expressions.  In the latter case, has_div is true and
1313
 * we need to take into account the extra denominator column.
1314
 */
1315
static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
1316
  unsigned n_div, int has_div, struct isl_mat *mat)
1317
1.66M
{
1318
1.66M
  int i;
1319
1.66M
  struct isl_mat *t;
1320
1.66M
  int e;
1321
1.66M
1322
1.66M
  if (mat->n_col >= mat->n_row)
1323
1.03M
    e = 0;
1324
632k
  else
1325
632k
    e = mat->n_row - mat->n_col;
1326
1.66M
  if (has_div)
1327
554k
    
for (i = 0; 554k
i < n;
++i137
)
1328
554k
      
isl_int_mul137
(q[i][0], q[i][0], mat->row[0][0]);
1329
1.66M
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1330
1.66M
  t = isl_mat_product(t, mat);
1331
1.66M
  if (!t)
1332
0
    return -1;
1333
4.47M
  
for (i = 0; 1.66M
i < n;
++i2.80M
) {
1334
2.80M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1335
2.80M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1336
2.80M
          q[i] + has_div + t->n_col + e, n_div);
1337
2.80M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1338
2.80M
  }
1339
1.66M
  isl_mat_free(t);
1340
1.66M
  return 0;
1341
1.66M
}
1342
1343
/* Replace the variables x in bset by x' given by x = M x', with
1344
 * M the matrix mat.
1345
 *
1346
 * If there are fewer variables x' then there are x, then we perform
1347
 * the transformation in place, which means that, in principle,
1348
 * this frees up some extra variables as the number
1349
 * of columns remains constant, but we would have to extend
1350
 * the div array too as the number of rows in this array is assumed
1351
 * to be equal to extra.
1352
 */
1353
__isl_give isl_basic_set *isl_basic_set_preimage(
1354
  __isl_take isl_basic_set *bset, __isl_take isl_mat *mat)
1355
554k
{
1356
554k
  struct isl_ctx *ctx;
1357
554k
1358
554k
  if (!bset || !mat)
1359
0
    goto error;
1360
554k
1361
554k
  ctx = bset->ctx;
1362
554k
  bset = isl_basic_set_cow(bset);
1363
554k
  if (!bset)
1364
0
    goto error;
1365
554k
1366
554k
  isl_assert(ctx, bset->dim->nparam == 0, goto error);
1367
554k
  isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1368
554k
  isl_assert(ctx, mat->n_col > 0, goto error);
1369
554k
1370
554k
  if (mat->n_col > mat->n_row) {
1371
100k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1372
100k
    if (!bset)
1373
0
      goto error;
1374
454k
  } else if (mat->n_col < mat->n_row) {
1375
210k
    bset->dim = isl_space_cow(bset->dim);
1376
210k
    if (!bset->dim)
1377
0
      goto error;
1378
210k
    bset->dim->n_out -= mat->n_row - mat->n_col;
1379
210k
  }
1380
554k
1381
554k
  if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1382
554k
      isl_mat_copy(mat)) < 0)
1383
0
    goto error;
1384
554k
1385
554k
  if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1386
554k
      isl_mat_copy(mat)) < 0)
1387
0
    goto error;
1388
554k
1389
554k
  if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
1390
0
    goto error2;
1391
554k
1392
554k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1393
554k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1394
554k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1395
554k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1396
554k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1397
554k
1398
554k
  bset = isl_basic_set_simplify(bset);
1399
554k
  bset = isl_basic_set_finalize(bset);
1400
554k
1401
554k
  return bset;
1402
0
error:
1403
0
  isl_mat_free(mat);
1404
0
error2:
1405
0
  isl_basic_set_free(bset);
1406
0
  return NULL;
1407
0
}
1408
1409
__isl_give isl_set *isl_set_preimage(
1410
  __isl_take isl_set *set, __isl_take isl_mat *mat)
1411
31.9k
{
1412
31.9k
  int i;
1413
31.9k
1414
31.9k
  set = isl_set_cow(set);
1415
31.9k
  if (!set)
1416
0
    goto error;
1417
31.9k
1418
90.8k
  
for (i = 0; 31.9k
i < set->n;
++i58.9k
) {
1419
58.9k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1420
58.9k
                isl_mat_copy(mat));
1421
58.9k
    if (!set->p[i])
1422
0
      goto error;
1423
58.9k
  }
1424
31.9k
  if (mat->n_col != mat->n_row) {
1425
13.6k
    set->dim = isl_space_cow(set->dim);
1426
13.6k
    if (!set->dim)
1427
0
      goto error;
1428
13.6k
    set->dim->n_out += mat->n_col;
1429
13.6k
    set->dim->n_out -= mat->n_row;
1430
13.6k
  }
1431
31.9k
  isl_mat_free(mat);
1432
31.9k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1433
31.9k
  return set;
1434
0
error:
1435
0
  isl_set_free(set);
1436
0
  isl_mat_free(mat);
1437
0
  return NULL;
1438
31.9k
}
1439
1440
/* Replace the variables x starting at "first_col" in the rows "rows"
1441
 * of some coefficient matrix by x' with x = M x' with M the matrix mat.
1442
 * That is, replace the corresponding coefficients c by c M.
1443
 */
1444
isl_stat isl_mat_sub_transform(isl_int **row, unsigned n_row,
1445
  unsigned first_col, __isl_take isl_mat *mat)
1446
222
{
1447
222
  int i;
1448
222
  isl_ctx *ctx;
1449
222
  isl_mat *t;
1450
222
1451
222
  if (!mat)
1452
0
    return isl_stat_error;
1453
222
  ctx = isl_mat_get_ctx(mat);
1454
222
  t = isl_mat_sub_alloc6(ctx, row, 0, n_row, first_col, mat->n_row);
1455
222
  t = isl_mat_product(t, mat);
1456
222
  if (!t)
1457
0
    return isl_stat_error;
1458
338
  
for (i = 0; 222
i < n_row;
++i116
)
1459
116
    isl_seq_swp_or_cpy(row[i] + first_col, t->row[i], t->n_col);
1460
222
  isl_mat_free(t);
1461
222
  return isl_stat_ok;
1462
222
}
1463
1464
void isl_mat_print_internal(__isl_keep isl_mat *mat, FILE *out, int indent)
1465
0
{
1466
0
  int i, j;
1467
0
1468
0
  if (!mat) {
1469
0
    fprintf(out, "%*snull mat\n", indent, "");
1470
0
    return;
1471
0
  }
1472
0
1473
0
  if (mat->n_row == 0)
1474
0
    fprintf(out, "%*s[]\n", indent, "");
1475
0
1476
0
  for (i = 0; i < mat->n_row; ++i) {
1477
0
    if (!i)
1478
0
      fprintf(out, "%*s[[", indent, "");
1479
0
    else
1480
0
      fprintf(out, "%*s[", indent+1, "");
1481
0
    for (j = 0; j < mat->n_col; ++j) {
1482
0
      if (j)
1483
0
          fprintf(out, ",");
1484
0
      isl_int_print(out, mat->row[i][j], 0);
1485
0
    }
1486
0
    if (i == mat->n_row-1)
1487
0
      fprintf(out, "]]\n");
1488
0
    else
1489
0
      fprintf(out, "]\n");
1490
0
  }
1491
0
}
1492
1493
void isl_mat_dump(__isl_keep isl_mat *mat)
1494
0
{
1495
0
  isl_mat_print_internal(mat, stderr, 0);
1496
0
}
1497
1498
__isl_give isl_mat *isl_mat_drop_cols(__isl_take isl_mat *mat,
1499
  unsigned col, unsigned n)
1500
84.6k
{
1501
84.6k
  int r;
1502
84.6k
1503
84.6k
  if (n == 0)
1504
22
    return mat;
1505
84.6k
1506
84.6k
  mat = isl_mat_cow(mat);
1507
84.6k
  if (check_col_range(mat, col, n) < 0)
1508
0
    return isl_mat_free(mat);
1509
84.6k
1510
84.6k
  if (col != mat->n_col-n) {
1511
66.6k
    for (r = 0; r < mat->n_row; 
++r52.2k
)
1512
52.2k
      isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1513
52.2k
          mat->n_col - col - n);
1514
14.4k
  }
1515
84.6k
  mat->n_col -= n;
1516
84.6k
  return mat;
1517
84.6k
}
1518
1519
__isl_give isl_mat *isl_mat_drop_rows(__isl_take isl_mat *mat,
1520
  unsigned row, unsigned n)
1521
120k
{
1522
120k
  int r;
1523
120k
1524
120k
  mat = isl_mat_cow(mat);
1525
120k
  if (check_row_range(mat, row, n) < 0)
1526
0
    return isl_mat_free(mat);
1527
120k
1528
374k
  
for (r = row; 120k
r+n < mat->n_row;
++r254k
)
1529
254k
    mat->row[r] = mat->row[r+n];
1530
120k
1531
120k
  mat->n_row -= n;
1532
120k
  return mat;
1533
120k
}
1534
1535
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1536
        unsigned col, unsigned n)
1537
87.0k
{
1538
87.0k
  isl_mat *ext;
1539
87.0k
1540
87.0k
  if (check_col_range(mat, col, 0) < 0)
1541
0
    return isl_mat_free(mat);
1542
87.0k
  if (n == 0)
1543
448
    return mat;
1544
86.5k
1545
86.5k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1546
86.5k
  if (!ext)
1547
0
    goto error;
1548
86.5k
1549
86.5k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1550
86.5k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1551
86.5k
        col + n, col, mat->n_col - col);
1552
86.5k
1553
86.5k
  isl_mat_free(mat);
1554
86.5k
  return ext;
1555
0
error:
1556
0
  isl_mat_free(mat);
1557
0
  return NULL;
1558
86.5k
}
1559
1560
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1561
  unsigned first, unsigned n)
1562
87.0k
{
1563
87.0k
  int i;
1564
87.0k
1565
87.0k
  if (!mat)
1566
0
    return NULL;
1567
87.0k
  mat = isl_mat_insert_cols(mat, first, n);
1568
87.0k
  if (!mat)
1569
0
    return NULL;
1570
87.0k
1571
91.2k
  
for (i = 0; 87.0k
i < mat->n_row;
++i4.20k
)
1572
4.20k
    isl_seq_clr(mat->row[i] + first, n);
1573
87.0k
1574
87.0k
  return mat;
1575
87.0k
}
1576
1577
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1578
4.66k
{
1579
4.66k
  if (!mat)
1580
0
    return NULL;
1581
4.66k
1582
4.66k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1583
4.66k
}
1584
1585
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1586
        unsigned row, unsigned n)
1587
5.10k
{
1588
5.10k
  isl_mat *ext;
1589
5.10k
1590
5.10k
  if (check_row_range(mat, row, 0) < 0)
1591
0
    return isl_mat_free(mat);
1592
5.10k
  if (n == 0)
1593
444
    return mat;
1594
4.65k
1595
4.65k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1596
4.65k
  if (!ext)
1597
0
    goto error;
1598
4.65k
1599
4.65k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1600
4.65k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1601
4.65k
        mat->n_row - row, 0, 0, mat->n_col);
1602
4.65k
1603
4.65k
  isl_mat_free(mat);
1604
4.65k
  return ext;
1605
0
error:
1606
0
  isl_mat_free(mat);
1607
0
  return NULL;
1608
4.65k
}
1609
1610
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1611
5.09k
{
1612
5.09k
  if (!mat)
1613
0
    return NULL;
1614
5.09k
1615
5.09k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1616
5.09k
}
1617
1618
__isl_give isl_mat *isl_mat_insert_zero_rows(__isl_take isl_mat *mat,
1619
  unsigned row, unsigned n)
1620
2
{
1621
2
  int i;
1622
2
1623
2
  mat = isl_mat_insert_rows(mat, row, n);
1624
2
  if (!mat)
1625
0
    return NULL;
1626
2
  
1627
4
  
for (i = 0; 2
i < n;
++i2
)
1628
2
    isl_seq_clr(mat->row[row + i], mat->n_col);
1629
2
1630
2
  return mat;
1631
2
}
1632
1633
__isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n)
1634
0
{
1635
0
  if (!mat)
1636
0
    return NULL;
1637
0
1638
0
  return isl_mat_insert_zero_rows(mat, mat->n_row, n);
1639
0
}
1640
1641
void isl_mat_col_submul(struct isl_mat *mat,
1642
      int dst_col, isl_int f, int src_col)
1643
7.00k
{
1644
7.00k
  int i;
1645
7.00k
1646
48.4k
  for (i = 0; i < mat->n_row; 
++i41.4k
)
1647
41.4k
    isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1648
7.00k
}
1649
1650
void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1651
2
{
1652
2
  int i;
1653
2
1654
2
  if (!mat)
1655
0
    return;
1656
2
1657
7
  
for (i = 0; 2
i < mat->n_row;
++i5
)
1658
5
    isl_int_add(mat->row[i][dst_col],
1659
2
          mat->row[i][dst_col], mat->row[i][src_col]);
1660
2
}
1661
1662
void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1663
14.4k
{
1664
14.4k
  int i;
1665
14.4k
1666
72.1k
  for (i = 0; i < mat->n_row; 
++i57.7k
)
1667
57.7k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1668
14.4k
}
1669
1670
/* Add "f" times column "src_col" to column "dst_col" of "mat" and
1671
 * return the result.
1672
 */
1673
__isl_give isl_mat *isl_mat_col_addmul(__isl_take isl_mat *mat, int dst_col,
1674
  isl_int f, int src_col)
1675
6
{
1676
6
  int i;
1677
6
1678
6
  if (check_col(mat, dst_col) < 0 || check_col(mat, src_col) < 0)
1679
0
    return isl_mat_free(mat);
1680
6
1681
14
  
for (i = 0; 6
i < mat->n_row;
++i8
) {
1682
8
    if (isl_int_is_zero(mat->row[i][src_col]))
1683
8
      
continue2
;
1684
6
    mat = isl_mat_cow(mat);
1685
6
    if (!mat)
1686
0
      return NULL;
1687
6
    isl_int_addmul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1688
6
  }
1689
6
1690
6
  return mat;
1691
6
}
1692
1693
/* Negate column "col" of "mat" and return the result.
1694
 */
1695
__isl_give isl_mat *isl_mat_col_neg(__isl_take isl_mat *mat, int col)
1696
4
{
1697
4
  int i;
1698
4
1699
4
  if (check_col(mat, col) < 0)
1700
0
    return isl_mat_free(mat);
1701
4
1702
10
  
for (i = 0; 4
i < mat->n_row;
++i6
) {
1703
6
    if (isl_int_is_zero(mat->row[i][col]))
1704
6
      
continue2
;
1705
4
    mat = isl_mat_cow(mat);
1706
4
    if (!mat)
1707
0
      return NULL;
1708
4
    isl_int_neg(mat->row[i][col], mat->row[i][col]);
1709
4
  }
1710
4
1711
4
  return mat;
1712
4
}
1713
1714
/* Negate row "row" of "mat" and return the result.
1715
 */
1716
__isl_give isl_mat *isl_mat_row_neg(__isl_take isl_mat *mat, int row)
1717
16
{
1718
16
  if (check_row(mat, row) < 0)
1719
0
    return isl_mat_free(mat);
1720
16
  if (isl_seq_first_non_zero(mat->row[row], mat->n_col) == -1)
1721
0
    return mat;
1722
16
  mat = isl_mat_cow(mat);
1723
16
  if (!mat)
1724
0
    return NULL;
1725
16
  isl_seq_neg(mat->row[row], mat->row[row], mat->n_col);
1726
16
  return mat;
1727
16
}
1728
1729
__isl_give isl_mat *isl_mat_unimodular_complete(__isl_take isl_mat *M, int row)
1730
23.6k
{
1731
23.6k
  int r;
1732
23.6k
  struct isl_mat *H = NULL, *Q = NULL;
1733
23.6k
1734
23.6k
  if (!M)
1735
0
    return NULL;
1736
23.6k
1737
23.6k
  isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1738
23.6k
  M->n_row = row;
1739
23.6k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1740
23.6k
  M->n_row = M->n_col;
1741
23.6k
  if (!H)
1742
0
    goto error;
1743
47.2k
  
for (r = 0; 23.6k
r < row;
++r23.6k
)
1744
23.6k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1745
98.1k
  
for (r = row; 23.6k
r < M->n_row;
++r74.5k
)
1746
74.5k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1747
23.6k
  isl_mat_free(H);
1748
23.6k
  isl_mat_free(Q);
1749
23.6k
  return M;
1750
0
error:
1751
0
  isl_mat_free(H);
1752
0
  isl_mat_free(Q);
1753
0
  isl_mat_free(M);
1754
0
  return NULL;
1755
23.6k
}
1756
1757
__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1758
  __isl_take isl_mat *bot)
1759
455
{
1760
455
  struct isl_mat *mat;
1761
455
1762
455
  if (!top || !bot)
1763
0
    goto error;
1764
455
1765
455
  isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1766
455
  if (top->n_row == 0) {
1767
402
    isl_mat_free(top);
1768
402
    return bot;
1769
402
  }
1770
53
  if (bot->n_row == 0) {
1771
0
    isl_mat_free(bot);
1772
0
    return top;
1773
0
  }
1774
53
1775
53
  mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1776
53
  if (!mat)
1777
0
    goto error;
1778
53
  isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1779
53
       0, 0, mat->n_col);
1780
53
  isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1781
53
       0, 0, mat->n_col);
1782
53
  isl_mat_free(top);
1783
53
  isl_mat_free(bot);
1784
53
  return mat;
1785
0
error:
1786
0
  isl_mat_free(top);
1787
0
  isl_mat_free(bot);
1788
0
  return NULL;
1789
53
}
1790
1791
isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1792
48.5k
{
1793
48.5k
  int i;
1794
48.5k
1795
48.5k
  if (!mat1 || !mat2)
1796
0
    return isl_bool_error;
1797
48.5k
1798
48.5k
  if (mat1->n_row != mat2->n_row)
1799
7.53k
    return isl_bool_false;
1800
41.0k
1801
41.0k
  if (mat1->n_col != mat2->n_col)
1802
0
    return isl_bool_false;
1803
41.0k
1804
75.6k
  
for (i = 0; 41.0k
i < mat1->n_row;
++i34.6k
)
1805
35.0k
    if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))
1806
467
      return isl_bool_false;
1807
41.0k
1808
41.0k
  
return isl_bool_true40.5k
;
1809
41.0k
}
1810
1811
__isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1812
0
{
1813
0
  struct isl_mat *mat;
1814
0
1815
0
  if (!vec)
1816
0
    return NULL;
1817
0
  mat = isl_mat_alloc(vec->ctx, 1, vec->size);
1818
0
  if (!mat)
1819
0
    goto error;
1820
0
1821
0
  isl_seq_cpy(mat->row[0], vec->el, vec->size);
1822
0
1823
0
  isl_vec_free(vec);
1824
0
  return mat;
1825
0
error:
1826
0
  isl_vec_free(vec);
1827
0
  return NULL;
1828
0
}
1829
1830
/* Return a copy of row "row" of "mat" as an isl_vec.
1831
 */
1832
__isl_give isl_vec *isl_mat_get_row(__isl_keep isl_mat *mat, unsigned row)
1833
24
{
1834
24
  isl_vec *v;
1835
24
1836
24
  if (!mat)
1837
0
    return NULL;
1838
24
  if (row >= mat->n_row)
1839
24
    
isl_die0
(mat->ctx, isl_error_invalid, "row out of range",
1840
24
      return NULL);
1841
24
1842
24
  v = isl_vec_alloc(isl_mat_get_ctx(mat), mat->n_col);
1843
24
  if (!v)
1844
0
    return NULL;
1845
24
  isl_seq_cpy(v->el, mat->row[row], mat->n_col);
1846
24
1847
24
  return v;
1848
24
}
1849
1850
__isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1851
  __isl_take isl_vec *bot)
1852
0
{
1853
0
  return isl_mat_concat(top, isl_mat_from_row_vec(bot));
1854
0
}
1855
1856
__isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1857
  unsigned dst_col, unsigned src_col, unsigned n)
1858
626
{
1859
626
  isl_mat *res;
1860
626
1861
626
  if (!mat)
1862
0
    return NULL;
1863
626
  if (n == 0 || dst_col == src_col)
1864
496
    return mat;
1865
130
1866
130
  res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1867
130
  if (!res)
1868
0
    goto error;
1869
130
1870
130
  if (dst_col < src_col) {
1871
130
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1872
130
         0, 0, dst_col);
1873
130
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1874
130
         dst_col, src_col, n);
1875
130
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1876
130
         dst_col + n, dst_col, src_col - dst_col);
1877
130
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1878
130
         src_col + n, src_col + n,
1879
130
         res->n_col - src_col - n);
1880
130
  } else {
1881
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1882
0
         0, 0, src_col);
1883
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1884
0
         src_col, src_col + n, dst_col - src_col);
1885
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1886
0
         dst_col, src_col, n);
1887
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1888
0
         dst_col + n, dst_col + n,
1889
0
         res->n_col - dst_col - n);
1890
0
  }
1891
130
  isl_mat_free(mat);
1892
130
1893
130
  return res;
1894
0
error:
1895
0
  isl_mat_free(mat);
1896
0
  return NULL;
1897
130
}
1898
1899
/* Return the gcd of the elements in row "row" of "mat" in *gcd.
1900
 * Return isl_stat_ok on success and isl_stat_error on failure.
1901
 */
1902
isl_stat isl_mat_row_gcd(__isl_keep isl_mat *mat, int row, isl_int *gcd)
1903
7
{
1904
7
  if (check_row(mat, row) < 0)
1905
0
    return isl_stat_error;
1906
7
1907
7
  isl_seq_gcd(mat->row[row], mat->n_col, gcd);
1908
7
1909
7
  return isl_stat_ok;
1910
7
}
1911
1912
void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1913
778
{
1914
778
  int i;
1915
778
  isl_int g;
1916
778
1917
778
  isl_int_set_si(*gcd, 0);
1918
778
  if (!mat)
1919
0
    return;
1920
778
1921
778
  isl_int_init(g);
1922
4.99k
  for (i = 0; i < mat->n_row; 
++i4.21k
) {
1923
4.21k
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1924
4.21k
    isl_int_gcd(*gcd, *gcd, g);
1925
4.21k
  }
1926
778
  isl_int_clear(g);
1927
778
}
1928
1929
/* Return the result of scaling "mat" by a factor of "m".
1930
 */
1931
__isl_give isl_mat *isl_mat_scale(__isl_take isl_mat *mat, isl_int m)
1932
0
{
1933
0
  int i;
1934
0
1935
0
  if (isl_int_is_one(m))
1936
0
    return mat;
1937
0
1938
0
  mat = isl_mat_cow(mat);
1939
0
  if (!mat)
1940
0
    return NULL;
1941
0
1942
0
  for (i = 0; i < mat->n_row; ++i)
1943
0
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
1944
0
1945
0
  return mat;
1946
0
}
1947
1948
__isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1949
778
{
1950
778
  int i;
1951
778
1952
778
  if (isl_int_is_one(m))
1953
778
    
return mat0
;
1954
778
1955
778
  mat = isl_mat_cow(mat);
1956
778
  if (!mat)
1957
0
    return NULL;
1958
778
1959
4.99k
  
for (i = 0; 778
i < mat->n_row;
++i4.21k
)
1960
4.21k
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1961
778
1962
778
  return mat;
1963
778
}
1964
1965
__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row,
1966
  isl_int m)
1967
10
{
1968
10
  if (isl_int_is_one(m))
1969
10
    
return mat0
;
1970
10
1971
10
  mat = isl_mat_cow(mat);
1972
10
  if (!mat)
1973
0
    return NULL;
1974
10
1975
10
  isl_seq_scale_down(mat->row[row], mat->row[row], m, mat->n_col);
1976
10
1977
10
  return mat;
1978
10
}
1979
1980
__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1981
778
{
1982
778
  isl_int gcd;
1983
778
1984
778
  if (!mat)
1985
0
    return NULL;
1986
778
1987
778
  isl_int_init(gcd);
1988
778
  isl_mat_gcd(mat, &gcd);
1989
778
  mat = isl_mat_scale_down(mat, gcd);
1990
778
  isl_int_clear(gcd);
1991
778
1992
778
  return mat;
1993
778
}
1994
1995
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1996
2.05k
{
1997
2.05k
  mat = isl_mat_cow(mat);
1998
2.05k
  if (!mat)
1999
0
    return NULL;
2000
2.05k
2001
2.05k
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
2002
2.05k
2003
2.05k
  return mat;
2004
2.05k
}
2005
2006
/* Number of initial non-zero columns.
2007
 */
2008
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
2009
1.26k
{
2010
1.26k
  int i;
2011
1.26k
2012
1.26k
  if (!mat)
2013
0
    return -1;
2014
1.26k
2015
2.07k
  
for (i = 0; 1.26k
i < mat->n_col;
++i813
)
2016
1.62k
    if (row_first_non_zero(mat->row, mat->n_row, i) < 0)
2017
812
      break;
2018
1.26k
2019
1.26k
  return i;
2020
1.26k
}
2021
2022
/* Return a basis for the space spanned by the rows of "mat".
2023
 * Any basis will do, so simply perform Gaussian elimination and
2024
 * remove the empty rows.
2025
 */
2026
__isl_give isl_mat *isl_mat_row_basis(__isl_take isl_mat *mat)
2027
0
{
2028
0
  return isl_mat_reverse_gauss(mat);
2029
0
}
2030
2031
/* Return rows that extend a basis of "mat1" to one
2032
 * that covers both "mat1" and "mat2".
2033
 * The Hermite normal form of the concatenation of the two matrices is
2034
 *
2035
 *                       [ Q1 ]
2036
 *  [ M1 ] = [ H1 0  0 ] [ Q2 ]
2037
 *  [ M2 ] = [ H2 H3 0 ] [ Q3 ]
2038
 *
2039
 * The number of columns in H1 and H3 determine the number of rows
2040
 * in Q1 and Q2.  Q1 is a basis for M1, while Q2 extends this basis
2041
 * to also cover M2.
2042
 */
2043
__isl_give isl_mat *isl_mat_row_basis_extension(
2044
  __isl_take isl_mat *mat1, __isl_take isl_mat *mat2)
2045
0
{
2046
0
  int n_row;
2047
0
  int r1, r, n1;
2048
0
  isl_mat *H, *Q;
2049
0
2050
0
  n1 = isl_mat_rows(mat1);
2051
0
  H = isl_mat_concat(mat1, mat2);
2052
0
  H = isl_mat_left_hermite(H, 0, NULL, &Q);
2053
0
  if (!H || !Q)
2054
0
    goto error;
2055
0
2056
0
  r1 = hermite_first_zero_col(H, 0, n1);
2057
0
  r = hermite_first_zero_col(H, r1, H->n_row);
2058
0
  n_row = isl_mat_rows(Q);
2059
0
  Q = isl_mat_drop_rows(Q, r, n_row - r);
2060
0
  Q = isl_mat_drop_rows(Q, 0, r1);
2061
0
2062
0
  isl_mat_free(H);
2063
0
  return Q;
2064
0
error:
2065
0
  isl_mat_free(H);
2066
0
  isl_mat_free(Q);
2067
0
  return NULL;
2068
0
}
2069
2070
/* Are the rows of "mat1" linearly independent of those of "mat2"?
2071
 * That is, is there no linear dependence among the combined rows
2072
 * that is not already present in either "mat1" or "mat2"?
2073
 * In other words, is the rank of "mat1" and "mat2" combined equal
2074
 * to the sum of the ranks of "mat1" and "mat2"?
2075
 */
2076
isl_bool isl_mat_has_linearly_independent_rows(__isl_keep isl_mat *mat1,
2077
  __isl_keep isl_mat *mat2)
2078
0
{
2079
0
  int r1, r2, r;
2080
0
  isl_mat *mat;
2081
0
2082
0
  r1 = isl_mat_rank(mat1);
2083
0
  if (r1 < 0)
2084
0
    return isl_bool_error;
2085
0
  if (r1 == 0)
2086
0
    return isl_bool_true;
2087
0
  r2 = isl_mat_rank(mat2);
2088
0
  if (r2 < 0)
2089
0
    return isl_bool_error;
2090
0
  if (r2 == 0)
2091
0
    return isl_bool_true;
2092
0
2093
0
  mat = isl_mat_concat(isl_mat_copy(mat1), isl_mat_copy(mat2));
2094
0
  r = isl_mat_rank(mat);
2095
0
  isl_mat_free(mat);
2096
0
  if (r < 0)
2097
0
    return isl_bool_error;
2098
0
  return r == r1 + r2;
2099
0
}