Coverage Report

Created: 2017-04-29 12:21

/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.04M
{
24
2.04M
  return mat ? mat->ctx : NULL;
25
2.04M
}
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
3.65M
{
53
3.65M
  int i;
54
3.65M
  struct isl_mat *mat;
55
3.65M
56
3.65M
  mat = isl_alloc_type(ctx, struct isl_mat);
57
3.65M
  if (!mat)
58
0
    return NULL;
59
3.65M
60
3.65M
  mat->row = NULL;
61
3.65M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
62
3.65M
  if (isl_blk_is_error(mat->block))
63
0
    goto error;
64
3.65M
  
mat->row = 3.65M
isl_alloc_array3.65M
(ctx, isl_int *, n_row);
65
3.65M
  if (
n_row && 3.65M
!mat->row2.65M
)
66
0
    goto error;
67
3.65M
68
18.4M
  
for (i = 0; 3.65M
i < n_row18.4M
;
++i14.7M
)
69
14.7M
    mat->row[i] = mat->block.data + i * n_col;
70
3.65M
71
3.65M
  mat->ctx = ctx;
72
3.65M
  isl_ctx_ref(ctx);
73
3.65M
  mat->ref = 1;
74
3.65M
  mat->n_row = n_row;
75
3.65M
  mat->n_col = n_col;
76
3.65M
  mat->max_col = n_col;
77
3.65M
  mat->flags = 0;
78
3.65M
79
3.65M
  return mat;
80
0
error:
81
0
  isl_blk_free(ctx, mat->block);
82
0
  free(mat);
83
0
  return NULL;
84
3.65M
}
85
86
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
87
  unsigned n_row, unsigned n_col)
88
70.0k
{
89
70.0k
  int i;
90
70.0k
  isl_int *old;
91
70.0k
  isl_int **row;
92
70.0k
93
70.0k
  if (!mat)
94
0
    return NULL;
95
70.0k
96
70.0k
  
if (70.0k
mat->max_col >= n_col && 70.0k
mat->n_row >= n_row67.3k
)
{5.26k
97
5.26k
    if (mat->n_col < n_col)
98
189
      mat->n_col = n_col;
99
5.26k
    return mat;
100
5.26k
  }
101
70.0k
102
64.8k
  
if (64.8k
mat->max_col < n_col64.8k
)
{2.72k
103
2.72k
    struct isl_mat *new_mat;
104
2.72k
105
2.72k
    if (n_row < mat->n_row)
106
2
      n_row = mat->n_row;
107
2.72k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
108
2.72k
    if (!new_mat)
109
0
      goto error;
110
21.1k
    
for (i = 0; 2.72k
i < mat->n_row21.1k
;
++i18.4k
)
111
18.4k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
112
2.72k
    isl_mat_free(mat);
113
2.72k
    return new_mat;
114
2.72k
  }
115
64.8k
116
62.0k
  mat = isl_mat_cow(mat);
117
62.0k
  if (!mat)
118
0
    goto error;
119
62.0k
120
62.0k
  old = mat->block.data;
121
62.0k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
122
62.0k
  if (isl_blk_is_error(mat->block))
123
0
    goto error;
124
62.0k
  
row = 62.0k
isl_realloc_array62.0k
(mat->ctx, mat->row, isl_int *, n_row);
125
62.0k
  if (
n_row && 62.0k
!row62.0k
)
126
0
    goto error;
127
62.0k
  mat->row = row;
128
62.0k
129
544k
  for (i = 0; 
i < mat->n_row544k
;
++i482k
)
130
482k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
131
191k
  for (i = mat->n_row; 
i < n_row191k
;
++i129k
)
132
129k
    mat->row[i] = mat->block.data + i * mat->max_col;
133
62.0k
  mat->n_row = n_row;
134
62.0k
  if (mat->n_col < n_col)
135
0
    mat->n_col = n_col;
136
62.0k
137
62.0k
  return mat;
138
0
error:
139
0
  isl_mat_free(mat);
140
0
  return NULL;
141
62.0k
}
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.76M
{
146
1.76M
  int i;
147
1.76M
  struct isl_mat *mat;
148
1.76M
149
1.76M
  mat = isl_alloc_type(ctx, struct isl_mat);
150
1.76M
  if (!mat)
151
0
    return NULL;
152
1.76M
  
mat->row = 1.76M
isl_alloc_array1.76M
(ctx, isl_int *, n_row);
153
1.76M
  if (
n_row && 1.76M
!mat->row945k
)
154
0
    goto error;
155
5.42M
  
for (i = 0; 1.76M
i < n_row5.42M
;
++i3.65M
)
156
3.65M
    mat->row[i] = row[first_row+i] + first_col;
157
1.76M
  mat->ctx = ctx;
158
1.76M
  isl_ctx_ref(ctx);
159
1.76M
  mat->ref = 1;
160
1.76M
  mat->n_row = n_row;
161
1.76M
  mat->n_col = n_col;
162
1.76M
  mat->block = isl_blk_empty();
163
1.76M
  mat->flags = ISL_MAT_BORROWED;
164
1.76M
  return mat;
165
0
error:
166
0
  free(mat);
167
0
  return NULL;
168
1.76M
}
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
349k
{
173
349k
  if (!mat)
174
0
    return NULL;
175
349k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
176
349k
          first_col, n_col);
