Coverage Report

Created: 2017-06-23 12:40

/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
2.47M
{
24
2.47M
  return mat ? mat->ctx : NULL;
25
2.47M
}
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
)
{0
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
4.06M
{
53
4.06M
  int i;
54
4.06M
  struct isl_mat *mat;
55
4.06M
56
4.06M
  mat = isl_alloc_type(ctx, struct isl_mat);
57
4.06M
  if (!mat)
58
0
    return NULL;
59
4.06M
60
4.06M
  mat->row = NULL;
61
4.06M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
62
4.06M
  if (isl_blk_is_error(mat->block))
63
0
    goto error;
64
4.06M
  
mat->row = 4.06M
isl_alloc_array4.06M
(ctx, isl_int *, n_row);
65
4.06M
  if (
n_row && 4.06M
!mat->row2.97M
)
66
0
    goto error;
67
4.06M
68
20.2M
  
for (i = 0; 4.06M
i < n_row20.2M
;
++i16.1M
)
69
16.1M
    mat->row[i] = mat->block.data + i * n_col;
70
4.06M
71
4.06M
  mat->ctx = ctx;
72
4.06M
  isl_ctx_ref(ctx);
73
4.06M
  mat->ref = 1;
74
4.06M
  mat->n_row = n_row;
75
4.06M
  mat->n_col = n_col;
76
4.06M
  mat->max_col = n_col;
77
4.06M
  mat->flags = 0;
78
4.06M
79
4.06M
  return mat;
80
0
error:
81
0
  isl_blk_free(ctx, mat->block);
82
0
  free(mat);
83
0
  return NULL;
84
4.06M
}
85
86
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
87
  unsigned n_row, unsigned n_col)
88
85.2k
{
89
85.2k
  int i;
90
85.2k
  isl_int *old;
91
85.2k
  isl_int **row;
92
85.2k
93
85.2k
  if (!mat)
94
0
    return NULL;
95
85.2k
96
85.2k
  
if (85.2k
mat->max_col >= n_col && 85.2k
mat->n_row >= n_row81.0k
)
{6.58k
97
6.58k
    if (mat->n_col < n_col)
98
199
      mat->n_col = n_col;
99
6.58k
    return mat;
100
6.58k
  }
101
85.2k
102
78.6k
  
if (78.6k
mat->max_col < n_col78.6k
)
{4.21k
103
4.21k
    struct isl_mat *new_mat;
104
4.21k
105
4.21k
    if (n_row < mat->n_row)
106
6
      n_row = mat->n_row;
107
4.21k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
108
4.21k
    if (!new_mat)
109
0
      goto error;
110
45.0k
    
for (i = 0; 4.21k
i < mat->n_row45.0k
;
++i40.8k
)
111
40.8k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
112
4.21k
    isl_mat_free(mat);
113
4.21k
    return new_mat;
114
4.21k
  }
115
78.6k
116
74.4k
  mat = isl_mat_cow(mat);
117
74.4k
  if (!mat)
118
0
    goto error;
119
74.4k
120
74.4k
  old = mat->block.data;
121
74.4k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
122
74.4k
  if (isl_blk_is_error(mat->block))
123
0
    goto error;
124
74.4k
  
row = 74.4k
isl_realloc_array74.4k
(mat->ctx, mat->row, isl_int *, n_row);
125
74.4k
  if (
n_row && 74.4k
!row74.4k
)
126
0
    goto error;
127
74.4k
  mat->row = row;
128
74.4k
129
702k
  for (i = 0; 
i < mat->n_row702k
;
++i628k
)
130
628k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
131
229k
  for (i = mat->n_row; 
i < n_row229k
;
++i155k
)
132
155k
    mat->row[i] = mat->block.data + i * mat->max_col;
133
74.4k
  mat->n_row = n_row;
134
74.4k
  if (mat->n_col < n_col)
135
0
    mat->n_col = n_col;
136
74.4k
137
74.4k
  return mat;
138
0
error:
139
0
  isl_mat_free(mat);
140
0
  return NULL;
141
74.4k
}
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
1.97M
{
146
1.97M
  int i;
147
1.97M
  struct isl_mat *mat;
148
1.97M
149
1.97M
  mat = isl_alloc_type(ctx, struct isl_mat);
150
1.97M
  if (!mat)
151
0
    return NULL;
152
1.97M
  
mat->row = 1.97M
isl_alloc_array1.97M
(ctx, isl_int *, n_row);
153
1.97M
  if (
n_row && 1.97M
!mat->row1.06M
)
154
0
    goto error;
155
5.94M
  
for (i = 0; 1.97M
i < n_row5.94M
;
++i3.97M
)
156
3.97M
    mat->row[i] = row[first_row+i] + first_col;
157
1.97M
  mat->ctx = ctx;
158
1.97M
  isl_ctx_ref(ctx);
159
1.97M
  mat->ref = 1;
160
1.97M
  mat->n_row = n_row;
161
1.97M
  mat->n_col = n_col;
162
1.97M
  mat->block = isl_blk_empty();
163
1.97M
  mat->flags = ISL_MAT_BORROWED;
164
1.97M
  return mat;
165
0
error:
166
0
  free(mat);
167
0
  return NULL;
168
1.97M
}
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
406k
{
173
406k
  if (!mat)
174
0
    return NULL;
175
406k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
176
406k
          first_col, n_col);
177
406k
}
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
243k
{
182
243k
  int i;
183
243k
184
763k
  for (i = 0; 
i < n_row763k
;
++i519k
)
185
519k
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
186
243k
}
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
91.8k
{
191
91.8k
  int i;
192
91.8k
193
286k
  for (i = 0; 
i < n_row286k
;
++i194k
)
194
194k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
195
91.8k
}
196
197
__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
198
1.51M
{
199
1.51M
  if (!mat)
200
0
    return NULL;
201
1.51M
202
1.51M
  mat->ref++;
203
1.51M
  return mat;
204
1.51M
}
205
206
struct isl_mat *isl_mat_dup(struct isl_mat *mat)
207
259k
{
208
259k
  int i;
209
259k
  struct isl_mat *mat2;
210
259k
211
259k
  if (!mat)
212
2.35k
    return NULL;
213
257k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
214
257k
  if (!mat2)
215
0
    return NULL;
216
785k
  
for (i = 0; 257k
i < mat->n_row785k
;
++i527k
)
217
527k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
218
257k
  return mat2;
219
257k
}
220
221
struct isl_mat *isl_mat_cow(struct isl_mat *mat)
222
2.51M
{
223
2.51M
  struct isl_mat *mat2;
224
2.51M
  if (!mat)
225
0
    return NULL;
226
2.51M
227
2.51M
  
if (2.51M
mat->ref == 1 && 2.51M
!2.48M
ISL_F_ISSET2.48M
(mat, ISL_MAT_BORROWED))
228
2.25M
    return mat;
229
2.51M
230
255k
  mat2 = isl_mat_dup(mat);
231
255k
  isl_mat_free(mat);
232
255k
  return mat2;
233
2.51M
}
234
235
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
236
8.59M
{
237
8.59M
  if (!mat)
238
1.04M
    return NULL;
239
8.59M
240
7.55M
  
if (7.55M
--mat->ref > 07.55M
)
241
1.51M
    return NULL;
242
7.55M
243
6.03M
  
if (6.03M
!6.03M
ISL_F_ISSET6.03M
(mat, ISL_MAT_BORROWED))
244
4.06M
    isl_blk_free(mat->ctx, mat->block);
245
6.03M
  isl_ctx_deref(mat->ctx);
246
6.03M
  free(mat->row);
247
6.03M
  free(mat);
248
6.03M
249
6.03M
  return NULL;
250
7.55M
}
251
252
int isl_mat_rows(__isl_keep isl_mat *mat)
253
88.8k
{
254
88.8k
  return mat ? 
mat->n_row88.8k
:
-10
;
255
88.8k
}
256
257
int isl_mat_cols(__isl_keep isl_mat *mat)
258
29.2k
{
259
29.2k
  return mat ? 
mat->n_col29.2k
:
-10
;
260
29.2k
}
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.35k
{
266
4.35k
  if (!mat)
267
0
    return isl_stat_error;
268
4.35k
  
if (4.35k
col < 0 || 4.35k
col >= mat->n_col4.35k
)
269
0
    isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
270
4.35k
      "column out of range", return isl_stat_error);
271
4.35k
  return isl_stat_ok;
272
4.35k
}
273
274
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
275
2.87k
{
276
2.87k
  if (!mat)
277
0
    return -1;
278
2.87k
  
if (2.87k
row < 0 || 2.87k
row >= mat->n_row2.87k
)
279
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
280
2.87k
      return -1);
281
2.87k
  
if (2.87k
check_col(mat, col) < 02.87k
)
282
0
    return -1;
283
2.87k
  
