Coverage Report

Created: 2018-12-09 11:54

/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.76M
{
27
3.76M
  return mat ? mat->ctx : NULL;
28
3.76M
}
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.99M
{
56
5.99M
  int i;
57
5.99M
  struct isl_mat *mat;
58
5.99M
59
5.99M
  mat = isl_alloc_type(ctx, struct isl_mat);
60
5.99M
  if (!mat)
61
0
    return NULL;
62
5.99M
63
5.99M
  mat->row = NULL;
64
5.99M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
65
5.99M
  if (isl_blk_is_error(mat->block))
66
0
    goto error;
67
5.99M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
68
5.99M
  if (n_row && 
!mat->row4.46M
)
69
0
    goto error;
70
5.99M
71
33.2M
  
for (i = 0; 5.99M
i < n_row;
++i27.2M
)
72
27.2M
    mat->row[i] = mat->block.data + i * n_col;
73
5.99M
74
5.99M
  mat->ctx = ctx;
75
5.99M
  isl_ctx_ref(ctx);
76
5.99M
  mat->ref = 1;
77
5.99M
  mat->n_row = n_row;
78
5.99M
  mat->n_col = n_col;
79
5.99M
  mat->max_col = n_col;
80
5.99M
  mat->flags = 0;
81
5.99M
82
5.99M
  return mat;
83
0
error:
84
0
  isl_blk_free(ctx, mat->block);
85
0
  free(mat);
86
0
  return NULL;
87
5.99M
}
88
89
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
90
  unsigned n_row, unsigned n_col)
91
105k
{
92
105k
  int i;
93
105k
  isl_int *old;
94
105k
  isl_int **row;
95
105k
96
105k
  if (!mat)
97
0
    return NULL;
98
105k
99
105k
  if (mat->max_col >= n_col && 
mat->n_row >= n_row101k
) {
100
7.92k
    if (mat->n_col < n_col)
101
42
      mat->n_col = n_col;
102
7.92k
    return mat;
103
7.92k
  }
104
97.2k
105
97.2k
  if (mat->max_col < n_col) {
106
3.87k
    struct isl_mat *new_mat;
107
3.87k
108
3.87k
    if (n_row < mat->n_row)
109
9
      n_row = mat->n_row;
110
3.87k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
111
3.87k
    if (!new_mat)
112
0
      goto error;
113
30.4k
    
for (i = 0; 3.87k
i < mat->n_row;
++i26.6k
)
114
26.6k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
115
3.87k
    isl_mat_free(mat);
116
3.87k
    return new_mat;
117
3.87k
  }
118
93.4k
119
93.4k
  mat = isl_mat_cow(mat);
120
93.4k
  if (!mat)
121
0
    goto error;
122
93.4k
123
93.4k
  old = mat->block.data;
124
93.4k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
125
93.4k
  if (isl_blk_is_error(mat->block))
126
0
    goto error;
127
93.4k
  row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
128
93.4k
  if (n_row && !row)
129
0
    goto error;
130
93.4k
  mat->row = row;
131
93.4k
132
984k
  for (i = 0; i < mat->n_row; 
++i891k
)
133
891k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
134
331k
  for (i = mat->n_row; i < n_row; 
++i237k
)
135
237k
    mat->row[i] = mat->block.data + i * mat->max_col;
136
93.4k
  mat->n_row = n_row;
137
93.4k
  if (mat->n_col < n_col)
138
0
    mat->n_col = n_col;
139
93.4k
140
93.4k
  return mat;
141
0
error:
142
0
  isl_mat_free(mat);
143
0
  return NULL;
144
93.4k
}
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.84M
{
149
2.84M
  int i;
150
2.84M
  struct isl_mat *mat;
151
2.84M
152
2.84M
  mat = isl_alloc_type(ctx, struct isl_mat);
153
2.84M
  if (!mat)
154
0
    return NULL;
155
2.84M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
156
2.84M
  if (n_row && 
!mat->row1.60M
)
157
0
    goto error;
158
9.83M
  
for (i = 0; 2.84M
i < n_row;
++i6.99M
)
159
6.99M
    mat->row[i] = row[first_row+i] + first_col;
160
2.84M
  mat->ctx = ctx;
161
2.84M
  isl_ctx_ref(ctx);
162
2.84M
  mat->ref = 1;
163
2.84M
  mat->n_row = n_row;
164
2.84M
  mat->n_col = n_col;
165
2.84M
  mat->block = isl_blk_empty();
166
2.84M
  mat->flags = ISL_MAT_BORROWED;
167
2.84M
  return mat;
168
0
error:
169
0
  free(mat);
170
0
  return NULL;
171
2.84M
}
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
663k
{
176
663k
  if (!mat)
177
0
    return NULL;
178
663k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
179
663k
          first_col, n_col);