177
349k
}
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
220k
{
182
220k
  int i;
183
220k
184
696k
  for (i = 0; 
i < n_row696k
;
++i475k
)
185
475k
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
186
220k
}
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
79.4k
{
191
79.4k
  int i;
192
79.4k
193
251k
  for (i = 0; 
i < n_row251k
;
++i171k
)
194
171k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
195
79.4k
}
196
197
__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
198
1.38M
{
199
1.38M
  if (!mat)
200
0
    return NULL;
201
1.38M
202
1.38M
  mat->ref++;
203
1.38M
  return mat;
204
1.38M
}
205
206
struct isl_mat *isl_mat_dup(struct isl_mat *mat)
207
232k
{
208
232k
  int i;
209
232k
  struct isl_mat *mat2;
210
232k
211
232k
  if (!mat)
212
1.96k
    return NULL;
213
230k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
214
230k
  if (!mat2)
215
0
    return NULL;
216
712k
  
for (i = 0; 230k
i < mat->n_row712k
;
++i481k
)
217
481k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
218
230k
  return mat2;
219
230k
}
220
221
struct isl_mat *isl_mat_cow(struct isl_mat *mat)
222
2.18M
{
223
2.18M
  struct isl_mat *mat2;
224
2.18M
  if (!mat)
225
0
    return NULL;
226
2.18M
227
2.18M
  
if (2.18M
mat->ref == 1 && 2.18M
!2.15M
ISL_F_ISSET2.15M
(mat, ISL_MAT_BORROWED))
228
1.95M
    return mat;
229
2.18M
230
228k
  mat2 = isl_mat_dup(mat);
231
228k
  isl_mat_free(mat);
232
228k
  return mat2;
233
2.18M
}
234
235
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
236
7.73M
{
237
7.73M
  if (!mat)
238
934k
    return NULL;
239
7.73M
240
6.80M
  
if (6.80M
--mat->ref > 06.80M
)
241
1.38M
    return NULL;
242
6.80M
243
5.41M
  
if (5.41M
!5.41M
ISL_F_ISSET5.41M
(mat, ISL_MAT_BORROWED))
244
3.65M
    isl_blk_free(mat->ctx, mat->block);
245
5.41M
  isl_ctx_deref(mat->ctx);
246
5.41M
  free(mat->row);
247
5.41M
  free(mat);
248
5.41M
249
5.41M
  return NULL;
250
6.80M
}
251
252
int isl_mat_rows(__isl_keep isl_mat *mat)
253
81.3k
{
254
81.3k
  return mat ? 
mat->n_row81.3k
:
-10
;
255
81.3k
}
256
257
int isl_mat_cols(__isl_keep isl_mat *mat)
258
24.6k
{
259
24.6k
  return mat ? 
mat->n_col24.6k
:
-10
;
260
24.6k
}
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.38k
{
266
4.38k
  if (!mat)
267
0
    return isl_stat_error;
268
4.38k
  
if (4.38k
col < 0 || 4.38k
col >= mat->n_col4.38k
)
269
0
    isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
270
4.38k
      "column out of range", return isl_stat_error);
271
4.38k
  return isl_stat_ok;
272
4.38k
}
273
274
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
275
2.90k
{
276
2.90k
  if (!mat)
277
0
    return -1;
278
2.90k
  
if (2.90k
row < 0 || 2.90k
row >= mat->n_row2.90k
)
279
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
280
2.90k
      return -1);
281
2.90k
  
if (2.90k
check_col(mat, col) < 02.90k
)
282
0
    return -1;
283
2.90k
  
isl_int_set2.90k
(*v, mat->row[row][col]);2.90k
284
2.90k
  return 0;
285
2.90k
}
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
567k
{
361
567k
  int i;
362
567k
  struct isl_mat *mat;
363
567k
364
567k
  mat = isl_mat_alloc(ctx, n_row, n_row);
365
567k
  if (!mat)
366
0
    return NULL;
367
2.90M
  
for (i = 0; 567k
i < n_row2.90M
;
++i2.34M
)
{2.34M
368
2.34M
    isl_seq_clr(mat->row[i], i);
369
2.34M
    isl_int_set(mat->row[i][i], d);
370
2.34M
    isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
371
2.34M
  }
372
567k
373
567k
  return mat;
374
567k
}
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
567k
{
394
567k
  if (!ctx)
395
0
    return NULL;
396
567k
  return isl_mat_diag(ctx, n_row, ctx->one);
397
567k
}
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
205k
{
426
205k
  int i;
427
205k
  struct isl_vec *prod;
428
205k
429
205k
  if (
!mat || 205k
!vec205k
)
430
0
    goto error;
431
205k
432
205k
  
isl_assert205k
(mat->ctx, mat->n_col == vec->size, goto error);205k
433
205k
434
205k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
435
205k
  if (!prod)
436
0
    goto error;
437
205k
438
1.47M
  
for (i = 0; 205k
i < prod->size1.47M
;
++i1.27M
)
439
1.27M
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
440
1.27M
          &prod->block.data[i]);
441
205k
  isl_mat_free(mat);
442
205k
  isl_vec_free(vec);
443
205k
  return prod;
444
0
error:
445
0
  isl_mat_free(mat);
446
0
  isl_vec_free(vec);
447
0
  return NULL;
448
205k
}
449
450
__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
451
  __isl_take isl_vec *vec)
452
245
{
453
245
  struct isl_mat *vec_mat;
454
245
  int i;
455
245
456
245
  if (
!mat || 245
!vec245
)
457
0
    goto error;
458
245
  vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
459
245
  if (!vec_mat)
460
0
    goto error;
461
1.99k
  
for (i = 0; 245
i < vec->size1.99k
;
++i1.75k
)
462
1.75k
    isl_int_set(vec_mat->row[i][0], vec->el[i]);
463
245
  vec_mat = isl_mat_inverse_product(mat, vec_mat);
464
245
  isl_vec_free(vec);
465
245
  if (!vec_mat)
466
0
    return NULL;
467
245
  vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
468
245
  if (vec)
469
1.99k
    
for (i = 0; 245
i < vec->size1.99k
;
++i1.75k
)
470
1.75k
      isl_int_set(vec->el[i], vec_mat->row[i][0]);
471
245
  isl_mat_free(vec_mat);
472
245
  return vec;
473
0
error:
474
0
  isl_mat_free(mat);
475
0
  isl_vec_free(vec);
476
0
  return NULL;
477
245
}
478
479
__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,
480
  __isl_take isl_mat *mat)
481
11.8k
{
482
11.8k
  int i, j;
483
11.8k
  struct isl_vec *prod;
484
11.8k
485
11.8k
  if (
!mat || 11.8k
!vec11.8k
)
486
0
    goto error;
487
11.8k
488
11.8k
  
isl_assert11.8k
(mat->ctx, mat->n_row == vec->size, goto error);11.8k
489
11.8k
490
11.8k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
491
11.8k
  if (!prod)
492
0
    goto error;
493
11.8k
494
60.7k
  
for (i = 0; 11.8k
i < prod->size60.7k
;
++i48.8k
)
{48.8k
495
48.8k
    isl_int_set_si(prod->el[i], 0);
496
370k
    for (j = 0; 
j < vec->size370k
;
++j322k
)
497
322k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
498
48.8k
  }
499
11.8k
  isl_mat_free(mat);
500
11.8k
  isl_vec_free(vec);
501
11.8k
  return prod;
502
0
error:
503
0
  isl_mat_free(mat);
504
0
  isl_vec_free(vec);
505
0
  return NULL;
506
11.8k
}
507
508
__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
509
  __isl_take isl_mat *right)
