Coverage Report

Created: 2018-04-24 22:41

/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.36M
{
27
3.36M
  return mat ? mat->ctx : NULL;
28
3.36M
}
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.92M
{
56
5.92M
  int i;
57
5.92M
  struct isl_mat *mat;
58
5.92M
59
5.92M
  mat = isl_alloc_type(ctx, struct isl_mat);
60
5.92M
  if (!mat)
61
0
    return NULL;
62
5.92M
63
5.92M
  mat->row = NULL;
64
5.92M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
65
5.92M
  if (isl_blk_is_error(mat->block))
66
0
    goto error;
67
5.92M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
68
5.92M
  if (n_row && 
!mat->row4.36M
)
69
0
    goto error;
70
5.92M
71
32.4M
  
for (i = 0; 5.92M
i < n_row;
++i26.5M
)
72
26.5M
    mat->row[i] = mat->block.data + i * n_col;
73
5.92M
74
5.92M
  mat->ctx = ctx;
75
5.92M
  isl_ctx_ref(ctx);
76
5.92M
  mat->ref = 1;
77
5.92M
  mat->n_row = n_row;
78
5.92M
  mat->n_col = n_col;
79
5.92M
  mat->max_col = n_col;
80
5.92M
  mat->flags = 0;
81
5.92M
82
5.92M
  return mat;
83
0
error:
84
0
  isl_blk_free(ctx, mat->block);
85
0
  free(mat);
86
0
  return NULL;
87
5.92M
}
88
89
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
90
  unsigned n_row, unsigned n_col)
91
106k
{
92
106k
  int i;
93
106k
  isl_int *old;
94
106k
  isl_int **row;
95
106k
96
106k
  if (!mat)
97
0
    return NULL;
98
106k
99
106k
  if (mat->max_col >= n_col && 
mat->n_row >= n_row102k
) {
100
7.98k
    if (mat->n_col < n_col)
101
34
      mat->n_col = n_col;
102
7.98k
    return mat;
103
7.98k
  }
104
98.1k
105
98.1k
  if (mat->max_col < n_col) {
106
3.72k
    struct isl_mat *new_mat;
107
3.72k
108
3.72k
    if (n_row < mat->n_row)
109
9
      n_row = mat->n_row;
110
3.72k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
111
3.72k
    if (!new_mat)
112
0
      goto error;
113
26.4k
    
for (i = 0; 3.72k
i < mat->n_row;
++i22.7k
)
114
22.7k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
115
3.72k
    isl_mat_free(mat);
116
3.72k
    return new_mat;
117
3.72k
  }
118
94.3k
119
94.3k
  mat = isl_mat_cow(mat);
120
94.3k
  if (!mat)
121
0
    goto error;
122
94.3k
123
94.3k
  old = mat->block.data;
124
94.3k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
125
94.3k
  if (isl_blk_is_error(mat->block))
126
0
    goto error;
127
94.3k
  row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
128
94.3k
  if (n_row && !row)
129
0
    goto error;
130
94.3k
  mat->row = row;
131
94.3k
132
1.00M
  for (i = 0; i < mat->n_row; 
++i908k
)
133
908k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
134
332k
  for (i = mat->n_row; i < n_row; 
++i238k
)
135
238k
    mat->row[i] = mat->block.data + i * mat->max_col;
136
94.3k
  mat->n_row = n_row;
137
94.3k
  if (mat->n_col < n_col)
138
0
    mat->n_col = n_col;
139
94.3k
140
94.3k
  return mat;
141
0
error:
142
0
  isl_mat_free(mat);
143
0
  return NULL;
144
94.3k
}
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.82M
{
149
2.82M
  int i;
150
2.82M
  struct isl_mat *mat;
151
2.82M
152
2.82M
  mat = isl_alloc_type(ctx, struct isl_mat);
153
2.82M
  if (!mat)
154
0
    return NULL;
155
2.82M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
156
2.82M
  if (n_row && 
!mat->row1.56M
)
157
0
    goto error;
158
9.82M
  
for (i = 0; 2.82M
i < n_row;
++i6.99M
)
159
6.99M
    mat->row[i] = row[first_row+i] + first_col;
160
2.82M
  mat->ctx = ctx;
161
2.82M
  isl_ctx_ref(ctx);
162
2.82M
  mat->ref = 1;
163
2.82M
  mat->n_row = n_row;
164
2.82M
  mat->n_col = n_col;
165
2.82M
  mat->block = isl_blk_empty();
166
2.82M
  mat->flags = ISL_MAT_BORROWED;
167
2.82M
  return mat;
168
0
error:
169
0
  free(mat);
170
0
  return NULL;
171
2.82M
}
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
625k
{
176
625k
  if (!mat)
177
0
    return NULL;
178
625k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
179
625k
          first_col, n_col);
180
625k
}
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
423k
{
185
423k
  int i;
186
423k
187
1.46M
  for (i = 0; i < n_row; 
++i1.04M
)
188
1.04M
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
189
423k
}
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
144k
{
194
144k
  int i;
195
144k
196
440k
  for (i = 0; i < n_row; 
++i296k
)
197
296k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
198
144k
}
199
200
__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
201
2.15M
{
202
2.15M
  if (!mat)
203
0
    return NULL;
204
2.15M
205
2.15M
  mat->ref++;
206
2.15M
  return mat;
207
2.15M
}
208
209
__isl_give isl_mat *isl_mat_dup(__isl_keep isl_mat *mat)
210
365k
{
211
365k
  int i;
212
365k
  struct isl_mat *mat2;
213
365k
214
365k
  if (!mat)
215
2.75k
    return NULL;
216
362k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
217
362k
  if (!mat2)
218
0
    return NULL;
219
1.17M
  
for (i = 0; 362k
i < mat->n_row;
++i813k
)
220
813k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
221
362k
  return mat2;
222
362k
}
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
359k
233
359k
  mat2 = isl_mat_dup(mat);
234
359k
  isl_mat_free(mat);
235
359k
  return mat2;
