Coverage Report

Created: 2017-03-28 09:59

/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
1.93M
{
24
1.93M
  return mat ? mat->ctx : NULL;
25
1.93M
}
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.33M
{
53
3.33M
  int i;
54
3.33M
  struct isl_mat *mat;
55
3.33M
56
3.33M
  mat = isl_alloc_type(ctx, struct isl_mat);
57
3.33M
  if (!mat)
58
0
    return NULL;
59
3.33M
60
3.33M
  mat->row = NULL;
61
3.33M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
62
3.33M
  if (isl_blk_is_error(mat->block))
63
0
    goto error;
64
3.33M
  
mat->row = 3.33M
isl_alloc_array3.33M
(ctx, isl_int *, n_row);
65
3.33M
  if (
n_row && 3.33M
!mat->row2.40M
)
66
0
    goto error;
67
3.33M
68
16.9M
  
for (i = 0; 3.33M
i < n_row16.9M
;
++i13.5M
)
69
13.5M
    mat->row[i] = mat->block.data + i * n_col;
70
3.33M
71
3.33M
  mat->ctx = ctx;
72
3.33M
  isl_ctx_ref(ctx);
73
3.33M
  mat->ref = 1;
74
3.33M
  mat->n_row = n_row;
75
3.33M
  mat->n_col = n_col;
76
3.33M
  mat->max_col = n_col;
77
3.33M
  mat->flags = 0;
78
3.33M
79
3.33M
  return mat;
80
0
error:
81
0
  isl_blk_free(ctx, mat->block);
82
0
  free(mat);
83
0
  return NULL;
84
3.33M
}
85
86
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
87
  unsigned n_row, unsigned n_col)
88
63.8k
{
89
63.8k
  int i;
90
63.8k
  isl_int *old;
91
63.8k
  isl_int **row;
92
63.8k
93
63.8k
  if (!mat)
94
0
    return NULL;
95
63.8k
96
63.8k
  
if (63.8k
mat->max_col >= n_col && 63.8k
mat->n_row >= n_row61.1k
)
{3.66k
97
3.66k
    if (mat->n_col < n_col)
98
189
      mat->n_col = n_col;
99
3.66k
    return mat;
100
3.66k
  }
101
63.8k
102
60.1k
  
if (60.1k
mat->max_col < n_col60.1k
)
{2.69k
103
2.69k
    struct isl_mat *new_mat;
104
2.69k
105
2.69k
    if (n_row < mat->n_row)
106
2
      n_row = mat->n_row;
107
2.69k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
108
2.69k
    if (!new_mat)
109
0
      goto error;
110
21.0k
    
for (i = 0; 2.69k
i < mat->n_row21.0k
;
++i18.3k
)
111
18.3k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
112
2.69k
    isl_mat_free(mat);
113
2.69k
    return new_mat;
114
2.69k
  }
115
60.1k
116
57.4k
  mat = isl_mat_cow(mat);
117
57.4k
  if (!mat)
118
0
    goto error;
119
57.4k
120
57.4k
  old = mat->block.data;
121
57.4k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
122
57.4k
  if (isl_blk_is_error(mat->block))
123
0
    goto error;
124
57.4k
  
row = 57.4k
isl_realloc_array57.4k
(mat->ctx, mat->row, isl_int *, n_row);
125
57.4k
  if (
n_row && 57.4k
!row57.4k
)
126
0
    goto error;
127
57.4k
  mat->row = row;
128
57.4k
129
503k
  for (i = 0; 
i < mat->n_row503k
;
++i446k
)
130
446k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
131
174k
  for (i = mat->n_row; 
i < n_row174k
;
++i117k
)
132
117k
    mat->row[i] = mat->block.data + i * mat->max_col;
133
57.4k
  mat->n_row = n_row;
134
57.4k
  if (mat->n_col < n_col)
135
0
    mat->n_col = n_col;
136
57.4k
137
57.4k
  return mat;
138
0
error:
139
0
  isl_mat_free(mat);
140
0
  return NULL;
141
57.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.59M
{
146
1.59M
  int i;
147
1.59M
  struct isl_mat *mat;
148
1.59M
149
1.59M
  mat = isl_alloc_type(ctx, struct isl_mat);
150
1.59M
  if (!mat)
151
0
    return NULL;
152
1.59M
  
mat->row = 1.59M
isl_alloc_array1.59M
(ctx, isl_int *, n_row);
153
1.59M
  if (
n_row && 1.59M
!mat->row832k
)
154
0
    goto error;
155
4.82M
  
for (i = 0; 1.59M
i < n_row4.82M
;
++i3.23M
)
156
3.23M
    mat->row[i] = row[first_row+i] + first_col;
157
1.59M
  mat->ctx = ctx;
158
1.59M
  isl_ctx_ref(ctx);
159
1.59M
  mat->ref = 1;
160
1.59M
  mat->n_row = n_row;
161
1.59M
  mat->n_col = n_col;
162
1.59M
  mat->block = isl_blk_empty();
163
1.59M
  mat->flags = ISL_MAT_BORROWED;
164
1.59M
  return mat;
165
0
error:
166
0
  free(mat);
167
0
  return NULL;
168
1.59M
}
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
292k
{
173
292k
  if (!mat)
174
0
    return NULL;
175
292k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
176
292k
          first_col, n_col);
177
292k
}
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
195k
{
182
195k
  int i;
183
195k
184
650k
  for (i = 0; 
i < n_row650k
;
++i454k
)
185
454k
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
186
195k
}
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
66.3k
{
191
66.3k
  int i;
192
66.3k
193
198k
  for (i = 0; 
i < n_row198k
;
++i132k
)
194
132k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
195
66.3k
}
196
197
struct isl_mat *isl_mat_copy(struct isl_mat *mat)
198
1.27M
{
199
1.27M
  if (!mat)
200
0
    return NULL;
201
1.27M
202
1.27M
  mat->ref++;
203
1.27M
  return mat;
204
1.27M
}
205
206
struct isl_mat *isl_mat_dup(struct isl_mat *mat)
207
207k
{
208
207k
  int i;
209
207k
  struct isl_mat *mat2;
210
207k
211
207k
  if (!mat)
212
1.15k
    return NULL;
213
206k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
214
206k
  if (!mat2)
215
0
    return NULL;
216
623k
  
for (i = 0; 206k
i < mat->n_row623k
;
++i416k
)
217
416k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
218
206k
  return mat2;
219
206k
}
220
221
struct isl_mat *isl_mat_cow(struct isl_mat *mat)
222
1.99M
{
223
1.99M
  struct isl_mat *mat2;
224
1.99M
  if (!mat)
225
0
    return NULL;
226
1.99M
227
1.99M
  
if (1.99M
mat->ref == 1 && 1.99M
!1.97M
ISL_F_ISSET1.97M
(mat, ISL_MAT_BORROWED))
228
1.79M
    return mat;
229
1.99M
230
205k
  mat2 = isl_mat_dup(mat);
231
205k
  isl_mat_free(mat);
232
205k
  return mat2;
233
1.99M
}
234
235
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
236
7.09M
{
237
7.09M
  if (!mat)
238
882k
    return NULL;
239
7.09M
240
6.20M
  
if (6.20M
--mat->ref > 06.20M
)
241
1.27M
    return NULL;
242
6.20M
243
4.92M
  
if (4.92M
!4.92M
ISL_F_ISSET4.92M
(mat, ISL_MAT_BORROWED))
244
3.33M
    isl_blk_free(mat->ctx, mat->block);
245
4.92M
  isl_ctx_deref(mat->ctx);
246
4.92M
  free(mat->row);
247
4.92M
  free(mat);
248
4.92M
249
4.92M
  return NULL;
250
6.20M
}
251
252
int isl_mat_rows(__isl_keep isl_mat *mat)
253
78.7k
{
254
78.7k
  return mat ? 
mat->n_row78.7k
:
-10
;
255
78.7k
}
256
257
int isl_mat_cols(__isl_keep isl_mat *mat)
258
22.6k
{
259
22.6k
  return mat ? 
mat->n_col22.6k
:
-10
;
260
22.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.30k
{
266
4.30k
  if (!mat)
267
0
    return isl_stat_error;
268
4.30k
  
if (4.30k
col < 0 || 4.30k
col >= mat->n_col4.30k
)
269
0
    isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
270
4.30k
      "column out of range", return isl_stat_error);
271
4.30k
  return isl_stat_ok;
272
4.30k
}
273
274
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
275
2.86k
{
276
2.86k
  if (!mat)
277
0
    return -1;
278
2.86k
  
if (2.86k
row < 0 || 2.86k
row >= mat->n_row2.86k
)
279
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
280
2.86k
      return -1);
281
2.86k
  
if (2.86k
check_col(mat, col) < 02.86k
)
282
0
    return -1;
283
2.86k
  
isl_int_set2.86k
(*v, mat->row[row][col]);2.86k
284
2.86k
  return 0;
