Coverage Report

Created: 2017-10-03 07:32

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/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 2014      Ecole Normale Superieure
4
 *
5
 * Use of this software is governed by the MIT license
6
 *
7
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
8
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
9
 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
10
 */
11
12
#include <isl_ctx_private.h>
13
#include <isl_map_private.h>
14
#include <isl/space.h>
15
#include <isl_seq.h>
16
#include <isl_mat_private.h>
17
#include <isl_vec_private.h>
18
#include <isl_space_private.h>
19
#include <isl_val_private.h>
20
#include <isl/deprecated/mat_int.h>
21
22
isl_ctx *isl_mat_get_ctx(__isl_keep isl_mat *mat)
23
3.37M
{
24
3.37M
  return mat ? mat->ctx : NULL;
25
3.37M
}
26
27
/* Return a hash value that digests "mat".
28
 */
29
uint32_t isl_mat_get_hash(__isl_keep isl_mat *mat)
30
0
{
31
0
  int i;
32
0
  uint32_t hash;
33
0
34
0
  if (!mat)
35
0
    return 0;
36
0
37
0
  
hash = 0
isl_hash_init0
();
38
0
  isl_hash_byte(hash, mat->n_row & 0xFF);
39
0
  isl_hash_byte(hash, mat->n_col & 0xFF);
40
0
  for (i = 0; 
i < mat->n_row0
;
++i0
) {
41
0
    uint32_t row_hash;
42
0
43
0
    row_hash = isl_seq_get_hash(mat->row[i], mat->n_col);
44
0
    isl_hash_hash(hash, row_hash);
45
0
  }
46
0
47
0
  return hash;
48
0
}
49
50
struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
51
  unsigned n_row, unsigned n_col)
52
5.80M
{
53
5.80M
  int i;
54
5.80M
  struct isl_mat *mat;
55
5.80M
56
5.80M
  mat = isl_alloc_type(ctx, struct isl_mat);
57
5.80M
  if (!mat)
58
0
    return NULL;
59
5.80M
60
5.80M
  mat->row = NULL;
61
5.80M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
62
5.80M
  if (isl_blk_is_error(mat->block))
63
0
    goto error;
64
5.80M
  
mat->row = 5.80M
isl_alloc_array5.80M
(ctx, isl_int *, n_row);
65
5.80M
  if (
n_row && 5.80M
!mat->row4.29M
)
66
0
    goto error;
67
5.80M
68
31.9M
  
for (i = 0; 5.80M
i < n_row31.9M
;
++i26.1M
)
69
26.1M
    mat->row[i] = mat->block.data + i * n_col;
70
5.80M
71
5.80M
  mat->ctx = ctx;
72
5.80M
  isl_ctx_ref(ctx);
73
5.80M
  mat->ref = 1;
74
5.80M
  mat->n_row = n_row;
75
5.80M
  mat->n_col = n_col;
76
5.80M
  mat->max_col = n_col;
77
5.80M
  mat->flags = 0;
78
5.80M
79
5.80M
  return mat;
80
0
error:
81
0
  isl_blk_free(ctx, mat->block);
82
0
  free(mat);
83
0
  return NULL;
84
5.80M
}
85
86
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
87
  unsigned n_row, unsigned n_col)
88
109k
{
89
109k
  int i;
90
109k
  isl_int *old;
91
109k
  isl_int **row;
92
109k
93
109k
  if (!mat)
94
0
    return NULL;
95
109k
96
109k
  
if (109k
mat->max_col >= n_col && 109k
mat->n_row >= n_row105k
) {
97
7.80k
    if (mat->n_col < n_col)
98
22
      mat->n_col = n_col;
99
7.80k
    return mat;
100
7.80k
  }
101
101k
102
101k
  
if (101k
mat->max_col < n_col101k
) {
103
3.60k
    struct isl_mat *new_mat;
104
3.60k
105
3.60k
    if (n_row < mat->n_row)
106
8
      n_row = mat->n_row;
107
3.60k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
108
3.60k
    if (!new_mat)
109
0
      goto error;
110
25.0k
    
for (i = 0; 3.60k
i < mat->n_row25.0k
;
++i21.4k
)
111
21.4k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
112
3.60k
    isl_mat_free(mat);
113
3.60k
    return new_mat;
114
3.60k
  }
115
98.1k
116
98.1k
  mat = isl_mat_cow(mat);
117
98.1k
  if (!mat)
118
0
    goto error;
119
98.1k
120
98.1k
  old = mat->block.data;
121
98.1k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
122
98.1k
  if (isl_blk_is_error(mat->block))
123
0
    goto error;
124
98.1k
  
row = 98.1k
isl_realloc_array98.1k
(mat->ctx, mat->row, isl_int *, n_row);
125
98.1k
  if (
n_row && 98.1k
!row98.1k
)
126
0
    goto error;
127
98.1k
  mat->row = row;
128
98.1k
129
1.02M
  for (i = 0; 
i < mat->n_row1.02M
;
++i931k
)
130
931k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
131
343k
  for (i = mat->n_row; 
i < n_row343k
;
++i245k
)
132
245k
    mat->row[i] = mat->block.data + i * mat->max_col;
133
98.1k
  mat->n_row = n_row;
134
98.1k
  if (mat->n_col < n_col)
135
0
    mat->n_col = n_col;
136
98.1k
137
98.1k
  return mat;
138
0
error:
139
0
  isl_mat_free(mat);
140
0
  return NULL;
141
109k
}
142
143
__isl_give isl_mat *isl_mat_sub_alloc6(isl_ctx *ctx, isl_int **row,
144
  unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
145
2.77M
{
146
2.77M
  int i;
147
2.77M
  struct isl_mat *mat;
148
2.77M
149
2.77M
  mat = isl_alloc_type(ctx, struct isl_mat);
150
2.77M
  if (!mat)
151
0
    return NULL;
152
2.77M
  
mat->row = 2.77M
isl_alloc_array2.77M
(ctx, isl_int *, n_row);
153
2.77M
  if (
n_row && 2.77M
!mat->row1.53M
)
154
0
    goto error;
155
9.64M
  
for (i = 0; 2.77M
i < n_row9.64M
;
++i6.86M
)
156
6.86M
    mat->row[i] = row[first_row+i] + first_col;
157
2.77M
  mat->ctx = ctx;
158
2.77M
  isl_ctx_ref(ctx);
159
2.77M
  mat->ref = 1;
160
2.77M
  mat->n_row = n_row;
161
2.77M
  mat->n_col = n_col;
162
2.77M
  mat->block = isl_blk_empty();
163
2.77M
  mat->flags = ISL_MAT_BORROWED;
164
2.77M
  return mat;
165
0
error:
166
0
  free(mat);
167
0
  return NULL;
168
2.77M
}
169
170
__isl_give isl_mat *isl_mat_sub_alloc(__isl_keep isl_mat *mat,
171
  unsigned first_row, unsigned n_row, unsigned first_col, unsigned n_col)
172
609k
{
173
609k
  if (!mat)
174
0
    return NULL;
175
609k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
176
609k
          first_col, n_col);
177
609k
}
178
179
void isl_mat_sub_copy(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
180
  unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
181
406k
{
182
406k
  int i;
183
406k
184
1.42M
  for (i = 0; 
i < n_row1.42M
;
++i1.01M
)
185
1.01M
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
186
406k
}
187
188
void isl_mat_sub_neg(struct isl_ctx *ctx, isl_int **dst, isl_int **src,
189
  unsigned n_row, unsigned dst_col, unsigned src_col, unsigned n_col)
190
140k
{
191
140k
  int i;
192
140k
193
425k
  for (i = 0; 
i < n_row425k
;
++i285k
)
194
285k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
195
140k
}
196
197
__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
198
2.11M
{
199
2.11M
  if (!mat)
200
0
    return NULL;
201
2.11M
202
2.11M
  mat->ref++;
203
2.11M
  return mat;
204
2.11M
}
205
206
__isl_give isl_mat *isl_mat_dup(__isl_keep isl_mat *mat)
207
357k
{
208
357k
  int i;
209
357k
  struct isl_mat *mat2;
210
357k
211
357k
  if (!mat)
212
2.85k
    return NULL;
213
354k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
214
354k
  if (!mat2)
215
0
    return NULL;
216
1.14M
  
for (i = 0; 354k
i < mat->n_row1.14M
;
++i794k
)
217
794k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
218
357k
  return mat2;
219
357k
}
220
221
__isl_give isl_mat *isl_mat_cow(__isl_take isl_mat *mat)
222
3.94M
{
223
3.94M
  struct isl_mat *mat2;
224
3.94M
  if (!mat)
225
0
    return NULL;
226
3.94M
227
3.94M
  
if (3.94M
mat->ref == 1 && 3.94M
!3.91M
ISL_F_ISSET3.91M
(mat, ISL_MAT_BORROWED))
228
3.59M
    return mat;
229
351k
230
351k
  mat2 = isl_mat_dup(mat);
231
351k
  isl_mat_free(mat);
232
351k
  return mat2;
233
351k
}
234
235
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
236
12.0M
{
237
12.0M
  if (!mat)
238
1.34M
    return NULL;
239
10.6M
240
10.6M
  
if (10.6M
--mat->ref > 010.6M
)
241
2.11M
    return NULL;
242
8.58M
243
8.58M
  
if (8.58M
!8.58M
ISL_F_ISSET8.58M
(mat, ISL_MAT_BORROWED))
244
5.80M
    isl_blk_free(mat->ctx, mat->block);
245
12.0M
  isl_ctx_deref(mat->ctx);
246
12.0M
  free(mat->row);
247
12.0M
  free(mat);
248
12.0M
249
12.0M
  return NULL;
250
12.0M
}
251
252
int isl_mat_rows(__isl_keep isl_mat *mat)
253
118k
{
254
118k
  return mat ? 
mat->n_row118k
:
-10
;
255
118k
}
256
257
int isl_mat_cols(__isl_keep isl_mat *mat)
258
42.3k
{
259
42.3k
  return mat ? 
mat->n_col42.3k
:
-10
;
260
42.3k
}
261
262
/* Check that "col" is a valid column position for "mat".
263
 */
264
static isl_stat check_col(__isl_keep isl_mat *mat, int col)
265
4.62k
{
266
4.62k
  if (!mat)
267
0
    return isl_stat_error;
268
4.62k
  
if (4.62k
col < 0 || 4.62k
col >= mat->n_col4.62k
)
269
0
    isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
270
4.62k
      "column out of range", return isl_stat_error);
271
4.62k
  return isl_stat_ok;
272
4.62k
}
273
274
/* Check that "row" is a valid row position for "mat".
275
 */