236
359k
}
237
238
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
239
12.2M
{
240
12.2M
  if (!mat)
241
1.34M
    return NULL;
242
10.9M
243
10.9M
  if (--mat->ref > 0)
244
2.15M
    return NULL;
245
8.75M
246
8.75M
  if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
247
8.75M
    
isl_blk_free(mat->ctx, mat->block)5.92M
;
248
8.75M
  isl_ctx_deref(mat->ctx);
249
8.75M
  free(mat->row);
250
8.75M
  free(mat);
251
8.75M
252
8.75M
  return NULL;
253
8.75M
}
254
255
int isl_mat_rows(__isl_keep isl_mat *mat)
256
475k
{
257
475k
  return mat ? mat->n_row : 
-10
;
258
475k
}
259
260
int isl_mat_cols(__isl_keep isl_mat *mat)
261
45.1k
{
262
45.1k
  return mat ? mat->n_col : 
-10
;
263
45.1k
}
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.09k
{
269
5.09k
  if (!mat)
270
0
    return isl_stat_error;
271
5.09k
  if (col < 0 || col >= mat->n_col)
272
5.09k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
273
5.09k
      "column out of range", return isl_stat_error);
274
5.09k
  return isl_stat_ok;
275
5.09k
}
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.10k
{
281
5.10k
  if (!mat)
282
0
    return isl_stat_error;
283
5.10k
  if (row < 0 || row >= mat->n_row)
284
5.10k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
285
5.10k
      "row out of range", return isl_stat_error);
286
5.10k
  return isl_stat_ok;
287
5.10k
}
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.46M
{
294
1.46M
  if (!mat)
295
0
    return isl_stat_error;
296
1.46M
  if (first + n > mat->n_col || first + n < first)
297
1.46M
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
298
1.46M
      "column position or range out of bounds",
299
1.46M
      return isl_stat_error);
300
1.46M
  return isl_stat_ok;
301
1.46M
}
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.25M
{
308
4.25M
  if (!mat)
309
0
    return isl_stat_error;
310
4.25M
  if (first + n > mat->n_row || first + n < first)
311
4.25M
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
312
4.25M
      "row position or range out of bounds",
313
4.25M
      return isl_stat_error);
314
4.25M
  return isl_stat_ok;
315
4.25M
}
316
317
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
318
3.48k
{
319
3.48k
  if (check_row(mat, row) < 0)
320
0
    return -1;
321
3.48k
  if (check_col(mat, col) < 0)
322
0
    return -1;
323
3.48k
  isl_int_set(*v, mat->row[row][col]);
324
3.48k
  return 0;
325
3.48k
}
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.52k
{
345
1.52k
  mat = isl_mat_cow(mat);
346
1.52k
  if (check_row(mat, row) < 0)
347
0
    return isl_mat_free(mat);
348
1.52k
  if (check_col(mat, col) < 0)
349
0
    return isl_mat_free(mat);
350
1.52k
  isl_int_set(mat->row[row][col], v);
351
1.52k
  return mat;
352
1.52k
}
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.45M
  
for (i = 0; 1.01M
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
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
151
{
405
151
  int i;
406
151
  isl_mat *mat;
407
151
408
151
  mat = isl_mat_alloc(ctx, n_row, n_col);
409
151
  if (!mat)
410
0
    return NULL;
411
192
  
for (i = 0; 151
i < n_row;
++i41
)
412
41
    isl_seq_clr(mat->row[i], n_col);
413
151
414
151
  return mat;
415
151
}
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
343k
{
451
343k
  int i;
452
343k
  struct isl_vec *prod;
453
343k
454
343k
  if (!mat || !vec)
455
0
    goto error;
456
343k
457
343k
  isl_assert(mat->ctx, mat->n_col == vec->size, goto error);
458
343k
459
343k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
460
343k
  if (!prod)
461
0
    goto error;
462
343k
463
2.66M
  
for (i = 0; 343k
i < prod->size;
++i2.32M
)
464
2.32M
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
465
2.32M
          &prod->block.data[i]);
466
343k
  isl_mat_free(mat);
467
343k
  isl_vec_free(vec);
468
343k
  return prod;
469
0
error:
470
0
  isl_mat_free(mat);
471
0
  isl_vec_free(vec);
472
0
  return NULL;
473
343k
}
474
475
__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
476
  __isl_take isl_vec *vec)
477
159
{
478
159
  struct isl_mat *vec_mat;
479
159
  int i;
480
159
481
159
  if (!mat || !vec)
482
0
    goto error;
483
159
  vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
484
159
  if (!vec_mat)
485
0
    goto error;
486
1.17k
  
for (i = 0; 159
i < vec->size;
++i1.02k
)
487
1.02k
    isl_int_set(vec_mat->row[i][0], vec->el[i]);
488
159
  vec_mat = isl_mat_inverse_product(mat, vec_mat);
489
159
  isl_vec_free(vec);
490
159
  if (!vec_mat)
491
0
    return NULL;
492
159
  vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
493
159
  if (vec)
494
1.17k
    
for (i = 0; 159
i < vec->size;
++i1.02k
)
495
1.02k
      isl_int_set(vec->el[i], vec_mat->row[i][0]);
496
159
  isl_mat_free(vec_mat);
497
159
  return vec;
498
0
error:
499
0
  isl_mat_free(mat);
500
0
  isl_vec_free(vec);
501
0
  return NULL;
502
159
}
503
504
__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,
505
  __isl_take isl_mat *mat)
506
19.8k
{
507
19.8k
  int i, j;
508
19.8k
  struct isl_vec *prod;
509
19.8k
510
19.8k
  if (!mat || !vec)
511
0
    goto error;
512
19.8k
513
19.8k
  isl_assert(mat->ctx, mat->n_row == vec->size, goto error);
514
19.8k
515
19.8k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
516
19.8k
  if (!prod)
517
0
    goto error;
518
19.8k
519
107k
  
for (i = 0; 19.8k
i < prod->size;
++i87.5k
) {
520
87.5k
    isl_int_set_si(prod->el[i], 0);
521
659k
    for (j = 0; j < vec->size; 
++j571k
)
522
571k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
523
87.5k
  }
524
19.8k
  isl_mat_free(mat);
525
19.8k
  isl_vec_free(vec);
526
19.8k
  return prod;
527
0
error:
528
0
  isl_mat_free(mat);
529
0
  isl_vec_free(vec);
530
0
  return NULL;
531
19.8k
}
532
533
__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
534
  __isl_take isl_mat *right)
