Coverage Report

Created: 2019-04-21 11:35

/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
1.24M
{
27
1.24M
  return mat ? mat->ctx : NULL;
28
1.24M
}
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
2.08M
{
56
2.08M
  int i;
57
2.08M
  struct isl_mat *mat;
58
2.08M
59
2.08M
  mat = isl_alloc_type(ctx, struct isl_mat);
60
2.08M
  if (!mat)
61
0
    return NULL;
62
2.08M
63
2.08M
  mat->row = NULL;
64
2.08M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
65
2.08M
  if (isl_blk_is_error(mat->block))
66
0
    goto error;
67
2.08M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
68
2.08M
  if (n_row && 
!mat->row1.50M
)
69
0
    goto error;
70
2.08M
71
10.9M
  
for (i = 0; 2.08M
i < n_row;
++i8.90M
)
72
8.90M
    mat->row[i] = mat->block.data + i * n_col;
73
2.08M
74
2.08M
  mat->ctx = ctx;
75
2.08M
  isl_ctx_ref(ctx);
76
2.08M
  mat->ref = 1;
77
2.08M
  mat->n_row = n_row;
78
2.08M
  mat->n_col = n_col;
79
2.08M
  mat->max_col = n_col;
80
2.08M
  mat->flags = 0;
81
2.08M
82
2.08M
  return mat;
83
0
error:
84
0
  isl_blk_free(ctx, mat->block);
85
0
  free(mat);
86
0
  return NULL;
87
2.08M
}
88
89
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
90
  unsigned n_row, unsigned n_col)
91
30.9k
{
92
30.9k
  int i;
93
30.9k
  isl_int *old;
94
30.9k
  isl_int **row;
95
30.9k
96
30.9k
  if (!mat)
97
0
    return NULL;
98
30.9k
99
30.9k
  if (mat->max_col >= n_col && 
mat->n_row >= n_row29.9k
) {
100
2.20k
    if (mat->n_col < n_col)
101
9
      mat->n_col = n_col;
102
2.20k
    return mat;
103
2.20k
  }
104
28.7k
105
28.7k
  if (mat->max_col < n_col) {
106
1.00k
    struct isl_mat *new_mat;
107
1.00k
108
1.00k
    if (n_row < mat->n_row)
109
2
      n_row = mat->n_row;
110
1.00k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
111
1.00k
    if (!new_mat)
112
0
      goto error;
113
6.78k
    
for (i = 0; 1.00k
i < mat->n_row;
++i5.78k
)
114
5.78k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
115
1.00k
    isl_mat_free(mat);
116
1.00k
    return new_mat;
117
1.00k
  }
118
27.6k
119
27.6k
  mat = isl_mat_cow(mat);
120
27.6k
  if (!mat)
121
0
    goto error;
122
27.6k
123
27.6k
  old = mat->block.data;
124
27.6k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
125
27.6k
  if (isl_blk_is_error(mat->block))
126
0
    goto error;
127
27.6k
  row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
128
27.6k
  if (n_row && !row)
129
0
    goto error;
130
27.6k
  mat->row = row;
131
27.6k
132
267k
  for (i = 0; i < mat->n_row; 
++i239k
)
133
239k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
134
90.8k
  for (i = mat->n_row; i < n_row; 
++i63.1k
)
135
63.1k
    mat->row[i] = mat->block.data + i * mat->max_col;
136
27.6k
  mat->n_row = n_row;
137
27.6k
  if (mat->n_col < n_col)
138
0
    mat->n_col = n_col;
139
27.6k
140
27.6k
  return mat;
141
0
error:
142
0
  isl_mat_free(mat);
143
0
  return NULL;
144
27.6k
}
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
1.01M
{
149
1.01M
  int i;
150
1.01M
  struct isl_mat *mat;
151
1.01M
152
1.01M
  mat = isl_alloc_type(ctx, struct isl_mat);
153
1.01M
  if (!mat)
154
0
    return NULL;
155
1.01M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
156
1.01M
  if (n_row && 
!mat->row551k
)
157
0
    goto error;
158
3.23M
  
for (i = 0; 1.01M
i < n_row;
++i2.21M
)
159
2.21M
    mat->row[i] = row[first_row+i] + first_col;
160
1.01M
  mat->ctx = ctx;
161
1.01M
  isl_ctx_ref(ctx);
162
1.01M
  mat->ref = 1;
163
1.01M
  mat->n_row = n_row;
164
1.01M
  mat->n_col = n_col;
165
1.01M
  mat->block = isl_blk_empty();
166
1.01M
  mat->flags = ISL_MAT_BORROWED;
167
1.01M
  return mat;
168
0
error:
169
0
  free(mat);
170
0
  return NULL;
171
1.01M
}
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
188k
{
176
188k
  if (!mat)
177
0
    return NULL;
178
188k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
179
188k
          first_col, n_col);
180
188k
}
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
126k
{
185
126k
  int i;
186
126k
187
390k
  for (i = 0; i < n_row; 
++i264k
)
188
264k
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
189
126k
}
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
43.7k
{
194
43.7k
  int i;
195
43.7k
196
131k
  for (i = 0; i < n_row; 
++i87.3k
)
197
87.3k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
198
43.7k
}
199
200
__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
201
808k
{
202
808k
  if (!mat)
203
0
    return NULL;
204
808k
205
808k
  mat->ref++;
206
808k
  return mat;
207
808k
}
208
209
__isl_give isl_mat *isl_mat_dup(__isl_keep isl_mat *mat)
210
124k
{
211
124k
  int i;
212
124k
  struct isl_mat *mat2;
213
124k
214
124k
  if (!mat)
215
783
    return NULL;
216
123k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
217
123k
  if (!mat2)
218
0
    return NULL;
219
398k
  
for (i = 0; 123k
i < mat->n_row;
++i274k
)
220
274k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
221
123k
  return mat2;
222
123k
}
223
224
__isl_give isl_mat *isl_mat_cow(__isl_take isl_mat *mat)
225
1.30M
{
226
1.30M
  struct isl_mat *mat2;
227
1.30M
  if (!mat)
228
0
    return NULL;
229
1.30M
230
1.30M
  if (mat->ref == 1 && 
!1.28M
ISL_F_ISSET1.28M
(mat, ISL_MAT_BORROWED))
231
1.30M
    
return mat1.18M
;
232
122k
233
122k
  mat2 = isl_mat_dup(mat);
234
122k
  isl_mat_free(mat);
235
122k
  return mat2;
236
122k
}
237
238
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
239
4.45M
{
240
4.45M
  if (!mat)
241
546k
    return NULL;
242
3.91M
243
3.91M
  if (--mat->ref > 0)
244
808k
    return NULL;
245
3.10M
246
3.10M
  if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
247
3.10M
    
isl_blk_free(mat->ctx, mat->block)2.08M
;
248
3.10M
  isl_ctx_deref(mat->ctx);
249
3.10M
  free(mat->row);
250
3.10M
  free(mat);
251
3.10M
252
3.10M
  return NULL;
253
3.10M
}
254
255
int isl_mat_rows(__isl_keep isl_mat *mat)
256
127k
{
257
127k
  return mat ? mat->n_row : 
-10
;
258
127k
}
259
260
int isl_mat_cols(__isl_keep isl_mat *mat)
261
16.2k
{
262
16.2k
  return mat ? mat->n_col : 
-10
;
263
16.2k
}
264
265
/* Check that "col" is a valid column position for "mat".
266
 */
