Coverage Report

Created: 2017-08-21 19:50

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