isl_int_set2.87k
(*v, mat->row[row][col]);2.87k
284
2.87k
  return 0;
285
2.87k
}
286
287
/* Extract the element at row "row", oolumn "col" of "mat".
288
 */
289
__isl_give isl_val *isl_mat_get_element_val(__isl_keep isl_mat *mat,
290
  int row, int col)
291
24
{
292
24
  isl_ctx *ctx;
293
24
294
24
  if (!mat)
295
0
    return NULL;
296
24
  ctx = isl_mat_get_ctx(mat);
297
24
  if (
row < 0 || 24
row >= mat->n_row24
)
298
0
    isl_die(ctx, isl_error_invalid, "row out of range",
299
24
      return NULL);
300
24
  
if (24
check_col(mat, col) < 024
)
301
0
    return NULL;
302
24
  return isl_val_int_from_isl_int(ctx, mat->row[row][col]);
303
24
}
304
305
__isl_give isl_mat *isl_mat_set_element(__isl_take isl_mat *mat,
306
  int row, int col, isl_int v)
307
1.39k
{
308
1.39k
  mat = isl_mat_cow(mat);
309
1.39k
  if (!mat)
310
0
    return NULL;
311
1.39k
  
if (1.39k
row < 0 || 1.39k
row >= mat->n_row1.39k
)
312
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
313
1.39k
      goto error);
314
1.39k
  
if (1.39k
check_col(mat, col) < 01.39k
)
315
0
    return isl_mat_free(mat);
316
1.39k
  
isl_int_set1.39k
(mat->row[row][col], v);1.39k
317
1.39k
  return mat;
318
0
error:
319
0
  isl_mat_free(mat);
320
0
  return NULL;
321
1.39k
}
322
323
__isl_give isl_mat *isl_mat_set_element_si(__isl_take isl_mat *mat,
324
  int row, int col, int v)
325
41
{
326
41
  mat = isl_mat_cow(mat);
327
41
  if (!mat)
328
0
    return NULL;
329
41
  
if (41
row < 0 || 41
row >= mat->n_row41
)
330
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
331
41
      goto error);
332
41
  
if (41
check_col(mat, col) < 041
)
333
0
    return isl_mat_free(mat);
334
41
  
isl_int_set_si41
(mat->row[row][col], v);41
335
41
  return mat;
336
0
error:
337
0
  isl_mat_free(mat);
338
0
  return NULL;
339
41
}
340
341
/* Replace the element at row "row", column "col" of "mat" by "v".
342
 */
343
__isl_give isl_mat *isl_mat_set_element_val(__isl_take isl_mat *mat,
344
  int row, int col, __isl_take isl_val *v)
345
0
{
346
0
  if (!v)
347
0
    return isl_mat_free(mat);
348
0
  
if (0
!isl_val_is_int(v)0
)
349
0
    isl_die(isl_val_get_ctx(v), isl_error_invalid,
350
0
      "expecting integer value", goto error);
351
0
  mat = isl_mat_set_element(mat, row, col, v->n);
352
0
  isl_val_free(v);
353
0
  return mat;
354
0
error:
355
0
  isl_val_free(v);
356
0
  return isl_mat_free(mat);
357
0
}
358
359
__isl_give isl_mat *isl_mat_diag(isl_ctx *ctx, unsigned n_row, isl_int d)
360
631k
{
361
631k
  int i;
362
631k
  struct isl_mat *mat;
363
631k
364
631k
  mat = isl_mat_alloc(ctx, n_row, n_row);
365
631k
  if (!mat)
366
0
    return NULL;
367
3.17M
  
for (i = 0; 631k
i < n_row3.17M
;
++i2.53M
)
{2.53M
368
2.53M
    isl_seq_clr(mat->row[i], i);
369
2.53M
    isl_int_set(mat->row[i][i], d);
370
2.53M
    isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
371
2.53M
  }
372
631k
373
631k
  return mat;
374
631k
}
375
376
/* Create an "n_row" by "n_col" matrix with zero elements.
377
 */
378
__isl_give isl_mat *isl_mat_zero(isl_ctx *ctx, unsigned n_row, unsigned n_col)
379
45
{
380
45
  int i;
381
45
  isl_mat *mat;
382
45
383
45
  mat = isl_mat_alloc(ctx, n_row, n_col);
384
45
  if (!mat)
385
0
    return NULL;
386
86
  
for (i = 0; 45
i < n_row86
;
++i41
)
387
41
    isl_seq_clr(mat->row[i], n_col);
388
45
389
45
  return mat;
390
45
}
391
392
__isl_give isl_mat *isl_mat_identity(isl_ctx *ctx, unsigned n_row)
393
631k
{
394
631k
  if (!ctx)
395
0
    return NULL;
396
631k
  return isl_mat_diag(ctx, n_row, ctx->one);
397
631k
}
398
399
/* Is "mat" a (possibly scaled) identity matrix?
400
 */
401
int isl_mat_is_scaled_identity(__isl_keep isl_mat *mat)
402
0
{
403
0
  int i;
404
0
405
0
  if (!mat)
406
0
    return -1;
407
0
  
if (0
mat->n_row != mat->n_col0
)
408
0
    return 0;
409
0
410
0
  
for (i = 0; 0
i < mat->n_row0
;
++i0
)
{0
411
0
    if (isl_seq_first_non_zero(mat->row[i], i) != -1)
412
0
      return 0;
413
0
    
if (0
isl_int_ne0
(mat->row[0][0], mat->row[i][i]))
414
0
      return 0;
415
0
    
if (0
isl_seq_first_non_zero(mat->row[i] + i + 1,0
416
0
              mat->n_col - (i + 1)) != -1)
417
0
      return 0;
418
0
  }
419
0
420
0
  return 1;
421
0
}
422
423
__isl_give isl_vec *isl_mat_vec_product(__isl_take isl_mat *mat,
424
  __isl_take isl_vec *vec)
425
228k
{
426
228k
  int i;
427
228k
  struct isl_vec *prod;
428
228k
429
228k
  if (
!mat || 228k
!vec228k
)
430
0
    goto error;
431
228k
432
228k
  
isl_assert228k
(mat->ctx, mat->n_col == vec->size, goto error);228k
433
228k
434
228k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
435
228k
  if (!prod)
436
0
    goto error;
437
228k
438
1.58M
  
for (i = 0; 228k
i < prod->size1.58M
;
++i1.35M
)
439
1.35M
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
440
1.35M
          &prod->block.data[i]);
441
228k
  isl_mat_free(mat);
442
228k
  isl_vec_free(vec);
443
228k
  return prod;
444
0
error:
445
0
  isl_mat_free(mat);
446
0
  isl_vec_free(vec);
447
0
  return NULL;
448
228k
}
449
450
__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
451
  __isl_take isl_vec *vec)
452
258
{
453
258
  struct isl_mat *vec_mat;
454
258
  int i;
455
258
456
258
  if (
!mat || 258
!vec258
)
457
0
    goto error;
458
258
  vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
459
258
  if (!vec_mat)
460
0
    goto error;
461
2.09k
  
for (i = 0; 258
i < vec->size2.09k
;
++i1.84k
)
462
1.84k
    isl_int_set(vec_mat->row[i][0], vec->el[i]);
463
258
  vec_mat = isl_mat_inverse_product(mat, vec_mat);
464
258
  isl_vec_free(vec);
465
258
  if (!vec_mat)
466
0
    return NULL;
467
258
  vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
468
258
  if (vec)
469
2.09k
    
for (i = 0; 258
i < vec->size2.09k
;
++i1.84k
)
470
1.84k
      isl_int_set(vec->el[i], vec_mat->row[i][0]);
471
258
  isl_mat_free(vec_mat);
472
258
  return vec;
473
0
error:
474
0
  isl_mat_free(mat);
475
0
  isl_vec_free(vec);
476
0
  return NULL;
477
258
}
478
479
__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,
480
  __isl_take isl_mat *mat)
481
14.3k
{
482
14.3k
  int i, j;
483
14.3k
  struct isl_vec *prod;
484
14.3k
485
14.3k
  if (
!mat || 14.3k
!vec14.3k
)
486
0
    goto error;
487
14.3k
488
14.3k
  
isl_assert14.3k
(mat->ctx, mat->n_row == vec->size, goto error);14.3k
489
14.3k
490
14.3k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
491
14.3k
  if (!prod)
492
0
    goto error;
493
14.3k
494
70.6k
  
for (i = 0; 14.3k
i < prod->size70.6k
;
++i56.3k
)
{56.3k
495
56.3k
    isl_int_set_si(prod->el[i], 0);
496
406k
    for (j = 0; 
j < vec->size406k
;
++j349k
)
497
349k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
498
56.3k
  }
499
14.3k
  isl_mat_free(mat);
500
14.3k
  isl_vec_free(vec);
501
14.3k
  return prod;
502
0
error:
503
0
  isl_mat_free(mat);
504
0
  isl_vec_free(vec);
505
0
  return NULL;
506
14.3k
}
507
508
__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
509
  __isl_take isl_mat *right)