285
2.86k
}
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
12
{
292
12
  isl_ctx *ctx;
293
12
294
12
  if (!mat)
295
0
    return NULL;
296
12
  ctx = isl_mat_get_ctx(mat);
297
12
  if (
row < 0 || 12
row >= mat->n_row12
)
298
0
    isl_die(ctx, isl_error_invalid, "row out of range",
299
12
      return NULL);
300
12
  
if (12
check_col(mat, col) < 012
)
301
0
    return NULL;
302
12
  return isl_val_int_from_isl_int(ctx, mat->row[row][col]);
303
12
}
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.36k
{
308
1.36k
  mat = isl_mat_cow(mat);
309
1.36k
  if (!mat)
310
0
    return NULL;
311
1.36k
  
if (1.36k
row < 0 || 1.36k
row >= mat->n_row1.36k
)
312
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
313
1.36k
      goto error);
314
1.36k
  
if (1.36k
check_col(mat, col) < 01.36k
)
315
0
    return isl_mat_free(mat);
316
1.36k
  
isl_int_set1.36k
(mat->row[row][col], v);1.36k
317
1.36k
  return mat;
318
0
error:
319
0
  isl_mat_free(mat);
320
0
  return NULL;
321
1.36k
}
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
522k
{
361
522k
  int i;
362
522k
  struct isl_mat *mat;
363
522k
364
522k
  mat = isl_mat_alloc(ctx, n_row, n_row);
365
522k
  if (!mat)
366
0
    return NULL;
367
2.68M
  
for (i = 0; 522k
i < n_row2.68M
;
++i2.16M
)
{2.16M
368
2.16M
    isl_seq_clr(mat->row[i], i);
369
2.16M
    isl_int_set(mat->row[i][i], d);
370
2.16M
    isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
371
2.16M
  }
372
522k
373
522k
  return mat;
374
522k
}
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
522k
{
394
522k
  if (!ctx)
395
0
    return NULL;
396
522k
  return isl_mat_diag(ctx, n_row, ctx->one);
397
522k
}
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
struct isl_vec *isl_mat_vec_product(struct isl_mat *mat, struct isl_vec *vec)
424
184k
{
425
184k
  int i;
426
184k
  struct isl_vec *prod;
427
184k
428
184k
  if (
!mat || 184k
!vec184k
)
429
0
    goto error;
430
184k
431
184k
  
isl_assert184k
(mat->ctx, mat->n_col == vec->size, goto error);184k
432
184k
433
184k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
434
184k
  if (!prod)
435
0
    goto error;
436
184k
437
1.33M
  
for (i = 0; 184k
i < prod->size1.33M
;
++i1.15M
)
438
1.15M
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
439
1.15M
          &prod->block.data[i]);
440
184k
  isl_mat_free(mat);
441
184k
  isl_vec_free(vec);
442
184k
  return prod;
443
0
error:
444
0
  isl_mat_free(mat);
445
0
  isl_vec_free(vec);
446
0
  return NULL;
447
184k
}
448
449
__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
450
  __isl_take isl_vec *vec)
451
243
{
452
243
  struct isl_mat *vec_mat;
453
243
  int i;
454
243
455
243
  if (
!mat || 243
!vec243
)
456
0
    goto error;
457
243
  vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
458
243
  if (!vec_mat)
459
0
    goto error;
460
1.97k
  
for (i = 0; 243
i < vec->size1.97k
;
++i1.73k
)
461
1.73k
    isl_int_set(vec_mat->row[i][0], vec->el[i]);
462
243
  vec_mat = isl_mat_inverse_product(mat, vec_mat);
463
243
  isl_vec_free(vec);
464
243
  if (!vec_mat)
465
0
    return NULL;
466
243
  vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
467
243
  if (vec)
468
1.97k
    
for (i = 0; 243
i < vec->size1.97k
;
++i1.73k
)
469
1.73k
      isl_int_set(vec->el[i], vec_mat->row[i][0]);
470
243
  isl_mat_free(vec_mat);
471
243
  return vec;
472
0
error:
473
0
  isl_mat_free(mat);
474
0
  isl_vec_free(vec);
475
0
  return NULL;
476
243
}
477
478
struct isl_vec *isl_vec_mat_product(struct isl_vec *vec, struct isl_mat *mat)
479
11.4k
{
480
11.4k
  int i, j;
481
11.4k
  struct isl_vec *prod;
482
11.4k
483
11.4k
  if (
!mat || 11.4k
!vec11.4k
)
484
0
    goto error;
485
11.4k
486
11.4k
  
isl_assert11.4k
(mat->ctx, mat->n_row == vec->size, goto error);11.4k
487
11.4k
488
11.4k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
489
11.4k
  if (!prod)
490
0
    goto error;
491
11.4k
492
58.7k
  
for (i = 0; 11.4k
i < prod->size58.7k
;
++i47.3k
)
{47.3k
493
47.3k
    isl_int_set_si(prod->el[i], 0);
494
361k
    for (j = 0; 
j < vec->size361k
;
++j313k
)
495
313k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
496
47.3k
  }
497
11.4k
  isl_mat_free(mat);
498
11.4k
  isl_vec_free(vec);
499
11.4k
  return prod;
500
0
error:
501
0
  isl_mat_free(mat);
502
0
  isl_vec_free(vec);
503
0
  return NULL;
504
11.4k
}
505
506
struct isl_mat *isl_mat_aff_direct_sum(struct isl_mat *left,
507
  struct isl_mat *right)
508
66.3k
{
509
66.3k
  int i;
510
66.3k
  struct isl_mat *sum;
511
66.3k
512
66.3k
  if (
!left || 66.3k
!right66.3k
)
513
0
    goto error;
514
66.3k
515
66.3k
  
isl_assert66.3k
(left->ctx, left->n_row == right->n_row, goto error);66.3k
516
66.3k
  
isl_assert66.3k
(left->ctx, left->n_row >= 1, goto error);66.3k
517
66.3k
  
isl_assert66.3k
(left->ctx, left->n_col >= 1, goto error);66.3k
518
66.3k
  
isl_assert66.3k
(left->ctx, right->n_col >= 1, goto error);66.3k
519
66.3k
  
isl_assert66.3k
(left->ctx,66.3k
520
66.3k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
521
66.3k
      goto error);
522
66.3k
  
isl_assert66.3k
(left->ctx,66.3k
523
66.3k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
524
66.3k
      goto error);
525
66.3k
526
66.3k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
527
66.3k
  if (!sum)
528
0
    goto error;
529
66.3k
  
isl_int_lcm66.3k
(sum->row[0][0], left->row[0][0], right->row[0][0]);66.3k
530
66.3k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
531
66.3k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
532
66.3k
533
66.3k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
534
462k
  for (i = 1; 
i < sum->n_row462k
;
++i396k
)
{396k
535
396k
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
536
396k
    isl_int_addmul(sum->row[i][0],
537
396k
        right->row[0][0], right->row[i][0]);
538
396k
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
539
396k
        left->n_col-1);
540
396k
    isl_seq_scale(sum->row[i]+left->n_col,
541
396k
        right->row[i]+1, right->row[0][0],
542
396k
        right->n_col-1);
543
396k
  }
544
66.3k
545
66.3k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
546
66.3k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
547
66.3k
  isl_mat_free(left);
548
66.3k
  isl_mat_free(right);
549
66.3k
  return sum;
550
0
error:
551
0
  isl_mat_free(left);
552
0
  isl_mat_free(right);
553
0
  return NULL;
554
66.3k
}
555
556
static void exchange(struct isl_mat *M, struct isl_mat **U,
557
  struct isl_mat **Q, unsigned row, unsigned i, unsigned j)
558
513k
{
559
513k
  int r;
560
3.07M
  for (r = row; 
r < M->n_row3.07M
;
++r2.56M
)
561
2.56M
    isl_int_swap(M->row[r][i], M->row[r][j]);
562
513k
  if (
U513k
)
{497k
563
4.10M
    for (r = 0; 
r < (*U)->n_row4.10M
;
++r3.60M
)
564
3.60M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
565
497k
  }
566
513k
  if (Q)
567
135k
    isl_mat_swap_rows(*Q, i, j);
568
513k
}
569
570
static void subtract(struct isl_mat *M, struct isl_mat **U,
571
  struct isl_mat **Q, unsigned row, unsigned i, unsigned j, isl_int m)
572
508k
{
573
508k
  int r;
574
1.83M
  for (r = row; 
r < M->n_row1.83M
;
++r1.32M
)
575
1.32M
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
576
508k
  if (
U508k
)
{489k
577
3.03M
    for (r = 0; 
r < (*U)->n_row3.03M
;
++r2.54M
)
578
2.54M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
579
489k
  }
580
508k
  if (
Q508k
)
{111k
581
690k
    for (r = 0; 
r < (*Q)->n_col690k
;
++r579k
)
582
579k
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
583
111k
  }
584
508k
}
585
586
static void oppose(struct isl_mat *M, struct isl_mat **U,
587
  struct isl_mat **Q, unsigned row, unsigned col)