535
144k
{
536
144k
  int i;
537
144k
  struct isl_mat *sum;
538
144k
539
144k
  if (!left || !right)
540
0
    goto error;
541
144k
542
144k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
543
144k
  isl_assert(left->ctx, left->n_row >= 1, goto error);
544
144k
  isl_assert(left->ctx, left->n_col >= 1, goto error);
545
144k
  isl_assert(left->ctx, right->n_col >= 1, goto error);
546
144k
  isl_assert(left->ctx,
547
144k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
548
144k
      goto error);
549
144k
  isl_assert(left->ctx,
550
144k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
551
144k
      goto error);
552
144k
553
144k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
554
144k
  if (!sum)
555
0
    goto error;
556
144k
  isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
557
144k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
558
144k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
559
144k
560
144k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
561
1.08M
  for (i = 1; i < sum->n_row; 
++i942k
) {
562
942k
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
563
942k
    isl_int_addmul(sum->row[i][0],
564
942k
        right->row[0][0], right->row[i][0]);
565
942k
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
566
942k
        left->n_col-1);
567
942k
    isl_seq_scale(sum->row[i]+left->n_col,
568
942k
        right->row[i]+1, right->row[0][0],
569
942k
        right->n_col-1);
570
942k
  }
571
144k
572
144k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
573
144k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
574
144k
  isl_mat_free(left);
575
144k
  isl_mat_free(right);
576
144k
  return sum;
577
0
error:
578
0
  isl_mat_free(left);
579
0
  isl_mat_free(right);
580
0
  return NULL;
581
144k
}
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
888k
{
586
888k
  int r;
587
6.94M
  for (r = row; r < M->n_row; 
++r6.05M
)
588
6.05M
    isl_int_swap(M->row[r][i], M->row[r][j]);
589
888k
  if (U) {
590
8.68M
    for (r = 0; r < (*U)->n_row; 
++r7.82M
)
591
7.82M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
592
860k
  }
593
888k
  if (Q)
594
332k
    isl_mat_swap_rows(*Q, i, j);
595
888k
}
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
3.05M
  for (r = row; r < M->n_row; 
++r2.32M
)
602
2.32M
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
603
732k
  if (U) {
604
4.54M
    for (r = 0; r < (*U)->n_row; 
++r3.85M
)
605
3.85M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
606
694k
  }
607
732k
  if (Q) {
608
1.34M
    for (r = 0; r < (*Q)->n_col; 
++r1.12M
)
609
1.12M
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
610
219k
  }
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
257k
{
616
257k
  int r;
617
1.09M
  for (r = row; r < M->n_row; 
++r840k
)
618
840k
    isl_int_neg(M->row[r][col], M->row[r][col]);
619
257k
  if (U) {
620
1.46M
    for (r = 0; r < (*U)->n_row; 
++r1.23M
)
621
1.23M
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
622
224k
  }
623
257k
  if (Q)
624
123k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
625
257k
}
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
444k
    *U = NULL;
647
479k
  if (Q)
648
235k
    *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
444k
    *U = isl_mat_identity(M->ctx, M->n_col);
656
444k
    if (!*U)
657
0
      goto error;
658
479k
  }
659
479k
  if (Q) {
660
235k
    *Q = isl_mat_identity(M->ctx, M->n_col);
661
235k
    if (!*Q)
662
0
      goto error;
663
479k
  }
664
479k
665
479k
  col = 0;
666
479k
  isl_int_init(c);
667
2.25M
  for (row = 0; row < M->n_row; 
++row1.77M
) {
668
1.77M
    int first, i, off;
669
1.77M
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
670
1.77M
    if (first == -1)
671
567k
      continue;
672
1.21M
    first += col;
673
1.21M
    if (first != col)
674
631k
      exchange(M, U, Q, row, first, col);
675
1.21M
    if (isl_int_is_neg(M->row[row][col]))
676
1.21M
      
oppose(M, U, Q, row, col)215k
;
677
1.21M
    first = col+1;
678
1.57M
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
679
1.57M
                   M->n_col-first)) != -1) {
680
369k
      first += off;
681
369k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
682
369k
      subtract(M, U, Q, row, col, first, c);
683
369k
      if (!isl_int_is_zero(M->row[row][first]))
684
369k
        
exchange(M, U, Q, row, first, col)30.9k
;
685
338k
      else
686
338k
        ++first;
687
369k
    }
688
5.01M
    for (i = 0; i < col; 
++i3.80M
) {
689
3.80M
      if (isl_int_is_zero(M->row[row][i]))
690
3.80M
        
continue3.65M
;
691
154k
      if (neg)
692
154k
        
isl_int_cdiv_q0
(c, M->row[row][i], M->row[row][col]);
693
154k
      else
694
154k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
695
154k
      if (isl_int_is_zero(c))
696
154k
        
continue52.8k
;
697
102k
      subtract(M, U, Q, row, col, i, c);
698
102k
    }
699
1.21M
    ++col;