180
663k
}
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
422k
{
185
422k
  int i;
186
422k
187
1.40M
  for (i = 0; i < n_row; 
++i982k
)
188
982k
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
189
422k
}
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
149k
{
194
149k
  int i;
195
149k
196
463k
  for (i = 0; i < n_row; 
++i313k
)
197
313k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
198
149k
}
199
200
__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
201
2.13M
{
202
2.13M
  if (!mat)
203
0
    return NULL;
204
2.13M
205
2.13M
  mat->ref++;
206
2.13M
  return mat;
207
2.13M
}
208
209
__isl_give isl_mat *isl_mat_dup(__isl_keep isl_mat *mat)
210
372k
{
211
372k
  int i;
212
372k
  struct isl_mat *mat2;
213
372k
214
372k
  if (!mat)
215
2.70k
    return NULL;
216
370k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
217
370k
  if (!mat2)
218
0
    return NULL;
219
1.20M
  
for (i = 0; 370k
i < mat->n_row;
++i833k
)
220
833k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
221
370k
  return mat2;
222
370k
}
223
224
__isl_give isl_mat *isl_mat_cow(__isl_take isl_mat *mat)
225
4.16M
{
226
4.16M
  struct isl_mat *mat2;
227
4.16M
  if (!mat)
228
0
    return NULL;
229
4.16M
230
4.16M
  if (mat->ref == 1 && 
!4.11M
ISL_F_ISSET4.11M
(mat, ISL_MAT_BORROWED))
231
4.16M
    
return mat3.79M
;
232
367k
233
367k
  mat2 = isl_mat_dup(mat);
234
367k
  isl_mat_free(mat);
235
367k
  return mat2;
236
367k
}
237
238
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
239
12.3M
{
240
12.3M
  if (!mat)
241
1.40M
    return NULL;
242
10.9M
243
10.9M
  if (--mat->ref > 0)
244
2.13M
    return NULL;
245
8.84M
246
8.84M
  if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
247
8.84M
    
isl_blk_free(mat->ctx, mat->block)5.99M
;
248
8.84M
  isl_ctx_deref(mat->ctx);
249
8.84M
  free(mat->row);
250
8.84M
  free(mat);
251
8.84M
252
8.84M
  return NULL;
253
8.84M
}
254
255
int isl_mat_rows(__isl_keep isl_mat *mat)
256
517k
{
257
517k
  return mat ? mat->n_row : 
-10
;
258
517k
}
259
260
int isl_mat_cols(__isl_keep isl_mat *mat)
261
56.5k
{
262
56.5k
  return mat ? mat->n_col : 
-10
;
263
56.5k
}
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
5.03k
{
269
5.03k
  if (!mat)
270
0
    return isl_stat_error;
271
5.03k
  if (col < 0 || col >= mat->n_col)
272
5.03k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
273
5.03k
      "column out of range", return isl_stat_error);
274
5.03k
  return isl_stat_ok;
275
5.03k
}
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
5.04k
{
281
5.04k
  if (!mat)
282
0
    return isl_stat_error;
283
5.04k
  if (row < 0 || row >= mat->n_row)
284
5.04k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
285
5.04k
      "row out of range", return isl_stat_error);
286
5.04k
  return isl_stat_ok;
287
5.04k
}
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.65M
{
294
1.65M
  if (!mat)
295
0
    return isl_stat_error;
296
1.65M
  if (first + n > mat->n_col || first + n < first)
297
1.65M
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
298
1.65M
      "column position or range out of bounds",
299
1.65M
      return isl_stat_error);
300
1.65M
  return isl_stat_ok;
301
1.65M
}
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.39M
{
308
4.39M
  if (!mat)
309
0
    return isl_stat_error;
310
4.39M
  if (first + n > mat->n_row || first + n < first)
311
4.39M
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
312
4.39M
      "row position or range out of bounds",
313
4.39M
      return isl_stat_error);
314
4.39M
  return isl_stat_ok;
315
4.39M
}
316
317
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
318
3.44k
{
319
3.44k
  if (check_row(mat, row) < 0)
320
0
    return -1;
321
3.44k
  if (check_col(mat, col) < 0)
322
0
    return -1;
323
3.44k
  isl_int_set(*v, mat->row[row][col]);
324
3.44k
  return 0;
325
3.44k
}
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.51k
{
345
1.51k
  mat = isl_mat_cow(mat);
346
1.51k
  if (check_row(mat, row) < 0)
347
0
    return isl_mat_free(mat);
348
1.51k
  if (check_col(mat, col) < 0)
349
0
    return isl_mat_free(mat);
350
1.51k
  isl_int_set(mat->row[row][col], v);
351
1.51k
  return mat;
352
1.51k
}
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
1.01M
{
386
1.01M
  int i;
387
1.01M
  struct isl_mat *mat;
388
1.01M
389
1.01M
  mat = isl_mat_alloc(ctx, n_row, n_row);
390
1.01M
  if (!mat)
391
0
    return NULL;
392
5.61M
  
for (i = 0; 1.01M
i < n_row;
++i4.59M
) {
393
4.59M
    isl_seq_clr(mat->row[i], i);
394
4.59M
    isl_int_set(mat->row[i][i], d);
395
4.59M
    isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
396
4.59M
  }
397
1.01M
398
1.01M
  return mat;
399
1.01M
}
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
1.01M
{
419
1.01M
  if (!ctx)
420
0
    return NULL;
421
1.01M
  return isl_mat_diag(ctx, n_row, ctx->one);
422
1.01M
}
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
335k
{
451
335k
  int i;
452
335k
  struct isl_vec *prod;
453
335k
454
335k
  if (!mat || !vec)
455
0
    goto error;
456
335k
457
335k
  isl_assert(mat->ctx, mat->n_col == vec->size, goto error);
458
335k
459
335k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
460
335k
  if (!prod)
461
0
    goto error;
462
335k
463
2.62M
  
for (i = 0; 335k
i < prod->size;
++i2.29M
)
464
2.29M
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
465
2.29M
          &prod->block.data[i]);
466
335k
  isl_mat_free(mat);
467
335k
  isl_vec_free(vec);
468
335k
  return prod;
469
0
error:
470
0
  isl_mat_free(mat);
471
0
  isl_vec_free(vec);
472
0
  return NULL;