510
91.8k
{
511
91.8k
  int i;
512
91.8k
  struct isl_mat *sum;
513
91.8k
514
91.8k
  if (
!left || 91.8k
!right91.8k
)
515
0
    goto error;
516
91.8k
517
91.8k
  
isl_assert91.8k
(left->ctx, left->n_row == right->n_row, goto error);91.8k
518
91.8k
  
isl_assert91.8k
(left->ctx, left->n_row >= 1, goto error);91.8k
519
91.8k
  
isl_assert91.8k
(left->ctx, left->n_col >= 1, goto error);91.8k
520
91.8k
  
isl_assert91.8k
(left->ctx, right->n_col >= 1, goto error);91.8k
521
91.8k
  
isl_assert91.8k
(left->ctx,91.8k
522
91.8k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
523
91.8k
      goto error);
524
91.8k
  
isl_assert91.8k
(left->ctx,91.8k
525
91.8k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
526
91.8k
      goto error);
527
91.8k
528
91.8k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
529
91.8k
  if (!sum)
530
0
    goto error;
531
91.8k
  
isl_int_lcm91.8k
(sum->row[0][0], left->row[0][0], right->row[0][0]);91.8k
532
91.8k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
533
91.8k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
534
91.8k
535
91.8k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
536
606k
  for (i = 1; 
i < sum->n_row606k
;
++i514k
)
{514k
537
514k
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
538
514k
    isl_int_addmul(sum->row[i][0],
539
514k
        right->row[0][0], right->row[i][0]);
540
514k
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
541
514k
        left->n_col-1);
542
514k
    isl_seq_scale(sum->row[i]+left->n_col,
543
514k
        right->row[i]+1, right->row[0][0],
544
514k
        right->n_col-1);
545
514k
  }
546
91.8k
547
91.8k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
548
91.8k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
549
91.8k
  isl_mat_free(left);
550
91.8k
  isl_mat_free(right);
551
91.8k
  return sum;
552
0
error:
553
0
  isl_mat_free(left);
554
0
  isl_mat_free(right);
555
0
  return NULL;
556
91.8k
}
557
558
static void exchange(struct isl_mat *M, struct isl_mat **U,
559
  struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
560
619k
{
561
619k
  int r;
562
3.66M
  for (r = row; 
r < M->n_row3.66M
;
++r3.04M
)
563
3.04M
    isl_int_swap(M->row[r][i], M->row[r][j]);
564
619k
  if (
U619k
)
{598k
565
4.84M
    for (r = 0; 
r < (*U)->n_row4.84M
;
++r4.24M
)
566
4.24M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
567
598k
  }
568
619k
  if (Q)
569
179k
    isl_mat_swap_rows(*Q, i, j);
570
619k
}
571
572
static void subtract(struct isl_mat *M, struct isl_mat **U,
573
  struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
574
653k
{
575
653k
  int r;
576
2.68M
  for (r = row; 
r < M->n_row2.68M
;
++r2.03M
)
577
2.03M
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
578
653k
  if (
U653k
)
{628k
579
4.14M
    for (r = 0; 
r < (*U)->n_row4.14M
;
++r3.51M
)
580
3.51M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
581
628k
  }
582
653k
  if (
Q653k
)
{178k
583
1.20M
    for (r = 0; 
r < (*Q)->n_col1.20M
;
++r1.02M
)
584
1.02M
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
585
178k
  }
586
653k
}
587
588
static void oppose(struct isl_mat *M, struct isl_mat **U,
589
  struct isl_mat **Q, unsigned row, unsigned col)
590
192k
{
591
192k
  int r;
592
718k
  for (r = row; 
r < M->n_row718k
;
++r525k
)
593
525k
    isl_int_neg(M->row[r][col], M->row[r][col]);
594
192k
  if (
U192k
)
{167k
595
1.03M
    for (r = 0; 
r < (*U)->n_row1.03M
;
++r872k
)
596
872k
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
597
167k
  }
598
192k
  if (Q)
599
85.9k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
600
192k
}
601
602
/* Given matrix M, compute
603
 *
604
 *    M U = H
605
 *    M   = H Q
606
 *
607
 * with U and Q unimodular matrices and H a matrix in column echelon form
608
 * such that on each echelon row the entries in the non-echelon column
609
 * are non-negative (if neg == 0) or non-positive (if neg == 1)
610
 * and strictly smaller (in absolute value) than the entries in the echelon
611
 * column.
612
 * If U or Q are NULL, then these matrices are not computed.
613
 */
614
__isl_give isl_mat *isl_mat_left_hermite(__isl_take isl_mat *M, int neg,
615
  __isl_give isl_mat **U, __isl_give isl_mat **Q)
616
317k
{
617
317k
  isl_int c;
618
317k
  int row, col;
619
317k
620
317k
  if (U)
621
289k
    *U = NULL;
622
317k
  if (Q)
623
140k
    *Q = NULL;
624
317k
  if (!M)
625
0
    goto error;
626
317k
  M = isl_mat_cow(M);
627
317k
  if (!M)
628
0
    goto error;
629
317k
  
if (317k
U317k
)
{289k
630
289k
    *U = isl_mat_identity(M->ctx, M->n_col);
631
289k
    if (!*U)
632
0
      goto error;
633
289k
  }
634
317k
  
if (317k
Q317k
)
{140k
635
140k
    *Q = isl_mat_identity(M->ctx, M->n_col);
636
140k
    if (!*Q)
637
0
      goto error;
638
140k
  }
639
317k
640
317k
  col = 0;
641
317k
  isl_int_init(c);
642
1.28M
  for (row = 0; 
row < M->n_row1.28M
;
++row965k
)
{965k
643
965k
    int first, i, off;
644
965k
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
645
965k
    if (first == -1)
646
266k
      continue;
647
698k
    first += col;
648
698k
    if (first != col)
649
354k
      exchange(M, U, Q, row, first, col);
650
698k
    if (isl_int_is_neg(M->row[row][col]))
651
155k
      oppose(M, U, Q, row, col);
652
698k
    first = col+1;
653
998k
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
654
299k
                   M->n_col-first)) != -1) {
655
299k
      first += off;
656
299k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
657
299k
      subtract(M, U, Q, row, col, first, c);
658
299k
      if (
!299k
isl_int_is_zero299k
(M->row[row][first]))
659
44.5k
        exchange(M, U, Q, row, first, col);
660
299k
      else
661
254k
        ++first;
662
299k
    }
663
2.39M
    for (i = 0; 
i < col2.39M
;
++i1.69M
)
{1.69M
664
1.69M
      if (isl_int_is_zero(M->row[row][i]))
665
1.57M
        continue;
666
122k
      
if (122k
neg122k
)
667
0
        isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
668
122k
      else
669
122k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
670
122k
      if (isl_int_is_zero(c))
671
26.0k
        continue;
672
96.1k
      subtract(M, U, Q, row, col, i, c);
673
96.1k
    }
674
698k
    ++col;
675
698k
  }
676
317k
  isl_int_clear(c);
677
317k
678
317k
  return M;
679
0
error:
680
0
  if (
Q0
)
{0
681
0
    isl_mat_free(*Q);
682
0
    *Q = NULL;
683
0
  }
684
0
  if (
U0
)
{0
685
0
    isl_mat_free(*U);
686
0
    *U = NULL;
687
0
  }
688
0
  isl_mat_free(M);
689
0
  return NULL;
690
317k
}
691
692
struct isl_mat *isl_mat_right_kernel(struct isl_mat *mat)
693
0
{
694
0
  int i, rank;
695
0
  struct isl_mat *U = NULL;
696
0
  struct isl_mat *K;
697
0
698
0
  mat = isl_mat_left_hermite(mat, 0, &U, NULL);
699
0
  if (
!mat || 0
!U0
)
700
0
    goto error;
701
0
702
0
  
for (i = 0, rank = 0; 0
rank < mat->n_col0
;
++rank0
)
{0
703
0
    while (
i < mat->n_row && 0
isl_int_is_zero0
(mat->row[i][rank]))
704
0
      ++i;
705
0
    if (i >= mat->n_row)
706
0
      break;
707
0
  }
708
0
  K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);
709
0
  if (!K)
710
0
    goto error;
711
0
  isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
712
0
  isl_mat_free(mat);
713
0
  isl_mat_free(U);
714
0
  return K;
715
0
error:
716
0
  isl_mat_free(mat);
717
0
  isl_mat_free(U);
718
0
  return NULL;