276
static isl_stat check_row(__isl_keep isl_mat *mat, int row)
277
4.63k
{
278
4.63k
  if (!mat)
279
0
    return isl_stat_error;
280
4.63k
  
if (4.63k
row < 0 || 4.63k
row >= mat->n_row4.63k
)
281
0
    isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
282
4.63k
      "row out of range", return isl_stat_error);
283
4.63k
  return isl_stat_ok;
284
4.63k
}
285
286
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
287
3.02k
{
288
3.02k
  if (check_row(mat, row) < 0)
289
0
    return -1;
290
3.02k
  
if (3.02k
check_col(mat, col) < 03.02k
)
291
0
    return -1;
292
3.02k
  
isl_int_set3.02k
(*v, mat->row[row][col]);
293
3.02k
  return 0;
294
3.02k
}
295
296
/* Extract the element at row "row", oolumn "col" of "mat".
297
 */
298
__isl_give isl_val *isl_mat_get_element_val(__isl_keep isl_mat *mat,
299
  int row, int col)
300
24
{
301
24
  isl_ctx *ctx;
302
24
303
24
  if (check_row(mat, row) < 0)
304
0
    return NULL;
305
24
  
if (24
check_col(mat, col) < 024
)
306
0
    return NULL;
307
24
  ctx = isl_mat_get_ctx(mat);
308
24
  return isl_val_int_from_isl_int(ctx, mat->row[row][col]);
309
24
}
310
311
__isl_give isl_mat *isl_mat_set_element(__isl_take isl_mat *mat,
312
  int row, int col, isl_int v)
313
1.51k
{
314
1.51k
  mat = isl_mat_cow(mat);
315
1.51k
  if (check_row(mat, row) < 0)
316
0
    return isl_mat_free(mat);
317
1.51k
  
if (1.51k
check_col(mat, col) < 01.51k
)
318
0
    return isl_mat_free(mat);
319
1.51k
  
isl_int_set1.51k
(mat->row[row][col], v);
320
1.51k
  return mat;
321
1.51k
}
322
323
__isl_give isl_mat *isl_mat_set_element_si(__isl_take isl_mat *mat,
324
  int row, int col, int v)
325
43
{
326
43
  mat = isl_mat_cow(mat);
327
43
  if (check_row(mat, row) < 0)
328
0
    return isl_mat_free(mat);
329
43
  
if (43
check_col(mat, col) < 043
)
330
0
    return isl_mat_free(mat);
331
43
  
isl_int_set_si43
(mat->row[row][col], v);
332
43
  return mat;
333
43
}
334
335
/* Replace the element at row "row", column "col" of "mat" by "v".
336
 */
337
__isl_give isl_mat *isl_mat_set_element_val(__isl_take isl_mat *mat,
338
  int row, int col, __isl_take isl_val *v)
339
0
{
340
0
  if (!v)
341
0
    return isl_mat_free(mat);
342
0
  
if (0
!isl_val_is_int(v)0
)
343
0
    isl_die(isl_val_get_ctx(v), isl_error_invalid,
344
0
      "expecting integer value", goto error);
345
0
  mat = isl_mat_set_element(mat, row, col, v->n);
346
0
  isl_val_free(v);
347
0
  return mat;
348
0
error:
349
0
  isl_val_free(v);
350
0
  return isl_mat_free(mat);
351
0
}
352
353
__isl_give isl_mat *isl_mat_diag(isl_ctx *ctx, unsigned n_row, isl_int d)
354
994k
{
355
994k
  int i;
356
994k
  struct isl_mat *mat;
357
994k
358
994k
  mat = isl_mat_alloc(ctx, n_row, n_row);
359
994k
  if (!mat)
360
0
    return NULL;
361
5.34M
  
for (i = 0; 994k
i < n_row5.34M
;
++i4.35M
) {
362
4.35M
    isl_seq_clr(mat->row[i], i);
363
4.35M
    isl_int_set(mat->row[i][i], d);
364
4.35M
    isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
365
4.35M
  }
366
994k
367
994k
  return mat;
368
994k
}
369
370
/* Create an "n_row" by "n_col" matrix with zero elements.
371
 */
372
__isl_give isl_mat *isl_mat_zero(isl_ctx *ctx, unsigned n_row, unsigned n_col)
373
151
{
374
151
  int i;
375
151
  isl_mat *mat;
376
151
377
151
  mat = isl_mat_alloc(ctx, n_row, n_col);
378
151
  if (!mat)
379
0
    return NULL;
380
192
  
for (i = 0; 151
i < n_row192
;
++i41
)
381
41
    isl_seq_clr(mat->row[i], n_col);
382
151
383
151
  return mat;
384
151
}
385
386
__isl_give isl_mat *isl_mat_identity(isl_ctx *ctx, unsigned n_row)
387
994k
{
388
994k
  if (!ctx)
389
0
    return NULL;
390
994k
  return isl_mat_diag(ctx, n_row, ctx->one);
391
994k
}
392
393
/* Is "mat" a (possibly scaled) identity matrix?
394
 */
395
int isl_mat_is_scaled_identity(__isl_keep isl_mat *mat)
396
0
{
397
0
  int i;
398
0
399
0
  if (!mat)
400
0
    return -1;
401
0
  
if (0
mat->n_row != mat->n_col0
)
402
0
    return 0;
403
0
404
0
  
for (i = 0; 0
i < mat->n_row0
;
++i0
) {
405
0
    if (isl_seq_first_non_zero(mat->row[i], i) != -1)
406
0
      return 0;
407
0
    
if (0
isl_int_ne0
(mat->row[0][0], mat->row[i][i]))
408
0
      return 0;
409
0
    
if (0
isl_seq_first_non_zero(mat->row[i] + i + 1,
410
0
              mat->n_col - (i + 1)) != -1)
411
0
      return 0;
412
0
  }
413
0
414
0
  return 1;
415
0
}
416
417
__isl_give isl_vec *isl_mat_vec_product(__isl_take isl_mat *mat,
418
  __isl_take isl_vec *vec)
419
337k
{
420
337k
  int i;
421
337k
  struct isl_vec *prod;
422
337k
423
337k
  if (
!mat || 337k
!vec337k
)
424
0
    goto error;
425
337k
426
337k
  
isl_assert337k
(mat->ctx, mat->n_col == vec->size, goto error);
427
337k
428
337k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
429
337k
  if (!prod)
430
0
    goto error;
431
337k
432
2.61M
  
for (i = 0; 337k
i < prod->size2.61M
;
++i2.27M
)
433
2.27M
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
434
2.27M
          &prod->block.data[i]);
435
337k
  isl_mat_free(mat);
436
337k
  isl_vec_free(vec);
437
337k
  return prod;
438
0
error:
439
0
  isl_mat_free(mat);
440
0
  isl_vec_free(vec);
441
0
  return NULL;
442
337k
}
443
444
__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
445
  __isl_take isl_vec *vec)
446
160
{
447
160
  struct isl_mat *vec_mat;
448
160
  int i;
449
160
450
160
  if (
!mat || 160
!vec160
)
451
0
    goto error;
452
160
  vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
453
160
  if (!vec_mat)
454
0
    goto error;
455
1.18k
  
for (i = 0; 160
i < vec->size1.18k
;
++i1.02k
)
456
1.02k
    isl_int_set(vec_mat->row[i][0], vec->el[i]);
457
160
  vec_mat = isl_mat_inverse_product(mat, vec_mat);
458
160
  isl_vec_free(vec);
459
160
  if (!vec_mat)
460
0
    return NULL;
461
160
  vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
462
160
  if (vec)
463
1.18k
    
for (i = 0; 160
i < vec->size1.18k
;
++i1.02k
)
464
1.02k
      isl_int_set(vec->el[i], vec_mat->row[i][0]);
465
160
  isl_mat_free(vec_mat);
466
160
  return vec;
467
0
error:
468
0
  isl_mat_free(mat);
469
0
  isl_vec_free(vec);
470
0
  return NULL;
471
160
}
472
473
__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,
474
  __isl_take isl_mat *mat)
475
19.5k
{
476
19.5k
  int i, j;
477
19.5k
  struct isl_vec *prod;
478
19.5k
479
19.5k
  if (
!mat || 19.5k
!vec19.5k
)
480
0
    goto error;
481
19.5k
482
19.5k
  
isl_assert19.5k
(mat->ctx, mat->n_row == vec->size, goto error);
483
19.5k
484
19.5k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
485
19.5k
  if (!prod)
486
0
    goto error;
487
19.5k
488
110k
  
for (i = 0; 19.5k
i < prod->size110k
;
++i90.5k
) {
489
90.5k
    isl_int_set_si(prod->el[i], 0);
490
698k
    for (j = 0; 
j < vec->size698k
;
++j608k
)
491
608k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
492
90.5k
  }
493
19.5k
  isl_mat_free(mat);
494
19.5k
  isl_vec_free(vec);
495
19.5k
  return prod;
496
0
error:
497
0
  isl_mat_free(mat);
498
0
  isl_vec_free(vec);
499
0
  return NULL;
500
19.5k
}
501
502
__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
503
  __isl_take isl_mat *right)