267
static isl_stat check_col(__isl_keep isl_mat *mat, int col)
268
4.01k
{
269
4.01k
  if (!mat)
270
0
    return isl_stat_error;
271
4.01k
  if (col < 0 || col >= mat->n_col)
272
4.01k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
273
4.01k
      "column out of range", return isl_stat_error);
274
4.01k
  return isl_stat_ok;
275
4.01k
}
276
277
/* Check that "row" is a valid row position for "mat".
278
 */
279
static isl_stat check_row(__isl_keep isl_mat *mat, int row)
280
4.02k
{
281
4.02k
  if (!mat)
282
0
    return isl_stat_error;
283
4.02k
  if (row < 0 || row >= mat->n_row)
284
4.02k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
285
4.02k
      "row out of range", return isl_stat_error);
286
4.02k
  return isl_stat_ok;
287
4.02k
}
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
559k
{
294
559k
  if (!mat)
295
0
    return isl_stat_error;
296
559k
  if (first + n > mat->n_col || first + n < first)
297
559k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
298
559k
      "column position or range out of bounds",
299
559k
      return isl_stat_error);
300
559k
  return isl_stat_ok;
301
559k
}
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
1.27M
{
308
1.27M
  if (!mat)
309
0
    return isl_stat_error;
310
1.27M
  if (first + n > mat->n_row || first + n < first)
311
1.27M
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
312
1.27M
      "row position or range out of bounds",
313
1.27M
      return isl_stat_error);
314
1.27M
  return isl_stat_ok;
315
1.27M
}
316
317
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
318
2.79k
{
319
2.79k
  if (check_row(mat, row) < 0)
320
0
    return -1;
321
2.79k
  if (check_col(mat, col) < 0)
322
0
    return -1;
323
2.79k
  isl_int_set(*v, mat->row[row][col]);
324
2.79k
  return 0;
325
2.79k
}
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.13k
{
345
1.13k
  mat = isl_mat_cow(mat);
346
1.13k
  if (check_row(mat, row) < 0)
347
0
    return isl_mat_free(mat);
348
1.13k
  if (check_col(mat, col) < 0)
349
0
    return isl_mat_free(mat);
350
1.13k
  isl_int_set(mat->row[row][col], v);
351
1.13k
  return mat;
352
1.13k
}
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
315k
{
386
315k
  int i;
387
315k
  struct isl_mat *mat;
388
315k
389
315k
  mat = isl_mat_alloc(ctx, n_row, n_row);
390
315k
  if (!mat)
391
0
    return NULL;
392
1.71M
  
for (i = 0; 315k
i < n_row;
++i1.40M
) {
393
1.40M
    isl_seq_clr(mat->row[i], i);
394
1.40M
    isl_int_set(mat->row[i][i], d);
395
1.40M
    isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
396
1.40M
  }
397
315k
398
315k
  return mat;
399
315k
}
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
315k
{
419
315k
  if (!ctx)
420
0
    return NULL;
421
315k
  return isl_mat_diag(ctx, n_row, ctx->one);
422
315k
}
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
104k
{
451
104k
  int i;
452
104k
  struct isl_vec *prod;
453
104k
454
104k
  if (!mat || !vec)
455
0
    goto error;
456
104k
457
104k
  isl_assert(mat->ctx, mat->n_col == vec->size, goto error);
458
104k
459
104k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
460
104k
  if (!prod)
461
0
    goto error;
462
104k
463
835k
  
for (i = 0; 104k
i < prod->size;
++i730k
)
464
730k
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
465
730k
          &prod->block.data[i]);
466
104k
  isl_mat_free(mat);
467
104k
  isl_vec_free(vec);
468
104k
  return prod;
469
0
error:
470
0
  isl_mat_free(mat);
471
0
  isl_vec_free(vec);
472
0
  return NULL;
473
104k
}
474
475
__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
476
  __isl_take isl_vec *vec)
477
36
{
478
36
  struct isl_mat *vec_mat;
479
36
  int i;
480
36
481
36
  if (!mat || !vec)
482
0
    goto error;
483
36
  vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
484
36
  if (!vec_mat)
485
0
    goto error;
486
252
  
for (i = 0; 36
i < vec->size;
++i216
)
487
216
    isl_int_set(vec_mat->row[i][0], vec->el[i]);
488
36
  vec_mat = isl_mat_inverse_product(mat, vec_mat);
489
36
  isl_vec_free(vec);
490
36
  if (!vec_mat)
491
0
    return NULL;
492
36
  vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
493
36
  if (vec)
494
252
    
for (i = 0; 36
i < vec->size;
++i216
)
495
216
      isl_int_set(vec->el[i], vec_mat->row[i][0]);
496
36
  isl_mat_free(vec_mat);
497
36
  return vec;
498
0
error:
499
0
  isl_mat_free(mat);
500
0
  isl_vec_free(vec);
501
0
  return NULL;
502
36
}
503
504
__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,
505
  __isl_take isl_mat *mat)