719
0
}
720
721
__isl_give isl_mat *isl_mat_lin_to_aff(__isl_take isl_mat *mat)
722
464k
{
723
464k
  int i;
724
464k
  struct isl_mat *mat2;
725
464k
726
464k
  if (!mat)
727
0
    return NULL;
728
464k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
729
464k
  if (!mat2)
730
0
    goto error;
731
464k
  
isl_int_set_si464k
(mat2->row[0][0], 1);464k
732
464k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
733
2.53M
  for (i = 0; 
i < mat->n_row2.53M
;
++i2.06M
)
{2.06M
734
2.06M
    isl_int_set_si(mat2->row[1+i][0], 0);
735
2.06M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
736
2.06M
  }
737
464k
  isl_mat_free(mat);
738
464k
  return mat2;
739
0
error:
740
0
  isl_mat_free(mat);
741
0
  return NULL;
742
464k
}
743
744
/* Given two matrices M1 and M2, return the block matrix
745
 *
746
 *  [ M1  0  ]
747
 *  [ 0   M2 ]
748
 */
749
__isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1,
750
  __isl_take isl_mat *mat2)
751
73.8k
{
752
73.8k
  int i;
753
73.8k
  isl_mat *mat;
754
73.8k
755
73.8k
  if (
!mat1 || 73.8k
!mat273.8k
)
756
0
    goto error;
757
73.8k
758
73.8k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
759
73.8k
               mat1->n_col + mat2->n_col);
760
73.8k
  if (!mat)
761
0
    goto error;
762
147k
  
for (i = 0; 73.8k
i < mat1->n_row147k
;
++i73.9k
)
{73.9k
763
73.9k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
764
73.9k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
765
73.9k
  }
766
404k
  for (i = 0; 
i < mat2->n_row404k
;
++i330k
)
{330k
767
330k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
768
330k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
769
330k
                mat2->row[i], mat2->n_col);
770
330k
  }
771
73.8k
  isl_mat_free(mat1);
772
73.8k
  isl_mat_free(mat2);
773
73.8k
  return mat;
774
0
error:
775
0
  isl_mat_free(mat1);
776
0
  isl_mat_free(mat2);
777
0
  return NULL;
778
73.8k
}
779
780
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
781
640k
{
782
640k
  int i;
783
640k
784
1.09M
  for (i = 0; 
i < n_row1.09M
;
++i453k
)
785
776k
    
if (776k
!776k
isl_int_is_zero776k
(row[i][col]))
786
322k
      return i;
787
317k
  return -1;
788
640k
}
789
790
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
791
316k
{
792
316k
  int i, min = row_first_non_zero(row, n_row, col);
793
316k
  if (min < 0)
794
0
    return -1;
795
773k
  
for (i = min + 1; 316k
i < n_row773k
;
++i457k
)
{457k
796
457k
    if (isl_int_is_zero(row[i][col]))
797
451k
      continue;
798
5.84k
    
if (5.84k
isl_int_abs_lt5.84k
(row[i][col], row[min][col]))
799
670
      min = i;
800
5.84k
  }
801
316k
  return min;
802
316k
}
803
804
static isl_stat inv_exchange(__isl_keep isl_mat **left,
805
  __isl_keep isl_mat **right, unsigned i, unsigned j)
806
1.11k
{
807
1.11k
  *left = isl_mat_swap_rows(*left, i, j);
808
1.11k
  *right = isl_mat_swap_rows(*right, i, j);
809
1.11k
810
1.11k
  if (
!*left || 1.11k
!*right1.11k
)
811
0
    return isl_stat_error;
812
1.11k
  return isl_stat_ok;
813
1.11k
}
814
815
static void inv_oppose(
816
  struct isl_mat *left, struct isl_mat *right, unsigned row)
817
602
{
818
602
  isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
819
602
  isl_seq_neg(right->row[row], right->row[row], right->n_col);
820
602
}
821
822
static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
823
  unsigned row, unsigned i, isl_int m)
824
6.16k
{
825
6.16k
  isl_int_neg(m, m);
826
6.16k
  isl_seq_combine(left->row[i]+row,
827
6.16k
      left->ctx->one, left->row[i]+row,
828
6.16k
      m, left->row[row]+row,
829
6.16k
      left->n_col-row);
830
6.16k
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
831
6.16k
      m, right->row[row], right->n_col);
832
6.16k
}
833
834
/* Compute inv(left)*right
835
 */
836
__isl_give isl_mat *isl_mat_inverse_product(__isl_take isl_mat *left,
837
  __isl_take isl_mat *right)
838
104k
{
839
104k
  int row;
840
104k
  isl_int a, b;
841
104k
842
104k
  if (
!left || 104k
!right104k
)
843
0
    goto error;
844
104k
845
104k
  
isl_assert104k
(left->ctx, left->n_row == left->n_col, goto error);104k
846
104k
  
isl_assert104k
(left->ctx, left->n_row == right->n_row, goto error);104k
847
104k
848
104k
  
if (104k
left->n_row == 0104k
)
{0
849
0
    isl_mat_free(left);
850
0
    return right;
851
0
  }
852
104k
853
104k
  left = isl_mat_cow(left);
854
104k
  right = isl_mat_cow(right);
855
104k
  if (
!left || 104k
!right104k
)
856
0
    goto error;
857
104k
858
104k
  
isl_int_init104k
(a);104k
859
104k
  isl_int_init(b);
860
420k
  for (row = 0; 
row < left->n_row420k
;
++row316k
)
{316k
861
316k
    int pivot, first, i, off;
862
316k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
863
316k
    if (
pivot < 0316k
)
{0
864
0
      isl_int_clear(a);
865
0
      isl_int_clear(b);
866
0
      isl_assert(left->ctx, pivot >= 0, goto error);
867
0
    }
868
316k
    pivot += row;
869
316k
    if (pivot != row)
870
792
      
if (792
inv_exchange(&left, &right, pivot, row) < 0792
)
871
0
        goto error;
872
316k
    
if (316k
isl_int_is_neg316k
(left->row[row][row]))
873
602
      inv_oppose(left, right, row);
874
316k
    first = row+1;
875
322k
    while ((off = row_first_non_zero(left->row+first,
876
6.16k
          left->n_row-first, row)) != -1) {
877
6.16k
      first += off;
878
6.16k
      isl_int_fdiv_q(a, left->row[first][row],
879
6.16k
          left->row[row][row]);
880
6.16k
      inv_subtract(left, right, row, first, a);
881
6.16k
      if (
!6.16k
isl_int_is_zero6.16k
(left->row[first][row]))
{319
882
319
        if (inv_exchange(&left, &right, row, first) < 0)
883
0
          goto error;
884
5.84k
      } else {
885
5.84k
        ++first;
886
5.84k
      }
887
6.16k
    }
888
774k
    
for (i = 0; 316k
i < row774k
;
++i458k
)
{458k
889
458k
      if (isl_int_is_zero(left->row[i][row]))
890
455k
        continue;
891
2.22k
      
isl_int_gcd2.22k
(a, left->row[row][row], left->row[i][row]);2.22k
892
2.22k
      isl_int_divexact(b, left->row[i][row], a);
893
2.22k
      isl_int_divexact(a, left->row[row][row], a);
894
2.22k
      isl_int_neg(b, b);
895
2.22k
      isl_seq_combine(left->row[i] + i,
896
2.22k
          a, left->row[i] + i,
897
2.22k
          b, left->row[row] + i,
898
2.22k
          left->n_col - i);
899
2.22k
      isl_seq_combine(right->row[i], a, right->row[i],
900
2.22k
          b, right->row[row], right->n_col);
901
2.22k
    }
902
316k
  }
903
104k
  
isl_int_clear104k
(b);104k
904
104k
905
104k
  isl_int_set(a, left->row[0][0]);
906
316k
  for (row = 1; 
row < left->n_row316k
;
++row211k
)
907
211k
    isl_int_lcm(a, a, left->row[row][row]);
908
104k
  if (
isl_int_is_zero104k
(a))
{0
909
0
    isl_int_clear(a);
910
0
    isl_assert(left->ctx, 0, goto error);
911
0
  }
912
420k
  
for (row = 0; 104k
row < left->n_row420k
;
++row316k
)
{316k
913
316k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
914
316k
    if (isl_int_is_one(left->row[row][row]))
915
309k
      continue;
916
7.03k
    isl_seq_scale(right->row[row], right->row[row],
917
7.03k
        left->row[row][row], right->n_col);
918
7.03k
  }
919
104k
  isl_int_clear(a);
920
104k
921
104k
  isl_mat_free(left);
922
104k
  return right;
923
0
error:
924
0
  isl_mat_free(left);
925
0
  isl_mat_free(right);
926
0
  return NULL;
927
104k
}
928
929
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
930
55.2k
{
931
55.2k
  int i;
932
55.2k
933
259k
  for (i = 0; 
i < mat->n_row259k
;
++i204k
)
934
204k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
935
55.2k
}
936
937
void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
938
  isl_int m1, unsigned src1, isl_int m2, unsigned src2)