504
140k
{
505
140k
  int i;
506
140k
  struct isl_mat *sum;
507
140k
508
140k
  if (
!left || 140k
!right140k
)
509
0
    goto error;
510
140k
511
140k
  
isl_assert140k
(left->ctx, left->n_row == right->n_row, goto error);
512
140k
  
isl_assert140k
(left->ctx, left->n_row >= 1, goto error);
513
140k
  
isl_assert140k
(left->ctx, left->n_col >= 1, goto error);
514
140k
  
isl_assert140k
(left->ctx, right->n_col >= 1, goto error);
515
140k
  
isl_assert140k
(left->ctx,
516
140k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
517
140k
      goto error);
518
140k
  
isl_assert140k
(left->ctx,
519
140k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
520
140k
      goto error);
521
140k
522
140k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
523
140k
  if (!sum)
524
0
    goto error;
525
140k
  
isl_int_lcm140k
(sum->row[0][0], left->row[0][0], right->row[0][0]);
526
140k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
527
140k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
528
140k
529
140k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
530
1.06M
  for (i = 1; 
i < sum->n_row1.06M
;
++i922k
) {
531
922k
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
532
922k
    isl_int_addmul(sum->row[i][0],
533
922k
        right->row[0][0], right->row[i][0]);
534
922k
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
535
922k
        left->n_col-1);
536
922k
    isl_seq_scale(sum->row[i]+left->n_col,
537
922k
        right->row[i]+1, right->row[0][0],
538
922k
        right->n_col-1);
539
922k
  }
540
140k
541
140k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
542
140k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
543
140k
  isl_mat_free(left);
544
140k
  isl_mat_free(right);
545
140k
  return sum;
546
0
error:
547
0
  isl_mat_free(left);
548
0
  isl_mat_free(right);
549
0
  return NULL;
550
140k
}
551
552
static void exchange(struct isl_mat *M, struct isl_mat **U,
553
  struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
554
872k
{
555
872k
  int r;
556
6.52M
  for (r = row; 
r < M->n_row6.52M
;
++r5.65M
)
557
5.65M
    isl_int_swap(M->row[r][i], M->row[r][j]);
558
872k
  if (
U872k
) {
559
8.28M
    for (r = 0; 
r < (*U)->n_row8.28M
;
++r7.44M
)
560
7.44M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
561
844k
  }
562
872k
  if (Q)
563
324k
    isl_mat_swap_rows(*Q, i, j);
564
872k
}
565
566
static void subtract(struct isl_mat *M, struct isl_mat **U,
567
  struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
568
722k
{
569
722k
  int r;
570
3.01M
  for (r = row; 
r < M->n_row3.01M
;
++r2.29M
)
571
2.29M
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
572
722k
  if (
U722k
) {
573
4.49M
    for (r = 0; 
r < (*U)->n_row4.49M
;
++r3.80M
)
574
3.80M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
575
684k
  }
576
722k
  if (
Q722k
) {
577
1.31M
    for (r = 0; 
r < (*Q)->n_col1.31M
;
++r1.10M
)
578
1.10M
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
579
214k
  }
580
722k
}
581
582
static void oppose(struct isl_mat *M, struct isl_mat **U,
583
  struct isl_mat **Q, unsigned row, unsigned col)
584
250k
{
585
250k
  int r;
586
1.07M
  for (r = row; 
r < M->n_row1.07M
;
++r823k
)
587
823k
    isl_int_neg(M->row[r][col], M->row[r][col]);
588
250k
  if (
U250k
) {
589
1.43M
    for (r = 0; 
r < (*U)->n_row1.43M
;
++r1.21M
)
590
1.21M
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
591
218k
  }
592
250k
  if (Q)
593
119k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
594
250k
}
595
596
/* Given matrix M, compute
597
 *
598
 *    M U = H
599
 *    M   = H Q
600
 *
601
 * with U and Q unimodular matrices and H a matrix in column echelon form
602
 * such that on each echelon row the entries in the non-echelon column
603
 * are non-negative (if neg == 0) or non-positive (if neg == 1)
604
 * and strictly smaller (in absolute value) than the entries in the echelon
605
 * column.
606
 * If U or Q are NULL, then these matrices are not computed.
607
 */
608
__isl_give isl_mat *isl_mat_left_hermite(__isl_take isl_mat *M, int neg,
609
  __isl_give isl_mat **U, __isl_give isl_mat **Q)
610
470k
{
611
470k
  isl_int c;
612
470k
  int row, col;
613
470k
614
470k
  if (U)
615
435k
    *U = NULL;
616
470k
  if (Q)
617
230k
    *Q = NULL;
618
470k
  if (!M)
619
0
    goto error;
620
470k
  M = isl_mat_cow(M);
621
470k
  if (!M)
622
0
    goto error;
623
470k
  
if (470k
U470k
) {
624
435k
    *U = isl_mat_identity(M->ctx, M->n_col);
625
435k
    if (!*U)
626
0
      goto error;
627
470k
  }
628
470k
  
if (470k
Q470k
) {
629
230k
    *Q = isl_mat_identity(M->ctx, M->n_col);
630
230k
    if (!*Q)
631
0
      goto error;
632
470k
  }
633
470k
634
470k
  col = 0;
635
470k
  isl_int_init(c);
636
2.20M
  for (row = 0; 
row < M->n_row2.20M
;
++row1.73M
) {
637
1.73M
    int first, i, off;
638
1.73M
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
639
1.73M
    if (first == -1)
640
556k
      continue;
641
1.18M
    first += col;
642
1.18M
    if (first != col)
643
617k
      exchange(M, U, Q, row, first, col);
644
1.18M
    if (isl_int_is_neg(M->row[row][col]))
645
208k
      oppose(M, U, Q, row, col);
646
1.18M
    first = col+1;
647
1.54M
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
648
1.18M
                   M->n_col-first)) != -1) {
649
360k
      first += off;
650
360k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
651
360k
      subtract(M, U, Q, row, col, first, c);
652
360k
      if (
!360k
isl_int_is_zero360k
(M->row[row][first]))
653
30.1k
        exchange(M, U, Q, row, first, col);
654
360k
      else
655
329k
        ++first;
656
360k
    }
657
4.65M
    for (i = 0; 
i < col4.65M
;
++i3.47M
) {
658
3.47M
      if (isl_int_is_zero(M->row[row][i]))
659
3.32M
        continue;
660
152k
      
if (152k
neg152k
)
661
0
        isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
662
152k
      else
663
152k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
664
152k
      if (isl_int_is_zero(c))
665
52.6k
        continue;
666
99.7k
      subtract(M, U, Q, row, col, i, c);
667
99.7k
    }
668
1.73M
    ++col;
669
1.73M
  }
670
470k
  isl_int_clear(c);
671
470k
672
470k
  return M;
673
0
error:
674
0
  if (
Q0
) {
675
0
    isl_mat_free(*Q);
676
0
    *Q = NULL;
677
0
  }
678
0
  if (
U0
) {
679
0
    isl_mat_free(*U);
680
0
    *U = NULL;
681
0
  }
682
0
  isl_mat_free(M);
683
0
  return NULL;
684
470k
}
685
686
/* Use row "row" of "mat" to eliminate column "col" from all other rows.
687
 */
688
static __isl_give isl_mat *eliminate(__isl_take isl_mat *mat, int row, int col)
689
1.19k
{
690
1.19k
  int k, nr, nc;
691
1.19k
  isl_ctx *ctx;
692
1.19k
693
1.19k
  if (!mat)
694
0
    return NULL;
695
1.19k
696
1.19k
  ctx = isl_mat_get_ctx(mat);
697
1.19k
  nr = isl_mat_rows(mat);
698
1.19k
  nc = isl_mat_cols(mat);
699
1.19k
700
3.30k
  for (k = 0; 
k < nr3.30k
;
++k2.11k
) {
701
2.11k
    if (k == row)
702
1.19k
      continue;
703
922
    
if (922
isl_int_is_zero922
(mat->row[k][col]))
704
922
      continue;
705
0
    mat = isl_mat_cow(mat);
706
0
    if (!mat)
707
0
      return NULL;
708
0
    isl_seq_elim(mat->row[k], mat->row[row], col, nc, NULL);
709
0
    isl_seq_normalize(ctx, mat->row[k], nc);
710
0
  }
711
1.19k
712
1.19k
  return mat;
713
1.19k
}
714
715
/* Perform Gaussian elimination on the rows of "mat", but start
716
 * from the final row and the final column.
717
 * Any zero rows that result from the elimination are removed.
718
 *
719
 * In particular, for each column from last to first,
720
 * look for the last row with a non-zero coefficient in that column,
721
 * move it last (but before other rows moved last in previous steps) and
722
 * use it to eliminate the column from the other rows.
723
 */