473
335k
}
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
22.7k
{
507
22.7k
  int i, j;
508
22.7k
  struct isl_vec *prod;
509
22.7k
510
22.7k
  if (!mat || !vec)
511
0
    goto error;
512
22.7k
513
22.7k
  isl_assert(mat->ctx, mat->n_row == vec->size, goto error);
514
22.7k
515
22.7k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
516
22.7k
  if (!prod)
517
0
    goto error;
518
22.7k
519
132k
  
for (i = 0; 22.7k
i < prod->size;
++i109k
) {
520
109k
    isl_int_set_si(prod->el[i], 0);
521
892k
    for (j = 0; j < vec->size; 
++j782k
)
522
782k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
523
109k
  }
524
22.7k
  isl_mat_free(mat);
525
22.7k
  isl_vec_free(vec);
526
22.7k
  return prod;
527
0
error:
528
0
  isl_mat_free(mat);
529
0
  isl_vec_free(vec);
530
0
  return NULL;
531
22.7k
}
532
533
__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
534
  __isl_take isl_mat *right)
535
149k
{
536
149k
  int i;
537
149k
  struct isl_mat *sum;
538
149k
539
149k
  if (!left || !right)
540
0
    goto error;
541
149k
542
149k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
543
149k
  isl_assert(left->ctx, left->n_row >= 1, goto error);
544
149k
  isl_assert(left->ctx, left->n_col >= 1, goto error);
545
149k
  isl_assert(left->ctx, right->n_col >= 1, goto error);
546
149k
  isl_assert(left->ctx,
547
149k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
548
149k
      goto error);
549
149k
  isl_assert(left->ctx,
550
149k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
551
149k
      goto error);
552
149k
553
149k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
554
149k
  if (!sum)
555
0
    goto error;
556
149k
  isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
557
149k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
558
149k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
559
149k
560
149k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
561
1.15M
  for (i = 1; i < sum->n_row; 
++i1.00M
) {
562
1.00M
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
563
1.00M
    isl_int_addmul(sum->row[i][0],
564
1.00M
        right->row[0][0], right->row[i][0]);
565
1.00M
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
566
1.00M
        left->n_col-1);
567
1.00M
    isl_seq_scale(sum->row[i]+left->n_col,
568
1.00M
        right->row[i]+1, right->row[0][0],
569
1.00M
        right->n_col-1);
570
1.00M
  }
571
149k
572
149k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
573
149k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
574
149k
  isl_mat_free(left);
575
149k
  isl_mat_free(right);
576
149k
  return sum;
577
0
error:
578
0
  isl_mat_free(left);
579
0
  isl_mat_free(right);
580
0
  return NULL;
581
149k
}
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
918k
{
586
918k
  int r;
587
6.92M
  for (r = row; r < M->n_row; 
++r6.00M
)
588
6.00M
    isl_int_swap(M->row[r][i], M->row[r][j]);
589
918k
  if (U) {
590
8.98M
    for (r = 0; r < (*U)->n_row; 
++r8.10M
)
591
8.10M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
592
879k
  }
593
918k
  if (Q)
594
343k
    isl_mat_swap_rows(*Q, i, j);
595
918k
}
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
732k
{
600
732k
  int r;
601
2.80M
  for (r = row; r < M->n_row; 
++r2.07M
)
602
2.07M
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
603
732k
  if (U) {
604
4.55M
    for (r = 0; r < (*U)->n_row; 
++r3.88M
)
605
3.88M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
606
674k
  }
607
732k
  if (Q) {
608
1.30M
    for (r = 0; r < (*Q)->n_col; 
++r1.10M
)
609
1.10M
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
610
196k
  }
611
732k
}
612
613
static void oppose(struct isl_mat *M, struct isl_mat **U,
614
  struct isl_mat **Q, unsigned row, unsigned col)
615
286k
{
616
286k
  int r;
617
1.23M
  for (r = row; r < M->n_row; 
++r949k
)
618
949k
    isl_int_neg(M->row[r][col], M->row[r][col]);
619
286k
  if (U) {
620
1.58M
    for (r = 0; r < (*U)->n_row; 
++r1.33M
)
621
1.33M
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
622
241k
  }
623
286k
  if (Q)
624
140k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
625
286k
}
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
491k
{
642
491k
  isl_int c;
643
491k
  int row, col;
644
491k
645
491k
  if (U)
646
444k
    *U = NULL;
647
491k
  if (Q)
648
238k
    *Q = NULL;
649
491k
  if (!M)
650
0
    goto error;
651
491k
  M = isl_mat_cow(M);
652
491k
  if (!M)
653
0
    goto error;
654
491k
  if (U) {
655
444k
    *U = isl_mat_identity(M->ctx, M->n_col);
656
444k
    if (!*U)
657
0
      goto error;
658
491k
  }
659
491k
  if (Q) {
660
238k
    *Q = isl_mat_identity(M->ctx, M->n_col);
661
238k
    if (!*Q)
662
0
      goto error;
663
491k
  }
664
491k
665
491k
  col = 0;
666
491k
  isl_int_init(c);
667
2.24M
  for (row = 0; row < M->n_row; 
++row1.75M
) {
668
1.75M
    int first, i, off;
669
1.75M
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
670
1.75M
    if (first == -1)
671
485k
      continue;
672
1.26M
    first += col;
673
1.26M
    if (first != col)
674
666k
      exchange(M, U, Q, row, first, col);
675
1.26M
    if (isl_int_is_neg(M->row[row][col]))
676
1.26M
      
oppose(M, U, Q, row, col)240k
;
677
1.26M
    first = col+1;
678
1.63M
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
679
1.63M
                   M->n_col-first)) != -1) {
680
368k
      first += off;
681
368k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
682
368k
      subtract(M, U, Q, row, col, first, c);
683
368k
      if (!isl_int_is_zero(M->row[row][first]))
684
368k
        
exchange(M, U, Q, row, first, col)23.2k
;
685
345k
      else
686
345k
        ++first;
687
368k
    }