588
152k
{
589
152k
  int r;
590
552k
  for (r = row; 
r < M->n_row552k
;
++r400k
)
591
400k
    isl_int_neg(M->row[r][col], M->row[r][col]);
592
152k
  if (
U152k
)
{133k
593
829k
    for (r = 0; 
r < (*U)->n_row829k
;
++r695k
)
594
695k
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
595
133k
  }
596
152k
  if (Q)
597
69.0k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
598
152k
}
599
600
/* Given matrix M, compute
601
 *
602
 *    M U = H
603
 *    M   = H Q
604
 *
605
 * with U and Q unimodular matrices and H a matrix in column echelon form
606
 * such that on each echelon row the entries in the non-echelon column
607
 * are non-negative (if neg == 0) or non-positive (if neg == 1)
608
 * and strictly smaller (in absolute value) than the entries in the echelon
609
 * column.
610
 * If U or Q are NULL, then these matrices are not computed.
611
 */
612
struct isl_mat *isl_mat_left_hermite(struct isl_mat *M, int neg,
613
  struct isl_mat **U, struct isl_mat **Q)
614
252k
{
615
252k
  isl_int c;
616
252k
  int row, col;
617
252k
618
252k
  if (U)
619
232k
    *U = NULL;
620
252k
  if (Q)
621
115k
    *Q = NULL;
622
252k
  if (!M)
623
0
    goto error;
624
252k
  M = isl_mat_cow(M);
625
252k
  if (!M)
626
0
    goto error;
627
252k
  
if (252k
U252k
)
{232k
628
232k
    *U = isl_mat_identity(M->ctx, M->n_col);
629
232k
    if (!*U)
630
0
      goto error;
631
232k
  }
632
252k
  
if (252k
Q252k
)
{115k
633
115k
    *Q = isl_mat_identity(M->ctx, M->n_col);
634
115k
    if (!*Q)
635
0
      goto error;
636
115k
  }
637
252k
638
252k
  col = 0;
639
252k
  isl_int_init(c);
640
1.05M
  for (row = 0; 
row < M->n_row1.05M
;
++row802k
)
{802k
641
802k
    int first, i, off;
642
802k
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
643
802k
    if (first == -1)
644
238k
      continue;
645
564k
    first += col;
646
564k
    if (first != col)
647
283k
      exchange(M, U, Q, row, first, col);
648
564k
    if (isl_int_is_neg(M->row[row][col]))
649
118k
      oppose(M, U, Q, row, col);
650
564k
    first = col+1;
651
758k
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
652
194k
                   M->n_col-first)) != -1) {
653
194k
      first += off;
654
194k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
655
194k
      subtract(M, U, Q, row, col, first, c);
656
194k
      if (
!194k
isl_int_is_zero194k
(M->row[row][first]))
657
12.1k
        exchange(M, U, Q, row, first, col);
658
194k
      else
659
182k
        ++first;
660
194k
    }
661
2.08M
    for (i = 0; 
i < col2.08M
;
++i1.52M
)
{1.52M
662
1.52M
      if (isl_int_is_zero(M->row[row][i]))
663
1.44M
        continue;
664
78.7k
      
if (78.7k
neg78.7k
)
665
0
        isl_int_cdiv_q(c, M->row[row][i], M->row[row][col]);
666
78.7k
      else
667
78.7k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
668
78.7k
      if (isl_int_is_zero(c))
669
21.1k
        continue;
670
57.5k
      subtract(M, U, Q, row, col, i, c);
671
57.5k
    }
672
564k
    ++col;
673
564k
  }
674
252k
  isl_int_clear(c);
675
252k
676
252k
  return M;
677
0
error:
678
0
  if (
Q0
)
{0
679
0
    isl_mat_free(*Q);
680
0
    *Q = NULL;
681
0
  }
682
0
  if (
U0
)
{0
683
0
    isl_mat_free(*U);
684
0
    *U = NULL;
685
0
  }
686
0
  isl_mat_free(M);
687
0
  return NULL;
688
252k
}
689
690
struct isl_mat *isl_mat_right_kernel(struct isl_mat *mat)
691
0
{
692
0
  int i, rank;
693
0
  struct isl_mat *U = NULL;
694
0
  struct isl_mat *K;
695
0
696
0
  mat = isl_mat_left_hermite(mat, 0, &U, NULL);
697
0
  if (
!mat || 0
!U0
)
698
0
    goto error;
699
0
700
0
  
for (i = 0, rank = 0; 0
rank < mat->n_col0
;
++rank0
)
{0
701
0
    while (
i < mat->n_row && 0
isl_int_is_zero0
(mat->row[i][rank]))
702
0
      ++i;
703
0
    if (i >= mat->n_row)
704
0
      break;
705
0
  }
706
0
  K = isl_mat_alloc(U->ctx, U->n_row, U->n_col - rank);
707
0
  if (!K)
708
0
    goto error;
709
0
  isl_mat_sub_copy(K->ctx, K->row, U->row, U->n_row, 0, rank, U->n_col-rank);
710
0
  isl_mat_free(mat);
711
0
  isl_mat_free(U);
712
0
  return K;
713
0
error:
714
0
  isl_mat_free(mat);
715
0
  isl_mat_free(U);
716
0
  return NULL;
717
0
}
718
719
struct isl_mat *isl_mat_lin_to_aff(struct isl_mat *mat)
720
355k
{
721
355k
  int i;
722
355k
  struct isl_mat *mat2;
723
355k
724
355k
  if (!mat)
725
0
    return NULL;
726
355k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
727
355k
  if (!mat2)
728
0
    goto error;
729
355k
  
isl_int_set_si355k
(mat2->row[0][0], 1);355k
730
355k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
731
2.03M
  for (i = 0; 
i < mat->n_row2.03M
;
++i1.67M
)
{1.67M
732
1.67M
    isl_int_set_si(mat2->row[1+i][0], 0);
733
1.67M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
734
1.67M
  }
735
355k
  isl_mat_free(mat);
736
355k
  return mat2;
737
0
error:
738
0
  isl_mat_free(mat);
739
0
  return NULL;
740
355k
}
741
742
/* Given two matrices M1 and M2, return the block matrix
743
 *
744
 *  [ M1  0  ]
745
 *  [ 0   M2 ]
746
 */
747
__isl_give isl_mat *isl_mat_diagonal(__isl_take isl_mat *mat1,
748
  __isl_take isl_mat *mat2)
749
59.5k
{
750
59.5k
  int i;
751
59.5k
  isl_mat *mat;
752
59.5k
753
59.5k
  if (
!mat1 || 59.5k
!mat259.5k
)
754
0
    goto error;
755
59.5k
756
59.5k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
757
59.5k
               mat1->n_col + mat2->n_col);
758
59.5k
  if (!mat)
759
0
    goto error;
760
119k
  
for (i = 0; 59.5k
i < mat1->n_row119k
;
++i59.6k
)
{59.6k
761
59.6k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
762
59.6k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
763
59.6k
  }
764
351k
  for (i = 0; 
i < mat2->n_row351k
;
++i291k
)
{291k
765
291k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
766
291k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
767
291k
                mat2->row[i], mat2->n_col);
768
291k
  }
769
59.5k
  isl_mat_free(mat1);
770
59.5k
  isl_mat_free(mat2);
771
59.5k
  return mat;
772
0
error:
773
0
  isl_mat_free(mat1);
774
0
  isl_mat_free(mat2);
775
0
  return NULL;
776
59.5k
}
777
778
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
779
445k
{
780
445k
  int i;
781
445k
782
742k
  for (i = 0; 
i < n_row742k
;
++i296k
)
783
521k
    
if (521k
!521k
isl_int_is_zero521k
(row[i][col]))
784
224k
      return i;
785
220k
  return -1;
786
445k
}
787
788
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
789
219k
{
790
219k
  int i, min = row_first_non_zero(row, n_row, col);
791
219k
  if (min < 0)
792
0
    return -1;
793
519k
  
for (i = min + 1; 219k
i < n_row519k
;
++i299k
)
{299k
794
299k
    if (isl_int_is_zero(row[i][col]))
795
295k
      continue;
796
4.39k
    
if (4.39k
isl_int_abs_lt4.39k
(row[i][col], row[min][col]))
797
662
      min = i;
798
4.39k
  }
799
219k
  return min;
800
219k
}
801
802
static isl_stat inv_exchange(__isl_keep isl_mat **left,
803
  __isl_keep isl_mat **right, unsigned i, unsigned j)
804
1.05k
{
805
1.05k
  *left = isl_mat_swap_rows(*left, i, j);
806
1.05k
  *right = isl_mat_swap_rows(*right, i, j);
807
1.05k
808
1.05k
  if (
!*left || 1.05k
!*right1.05k
)
809
0
    return isl_stat_error;
810
1.05k
  return isl_stat_ok;
811
1.05k
}
812
813
static void inv_oppose(
814
  struct isl_mat *left, struct isl_mat *right, unsigned row)