506
5.67k
{
507
5.67k
  int i, j;
508
5.67k
  struct isl_vec *prod;
509
5.67k
510
5.67k
  if (!mat || !vec)
511
0
    goto error;
512
5.67k
513
5.67k
  isl_assert(mat->ctx, mat->n_row == vec->size, goto error);
514
5.67k
515
5.67k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
516
5.67k
  if (!prod)
517
0
    goto error;
518
5.67k
519
33.9k
  
for (i = 0; 5.67k
i < prod->size;
++i28.2k
) {
520
28.2k
    isl_int_set_si(prod->el[i], 0);
521
229k
    for (j = 0; j < vec->size; 
++j201k
)
522
201k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
523
28.2k
  }
524
5.67k
  isl_mat_free(mat);
525
5.67k
  isl_vec_free(vec);
526
5.67k
  return prod;
527
0
error:
528
0
  isl_mat_free(mat);
529
0
  isl_vec_free(vec);
530
0
  return NULL;
531
5.67k
}
532
533
__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
534
  __isl_take isl_mat *right)
535
43.7k
{
536
43.7k
  int i;
537
43.7k
  struct isl_mat *sum;
538
43.7k
539
43.7k
  if (!left || !right)
540
0
    goto error;
541
43.7k
542
43.7k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
543
43.7k
  isl_assert(left->ctx, left->n_row >= 1, goto error);
544
43.7k
  isl_assert(left->ctx, left->n_col >= 1, goto error);
545
43.7k
  isl_assert(left->ctx, right->n_col >= 1, goto error);
546
43.7k
  isl_assert(left->ctx,
547
43.7k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
548
43.7k
      goto error);
549
43.7k
  isl_assert(left->ctx,
550
43.7k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
551
43.7k
      goto error);
552
43.7k
553
43.7k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
554
43.7k
  if (!sum)
555
0
    goto error;
556
43.7k
  isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
557
43.7k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
558
43.7k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
559
43.7k
560
43.7k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
561
334k
  for (i = 1; i < sum->n_row; 
++i290k
) {
562
290k
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
563
290k
    isl_int_addmul(sum->row[i][0],
564
290k
        right->row[0][0], right->row[i][0]);
565
290k
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
566
290k
        left->n_col-1);
567
290k
    isl_seq_scale(sum->row[i]+left->n_col,
568
290k
        right->row[i]+1, right->row[0][0],
569
290k
        right->n_col-1);
570
290k
  }
571
43.7k
572
43.7k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
573
43.7k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
574
43.7k
  isl_mat_free(left);
575
43.7k
  isl_mat_free(right);
576
43.7k
  return sum;
577
0
error:
578
0
  isl_mat_free(left);
579
0
  isl_mat_free(right);
580
0
  return NULL;
581
43.7k
}
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
404k
{
586
404k
  int r;
587
2.45M
  for (r = row; r < M->n_row; 
++r2.05M
)
588
2.05M
    isl_int_swap(M->row[r][i], M->row[r][j]);
589
404k
  if (U) {
590
3.38M
    for (r = 0; r < (*U)->n_row; 
++r2.98M
)
591
2.98M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
592
397k
  }
593
404k
  if (Q)
594
89.5k
    isl_mat_swap_rows(*Q, i, j);
595
404k
}
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
388k
{
600
388k
  int r;
601
1.25M
  for (r = row; r < M->n_row; 
++r870k
)
602
870k
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
603
388k
  if (U) {
604
2.25M
    for (r = 0; r < (*U)->n_row; 
++r1.87M
)
605
1.87M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
606
377k
  }
607
388k
  if (Q) {
608
316k
    for (r = 0; r < (*Q)->n_col; 
++r266k
)
609
266k
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
610
49.7k
  }
611
388k
}
612
613
static void oppose(struct isl_mat *M, struct isl_mat **U,
614
  struct isl_mat **Q, unsigned row, unsigned col)
615
98.5k
{
616
98.5k
  int r;
617
360k
  for (r = row; r < M->n_row; 
++r261k
)
618
261k
    isl_int_neg(M->row[r][col], M->row[r][col]);
619
98.5k
  if (U) {
620
573k
    for (r = 0; r < (*U)->n_row; 
++r483k
)
621
483k
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
622
89.5k
  }
623
98.5k
  if (Q)
624
38.3k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
625
98.5k
}
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
142k
{
642
142k
  isl_int c;
643
142k
  int row, col;
644
142k
645
142k
  if (U)
646
133k
    *U = NULL;
647
142k
  if (Q)
648
64.8k
    *Q = NULL;
649
142k
  if (!M)
650
0
    goto error;
651
142k
  M = isl_mat_cow(M);
652
142k
  if (!M)
653
0
    goto error;
654
142k
  if (U) {
655
133k
    *U = isl_mat_identity(M->ctx, M->n_col);
656
133k
    if (!*U)
657
0
      goto error;
658
142k
  }
659
142k
  if (Q) {
660
64.8k
    *Q = isl_mat_identity(M->ctx, M->n_col);
661
64.8k
    if (!*Q)
662
0
      goto error;
663
142k
  }
664
142k
665
142k
  col = 0;
666
142k
  isl_int_init(c);
667
624k
  for (row = 0; row < M->n_row; 
++row481k
) {
668
481k
    int first, i, off;
669
481k
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
670
481k
    if (first == -1)
671
133k
      continue;
672
348k
    first += col;
673
348k
    if (first != col)
674
183k
      exchange(M, U, Q, row, first, col);
675
348k
    if (isl_int_is_neg(M->row[row][col]))
676
348k
      
oppose(M, U, Q, row, col)67.8k
;
677
348k
    first = col+1;
678
450k
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
679
450k
                   M->n_col-first)) != -1) {
680
102k
      first += off;
681
102k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
682
102k
      subtract(M, U, Q, row, col, first, c);
683
102k
      if (!isl_int_is_zero(M->row[row][first]))
684
102k
        
exchange(M, U, Q, row, first, col)5.73k
;
685
96.8k
      else
686
96.8k
        ++first;
687
102k
    }
688
1.50M
    for (i = 0; i < col; 
++i1.15M
) {
689
1.15M
      if (isl_int_is_zero(M->row[row][i]))
690
1.15M
        
continue1.11M
;
691
40.8k
      if (neg)
692
40.8k
        
isl_int_cdiv_q0
(c, M->row[row][i], M->row[row][col]);
693
40.8k
      else
694
40.8k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
695
40.8k
      if (isl_int_is_zero(c))
696
40.8k
        
continue9.96k
;
697
30.8k
      subtract(M, U, Q, row, col, i, c);
698
30.8k
    }
699
348k
    ++col;
700
348k
  }
701
142k
  isl_int_clear(c);
702
142k
703
142k
  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