688
5.24M
    for (i = 0; i < col; 
++i3.98M
) {
689
3.98M
      if (isl_int_is_zero(M->row[row][i]))
690
3.98M
        
continue3.83M
;
691
144k
      if (neg)
692
144k
        
isl_int_cdiv_q0
(c, M->row[row][i], M->row[row][col]);
693
144k
      else
694
144k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
695
144k
      if (isl_int_is_zero(c))
696
144k
        
continue45.2k
;
697
99.5k
      subtract(M, U, Q, row, col, i, c);
698
99.5k
    }
699
1.26M
    ++col;
700
1.26M
  }
701
491k
  isl_int_clear(c);
702
491k
703
491k
  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
491k
}
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.19k
{
721
1.19k
  int k, nr, nc;
722
1.19k
  isl_ctx *ctx;
723
1.19k
724
1.19k
  if (!mat)
725
0
    return NULL;
726
1.19k
727
1.19k
  ctx = isl_mat_get_ctx(mat);
728
1.19k
  nr = isl_mat_rows(mat);
729
1.19k
  nc = isl_mat_cols(mat);
730
1.19k
731
3.31k
  for (k = 0; k < nr; 
++k2.12k
) {
732
2.12k
    if (k == row)
733
1.19k
      continue;
734
922
    if (isl_int_is_zero(mat->row[k][col]))
735
922
      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.19k
743
1.19k
  return mat;
744
1.19k
}
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.27k
{
757
1.27k
  int k, row, last, nr, nc;
758
1.27k
759
1.27k
  if (!mat)
760
0
    return NULL;
761
1.27k
762
1.27k
  nr = isl_mat_rows(mat);
763
1.27k
  nc = isl_mat_cols(mat);
764
1.27k
765
1.27k
  last = nc - 1;
766
2.46k
  for (row = nr - 1; row >= 0; 
--row1.19k
) {
767
1.58k
    for (; last >= 0; 
--last391
) {
768
2.06k
      for (k = row; k >= 0; 
--k471
)
769
1.66k
        if (!isl_int_is_zero(mat->row[k][last]))
770
1.66k
          
break1.19k
;
771
1.58k
      if (k >= 0)
772
1.19k
        break;
773
1.58k
    }
774
1.19k
    if (last < 0)
775
0
      break;
776
1.19k
    if (k != row)
777
0
      mat = isl_mat_swap_rows(mat, k, row);
778
1.19k
    if (!mat)
779
0
      return NULL;
780
1.19k
    if (isl_int_is_neg(mat->row[row][last]))
781
1.19k
      
mat = isl_mat_row_neg(mat, row)4
;
782
1.19k
    mat = eliminate(mat, row, last);
783
1.19k
    if (!mat)
784
0
      return NULL;
785
1.19k
  }
786
1.27k
  mat = isl_mat_drop_rows(mat, 0, row + 1);
787
1.27k
788
1.27k
  return mat;
789
1.27k
}
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.27k
{
796
1.27k
  int i, nr, nc;
797
1.27k
798
1.27k
  if (!mat)
799
0
    return NULL;
800
1.27k
801
1.27k
  nr = isl_mat_rows(mat);
802
1.27k
  nc = isl_mat_cols(mat);
803
1.27k
804
2.46k
  for (i = 0; i < nr; 
++i1.19k
) {
805
1.19k
    int pos;
806
1.19k
807
1.19k
    pos = isl_seq_first_non_zero(mat->row[i], nc);
808
1.19k
    if (pos < 0)
809
0
      continue;
810
1.19k
    if (isl_int_is_nonneg(mat->row[i][pos]))
811
1.19k
      
continue1.18k
;
812
12
    mat = isl_mat_row_neg(mat, i);
813
12
    if (!mat)
814
0
      return NULL;
815
12
  }
816
1.27k
817
1.27k
  return mat;
818
1.27k
}
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
709k
{
888
709k
  int i;
889
709k
  struct isl_mat *mat2;
890
709k
891
709k
  if (!mat)
892
0
    return NULL;
893
709k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
894
709k
  if (!mat2)
895
0
    goto error;
896
709k
  isl_int_set_si(mat2->row[0][0], 1);
897
709k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
898
4.38M
  for (i = 0; i < mat->n_row; 
++i3.67M
) {
899
3.67M
    isl_int_set_si(mat2->row[1+i][0], 0);
900
3.67M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
901
3.67M
  }
902
709k
  isl_mat_free(mat);
903
709k
  return mat2;
904
0
error:
905
0
  isl_mat_free(mat);
906
0
  return NULL;
907
709k
}
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
137k
{
917
137k
  int i;
918
137k
  isl_mat *mat;
919
137k
920
137k
  if (!mat1 || !mat2)
921
0
    goto error;
922
137k
923
137k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
924
137k
               mat1->n_col + mat2->n_col);
925
137k
  if (!mat)
926
0
    goto error;
927
274k
  
for (i = 0; 137k
i < mat1->n_row;
++i137k
) {
928
137k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
929
137k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
930
137k
  }
931
882k
  for (i = 0; i < mat2->n_row; 
++i745k
) {
932
745k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
933
745k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
934
745k
                mat2->row[i], mat2->n_col);
935
745k
  }
936
137k
  isl_mat_free(mat1);
937
137k
  isl_mat_free(mat2);
938
137k
  return mat;
939
0
error:
940
0
  isl_mat_free(mat1);
941
0
  isl_mat_free(mat2);
942
0
  return NULL;