510
79.4k
{
511
79.4k
  int i;
512
79.4k
  struct isl_mat *sum;
513
79.4k
514
79.4k
  if (
!left || 79.4k
!right79.4k
)
515
0
    goto error;
516
79.4k
517
79.4k
  
isl_assert79.4k
(left->ctx, left->n_row == right->n_row, goto error);79.4k
518
79.4k
  
isl_assert79.4k
(left->ctx, left->n_row >= 1, goto error);79.4k
519
79.4k
  
isl_assert79.4k
(left->ctx, left->n_col >= 1, goto error);79.4k
520
79.4k
  
isl_assert79.4k
(left->ctx, right->n_col >= 1, goto error);79.4k
521
79.4k
  
isl_assert79.4k
(left->ctx,79.4k
522
79.4k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
523
79.4k
      goto error);
524
79.4k
  
isl_assert79.4k
(left->ctx,79.4k
525
79.4k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
526
79.4k
      goto error);
527
79.4k
528
79.4k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
529
79.4k
  if (!sum)
530
0
    goto error;
531
79.4k
  
isl_int_lcm79.4k
(sum->row[0][0], left->row[0][0], right->row[0][0]);79.4k
532
79.4k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
533
79.4k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
534
79.4k
535
79.4k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
536
548k
  for (i = 1; 
i < sum->n_row548k
;
++i468k
)
{468k
537
468k
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
538
468k
    isl_int_addmul(sum->row[i][0],
539
468k
        right->row[0][0], right->row[i][0]);
540
468k
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
541
468k
        left->n_col-1);
542
468k
    isl_seq_scale(sum->row[i]+left->n_col,
543
468k
        right->row[i]+1, right->row[0][0],
544
468k
        right->n_col-1);
545
468k
  }
546
79.4k
547
79.4k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
548
79.4k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
549
79.4k
  isl_mat_free(left);
550
79.4k
  isl_mat_free(right);
551
79.4k
  return sum;
552
0
error:
553
0
  isl_mat_free(left);
554
0
  isl_mat_free(right);
555
0
  return NULL;
556
79.4k
}
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
551k
{
561
551k
  int r;
562
3.23M
  for (r = row; 
r < M->n_row3.23M
;
++r2.68M
)
563
2.68M
    isl_int_swap(M->row[r][i], M->row[r][j]);
564
551k
  if (
U551k
)
{534k
565
4.37M
    for (r = 0; 
r < (*U)->n_row4.37M
;
++r3.84M
)
566
3.84M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
567
534k
  }
568
551k
  if (Q)
569
149k
    isl_mat_swap_rows(*Q, i, j);
570
551k
}
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
559k
{
575
559k
  int r;
576
2.02M
  for (r = row; 
r < M->n_row2.02M
;
++r1.46M
)
577
1.46M
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
578
559k
  if (
U559k
)
{538k
579
3.44M
    for (r = 0; 
r < (*U)->n_row3.44M
;
++r2.90M
)
580
2.90M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
581
538k
  }
582
559k
  if (
Q559k
)
{126k
583
826k
    for (r = 0; 
r < (*Q)->n_col826k
;
++r699k
)
584
699k
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
585
126k
  }
586
559k
}
587
588
static void oppose(struct isl_mat *M, struct isl_mat **U,
589
  struct isl_mat **Q, unsigned row, unsigned col)
590
170k
{
591
170k
  int r;
592
627k
  for (r = row; 
r < M->n_row627k
;
++r456k
)
593
456k
    isl_int_neg(M->row[r][col], M->row[r][col]);
594
170k
  if (
U170k
)
{149k
595
942k
    for (r = 0; 
r < (*U)->n_row942k
;
++r793k
)
596
793k
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
597
149k
  }
598
170k
  if (Q)
599
74.6k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
600
170k
}
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
281k
{
617
281k
  isl_int c;
618
281k
  int row, col;
619
281k
620
281k
  if (U)
621
257k
    *U = NULL;
622
281k
  if (Q)
623
125k
    *Q = NULL;
624
281k
  if (!M)
625
0
    goto error;
626
281k
  M = isl_mat_cow(M);
627
281k
  if (!M)
628
0
    goto error;
629
281k
  
if (281k
U281k
)
{257k
630
257k
    *U = isl_mat_identity(M->ctx, M->n_col);
631
257k
    if (!*U)
632
0
      goto error;
633
257k
  }
634
281k
  
if (281k
Q281k
)
{125k
635
125k
    *Q = isl_mat_identity(M->ctx, M->n_col);
636
125k
    if (!*Q)
637
0
      goto error;
638
125k
  }
639
281k
640
281k
  col = 0;
641
281k
  isl_int_init(c);
642
1.16M
  for (row = 0; 
row < M->n_row1.16M
;
++row879k
)
{879k
643
879k
    int first, i, off;
644
879k
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
645
879k
    if (first == -1)
646
247k
      continue;
647
632k
    first += col;
648
632k
    if (first != col)
649
320k
      exchange(M, U, Q, row, first, col);
650
632k
    if (isl_int_is_neg(M->row[row][col]))
651
136k
      oppose(M, U, Q, row, col);
652
632k
    first = col+1;
653
864k
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
654
232k
                   M->n_col-first)) != -1) {
655
232k
      first += off;
656
232k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
657
232k
      subtract(M, U, Q, row, col, first, c);
658
232k
      if (
!232k
isl_int_is_zero232k
(M->row[row][first]))
659
13.2k
        exchange(M, U, Q, row, first, col);
660
232k
      else
661
218k
        ++first;
662
232k
    }
663
2.24M
    for (i = 0; 
i < col2.24M
;
++i1.61M
)
{1.61M
664
1.61M
      if (isl_int_is_zero(M->row[row][i]))
665
1.52M
        continue;
666
93.5k
      
if (93.5k
neg93.5k
)
667
0
        isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
668
93.5k
      else
669
93.5k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
670
93.5k
      if (isl_int_is_zero(c))
671
22.5k
        continue;
672
71.0k
      subtract(M, U, Q, row, col, i, c);
673
71.0k
    }
674
632k
    ++col;
675
632k
  }
676
281k
  isl_int_clear(c);
677
281k
678
281k
  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