815
597
{
816
597
  isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
817
597
  isl_seq_neg(right->row[row], right->row[row], right->n_col);
818
597
}
819
820
static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
821
  unsigned row, unsigned i, isl_int m)
822
4.68k
{
823
4.68k
  isl_int_neg(m, m);
824
4.68k
  isl_seq_combine(left->row[i]+row,
825
4.68k
      left->ctx->one, left->row[i]+row,
826
4.68k
      m, left->row[row]+row,
827
4.68k
      left->n_col-row);
828
4.68k
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
829
4.68k
      m, right->row[row], right->n_col);
830
4.68k
}
831
832
/* Compute inv(left)*right
833
 */
834
struct isl_mat *isl_mat_inverse_product(struct isl_mat *left,
835
  struct isl_mat *right)
836
75.0k
{
837
75.0k
  int row;
838
75.0k
  isl_int a, b;
839
75.0k
840
75.0k
  if (
!left || 75.0k
!right75.0k
)
841
0
    goto error;
842
75.0k
843
75.0k
  
isl_assert75.0k
(left->ctx, left->n_row == left->n_col, goto error);75.0k
844
75.0k
  
isl_assert75.0k
(left->ctx, left->n_row == right->n_row, goto error);75.0k
845
75.0k
846
75.0k
  
if (75.0k
left->n_row == 075.0k
)
{0
847
0
    isl_mat_free(left);
848
0
    return right;
849
0
  }
850
75.0k
851
75.0k
  left = isl_mat_cow(left);
852
75.0k
  right = isl_mat_cow(right);
853
75.0k
  if (
!left || 75.0k
!right75.0k
)
854
0
    goto error;
855
75.0k
856
75.0k
  
isl_int_init75.0k
(a);75.0k
857
75.0k
  isl_int_init(b);
858
294k
  for (row = 0; 
row < left->n_row294k
;
++row219k
)
{219k
859
219k
    int pivot, first, i, off;
860
219k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
861
219k
    if (
pivot < 0219k
)
{0
862
0
      isl_int_clear(a);
863
0
      isl_int_clear(b);
864
0
      isl_assert(left->ctx, pivot >= 0, goto error);
865
0
    }
866
219k
    pivot += row;
867
219k
    if (pivot != row)
868
769
      
if (769
inv_exchange(&left, &right, pivot, row) < 0769
)
869
0
        goto error;
870
219k
    
if (219k
isl_int_is_neg219k
(left->row[row][row]))
871
597
      inv_oppose(left, right, row);
872
219k
    first = row+1;
873
224k
    while ((off = row_first_non_zero(left->row+first,
874
4.68k
          left->n_row-first, row)) != -1) {
875
4.68k
      first += off;
876
4.68k
      isl_int_fdiv_q(a, left->row[first][row],
877
4.68k
          left->row[row][row]);
878
4.68k
      inv_subtract(left, right, row, first, a);
879
4.68k
      if (
!4.68k
isl_int_is_zero4.68k
(left->row[first][row]))
{287
880
287
        if (inv_exchange(&left, &right, row, first) < 0)
881
0
          goto error;
882
4.39k
      } else {
883
4.39k
        ++first;
884
4.39k
      }
885
4.68k
    }
886
520k
    
for (i = 0; 219k
i < row520k
;
++i300k
)
{300k
887
300k
      if (isl_int_is_zero(left->row[i][row]))
888
298k
        continue;
889
2.17k
      
isl_int_gcd2.17k
(a, left->row[row][row], left->row[i][row]);2.17k
890
2.17k
      isl_int_divexact(b, left->row[i][row], a);
891
2.17k
      isl_int_divexact(a, left->row[row][row], a);
892
2.17k
      isl_int_neg(b, b);
893
2.17k
      isl_seq_combine(left->row[i] + i,
894
2.17k
          a, left->row[i] + i,
895
2.17k
          b, left->row[row] + i,
896
2.17k
          left->n_col - i);
897
2.17k
      isl_seq_combine(right->row[i], a, right->row[i],
898
2.17k
          b, right->row[row], right->n_col);
899
2.17k
    }
900
219k
  }
901
75.0k
  
isl_int_clear75.0k
(b);75.0k
902
75.0k
903
75.0k
  isl_int_set(a, left->row[0][0]);
904
219k
  for (row = 1; 
row < left->n_row219k
;
++row144k
)
905
144k
    isl_int_lcm(a, a, left->row[row][row]);
906
75.0k
  if (
isl_int_is_zero75.0k
(a))
{0
907
0
    isl_int_clear(a);
908
0
    isl_assert(left->ctx, 0, goto error);
909
0
  }
910
294k
  
for (row = 0; 75.0k
row < left->n_row294k
;
++row219k
)
{219k
911
219k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
912
219k
    if (isl_int_is_one(left->row[row][row]))
913
215k
      continue;
914
4.01k
    isl_seq_scale(right->row[row], right->row[row],
915
4.01k
        left->row[row][row], right->n_col);
916
4.01k
  }
917
75.0k
  isl_int_clear(a);
918
75.0k
919
75.0k
  isl_mat_free(left);
920
75.0k
  return right;
921
0
error:
922
0
  isl_mat_free(left);
923
0
  isl_mat_free(right);
924
0
  return NULL;
925
75.0k
}
926
927
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
928
54.4k
{
929
54.4k
  int i;
930
54.4k
931
254k
  for (i = 0; 
i < mat->n_row254k
;
++i200k
)
932
200k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
933
54.4k
}
934
935
void isl_mat_col_combine(struct isl_mat *mat, unsigned dst,
936
  isl_int m1, unsigned src1, isl_int m2, unsigned src2)
937
155k
{
938
155k
  int i;
939
155k
  isl_int tmp;
940
155k
941
155k
  isl_int_init(tmp);
942
694k
  for (i = 0; 
i < mat->n_row694k
;
++i538k
)
{538k
943
538k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
944
538k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
945
538k
    isl_int_set(mat->row[i][dst], tmp);
946
538k
  }
947
155k
  isl_int_clear(tmp);
948
155k
}
949
950
struct isl_mat *isl_mat_right_inverse(struct isl_mat *mat)
951
50.6k
{
952
50.6k
  struct isl_mat *inv;
953
50.6k
  int row;
954
50.6k
  isl_int a, b;
955
50.6k
956
50.6k
  mat = isl_mat_cow(mat);
957
50.6k
  if (!mat)
958
0
    return NULL;
959
50.6k
960
50.6k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
961
50.6k
  inv = isl_mat_cow(inv);
962
50.6k
  if (!inv)
963
0
    goto error;
964
50.6k
965
50.6k
  
isl_int_init50.6k
(a);50.6k
966
50.6k
  isl_int_init(b);
967
196k
  for (row = 0; 
row < mat->n_row196k
;
++row145k
)
{145k
968
145k
    int pivot, first, i, off;
969
145k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
970
145k
    if (
pivot < 0145k
)
{0
971
0
      isl_int_clear(a);
972
0
      isl_int_clear(b);
973
0
      isl_assert(mat->ctx, pivot >= 0, goto error);
974
0
    }
975
145k
    pivot += row;
976
145k
    if (pivot != row)
977
21.3k
      exchange(mat, &inv, NULL, row, pivot, row);
978
145k
    if (isl_int_is_neg(mat->row[row][row]))
979
33.3k
      oppose(mat, &inv, NULL, row, row);
980
145k
    first = row+1;
981
401k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
982
256k
                mat->n_col-first)) != -1) {
983
256k
      first += off;
984
256k
      isl_int_fdiv_q(a, mat->row[row][first],
985
256k
                mat->row[row][row]);
986
256k
      subtract(mat, &inv, NULL, row, row, first, a);
987
256k
      if (
!256k
isl_int_is_zero256k
(mat->row[row][first]))
988
195k
        exchange(mat, &inv, NULL, row, row, first);
989
256k
      else
990
60.4k
        ++first;
991
256k
    }
992
309k
    for (i = 0; 
i < row309k
;
++i164k
)
{164k
993
164k
      if (isl_int_is_zero(mat->row[row][i]))
994
86.6k
        continue;
995
77.7k
      
isl_int_gcd77.7k
(a, mat->row[row][row], mat->row[row][i]);77.7k
996
77.7k
      isl_int_divexact(b, mat->row[row][i], a);
997
77.7k
      isl_int_divexact(a, mat->row[row][row], a);
998
77.7k
      isl_int_neg(a, a);
999
77.7k
      isl_mat_col_combine(mat, i, a, i, b, row);
1000
77.7k
      isl_mat_col_combine(inv, i, a, i, b, row);
1001
77.7k
    }
1002
145k
  }
1003
50.6k
  
isl_int_clear50.6k
(b);50.6k
1004
50.6k
1005
50.6k
  isl_int_set(a, mat->row[0][0]);
1006
145k
  for (row = 1; 
row < mat->n_row145k
;
++row94.9k
)
1007
94.9k
    isl_int_lcm(a, a, mat->row[row][row]);