943
137k
}
944
945
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
946
1.02M
{
947
1.02M
  int i;
948
1.02M
949
1.74M
  for (i = 0; i < n_row; 
++i719k
)
950
1.23M
    if (!isl_int_is_zero(row[i][col]))
951
1.23M
      
return i515k
;
952
1.02M
  
return -1511k
;
953
1.02M
}
954
955
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
956
510k
{
957
510k
  int i, min = row_first_non_zero(row, n_row, col);
958
510k
  if (min < 0)
959
0
    return -1;
960
1.23M
  
for (i = min + 1; 510k
i < n_row;
++i720k
) {
961
720k
    if (isl_int_is_zero(row[i][col]))
962
720k
      
continue716k
;
963
4.02k
    if (isl_int_abs_lt(row[i][col], row[min][col]))
964
4.02k
      
min = i63
;
965
4.02k
  }
966
510k
  return min;
967
510k
}
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.02k
{
990
4.02k
  isl_int_neg(m, m);
991
4.02k
  isl_seq_combine(left->row[i]+row,
992
4.02k
      left->ctx->one, left->row[i]+row,
993
4.02k
      m, left->row[row]+row,
994
4.02k
      left->n_col-row);
995
4.02k
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
996
4.02k
      m, right->row[row], right->n_col);
997
4.02k
}
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
170k
{
1004
170k
  int row;
1005
170k
  isl_int a, b;
1006
170k
1007
170k
  if (!left || !right)
1008
0
    goto error;
1009
170k
1010
170k
  isl_assert(left->ctx, left->n_row == left->n_col, goto error);
1011
170k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
1012
170k
1013
170k
  if (left->n_row == 0) {
1014
0
    isl_mat_free(left);
1015
0
    return right;
1016
0
  }
1017
170k
1018
170k
  left = isl_mat_cow(left);
1019
170k
  right = isl_mat_cow(right);
1020
170k
  if (!left || !right)
1021
0
    goto error;
1022
170k
1023
170k
  isl_int_init(a);
1024
170k
  isl_int_init(b);
1025
681k
  for (row = 0; row < left->n_row; 
++row510k
) {
1026
510k
    int pivot, first, i, off;
1027
510k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
1028
510k
    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
510k
    pivot += row;
1034
510k
    if (pivot != row)
1035
438
      if (inv_exchange(&left, &right, pivot, row) < 0)
1036
0
        goto error;
1037
510k
    if (isl_int_is_neg(left->row[row][row]))
1038
510k
      
inv_oppose(left, right, row)133
;
1039
510k
    first = row+1;
1040
514k
    while ((off = row_first_non_zero(left->row+first,
1041
514k
          left->n_row-first, row)) != -1) {
1042
4.02k
      first += off;
1043
4.02k
      isl_int_fdiv_q(a, left->row[first][row],
1044
4.02k
          left->row[row][row]);
1045
4.02k
      inv_subtract(left, right, row, first, a);
1046
4.02k
      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.02k
      } else {
1050
4.02k
        ++first;
1051
4.02k
      }
1052
4.02k
    }
1053
1.23M
    
for (i = 0; 510k
i < row;
++i721k
) {
1054
721k
      if (isl_int_is_zero(left->row[i][row]))
1055
721k
        
continue721k
;
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
510k
  }
1068
170k
  isl_int_clear(b);
1069
170k
1070
170k
  isl_int_set(a, left->row[0][0]);
1071
510k
  for (row = 1; row < left->n_row; 
++row340k
)
1072
340k
    isl_int_lcm(a, a, left->row[row][row]);
1073
170k
  if (isl_int_is_zero(a)){
1074
0
    isl_int_clear(a);
1075
0
    isl_assert(left->ctx, 0, goto error);
1076
0
  }
1077
681k
  
for (row = 0; 170k
row < left->n_row;
++row510k
) {
1078
510k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
1079
510k
    if (isl_int_is_one(left->row[row][row]))
1080
510k
      
continue502k
;
1081
8.33k
    isl_seq_scale(right->row[row], right->row[row],
1082
8.33k
        left->row[row][row], right->n_col);
1083
8.33k
  }
1084
170k
  isl_int_clear(a);
1085
170k
1086
170k
  isl_mat_free(left);
1087
170k
  return right;
1088
0
error:
1089
0
  isl_mat_free(left);
1090
0
  isl_mat_free(right);
1091
0
  return NULL;
1092
170k
}
1093
1094
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
1095
60.1k
{
1096
60.1k
  int i;
1097
60.1k
1098
302k
  for (i = 0; i < mat->n_row; 
++i241k
)
1099
241k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
1100
60.1k
}
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
170k
{
1105
170k
  int i;
1106
170k
  isl_int tmp;
1107
170k
1108
170k
  isl_int_init(tmp);
1109
792k
  for (i = 0; i < mat->n_row; 
++i622k
) {
1110
622k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
1111
622k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
1112
622k
    isl_int_set(mat->row[i][dst], tmp);
1113
622k
  }
1114
170k
  isl_int_clear(tmp);
1115
170k
}
1116
1117
__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
1118
66.0k
{
1119
66.0k
  struct isl_mat *inv;
1120
66.0k
  int row;
1121
66.0k
  isl_int a, b;
1122
66.0k
1123
66.0k
  mat = isl_mat_cow(mat);
1124
66.0k
  if (!mat)
1125
0
    return NULL;
1126
66.0k
1127
66.0k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
1128
66.0k
  inv = isl_mat_cow(inv);
1129
66.0k
  if (!inv)
1130
0
    goto error;
1131
66.0k
1132
66.0k
  isl_int_init(a);
1133
66.0k
  isl_int_init(b);
1134
281k
  for (row = 0; row < mat->n_row; 
++row215k
) {
1135
215k
    int pivot, first, i, off;
1136
215k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
1137
215k
    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
215k
    pivot += row;
1143
215k
    if (pivot != row)
1144
33.2k
      exchange(mat, &inv, NULL, row, pivot, row);
1145
215k
    if (isl_int_is_neg(mat->row[row][row]))
1146
215k
      
oppose(mat, &inv, NULL, row, row)46.2k
;
1147
215k
    first = row+1;
1148
480k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
1149
480k
                mat->n_col-first)) != -1) {
1150
264k
      first += off;
1151
264k
      isl_int_fdiv_q(a, mat->row[row][first],
1152
264k
                mat->row[row][row]);
1153
264k
      subtract(mat, &inv, NULL, row, row, first, a);
1154
264k
      if (!isl_int_is_zero(mat->row[row][first]))
1155
264k
        
exchange(mat, &inv, NULL, row, row, first)196k
;
1156
67.9k
      else
1157
67.9k
        ++first;
1158
264k
    }