939
157k
{
940
157k
  int i;
941
157k
  isl_int tmp;
942
157k
943
157k
  isl_int_init(tmp);
944
705k
  for (i = 0; 
i < mat->n_row705k
;
++i547k
)
{547k
945
547k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
946
547k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
947
547k
    isl_int_set(mat->row[i][dst], tmp);
948
547k
  }
949
157k
  isl_int_clear(tmp);
950
157k
}
951
952
__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
953
54.8k
{
954
54.8k
  struct isl_mat *inv;
955
54.8k
  int row;
956
54.8k
  isl_int a, b;
957
54.8k
958
54.8k
  mat = isl_mat_cow(mat);
959
54.8k
  if (!mat)
960
0
    return NULL;
961
54.8k
962
54.8k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
963
54.8k
  inv = isl_mat_cow(inv);
964
54.8k
  if (!inv)
965
0
    goto error;
966
54.8k
967
54.8k
  
isl_int_init54.8k
(a);54.8k
968
54.8k
  isl_int_init(b);
969
210k
  for (row = 0; 
row < mat->n_row210k
;
++row155k
)
{155k
970
155k
    int pivot, first, i, off;
971
155k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
972
155k
    if (
pivot < 0155k
)
{0
973
0
      isl_int_clear(a);
974
0
      isl_int_clear(b);
975
0
      isl_assert(mat->ctx, pivot >= 0, goto error);
976
0
    }
977
155k
    pivot += row;
978
155k
    if (pivot != row)
979
23.8k
      exchange(mat, &inv, NULL, row, pivot, row);
980
155k
    if (isl_int_is_neg(mat->row[row][row]))
981
36.8k
      oppose(mat, &inv, NULL, row, row);
982
155k
    first = row+1;
983
413k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
984
257k
                mat->n_col-first)) != -1) {
985
257k
      first += off;
986
257k
      isl_int_fdiv_q(a, mat->row[row][first],
987
257k
                mat->row[row][row]);
988
257k
      subtract(mat, &inv, NULL, row, row, first, a);
989
257k
      if (
!257k
isl_int_is_zero257k
(mat->row[row][first]))
990
196k
        exchange(mat, &inv, NULL, row, row, first);
991
257k
      else
992
61.7k
        ++first;
993
257k
    }
994
328k
    for (i = 0; 
i < row328k
;
++i172k
)
{172k
995
172k
      if (isl_int_is_zero(mat->row[row][i]))
996
93.9k
        continue;
997
78.8k
      
isl_int_gcd78.8k
(a, mat->row[row][row], mat->row[row][i]);78.8k
998
78.8k
      isl_int_divexact(b, mat->row[row][i], a);
999
78.8k
      isl_int_divexact(a, mat->row[row][row], a);
1000
78.8k
      isl_int_neg(a, a);
1001
78.8k
      isl_mat_col_combine(mat, i, a, i, b, row);
1002
78.8k
      isl_mat_col_combine(inv, i, a, i, b, row);
1003
78.8k
    }
1004
155k
  }
1005
54.8k
  
isl_int_clear54.8k
(b);54.8k
1006
54.8k
1007
54.8k
  isl_int_set(a, mat->row[0][0]);
1008
155k
  for (row = 1; 
row < mat->n_row155k
;
++row100k
)
1009
100k
    isl_int_lcm(a, a, mat->row[row][row]);
1010
54.8k
  if (
isl_int_is_zero54.8k
(a))
{0
1011
0
    isl_int_clear(a);
1012
0
    goto error;
1013
0
  }
1014
210k
  
for (row = 0; 54.8k
row < mat->n_row210k
;
++row155k
)
{155k
1015
155k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1016
155k
    if (isl_int_is_one(mat->row[row][row]))
1017
100k
      continue;
1018
55.2k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1019
55.2k
  }
1020
54.8k
  isl_int_clear(a);
1021
54.8k
1022
54.8k
  isl_mat_free(mat);
1023
54.8k
1024
54.8k
  return inv;
1025
0
error:
1026
0
  isl_mat_free(mat);
1027
0
  isl_mat_free(inv);
1028
0
  return NULL;
1029
54.8k
}
1030
1031
__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1032
3.71k
{
1033
3.71k
  struct isl_mat *transpose = NULL;
1034
3.71k
  int i, j;
1035
3.71k
1036
3.71k
  if (!mat)
1037
0
    return NULL;
1038
3.71k
1039
3.71k
  
if (3.71k
mat->n_col == mat->n_row3.71k
)
{3.71k
1040
3.71k
    mat = isl_mat_cow(mat);
1041
3.71k
    if (!mat)
1042
0
      return NULL;
1043
15.2k
    
for (i = 0; 3.71k
i < mat->n_row15.2k
;
++i11.5k
)
1044
33.5k
      
for (j = i + 1; 11.5k
j < mat->n_col33.5k
;
++j22.0k
)
1045
22.0k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1046
3.71k
    return mat;
1047
3.71k
  }
1048
0
  transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
1049
0
  if (!transpose)
1050
0
    goto error;
1051
0
  
for (i = 0; 0
i < mat->n_row0
;
++i0
)
1052
0
    
for (j = 0; 0
j < mat->n_col0
;
++j0
)
1053
0
      isl_int_set(transpose->row[j][i], mat->row[i][j]);
1054
0
  isl_mat_free(mat);
1055
0
  return transpose;
1056
0
error:
1057
0
  isl_mat_free(mat);
1058
0
  return NULL;
1059
0
}
1060
1061
__isl_give isl_mat *isl_mat_swap_cols(__isl_take isl_mat *mat,
1062
  unsigned i, unsigned j)
1063
395k
{
1064
395k
  int r;
1065
395k
1066
395k
  mat = isl_mat_cow(mat);
1067
395k
  if (!mat)
1068
0
    return NULL;
1069
395k
  
isl_assert395k
(mat->ctx, i < mat->n_col, goto error);395k
1070
395k
  
isl_assert395k
(mat->ctx, j < mat->n_col, goto error);395k
1071
395k
1072
9.55M
  
for (r = 0; 395k
r < mat->n_row9.55M
;
++r9.16M
)
1073
9.16M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1074
395k
  return mat;
1075
0
error:
1076
0
  isl_mat_free(mat);
1077
0
  return NULL;
1078
395k
}
1079
1080
__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1081
  unsigned i, unsigned j)
1082
1.25M
{
1083
1.25M
  isl_int *t;
1084
1.25M
1085
1.25M
  if (!mat)
1086
0
    return NULL;
1087
1.25M
  mat = isl_mat_cow(mat);
1088
1.25M
  if (!mat)
1089
0
    return NULL;
1090
1.25M
  t = mat->row[i];
1091
1.25M
  mat->row[i] = mat->row[j];
1092
1.25M
  mat->row[j] = t;
1093
1.25M
  return mat;
1094
1.25M
}
1095
1096
/* Calculate the product of two matrices.
1097
 *
1098
 * This function is optimized for operand matrices that contain many zeros and
1099
 * skips multiplications where we know one of the operands is zero.
1100
 */
1101
__isl_give isl_mat *isl_mat_product(__isl_take isl_mat *left,
1102
  __isl_take isl_mat *right)
1103
1.43M
{
1104
1.43M
  int i, j, k;
1105
1.43M
  struct isl_mat *prod;
1106
1.43M
1107
1.43M
  if (
!left || 1.43M
!right1.43M
)
1108
0
    goto error;
1109
1.43M
  
isl_assert1.43M
(left->ctx, left->n_col == right->n_row, goto error);1.43M
1110
1.43M
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1111
1.43M
  if (!prod)
1112
0
    goto error;
1113
1.43M
  
if (1.43M
left->n_col == 01.43M
)
{231
1114
359
    for (i = 0; 
i < prod->n_row359
;
++i128
)
1115
128
      isl_seq_clr(prod->row[i], prod->n_col);
1116
231
    isl_mat_free(left);
1117
231
    isl_mat_free(right);
1118
231
    return prod;
1119
231
  }
1120
4.12M
  
for (i = 0; 1.43M
i < prod->n_row4.12M
;
++i2.68M
)
{2.68M
1121
19.2M
    for (j = 0; 
j < prod->n_col19.2M
;
++j16.5M
)
1122
16.5M
      isl_int_mul(prod->row[i][j],
1123
2.68M
            left->row[i][0], right->row[0][j]);
1124
22.1M
    for (k = 1; 
k < left->n_col22.1M
;
++k19.4M
)
{19.4M
1125
19.4M
      if (isl_int_is_zero(left->row[i][k]))
1126
16.2M
        continue;
1127
25.7M
      
for (j = 0; 3.18M
j < prod->n_col25.7M
;
++j22.5M
)
1128
22.5M
        isl_int_addmul(prod->row[i][j],
1129
3.18M
              left->row[i][k], right->row[k][j]);
1130
3.18M
    }
1131
2.68M
  }
1132
1.43M
  isl_mat_free(left);
1133
1.43M
  isl_mat_free(right);
1134
1.43M
  return prod;
1135
0
error:
1136
0
  isl_mat_free(left);
1137
0
  isl_mat_free(right);
1138
0
  return NULL;
1139
1.43M
}
1140
1141
/* Replace the variables x in the rows q by x' given by x = M x',
1142
 * with M the matrix mat.
1143
 *
1144
 * If the number of new variables is greater than the original
1145
 * number of variables, then the rows q have already been
1146
 * preextended.  If the new number is smaller, then the coefficients
1147
 * of the divs, which are not changed, need to be shifted down.
1148
 * The row q may be the equalities, the inequalities or the
1149
 * div expressions.  In the latter case, has_div is true and
1150
 * we need to take into account the extra denominator column.
1151
 */