700
1.21M
  }
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.22k
{
721
1.22k
  int k, nr, nc;
722
1.22k
  isl_ctx *ctx;
723
1.22k
724
1.22k
  if (!mat)
725
0
    return NULL;
726
1.22k
727
1.22k
  ctx = isl_mat_get_ctx(mat);
728
1.22k
  nr = isl_mat_rows(mat);
729
1.22k
  nc = isl_mat_cols(mat);
730
1.22k
731
3.38k
  for (k = 0; k < nr; 
++k2.16k
) {
732
2.16k
    if (k == row)
733
1.22k
      continue;
734
942
    if (isl_int_is_zero(mat->row[k][col]))
735
942
      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.22k
743
1.22k
  return mat;
744
1.22k
}
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.28k
{
757
1.28k
  int k, row, last, nr, nc;
758
1.28k
759
1.28k
  if (!mat)
760
0
    return NULL;
761
1.28k
762
1.28k
  nr = isl_mat_rows(mat);
763
1.28k
  nc = isl_mat_cols(mat);
764
1.28k
765
1.28k
  last = nc - 1;
766
2.51k
  for (row = nr - 1; row >= 0; 
--row1.22k
) {
767
1.62k
    for (; last >= 0; 
--last399
) {
768
2.10k
      for (k = row; k >= 0; 
--k479
)
769
1.70k
        if (!isl_int_is_zero(mat->row[k][last]))
770
1.70k
          
break1.22k
;
771
1.62k
      if (k >= 0)
772
1.22k
        break;
773
1.62k
    }
774
1.22k
    if (last < 0)
775
0
      break;
776
1.22k
    if (k != row)
777
0
      mat = isl_mat_swap_rows(mat, k, row);
778
1.22k
    if (!mat)
779
0
      return NULL;
780
1.22k
    if (isl_int_is_neg(mat->row[row][last]))
781
1.22k
      
mat = isl_mat_row_neg(mat, row)4
;
782
1.22k
    mat = eliminate(mat, row, last);
783
1.22k
    if (!mat)
784
0
      return NULL;
785
1.22k
  }
786
1.28k
  mat = isl_mat_drop_rows(mat, 0, row + 1);
787
1.28k
788
1.28k
  return mat;
789
1.28k
}
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.28k
{
796
1.28k
  int i, nr, nc;
797
1.28k
798
1.28k
  if (!mat)
799
0
    return NULL;
800
1.28k
801
1.28k
  nr = isl_mat_rows(mat);
802
1.28k
  nc = isl_mat_cols(mat);
803
1.28k
804
2.51k
  for (i = 0; i < nr; 
++i1.22k
) {
805
1.22k
    int pos;
806
1.22k
807
1.22k
    pos = isl_seq_first_non_zero(mat->row[i], nc);
808
1.22k
    if (pos < 0)
809
0
      continue;
810
1.22k
    if (isl_int_is_nonneg(mat->row[i][pos]))
811
1.22k
      
continue1.21k
;
812
12
    mat = isl_mat_row_neg(mat, i);
813
12
    if (!mat)
814
0
      return NULL;
815
12
  }
816
1.28k
817
1.28k
  return mat;
818
1.28k
}
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
690k
{
888
690k
  int i;
889
690k
  struct isl_mat *mat2;
890
690k
891
690k
  if (!mat)
892
0
    return NULL;
893
690k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
894
690k
  if (!mat2)
895
0
    goto error;
896
690k
  isl_int_set_si(mat2->row[0][0], 1);
897
690k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
898
4.22M
  for (i = 0; i < mat->n_row; 
++i3.53M
) {
899
3.53M
    isl_int_set_si(mat2->row[1+i][0], 0);
900
3.53M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
901
3.53M
  }
902
690k
  isl_mat_free(mat);
903
690k
  return mat2;
904
0
error:
905
0
  isl_mat_free(mat);
906
0
  return NULL;
907
690k
}
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
131k
{
917
131k
  int i;
918
131k
  isl_mat *mat;
919
131k
920
131k
  if (!mat1 || !mat2)
921
0
    goto error;
922
131k
923
131k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
924
131k
               mat1->n_col + mat2->n_col);
925
131k
  if (!mat)
926
0
    goto error;
927
262k
  
for (i = 0; 131k
i < mat1->n_row;
++i131k
) {
928
131k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
929
131k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
930
131k
  }
931
837k
  for (i = 0; i < mat2->n_row; 
++i705k
) {
932
705k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
933
705k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
934
705k
                mat2->row[i], mat2->n_col);
935
705k
  }
936
131k
  isl_mat_free(mat1);
937
131k
  isl_mat_free(mat2);
938
131k
  return mat;
939
0
error:
940
0
  isl_mat_free(mat1);
941
0
  isl_mat_free(mat2);
942
0
  return NULL;
943
131k
}
944
945
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
946
957k
{
947
957k
  int i;
948
957k
949
1.62M
  for (i = 0; i < n_row; 
++i665k
)
950
1.14M
    if (!isl_int_is_zero(row[i][col]))
951
1.14M
      
return i480k
;
952
957k
  
return -1476k
;
953
957k
}
954
955
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
956
475k
{
957
475k
  int i, min = row_first_non_zero(row, n_row, col);
958
475k
  if (min < 0)
959
0
    return -1;
960
1.14M
  
for (i = min + 1; 475k
i < n_row;
++i667k
) {
961
667k
    if (isl_int_is_zero(row[i][col]))
962
667k
      
continue663k
;
963
3.58k
    if (isl_int_abs_lt(row[i][col], row[min][col]))
964
3.58k
      
min = i49
;
965
3.58k
  }
966
475k
  return min;
967
475k
}
968
969
static isl_stat inv_exchange(__isl_keep isl_mat **left,
970
  __isl_keep isl_mat **right, unsigned i, unsigned j)
971
340
{
972
340
  *left = isl_mat_swap_rows(*left, i, j);
973
340
  *right = isl_mat_swap_rows(*right, i, j);
974
340
975
340
  if (!*left || !*right)
976
0
    return isl_stat_error;
977
340
  return isl_stat_ok;
978
340
}
979
980
static void inv_oppose(
981
  struct isl_mat *left, struct isl_mat *right, unsigned row)
982
105
{
983
105
  isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
984
105
  isl_seq_neg(right->row[row], right->row[row], right->n_col);
985
105
}
986
987
static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
988
  unsigned row, unsigned i, isl_int m)
989
3.58k
{
990
3.58k
  isl_int_neg(m, m);
991
3.58k
  isl_seq_combine(left->row[i]+row,
992
3.58k
      left->ctx->one, left->row[i]+row,
993
3.58k
      m, left->row[row]+row,
994
3.58k
      left->n_col-row);
995
3.58k
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
996
3.58k
      m, right->row[row], right->n_col);