1159
574k
    for (i = 0; i < row; 
++i358k
) {
1160
358k
      if (isl_int_is_zero(mat->row[row][i]))
1161
358k
        
continue273k
;
1162
85.0k
      isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
1163
85.0k
      isl_int_divexact(b, mat->row[row][i], a);
1164
85.0k
      isl_int_divexact(a, mat->row[row][row], a);
1165
85.0k
      isl_int_neg(a, a);
1166
85.0k
      isl_mat_col_combine(mat, i, a, i, b, row);
1167
85.0k
      isl_mat_col_combine(inv, i, a, i, b, row);
1168
85.0k
    }
1169
215k
  }
1170
66.0k
  isl_int_clear(b);
1171
66.0k
1172
66.0k
  isl_int_set(a, mat->row[0][0]);
1173
215k
  for (row = 1; row < mat->n_row; 
++row149k
)
1174
149k
    isl_int_lcm(a, a, mat->row[row][row]);
1175
66.0k
  if (isl_int_is_zero(a)){
1176
0
    isl_int_clear(a);
1177
0
    goto error;
1178
0
  }
1179
281k
  
for (row = 0; 66.0k
row < mat->n_row;
++row215k
) {
1180
215k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1181
215k
    if (isl_int_is_one(mat->row[row][row]))
1182
215k
      
continue155k
;
1183
60.1k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1184
60.1k
  }
1185
66.0k
  isl_int_clear(a);
1186
66.0k
1187
66.0k
  isl_mat_free(mat);
1188
66.0k
1189
66.0k
  return inv;
1190
0
error:
1191
0
  isl_mat_free(mat);
1192
0
  isl_mat_free(inv);
1193
0
  return NULL;
1194
66.0k
}
1195
1196
__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1197
4.59k
{
1198
4.59k
  struct isl_mat *transpose = NULL;
1199
4.59k
  int i, j;
1200
4.59k
1201
4.59k
  if (!mat)
1202
0
    return NULL;
1203
4.59k
1204
4.59k
  if (mat->n_col == mat->n_row) {
1205
4.59k
    mat = isl_mat_cow(mat);
1206
4.59k
    if (!mat)
1207
0
      return NULL;
1208
26.6k
    
for (i = 0; 4.59k
i < mat->n_row;
++i22.0k
)
1209
83.9k
      
for (j = i + 1; 22.0k
j < mat->n_col;
++j61.9k
)
1210
61.9k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1211
4.59k
    return mat;
1212
4.59k
  }
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
741k
{
1229
741k
  int r;
1230
741k
1231
741k
  mat = isl_mat_cow(mat);
1232
741k
  if (check_col_range(mat, i, 1) < 0 ||
1233
741k
      check_col_range(mat, j, 1) < 0)
1234
0
    return isl_mat_free(mat);
1235
741k
1236
21.2M
  
for (r = 0; 741k
r < mat->n_row;
++r20.5M
)
1237
20.5M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1238
741k
  return mat;
1239
741k
}
1240
1241
__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1242
  unsigned i, unsigned j)
1243
2.13M
{
1244
2.13M
  isl_int *t;
1245
2.13M
1246
2.13M
  if (!mat)
1247
0
    return NULL;
1248
2.13M
  mat = isl_mat_cow(mat);
1249
2.13M
  if (check_row_range(mat, i, 1) < 0 ||
1250
2.13M
      check_row_range(mat, j, 1) < 0)
1251
0
    return isl_mat_free(mat);
1252
2.13M
1253
2.13M
  t = mat->row[i];
1254
2.13M
  mat->row[i] = mat->row[j];
1255
2.13M
  mat->row[j] = t;
1256
2.13M
  return mat;
1257
2.13M
}
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
2.02M
{
1267
2.02M
  int i, j, k;
1268
2.02M
  struct isl_mat *prod;
1269
2.02M
1270
2.02M
  if (!left || !right)
1271
0
    goto error;
1272
2.02M
  isl_assert(left->ctx, left->n_col == right->n_row, goto error);
1273
2.02M
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1274
2.02M
  if (!prod)
1275
0
    goto error;
1276
2.02M
  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.73M
  
for (i = 0; 2.02M
i < prod->n_row;
++i4.70M
) {
1284
39.1M
    for (j = 0; j < prod->n_col; 
++j34.4M
)
1285
34.4M
      isl_int_mul(prod->row[i][j],
1286
4.70M
            left->row[i][0], right->row[0][j]);
1287
45.1M
    for (k = 1; k < left->n_col; 
++k40.4M
) {
1288
40.4M
      if (isl_int_is_zero(left->row[i][k]))
1289
40.4M
        
continue35.4M
;
1290
47.4M
      
for (j = 0; 5.01M
j < prod->n_col;
++j42.4M
)
1291
42.4M
        isl_int_addmul(prod->row[i][j],
1292
5.01M
              left->row[i][k], right->row[k][j]);
1293
5.01M
    }
1294
4.70M
  }
1295
2.02M
  isl_mat_free(left);
1296
2.02M
  isl_mat_free(right);
1297
2.02M
  return prod;
1298
0
error:
1299
0
  isl_mat_free(left);
1300
0
  isl_mat_free(right);
1301
0
  return NULL;
1302
2.02M
}
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.69M
{
1318
1.69M
  int i;
1319
1.69M
  struct isl_mat *t;
1320
1.69M
  int e;
1321
1.69M
1322
1.69M
  if (mat->n_col >= mat->n_row)
1323
1.04M
    e = 0;
1324
644k
  else
1325
644k
    e = mat->n_row - mat->n_col;
1326
1.69M
  if (has_div)
1327
564k
    
for (i = 0; 563k
i < n;
++i147
)
1328
563k
      
isl_int_mul147
(q[i][0], q[i][0], mat->row[0][0]);
1329
1.69M
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1330
1.69M
  t = isl_mat_product(t, mat);
1331
1.69M
  if (!t)
1332
0
    return -1;
1333
4.58M
  
for (i = 0; 1.69M
i < n;
++i2.89M
) {
1334
2.89M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1335
2.89M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1336
2.89M
          q[i] + has_div + t->n_col + e, n_div);
1337
2.89M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1338
2.89M
  }