281k
}
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
406k
{
723
406k
  int i;
724
406k
  struct isl_mat *mat2;
725
406k
726
406k
  if (!mat)
727
0
    return NULL;
728
406k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
729
406k
  if (!mat2)
730
0
    goto error;
731
406k
  
isl_int_set_si406k
(mat2->row[0][0], 1);406k
732
406k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
733
2.31M
  for (i = 0; 
i < mat->n_row2.31M
;
++i1.90M
)
{1.90M
734
1.90M
    isl_int_set_si(mat2->row[1+i][0], 0);
735
1.90M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
736
1.90M
  }
737
406k
  isl_mat_free(mat);
738
406k
  return mat2;
739
0
error:
740
0
  isl_mat_free(mat);
741
0
  return NULL;
742
406k
}
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
66.8k
{
752
66.8k
  int i;
753
66.8k
  isl_mat *mat;
754
66.8k
755
66.8k
  if (
!mat1 || 66.8k
!mat266.8k
)
756
0
    goto error;
757
66.8k
758
66.8k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
759
66.8k
               mat1->n_col + mat2->n_col);
760
66.8k
  if (!mat)
761
0
    goto error;
762
133k
  
for (i = 0; 66.8k
i < mat1->n_row133k
;
++i66.9k
)
{66.9k
763
66.9k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
764
66.9k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
765
66.9k
  }
766
379k
  for (i = 0; 
i < mat2->n_row379k
;
++i312k
)
{312k
767
312k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
768
312k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
769
312k
                mat2->row[i], mat2->n_col);
770
312k
  }
771
66.8k
  isl_mat_free(mat1);
772
66.8k
  isl_mat_free(mat2);
773
66.8k
  return mat;
774
0
error:
775
0
  isl_mat_free(mat1);
776
0
  isl_mat_free(mat2);
777
0
  return NULL;
778
66.8k
}
779
780
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
781
559k
{
782
559k
  int i;
783
559k
784
968k
  for (i = 0; 
i < n_row968k
;
++i409k
)
785
691k
    
if (691k
!691k
isl_int_is_zero691k
(row[i][col]))
786
282k
      return i;
787
276k
  return -1;
788
559k
}
789
790
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
791
276k
{
792
276k
  int i, min = row_first_non_zero(row, n_row, col);
793
276k
  if (min < 0)
794
0
    return -1;
795
689k
  
for (i = min + 1; 276k
i < n_row689k
;
++i413k
)
{413k
796
413k
    if (isl_int_is_zero(row[i][col]))
797
407k
      continue;
798
5.48k
    
if (5.48k
isl_int_abs_lt5.48k
(row[i][col], row[min][col]))
799
666
      min = i;
800
5.48k
  }
801
276k
  return min;
802
276k
}
803
804
static isl_stat inv_exchange(__isl_keep isl_mat **left,
805
  __isl_keep isl_mat **right, unsigned i, unsigned j)
806
1.06k
{
807
1.06k
  *left = isl_mat_swap_rows(*left, i, j);
808
1.06k
  *right = isl_mat_swap_rows(*right, i, j);
809
1.06k
810
1.06k
  if (
!*left || 1.06k
!*right1.06k
)
811
0
    return isl_stat_error;
812
1.06k
  return isl_stat_ok;
813
1.06k
}
814
815
static void inv_oppose(
816
  struct isl_mat *left, struct isl_mat *right, unsigned row)
817
601
{
818
601
  isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
819
601
  isl_seq_neg(right->row[row], right->row[row], right->n_col);
820
601
}
821
822
static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
823
  unsigned row, unsigned i, isl_int m)
824
5.77k
{
825
5.77k
  isl_int_neg(m, m);
826
5.77k
  isl_seq_combine(left->row[i]+row,
827
5.77k
      left->ctx->one, left->row[i]+row,
828
5.77k
      m, left->row[row]+row,
829
5.77k
      left->n_col-row);
830
5.77k
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
831
5.77k
      m, right->row[row], right->n_col);