142k
}
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
915
{
721
915
  int k, nr, nc;
722
915
  isl_ctx *ctx;
723
915
724
915
  if (!mat)
725
0
    return NULL;
726
915
727
915
  ctx = isl_mat_get_ctx(mat);
728
915
  nr = isl_mat_rows(mat);
729
915
  nc = isl_mat_cols(mat);
730
915
731
2.43k
  for (k = 0; k < nr; 
++k1.51k
) {
732
1.51k
    if (k == row)
733
915
      continue;
734
600
    if (isl_int_is_zero(mat->row[k][col]))
735
600
      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
915
743
915
  return mat;
744
915
}
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.07k
{
757
1.07k
  int k, row, last, nr, nc;
758
1.07k
759
1.07k
  if (!mat)
760
0
    return NULL;
761
1.07k
762
1.07k
  nr = isl_mat_rows(mat);
763
1.07k
  nc = isl_mat_cols(mat);
764
1.07k
765
1.07k
  last = nc - 1;
766
1.98k
  for (row = nr - 1; row >= 0; 
--row915
) {
767
1.17k
    for (; last >= 0; 
--last264
) {
768
1.48k
      for (k = row; k >= 0; 
--k308
)
769
1.22k
        if (!isl_int_is_zero(mat->row[k][last]))
770
1.22k
          
break915
;
771
1.17k
      if (k >= 0)
772
915
        break;
773
1.17k
    }
774
915
    if (last < 0)
775
0
      break;
776
915
    if (k != row)
777
0
      mat = isl_mat_swap_rows(mat, k, row);
778
915
    if (!mat)
779
0
      return NULL;
780
915
    if (isl_int_is_neg(mat->row[row][last]))
781
915
      
mat = isl_mat_row_neg(mat, row)4
;
782
915
    mat = eliminate(mat, row, last);
783
915
    if (!mat)
784
0
      return NULL;
785
915
  }
786
1.07k
  mat = isl_mat_drop_rows(mat, 0, row + 1);
787
1.07k
788
1.07k
  return mat;
789
1.07k
}
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.07k
{
796
1.07k
  int i, nr, nc;
797
1.07k
798
1.07k
  if (!mat)
799
0
    return NULL;
800
1.07k
801
1.07k
  nr = isl_mat_rows(mat);
802
1.07k
  nc = isl_mat_cols(mat);
803
1.07k
804
1.98k
  for (i = 0; i < nr; 
++i915
) {
805
915
    int pos;
806
915
807
915
    pos = isl_seq_first_non_zero(mat->row[i], nc);
808
915
    if (pos < 0)
809
0
      continue;
810
915
    if (isl_int_is_nonneg(mat->row[i][pos]))
811
915
      
continue905
;
812
10
    mat = isl_mat_row_neg(mat, i);
813
10
    if (!mat)
814
0
      return NULL;
815
10
  }
816
1.07k
817
1.07k
  return mat;
818
1.07k
}
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
211k
{
888
211k
  int i;
889
211k
  struct isl_mat *mat2;
890
211k
891
211k
  if (!mat)
892
0
    return NULL;
893
211k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
894
211k
  if (!mat2)
895
0
    goto error;
896
211k
  isl_int_set_si(mat2->row[0][0], 1);
897
211k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
898
1.31M
  for (i = 0; i < mat->n_row; 
++i1.09M
) {
899
1.09M
    isl_int_set_si(mat2->row[1+i][0], 0);
900
1.09M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
901
1.09M
  }
902
211k
  isl_mat_free(mat);
903
211k
  return mat2;
904
0
error:
905
0
  isl_mat_free(mat);
906
0
  return NULL;
907
211k
}
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
37.4k
{
917
37.4k
  int i;
918
37.4k
  isl_mat *mat;
919
37.4k
920
37.4k
  if (!mat1 || !mat2)
921
0
    goto error;
922
37.4k
923
37.4k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
924
37.4k
               mat1->n_col + mat2->n_col);
925
37.4k
  if (!mat)
926
0
    goto error;
927
74.9k
  
for (i = 0; 37.4k
i < mat1->n_row;
++i37.4k
) {
928
37.4k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
929
37.4k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
930
37.4k
  }
931
249k
  for (i = 0; i < mat2->n_row; 
++i212k
) {
932
212k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
933
212k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
934
212k
                mat2->row[i], mat2->n_col);
935
212k
  }
936
37.4k
  isl_mat_free(mat1);
937
37.4k
  isl_mat_free(mat2);
938
37.4k
  return mat;
939
0
error:
940
0
  isl_mat_free(mat1);
941
0
  isl_mat_free(mat2);
942
0
  return NULL;
943
37.4k
}
944
945
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
946
282k
{
947
282k
  int i;
948
282k
949
473k
  for (i = 0; i < n_row; 
++i190k
)
950
332k
    if (!isl_int_is_zero(row[i][col]))
951
332k
      
return i141k
;
952
282k
  
return -1141k
;
953
282k
}
954
955
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
956
140k
{
957
140k
  int i, min = row_first_non_zero(row, n_row, col);
958
140k
  if (min < 0)
959
0
    return -1;
960
331k
  
for (i = min + 1; 140k
i < n_row;
++i190k
) {
961
190k
    if (isl_int_is_zero(row[i][col]))
962
190k
      
continue190k
;
963
596
    if (isl_int_abs_lt(row[i][col], row[min][col]))
964
596
      
min = i10
;
965
596
  }
966
140k
  return min;
967
140k
}
968
969
static isl_stat inv_exchange(__isl_keep isl_mat **left,
970
  __isl_keep isl_mat **right, unsigned i, unsigned j)
971
59
{
972
59
  *left = isl_mat_swap_rows(*left, i, j);
973
59
  *right = isl_mat_swap_rows(*right, i, j);
974
59
975
59
  if (!*left || !*right)
976
0
    return isl_stat_error;
977
59
  return isl_stat_ok;
978
59
}
979
980
static void inv_oppose(
981
  struct isl_mat *left, struct isl_mat *right, unsigned row)
982
30
{
983
30
  isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
984
30
  isl_seq_neg(right->row[row], right->row[row], right->n_col);
985
30
}
986
987
static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
988
  unsigned row, unsigned i, isl_int m)
989
598
{
990
598
  isl_int_neg(m, m);
991
598
  isl_seq_combine(left->row[i]+row,
992
598
      left->ctx->one, left->row[i]+row,
993
598
      m, left->row[row]+row,
994
598
      left->n_col-row);
995
598
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
996
598
      m, right->row[row], right->n_col);