997
3.58k
}
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
159k
{
1004
159k
  int row;
1005
159k
  isl_int a, b;
1006
159k
1007
159k
  if (!left || !right)
1008
0
    goto error;
1009
159k
1010
159k
  isl_assert(left->ctx, left->n_row == left->n_col, goto error);
1011
159k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
1012
159k
1013
159k
  if (left->n_row == 0) {
1014
0
    isl_mat_free(left);
1015
0
    return right;
1016
0
  }
1017
159k
1018
159k
  left = isl_mat_cow(left);
1019
159k
  right = isl_mat_cow(right);
1020
159k
  if (!left || !right)
1021
0
    goto error;
1022
159k
1023
159k
  isl_int_init(a);
1024
159k
  isl_int_init(b);
1025
635k
  for (row = 0; row < left->n_row; 
++row475k
) {
1026
475k
    int pivot, first, i, off;
1027
475k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
1028
475k
    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
475k
    pivot += row;
1034
475k
    if (pivot != row)
1035
338
      if (inv_exchange(&left, &right, pivot, row) < 0)
1036
0
        goto error;
1037
475k
    if (isl_int_is_neg(left->row[row][row]))
1038
475k
      
inv_oppose(left, right, row)105
;
1039
475k
    first = row+1;
1040
479k
    while ((off = row_first_non_zero(left->row+first,
1041
479k
          left->n_row-first, row)) != -1) {
1042
3.58k
      first += off;
1043
3.58k
      isl_int_fdiv_q(a, left->row[first][row],
1044
3.58k
          left->row[row][row]);
1045
3.58k
      inv_subtract(left, right, row, first, a);
1046
3.58k
      if (!isl_int_is_zero(left->row[first][row])) {
1047
2
        if (inv_exchange(&left, &right, row, first) < 0)
1048
0
          goto error;
1049
3.58k
      } else {
1050
3.58k
        ++first;
1051
3.58k
      }
1052
3.58k
    }
1053
1.14M
    
for (i = 0; 475k
i < row;
++i668k
) {
1054
668k
      if (isl_int_is_zero(left->row[i][row]))
1055
668k
        
continue667k
;
1056
346
      isl_int_gcd(a, left->row[row][row], left->row[i][row]);
1057
346
      isl_int_divexact(b, left->row[i][row], a);
1058
346
      isl_int_divexact(a, left->row[row][row], a);
1059
346
      isl_int_neg(b, b);
1060
346
      isl_seq_combine(left->row[i] + i,
1061
346
          a, left->row[i] + i,
1062
346
          b, left->row[row] + i,
1063
346
          left->n_col - i);
1064
346
      isl_seq_combine(right->row[i], a, right->row[i],
1065
346
          b, right->row[row], right->n_col);
1066
346
    }
1067
475k
  }
1068
159k
  isl_int_clear(b);
1069
159k
1070
159k
  isl_int_set(a, left->row[0][0]);
1071
475k
  for (row = 1; row < left->n_row; 
++row316k
)
1072
316k
    isl_int_lcm(a, a, left->row[row][row]);
1073
159k
  if (isl_int_is_zero(a)){
1074
0
    isl_int_clear(a);
1075
0
    isl_assert(left->ctx, 0, goto error);
1076
0
  }
1077
635k
  
for (row = 0; 159k
row < left->n_row;
++row475k
) {
1078
475k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
1079
475k
    if (isl_int_is_one(left->row[row][row]))
1080
475k
      
continue469k
;
1081
6.84k
    isl_seq_scale(right->row[row], right->row[row],
1082
6.84k
        left->row[row][row], right->n_col);
1083
6.84k
  }
1084
159k
  isl_int_clear(a);
1085
159k
1086
159k
  isl_mat_free(left);
1087
159k
  return right;
1088
0
error:
1089
0
  isl_mat_free(left);
1090
0
  isl_mat_free(right);
1091
0
  return NULL;
1092
159k
}
1093
1094
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
1095
58.4k
{
1096
58.4k
  int i;
1097
58.4k
1098
285k
  for (i = 0; i < mat->n_row; 
++i226k
)
1099
226k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
1100
58.4k
}
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
773k
  for (i = 0; i < mat->n_row; 
++i605k
) {
1110
605k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
1111
605k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
1112
605k
    isl_int_set(mat->row[i][dst], tmp);
1113
605k
  }
1114
168k
  isl_int_clear(tmp);
1115
168k
}
1116
1117
__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
1118
61.1k
{
1119
61.1k
  struct isl_mat *inv;
1120
61.1k
  int row;
1121
61.1k
  isl_int a, b;
1122
61.1k
1123
61.1k
  mat = isl_mat_cow(mat);
1124
61.1k
  if (!mat)
1125
0
    return NULL;
1126
61.1k
1127
61.1k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
1128
61.1k
  inv = isl_mat_cow(inv);
1129
61.1k
  if (!inv)
1130
0
    goto error;
1131
61.1k
1132
61.1k
  isl_int_init(a);
1133
61.1k
  isl_int_init(b);
1134
248k
  for (row = 0; row < mat->n_row; 
++row187k
) {
1135
187k
    int pivot, first, i, off;
1136
187k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
1137
187k
    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
187k
    pivot += row;
1143
187k
    if (pivot != row)
1144
30.2k
      exchange(mat, &inv, NULL, row, pivot, row);
1145
187k
    if (isl_int_is_neg(mat->row[row][row]))
1146
187k
      
oppose(mat, &inv, NULL, row, row)42.5k
;
1147
187k
    first = row+1;
1148
448k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
1149
448k
                mat->n_col-first)) != -1) {
1150
260k
      first += off;
1151
260k
      isl_int_fdiv_q(a, mat->row[row][first],
1152
260k
                mat->row[row][row]);
1153
260k
      subtract(mat, &inv, NULL, row, row, first, a);
1154
260k
      if (!isl_int_is_zero(mat->row[row][first]))
1155
260k
        
exchange(mat, &inv, NULL, row, row, first)196k
;
1156
64.7k
      else
1157
64.7k
        ++first;
1158
260k
    }