724
__isl_give isl_mat *isl_mat_reverse_gauss(__isl_take isl_mat *mat)
725
1.00k
{
726
1.00k
  int k, row, last, nr, nc;
727
1.00k
728
1.00k
  if (!mat)
729
0
    return NULL;
730
1.00k
731
1.00k
  nr = isl_mat_rows(mat);
732
1.00k
  nc = isl_mat_cols(mat);
733
1.00k
734
1.00k
  last = nc - 1;
735
2.19k
  for (row = nr - 1; 
row >= 02.19k
;
--row1.19k
) {
736
1.58k
    for (; 
last >= 01.58k
;
--last389
) {
737
2.05k
      for (k = row; 
k >= 02.05k
;
--k469
)
738
1.66k
        
if (1.66k
!1.66k
isl_int_is_zero1.66k
(mat->row[k][last]))
739
1.19k
          break;
740
1.58k
      if (k >= 0)
741
1.19k
        break;
742
1.58k
    }
743
1.19k
    if (last < 0)
744
0
      break;
745
1.19k
    
if (1.19k
k != row1.19k
)
746
0
      mat = isl_mat_swap_rows(mat, k, row);
747
1.19k
    if (!mat)
748
0
      return NULL;
749
1.19k
    
if (1.19k
isl_int_is_neg1.19k
(mat->row[row][last]))
750
4
      mat = isl_mat_row_neg(mat, row);
751
1.19k
    mat = eliminate(mat, row, last);
752
1.19k
    if (!mat)
753
0
      return NULL;
754
1.19k
  }
755
1.00k
  mat = isl_mat_drop_rows(mat, 0, row + 1);
756
1.00k
757
1.00k
  return mat;
758
1.00k
}
759
760
/* Negate the lexicographically negative rows of "mat" such that
761
 * all rows in the result are lexicographically non-negative.
762
 */
763
__isl_give isl_mat *isl_mat_lexnonneg_rows(__isl_take isl_mat *mat)
764
1.00k
{
765
1.00k
  int i, nr, nc;
766
1.00k
767
1.00k
  if (!mat)
768
0
    return NULL;
769
1.00k
770
1.00k
  nr = isl_mat_rows(mat);
771
1.00k
  nc = isl_mat_cols(mat);
772
1.00k
773
2.19k
  for (i = 0; 
i < nr2.19k
;
++i1.19k
) {
774
1.19k
    int pos;
775
1.19k
776
1.19k
    pos = isl_seq_first_non_zero(mat->row[i], nc);
777
1.19k
    if (pos < 0)
778
0
      continue;
779
1.19k
    
if (1.19k
isl_int_is_nonneg1.19k
(mat->row[i][pos]))
780
1.17k
      continue;
781
14
    mat = isl_mat_row_neg(mat, i);
782
14
    if (!mat)
783
0
      return NULL;
784
1.19k
  }
785
1.00k
786
1.00k
  return mat;
787
1.00k
}
788
789
struct isl_mat *isl_mat_right_kernel(struct isl_mat *mat)
790
0
{
791
0
  int i, rank;
792
0
  struct isl_mat *U = NULL;
793
0
  struct isl_mat *K;
794
0
795
0
  mat = isl_mat_left_hermite(mat, 0, &U, NULL);
796
0
  if (
!mat || 0
!U0
)
797
0
    goto error;
798
0
799
0
  
for (i = 0, rank = 0; 0
rank < mat->n_col0
;
++rank0
) {
800
0
    while (
i < mat->n_row && 0
isl_int_is_zero0
(mat->row[i][rank]))
801
0
      ++i;
802
0
    if (i >= mat->n_row)
803
0
      break;
804
0
  }
805
0
  K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);
806
0
  if (!K)
807
0
    goto error;
808
0
  isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
809
0
  isl_mat_free(mat);
810
0
  isl_mat_free(U);
811
0
  return K;
812
0
error:
813
0
  isl_mat_free(mat);
814
0
  isl_mat_free(U);
815
0
  return NULL;
816
0
}
817
818
__isl_give isl_mat *isl_mat_lin_to_aff(__isl_take isl_mat *mat)
819
676k
{
820
676k
  int i;
821
676k
  struct isl_mat *mat2;
822
676k
823
676k
  if (!mat)
824
0
    return NULL;
825
676k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
826
676k
  if (!mat2)
827
0
    goto error;
828
676k
  
isl_int_set_si676k
(mat2->row[0][0], 1);
829
676k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
830
4.13M
  for (i = 0; 
i < mat->n_row4.13M
;
++i3.46M
) {
831
3.46M
    isl_int_set_si(mat2->row[1+i][0], 0);
832
3.46M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
833
3.46M
  }
834
676k
  isl_mat_free(mat);
835
676k
  return mat2;
836
0
error:
837
0
  isl_mat_free(mat);
838
0
  return NULL;
839
676k
}
840
841
/* Given two matrices M1 and M2, return the block matrix
842
 *
843
 *  [ M1  0  ]
844
 *  [ 0   M2 ]
845
 */
846
__isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1,
847
  __isl_take isl_mat *mat2)
848
127k
{
849
127k
  int i;
850
127k
  isl_mat *mat;
851
127k
852
127k
  if (
!mat1 || 127k
!mat2127k
)
853
0
    goto error;
854
127k
855
127k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
856
127k
               mat1->n_col + mat2->n_col);
857
127k
  if (!mat)
858
0
    goto error;
859
255k
  
for (i = 0; 127k
i < mat1->n_row255k
;
++i127k
) {
860
127k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
861
127k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
862
127k
  }
863
814k
  for (i = 0; 
i < mat2->n_row814k
;
++i686k
) {
864
686k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
865
686k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
866
686k
                mat2->row[i], mat2->n_col);
867
686k
  }
868
127k
  isl_mat_free(mat1);
869
127k
  isl_mat_free(mat2);
870
127k
  return mat;
871
0
error:
872
0
  isl_mat_free(mat1);
873
0
  isl_mat_free(mat2);
874
0
  return NULL;
875
127k
}
876
877
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
878
925k
{
879
925k
  int i;
880
925k
881
1.56M
  for (i = 0; 
i < n_row1.56M
;
++i635k
)
882
1.09M
    
if (1.09M
!1.09M
isl_int_is_zero1.09M
(row[i][col]))
883
464k
      return i;
884
461k
  return -1;
885
925k
}
886
887
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
888
460k
{
889
460k
  int i, min = row_first_non_zero(row, n_row, col);
890
460k
  if (min < 0)
891
0
    return -1;
892
1.09M
  
for (i = min + 1; 460k
i < n_row1.09M
;
++i636k
) {
893
636k
    if (isl_int_is_zero(row[i][col]))
894
632k
      continue;
895
3.55k
    
if (3.55k
isl_int_abs_lt3.55k
(row[i][col], row[min][col]))
896
50
      min = i;
897
636k
  }
898
460k
  return min;
899
460k
}
900
901
static isl_stat inv_exchange(__isl_keep isl_mat **left,
902
  __isl_keep isl_mat **right, unsigned i, unsigned j)
903
341
{
904
341
  *left = isl_mat_swap_rows(*left, i, j);
905
341
  *right = isl_mat_swap_rows(*right, i, j);
906
341
907
341
  if (
!*left || 341
!*right341
)
908
0
    return isl_stat_error;
909
341
  return isl_stat_ok;
910
341
}
911
912
static void inv_oppose(
913
  struct isl_mat *left, struct isl_mat *right, unsigned row)
914
106
{
915
106
  isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
916
106
  isl_seq_neg(right->row[row], right->row[row], right->n_col);
917
106
}
918
919
static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
920
  unsigned row, unsigned i, isl_int m)
921
3.55k
{
922
3.55k
  isl_int_neg(m, m);
923
3.55k
  isl_seq_combine(left->row[i]+row,
924
3.55k
      left->ctx->one, left->row[i]+row,
925
3.55k
      m, left->row[row]+row,
926
3.55k
      left->n_col-row);
927
3.55k
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
928
3.55k
      m, right->row[row], right->n_col);
929
3.55k
}
930
931
/* Compute inv(left)*right
932
 */
933
__isl_give isl_mat *isl_mat_inverse_product(__isl_take isl_mat *left,
934
  __isl_take isl_mat *right)