997
598
}
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
47.8k
{
1004
47.8k
  int row;
1005
47.8k
  isl_int a, b;
1006
47.8k
1007
47.8k
  if (!left || !right)
1008
0
    goto error;
1009
47.8k
1010
47.8k
  isl_assert(left->ctx, left->n_row == left->n_col, goto error);
1011
47.8k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
1012
47.8k
1013
47.8k
  if (left->n_row == 0) {
1014
0
    isl_mat_free(left);
1015
0
    return right;
1016
0
  }
1017
47.8k
1018
47.8k
  left = isl_mat_cow(left);
1019
47.8k
  right = isl_mat_cow(right);
1020
47.8k
  if (!left || !right)
1021
0
    goto error;
1022
47.8k
1023
47.8k
  isl_int_init(a);
1024
47.8k
  isl_int_init(b);
1025
188k
  for (row = 0; row < left->n_row; 
++row140k
) {
1026
140k
    int pivot, first, i, off;
1027
140k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
1028
140k
    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
140k
    pivot += row;
1034
140k
    if (pivot != row)
1035
57
      if (inv_exchange(&left, &right, pivot, row) < 0)
1036
0
        goto error;
1037
140k
    if (isl_int_is_neg(left->row[row][row]))
1038
140k
      
inv_oppose(left, right, row)30
;
1039
140k
    first = row+1;
1040
141k
    while ((off = row_first_non_zero(left->row+first,
1041
141k
          left->n_row-first, row)) != -1) {
1042
598
      first += off;
1043
598
      isl_int_fdiv_q(a, left->row[first][row],
1044
598
          left->row[row][row]);
1045
598
      inv_subtract(left, right, row, first, a);
1046
598
      if (!isl_int_is_zero(left->row[first][row])) {
1047
2
        if (inv_exchange(&left, &right, row, first) < 0)
1048
0
          goto error;
1049
596
      } else {
1050
596
        ++first;
1051
596
      }
1052
598
    }
1053
331k
    
for (i = 0; 140k
i < row;
++i190k
) {
1054
190k
      if (isl_int_is_zero(left->row[i][row]))
1055
190k
        
continue190k
;
1056
77
      isl_int_gcd(a, left->row[row][row], left->row[i][row]);
1057
77
      isl_int_divexact(b, left->row[i][row], a);
1058
77
      isl_int_divexact(a, left->row[row][row], a);
1059
77
      isl_int_neg(b, b);
1060
77
      isl_seq_combine(left->row[i] + i,
1061
77
          a, left->row[i] + i,
1062
77
          b, left->row[row] + i,
1063
77
          left->n_col - i);
1064
77
      isl_seq_combine(right->row[i], a, right->row[i],
1065
77
          b, right->row[row], right->n_col);
1066
77
    }
1067
140k
  }
1068
47.8k
  isl_int_clear(b);
1069
47.8k
1070
47.8k
  isl_int_set(a, left->row[0][0]);
1071
140k
  for (row = 1; row < left->n_row; 
++row92.5k
)
1072
92.5k
    isl_int_lcm(a, a, left->row[row][row]);
1073
47.8k
  if (isl_int_is_zero(a)){
1074
0
    isl_int_clear(a);
1075
0
    isl_assert(left->ctx, 0, goto error);
1076
0
  }
1077
188k
  
for (row = 0; 47.8k
row < left->n_row;
++row140k
) {
1078
140k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
1079
140k
    if (isl_int_is_one(left->row[row][row]))
1080
140k
      
continue138k
;
1081
1.45k
    isl_seq_scale(right->row[row], right->row[row],
1082
1.45k
        left->row[row][row], right->n_col);
1083
1.45k
  }
1084
47.8k
  isl_int_clear(a);
1085
47.8k
1086
47.8k
  isl_mat_free(left);
1087
47.8k
  return right;
1088
0
error:
1089
0
  isl_mat_free(left);
1090
0
  isl_mat_free(right);
1091
0
  return NULL;
1092
47.8k
}
1093
1094
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
1095
53.4k
{
1096
53.4k
  int i;
1097
53.4k
1098
248k
  for (i = 0; i < mat->n_row; 
++i195k
)
1099
195k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
1100
53.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
155k
{
1105
155k
  int i;
1106
155k
  isl_int tmp;
1107
155k
1108
155k
  isl_int_init(tmp);
1109
694k
  for (i = 0; i < mat->n_row; 
++i538k
) {
1110
538k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
1111
538k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
1112
538k
    isl_int_set(mat->row[i][dst], tmp);
1113
538k
  }
1114
155k
  isl_int_clear(tmp);
1115
155k
}
1116
1117
__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
1118
46.4k
{
1119
46.4k
  struct isl_mat *inv;
1120
46.4k
  int row;
1121
46.4k
  isl_int a, b;
1122
46.4k
1123
46.4k
  mat = isl_mat_cow(mat);
1124
46.4k
  if (!mat)
1125
0
    return NULL;
1126
46.4k
1127
46.4k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
1128
46.4k
  inv = isl_mat_cow(inv);
1129
46.4k
  if (!inv)
1130
0
    goto error;
1131
46.4k
1132
46.4k
  isl_int_init(a);
1133
46.4k
  isl_int_init(b);
1134
178k
  for (row = 0; row < mat->n_row; 
++row132k
) {
1135
132k
    int pivot, first, i, off;
1136
132k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
1137
132k
    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
132k
    pivot += row;
1143
132k
    if (pivot != row)
1144
20.0k
      exchange(mat, &inv, NULL, row, pivot, row);
1145
132k
    if (isl_int_is_neg(mat->row[row][row]))
1146
132k
      
oppose(mat, &inv, NULL, row, row)30.6k
;
1147
132k
    first = row+1;
1148
387k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
1149
387k
                mat->n_col-first)) != -1) {
1150
255k
      first += off;
1151
255k
      isl_int_fdiv_q(a, mat->row[row][first],
1152
255k
                mat->row[row][row]);
1153
255k
      subtract(mat, &inv, NULL, row, row, first, a);
1154
255k
      if (!isl_int_is_zero(mat->row[row][first]))
1155
255k
        
exchange(mat, &inv, NULL, row, row, first)195k
;
1156
59.2k
      else
1157
59.2k
        ++first;
1158
255k
    }