1159
462k
    for (i = 0; i < row; 
++i275k
) {
1160
275k
      if (isl_int_is_zero(mat->row[row][i]))
1161
275k
        
continue190k
;
1162
84.2k
      isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
1163
84.2k
      isl_int_divexact(b, mat->row[row][i], a);
1164
84.2k
      isl_int_divexact(a, mat->row[row][row], a);
1165
84.2k
      isl_int_neg(a, a);
1166
84.2k
      isl_mat_col_combine(mat, i, a, i, b, row);
1167
84.2k
      isl_mat_col_combine(inv, i, a, i, b, row);
1168
84.2k
    }
1169
187k
  }
1170
61.1k
  isl_int_clear(b);
1171
61.1k
1172
61.1k
  isl_int_set(a, mat->row[0][0]);
1173
187k
  for (row = 1; row < mat->n_row; 
++row126k
)
1174
126k
    isl_int_lcm(a, a, mat->row[row][row]);
1175
61.1k
  if (isl_int_is_zero(a)){
1176
0
    isl_int_clear(a);
1177
0
    goto error;
1178
0
  }
1179
248k
  
for (row = 0; 61.1k
row < mat->n_row;
++row187k
) {
1180
187k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1181
187k
    if (isl_int_is_one(mat->row[row][row]))
1182
187k
      
continue128k
;
1183
58.4k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1184
58.4k
  }
1185
61.1k
  isl_int_clear(a);
1186
61.1k
1187
61.1k
  isl_mat_free(mat);
1188
61.1k
1189
61.1k
  return inv;
1190
0
error:
1191
0
  isl_mat_free(mat);
1192
0
  isl_mat_free(inv);
1193
0
  return NULL;
1194
61.1k
}
1195
1196
__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1197
3.77k
{
1198
3.77k
  struct isl_mat *transpose = NULL;
1199
3.77k
  int i, j;
1200
3.77k
1201
3.77k
  if (!mat)
1202
0
    return NULL;
1203
3.77k
1204
3.77k
  if (mat->n_col == mat->n_row) {
1205
3.77k
    mat = isl_mat_cow(mat);
1206
3.77k
    if (!mat)
1207
0
      return NULL;
1208
19.5k
    
for (i = 0; 3.77k
i < mat->n_row;
++i15.7k
)
1209
55.9k
      
for (j = i + 1; 15.7k
j < mat->n_col;
++j40.1k
)
1210
40.1k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1211
3.77k
    return mat;
1212
3.77k
  }
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
644k
{
1229
644k
  int r;
1230
644k
1231
644k
  mat = isl_mat_cow(mat);
1232
644k
  if (check_col_range(mat, i, 1) < 0 ||
1233
644k
      check_col_range(mat, j, 1) < 0)
1234
0
    return isl_mat_free(mat);
1235
644k
1236
19.4M
  
for (r = 0; 644k
r < mat->n_row;
++r18.7M
)
1237
18.7M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1238
644k
  return mat;
1239
644k
}
1240
1241
__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1242
  unsigned i, unsigned j)
1243
2.06M
{
1244
2.06M
  isl_int *t;
1245
2.06M
1246
2.06M
  if (!mat)
1247
0
    return NULL;
1248
2.06M
  mat = isl_mat_cow(mat);
1249
2.06M
  if (check_row_range(mat, i, 1) < 0 ||
1250
2.06M
      check_row_range(mat, j, 1) < 0)
1251
0
    return isl_mat_free(mat);
1252
2.06M
1253
2.06M
  t = mat->row[i];
1254
2.06M
  mat->row[i] = mat->row[j];
1255
2.06M
  mat->row[j] = t;
1256
2.06M
  return mat;
1257
2.06M
}
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.03M
{
1267
2.03M
  int i, j, k;
1268
2.03M
  struct isl_mat *prod;
1269
2.03M
1270
2.03M
  if (!left || !right)
1271
0
    goto error;
1272
2.03M
  isl_assert(left->ctx, left->n_col == right->n_row, goto error);
1273
2.03M
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1274
2.03M
  if (!prod)
1275
0
    goto error;
1276
2.03M
  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.86M
  
for (i = 0; 2.03M
i < prod->n_row;
++i4.83M
) {
1284
40.6M
    for (j = 0; j < prod->n_col; 
++j35.8M
)
1285
35.8M
      isl_int_mul(prod->row[i][j],
1286
4.83M
            left->row[i][0], right->row[0][j]);
1287
46.5M
    for (k = 1; k < left->n_col; 
++k41.7M
) {
1288
41.7M
      if (isl_int_is_zero(left->row[i][k]))
1289
41.7M
        
continue36.3M
;
1290
49.8M
      
for (j = 0; 5.33M
j < prod->n_col;
++j44.5M
)
1291
44.5M
        isl_int_addmul(prod->row[i][j],
1292
5.33M
              left->row[i][k], right->row[k][j]);
1293
5.33M
    }
1294
4.83M
  }
1295
2.03M
  isl_mat_free(left);
1296
2.03M
  isl_mat_free(right);
1297
2.03M
  return prod;
1298
0
error:
1299
0
  isl_mat_free(left);
1300
0
  isl_mat_free(right);
1301
0
  return NULL;
1302
2.03M
}
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.71M
{
1318
1.71M
  int i;
1319
1.71M
  struct isl_mat *t;
1320
1.71M
  int e;
1321
1.71M
1322
1.71M
  if (mat->n_col >= mat->n_row)
1323
1.07M
    e = 0;
1324
644k
  else
1325
644k
    e = mat->n_row - mat->n_col;
1326
1.71M
  if (has_div)
1327
572k
    
for (i = 0; 572k
i < n;
++i141
)
1328
572k
      
isl_int_mul141
(q[i][0], q[i][0], mat->row[0][0]);
1329
1.71M
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1330
1.71M
  t = isl_mat_product(t, mat);
1331
1.71M
  if (!t)
1332
0
    return -1;
1333
4.80M
  
for (i = 0; 1.71M
i < n;
++i3.09M
) {
1334
3.09M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1335
3.09M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1336
3.09M
          q[i] + has_div + t->n_col + e, n_div);
1337
3.09M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1338
3.09M
  }
1339
1.71M
  isl_mat_free(t);
1340
1.71M
  return 0;