832
5.77k
}
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
89.8k
{
839
89.8k
  int row;
840
89.8k
  isl_int a, b;
841
89.8k
842
89.8k
  if (
!left || 89.8k
!right89.8k
)
843
0
    goto error;
844
89.8k
845
89.8k
  
isl_assert89.8k
(left->ctx, left->n_row == left->n_col, goto error);89.8k
846
89.8k
  
isl_assert89.8k
(left->ctx, left->n_row == right->n_row, goto error);89.8k
847
89.8k
848
89.8k
  
if (89.8k
left->n_row == 089.8k
)
{0
849
0
    isl_mat_free(left);
850
0
    return right;
851
0
  }
852
89.8k
853
89.8k
  left = isl_mat_cow(left);
854
89.8k
  right = isl_mat_cow(right);
855
89.8k
  if (
!left || 89.8k
!right89.8k
)
856
0
    goto error;
857
89.8k
858
89.8k
  
isl_int_init89.8k
(a);89.8k
859
89.8k
  isl_int_init(b);
860
365k
  for (row = 0; 
row < left->n_row365k
;
++row276k
)
{276k
861
276k
    int pivot, first, i, off;
862
276k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
863
276k
    if (
pivot < 0276k
)
{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
276k
    pivot += row;
869
276k
    if (pivot != row)
870
775
      
if (775
inv_exchange(&left, &right, pivot, row) < 0775
)
871
0
        goto error;
872
276k
    
if (276k
isl_int_is_neg276k
(left->row[row][row]))
873
601
      inv_oppose(left, right, row);
874
276k
    first = row+1;
875
281k
    while ((off = row_first_non_zero(left->row+first,
876
5.77k
          left->n_row-first, row)) != -1) {
877
5.77k
      first += off;
878
5.77k
      isl_int_fdiv_q(a, left->row[first][row],
879
5.77k
          left->row[row][row]);
880
5.77k
      inv_subtract(left, right, row, first, a);
881
5.77k
      if (
!5.77k
isl_int_is_zero5.77k
(left->row[first][row]))
{287
882
287
        if (inv_exchange(&left, &right, row, first) < 0)
883
0
          goto error;
884
5.48k
      } else {
885
5.48k
        ++first;
886
5.48k
      }
887
5.77k
    }
888
690k
    
for (i = 0; 276k
i < row690k
;
++i414k
)
{414k
889
414k
      if (isl_int_is_zero(left->row[i][row]))
890
411k
        continue;
891
2.18k
      
isl_int_gcd2.18k
(a, left->row[row][row], left->row[i][row]);2.18k
892
2.18k
      isl_int_divexact(b, left->row[i][row], a);
893
2.18k
      isl_int_divexact(a, left->row[row][row], a);
894
2.18k
      isl_int_neg(b, b);
895
2.18k
      isl_seq_combine(left->row[i] + i,
896
2.18k
          a, left->row[i] + i,
897
2.18k
          b, left->row[row] + i,
898
2.18k
          left->n_col - i);
899
2.18k
      isl_seq_combine(right->row[i], a, right->row[i],
900
2.18k
          b, right->row[row], right->n_col);
901
2.18k
    }
902
276k
  }
903
89.8k
  
isl_int_clear89.8k
(b);89.8k
904
89.8k
905
89.8k
  isl_int_set(a, left->row[0][0]);
906
276k
  for (row = 1; 
row < left->n_row276k
;
++row186k
)
907
186k
    isl_int_lcm(a, a, left->row[row][row]);
908
89.8k
  if (
isl_int_is_zero89.8k
(a))
{0
909
0
    isl_int_clear(a);
910
0
    isl_assert(left->ctx, 0, goto error);
911
0
  }
912
365k
  
for (row = 0; 89.8k
row < left->n_row365k
;
++row276k
)
{276k
913
276k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
914
276k
    if (isl_int_is_one(left->row[row][row]))
915
270k
      continue;
916
5.92k
    isl_seq_scale(right->row[row], right->row[row],
917
5.92k
        left->row[row][row], right->n_col);
918
5.92k
  }
919
89.8k
  isl_int_clear(a);
920
89.8k
921
89.8k
  isl_mat_free(left);
922
89.8k
  return right;
923
0
error:
924
0
  isl_mat_free(left);
925
0
  isl_mat_free(right);
926
0
  return NULL;
927
89.8k
}
928
929
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
930
54.8k
{
931
54.8k
  int i;
932
54.8k
933
256k
  for (i = 0; 
i < mat->n_row256k
;
++i201k
)
934
201k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
935
54.8k
}
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
155k
{
940
155k
  int i;
941
155k
  isl_int tmp;
942
155k
943
155k
  isl_int_init(tmp);
944
695k
  for (i = 0; 
i < mat->n_row695k
;
++i539k
)
{539k
945
539k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
946
539k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
947
539k
    isl_int_set(mat->row[i][dst], tmp);
948
539k
  }
949
155k
  isl_int_clear(tmp);
950
155k
}
951
952
__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
953
52.0k
{
954
52.0k
  struct isl_mat *inv;
955
52.0k
  int row;
956
52.0k
  isl_int a, b;
957
52.0k
958
52.0k
  mat = isl_mat_cow(mat);
959
52.0k
  if (!mat)
960
0
    return NULL;
961
52.0k
962
52.0k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
963
52.0k
  inv = isl_mat_cow(inv);
964
52.0k
  if (!inv)
965
0
    goto error;
966
52.0k
967
52.0k
  
isl_int_init52.0k
(a);52.0k
968
52.0k
  isl_int_init(b);
969
201k
  for (row = 0; 
row < mat->n_row201k
;
++row149k
)
{149k
970
149k
    int pivot, first, i, off;
971
149k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
972
149k
    if (
pivot < 0149k
)
{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
149k
    pivot += row;
978
149k
    if (pivot != row)
979
21.7k
      exchange(mat, &inv, NULL, row, pivot, row);
980
149k
    if (isl_int_is_neg(mat->row[row][row]))
981
34.3k
      oppose(mat, &inv, NULL, row, row);
982
149k
    first = row+1;
983
405k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
984
256k
                mat->n_col-first)) != -1) {
985
256k
      first += off;
986
256k
      isl_int_fdiv_q(a, mat->row[row][first],
987
256k
                mat->row[row][row]);
988
256k
      subtract(mat, &inv, NULL, row, row, first, a);
989
256k
      if (
!256k
isl_int_is_zero256k
(mat->row[row][first]))
990
195k
        exchange(mat, &inv, NULL, row, row, first);
991
256k
      else
992
60.6k
        ++first;
993
256k
    }
994
316k
    for (i = 0; 
i < row316k
;
++i167k
)
{167k
995
167k
      if (isl_int_is_zero(mat->row[row][i]))
996
89.8k
        continue;
997
77.8k
      
isl_int_gcd77.8k
(a, mat->row[row][row], mat->row[row][i]);77.8k
998
77.8k
      isl_int_divexact(b, mat->row[row][i], a);
999
77.8k
      isl_int_divexact(a, mat->row[row][row], a);
1000
77.8k
      isl_int_neg(a, a);
1001
77.8k
      isl_mat_col_combine(mat, i, a, i, b, row);
1002
77.8k
      isl_mat_col_combine(inv, i, a, i, b, row);
1003
77.8k
    }
1004
149k
  }
1005
52.0k
  
isl_int_clear52.0k
(b);52.0k
1006
52.0k
1007
52.0k
  isl_int_set(a, mat->row[0][0]);
1008
149k
  for (row = 1; 
row < mat->n_row149k
;
++row96.9k
)
1009
96.9k
    isl_int_lcm(a, a, mat->row[row][row]);
1010
52.0k
  if (
isl_int_is_zero52.0k
(a))
{0
1011
0
    isl_int_clear(a);
1012
0
    goto error;
1013
0
  }
1014
201k
  
for (row = 0; 52.0k
row < mat->n_row201k
;
++row149k
)
{149k
1015
149k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1016
149k
    if (isl_int_is_one(mat->row[row][row]))
1017
94.2k
      continue;
1018
54.8k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1019
54.8k
  }
1020
52.0k
  isl_int_clear(a);
1021
52.0k
1022
52.0k
  isl_mat_free(mat);
1023
52.0k
1024
52.0k
  return inv;
1025
0
error:
1026
0
  isl_mat_free(mat);
1027
0
  isl_mat_free(inv);
1028
0
  return NULL;