1008
50.6k
  if (
isl_int_is_zero50.6k
(a))
{0
1009
0
    isl_int_clear(a);
1010
0
    goto error;
1011
0
  }
1012
196k
  
for (row = 0; 50.6k
row < mat->n_row196k
;
++row145k
)
{145k
1013
145k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1014
145k
    if (isl_int_is_one(mat->row[row][row]))
1015
91.1k
      continue;
1016
54.4k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1017
54.4k
  }
1018
50.6k
  isl_int_clear(a);
1019
50.6k
1020
50.6k
  isl_mat_free(mat);
1021
50.6k
1022
50.6k
  return inv;
1023
0
error:
1024
0
  isl_mat_free(mat);
1025
0
  isl_mat_free(inv);
1026
0
  return NULL;
1027
50.6k
}
1028
1029
struct isl_mat *isl_mat_transpose(struct isl_mat *mat)
1030
3.30k
{
1031
3.30k
  struct isl_mat *transpose = NULL;
1032
3.30k
  int i, j;
1033
3.30k
1034
3.30k
  if (!mat)
1035
0
    return NULL;
1036
3.30k
1037
3.30k
  
if (3.30k
mat->n_col == mat->n_row3.30k
)
{3.30k
1038
3.30k
    mat = isl_mat_cow(mat);
1039
3.30k
    if (!mat)
1040
0
      return NULL;
1041
13.4k
    
for (i = 0; 3.30k
i < mat->n_row13.4k
;
++i10.1k
)
1042
30.2k
      
for (j = i + 1; 10.1k
j < mat->n_col30.2k
;
++j20.1k
)
1043
20.1k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1044
3.30k
    return mat;
1045
3.30k
  }
1046
0
  transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
1047
0
  if (!transpose)
1048
0
    goto error;
1049
0
  
for (i = 0; 0
i < mat->n_row0
;
++i0
)
1050
0
    
for (j = 0; 0
j < mat->n_col0
;
++j0
)
1051
0
      isl_int_set(transpose->row[j][i], mat->row[i][j]);
1052
0
  isl_mat_free(mat);
1053
0
  return transpose;
1054
0
error:
1055
0
  isl_mat_free(mat);
1056
0
  return NULL;
1057
0
}
1058
1059
struct isl_mat *isl_mat_swap_cols(struct isl_mat *mat, unsigned i, unsigned j)
1060
324k
{
1061
324k
  int r;
1062
324k
1063
324k
  mat = isl_mat_cow(mat);
1064
324k
  if (!mat)
1065
0
    return NULL;
1066
324k
  
isl_assert324k
(mat->ctx, i < mat->n_col, goto error);324k
1067
324k
  
isl_assert324k
(mat->ctx, j < mat->n_col, goto error);324k
1068
324k
1069
8.20M
  
for (r = 0; 324k
r < mat->n_row8.20M
;
++r7.87M
)
1070
7.87M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1071
324k
  return mat;
1072
0
error:
1073
0
  isl_mat_free(mat);
1074
0
  return NULL;
1075
324k
}
1076
1077
struct isl_mat *isl_mat_swap_rows(struct isl_mat *mat, unsigned i, unsigned j)
1078
978k
{
1079
978k
  isl_int *t;
1080
978k
1081
978k
  if (!mat)
1082
0
    return NULL;
1083
978k
  mat = isl_mat_cow(mat);
1084
978k
  if (!mat)
1085
0
    return NULL;
1086
978k
  t = mat->row[i];
1087
978k
  mat->row[i] = mat->row[j];
1088
978k
  mat->row[j] = t;
1089
978k
  return mat;
1090
978k
}
1091
1092
/* Calculate the product of two matrices.
1093
 *
1094
 * This function is optimized for operand matrices that contain many zeros and
1095
 * skips multiplications where we know one of the operands is zero.
1096
 */
1097
__isl_give isl_mat *isl_mat_product(__isl_take isl_mat *left,
1098
  __isl_take isl_mat *right)
1099
1.19M
{
1100
1.19M
  int i, j, k;
1101
1.19M
  struct isl_mat *prod;
1102
1.19M
1103
1.19M
  if (
!left || 1.19M
!right1.19M
)
1104
0
    goto error;
1105
1.19M
  
isl_assert1.19M
(left->ctx, left->n_col == right->n_row, goto error);1.19M
1106
1.19M
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1107
1.19M
  if (!prod)
1108
0
    goto error;
1109
1.19M
  
if (1.19M
left->n_col == 01.19M
)
{231
1110
363
    for (i = 0; 
i < prod->n_row363
;
++i132
)
1111
132
      isl_seq_clr(prod->row[i], prod->n_col);
1112
231
    isl_mat_free(left);
1113
231
    isl_mat_free(right);
1114
231
    return prod;
1115
231
  }
1116
3.46M
  
for (i = 0; 1.19M
i < prod->n_row3.46M
;
++i2.26M
)
{2.26M
1117
17.4M
    for (j = 0; 
j < prod->n_col17.4M
;
++j15.1M
)
1118
15.1M
      isl_int_mul(prod->row[i][j],
1119
2.26M
            left->row[i][0], right->row[0][j]);
1120
19.8M
    for (k = 1; 
k < left->n_col19.8M
;
++k17.5M
)
{17.5M
1121
17.5M
      if (isl_int_is_zero(left->row[i][k]))
1122
14.9M
        continue;
1123
22.7M
      
for (j = 0; 2.66M
j < prod->n_col22.7M
;
++j20.1M
)
1124
20.1M
        isl_int_addmul(prod->row[i][j],
1125
2.66M
              left->row[i][k], right->row[k][j]);
1126
2.66M
    }
1127
2.26M
  }
1128
1.19M
  isl_mat_free(left);
1129
1.19M
  isl_mat_free(right);
1130
1.19M
  return prod;
1131
0
error:
1132
0
  isl_mat_free(left);
1133
0
  isl_mat_free(right);
1134
0
  return NULL;
1135
1.19M
}
1136
1137
/* Replace the variables x in the rows q by x' given by x = M x',
1138
 * with M the matrix mat.
1139
 *
1140
 * If the number of new variables is greater than the original
1141
 * number of variables, then the rows q have already been
1142
 * preextended.  If the new number is smaller, then the coefficients
1143
 * of the divs, which are not changed, need to be shifted down.
1144
 * The row q may be the equalities, the inequalities or the
1145
 * div expressions.  In the latter case, has_div is true and
1146
 * we need to take into account the extra denominator column.
1147
 */
1148
static int preimage(struct isl_ctx *ctx, isl_int **q, unsigned n,
1149
  unsigned n_div, int has_div, struct isl_mat *mat)
1150
1.04M
{
1151
1.04M
  int i;
1152
1.04M
  struct isl_mat *t;
1153
1.04M
  int e;
1154
1.04M
1155
1.04M
  if (mat->n_col >= mat->n_row)
1156
689k
    e = 0;
1157
1.04M
  else
1158
351k
    e = mat->n_row - mat->n_col;
1159
1.04M
  if (has_div)
1160
346k
    
for (i = 0; 346k
i < n346k
;
++i36
)
1161
36
      isl_int_mul(q[i][0], q[i][0], mat->row[0][0]);
1162
1.04M
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1163
1.04M
  t = isl_mat_product(t, mat);
1164
1.04M
  if (!t)
1165
0
    return -1;
1166
2.55M
  
for (i = 0; 1.04M
i < n2.55M
;
++i1.51M
)
{1.51M
1167
1.51M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1168
1.51M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1169
1.51M
          q[i] + has_div + t->n_col + e, n_div);
1170
1.51M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1171
1.51M
  }
1172
1.04M
  isl_mat_free(t);
1173
1.04M
  return 0;
1174
1.04M
}
1175
1176
/* Replace the variables x in bset by x' given by x = M x', with
1177
 * M the matrix mat.
1178
 *
1179
 * If there are fewer variables x' then there are x, then we perform
1180
 * the transformation in place, which means that, in principle,
1181
 * this frees up some extra variables as the number
1182
 * of columns remains constant, but we would have to extend
1183
 * the div array too as the number of rows in this array is assumed
1184
 * to be equal to extra.
1185
 */
1186
struct isl_basic_set *isl_basic_set_preimage(struct isl_basic_set *bset,
1187
  struct isl_mat *mat)
1188
346k
{
1189
346k
  struct isl_ctx *ctx;
1190
346k
1191
346k
  if (
!bset || 346k
!mat346k
)
1192
0
    goto error;
1193
346k
1194
346k
  ctx = bset->ctx;
1195
346k
  bset = isl_basic_set_cow(bset);
1196
346k
  if (!bset)
1197
0
    goto error;
1198
346k
1199
346k
  
isl_assert346k
(ctx, bset->dim->nparam == 0, goto error);346k
1200
346k
  
isl_assert346k
(ctx, 1+bset->dim->n_out == mat->n_row, goto error);346k
1201
346k
  
isl_assert346k
(ctx, mat->n_col > 0, goto error);346k
1202
346k
1203
346k
  
if (346k
mat->n_col > mat->n_row346k
)
{59.6k
1204
59.6k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1205
59.6k
    if (!bset)
1206
0
      goto error;
1207
287k
  } else 
if (287k
mat->n_col < mat->n_row287k
)
{117k
1208
117k
    bset->dim = isl_space_cow(bset->dim);
1209
117k
    if (!bset->dim)
1210
0
      goto error;
1211
117k
    bset->dim->n_out -= mat->n_row - mat->n_col;
1212
117k
  }
1213
346k
1214
346k
  
if (346k
preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,346k
1215
346k
      isl_mat_copy(mat)) < 0)