935
155k
{
936
155k
  int row;
937
155k
  isl_int a, b;
938
155k
939
155k
  if (
!left || 155k
!right155k
)
940
0
    goto error;
941
155k
942
155k
  
isl_assert155k
(left->ctx, left->n_row == left->n_col, goto error);
943
155k
  
isl_assert155k
(left->ctx, left->n_row == right->n_row, goto error);
944
155k
945
155k
  
if (155k
left->n_row == 0155k
) {
946
0
    isl_mat_free(left);
947
0
    return right;
948
0
  }
949
155k
950
155k
  left = isl_mat_cow(left);
951
155k
  right = isl_mat_cow(right);
952
155k
  if (
!left || 155k
!right155k
)
953
0
    goto error;
954
155k
955
155k
  
isl_int_init155k
(a);
956
155k
  isl_int_init(b);
957
616k
  for (row = 0; 
row < left->n_row616k
;
++row460k
) {
958
460k
    int pivot, first, i, off;
959
460k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
960
460k
    if (
pivot < 0460k
) {
961
0
      isl_int_clear(a);
962
0
      isl_int_clear(b);
963
0
      isl_assert(left->ctx, pivot >= 0, goto error);
964
0
    }
965
460k
    pivot += row;
966
460k
    if (pivot != row)
967
339
      
if (339
inv_exchange(&left, &right, pivot, row) < 0339
)
968
0
        goto error;
969
460k
    
if (460k
isl_int_is_neg460k
(left->row[row][row]))
970
106
      inv_oppose(left, right, row);
971
460k
    first = row+1;
972
463k
    while ((off = row_first_non_zero(left->row+first,
973
460k
          left->n_row-first, row)) != -1) {
974
3.55k
      first += off;
975
3.55k
      isl_int_fdiv_q(a, left->row[first][row],
976
3.55k
          left->row[row][row]);
977
3.55k
      inv_subtract(left, right, row, first, a);
978
3.55k
      if (
!3.55k
isl_int_is_zero3.55k
(left->row[first][row])) {
979
2
        if (inv_exchange(&left, &right, row, first) < 0)
980
0
          goto error;
981
3.55k
      } else {
982
3.55k
        ++first;
983
3.55k
      }
984
3.55k
    }
985
1.09M
    
for (i = 0; 460k
i < row1.09M
;
++i637k
) {
986
637k
      if (isl_int_is_zero(left->row[i][row]))
987
636k
        continue;
988
347
      
isl_int_gcd347
(a, left->row[row][row], left->row[i][row]);
989
347
      isl_int_divexact(b, left->row[i][row], a);
990
347
      isl_int_divexact(a, left->row[row][row], a);
991
347
      isl_int_neg(b, b);
992
637k
      isl_seq_combine(left->row[i] + i,
993
637k
          a, left->row[i] + i,
994
637k
          b, left->row[row] + i,
995
637k
          left->n_col - i);
996
637k
      isl_seq_combine(right->row[i], a, right->row[i],
997
637k
          b, right->row[row], right->n_col);
998
637k
    }
999
460k
  }
1000
155k
  
isl_int_clear155k
(b);
1001
155k
1002
155k
  isl_int_set(a, left->row[0][0]);
1003
460k
  for (row = 1; 
row < left->n_row460k
;
++row304k
)
1004
304k
    isl_int_lcm(a, a, left->row[row][row]);
1005
155k
  if (
isl_int_is_zero155k
(a)){
1006
0
    isl_int_clear(a);
1007
0
    isl_assert(left->ctx, 0, goto error);
1008
0
  }
1009
616k
  
for (row = 0; 155k
row < left->n_row616k
;
++row460k
) {
1010
460k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
1011
460k
    if (isl_int_is_one(left->row[row][row]))
1012
453k
      continue;
1013
6.93k
    isl_seq_scale(right->row[row], right->row[row],
1014
6.93k
        left->row[row][row], right->n_col);
1015
6.93k
  }
1016
155k
  isl_int_clear(a);
1017
155k
1018
155k
  isl_mat_free(left);
1019
155k
  return right;
1020
0
error:
1021
0
  isl_mat_free(left);
1022
0
  isl_mat_free(right);
1023
0
  return NULL;
1024
155k
}
1025
1026
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
1027
58.5k
{
1028
58.5k
  int i;
1029
58.5k
1030
287k
  for (i = 0; 
i < mat->n_row287k
;
++i228k
)
1031
228k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
1032
58.5k
}
1033
1034
void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
1035
  isl_int m1, unsigned src1, isl_int m2, unsigned src2)
1036
170k
{
1037
170k
  int i;
1038
170k
  isl_int tmp;
1039
170k
1040
170k
  isl_int_init(tmp);
1041
791k
  for (i = 0; 
i < mat->n_row791k
;
++i620k
) {
1042
620k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
1043
620k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
1044
620k
    isl_int_set(mat->row[i][dst], tmp);
1045
620k
  }
1046
170k
  isl_int_clear(tmp);
1047
170k
}
1048
1049
__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
1050
60.6k
{
1051
60.6k
  struct isl_mat *inv;
1052
60.6k
  int row;
1053
60.6k
  isl_int a, b;
1054
60.6k
1055
60.6k
  mat = isl_mat_cow(mat);
1056
60.6k
  if (!mat)
1057
0
    return NULL;
1058
60.6k
1059
60.6k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
1060
60.6k
  inv = isl_mat_cow(inv);
1061
60.6k
  if (!inv)
1062
0
    goto error;
1063
60.6k
1064
60.6k
  
isl_int_init60.6k
(a);
1065
60.6k
  isl_int_init(b);
1066
246k
  for (row = 0; 
row < mat->n_row246k
;
++row185k
) {
1067
185k
    int pivot, first, i, off;
1068
185k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
1069
185k
    if (
pivot < 0185k
) {
1070
0
      isl_int_clear(a);
1071
0
      isl_int_clear(b);
1072
0
      isl_assert(mat->ctx, pivot >= 0, goto error);
1073
0
    }
1074
185k
    pivot += row;
1075
185k
    if (pivot != row)
1076
29.1k
      exchange(mat, &inv, NULL, row, pivot, row);
1077
185k
    if (isl_int_is_neg(mat->row[row][row]))
1078
42.1k
      oppose(mat, &inv, NULL, row, row);
1079
185k
    first = row+1;
1080
447k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
1081
185k
                mat->n_col-first)) != -1) {
1082
262k
      first += off;
1083
262k
      isl_int_fdiv_q(a, mat->row[row][first],
1084
262k
                mat->row[row][row]);
1085
262k
      subtract(mat, &inv, NULL, row, row, first, a);
1086
262k
      if (
!262k
isl_int_is_zero262k
(mat->row[row][first]))
1087
196k
        exchange(mat, &inv, NULL, row, row, first);
1088
262k
      else
1089
66.0k
        ++first;
1090
262k
    }
1091
458k
    for (i = 0; 
i < row458k
;
++i273k
) {
1092
273k
      if (isl_int_is_zero(mat->row[row][i]))
1093
187k
        continue;
1094
85.3k
      
isl_int_gcd85.3k
(a, mat->row[row][row], mat->row[row][i]);
1095
85.3k
      isl_int_divexact(b, mat->row[row][i], a);
1096
85.3k
      isl_int_divexact(a, mat->row[row][row], a);
1097
85.3k
      isl_int_neg(a, a);
1098
273k
      isl_mat_col_combine(mat, i, a, i, b, row);
1099
273k
      isl_mat_col_combine(inv, i, a, i, b, row);
1100
273k
    }
1101
185k
  }
1102
60.6k
  
isl_int_clear60.6k
(b);
1103
60.6k
1104
60.6k
  isl_int_set(a, mat->row[0][0]);
1105
185k
  for (row = 1; 
row < mat->n_row185k
;
++row125k
)
1106
125k
    isl_int_lcm(a, a, mat->row[row][row]);
1107
60.6k
  if (
isl_int_is_zero60.6k
(a)){
1108
0
    isl_int_clear(a);
1109
0
    goto error;
1110
0
  }
1111
246k
  
for (row = 0; 60.6k
row < mat->n_row246k
;
++row185k
) {
1112
185k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1113
185k
    if (isl_int_is_one(mat->row[row][row]))
1114
127k
      continue;
1115
58.5k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1116
58.5k
  }
1117
60.6k
  isl_int_clear(a);
1118
60.6k
1119
60.6k
  isl_mat_free(mat);
1120
60.6k
1121
60.6k
  return inv;
1122
0
error:
1123
0
  isl_mat_free(mat);
1124
0
  isl_mat_free(inv);
1125
0
  return NULL;
1126
60.6k
}
1127
1128
__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1129
3.46k
{
1130
3.46k
  struct isl_mat *transpose = NULL;
1131
3.46k
  int i, j;
1132
3.46k
1133
3.46k
  if (!mat)
1134
0
    return NULL;
1135
3.46k
1136
3.46k
  
if (3.46k
mat->n_col == mat->n_row3.46k
) {
1137
3.46k
    mat = isl_mat_cow(mat);
1138
3.46k
    if (!mat)
1139
0
      return NULL;
1140
18.6k
    
for (i = 0; 3.46k
i < mat->n_row18.6k
;
++i15.2k
)
1141
54.7k
      
for (j = i + 1; 15.2k
j < mat->n_col54.7k
;
++j39.4k
)
1142
39.4k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1143
3.46k
    return mat;
1144
3.46k
  }
1145
0
  transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
1146
0
  if (!transpose)
1147
0
    goto error;
1148
0
  
for (i = 0; 0
i < mat->n_row0
;
++i0
)
1149
0
    
for (j = 0; 0
j < mat->n_col0
;
++j0
)
1150
0
      isl_int_set(transpose->row[j][i], mat->row[i][j]);
1151
0
  isl_mat_free(mat);
1152
0
  return transpose;
1153
0
error:
1154
0
  isl_mat_free(mat);
1155
0
  return NULL;
1156
3.46k
}
1157
1158
__isl_give isl_mat *isl_mat_swap_cols(__isl_take isl_mat *mat,
1159
  unsigned i, unsigned j)
1160
647k
{
1161
647k
  int r;
1162
647k
1163
647k
  mat = isl_mat_cow(mat);
1164
647k
  if (!mat)
1165
0
    return NULL;
1166
647k
  
isl_assert647k
(mat->ctx, i < mat->n_col, goto error);
1167
647k
  
isl_assert647k
(mat->ctx, j < mat->n_col, goto error);
1168
647k
1169
18.4M
  
for (r = 0; 647k
r < mat->n_row18.4M
;
++r17.8M
)
1170
17.8M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1171
647k
  return mat;
1172
0
error:
1173
0
  isl_mat_free(mat);
1174
0
  return NULL;
1175
647k
}
1176
1177
__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1178
  unsigned i, unsigned j)
1179
2.08M
{
1180
2.08M
  isl_int *t;
1181
2.08M
1182
2.08M
  if (!mat)
1183
0
    return NULL;
1184
2.08M
  mat = isl_mat_cow(mat);
1185
2.08M
  if (!mat)
1186
0
    return NULL;
1187
2.08M
  t = mat->row[i];
1188
2.08M
  mat->row[i] = mat->row[j];
1189
2.08M
  mat->row[j] = t;
1190
2.08M
  return mat;
1191
2.08M
}
1192
1193
/* Calculate the product of two matrices.
1194
 *
1195
 * This function is optimized for operand matrices that contain many zeros and
1196
 * skips multiplications where we know one of the operands is zero.
1197
 */