1029
52.0k
}
1030
1031
__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1032
3.62k
{
1033
3.62k
  struct isl_mat *transpose = NULL;
1034
3.62k
  int i, j;
1035
3.62k
1036
3.62k
  if (!mat)
1037
0
    return NULL;
1038
3.62k
1039
3.62k
  
if (3.62k
mat->n_col == mat->n_row3.62k
)
{3.62k
1040
3.62k
    mat = isl_mat_cow(mat);
1041
3.62k
    if (!mat)
1042
0
      return NULL;
1043
14.8k
    
for (i = 0; 3.62k
i < mat->n_row14.8k
;
++i11.2k
)
1044
33.0k
      
for (j = i + 1; 11.2k
j < mat->n_col33.0k
;
++j21.8k
)
1045
21.8k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1046
3.62k
    return mat;
1047
3.62k
  }
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
352k
{
1064
352k
  int r;
1065
352k
1066
352k
  mat = isl_mat_cow(mat);
1067
352k
  if (!mat)
1068
0
    return NULL;
1069
352k
  
isl_assert352k
(mat->ctx, i < mat->n_col, goto error);352k
1070
352k
  
isl_assert352k
(mat->ctx, j < mat->n_col, goto error);352k
1071
352k
1072
8.59M
  
for (r = 0; 352k
r < mat->n_row8.59M
;
++r8.24M
)
1073
8.24M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1074
352k
  return mat;
1075
0
error:
1076
0
  isl_mat_free(mat);
1077
0
  return NULL;
1078
352k
}
1079
1080
__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1081
  unsigned i, unsigned j)
1082
1.06M
{
1083
1.06M
  isl_int *t;
1084
1.06M
1085
1.06M
  if (!mat)
1086
0
    return NULL;
1087
1.06M
  mat = isl_mat_cow(mat);
1088
1.06M
  if (!mat)
1089
0
    return NULL;
1090
1.06M
  t = mat->row[i];
1091
1.06M
  mat->row[i] = mat->row[j];
1092
1.06M
  mat->row[j] = t;
1093
1.06M
  return mat;
1094
1.06M
}
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.30M
{
1104
1.30M
  int i, j, k;
1105
1.30M
  struct isl_mat *prod;
1106
1.30M
1107
1.30M
  if (
!left || 1.30M
!right1.30M
)
1108
0
    goto error;
1109
1.30M
  
isl_assert1.30M
(left->ctx, left->n_col == right->n_row, goto error);1.30M
1110
1.30M
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1111
1.30M
  if (!prod)
1112
0
    goto error;
1113
1.30M
  
if (1.30M
left->n_col == 01.30M
)
{231
1114
363
    for (i = 0; 
i < prod->n_row363
;
++i132
)
1115
132
      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
3.79M
  
for (i = 0; 1.30M
i < prod->n_row3.79M
;
++i2.49M
)
{2.49M
1121
18.3M
    for (j = 0; 
j < prod->n_col18.3M
;
++j15.8M
)
1122
15.8M
      isl_int_mul(prod->row[i][j],
1123
2.49M
            left->row[i][0], right->row[0][j]);
1124
21.1M
    for (k = 1; 
k < left->n_col21.1M
;
++k18.7M
)
{18.7M
1125
18.7M
      if (isl_int_is_zero(left->row[i][k]))
1126
15.8M
        continue;
1127
23.8M
      
for (j = 0; 2.90M
j < prod->n_col23.8M
;
++j20.9M
)
1128
20.9M
        isl_int_addmul(prod->row[i][j],
1129
2.90M
              left->row[i][k], right->row[k][j]);
1130
2.90M
    }
1131
2.49M
  }
1132
1.30M
  isl_mat_free(left);
1133
1.30M
  isl_mat_free(right);
1134
1.30M
  return prod;
1135
0
error:
1136
0
  isl_mat_free(left);
1137
0
  isl_mat_free(right);
1138
0
  return NULL;
1139
1.30M
}
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.12M
{
1155
1.12M
  int i;
1156
1.12M
  struct isl_mat *t;
1157
1.12M
  int e;
1158
1.12M
1159
1.12M
  if (mat->n_col >= mat->n_row)
1160
731k
    e = 0;
1161
1.12M
  else
1162
393k
    e = mat->n_row - mat->n_col;
1163
1.12M
  if (has_div)
1164
375k
    
for (i = 0; 374k
i < n375k
;
++i38
)
1165
38
      isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
1166
1.12M
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1167
1.12M
  t = isl_mat_product(t, mat);
1168
1.12M
  if (!t)
1169
0
    return -1;
1170
2.76M
  
for (i = 0; 1.12M
i < n2.76M
;
++i1.63M
)
{1.63M
1171
1.63M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1172
1.63M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1173
1.63M
          q[i] + has_div + t->n_col + e, n_div);
1174
1.63M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1175
1.63M
  }
1176
1.12M
  isl_mat_free(t);
1177
1.12M
  return 0;
1178
1.12M
}
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
374k
{
1193
374k
  struct isl_ctx *ctx;
1194
374k
1195
374k
  if (
!bset || 374k
!mat374k
)
1196
0
    goto error;
1197
374k
1198
374k
  ctx = bset->ctx;
1199
374k
  bset = isl_basic_set_cow(bset);
1200
374k
  if (!bset)
1201
0
    goto error;
1202
374k
1203
374k
  
isl_assert374k
(ctx, bset->dim->nparam == 0, goto error);374k
1204
374k
  
isl_assert374k
(ctx, 1+bset->dim->n_out == mat->n_row, goto error);374k
1205
374k
  
isl_assert374k
(ctx, mat->n_col > 0, goto error);374k
1206
374k
1207
374k
  
if (374k
mat->n_col > mat->n_row374k
)
{62.9k
1208
62.9k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1209
62.9k
    if (!bset)
1210
0
      goto error;
1211
311k
  } else 
if (311k
mat->n_col < mat->n_row311k
)
{131k
1212
131k
    bset->dim = isl_space_cow(bset->dim);
1213
131k
    if (!bset->dim)
1214
0
      goto error;
1215
131k
    bset->dim->n_out -= mat->n_row - mat->n_col;
1216
131k
  }
1217
374k
1218
374k
  
if (374k
preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,374k
1219
374k
      isl_mat_copy(mat)) < 0)
1220
0
    goto error;
1221
374k
1222
374k
  
if (374k
preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,374k
1223
374k
      isl_mat_copy(mat)) < 0)
1224
0
    goto error;
1225
374k
1226
374k
  
if (374k
preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0374k
)
1227
0
    goto error2;
1228
374k
1229
374k
  
ISL_F_CLR374k
(bset, ISL_BASIC_SET_NO_IMPLICIT);374k
1230
374k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1231
374k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1232
374k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1233
374k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1234
374k
1235
374k
  bset = isl_basic_set_simplify(bset);