1216
0
    goto error;
1217
346k
1218
346k
  
if (346k
preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,346k
1219
346k
      isl_mat_copy(mat)) < 0)
1220
0
    goto error;
1221
346k
1222
346k
  
if (346k
preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0346k
)
1223
0
    goto error2;
1224
346k
1225
346k
  
ISL_F_CLR346k
(bset, ISL_BASIC_SET_NO_IMPLICIT);346k
1226
346k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1227
346k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1228
346k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1229
346k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1230
346k
1231
346k
  bset = isl_basic_set_simplify(bset);
1232
346k
  bset = isl_basic_set_finalize(bset);
1233
346k
1234
346k
  return bset;
1235
0
error:
1236
0
  isl_mat_free(mat);
1237
0
error2:
1238
0
  isl_basic_set_free(bset);
1239
0
  return NULL;
1240
0
}
1241
1242
struct isl_set *isl_set_preimage(struct isl_set *set, struct isl_mat *mat)
1243
28.2k
{
1244
28.2k
  int i;
1245
28.2k
1246
28.2k
  set = isl_set_cow(set);
1247
28.2k
  if (!set)
1248
0
    goto error;
1249
28.2k
1250
83.5k
  
for (i = 0; 28.2k
i < set->n83.5k
;
++i55.2k
)
{55.2k
1251
55.2k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1252
55.2k
                isl_mat_copy(mat));
1253
55.2k
    if (!set->p[i])
1254
0
      goto error;
1255
55.2k
  }
1256
28.2k
  
if (28.2k
mat->n_col != mat->n_row28.2k
)
{13.6k
1257
13.6k
    set->dim = isl_space_cow(set->dim);
1258
13.6k
    if (!set->dim)
1259
0
      goto error;
1260
13.6k
    set->dim->n_out += mat->n_col;
1261
13.6k
    set->dim->n_out -= mat->n_row;
1262
13.6k
  }
1263
28.2k
  isl_mat_free(mat);
1264
28.2k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1265
28.2k
  return set;
1266
0
error:
1267
0
  isl_set_free(set);
1268
0
  isl_mat_free(mat);
1269
0
  return NULL;
1270
28.2k
}
1271
1272
/* Replace the variables x starting at pos in the rows q
1273
 * by x' with x = M x' with M the matrix mat.
1274
 * That is, replace the corresponding coefficients c by c M.
1275
 */
1276
static int transform(isl_ctx *ctx, isl_int **q, unsigned n,
1277
  unsigned pos, __isl_take isl_mat *mat)
1278
4.08k
{
1279
4.08k
  int i;
1280
4.08k
  isl_mat *t;
1281
4.08k
1282
4.08k
  t = isl_mat_sub_alloc6(ctx, q, 0, n, pos, mat->n_row);
1283
4.08k
  t = isl_mat_product(t, mat);
1284
4.08k
  if (!t)
1285
0
    return -1;
1286
8.16k
  
for (i = 0; 4.08k
i < n8.16k
;
++i4.07k
)
1287
4.07k
    isl_seq_swp_or_cpy(q[i] + pos, t->row[i], t->n_col);
1288
4.08k
  isl_mat_free(t);
1289
4.08k
  return 0;
1290
4.08k
}
1291
1292
/* Replace the variables x of type "type" starting at "first" in "bmap"
1293
 * by x' with x = M x' with M the matrix trans.
1294
 * That is, replace the corresponding coefficients c by c M.
1295
 *
1296
 * The transformation matrix should be a square matrix.
1297
 */
1298
__isl_give isl_basic_map *isl_basic_map_transform_dims(
1299
  __isl_take isl_basic_map *bmap, enum isl_dim_type type, unsigned first,
1300
  __isl_take isl_mat *trans)
1301
1.36k
{
1302
1.36k
  isl_ctx *ctx;
1303
1.36k
  unsigned pos;
1304
1.36k
1305
1.36k
  bmap = isl_basic_map_cow(bmap);
1306
1.36k
  if (
!bmap || 1.36k
!trans1.36k
)
1307
0
    goto error;
1308
1.36k
1309
1.36k
  ctx = isl_basic_map_get_ctx(bmap);
1310
1.36k
  if (trans->n_row != trans->n_col)
1311
0
    isl_die(trans->ctx, isl_error_invalid,
1312
1.36k
      "expecting square transformation matrix", goto error);
1313
1.36k
  
if (1.36k
first + trans->n_row > isl_basic_map_dim(bmap, type)1.36k
)
1314
0
    isl_die(trans->ctx, isl_error_invalid,
1315
1.36k
      "oversized transformation matrix", goto error);
1316
1.36k
1317
1.36k
  pos = isl_basic_map_offset(bmap, type) + first;
1318
1.36k
1319
1.36k
  if (transform(ctx, bmap->eq, bmap->n_eq, pos, isl_mat_copy(trans)) < 0)
1320
0
    goto error;
1321
1.36k
  
if (1.36k
transform(ctx, bmap->ineq, bmap->n_ineq, pos,1.36k
1322
1.36k
          isl_mat_copy(trans)) < 0)
1323
0
    goto error;
1324
1.36k
  
if (1.36k
transform(ctx, bmap->div, bmap->n_div, 1 + pos,1.36k
1325
1.36k
          isl_mat_copy(trans)) < 0)
1326
0
    goto error;
1327
1.36k
1328
1.36k
  
ISL_F_CLR1.36k
(bmap, ISL_BASIC_MAP_NORMALIZED);1.36k
1329
1.36k
  ISL_F_CLR(bmap, ISL_BASIC_MAP_NORMALIZED_DIVS);
1330
1.36k
1331
1.36k
  isl_mat_free(trans);
1332
1.36k
  return bmap;
1333
0
error:
1334
0
  isl_mat_free(trans);
1335
0
  isl_basic_map_free(bmap);
1336
0
  return NULL;
1337
1.36k
}
1338
1339
/* Replace the variables x of type "type" starting at "first" in "bset"
1340
 * by x' with x = M x' with M the matrix trans.
1341
 * That is, replace the corresponding coefficients c by c M.
1342
 *
1343
 * The transformation matrix should be a square matrix.
1344
 */
1345
__isl_give isl_basic_set *isl_basic_set_transform_dims(
1346
  __isl_take isl_basic_set *bset, enum isl_dim_type type, unsigned first,
1347
  __isl_take isl_mat *trans)
1348
1.31k
{
1349
1.31k
  return isl_basic_map_transform_dims(bset, type, first, trans);
1350
1.31k
}
1351
1352
void isl_mat_print_internal(__isl_keep isl_mat *mat, FILE *out, int indent)
1353
0
{
1354
0
  int i, j;
1355
0
1356
0
  if (
!mat0
)
{0
1357
0
    fprintf(out, "%*snull mat\n", indent, "");
1358
0
    return;
1359
0
  }
1360
0
1361
0
  
if (0
mat->n_row == 00
)
1362
0
    fprintf(out, "%*s[]\n", indent, "");
1363
0
1364
0
  for (i = 0; 
i < mat->n_row0
;
++i0
)
{0
1365
0
    if (!i)
1366
0
      fprintf(out, "%*s[[", indent, "");
1367
0
    else
1368
0
      fprintf(out, "%*s[", indent+1, "");
1369
0
    for (j = 0; 
j < mat->n_col0
;
++j0
)
{0
1370
0
      if (j)
1371
0
          fprintf(out, ",");
1372
0
      isl_int_print(out, mat->row[i][j], 0);
1373
0
    }
1374
0
    if (i == mat->n_row-1)
1375
0
      fprintf(out, "]]\n");
1376
0
    else
1377
0
      fprintf(out, "]\n");
1378
0
  }
1379
0
}
1380
1381
void isl_mat_dump(__isl_keep isl_mat *mat)
1382
0
{
1383
0
  isl_mat_print_internal(mat, stderr, 0);
1384
0
}
1385
1386
struct isl_mat *isl_mat_drop_cols(struct isl_mat *mat, unsigned col, unsigned n)
1387
54.8k
{
1388
54.8k
  int r;
1389
54.8k
1390
54.8k
  if (n == 0)
1391
2
    return mat;
1392
54.8k
1393
54.8k
  mat = isl_mat_cow(mat);
1394
54.8k
  if (!mat)
1395
0
    return NULL;
1396
54.8k
1397
54.8k
  
if (54.8k
col != mat->n_col-n54.8k
)
{14.0k
1398
64.1k
    for (r = 0; 
r < mat->n_row64.1k
;
++r50.0k
)
1399
50.0k
      isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1400
50.0k
          mat->n_col - col - n);
1401
14.0k
  }