1159
274k
    for (i = 0; i < row; 
++i141k
) {
1160
141k
      if (isl_int_is_zero(mat->row[row][i]))
1161
141k
        
continue63.7k
;
1162
77.9k
      isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
1163
77.9k
      isl_int_divexact(b, mat->row[row][i], a);
1164
77.9k
      isl_int_divexact(a, mat->row[row][row], a);
1165
77.9k
      isl_int_neg(a, a);
1166
77.9k
      isl_mat_col_combine(mat, i, a, i, b, row);
1167
77.9k
      isl_mat_col_combine(inv, i, a, i, b, row);
1168
77.9k
    }
1169
132k
  }
1170
46.4k
  isl_int_clear(b);
1171
46.4k
1172
46.4k
  isl_int_set(a, mat->row[0][0]);
1173
132k
  for (row = 1; row < mat->n_row; 
++row85.9k
)
1174
85.9k
    isl_int_lcm(a, a, mat->row[row][row]);
1175
46.4k
  if (isl_int_is_zero(a)){
1176
0
    isl_int_clear(a);
1177
0
    goto error;
1178
0
  }
1179
178k
  
for (row = 0; 46.4k
row < mat->n_row;
++row132k
) {
1180
132k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1181
132k
    if (isl_int_is_one(mat->row[row][row]))
1182
132k
      
continue78.9k
;
1183
53.4k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1184
53.4k
  }
1185
46.4k
  isl_int_clear(a);
1186
46.4k
1187
46.4k
  isl_mat_free(mat);
1188
46.4k
1189
46.4k
  return inv;
1190
0
error:
1191
0
  isl_mat_free(mat);
1192
0
  isl_mat_free(inv);
1193
0
  return NULL;
1194
46.4k
}
1195
1196
__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1197
1.72k
{
1198
1.72k
  struct isl_mat *transpose = NULL;
1199
1.72k
  int i, j;
1200
1.72k
1201
1.72k
  if (!mat)
1202
0
    return NULL;
1203
1.72k
1204
1.72k
  if (mat->n_col == mat->n_row) {
1205
1.72k
    mat = isl_mat_cow(mat);
1206
1.72k
    if (!mat)
1207
0
      return NULL;
1208
6.59k
    
for (i = 0; 1.72k
i < mat->n_row;
++i4.87k
)
1209
14.3k
      
for (j = i + 1; 4.87k
j < mat->n_col;
++j9.45k
)
1210
9.45k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1211
1.72k
    return mat;
1212
1.72k
  }
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
247k
{
1229
247k
  int r;
1230
247k
1231
247k
  mat = isl_mat_cow(mat);
1232
247k
  if (check_col_range(mat, i, 1) < 0 ||
1233
247k
      check_col_range(mat, j, 1) < 0)
1234
0
    return isl_mat_free(mat);
1235
247k
1236
6.31M
  
for (r = 0; 247k
r < mat->n_row;
++r6.06M
)
1237
6.06M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1238
247k
  return mat;
1239
247k
}
1240
1241
__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1242
  unsigned i, unsigned j)
1243
615k
{
1244
615k
  isl_int *t;
1245
615k
1246
615k
  if (!mat)
1247
0
    return NULL;
1248
615k
  mat = isl_mat_cow(mat);
1249
615k
  if (check_row_range(mat, i, 1) < 0 ||
1250
615k
      check_row_range(mat, j, 1) < 0)
1251
0
    return isl_mat_free(mat);
1252
615k
1253
615k
  t = mat->row[i];
1254
615k
  mat->row[i] = mat->row[j];
1255
615k
  mat->row[j] = t;
1256
615k
  return mat;
1257
615k
}
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
772k
{
1267
772k
  int i, j, k;
1268
772k
  struct isl_mat *prod;
1269
772k
1270
772k
  if (!left || !right)
1271
0
    goto error;
1272
772k
  isl_assert(left->ctx, left->n_col == right->n_row, goto error);
1273
772k
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1274
772k
  if (!prod)
1275
0
    goto error;
1276
772k
  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
2.34M
  
for (i = 0; 772k
i < prod->n_row;
++i1.56M
) {
1284
12.9M
    for (j = 0; j < prod->n_col; 
++j11.3M
)
1285
11.3M
      isl_int_mul(prod->row[i][j],
1286
1.56M
            left->row[i][0], right->row[0][j]);
1287
14.6M
    for (k = 1; k < left->n_col; 
++k13.0M
) {
1288
13.0M
      if (isl_int_is_zero(left->row[i][k]))
1289
13.0M
        
continue11.2M
;
1290
15.9M
      
for (j = 0; 1.78M
j < prod->n_col;
++j14.1M
)
1291
14.1M
        isl_int_addmul(prod->row[i][j],
1292
1.78M
              left->row[i][k], right->row[k][j]);
1293
1.78M
    }
1294
1.56M
  }
1295
772k
  isl_mat_free(left);
1296
772k
  isl_mat_free(right);
1297
772k
  return prod;
1298
0
error:
1299
0
  isl_mat_free(left);
1300
0
  isl_mat_free(right);
1301
0
  return NULL;
1302
772k
}
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
680k
{
1318
680k
  int i;
1319
680k
  struct isl_mat *t;
1320
680k
  int e;
1321
680k
1322
680k
  if (mat->n_col >= mat->n_row)
1323
429k
    e = 0;
1324
250k
  else
1325
250k
    e = mat->n_row - mat->n_col;
1326
680k
  if (has_div)
1327
226k
    
for (i = 0; 226k
i < n;
++i27
)
1328
226k
      
isl_int_mul27
(q[i][0], q[i][0], mat->row[0][0]);
1329
680k
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1330
680k
  t = isl_mat_product(t, mat);
1331
680k
  if (!t)
1332
0
    return -1;
1333
1.73M
  
for (i = 0; 680k
i < n;
++i1.05M
) {
1334
1.05M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1335
1.05M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1336
1.05M
          q[i] + has_div + t->n_col + e, n_div);
1337
1.05M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1338
1.05M
  }
1339
680k
  isl_mat_free(t);
1340
680k
  return 0;
1341
680k
}
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
226k
{
1356
226k
  struct isl_ctx *ctx;
1357
226k
1358
226k
  if (!bset || !mat)
1359
0
    goto error;
1360
226k
1361
226k
  ctx = bset->ctx;
1362
226k
  bset = isl_basic_set_cow(bset);
1363
226k
  if (!bset)
1364
0
    goto error;
1365
226k
1366
226k
  isl_assert(ctx, bset->dim->nparam == 0, goto error);
1367
226k
  isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1368
226k
  isl_assert(ctx, mat->n_col > 0, goto error);
1369
226k
1370
226k
  if (mat->n_col > mat->n_row) {
1371
39.9k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1372
39.9k
    if (!bset)
1373
0
      goto error;
1374
186k
  } else if (mat->n_col < mat->n_row) {
1375
83.5k
    bset->dim = isl_space_cow(bset->dim);
1376
83.5k
    if (!bset->dim)
1377
0
      goto error;
1378
83.5k
    bset->dim->n_out -= mat->n_row - mat->n_col;
1379
83.5k
  }
1380
226k
1381
226k
  if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1382
226k
      isl_mat_copy(mat)) < 0)