1152
static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
1153
  unsigned n_div, int has_div, struct isl_mat *mat)
1154
1.23M
{
1155
1.23M
  int i;
1156
1.23M
  struct isl_mat *t;
1157
1.23M
  int e;
1158
1.23M
1159
1.23M
  if (mat->n_col >= mat->n_row)
1160
802k
    e = 0;
1161
1.23M
  else
1162
436k
    e = mat->n_row - mat->n_col;
1163
1.23M
  if (has_div)
1164
412k
    
for (i = 0; 412k
i < n412k
;
++i111
)
1165
111
      isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
1166
1.23M
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1167
1.23M
  t = isl_mat_product(t, mat);
1168
1.23M
  if (!t)
1169
0
    return -1;
1170
2.99M
  
for (i = 0; 1.23M
i < n2.99M
;
++i1.75M
)
{1.75M
1171
1.75M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1172
1.75M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1173
1.75M
          q[i] + has_div + t->n_col + e, n_div);
1174
1.75M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1175
1.75M
  }
1176
1.23M
  isl_mat_free(t);
1177
1.23M
  return 0;
1178
1.23M
}
1179
1180
/* Replace the variables x in bset by x' given by x = M x', with
1181
 * M the matrix mat.
1182
 *
1183
 * If there are fewer variables x' then there are x, then we perform
1184
 * the transformation in place, which means that, in principle,
1185
 * this frees up some extra variables as the number
1186
 * of columns remains constant, but we would have to extend
1187
 * the div array too as the number of rows in this array is assumed
1188
 * to be equal to extra.
1189
 */
1190
struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
1191
  struct isl_mat *mat)
1192
412k
{
1193
412k
  struct isl_ctx *ctx;
1194
412k
1195
412k
  if (
!bset || 412k
!mat412k
)
1196
0
    goto error;
1197
412k
1198
412k
  ctx = bset->ctx;
1199
412k
  bset = isl_basic_set_cow(bset);
1200
412k
  if (!bset)
1201
0
    goto error;
1202
412k
1203
412k
  
isl_assert412k
(ctx, bset->dim->nparam == 0, goto error);412k
1204
412k
  
isl_assert412k
(ctx, 1+bset->dim->n_out == mat->n_row, goto error);412k
1205
412k
  
isl_assert412k
(ctx, mat->n_col > 0, goto error);412k
1206
412k
1207
412k
  
if (412k
mat->n_col > mat->n_row412k
)
{69.3k
1208
69.3k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1209
69.3k
    if (!bset)
1210
0
      goto error;
1211
343k
  } else 
if (343k
mat->n_col < mat->n_row343k
)
{145k
1212
145k
    bset->dim = isl_space_cow(bset->dim);
1213
145k
    if (!bset->dim)
1214
0
      goto error;
1215
145k
    bset->dim->n_out -= mat->n_row - mat->n_col;
1216
145k
  }
1217
412k
1218
412k
  
if (412k
preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,412k
1219
412k
      isl_mat_copy(mat)) < 0)
1220
0
    goto error;
1221
412k
1222
412k
  
if (412k
preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,412k
1223
412k
      isl_mat_copy(mat)) < 0)
1224
0
    goto error;
1225
412k
1226
412k
  
if (412k
preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0412k
)
1227
0
    goto error2;
1228
412k
1229
412k
  
ISL_F_CLR412k
(bset, ISL_BASIC_SET_NO_IMPLICIT);412k
1230
412k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1231
412k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1232
412k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1233
412k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1234
412k
1235
412k
  bset = isl_basic_set_simplify(bset);
1236
412k
  bset = isl_basic_set_finalize(bset);
1237
412k
1238
412k
  return bset;
1239
0
error:
1240
0
  isl_mat_free(mat);
1241
0
error2:
1242
0
  isl_basic_set_free(bset);
1243
0
  return NULL;
1244
0
}
1245
1246
struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
1247
29.0k
{
1248
29.0k
  int i;
1249
29.0k
1250
29.0k
  set = isl_set_cow(set);
1251
29.0k
  if (!set)
1252
0
    goto error;
1253
29.0k
1254
85.0k
  
for (i = 0; 29.0k
i < set->n85.0k
;
++i56.0k
)
{56.0k
1255
56.0k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1256
56.0k
                isl_mat_copy(mat));
1257
56.0k
    if (!set->p[i])
1258
0
      goto error;
1259
56.0k
  }
1260
29.0k
  
if (29.0k
mat->n_col != mat->n_row29.0k
)
{13.6k
1261
13.6k
    set->dim = isl_space_cow(set->dim);
1262
13.6k
    if (!set->dim)
1263
0
      goto error;
1264
13.6k
    set->dim->n_out += mat->n_col;
1265
13.6k
    set->dim->n_out -= mat->n_row;
1266
13.6k
  }
1267
29.0k
  isl_mat_free(mat);
1268
29.0k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1269
29.0k
  return set;
1270
0
error:
1271
0
  isl_set_free(set);
1272
0
  isl_mat_free(mat);
1273
0
  return NULL;
1274
29.0k
}
1275
1276
/* Replace the variables x starting at "first_col" in the rows "rows"
1277
 * of some coefficient matrix by x' with x = M x' with M the matrix mat.
1278
 * That is, replace the corresponding coefficients c by c M.
1279
 */
1280
isl_stat isl_mat_sub_transform(isl_int **row, unsigned n_row,
1281
  unsigned first_col, __isl_take isl_mat *mat)
1282
4.14k
{
1283
4.14k
  int i;
1284
4.14k
  isl_ctx *ctx;
1285
4.14k
  isl_mat *t;
1286
4.14k
1287
4.14k
  if (!mat)
1288
0
    return isl_stat_error;
1289
4.14k
  ctx = isl_mat_get_ctx(mat);
1290
4.14k
  t = isl_mat_sub_alloc6(ctx, row, 0, n_row, first_col, mat->n_row);
1291
4.14k
  t = isl_mat_product(t, mat);
1292
4.14k
  if (!t)
1293
0
    return isl_stat_error;
1294
8.22k
  
for (i = 0; 4.14k
i < n_row8.22k
;
++i4.08k
)
1295
4.08k
    isl_seq_swp_or_cpy(row[i] + first_col, t->row[i], t->n_col);
1296
4.14k
  isl_mat_free(t);
1297
4.14k
  return isl_stat_ok;
1298
4.14k
}
1299
1300
void isl_mat_print_internal(__isl_keep isl_mat *mat, FILE *out, int indent)
1301
0
{
1302
0
  int i, j;
1303
0
1304
0
  if (
!mat0
)
{0
1305
0
    fprintf(out, "%*snull mat\n", indent, "");
1306
0
    return;
1307
0
  }
1308
0
1309
0
  
if (0
mat->n_row == 00
)
1310
0
    fprintf(out, "%*s[]\n", indent, "");
1311
0
1312
0
  for (i = 0; 
i < mat->n_row0
;
++i0
)
{0
1313
0
    if (!i)
1314
0
      fprintf(out, "%*s[[", indent, "");
1315
0
    else
1316
0
      fprintf(out, "%*s[", indent+1, "");
1317
0
    for (j = 0; 
j < mat->n_col0
;
++j0
)
{0
1318
0
      if (j)
1319
0
          fprintf(out, ",");
1320
0
      isl_int_print(out, mat->row[i][j], 0);
1321
0
    }
1322
0
    if (i == mat->n_row-1)
1323
0
      fprintf(out, "]]\n");
1324
0
    else
1325
0
      fprintf(out, "]\n");
1326
0
  }
1327
0
}
1328
1329
void isl_mat_dump(__isl_keep isl_mat *mat)
1330
0
{
1331
0
  isl_mat_print_internal(mat, stderr, 0);
1332
0
}
1333
1334
struct isl_mat *isl_mat_drop_cols(struct isl_mat *mat, unsigned col, unsigned n)
1335
62.6k
{
1336
62.6k
  int r;
1337
62.6k
1338
62.6k
  if (n == 0)
1339
2
    return mat;
1340
62.6k
1341
62.6k
  mat = isl_mat_cow(mat);
1342
62.6k
  if (!mat)
1343
0
    return NULL;
1344
62.6k
1345
62.6k
  
if (62.6k
col != mat->n_col-n62.6k
)
{14.2k
1346
66.3k
    for (r = 0; 
r < mat->n_row66.3k
;
++r52.0k
)
1347
52.0k
      isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1348
52.0k
          mat->n_col - col - n);
1349
14.2k
  }