1236
374k
  bset = isl_basic_set_finalize(bset);
1237
374k
1238
374k
  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
28.2k
{
1248
28.2k
  int i;
1249
28.2k
1250
28.2k
  set = isl_set_cow(set);
1251
28.2k
  if (!set)
1252
0
    goto error;
1253
28.2k
1254
83.4k
  
for (i = 0; 28.2k
i < set->n83.4k
;
++i55.2k
)
{55.2k
1255
55.2k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1256
55.2k
                isl_mat_copy(mat));
1257
55.2k
    if (!set->p[i])
1258
0
      goto error;
1259
55.2k
  }
1260
28.2k
  
if (28.2k
mat->n_col != mat->n_row28.2k
)
{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
28.2k
  isl_mat_free(mat);
1268
28.2k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1269
28.2k
  return set;
1270
0
error:
1271
0
  isl_set_free(set);
1272
0
  isl_mat_free(mat);
1273
0
  return NULL;
1274
28.2k
}
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.09k
{
1283
4.09k
  int i;
1284
4.09k
  isl_ctx *ctx;
1285
4.09k
  isl_mat *t;
1286
4.09k
1287
4.09k
  if (!mat)
1288
0
    return isl_stat_error;
1289
4.09k
  ctx = isl_mat_get_ctx(mat);
1290
4.09k
  t = isl_mat_sub_alloc6(ctx, row, 0, n_row, first_col, mat->n_row);
1291
4.09k
  t = isl_mat_product(t, mat);
1292
4.09k
  if (!t)
1293
0
    return isl_stat_error;
1294
8.18k
  
for (i = 0; 4.09k
i < n_row8.18k
;
++i4.08k
)
1295
4.08k
    isl_seq_swp_or_cpy(row[i] + first_col, t->row[i], t->n_col);
1296
4.09k
  isl_mat_free(t);
1297
4.09k
  return isl_stat_ok;
1298
4.09k
}
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
57.3k
{
1336
57.3k
  int r;
1337
57.3k
1338
57.3k
  if (n == 0)
1339
2
    return mat;
1340
57.3k
1341
57.3k
  mat = isl_mat_cow(mat);
1342
57.3k
  if (!mat)
1343
0
    return NULL;
1344
57.3k
1345
57.3k
  
if (57.3k
col != mat->n_col-n57.3k
)
{14.0k
1346
64.1k
    for (r = 0; 
r < mat->n_row64.1k
;
++r50.0k
)
1347
50.0k
      isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1348
50.0k
          mat->n_col - col - n);
1349
14.0k
  }
1350
57.3k
  mat->n_col -= n;
1351
57.3k
  return mat;
1352
57.3k
}
1353
1354
struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n)
1355
70.5k
{
1356
70.5k
  int r;
1357
70.5k
1358
70.5k
  mat = isl_mat_cow(mat);
1359
70.5k
  if (!mat)
1360
0
    return NULL;
1361
70.5k
1362
170k
  
for (r = row; 70.5k
r+n < mat->n_row170k
;
++r99.6k
)
1363
99.6k
    mat->row[r] = mat->row[r+n];
1364
70.5k
1365
70.5k
  mat->n_row -= n;
1366
70.5k
  return mat;
1367
70.5k
}
1368
1369
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1370
        unsigned col, unsigned n)
1371
46.0k
{
1372
46.0k
  isl_mat *ext;
1373
46.0k
1374
46.0k
  if (!mat)
1375
0
    return NULL;
1376
46.0k
  
if (46.0k
n == 046.0k
)
1377
87
    return mat;
1378
46.0k
1379
45.9k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1380
45.9k
  if (!ext)
1381
0
    goto error;
1382
45.9k
1383
45.9k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1384
45.9k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1385
45.9k
        col + n, col, mat->n_col - col);
1386
45.9k
1387
45.9k
  isl_mat_free(mat);
1388
45.9k
  return ext;
1389
0
error:
1390
0
  isl_mat_free(mat);
1391
0
  return NULL;
1392
45.9k
}
1393
1394
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1395
  unsigned first, unsigned n)
1396
46.0k
{
1397
46.0k
  int i;
1398
46.0k
1399
46.0k
  if (!mat)
1400
0
    return NULL;
1401
46.0k
  mat = isl_mat_insert_cols(mat, first, n);
1402
46.0k
  if (!mat)
1403
0
    return NULL;
1404
46.0k
1405
49.0k
  
for (i = 0; 46.0k
i < mat->n_row49.0k
;
++i2.95k
)
1406
2.95k
    isl_seq_clr(mat->row[i] + first, n);
1407
46.0k
1408
46.0k
  return mat;
1409
46.0k
}
1410
1411
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1412
3.72k
{
1413
3.72k
  if (!mat)
1414
0
    return NULL;
1415
3.72k
1416
3.72k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1417
3.72k
}
1418
1419
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1420
        unsigned row, unsigned n)
1421
4.14k
{
1422
4.14k
  isl_mat *ext;
1423
4.14k
1424
4.14k
  if (!mat)
1425
0
    return NULL;
1426
4.14k
  
if (4.14k
n == 04.14k
)
1427
85
    return mat;
1428
4.14k
1429
4.06k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1430
4.06k
  if (!ext)
1431
0
    goto error;
1432
4.06k
1433
4.06k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1434
4.06k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1435
4.06k
        mat->n_row - row, 0, 0, mat->n_col);
1436
4.06k
1437
4.06k
  isl_mat_free(mat);
1438
4.06k
  return ext;
1439
0
error:
1440
0
  isl_mat_free(mat);
1441
0
  return NULL;
1442
4.06k
}
1443
1444
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1445
4.14k
{
1446
4.14k
  if (!mat)
1447
0
    return NULL;
1448
4.14k
1449
4.14k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1450
4.14k
}
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
2.42k
{
1478
2.42k
  int i;
1479
2.42k
1480
11.2k
  for (i = 0; 
i < mat->n_row11.2k
;
++i8.81k
)
1481
8.81k
    isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1482
2.42k
}
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
7.46k
{
1498
7.46k
  int i;
1499
7.46k
1500
26.9k
  for (i = 0; 
i < mat->n_row26.9k
;
++i19.4k
)
1501
19.4k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1502
7.46k
}
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
12.4k
{
1550
12.4k
  int r;
1551
12.4k
  struct isl_mat *H = NULL, *Q = NULL;
1552
12.4k
1553
12.4k
  if (!M)
1554
0
    return NULL;
1555
12.4k
1556
12.4k
  
isl_assert12.4k
(M->ctx, M->n_row == M->n_col, goto error);12.4k
1557
12.4k
  M->n_row = row;
1558
12.4k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1559
12.4k
  M->n_row = M->n_col;
1560
12.4k
  if (!H)
1561
0
    goto error;
1562
24.9k
  
for (r = 0; 12.4k
r < row24.9k
;
++r12.4k
)
1563
12.4k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1564
36.7k
  
for (r = row; 12.4k
r < M->n_row36.7k
;
++r24.2k
)
1565
24.2k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1566
12.4k
  isl_mat_free(H);
1567
12.4k
  isl_mat_free(Q);
1568
12.4k
  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
12.4k
}
1575
1576
__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1577
  __isl_take isl_mat *bot)