1383
0
    goto error;
1384
226k
1385
226k
  if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1386
226k
      isl_mat_copy(mat)) < 0)
1387
0
    goto error;
1388
226k
1389
226k
  if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
1390
0
    goto error2;
1391
226k
1392
226k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1393
226k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1394
226k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1395
226k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1396
226k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1397
226k
1398
226k
  bset = isl_basic_set_simplify(bset);
1399
226k
  bset = isl_basic_set_finalize(bset);
1400
226k
1401
226k
  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
28.8k
{
1412
28.8k
  int i;
1413
28.8k
1414
28.8k
  set = isl_set_cow(set);
1415
28.8k
  if (!set)
1416
0
    goto error;
1417
28.8k
1418
84.5k
  
for (i = 0; 28.8k
i < set->n;
++i55.7k
) {
1419
55.7k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1420
55.7k
                isl_mat_copy(mat));
1421
55.7k
    if (!set->p[i])
1422
0
      goto error;
1423
55.7k
  }
1424
28.8k
  if (mat->n_col != mat->n_row) {
1425
13.5k
    set->dim = isl_space_cow(set->dim);
1426
13.5k
    if (!set->dim)
1427
0
      goto error;
1428
13.5k
    set->dim->n_out += mat->n_col;
1429
13.5k
    set->dim->n_out -= mat->n_row;
1430
13.5k
  }
1431
28.8k
  isl_mat_free(mat);
1432
28.8k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1433
28.8k
  return set;
1434
0
error:
1435
0
  isl_set_free(set);
1436
0
  isl_mat_free(mat);
1437
0
  return NULL;
1438
28.8k
}
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
34.8k
{
1501
34.8k
  int r;
1502
34.8k
1503
34.8k
  if (n == 0)
1504
22
    return mat;
1505
34.8k
1506
34.8k
  mat = isl_mat_cow(mat);
1507
34.8k
  if (check_col_range(mat, col, n) < 0)
1508
0
    return isl_mat_free(mat);
1509
34.8k
1510
34.8k
  if (col != mat->n_col-n) {
1511
62.3k
    for (r = 0; r < mat->n_row; 
++r48.6k
)
1512
48.6k
      isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1513
48.6k
          mat->n_col - col - n);
1514
13.7k
  }
1515
34.8k
  mat->n_col -= n;
1516
34.8k
  return mat;
1517
34.8k
}
1518
1519
__isl_give isl_mat *isl_mat_drop_rows(__isl_take isl_mat *mat,
1520
  unsigned row, unsigned n)
1521
45.2k
{
1522
45.2k
  int r;
1523
45.2k
1524
45.2k
  mat = isl_mat_cow(mat);
1525
45.2k
  if (check_row_range(mat, row, n) < 0)
1526
0
    return isl_mat_free(mat);
1527
45.2k
1528
125k
  
for (r = row; 45.2k
r+n < mat->n_row;
++r79.8k
)
1529
79.8k
    mat->row[r] = mat->row[r+n];
1530
45.2k
1531
45.2k
  mat->n_row -= n;
1532
45.2k
  return mat;
1533
45.2k
}
1534
1535
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1536
        unsigned col, unsigned n)
1537
29.6k
{
1538
29.6k
  isl_mat *ext;
1539
29.6k
1540
29.6k
  if (check_col_range(mat, col, 0) < 0)
1541
0
    return isl_mat_free(mat);
1542
29.6k
  if (n == 0)
1543
37
    return mat;
1544
29.5k
1545
29.5k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1546
29.5k
  if (!ext)
1547
0
    goto error;
1548
29.5k
1549
29.5k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1550
29.5k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1551
29.5k
        col + n, col, mat->n_col - col);
1552
29.5k
1553
29.5k
  isl_mat_free(mat);
1554
29.5k
  return ext;
1555
0
error:
1556
0
  isl_mat_free(mat);
1557
0
  return NULL;
1558
29.5k
}
1559
1560
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1561
  unsigned first, unsigned n)
1562
29.6k
{
1563
29.6k
  int i;
1564
29.6k
1565
29.6k
  if (!mat)
1566
0
    return NULL;
1567
29.6k
  mat = isl_mat_insert_cols(mat, first, n);
1568
29.6k
  if (!mat)
1569
0
    return NULL;
1570
29.6k
1571
30.6k
  
for (i = 0; 29.6k
i < mat->n_row;
++i1.01k
)
1572
1.01k
    isl_seq_clr(mat->row[i] + first, n);
1573
29.6k
1574
29.6k
  return mat;
1575
29.6k
}
1576
1577
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1578
1.66k
{
1579
1.66k
  if (!mat)
1580
0
    return NULL;
1581
1.66k
1582
1.66k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1583
1.66k
}
1584
1585
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1586
        unsigned row, unsigned n)
1587
2.00k
{
1588
2.00k
  isl_mat *ext;
1589
2.00k
1590
2.00k
  if (check_row_range(mat, row, 0) < 0)
1591
0
    return isl_mat_free(mat);
1592
2.00k
  if (n == 0)
1593
33
    return mat;
1594
1.97k
1595
1.97k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1596
1.97k
  if (!ext)
1597
0
    goto error;
1598
1.97k
1599
1.97k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1600
1.97k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1601
1.97k
        mat->n_row - row, 0, 0, mat->n_col);
1602
1.97k
1603
1.97k
  isl_mat_free(mat);
1604
1.97k
  return ext;
1605
0
error:
1606
0
  isl_mat_free(mat);
1607
0
  return NULL;