1350
62.6k
  mat->n_col -= n;
1351
62.6k
  return mat;
1352
62.6k
}
1353
1354
struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n)
1355
78.6k
{
1356
78.6k
  int r;
1357
78.6k
1358
78.6k
  mat = isl_mat_cow(mat);
1359
78.6k
  if (!mat)
1360
0
    return NULL;
1361
78.6k
1362
185k
  
for (r = row; 78.6k
r+n < mat->n_row185k
;
++r106k
)
1363
106k
    mat->row[r] = mat->row[r+n];
1364
78.6k
1365
78.6k
  mat->n_row -= n;
1366
78.6k
  return mat;
1367
78.6k
}
1368
1369
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1370
        unsigned col, unsigned n)
1371
51.5k
{
1372
51.5k
  isl_mat *ext;
1373
51.5k
1374
51.5k
  if (!mat)
1375
0
    return NULL;
1376
51.5k
  
if (51.5k
n == 051.5k
)
1377
87
    return mat;
1378
51.5k
1379
51.4k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1380
51.4k
  if (!ext)
1381
0
    goto error;
1382
51.4k
1383
51.4k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1384
51.4k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1385
51.4k
        col + n, col, mat->n_col - col);
1386
51.4k
1387
51.4k
  isl_mat_free(mat);
1388
51.4k
  return ext;
1389
0
error:
1390
0
  isl_mat_free(mat);
1391
0
  return NULL;
1392
51.4k
}
1393
1394
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1395
  unsigned first, unsigned n)
1396
51.5k
{
1397
51.5k
  int i;
1398
51.5k
1399
51.5k
  if (!mat)
1400
0
    return NULL;
1401
51.5k
  mat = isl_mat_insert_cols(mat, first, n);
1402
51.5k
  if (!mat)
1403
0
    return NULL;
1404
51.5k
1405
54.9k
  
for (i = 0; 51.5k
i < mat->n_row54.9k
;
++i3.36k
)
1406
3.36k
    isl_seq_clr(mat->row[i] + first, n);
1407
51.5k
1408
51.5k
  return mat;
1409
51.5k
}
1410
1411
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1412
3.95k
{
1413
3.95k
  if (!mat)
1414
0
    return NULL;
1415
3.95k
1416
3.95k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1417
3.95k
}
1418
1419
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1420
        unsigned row, unsigned n)
1421
4.37k
{
1422
4.37k
  isl_mat *ext;
1423
4.37k
1424
4.37k
  if (!mat)
1425
0
    return NULL;
1426
4.37k
  
if (4.37k
n == 04.37k
)
1427
85
    return mat;
1428
4.37k
1429
4.28k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1430
4.28k
  if (!ext)
1431
0
    goto error;
1432
4.28k
1433
4.28k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1434
4.28k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1435
4.28k
        mat->n_row - row, 0, 0, mat->n_col);
1436
4.28k
1437
4.28k
  isl_mat_free(mat);
1438
4.28k
  return ext;
1439
0
error:
1440
0
  isl_mat_free(mat);
1441
0
  return NULL;
1442
4.28k
}
1443
1444
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1445
4.37k
{
1446
4.37k
  if (!mat)
1447
0
    return NULL;
1448
4.37k
1449
4.37k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1450
4.37k
}
1451
1452
__isl_give isl_mat *isl_mat_insert_zero_rows(__isl_take isl_mat *mat,
1453
  unsigned row, unsigned n)
1454
0
{
1455
0
  int i;
1456
0
1457
0
  mat = isl_mat_insert_rows(mat, row, n);
1458
0
  if (!mat)
1459
0
    return NULL;
1460
0
  
1461
0
  
for (i = 0; 0
i < n0
;
++i0
)
1462
0
    isl_seq_clr(mat->row[row + i], mat->n_col);
1463
0
1464
0
  return mat;
1465
0
}
1466
1467
__isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n)
1468
0
{
1469
0
  if (!mat)
1470
0
    return NULL;
1471
0
1472
0
  return isl_mat_insert_zero_rows(mat, mat->n_row, n);
1473
0
}
1474
1475
void isl_mat_col_submul(struct isl_mat *mat,
1476
      int dst_col, isl_int f, int src_col)
1477
3.73k
{
1478
3.73k
  int i;
1479
3.73k
1480
16.5k
  for (i = 0; 
i < mat->n_row16.5k
;
++i12.8k
)
1481
12.8k
    isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1482
3.73k
}
1483
1484
void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1485
2
{
1486
2
  int i;
1487
2
1488
2
  if (!mat)
1489
0
    return;
1490
2
1491
7
  
for (i = 0; 2
i < mat->n_row7
;
++i5
)
1492
5
    isl_int_add(mat->row[i][dst_col],
1493
2
          mat->row[i][dst_col], mat->row[i][src_col]);
1494
2
}
1495
1496
void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1497
9.15k
{
1498
9.15k
  int i;
1499
9.15k
1500
31.9k
  for (i = 0; 
i < mat->n_row31.9k
;
++i22.7k
)
1501
22.7k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1502
9.15k
}
1503
1504
/* Add "f" times column "src_col" to column "dst_col" of "mat" and
1505
 * return the result.
1506
 */
1507
__isl_give isl_mat *isl_mat_col_addmul(__isl_take isl_mat *mat, int dst_col,
1508
  isl_int f, int src_col)
1509
6
{
1510
6
  int i;
1511
6
1512
6
  if (
check_col(mat, dst_col) < 0 || 6
check_col(mat, src_col) < 06
)
1513
0
    return isl_mat_free(mat);
1514
6
1515
14
  
for (i = 0; 6
i < mat->n_row14
;
++i8
)
{8
1516
8
    if (isl_int_is_zero(mat->row[i][src_col]))
1517
2
      continue;
1518
6
    mat = isl_mat_cow(mat);
1519
6
    if (!mat)
1520
0
      return NULL;
1521
6
    
isl_int_addmul6
(mat->row[i][dst_col], f, mat->row[i][src_col]);6
1522
6
  }
1523
6
1524
6
  return mat;
1525
6
}
1526
1527
/* Negate column "col" of "mat" and return the result.
1528
 */
1529
__isl_give isl_mat *isl_mat_col_neg(__isl_take isl_mat *mat, int col)
1530
4
{
1531
4
  int i;
1532
4
1533
4
  if (check_col(mat, col) < 0)
1534
0
    return isl_mat_free(mat);
1535
4
1536
10
  
for (i = 0; 4
i < mat->n_row10
;
++i6
)
{6
1537
6
    if (isl_int_is_zero(mat->row[i][col]))
1538
2
      continue;
1539
4
    mat = isl_mat_cow(mat);
1540
4
    if (!mat)
1541
0
      return NULL;
1542
4
    
isl_int_neg4
(mat->row[i][col], mat->row[i][col]);4
1543
4
  }
1544
4
1545
4
  return mat;
1546
4
}
1547
1548
__isl_give isl_mat *isl_mat_unimodular_complete(__isl_take isl_mat *M, int row)
1549
14.5k
{
1550
14.5k
  int r;
1551
14.5k
  struct isl_mat *H = NULL, *Q = NULL;
1552
14.5k
1553
14.5k
  if (!M)
1554
0
    return NULL;
1555
14.5k
1556
14.5k
  
isl_assert14.5k
(M->ctx, M->n_row == M->n_col, goto error);14.5k
1557
14.5k
  M->n_row = row;
1558
14.5k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1559
14.5k
  M->n_row = M->n_col;
1560
14.5k
  if (!H)
1561
0
    goto error;
1562
29.0k
  
for (r = 0; 14.5k
r < row29.0k
;
++r14.5k
)
1563
14.5k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1564
40.9k
  
for (r = row; 14.5k
r < M->n_row40.9k
;
++r26.3k
)
1565
26.3k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1566
14.5k
  isl_mat_free(H);
1567
14.5k
  isl_mat_free(Q);
1568
14.5k
  return M;
1569
0
error:
1570
0
  isl_mat_free(H);
1571
0
  isl_mat_free(Q);
1572
0
  isl_mat_free(M);
1573
0
  return NULL;
1574
14.5k
}
1575
1576
__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1577
  __isl_take isl_mat *bot)