1341
1.71M
}
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
572k
{
1356
572k
  struct isl_ctx *ctx;
1357
572k
1358
572k
  if (!bset || !mat)
1359
0
    goto error;
1360
572k
1361
572k
  ctx = bset->ctx;
1362
572k
  bset = isl_basic_set_cow(bset);
1363
572k
  if (!bset)
1364
0
    goto error;
1365
572k
1366
572k
  isl_assert(ctx, bset->dim->nparam == 0, goto error);
1367
572k
  isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1368
572k
  isl_assert(ctx, mat->n_col > 0, goto error);
1369
572k
1370
572k
  if (mat->n_col > mat->n_row) {
1371
104k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1372
104k
    if (!bset)
1373
0
      goto error;
1374
468k
  } 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
572k
1381
572k
  if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1382
572k
      isl_mat_copy(mat)) < 0)
1383
0
    goto error;
1384
572k
1385
572k
  if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1386
572k
      isl_mat_copy(mat)) < 0)
1387
0
    goto error;
1388
572k
1389
572k
  if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
1390
0
    goto error2;
1391
572k
1392
572k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1393
572k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1394
572k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1395
572k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1396
572k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1397
572k
1398
572k
  bset = isl_basic_set_simplify(bset);
1399
572k
  bset = isl_basic_set_finalize(bset);
1400
572k
1401
572k
  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.1k
{
1412
32.1k
  int i;
1413
32.1k
1414
32.1k
  set = isl_set_cow(set);
1415
32.1k
  if (!set)
1416
0
    goto error;
1417
32.1k
1418
91.2k
  
for (i = 0; 32.1k
i < set->n;
++i59.1k
) {
1419
59.1k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1420
59.1k
                isl_mat_copy(mat));
1421
59.1k
    if (!set->p[i])
1422
0
      goto error;
1423
59.1k
  }
1424
32.1k
  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.1k
  isl_mat_free(mat);
1432
32.1k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1433
32.1k
  return set;
1434
0
error:
1435
0
  isl_set_free(set);
1436
0
  isl_mat_free(mat);
1437
0
  return NULL;
1438
32.1k
}
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
87.6k
{
1501
87.6k
  int r;
1502
87.6k
1503
87.6k
  if (n == 0)
1504
22
    return mat;
1505
87.6k
1506
87.6k
  mat = isl_mat_cow(mat);
1507
87.6k
  if (check_col_range(mat, col, n) < 0)
1508
0
    return isl_mat_free(mat);
1509
87.6k
1510
87.6k
  if (col != mat->n_col-n) {
1511
66.4k
    for (r = 0; r < mat->n_row; 
++r51.9k
)
1512
51.9k
      isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1513
51.9k
          mat->n_col - col - n);
1514
14.4k
  }
1515
87.6k
  mat->n_col -= n;
1516
87.6k
  return mat;
1517
87.6k
}
1518
1519
__isl_give isl_mat *isl_mat_drop_rows(__isl_take isl_mat *mat,
1520
  unsigned row, unsigned n)
1521
118k
{
1522
118k
  int r;
1523
118k
1524
118k
  mat = isl_mat_cow(mat);
1525
118k
  if (check_row_range(mat, row, n) < 0)
1526
0
    return isl_mat_free(mat);
1527
118k
1528
345k
  
for (r = row; 118k
r+n < mat->n_row;
++r227k
)
1529
227k
    mat->row[r] = mat->row[r+n];
1530
118k
1531
118k
  mat->n_row -= n;
1532
118k
  return mat;
1533
118k
}
1534
1535
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1536
        unsigned col, unsigned n)
1537
85.7k
{
1538
85.7k
  isl_mat *ext;
1539
85.7k
1540
85.7k
  if (check_col_range(mat, col, 0) < 0)
1541
0
    return isl_mat_free(mat);
1542
85.7k
  if (n == 0)
1543
448
    return mat;
1544
85.2k
1545
85.2k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1546
85.2k
  if (!ext)
1547
0
    goto error;
1548
85.2k
1549
85.2k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1550
85.2k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1551
85.2k
        col + n, col, mat->n_col - col);
1552
85.2k
1553
85.2k
  isl_mat_free(mat);
1554
85.2k
  return ext;
1555
0
error:
1556
0
  isl_mat_free(mat);
1557
0
  return NULL;
1558
85.2k
}
1559
1560
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1561
  unsigned first, unsigned n)
1562
85.7k
{
1563
85.7k
  int i;
1564
85.7k
1565
85.7k
  if (!mat)
1566
0
    return NULL;
1567
85.7k
  mat = isl_mat_insert_cols(mat, first, n);
1568
85.7k
  if (!mat)
1569
0
    return NULL;
1570
85.7k
1571
89.6k
  
for (i = 0; 85.7k
i < mat->n_row;
++i3.92k
)
1572
3.92k
    isl_seq_clr(mat->row[i] + first, n);
1573
85.7k
1574
85.7k
  return mat;
1575
85.7k
}
1576
1577
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1578
4.62k
{
1579
4.62k
  if (!mat)
1580
0
    return NULL;
1581
4.62k
1582
4.62k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1583
4.62k
}
1584
1585
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1586
        unsigned row, unsigned n)
1587
5.07k
{
1588
5.07k
  isl_mat *ext;
1589
5.07k
1590
5.07k
  if (check_row_range(mat, row, 0) < 0)
1591
0
    return isl_mat_free(mat);
1592
5.07k
  if (n == 0)
1593
444
    return mat;
1594
4.62k
1595
4.62k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1596
4.62k
  if (!ext)
1597
0
    goto error;
1598
4.62k
1599
4.62k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1600
4.62k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1601
4.62k
        mat->n_row - row, 0, 0, mat->n_col);
1602
4.62k
1603
4.62k
  isl_mat_free(mat);
1604
4.62k
  return ext;
1605
0
error:
1606
0
  isl_mat_free(mat);
1607
0
  return NULL;