1198
__isl_give isl_mat *isl_mat_product(__isl_take isl_mat *left,
1199
  __isl_take isl_mat *right)
1200
1.99M
{
1201
1.99M
  int i, j, k;
1202
1.99M
  struct isl_mat *prod;
1203
1.99M
1204
1.99M
  if (
!left || 1.99M
!right1.99M
)
1205
0
    goto error;
1206
1.99M
  
isl_assert1.99M
(left->ctx, left->n_col == right->n_row, goto error);
1207
1.99M
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1208
1.99M
  if (!prod)
1209
0
    goto error;
1210
1.99M
  
if (1.99M
left->n_col == 01.99M
) {
1211
35
    for (i = 0; 
i < prod->n_row35
;
++i8
)
1212
8
      isl_seq_clr(prod->row[i], prod->n_col);
1213
27
    isl_mat_free(left);
1214
27
    isl_mat_free(right);
1215
27
    return prod;
1216
27
  }
1217
6.76M
  
for (i = 0; 1.99M
i < prod->n_row6.76M
;
++i4.76M
) {
1218
39.2M
    for (j = 0; 
j < prod->n_col39.2M
;
++j34.4M
)
1219
34.4M
      isl_int_mul(prod->row[i][j],
1220
4.76M
            left->row[i][0], right->row[0][j]);
1221
44.3M
    for (k = 1; 
k < left->n_col44.3M
;
++k39.6M
) {
1222
39.6M
      if (isl_int_is_zero(left->row[i][k]))
1223
34.3M
        continue;
1224
48.4M
      
for (j = 0; 5.27M
j < prod->n_col48.4M
;
++j43.2M
)
1225
43.2M
        isl_int_addmul(prod->row[i][j],
1226
39.6M
              left->row[i][k], right->row[k][j]);
1227
39.6M
    }
1228
4.76M
  }
1229
1.99M
  isl_mat_free(left);
1230
1.99M
  isl_mat_free(right);
1231
1.99M
  return prod;
1232
0
error:
1233
0
  isl_mat_free(left);
1234
0
  isl_mat_free(right);
1235
0
  return NULL;
1236
1.99M
}
1237
1238
/* Replace the variables x in the rows q by x' given by x = M x',
1239
 * with M the matrix mat.
1240
 *
1241
 * If the number of new variables is greater than the original
1242
 * number of variables, then the rows q have already been
1243
 * preextended.  If the new number is smaller, then the coefficients
1244
 * of the divs, which are not changed, need to be shifted down.
1245
 * The row q may be the equalities, the inequalities or the
1246
 * div expressions.  In the latter case, has_div is true and
1247
 * we need to take into account the extra denominator column.
1248
 */
1249
static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
1250
  unsigned n_div, int has_div, struct isl_mat *mat)
1251
1.69M
{
1252
1.69M
  int i;
1253
1.69M
  struct isl_mat *t;
1254
1.69M
  int e;
1255
1.69M
1256
1.69M
  if (mat->n_col >= mat->n_row)
1257
1.05M
    e = 0;
1258
1.69M
  else
1259
632k
    e = mat->n_row - mat->n_col;
1260
1.69M
  if (has_div)
1261
564k
    
for (i = 0; 564k
i < n564k
;
++i130
)
1262
130
      isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
1263
1.69M
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1264
1.69M
  t = isl_mat_product(t, mat);
1265
1.69M
  if (!t)
1266
0
    return -1;
1267
4.75M
  
for (i = 0; 1.69M
i < n4.75M
;
++i3.06M
) {
1268
3.06M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1269
3.06M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1270
3.06M
          q[i] + has_div + t->n_col + e, n_div);
1271
3.06M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1272
3.06M
  }
1273
1.69M
  isl_mat_free(t);
1274
1.69M
  return 0;
1275
1.69M
}
1276
1277
/* Replace the variables x in bset by x' given by x = M x', with
1278
 * M the matrix mat.
1279
 *
1280
 * If there are fewer variables x' then there are x, then we perform
1281
 * the transformation in place, which means that, in principle,
1282
 * this frees up some extra variables as the number
1283
 * of columns remains constant, but we would have to extend
1284
 * the div array too as the number of rows in this array is assumed
1285
 * to be equal to extra.
1286
 */
1287
__isl_give isl_basic_set *isl_basic_set_preimage(
1288
  __isl_take isl_basic_set *bset, __isl_take isl_mat *mat)
1289
564k
{
1290
564k
  struct isl_ctx *ctx;
1291
564k
1292
564k
  if (
!bset || 564k
!mat564k
)
1293
0
    goto error;
1294
564k
1295
564k
  ctx = bset->ctx;
1296
564k
  bset = isl_basic_set_cow(bset);
1297
564k
  if (!bset)
1298
0
    goto error;
1299
564k
1300
564k
  
isl_assert564k
(ctx, bset->dim->nparam == 0, goto error);
1301
564k
  
isl_assert564k
(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1302
564k
  
isl_assert564k
(ctx, mat->n_col > 0, goto error);
1303
564k
1304
564k
  
if (564k
mat->n_col > mat->n_row564k
) {
1305
102k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1306
102k
    if (!bset)
1307
0
      goto error;
1308
461k
  } else 
if (461k
mat->n_col < mat->n_row461k
) {
1309
210k
    bset->dim = isl_space_cow(bset->dim);
1310
210k
    if (!bset->dim)
1311
0
      goto error;
1312
210k
    bset->dim->n_out -= mat->n_row - mat->n_col;
1313
210k
  }
1314
564k
1315
564k
  
if (564k
preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1316
564k
      isl_mat_copy(mat)) < 0)
1317
0
    goto error;
1318
564k
1319
564k
  
if (564k
preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1320
564k
      isl_mat_copy(mat)) < 0)
1321
0
    goto error;
1322
564k
1323
564k
  
if (564k
preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0564k
)
1324
0
    goto error2;
1325
564k
1326
564k
  
ISL_F_CLR564k
(bset, ISL_BASIC_SET_NO_IMPLICIT);
1327
564k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1328
564k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1329
564k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1330
564k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1331
564k
1332
564k
  bset = isl_basic_set_simplify(bset);
1333
564k
  bset = isl_basic_set_finalize(bset);
1334
564k
1335
564k
  return bset;
1336
0
error:
1337
0
  isl_mat_free(mat);
1338
0
error2:
1339
0
  isl_basic_set_free(bset);
1340
0
  return NULL;
1341
564k
}
1342
1343
__isl_give isl_set *isl_set_preimage(
1344
  __isl_take isl_set *set, __isl_take isl_mat *mat)
1345
32.2k
{
1346
32.2k
  int i;
1347
32.2k
1348
32.2k
  set = isl_set_cow(set);
1349
32.2k
  if (!set)
1350
0
    goto error;
1351
32.2k
1352
91.4k
  
for (i = 0; 32.2k
i < set->n91.4k
;
++i59.2k
) {
1353
59.2k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1354
59.2k
                isl_mat_copy(mat));
1355
59.2k
    if (!set->p[i])
1356
0
      goto error;
1357
59.2k
  }
1358
32.2k
  
if (32.2k
mat->n_col != mat->n_row32.2k
) {
1359
13.6k
    set->dim = isl_space_cow(set->dim);
1360
13.6k
    if (!set->dim)
1361
0
      goto error;
1362
13.6k
    set->dim->n_out += mat->n_col;
1363
13.6k
    set->dim->n_out -= mat->n_row;
1364
13.6k
  }
1365
32.2k
  isl_mat_free(mat);
1366
32.2k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1367
32.2k
  return set;
1368
0
error:
1369
0
  isl_set_free(set);
1370
0
  isl_mat_free(mat);
1371
0
  return NULL;
1372
32.2k
}
1373
1374
/* Replace the variables x starting at "first_col" in the rows "rows"
1375
 * of some coefficient matrix by x' with x = M x' with M the matrix mat.
1376
 * That is, replace the corresponding coefficients c by c M.
1377
 */
1378
isl_stat isl_mat_sub_transform(isl_int **row, unsigned n_row,
1379
  unsigned first_col, __isl_take isl_mat *mat)
1380
144
{
1381
144
  int i;
1382
144
  isl_ctx *ctx;
1383
144
  isl_mat *t;
1384
144
1385
144
  if (!mat)
1386
0
    return isl_stat_error;
1387
144
  ctx = isl_mat_get_ctx(mat);
1388
144
  t = isl_mat_sub_alloc6(ctx, row, 0, n_row, first_col, mat->n_row);
1389
144
  t = isl_mat_product(t, mat);
1390
144
  if (!t)
1391
0
    return isl_stat_error;
1392
220
  
for (i = 0; 144
i < n_row220
;
++i76
)
1393
76
    isl_seq_swp_or_cpy(row[i] + first_col, t->row[i], t->n_col);
1394
144
  isl_mat_free(t);
1395
144
  return isl_stat_ok;
1396
144
}
1397
1398
void isl_mat_print_internal(__isl_keep isl_mat *mat, FILE *out, int indent)
1399
0
{
1400
0
  int i, j;
1401
0
1402
0
  if (
!mat0
) {
1403
0
    fprintf(out, "%*snull mat\n", indent, "");
1404
0
    return;
1405
0
  }
1406
0
1407
0
  
if (0
mat->n_row == 00
)
1408
0
    fprintf(out, "%*s[]\n", indent, "");
1409
0
1410
0
  for (i = 0; 
i < mat->n_row0
;
++i0
) {
1411
0
    if (!i)
1412
0
      fprintf(out, "%*s[[", indent, "");
1413
0
    else
1414
0
      fprintf(out, "%*s[", indent+1, "");
1415
0
    for (j = 0; 
j < mat->n_col0
;
++j0
) {
1416
0
      if (j)
1417
0
          fprintf(out, ",");
1418
0
      isl_int_print(out, mat->row[i][j], 0);
1419
0
    }
1420
0
    if (i == mat->n_row-1)
1421
0
      fprintf(out, "]]\n");
1422
0
    else
1423
0
      fprintf(out, "]\n");
1424
0
  }
1425
0
}
1426
1427
void isl_mat_dump(__isl_keep isl_mat *mat)
1428
0
{
1429
0
  isl_mat_print_internal(mat, stderr, 0);
1430
0
}
1431
1432
__isl_give isl_mat *isl_mat_drop_cols(__isl_take isl_mat *mat,
1433
  unsigned col, unsigned n)