1339
1.69M
  isl_mat_free(t);
1340
1.69M
  return 0;
1341
1.69M
}
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
563k
{
1356
563k
  struct isl_ctx *ctx;
1357
563k
1358
563k
  if (!bset || !mat)
1359
0
    goto error;
1360
563k
1361
563k
  ctx = bset->ctx;
1362
563k
  bset = isl_basic_set_cow(bset);
1363
563k
  if (!bset)
1364
0
    goto error;
1365
563k
1366
563k
  isl_assert(ctx, bset->dim->nparam == 0, goto error);
1367
563k
  isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1368
563k
  isl_assert(ctx, mat->n_col > 0, goto error);
1369
563k
1370
563k
  if (mat->n_col > mat->n_row) {
1371
101k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1372
101k
    if (!bset)
1373
0
      goto error;
1374
462k
  } else if (mat->n_col < mat->n_row) {
1375
214k
    bset->dim = isl_space_cow(bset->dim);
1376
214k
    if (!bset->dim)
1377
0
      goto error;
1378
214k
    bset->dim->n_out -= mat->n_row - mat->n_col;
1379
214k
  }
1380
563k
1381
563k
  if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1382
563k
      isl_mat_copy(mat)) < 0)
1383
0
    goto error;
1384
563k
1385
563k
  if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1386
563k
      isl_mat_copy(mat)) < 0)
1387
0
    goto error;
1388
563k
1389
563k
  if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
1390
0
    goto error2;
1391
563k
1392
563k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1393
563k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1394
563k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1395
563k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1396
563k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1397
563k
1398
563k
  bset = isl_basic_set_simplify(bset);
1399
563k
  bset = isl_basic_set_finalize(bset);
1400
563k
1401
563k
  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
32.3k
{
1412
32.3k
  int i;
1413
32.3k
1414
32.3k
  set = isl_set_cow(set);
1415
32.3k
  if (!set)
1416
0
    goto error;
1417
32.3k
1418
91.7k
  
for (i = 0; 32.3k
i < set->n;
++i59.3k
) {
1419
59.3k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1420
59.3k
                isl_mat_copy(mat));
1421
59.3k
    if (!set->p[i])
1422
0
      goto error;
1423
59.3k
  }
1424
32.3k
  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
32.3k
  isl_mat_free(mat);
1432
32.3k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1433
32.3k
  return set;
1434
0
error:
1435
0
  isl_set_free(set);
1436
0
  isl_mat_free(mat);
1437
0
  return NULL;
1438
32.3k
}
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
85.8k
{
1501
85.8k
  int r;
1502
85.8k
1503
85.8k
  if (n == 0)
1504
22
    return mat;
1505
85.8k
1506
85.8k
  mat = isl_mat_cow(mat);
1507
85.8k
  if (check_col_range(mat, col, n) < 0)
1508
0
    return isl_mat_free(mat);
1509
85.8k
1510
85.8k
  if (col != mat->n_col-n) {
1511
66.7k
    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
85.8k
  mat->n_col -= n;
1516
85.8k
  return mat;
1517
85.8k
}
1518
1519
__isl_give isl_mat *isl_mat_drop_rows(__isl_take isl_mat *mat,
1520
  unsigned row, unsigned n)
1521
122k
{
1522
122k
  int r;
1523
122k
1524
122k
  mat = isl_mat_cow(mat);
1525
122k
  if (check_row_range(mat, row, n) < 0)
1526
0
    return isl_mat_free(mat);
1527
122k
1528
384k
  
for (r = row; 122k
r+n < mat->n_row;
++r261k
)
1529
261k
    mat->row[r] = mat->row[r+n];
1530
122k
1531
122k
  mat->n_row -= n;
1532
122k
  return mat;
1533
122k
}
1534
1535
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1536
        unsigned col, unsigned n)
1537
88.1k
{
1538
88.1k
  isl_mat *ext;
1539
88.1k
1540
88.1k
  if (check_col_range(mat, col, 0) < 0)
1541
0
    return isl_mat_free(mat);
1542
88.1k
  if (n == 0)
1543
504
    return mat;
1544
87.6k
1545
87.6k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1546
87.6k
  if (!ext)
1547
0
    goto error;
1548
87.6k
1549
87.6k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1550
87.6k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1551
87.6k
        col + n, col, mat->n_col - col);
1552
87.6k
1553
87.6k
  isl_mat_free(mat);
1554
87.6k
  return ext;
1555
0
error:
1556
0
  isl_mat_free(mat);
1557
0
  return NULL;
1558
87.6k
}
1559
1560
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1561
  unsigned first, unsigned n)