1608
1.97k
}
1609
1610
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1611
2.00k
{
1612
2.00k
  if (!mat)
1613
0
    return NULL;
1614
2.00k
1615
2.00k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1616
2.00k
}
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
1.51k
{
1644
1.51k
  int i;
1645
1.51k
1646
7.31k
  for (i = 0; i < mat->n_row; 
++i5.79k
)
1647
5.79k
    isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1648
1.51k
}
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
2.94k
{
1664
2.94k
  int i;
1665
2.94k
1666
10.3k
  for (i = 0; i < mat->n_row; 
++i7.37k
)
1667
7.37k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1668
2.94k
}
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
14
{
1718
14
  if (check_row(mat, row) < 0)
1719
0
    return isl_mat_free(mat);
1720
14
  if (isl_seq_first_non_zero(mat->row[row], mat->n_col) == -1)
1721
0
    return mat;
1722
14
  mat = isl_mat_cow(mat);
1723
14
  if (!mat)
1724
0
    return NULL;
1725
14
  isl_seq_neg(mat->row[row], mat->row[row], mat->n_col);
1726
14
  return mat;
1727
14
}
1728
1729
__isl_give isl_mat *isl_mat_unimodular_complete(__isl_take isl_mat *M, int row)
1730
4.82k
{
1731
4.82k
  int r;
1732
4.82k
  struct isl_mat *H = NULL, *Q = NULL;
1733
4.82k
1734
4.82k
  if (!M)
1735
0
    return NULL;
1736
4.82k
1737
4.82k
  isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1738
4.82k
  M->n_row = row;
1739
4.82k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1740
4.82k
  M->n_row = M->n_col;
1741
4.82k
  if (!H)
1742
0
    goto error;
1743
9.65k
  
for (r = 0; 4.82k
r < row;
++r4.82k
)
1744
4.82k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1745
14.3k
  
for (r = row; 4.82k
r < M->n_row;
++r9.53k
)
1746
9.53k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1747
4.82k
  isl_mat_free(H);
1748
4.82k
  isl_mat_free(Q);
1749
4.82k
  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
4.82k
}
1756
1757
__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1758
  __isl_take isl_mat *bot)
1759
129
{
1760
129
  struct isl_mat *mat;
1761
129
1762
129
  if (!top || !bot)
1763
0
    goto error;
1764
129
1765
129
  isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1766
129
  if (top->n_row == 0) {
1767
117
    isl_mat_free(top);
1768
117
    return bot;
1769
117
  }
1770
12
  if (bot->n_row == 0) {
1771
0
    isl_mat_free(bot);
1772
0
    return top;
1773
0
  }
1774
12
1775
12
  mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1776
12
  if (!mat)
1777
0
    goto error;
1778
12
  isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1779
12
       0, 0, mat->n_col);
1780
12
  isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1781
12
       0, 0, mat->n_col);
1782
12
  isl_mat_free(top);
1783
12
  isl_mat_free(bot);
1784
12
  return mat;
1785
0
error:
1786
0
  isl_mat_free(top);
1787
0
  isl_mat_free(bot);
1788
0
  return NULL;
1789
12
}
1790
1791
isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1792
2.11k
{
1793
2.11k
  int i;
1794
2.11k
1795
2.11k
  if (!mat1 || !mat2)
1796
0
    return isl_bool_error;
1797
2.11k
1798
2.11k
  if (mat1->n_row != mat2->n_row)
1799
118
    return isl_bool_false;
1800
2.00k
1801
2.00k
  if (mat1->n_col != mat2->n_col)
1802
0
    return isl_bool_false;
1803
2.00k
1804
3.24k
  
for (i = 0; 2.00k
i < mat1->n_row;
++i1.24k
)
1805
1.30k
    if (!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col))
1806
53
      return isl_bool_false;
1807
2.00k
1808
2.00k
  
return isl_bool_true1.94k
;
1809
2.00k
}
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
151
{
1859
151
  isl_mat *res;
1860
151
1861
151
  if (!mat)
1862
0
    return NULL;
1863
151
  if (n == 0 || dst_col == src_col)
1864
121
    return mat;
1865
30
1866
30
  res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1867
30
  if (!res)
1868
0
    goto error;
1869
30
1870
30
  if (dst_col < src_col) {
1871
30
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1872
30
         0, 0, dst_col);
1873
30
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1874
30
         dst_col, src_col, n);
1875
30
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1876
30
         dst_col + n, dst_col, src_col - dst_col);
1877
30
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1878
30
         src_col + n, src_col + n,
1879
30
         res->n_col - src_col - n);
1880
30
  } 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
30
  isl_mat_free(mat);
1892
30
1893
30
  return res;
1894
0
error:
1895
0
  isl_mat_free(mat);
1896
0
  return NULL;
1897
30
}
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
81
{
1914
81
  int i;
1915
81
  isl_int g;
1916
81
1917
81
  isl_int_set_si(*gcd, 0);
1918
81
  if (!mat)
1919
0
    return;
1920
81
1921
81
  isl_int_init(g);
1922
443
  for (i = 0; i < mat->n_row; 
++i362
) {
1923
362
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1924
362
    isl_int_gcd(*gcd, *gcd, g);
1925
362
  }
1926
81
  isl_int_clear(g);
1927
81
}
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
81
{
1950
81
  int i;
1951
81
1952
81
  if (isl_int_is_one(m))
1953
81
    
return mat0
;
1954
81
1955
81
  mat = isl_mat_cow(mat);
1956
81
  if (!mat)
1957
0
    return NULL;
1958
81
1959
443
  
for (i = 0; 81
i < mat->n_row;
++i362
)
1960
362
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1961
81
1962
81
  return mat;
1963
81
}
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
81
{
1982
81
  isl_int gcd;
1983
81
1984
81
  if (!mat)
1985
0
    return NULL;
1986
81
1987
81
  isl_int_init(gcd);
1988
81
  isl_mat_gcd(mat, &gcd);
1989
81
  mat = isl_mat_scale_down(mat, gcd);
1990
81
  isl_int_clear(gcd);
1991
81
1992
81
  return mat;
1993
81
}
1994
1995
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1996
457
{
1997
457
  mat = isl_mat_cow(mat);
1998
457
  if (!mat)
1999
0
    return NULL;
2000
457
2001
457
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
2002
457
2003
457
  return mat;
2004
457
}
2005
2006
/* Number of initial non-zero columns.
2007
 */
2008
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
2009
1.07k
{
2010
1.07k
  int i;
2011
1.07k
2012
1.07k
  if (!mat)
2013
0
    return -1;
2014
1.07k
2015
1.69k
  
for (i = 0; 1.07k
i < mat->n_col;
++i625
)
2016
1.28k
    if (row_first_non_zero(mat->row, mat->n_row, i) < 0)
2017
659
      break;
2018
1.07k
2019
1.07k
  return i;
2020
1.07k
}
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
}