1434
85.8k
{
1435
85.8k
  int r;
1436
85.8k
1437
85.8k
  if (n == 0)
1438
2
    return mat;
1439
85.8k
1440
85.8k
  mat = isl_mat_cow(mat);
1441
85.8k
  if (!mat)
1442
0
    return NULL;
1443
85.8k
1444
85.8k
  
if (85.8k
col != mat->n_col-n85.8k
) {
1445
66.1k
    for (r = 0; 
r < mat->n_row66.1k
;
++r51.7k
)
1446
51.7k
      isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1447
51.7k
          mat->n_col - col - n);
1448
14.3k
  }
1449
85.8k
  mat->n_col -= n;
1450
85.8k
  return mat;
1451
85.8k
}
1452
1453
__isl_give isl_mat *isl_mat_drop_rows(__isl_take isl_mat *mat,
1454
  unsigned row, unsigned n)
1455
115k
{
1456
115k
  int r;
1457
115k
1458
115k
  mat = isl_mat_cow(mat);
1459
115k
  if (!mat)
1460
0
    return NULL;
1461
115k
1462
339k
  
for (r = row; 115k
r+n < mat->n_row339k
;
++r223k
)
1463
223k
    mat->row[r] = mat->row[r+n];
1464
115k
1465
115k
  mat->n_row -= n;
1466
115k
  return mat;
1467
115k
}
1468
1469
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1470
        unsigned col, unsigned n)
1471
79.0k
{
1472
79.0k
  isl_mat *ext;
1473
79.0k
1474
79.0k
  if (!mat)
1475
0
    return NULL;
1476
79.0k
  
if (79.0k
n == 079.0k
)
1477
435
    return mat;
1478
78.6k
1479
78.6k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1480
78.6k
  if (!ext)
1481
0
    goto error;
1482
78.6k
1483
78.6k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1484
78.6k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1485
78.6k
        col + n, col, mat->n_col - col);
1486
78.6k
1487
78.6k
  isl_mat_free(mat);
1488
78.6k
  return ext;
1489
0
error:
1490
0
  isl_mat_free(mat);
1491
0
  return NULL;
1492
79.0k
}
1493
1494
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1495
  unsigned first, unsigned n)
1496
79.0k
{
1497
79.0k
  int i;
1498
79.0k
1499
79.0k
  if (!mat)
1500
0
    return NULL;
1501
79.0k
  mat = isl_mat_insert_cols(mat, first, n);
1502
79.0k
  if (!mat)
1503
0
    return NULL;
1504
79.0k
1505
82.9k
  
for (i = 0; 79.0k
i < mat->n_row82.9k
;
++i3.85k
)
1506
3.85k
    isl_seq_clr(mat->row[i] + first, n);
1507
79.0k
1508
79.0k
  return mat;
1509
79.0k
}
1510
1511
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1512
4.50k
{
1513
4.50k
  if (!mat)
1514
0
    return NULL;
1515
4.50k
1516
4.50k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1517
4.50k
}
1518
1519
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1520
        unsigned row, unsigned n)
1521
4.95k
{
1522
4.95k
  isl_mat *ext;
1523
4.95k
1524
4.95k
  if (!mat)
1525
0
    return NULL;
1526
4.95k
  
if (4.95k
n == 04.95k
)
1527
433
    return mat;
1528
4.52k
1529
4.52k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1530
4.52k
  if (!ext)
1531
0
    goto error;
1532
4.52k
1533
4.52k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1534
4.52k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1535
4.52k
        mat->n_row - row, 0, 0, mat->n_col);
1536
4.52k
1537
4.52k
  isl_mat_free(mat);
1538
4.52k
  return ext;
1539
0
error:
1540
0
  isl_mat_free(mat);
1541
0
  return NULL;
1542
4.95k
}
1543
1544
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1545
4.95k
{
1546
4.95k
  if (!mat)
1547
0
    return NULL;
1548
4.95k
1549
4.95k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1550
4.95k
}
1551
1552
__isl_give isl_mat *isl_mat_insert_zero_rows(__isl_take isl_mat *mat,
1553
  unsigned row, unsigned n)
1554
2
{
1555
2
  int i;
1556
2
1557
2
  mat = isl_mat_insert_rows(mat, row, n);
1558
2
  if (!mat)
1559
0
    return NULL;
1560
2
  
1561
4
  
for (i = 0; 2
i < n4
;
++i2
)
1562
2
    isl_seq_clr(mat->row[row + i], mat->n_col);
1563
2
1564
2
  return mat;
1565
2
}
1566
1567
__isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n)
1568
0
{
1569
0
  if (!mat)
1570
0
    return NULL;
1571
0
1572
0
  return isl_mat_insert_zero_rows(mat, mat->n_row, n);
1573
0
}
1574
1575
void isl_mat_col_submul(struct isl_mat *mat,
1576
      int dst_col, isl_int f, int src_col)
1577
5.07k
{
1578
5.07k
  int i;
1579
5.07k
1580
32.5k
  for (i = 0; 
i < mat->n_row32.5k
;
++i27.4k
)
1581
27.4k
    isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1582
5.07k
}
1583
1584
void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1585
2
{
1586
2
  int i;
1587
2
1588
2
  if (!mat)
1589
0
    return;
1590
2
1591
7
  
for (i = 0; 2
i < mat->n_row7
;
++i5
)
1592
5
    isl_int_add(mat->row[i][dst_col],
1593
2
          mat->row[i][dst_col], mat->row[i][src_col]);
1594
2
}
1595
1596
void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1597
10.8k
{
1598
10.8k
  int i;
1599
10.8k
1600
48.4k
  for (i = 0; 
i < mat->n_row48.4k
;
++i37.6k
)
1601
37.6k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1602
10.8k
}
1603
1604
/* Add "f" times column "src_col" to column "dst_col" of "mat" and
1605
 * return the result.
1606
 */
1607
__isl_give isl_mat *isl_mat_col_addmul(__isl_take isl_mat *mat, int dst_col,
1608
  isl_int f, int src_col)
1609
6
{
1610
6
  int i;
1611
6
1612
6
  if (
check_col(mat, dst_col) < 0 || 6
check_col(mat, src_col) < 06
)
1613
0
    return isl_mat_free(mat);
1614
6
1615
14
  
for (i = 0; 6
i < mat->n_row14
;
++i8
) {
1616
8
    if (isl_int_is_zero(mat->row[i][src_col]))
1617
2
      continue;
1618
6
    mat = isl_mat_cow(mat);
1619
6
    if (!mat)
1620
0
      return NULL;
1621
6
    
isl_int_addmul6
(mat->row[i][dst_col], f, mat->row[i][src_col]);
1622
6
  }
1623
6
1624
6
  return mat;
1625
6
}
1626
1627
/* Negate column "col" of "mat" and return the result.
1628
 */
1629
__isl_give isl_mat *isl_mat_col_neg(__isl_take isl_mat *mat, int col)
1630
4
{
1631
4
  int i;
1632
4
1633
4
  if (check_col(mat, col) < 0)
1634
0
    return isl_mat_free(mat);
1635
4
1636
10
  
for (i = 0; 4
i < mat->n_row10
;
++i6
) {
1637
6
    if (isl_int_is_zero(mat->row[i][col]))
1638
2
      continue;
1639
4
    mat = isl_mat_cow(mat);
1640
4
    if (!mat)
1641
0
      return NULL;
1642
4
    
isl_int_neg4
(mat->row[i][col], mat->row[i][col]);
1643
4
  }
1644
4
1645
4
  return mat;
1646
4
}
1647
1648
/* Negate row "row" of "mat" and return the result.
1649
 */
1650
__isl_give isl_mat *isl_mat_row_neg(__isl_take isl_mat *mat, int row)
1651
18
{
1652
18
  if (check_row(mat, row) < 0)
1653
0
    return isl_mat_free(mat);
1654
18
  
if (18
isl_seq_first_non_zero(mat->row[row], mat->n_col) == -118
)
1655
0
    return mat;
1656
18
  mat = isl_mat_cow(mat);
1657
18
  if (!mat)
1658
0
    return NULL;
1659
18
  isl_seq_neg(mat->row[row], mat->row[row], mat->n_col);
1660
18
  return mat;
1661
18
}
1662
1663
__isl_give isl_mat *isl_mat_unimodular_complete(__isl_take isl_mat *M, int row)
1664
17.6k
{
1665
17.6k
  int r;
1666
17.6k
  struct isl_mat *H = NULL, *Q = NULL;
1667
17.6k
1668
17.6k
  if (!M)
1669
0
    return NULL;
1670
17.6k
1671
17.6k
  
isl_assert17.6k
(M->ctx, M->n_row == M->n_col, goto error);
1672
17.6k
  M->n_row = row;
1673
17.6k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1674
17.6k
  M->n_row = M->n_col;
1675
17.6k
  if (!H)
1676
0
    goto error;
1677
35.3k
  
for (r = 0; 17.6k
r < row35.3k
;
++r17.6k
)
1678
17.6k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1679
65.4k
  
for (r = row; 17.6k
r < M->n_row65.4k
;
++r47.7k
)
1680
47.7k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1681
17.6k
  isl_mat_free(H);
1682
17.6k
  isl_mat_free(Q);
1683
17.6k
  return M;
1684
0
error:
1685
0
  isl_mat_free(H);
1686
0
  isl_mat_free(Q);
1687
0
  isl_mat_free(M);
1688
0
  return NULL;
1689
17.6k
}
1690
1691
__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1692
  __isl_take isl_mat *bot)