1562
88.1k
{
1563
88.1k
  int i;
1564
88.1k
1565
88.1k
  if (!mat)
1566
0
    return NULL;
1567
88.1k
  mat = isl_mat_insert_cols(mat, first, n);
1568
88.1k
  if (!mat)
1569
0
    return NULL;
1570
88.1k
1571
92.5k
  
for (i = 0; 88.1k
i < mat->n_row;
++i4.36k
)
1572
4.36k
    isl_seq_clr(mat->row[i] + first, n);
1573
88.1k
1574
88.1k
  return mat;
1575
88.1k
}
1576
1577
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1578
4.82k
{
1579
4.82k
  if (!mat)
1580
0
    return NULL;
1581
4.82k
1582
4.82k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1583
4.82k
}
1584
1585
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1586
        unsigned row, unsigned n)
1587
5.27k
{
1588
5.27k
  isl_mat *ext;
1589
5.27k
1590
5.27k
  if (check_row_range(mat, row, 0) < 0)
1591
0
    return isl_mat_free(mat);
1592
5.27k
  if (n == 0)
1593
500
    return mat;
1594
4.77k
1595
4.77k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1596
4.77k
  if (!ext)
1597
0
    goto error;
1598
4.77k
1599
4.77k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1600
4.77k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1601
4.77k
        mat->n_row - row, 0, 0, mat->n_col);
1602
4.77k
1603
4.77k
  isl_mat_free(mat);
1604
4.77k
  return ext;
1605
0
error:
1606
0
  isl_mat_free(mat);
1607
0
  return NULL;
1608
4.77k
}
1609
1610
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1611
5.26k
{
1612
5.26k
  if (!mat)
1613
0
    return NULL;
1614
5.26k
1615
5.26k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1616
5.26k
}
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.7k
{
1664
14.7k
  int i;
1665
14.7k
1666
74.7k
  for (i = 0; i < mat->n_row; 
++i59.9k
)
1667
59.9k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1668
14.7k
}
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.9k
{
1731
23.9k
  int r;
1732
23.9k
  struct isl_mat *H = NULL, *Q = NULL;
1733
23.9k
1734
23.9k
  if (!M)
1735
0
    return NULL;
1736
23.9k
1737
23.9k
  isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1738
23.9k
  M->n_row = row;
1739
23.9k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1740
23.9k
  M->n_row = M->n_col;
1741
23.9k
  if (!H)
1742
0
    goto error;
1743
47.9k
  
for (r = 0; 23.9k
r < row;
++r23.9k
)
1744
23.9k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1745
101k
  
for (r = row; 23.9k
r < M->n_row;
++r77.4k
)
1746
77.4k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1747
23.9k
  isl_mat_free(H);
1748
23.9k
  isl_mat_free(Q);
1749
23.9k
  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.9k
}
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
53.6k
{
1793
53.6k
  int i;
1794
53.6k
1795
53.6k
  if (!mat1 || !mat2)
1796
0
    return isl_bool_error;
1797
53.6k
1798
53.6k
  if (mat1->n_row != mat2->n_row)
1799
8.14k
    return isl_bool_false;
1800
45.4k
1801
45.4k
  if (mat1->n_col != mat2->n_col)
1802
0
    return isl_bool_false;
1803
45.4k
1804
84.5k
  
for (i = 0; 45.4k
i < mat1->n_row;
++i39.0k
)
1805
39.5k
    if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))
1806
467
      return isl_bool_false;
1807
45.4k
1808
45.4k
  
return isl_bool_true45.0k
;
1809
45.4k
}
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
782
{
1914
782
  int i;
1915
782
  isl_int g;
1916
782
1917
782
  isl_int_set_si(*gcd, 0);
1918
782
  if (!mat)
1919
0
    return;
1920
782
1921
782
  isl_int_init(g);
1922
5.01k
  for (i = 0; i < mat->n_row; 
++i4.23k
) {
1923
4.23k
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1924
4.23k
    isl_int_gcd(*gcd, *gcd, g);
1925
4.23k
  }
1926
782
  isl_int_clear(g);
1927
782
}
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
782
{
1950
782
  int i;
1951
782
1952
782
  if (isl_int_is_one(m))
1953
782
    
return mat0
;
1954
782
1955
782
  mat = isl_mat_cow(mat);
1956
782
  if (!mat)
1957
0
    return NULL;
1958
782
1959
5.01k
  
for (i = 0; 782
i < mat->n_row;
++i4.23k
)
1960
4.23k
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1961
782
1962
782
  return mat;
1963
782
}
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
782
{
1982
782
  isl_int gcd;
1983
782
1984
782
  if (!mat)
1985
0
    return NULL;
1986
782
1987
782
  isl_int_init(gcd);
1988
782
  isl_mat_gcd(mat, &gcd);
1989
782
  mat = isl_mat_scale_down(mat, gcd);
1990
782
  isl_int_clear(gcd);
1991
782
1992
782
  return mat;
1993
782
}
1994
1995
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1996
2.09k
{
1997
2.09k
  mat = isl_mat_cow(mat);
1998
2.09k
  if (!mat)
1999
0
    return NULL;
2000
2.09k
2001
2.09k
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
2002
2.09k
2003
2.09k
  return mat;
2004
2.09k
}
2005
2006
/* Number of initial non-zero columns.
2007
 */
2008
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
2009
1.27k
{
2010
1.27k
  int i;
2011
1.27k
2012
1.27k
  if (!mat)
2013
0
    return -1;
2014
1.27k
2015
2.09k
  
for (i = 0; 1.27k
i < mat->n_col;
++i821
)
2016
1.63k
    if (row_first_non_zero(mat->row, mat->n_row, i) < 0)
2017
817
      break;
2018
1.27k
2019
1.27k
  return i;
2020
1.27k
}
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
}