1578
371
{
1579
371
  struct isl_mat *mat;
1580
371
1581
371
  if (
!top || 371
!bot371
)
1582
0
    goto error;
1583
371
1584
371
  
isl_assert371
(top->ctx, top->n_col == bot->n_col, goto error);371
1585
371
  
if (371
top->n_row == 0371
)
{169
1586
169
    isl_mat_free(top);
1587
169
    return bot;
1588
169
  }
1589
202
  
if (202
bot->n_row == 0202
)
{0
1590
0
    isl_mat_free(bot);
1591
0
    return top;
1592
0
  }
1593
202
1594
202
  mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1595
202
  if (!mat)
1596
0
    goto error;
1597
202
  isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1598
202
       0, 0, mat->n_col);
1599
202
  isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1600
202
       0, 0, mat->n_col);
1601
202
  isl_mat_free(top);
1602
202
  isl_mat_free(bot);
1603
202
  return mat;
1604
0
error:
1605
0
  isl_mat_free(top);
1606
0
  isl_mat_free(bot);
1607
0
  return NULL;
1608
202
}
1609
1610
isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1611
5.54k
{
1612
5.54k
  int i;
1613
5.54k
1614
5.54k
  if (
!mat1 || 5.54k
!mat25.54k
)
1615
0
    return isl_bool_error;
1616
5.54k
1617
5.54k
  
if (5.54k
mat1->n_row != mat2->n_row5.54k
)
1618
789
    return isl_bool_false;
1619
5.54k
1620
4.75k
  
if (4.75k
mat1->n_col != mat2->n_col4.75k
)
1621
0
    return isl_bool_false;
1622
4.75k
1623
8.57k
  
for (i = 0; 4.75k
i < mat1->n_row8.57k
;
++i3.82k
)
1624
4.05k
    
if (4.05k
!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col)4.05k
)
1625
238
      return isl_bool_false;
1626
4.75k
1627
4.51k
  return isl_bool_true;
1628
4.75k
}
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
423
{
1678
423
  isl_mat *res;
1679
423
1680
423
  if (!mat)
1681
0
    return NULL;
1682
423
  
if (423
n == 0 || 423
dst_col == src_col423
)
1683
415
    return mat;
1684
423
1685
8
  res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1686
8
  if (!res)
1687
0
    goto error;
1688
8
1689
8
  
if (8
dst_col < src_col8
)
{8
1690
8
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1691
8
         0, 0, dst_col);
1692
8
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1693
8
         dst_col, src_col, n);
1694
8
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1695
8
         dst_col + n, dst_col, src_col - dst_col);
1696
8
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1697
8
         src_col + n, src_col + n,
1698
8
         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
8
  isl_mat_free(mat);
1711
8
1712
8
  return res;
1713
0
error:
1714
0
  isl_mat_free(mat);
1715
0
  return NULL;
1716
8
}
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
713
{
1736
713
  int i;
1737
713
  isl_int g;
1738
713
1739
713
  isl_int_set_si(*gcd, 0);
1740
713
  if (!mat)
1741
0
    return;
1742
713
1743
713
  
isl_int_init713
(g);713
1744
4.60k
  for (i = 0; 
i < mat->n_row4.60k
;
++i3.89k
)
{3.89k
1745
3.89k
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1746
3.89k
    isl_int_gcd(*gcd, *gcd, g);
1747
3.89k
  }
1748
713
  isl_int_clear(g);
1749
713
}
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
713
{
1772
713
  int i;
1773
713
1774
713
  if (isl_int_is_one(m))
1775
0
    return mat;
1776
713
1777
713
  mat = isl_mat_cow(mat);
1778
713
  if (!mat)
1779
0
    return NULL;
1780
713
1781
4.60k
  
for (i = 0; 713
i < mat->n_row4.60k
;
++i3.89k
)
1782
3.89k
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1783
713
1784
713
  return mat;
1785
713
}
1786
1787
__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row,
1788
  isl_int m)
1789
14
{
1790
14
  if (isl_int_is_one(m))
1791
0
    return mat;
1792
14
1793
14
  mat = isl_mat_cow(mat);
1794
14
  if (!mat)
1795
0
    return NULL;
1796
14
1797
14
  isl_seq_scale_down(mat->row[row], mat->row[row], m, mat->n_col);
1798
14
1799
14
  return mat;
1800
14
}
1801
1802
__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1803
713
{
1804
713
  isl_int gcd;
1805
713
1806
713
  if (!mat)
1807
0
    return NULL;
1808
713
1809
713
  
isl_int_init713
(gcd);713
1810
713
  isl_mat_gcd(mat, &gcd);
1811
713
  mat = isl_mat_scale_down(mat, gcd);
1812
713
  isl_int_clear(gcd);
1813
713
1814
713
  return mat;
1815
713
}
1816
1817
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1818
1.10k
{
1819
1.10k
  mat = isl_mat_cow(mat);
1820
1.10k
  if (!mat)
1821
0
    return NULL;
1822
1.10k
1823
1.10k
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
1824
1.10k
1825
1.10k
  return mat;
1826
1.10k
}
1827
1828
/* Number of initial non-zero columns.
1829
 */
1830
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
1831
952
{
1832
952
  int i;
1833
952
1834
952
  if (!mat)
1835
0
    return -1;
1836
952
1837
1.37k
  
for (i = 0; 952
i < mat->n_col1.37k
;
++i423
)
1838
1.18k
    
if (1.18k
row_first_non_zero(mat->row, mat->n_row, i) < 01.18k
)
1839
766
      break;
1840
952
1841
952
  return i;
1842
952
}