1693
420
{
1694
420
  struct isl_mat *mat;
1695
420
1696
420
  if (
!top || 420
!bot420
)
1697
0
    goto error;
1698
420
1699
420
  
isl_assert420
(top->ctx, top->n_col == bot->n_col, goto error);
1700
420
  
if (420
top->n_row == 0420
) {
1701
380
    isl_mat_free(top);
1702
380
    return bot;
1703
380
  }
1704
40
  
if (40
bot->n_row == 040
) {
1705
0
    isl_mat_free(bot);
1706
0
    return top;
1707
0
  }
1708
40
1709
40
  mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1710
40
  if (!mat)
1711
0
    goto error;
1712
40
  isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1713
40
       0, 0, mat->n_col);
1714
40
  isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1715
40
       0, 0, mat->n_col);
1716
40
  isl_mat_free(top);
1717
40
  isl_mat_free(bot);
1718
40
  return mat;
1719
0
error:
1720
0
  isl_mat_free(top);
1721
0
  isl_mat_free(bot);
1722
0
  return NULL;
1723
420
}
1724
1725
isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1726
42.9k
{
1727
42.9k
  int i;
1728
42.9k
1729
42.9k
  if (
!mat1 || 42.9k
!mat242.9k
)
1730
0
    return isl_bool_error;
1731
42.9k
1732
42.9k
  
if (42.9k
mat1->n_row != mat2->n_row42.9k
)
1733
6.90k
    return isl_bool_false;
1734
36.0k
1735
36.0k
  
if (36.0k
mat1->n_col != mat2->n_col36.0k
)
1736
0
    return isl_bool_false;
1737
36.0k
1738
70.0k
  
for (i = 0; 36.0k
i < mat1->n_row70.0k
;
++i34.0k
)
1739
34.1k
    
if (34.1k
!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col)34.1k
)
1740
190
      return isl_bool_false;
1741
36.0k
1742
35.8k
  return isl_bool_true;
1743
42.9k
}
1744
1745
__isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1746
0
{
1747
0
  struct isl_mat *mat;
1748
0
1749
0
  if (!vec)
1750
0
    return NULL;
1751
0
  mat = isl_mat_alloc(vec->ctx, 1, vec->size);
1752
0
  if (!mat)
1753
0
    goto error;
1754
0
1755
0
  isl_seq_cpy(mat->row[0], vec->el, vec->size);
1756
0
1757
0
  isl_vec_free(vec);
1758
0
  return mat;
1759
0
error:
1760
0
  isl_vec_free(vec);
1761
0
  return NULL;
1762
0
}
1763
1764
/* Return a copy of row "row" of "mat" as an isl_vec.
1765
 */
1766
__isl_give isl_vec *isl_mat_get_row(__isl_keep isl_mat *mat, unsigned row)
1767
24
{
1768
24
  isl_vec *v;
1769
24
1770
24
  if (!mat)
1771
0
    return NULL;
1772
24
  
if (24
row >= mat->n_row24
)
1773
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
1774
24
      return NULL);
1775
24
1776
24
  v = isl_vec_alloc(isl_mat_get_ctx(mat), mat->n_col);
1777
24
  if (!v)
1778
0
    return NULL;
1779
24
  isl_seq_cpy(v->el, mat->row[row], mat->n_col);
1780
24
1781
24
  return v;
1782
24
}
1783
1784
__isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1785
  __isl_take isl_vec *bot)
1786
0
{
1787
0
  return isl_mat_concat(top, isl_mat_from_row_vec(bot));
1788
0
}
1789
1790
__isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1791
  unsigned dst_col, unsigned src_col, unsigned n)
1792
567
{
1793
567
  isl_mat *res;
1794
567
1795
567
  if (!mat)
1796
0
    return NULL;
1797
567
  
if (567
n == 0 || 567
dst_col == src_col567
)
1798
447
    return mat;
1799
120
1800
120
  res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1801
120
  if (!res)
1802
0
    goto error;
1803
120
1804
120
  
if (120
dst_col < src_col120
) {
1805
120
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1806
120
         0, 0, dst_col);
1807
120
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1808
120
         dst_col, src_col, n);
1809
120
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1810
120
         dst_col + n, dst_col, src_col - dst_col);
1811
120
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1812
120
         src_col + n, src_col + n,
1813
120
         res->n_col - src_col - n);
1814
120
  } else {
1815
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1816
0
         0, 0, src_col);
1817
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1818
0
         src_col, src_col + n, dst_col - src_col);
1819
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1820
0
         dst_col, src_col, n);
1821
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1822
0
         dst_col + n, dst_col + n,
1823
0
         res->n_col - dst_col - n);
1824
0
  }
1825
120
  isl_mat_free(mat);
1826
120
1827
120
  return res;
1828
0
error:
1829
0
  isl_mat_free(mat);
1830
0
  return NULL;
1831
567
}
1832
1833
/* Return the gcd of the elements in row "row" of "mat" in *gcd.
1834
 * Return isl_stat_ok on success and isl_stat_error on failure.
1835
 */
1836
isl_stat isl_mat_row_gcd(__isl_keep isl_mat *mat, int row, isl_int *gcd)
1837
7
{
1838
7
  if (check_row(mat, row) < 0)
1839
0
    return isl_stat_error;
1840
7
1841
7
  isl_seq_gcd(mat->row[row], mat->n_col, gcd);
1842
7
1843
7
  return isl_stat_ok;
1844
7
}
1845
1846
void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1847
771
{
1848
771
  int i;
1849
771
  isl_int g;
1850
771
1851
771
  isl_int_set_si(*gcd, 0);
1852
771
  if (!mat)
1853
0
    return;
1854
771
1855
771
  
isl_int_init771
(g);
1856
4.72k
  for (i = 0; 
i < mat->n_row4.72k
;
++i3.95k
) {
1857
3.95k
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1858
3.95k
    isl_int_gcd(*gcd, *gcd, g);
1859
3.95k
  }
1860
771
  isl_int_clear(g);
1861
771
}
1862
1863
/* Return the result of scaling "mat" by a factor of "m".
1864
 */
1865
__isl_give isl_mat *isl_mat_scale(__isl_take isl_mat *mat, isl_int m)
1866
0
{
1867
0
  int i;
1868
0
1869
0
  if (isl_int_is_one(m))
1870
0
    return mat;
1871
0
1872
0
  mat = isl_mat_cow(mat);
1873
0
  if (!mat)
1874
0
    return NULL;
1875
0
1876
0
  
for (i = 0; 0
i < mat->n_row0
;
++i0
)
1877
0
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
1878
0
1879
0
  return mat;
1880
0
}
1881
1882
__isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1883
771
{
1884
771
  int i;
1885
771
1886
771
  if (isl_int_is_one(m))
1887
0
    return mat;
1888
771
1889
771
  mat = isl_mat_cow(mat);
1890
771
  if (!mat)
1891
0
    return NULL;
1892
771
1893
4.72k
  
for (i = 0; 771
i < mat->n_row4.72k
;
++i3.95k
)
1894
3.95k
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1895
771
1896
771
  return mat;
1897
771
}
1898
1899
__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row,
1900
  isl_int m)
1901
10
{
1902
10
  if (isl_int_is_one(m))
1903
0
    return mat;
1904
10
1905
10
  mat = isl_mat_cow(mat);
1906
10
  if (!mat)
1907
0
    return NULL;
1908
10
1909
10
  isl_seq_scale_down(mat->row[row], mat->row[row], m, mat->n_col);
1910
10
1911
10
  return mat;
1912
10
}
1913
1914
__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1915
771
{
1916
771
  isl_int gcd;
1917
771
1918
771
  if (!mat)
1919
0
    return NULL;
1920
771
1921
771
  
isl_int_init771
(gcd);
1922
771
  isl_mat_gcd(mat, &gcd);
1923
771
  mat = isl_mat_scale_down(mat, gcd);
1924
771
  isl_int_clear(gcd);
1925
771
1926
771
  return mat;
1927
771
}
1928
1929
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1930
1.63k
{
1931
1.63k
  mat = isl_mat_cow(mat);
1932
1.63k
  if (!mat)
1933
0
    return NULL;
1934
1.63k
1935
1.63k
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
1936
1.63k
1937
1.63k
  return mat;
1938
1.63k
}
1939
1940
/* Number of initial non-zero columns.
1941
 */
1942
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
1943
1.00k
{
1944
1.00k
  int i;
1945
1.00k
1946
1.00k
  if (!mat)
1947
0
    return -1;
1948
1.00k
1949
1.46k
  
for (i = 0; 1.00k
i < mat->n_col1.46k
;
++i456
)
1950
1.26k
    
if (1.26k
row_first_non_zero(mat->row, mat->n_row, i) < 01.26k
)
1951
812
      break;
1952
1.00k
1953
1.00k
  return i;
1954
1.00k
}