1402
54.8k
  mat->n_col -= n;
1403
54.8k
  return mat;
1404
54.8k
}
1405
1406
struct isl_mat *isl_mat_drop_rows(struct isl_mat *mat, unsigned row, unsigned n)
1407
65.3k
{
1408
65.3k
  int r;
1409
65.3k
1410
65.3k
  mat = isl_mat_cow(mat);
1411
65.3k
  if (!mat)
1412
0
    return NULL;
1413
65.3k
1414
155k
  
for (r = row; 65.3k
r+n < mat->n_row155k
;
++r90.3k
)
1415
90.3k
    mat->row[r] = mat->row[r+n];
1416
65.3k
1417
65.3k
  mat->n_row -= n;
1418
65.3k
  return mat;
1419
65.3k
}
1420
1421
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1422
        unsigned col, unsigned n)
1423
37.6k
{
1424
37.6k
  isl_mat *ext;
1425
37.6k
1426
37.6k
  if (!mat)
1427
0
    return NULL;
1428
37.6k
  
if (37.6k
n == 037.6k
)
1429
87
    return mat;
1430
37.6k
1431
37.5k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1432
37.5k
  if (!ext)
1433
0
    goto error;
1434
37.5k
1435
37.5k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1436
37.5k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1437
37.5k
        col + n, col, mat->n_col - col);
1438
37.5k
1439
37.5k
  isl_mat_free(mat);
1440
37.5k
  return ext;
1441
0
error:
1442
0
  isl_mat_free(mat);
1443
0
  return NULL;
1444
37.5k
}
1445
1446
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1447
  unsigned first, unsigned n)
1448
37.6k
{
1449
37.6k
  int i;
1450
37.6k
1451
37.6k
  if (!mat)
1452
0
    return NULL;
1453
37.6k
  mat = isl_mat_insert_cols(mat, first, n);
1454
37.6k
  if (!mat)
1455
0
    return NULL;
1456
37.6k
1457
40.5k
  
for (i = 0; 37.6k
i < mat->n_row40.5k
;
++i2.85k
)
1458
2.85k
    isl_seq_clr(mat->row[i] + first, n);
1459
37.6k
1460
37.6k
  return mat;
1461
37.6k
}
1462
1463
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1464
3.72k
{
1465
3.72k
  if (!mat)
1466
0
    return NULL;
1467
3.72k
1468
3.72k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1469
3.72k
}
1470
1471
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1472
        unsigned row, unsigned n)
1473
4.13k
{
1474
4.13k
  isl_mat *ext;
1475
4.13k
1476
4.13k
  if (!mat)
1477
0
    return NULL;
1478
4.13k
  
if (4.13k
n == 04.13k
)
1479
85
    return mat;
1480
4.13k
1481
4.04k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1482
4.04k
  if (!ext)
1483
0
    goto error;
1484
4.04k
1485
4.04k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1486
4.04k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1487
4.04k
        mat->n_row - row, 0, 0, mat->n_col);
1488
4.04k
1489
4.04k
  isl_mat_free(mat);
1490
4.04k
  return ext;
1491
0
error:
1492
0
  isl_mat_free(mat);
1493
0
  return NULL;
1494
4.04k
}
1495
1496
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1497
4.12k
{
1498
4.12k
  if (!mat)
1499
0
    return NULL;
1500
4.12k
1501
4.12k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1502
4.12k
}
1503
1504
__isl_give isl_mat *isl_mat_insert_zero_rows(__isl_take isl_mat *mat,
1505
  unsigned row, unsigned n)
1506
0
{
1507
0
  int i;
1508
0
1509
0
  mat = isl_mat_insert_rows(mat, row, n);
1510
0
  if (!mat)
1511
0
    return NULL;
1512
0
  
1513
0
  
for (i = 0; 0
i < n0
;
++i0
)
1514
0
    isl_seq_clr(mat->row[row + i], mat->n_col);
1515
0
1516
0
  return mat;
1517
0
}
1518
1519
__isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n)
1520
0
{
1521
0
  if (!mat)
1522
0
    return NULL;
1523
0
1524
0
  return isl_mat_insert_zero_rows(mat, mat->n_row, n);
1525
0
}
1526
1527
void isl_mat_col_submul(struct isl_mat *mat,
1528
      int dst_col, isl_int f, int src_col)
1529
2.03k
{
1530
2.03k
  int i;
1531
2.03k
1532
9.59k
  for (i = 0; 
i < mat->n_row9.59k
;
++i7.55k
)
1533
7.55k
    isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1534
2.03k
}
1535
1536
void isl_mat_col_add(__isl_keep isl_mat *mat, int dst_col, int src_col)
1537
2
{
1538
2
  int i;
1539
2
1540
2
  if (!mat)
1541
0
    return;
1542
2
1543
7
  
for (i = 0; 2
i < mat->n_row7
;
++i5
)
1544
5
    isl_int_add(mat->row[i][dst_col],
1545
2
          mat->row[i][dst_col], mat->row[i][src_col]);
1546
2
}
1547
1548
void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1549
6.31k
{
1550
6.31k
  int i;
1551
6.31k
1552
23.3k
  for (i = 0; 
i < mat->n_row23.3k
;
++i16.9k
)
1553
16.9k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1554
6.31k
}
1555
1556
/* Add "f" times column "src_col" to column "dst_col" of "mat" and
1557
 * return the result.
1558
 */
1559
__isl_give isl_mat *isl_mat_col_addmul(__isl_take isl_mat *mat, int dst_col,
1560
  isl_int f, int src_col)
1561
6
{
1562
6
  int i;
1563
6
1564
6
  if (
check_col(mat, dst_col) < 0 || 6
check_col(mat, src_col) < 06
)
1565
0
    return isl_mat_free(mat);
1566
6
1567
14
  
for (i = 0; 6
i < mat->n_row14
;
++i8
)
{8
1568
8
    if (isl_int_is_zero(mat->row[i][src_col]))
1569
2
      continue;
1570
6
    mat = isl_mat_cow(mat);
1571
6
    if (!mat)
1572
0
      return NULL;
1573
6
    
isl_int_addmul6
(mat->row[i][dst_col], f, mat->row[i][src_col]);6
1574
6
  }
1575
6
1576
6
  return mat;
1577
6
}
1578
1579
/* Negate column "col" of "mat" and return the result.
1580
 */
1581
__isl_give isl_mat *isl_mat_col_neg(__isl_take isl_mat *mat, int col)
1582
4
{
1583
4
  int i;
1584
4
1585
4
  if (check_col(mat, col) < 0)
1586
0
    return isl_mat_free(mat);
1587
4
1588
10
  
for (i = 0; 4
i < mat->n_row10
;
++i6
)
{6
1589
6
    if (isl_int_is_zero(mat->row[i][col]))
1590
2
      continue;
1591
4
    mat = isl_mat_cow(mat);
1592
4
    if (!mat)
1593
0
      return NULL;
1594
4
    
isl_int_neg4
(mat->row[i][col], mat->row[i][col]);4
1595
4
  }
1596
4
1597
4
  return mat;
1598
4
}
1599
1600
struct isl_mat *isl_mat_unimodular_complete(struct isl_mat *M, int row)
1601
10.6k
{
1602
10.6k
  int r;
1603
10.6k
  struct isl_mat *H = NULL, *Q = NULL;
1604
10.6k
1605
10.6k
  if (!M)
1606
0
    return NULL;
1607
10.6k
1608
10.6k
  
isl_assert10.6k
(M->ctx, M->n_row == M->n_col, goto error);10.6k
1609
10.6k
  M->n_row = row;
1610
10.6k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1611
10.6k
  M->n_row = M->n_col;
1612
10.6k
  if (!H)
1613
0
    goto error;
1614
21.3k
  
for (r = 0; 10.6k
r < row21.3k
;
++r10.6k
)
1615
10.6k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1616
32.5k
  
for (r = row; 10.6k
r < M->n_row32.5k
;
++r21.8k
)
1617
21.8k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1618
10.6k
  isl_mat_free(H);
1619
10.6k
  isl_mat_free(Q);
1620
10.6k
  return M;
1621
0
error:
1622
0
  isl_mat_free(H);
1623
0
  isl_mat_free(Q);
1624
0
  isl_mat_free(M);
1625
0
  return NULL;
1626
10.6k
}
1627
1628
__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1629
  __isl_take isl_mat *bot)