1608
4.62k
}
1609
1610
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1611
5.06k
{
1612
5.06k
  if (!mat)
1613
0
    return NULL;
1614
5.06k
1615
5.06k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1616
5.06k
}
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
5.31k
{
1644
5.31k
  int i;
1645
5.31k
1646
33.7k
  for (i = 0; i < mat->n_row; 
++i28.4k
)
1647
28.4k
    isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1648
5.31k
}
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
11.1k
{
1664
11.1k
  int i;
1665
11.1k
1666
49.9k
  for (i = 0; i < mat->n_row; 
++i38.7k
)
1667
38.7k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1668
11.1k
}
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
18.3k
{
1731
18.3k
  int r;
1732
18.3k
  struct isl_mat *H = NULL, *Q = NULL;
1733
18.3k
1734
18.3k
  if (!M)
1735
0
    return NULL;
1736
18.3k
1737
18.3k
  isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1738
18.3k
  M->n_row = row;
1739
18.3k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1740
18.3k
  M->n_row = M->n_col;
1741
18.3k
  if (!H)
1742
0
    goto error;
1743
36.7k
  
for (r = 0; 18.3k
r < row;
++r18.3k
)
1744
18.3k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1745
67.4k
  
for (r = row; 18.3k
r < M->n_row;
++r49.0k
)
1746
49.0k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1747
18.3k
  isl_mat_free(H);
1748
18.3k
  isl_mat_free(Q);
1749
18.3k
  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
18.3k
}
1756
1757
__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1758
  __isl_take isl_mat *bot)
1759
432
{
1760
432
  struct isl_mat *mat;
1761
432
1762
432
  if (!top || !bot)
1763
0
    goto error;
1764
432
1765
432
  isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1766
432
  if (top->n_row == 0) {
1767
392
    isl_mat_free(top);
1768
392
    return bot;
1769
392
  }
1770
40
  if (bot->n_row == 0) {
1771
0
    isl_mat_free(bot);
1772
0
    return top;
1773
0
  }
1774
40
1775
40
  mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1776
40
  if (!mat)
1777
0
    goto error;
1778
40
  isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1779
40
       0, 0, mat->n_col);
1780
40
  isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1781
40
       0, 0, mat->n_col);
1782
40
  isl_mat_free(top);
1783
40
  isl_mat_free(bot);
1784
40
  return mat;
1785
0
error:
1786
0
  isl_mat_free(top);
1787
0
  isl_mat_free(bot);
1788
0
  return NULL;
1789
40
}
1790
1791
isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1792
43.7k
{
1793
43.7k
  int i;
1794
43.7k
1795
43.7k
  if (!mat1 || !mat2)
1796
0
    return isl_bool_error;
1797
43.7k
1798
43.7k
  if (mat1->n_row != mat2->n_row)
1799
6.93k
    return isl_bool_false;
1800
36.7k
1801
36.7k
  if (mat1->n_col != mat2->n_col)
1802
0
    return isl_bool_false;
1803
36.7k
1804
71.2k
  
for (i = 0; 36.7k
i < mat1->n_row;
++i34.5k
)
1805
34.6k
    if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))
1806
189
      return isl_bool_false;
1807
36.7k
1808
36.7k
  
return isl_bool_true36.6k
;
1809
36.7k
}
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
604
{
1859
604
  isl_mat *res;
1860
604
1861
604
  if (!mat)
1862
0
    return NULL;
1863
604
  if (n == 0 || dst_col == src_col)
1864
472
    return mat;
1865
132
1866
132
  res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1867
132
  if (!res)
1868
0
    goto error;
1869
132
1870
132
  if (dst_col < src_col) {
1871
132
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1872
132
         0, 0, dst_col);
1873
132
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1874
132
         dst_col, src_col, n);
1875
132
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1876
132
         dst_col + n, dst_col, src_col - dst_col);
1877
132
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1878
132
         src_col + n, src_col + n,
1879
132
         res->n_col - src_col - n);
1880
132
  } 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
132
  isl_mat_free(mat);
1892
132
1893
132
  return res;
1894
0
error:
1895
0
  isl_mat_free(mat);
1896
0
  return NULL;
1897
132
}
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
732
{
1914
732
  int i;
1915
732
  isl_int g;
1916
732
1917
732
  isl_int_set_si(*gcd, 0);
1918
732
  if (!mat)
1919
0
    return;
1920
732
1921
732
  isl_int_init(g);
1922
4.65k
  for (i = 0; i < mat->n_row; 
++i3.92k
) {
1923
3.92k
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1924
3.92k
    isl_int_gcd(*gcd, *gcd, g);
1925
3.92k
  }
1926
732
  isl_int_clear(g);
1927
732
}
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
732
{
1950
732
  int i;
1951
732
1952
732
  if (isl_int_is_one(m))
1953
732
    
return mat0
;
1954
732
1955
732
  mat = isl_mat_cow(mat);
1956
732
  if (!mat)
1957
0
    return NULL;
1958
732
1959
4.65k
  
for (i = 0; 732
i < mat->n_row;
++i3.92k
)
1960
3.92k
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1961
732
1962
732
  return mat;
1963
732
}
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
732
{
1982
732
  isl_int gcd;
1983
732
1984
732
  if (!mat)
1985
0
    return NULL;
1986
732
1987
732
  isl_int_init(gcd);
1988
732
  isl_mat_gcd(mat, &gcd);
1989
732
  mat = isl_mat_scale_down(mat, gcd);
1990
732
  isl_int_clear(gcd);
1991
732
1992
732
  return mat;
1993
732
}
1994
1995
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1996
1.89k
{
1997
1.89k
  mat = isl_mat_cow(mat);
1998
1.89k
  if (!mat)
1999
0
    return NULL;
2000
1.89k
2001
1.89k
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
2002
1.89k
2003
1.89k
  return mat;
2004
1.89k
}
2005
2006
/* Number of initial non-zero columns.
2007
 */
2008
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
2009
1.28k
{
2010
1.28k
  int i;
2011
1.28k
2012
1.28k
  if (!mat)
2013
0
    return -1;
2014
1.28k
2015
2.12k
  
for (i = 0; 1.28k
i < mat->n_col;
++i833
)
2016
1.66k
    if (row_first_non_zero(mat->row, mat->n_row, i) < 0)
2017
832
      break;
2018
1.28k
2019
1.28k
  return i;
2020
1.28k
}
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
}