1578
473
{
1579
473
  struct isl_mat *mat;
1580
473
1581
473
  if (
!top || 473
!bot473
)
1582
0
    goto error;
1583
473
1584
473
  
isl_assert473
(top->ctx, top->n_col == bot->n_col, goto error);473
1585
473
  
if (473
top->n_row == 0473
)
{267
1586
267
    isl_mat_free(top);
1587
267
    return bot;
1588
267
  }
1589
206
  
if (206
bot->n_row == 0206
)
{0
1590
0
    isl_mat_free(bot);
1591
0
    return top;
1592
0
  }
1593
206
1594
206
  mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1595
206
  if (!mat)
1596
0
    goto error;
1597
206
  isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1598
206
       0, 0, mat->n_col);
1599
206
  isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1600
206
       0, 0, mat->n_col);
1601
206
  isl_mat_free(top);
1602
206
  isl_mat_free(bot);
1603
206
  return mat;
1604
0
error:
1605
0
  isl_mat_free(top);
1606
0
  isl_mat_free(bot);
1607
0
  return NULL;
1608
206
}
1609
1610
isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1611
6.13k
{
1612
6.13k
  int i;
1613
6.13k
1614
6.13k
  if (
!mat1 || 6.13k
!mat26.13k
)
1615
0
    return isl_bool_error;
1616
6.13k
1617
6.13k
  
if (6.13k
mat1->n_row != mat2->n_row6.13k
)
1618
921
    return isl_bool_false;
1619
6.13k
1620
5.21k
  
if (5.21k
mat1->n_col != mat2->n_col5.21k
)
1621
0
    return isl_bool_false;
1622
5.21k
1623
9.03k
  
for (i = 0; 5.21k
i < mat1->n_row9.03k
;
++i3.82k
)
1624
4.13k
    
if (4.13k
!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col)4.13k
)
1625
313
      return isl_bool_false;
1626
5.21k
1627
4.90k
  return isl_bool_true;
1628
5.21k
}
1629
1630
__isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1631
0
{
1632
0
  struct isl_mat *mat;
1633
0
1634
0
  if (!vec)
1635
0
    return NULL;
1636
0
  mat = isl_mat_alloc(vec->ctx, 1, vec->size);
1637
0
  if (!mat)
1638
0
    goto error;
1639
0
1640
0
  isl_seq_cpy(mat->row[0], vec->el, vec->size);
1641
0
1642
0
  isl_vec_free(vec);
1643
0
  return mat;
1644
0
error:
1645
0
  isl_vec_free(vec);
1646
0
  return NULL;
1647
0
}
1648
1649
/* Return a copy of row "row" of "mat" as an isl_vec.
1650
 */
1651
__isl_give isl_vec *isl_mat_get_row(__isl_keep isl_mat *mat, unsigned row)
1652
24
{
1653
24
  isl_vec *v;
1654
24
1655
24
  if (!mat)
1656
0
    return NULL;
1657
24
  
if (24
row >= mat->n_row24
)
1658
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
1659
24
      return NULL);
1660
24
1661
24
  v = isl_vec_alloc(isl_mat_get_ctx(mat), mat->n_col);
1662
24
  if (!v)
1663
0
    return NULL;
1664
24
  isl_seq_cpy(v->el, mat->row[row], mat->n_col);
1665
24
1666
24
  return v;
1667
24
}
1668
1669
__isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1670
  __isl_take isl_vec *bot)
1671
0
{
1672
0
  return isl_mat_concat(top, isl_mat_from_row_vec(bot));
1673
0
}
1674
1675
__isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1676
  unsigned dst_col, unsigned src_col, unsigned n)
1677
606
{
1678
606
  isl_mat *res;
1679
606
1680
606
  if (!mat)
1681
0
    return NULL;
1682
606
  
if (606
n == 0 || 606
dst_col == src_col606
)
1683
506
    return mat;
1684
606
1685
100
  res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1686
100
  if (!res)
1687
0
    goto error;
1688
100
1689
100
  
if (100
dst_col < src_col100
)
{100
1690
100
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1691
100
         0, 0, dst_col);
1692
100
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1693
100
         dst_col, src_col, n);
1694
100
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1695
100
         dst_col + n, dst_col, src_col - dst_col);
1696
100
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1697
100
         src_col + n, src_col + n,
1698
100
         res->n_col - src_col - n);
1699
0
  } else {
1700
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1701
0
         0, 0, src_col);
1702
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1703
0
         src_col, src_col + n, dst_col - src_col);
1704
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1705
0
         dst_col, src_col, n);
1706
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1707
0
         dst_col + n, dst_col + n,
1708
0
         res->n_col - dst_col - n);
1709
0
  }
1710
100
  isl_mat_free(mat);
1711
100
1712
100
  return res;
1713
0
error:
1714
0
  isl_mat_free(mat);
1715
0
  return NULL;
1716
100
}
1717
1718
/* Return the gcd of the elements in row "row" of "mat" in *gcd.
1719
 * Return isl_stat_ok on success and isl_stat_error on failure.
1720
 */
1721
isl_stat isl_mat_row_gcd(__isl_keep isl_mat *mat, int row, isl_int *gcd)
1722
7
{
1723
7
  if (!mat)
1724
0
    return isl_stat_error;
1725
7
1726
7
  
if (7
row < 0 || 7
row >= mat->n_row7
)
1727
0
    isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
1728
7
      "row out of range", return isl_stat_error);
1729
7
  isl_seq_gcd(mat->row[row], mat->n_col, gcd);
1730
7
1731
7
  return isl_stat_ok;
1732
7
}
1733
1734
void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1735
779
{
1736
779
  int i;
1737
779
  isl_int g;
1738
779
1739
779
  isl_int_set_si(*gcd, 0);
1740
779
  if (!mat)
1741
0
    return;
1742
779
1743
779
  
isl_int_init779
(g);779
1744
4.90k
  for (i = 0; 
i < mat->n_row4.90k
;
++i4.12k
)
{4.12k
1745
4.12k
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1746
4.12k
    isl_int_gcd(*gcd, *gcd, g);
1747
4.12k
  }
1748
779
  isl_int_clear(g);
1749
779
}
1750
1751
/* Return the result of scaling "mat" by a factor of "m".
1752
 */
1753
__isl_give isl_mat *isl_mat_scale(__isl_take isl_mat *mat, isl_int m)
1754
0
{
1755
0
  int i;
1756
0
1757
0
  if (isl_int_is_one(m))
1758
0
    return mat;
1759
0
1760
0
  mat = isl_mat_cow(mat);
1761
0
  if (!mat)
1762
0
    return NULL;
1763
0
1764
0
  
for (i = 0; 0
i < mat->n_row0
;
++i0
)
1765
0
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
1766
0
1767
0
  return mat;
1768
0
}
1769
1770
__isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1771
779
{
1772
779
  int i;
1773
779
1774
779
  if (isl_int_is_one(m))
1775
0
    return mat;
1776
779
1777
779
  mat = isl_mat_cow(mat);
1778
779
  if (!mat)
1779
0
    return NULL;
1780
779
1781
4.90k
  
for (i = 0; 779
i < mat->n_row4.90k
;
++i4.12k
)
1782
4.12k
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1783
779
1784
779
  return mat;
1785
779
}
1786
1787
__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row,
1788
  isl_int m)
1789
10
{
1790
10
  if (isl_int_is_one(m))
1791
0
    return mat;
1792
10
1793
10
  mat = isl_mat_cow(mat);
1794
10
  if (!mat)
1795
0
    return NULL;
1796
10
1797
10
  isl_seq_scale_down(mat->row[row], mat->row[row], m, mat->n_col);
1798
10
1799
10
  return mat;
1800
10
}
1801
1802
__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1803
779
{
1804
779
  isl_int gcd;
1805
779
1806
779
  if (!mat)
1807
0
    return NULL;
1808
779
1809
779
  
isl_int_init779
(gcd);779
1810
779
  isl_mat_gcd(mat, &gcd);
1811
779
  mat = isl_mat_scale_down(mat, gcd);
1812
779
  isl_int_clear(gcd);
1813
779
1814
779
  return mat;
1815
779
}
1816
1817
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1818
1.28k
{
1819
1.28k
  mat = isl_mat_cow(mat);
1820
1.28k
  if (!mat)
1821
0
    return NULL;
1822
1.28k
1823
1.28k
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
1824
1.28k
1825
1.28k
  return mat;
1826
1.28k
}
1827
1828
/* Number of initial non-zero columns.
1829
 */
1830
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
1831
951
{
1832
951
  int i;
1833
951
1834
951
  if (!mat)
1835
0
    return -1;
1836
951
1837
1.37k
  
for (i = 0; 951
i < mat->n_col1.37k
;
++i419
)
1838
1.18k
    
if (1.18k
row_first_non_zero(mat->row, mat->n_row, i) < 01.18k
)
1839
768
      break;
1840
951
1841
951
  return i;
1842
951
}