1630
368
{
1631
368
  struct isl_mat *mat;
1632
368
1633
368
  if (
!top || 368
!bot368
)
1634
0
    goto error;
1635
368
1636
368
  
isl_assert368
(top->ctx, top->n_col == bot->n_col, goto error);368
1637
368
  
if (368
top->n_row == 0368
)
{168
1638
168
    isl_mat_free(top);
1639
168
    return bot;
1640
168
  }
1641
200
  
if (200
bot->n_row == 0200
)
{0
1642
0
    isl_mat_free(bot);
1643
0
    return top;
1644
0
  }
1645
200
1646
200
  mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1647
200
  if (!mat)
1648
0
    goto error;
1649
200
  isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1650
200
       0, 0, mat->n_col);
1651
200
  isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1652
200
       0, 0, mat->n_col);
1653
200
  isl_mat_free(top);
1654
200
  isl_mat_free(bot);
1655
200
  return mat;
1656
0
error:
1657
0
  isl_mat_free(top);
1658
0
  isl_mat_free(bot);
1659
0
  return NULL;
1660
200
}
1661
1662
isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1663
5.36k
{
1664
5.36k
  int i;
1665
5.36k
1666
5.36k
  if (
!mat1 || 5.36k
!mat25.36k
)
1667
0
    return isl_bool_error;
1668
5.36k
1669
5.36k
  
if (5.36k
mat1->n_row != mat2->n_row5.36k
)
1670
781
    return isl_bool_false;
1671
5.36k
1672
4.58k
  
if (4.58k
mat1->n_col != mat2->n_col4.58k
)
1673
0
    return isl_bool_false;
1674
4.58k
1675
8.39k
  
for (i = 0; 4.58k
i < mat1->n_row8.39k
;
++i3.81k
)
1676
4.04k
    
if (4.04k
!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col)4.04k
)
1677
236
      return isl_bool_false;
1678
4.58k
1679
4.34k
  return isl_bool_true;
1680
4.58k
}
1681
1682
__isl_give isl_mat *isl_mat_from_row_vec(__isl_take isl_vec *vec)
1683
0
{
1684
0
  struct isl_mat *mat;
1685
0
1686
0
  if (!vec)
1687
0
    return NULL;
1688
0
  mat = isl_mat_alloc(vec->ctx, 1, vec->size);
1689
0
  if (!mat)
1690
0
    goto error;
1691
0
1692
0
  isl_seq_cpy(mat->row[0], vec->el, vec->size);
1693
0
1694
0
  isl_vec_free(vec);
1695
0
  return mat;
1696
0
error:
1697
0
  isl_vec_free(vec);
1698
0
  return NULL;
1699
0
}
1700
1701
/* Return a copy of row "row" of "mat" as an isl_vec.
1702
 */
1703
__isl_give isl_vec *isl_mat_get_row(__isl_keep isl_mat *mat, unsigned row)
1704
12
{
1705
12
  isl_vec *v;
1706
12
1707
12
  if (!mat)
1708
0
    return NULL;
1709
12
  
if (12
row >= mat->n_row12
)
1710
0
    isl_die(mat->ctx, isl_error_invalid, "row out of range",
1711
12
      return NULL);
1712
12
1713
12
  v = isl_vec_alloc(isl_mat_get_ctx(mat), mat->n_col);
1714
12
  if (!v)
1715
0
    return NULL;
1716
12
  isl_seq_cpy(v->el, mat->row[row], mat->n_col);
1717
12
1718
12
  return v;
1719
12
}
1720
1721
__isl_give isl_mat *isl_mat_vec_concat(__isl_take isl_mat *top,
1722
  __isl_take isl_vec *bot)
1723
0
{
1724
0
  return isl_mat_concat(top, isl_mat_from_row_vec(bot));
1725
0
}
1726
1727
__isl_give isl_mat *isl_mat_move_cols(__isl_take isl_mat *mat,
1728
  unsigned dst_col, unsigned src_col, unsigned n)
1729
388
{
1730
388
  isl_mat *res;
1731
388
1732
388
  if (!mat)
1733
0
    return NULL;
1734
388
  
if (388
n == 0 || 388
dst_col == src_col388
)
1735
380
    return mat;
1736
388
1737
8
  res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1738
8
  if (!res)
1739
0
    goto error;
1740
8
1741
8
  
if (8
dst_col < src_col8
)
{8
1742
8
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1743
8
         0, 0, dst_col);
1744
8
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1745
8
         dst_col, src_col, n);
1746
8
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1747
8
         dst_col + n, dst_col, src_col - dst_col);
1748
8
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1749
8
         src_col + n, src_col + n,
1750
8
         res->n_col - src_col - n);
1751
0
  } else {
1752
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1753
0
         0, 0, src_col);
1754
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1755
0
         src_col, src_col + n, dst_col - src_col);
1756
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1757
0
         dst_col, src_col, n);
1758
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1759
0
         dst_col + n, dst_col + n,
1760
0
         res->n_col - dst_col - n);
1761
0
  }
1762
8
  isl_mat_free(mat);
1763
8
1764
8
  return res;
1765
0
error:
1766
0
  isl_mat_free(mat);
1767
0
  return NULL;
1768
8
}
1769
1770
/* Return the gcd of the elements in row "row" of "mat" in *gcd.
1771
 * Return isl_stat_ok on success and isl_stat_error on failure.
1772
 */
1773
isl_stat isl_mat_row_gcd(__isl_keep isl_mat *mat, int row, isl_int *gcd)
1774
7
{
1775
7
  if (!mat)
1776
0
    return isl_stat_error;
1777
7
1778
7
  
if (7
row < 0 || 7
row >= mat->n_row7
)
1779
0
    isl_die(isl_mat_get_ctx(mat), isl_error_invalid,
1780
7
      "row out of range", return isl_stat_error);
1781
7
  isl_seq_gcd(mat->row[row], mat->n_col, gcd);
1782
7
1783
7
  return isl_stat_ok;
1784
7
}
1785
1786
void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1787
411
{
1788
411
  int i;
1789
411
  isl_int g;
1790
411
1791
411
  isl_int_set_si(*gcd, 0);
1792
411
  if (!mat)
1793
0
    return;
1794
411
1795
411
  
isl_int_init411
(g);411
1796
2.61k
  for (i = 0; 
i < mat->n_row2.61k
;
++i2.20k
)
{2.20k
1797
2.20k
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1798
2.20k
    isl_int_gcd(*gcd, *gcd, g);
1799
2.20k
  }
1800
411
  isl_int_clear(g);
1801
411
}
1802
1803
/* Return the result of scaling "mat" by a factor of "m".
1804
 */
1805
__isl_give isl_mat *isl_mat_scale(__isl_take isl_mat *mat, isl_int m)
1806
0
{
1807
0
  int i;
1808
0
1809
0
  if (isl_int_is_one(m))
1810
0
    return mat;
1811
0
1812
0
  mat = isl_mat_cow(mat);
1813
0
  if (!mat)
1814
0
    return NULL;
1815
0
1816
0
  
for (i = 0; 0
i < mat->n_row0
;
++i0
)
1817
0
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
1818
0
1819
0
  return mat;
1820
0
}
1821
1822
__isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1823
411
{
1824
411
  int i;
1825
411
1826
411
  if (isl_int_is_one(m))
1827
0
    return mat;
1828
411
1829
411
  mat = isl_mat_cow(mat);
1830
411
  if (!mat)
1831
0
    return NULL;
1832
411
1833
2.61k
  
for (i = 0; 411
i < mat->n_row2.61k
;
++i2.20k
)
1834
2.20k
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1835
411
1836
411
  return mat;
1837
411
}
1838
1839
__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row,
1840
  isl_int m)
1841
14
{
1842
14
  if (isl_int_is_one(m))
1843
0
    return mat;
1844
14
1845
14
  mat = isl_mat_cow(mat);
1846
14
  if (!mat)
1847
0
    return NULL;
1848
14
1849
14
  isl_seq_scale_down(mat->row[row], mat->row[row], m, mat->n_col);
1850
14
1851
14
  return mat;
1852
14
}
1853
1854
__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1855
411
{
1856
411
  isl_int gcd;
1857
411
1858
411
  if (!mat)
1859
0
    return NULL;
1860
411
1861
411
  
isl_int_init411
(gcd);411
1862
411
  isl_mat_gcd(mat, &gcd);
1863
411
  mat = isl_mat_scale_down(mat, gcd);
1864
411
  isl_int_clear(gcd);
1865
411
1866
411
  return mat;
1867
411
}
1868
1869
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1870
1.06k
{
1871
1.06k
  mat = isl_mat_cow(mat);
1872
1.06k
  if (!mat)
1873
0
    return NULL;
1874
1.06k
1875
1.06k
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
1876
1.06k
1877
1.06k
  return mat;
1878
1.06k
}
1879
1880
/* Number of initial non-zero columns.
1881
 */
1882
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
1883
933
{
1884
933
  int i;
1885
933
1886
933
  if (!mat)
1887
0
    return -1;
1888
933
1889
1.35k
  
for (i = 0; 933
i < mat->n_col1.35k
;
++i419
)
1890
1.16k
    
if (1.16k
row_first_non_zero(mat->row, mat->n_row, i) < 01.16k
)
1891
747
      break;
1892
933
1893
933
  return i;
1894
933
}