Coverage Report

Created: 2017-11-23 03:11

/Users/buildslave/jenkins/workspace/clang-stage2-coverage-R/llvm/tools/polly/lib/External/isl/isl_mat.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2008-2009 Katholieke Universiteit Leuven
3
 * Copyright 2014      Ecole Normale Superieure
4
 *
5
 * Use of this software is governed by the MIT license
6
 *
7
 * Written by Sven Verdoolaege, K.U.Leuven, Departement
8
 * Computerwetenschappen, Celestijnenlaan 200A, B-3001 Leuven, Belgium
9
 * and Ecole Normale Superieure, 45 rue d'Ulm, 75230 Paris, France
10
 */
11
12
#include <isl_ctx_private.h>
13
#include <isl_map_private.h>
14
#include <isl/space.h>
15
#include <isl_seq.h>
16
#include <isl_mat_private.h>
17
#include <isl_vec_private.h>
18
#include <isl_space_private.h>
19
#include <isl_val_private.h>
20
#include <isl/deprecated/mat_int.h>
21
22
isl_ctx *isl_mat_get_ctx(__isl_keep isl_mat *mat)
23
3.37M
{
24
3.37M
  return mat ? mat->ctx : NULL;
25
3.37M
}
26
27
/* Return a hash value that digests "mat".
28
 */
29
uint32_t isl_mat_get_hash(__isl_keep isl_mat *mat)
30
0
{
31
0
  int i;
32
0
  uint32_t hash;
33
0
34
0
  if (!mat)
35
0
    return 0;
36
0
37
0
  hash = isl_hash_init();
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_row; ++i) {
41
0
    uint32_t row_hash;
42
0
43
0
    row_hash = isl_seq_get_hash(mat->row[i], mat->n_col);
44
0
    isl_hash_hash(hash, row_hash);
45
0
  }
46
0
47
0
  return hash;
48
0
}
49
50
struct isl_mat *isl_mat_alloc(struct isl_ctx *ctx,
51
  unsigned n_row, unsigned n_col)
52
5.89M
{
53
5.89M
  int i;
54
5.89M
  struct isl_mat *mat;
55
5.89M
56
5.89M
  mat = isl_alloc_type(ctx, struct isl_mat);
57
5.89M
  if (!mat)
58
0
    return NULL;
59
5.89M
60
5.89M
  mat->row = NULL;
61
5.89M
  mat->block = isl_blk_alloc(ctx, n_row * n_col);
62
5.89M
  if (isl_blk_is_error(mat->block))
63
0
    goto error;
64
5.89M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
65
5.89M
  if (n_row && 
!mat->row4.36M
)
66
0
    goto error;
67
5.89M
68
32.2M
  
for (i = 0; 5.89M
i < n_row;
++i26.3M
)
69
26.3M
    mat->row[i] = mat->block.data + i * n_col;
70
5.89M
71
5.89M
  mat->ctx = ctx;
72
5.89M
  isl_ctx_ref(ctx);
73
5.89M
  mat->ref = 1;
74
5.89M
  mat->n_row = n_row;
75
5.89M
  mat->n_col = n_col;
76
5.89M
  mat->max_col = n_col;
77
5.89M
  mat->flags = 0;
78
5.89M
79
5.89M
  return mat;
80
0
error:
81
0
  isl_blk_free(ctx, mat->block);
82
0
  free(mat);
83
0
  return NULL;
84
5.89M
}
85
86
struct isl_mat *isl_mat_extend(struct isl_mat *mat,
87
  unsigned n_row, unsigned n_col)
88
110k
{
89
110k
  int i;
90
110k
  isl_int *old;
91
110k
  isl_int **row;
92
110k
93
110k
  if (!mat)
94
0
    return NULL;
95
110k
96
110k
  if (mat->max_col >= n_col && 
mat->n_row >= n_row106k
) {
97
8.22k
    if (mat->n_col < n_col)
98
22
      mat->n_col = n_col;
99
8.22k
    return mat;
100
8.22k
  }
101
102k
102
102k
  if (mat->max_col < n_col) {
103
3.67k
    struct isl_mat *new_mat;
104
3.67k
105
3.67k
    if (n_row < mat->n_row)
106
8
      n_row = mat->n_row;
107
3.67k
    new_mat = isl_mat_alloc(mat->ctx, n_row, n_col);
108
3.67k
    if (!new_mat)
109
0
      goto error;
110
25.7k
    
for (i = 0; 3.67k
i < mat->n_row;
++i22.0k
)
111
22.0k
      isl_seq_cpy(new_mat->row[i], mat->row[i], mat->n_col);
112
3.67k
    isl_mat_free(mat);
113
3.67k
    return new_mat;
114
3.67k
  }
115
98.6k
116
98.6k
  mat = isl_mat_cow(mat);
117
98.6k
  if (!mat)
118
0
    goto error;
119
98.6k
120
98.6k
  old = mat->block.data;
121
98.6k
  mat->block = isl_blk_extend(mat->ctx, mat->block, n_row * mat->max_col);
122
98.6k
  if (isl_blk_is_error(mat->block))
123
0
    goto error;
124
98.6k
  row = isl_realloc_array(mat->ctx, mat->row, isl_int *, n_row);
125
98.6k
  if (n_row && !row)
126
0
    goto error;
127
98.6k
  mat->row = row;
128
98.6k
129
1.03M
  for (i = 0; i < mat->n_row; 
++i932k
)
130
932k
    mat->row[i] = mat->block.data + (mat->row[i] - old);
131
343k
  for (i = mat->n_row; i < n_row; 
++i244k
)
132
244k
    mat->row[i] = mat->block.data + i * mat->max_col;
133
98.6k
  mat->n_row = n_row;
134
98.6k
  if (mat->n_col < n_col)
135
0
    mat->n_col = n_col;
136
98.6k
137
98.6k
  return mat;
138
0
error:
139
0
  isl_mat_free(mat);
140
0
  return NULL;
141
110k
}
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.82M
{
146
2.82M
  int i;
147
2.82M
  struct isl_mat *mat;
148
2.82M
149
2.82M
  mat = isl_alloc_type(ctx, struct isl_mat);
150
2.82M
  if (!mat)
151
0
    return NULL;
152
2.82M
  mat->row = isl_alloc_array(ctx, isl_int *, n_row);
153
2.82M
  if (n_row && 
!mat->row1.56M
)
154
0
    goto error;
155
9.78M
  
for (i = 0; 2.82M
i < n_row;
++i6.95M
)
156
6.95M
    mat->row[i] = row[first_row+i] + first_col;
157
2.82M
  mat->ctx = ctx;
158
2.82M
  isl_ctx_ref(ctx);
159
2.82M
  mat->ref = 1;
160
2.82M
  mat->n_row = n_row;
161
2.82M
  mat->n_col = n_col;
162
2.82M
  mat->block = isl_blk_empty();
163
2.82M
  mat->flags = ISL_MAT_BORROWED;
164
2.82M
  return mat;
165
0
error:
166
0
  free(mat);
167
0
  return NULL;
168
2.82M
}
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
622k
{
173
622k
  if (!mat)
174
0
    return NULL;
175
622k
  return isl_mat_sub_alloc6(mat->ctx, mat->row, first_row, n_row,
176
622k
          first_col, n_col);
177
622k
}
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
415k
{
182
415k
  int i;
183
415k
184
1.43M
  for (i = 0; i < n_row; 
++i1.02M
)
185
1.02M
    isl_seq_cpy(dst[i]+dst_col, src[i]+src_col, n_col);
186
415k
}
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
143k
{
191
143k
  int i;
192
143k
193
437k
  for (i = 0; i < n_row; 
++i293k
)
194
293k
    isl_seq_neg(dst[i]+dst_col, src[i]+src_col, n_col);
195
143k
}
196
197
__isl_give isl_mat *isl_mat_copy(__isl_keep isl_mat *mat)
198
2.14M
{
199
2.14M
  if (!mat)
200
0
    return NULL;
201
2.14M
202
2.14M
  mat->ref++;
203
2.14M
  return mat;
204
2.14M
}
205
206
__isl_give isl_mat *isl_mat_dup(__isl_keep isl_mat *mat)
207
364k
{
208
364k
  int i;
209
364k
  struct isl_mat *mat2;
210
364k
211
364k
  if (!mat)
212
2.96k
    return NULL;
213
361k
  mat2 = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
214
361k
  if (!mat2)
215
0
    return NULL;
216
1.16M
  
for (i = 0; 361k
i < mat->n_row;
++i807k
)
217
807k
    isl_seq_cpy(mat2->row[i], mat->row[i], mat->n_col);
218
364k
  return mat2;
219
364k
}
220
221
__isl_give isl_mat *isl_mat_cow(__isl_take isl_mat *mat)
222
3.96M
{
223
3.96M
  struct isl_mat *mat2;
224
3.96M
  if (!mat)
225
0
    return NULL;
226
3.96M
227
3.96M
  if (mat->ref == 1 && 
!3.92M
ISL_F_ISSET3.92M
(mat, ISL_MAT_BORROWED))
228
3.96M
    
return mat3.60M
;
229
358k
230
358k
  mat2 = isl_mat_dup(mat);
231
358k
  isl_mat_free(mat);
232
358k
  return mat2;
233
358k
}
234
235
__isl_null isl_mat *isl_mat_free(__isl_take isl_mat *mat)
236
12.2M
{
237
12.2M
  if (!mat)
238
1.35M
    return NULL;
239
10.8M
240
10.8M
  if (--mat->ref > 0)
241
2.14M
    return NULL;
242
8.72M
243
8.72M
  if (!ISL_F_ISSET(mat, ISL_MAT_BORROWED))
244
8.72M
    
isl_blk_free(mat->ctx, mat->block)5.89M
;
245
12.2M
  isl_ctx_deref(mat->ctx);
246
12.2M
  free(mat->row);
247
12.2M
  free(mat);
248
12.2M
249
12.2M
  return NULL;
250
12.2M
}
251
252
int isl_mat_rows(__isl_keep isl_mat *mat)
253
119k
{
254
119k
  return mat ? mat->n_row : 
-10
;
255
119k
}
256
257
int isl_mat_cols(__isl_keep isl_mat *mat)
258
44.1k
{
259
44.1k
  return mat ? mat->n_col : 
-10
;
260
44.1k
}
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.64k
{
266
4.64k
  if (!mat)
267
0
    return isl_stat_error;
268
4.64k
  if (col < 0 || col >= mat->n_col)
269
4.64k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
270
4.64k
      "column out of range", return isl_stat_error);
271
4.64k
  return isl_stat_ok;
272
4.64k
}
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.65k
{
278
4.65k
  if (!mat)
279
0
    return isl_stat_error;
280
4.65k
  if (row < 0 || row >= mat->n_row)
281
4.65k
    
isl_die0
(isl_mat_get_ctx(mat), isl_error_invalid,
282
4.65k
      "row out of range", return isl_stat_error);
283
4.65k
  return isl_stat_ok;
284
4.65k
}
285
286
int isl_mat_get_element(__isl_keep isl_mat *mat, int row, int col, isl_int *v)
287
3.03k
{
288
3.03k
  if (check_row(mat, row) < 0)
289
0
    return -1;
290
3.03k
  if (check_col(mat, col) < 0)
291
0
    return -1;
292
3.03k
  isl_int_set(*v, mat->row[row][col]);
293
3.03k
  return 0;
294
3.03k
}
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 (check_col(mat, col) < 0)
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.52k
{
314
1.52k
  mat = isl_mat_cow(mat);
315
1.52k
  if (check_row(mat, row) < 0)
316
0
    return isl_mat_free(mat);
317
1.52k
  if (check_col(mat, col) < 0)
318
0
    return isl_mat_free(mat);
319
1.52k
  isl_int_set(mat->row[row][col], v);
320
1.52k
  return mat;
321
1.52k
}
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 (check_col(mat, col) < 0)
330
0
    return isl_mat_free(mat);
331
43
  isl_int_set_si(mat->row[row][col], v);
332
43
  return mat;
333
43
}
334
335
/* Replace the element at row "row", column "col" of "mat" by "v".
336
 */
337
__isl_give isl_mat *isl_mat_set_element_val(__isl_take isl_mat *mat,
338
  int row, int col, __isl_take isl_val *v)
339
0
{
340
0
  if (!v)
341
0
    return isl_mat_free(mat);
342
0
  if (!isl_val_is_int(v))
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
1.00M
{
355
1.00M
  int i;
356
1.00M
  struct isl_mat *mat;
357
1.00M
358
1.00M
  mat = isl_mat_alloc(ctx, n_row, n_row);
359
1.00M
  if (!mat)
360
0
    return NULL;
361
5.39M
  
for (i = 0; 1.00M
i < n_row;
++i4.39M
) {
362
4.39M
    isl_seq_clr(mat->row[i], i);
363
4.39M
    isl_int_set(mat->row[i][i], d);
364
4.39M
    isl_seq_clr(mat->row[i]+i+1, n_row-(i+1));
365
4.39M
  }
366
1.00M
367
1.00M
  return mat;
368
1.00M
}
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_row;
++i41
)
381
151
    
isl_seq_clr(mat->row[i], n_col)41
;
382
151
383
151
  return mat;
384
151
}
385
386
__isl_give isl_mat *isl_mat_identity(isl_ctx *ctx, unsigned n_row)
387
1.00M
{
388
1.00M
  if (!ctx)
389
0
    return NULL;
390
1.00M
  return isl_mat_diag(ctx, n_row, ctx->one);
391
1.00M
}
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 (mat->n_row != mat->n_col)
402
0
    return 0;
403
0
404
0
  for (i = 0; i < mat->n_row; ++i) {
405
0
    if (isl_seq_first_non_zero(mat->row[i], i) != -1)
406
0
      return 0;
407
0
    if (isl_int_ne(mat->row[0][0], mat->row[i][i]))
408
0
      return 0;
409
0
    if (isl_seq_first_non_zero(mat->row[i] + i + 1,
410
0
              mat->n_col - (i + 1)) != -1)
411
0
      return 0;
412
0
  }
413
0
414
0
  return 1;
415
0
}
416
417
__isl_give isl_vec *isl_mat_vec_product(__isl_take isl_mat *mat,
418
  __isl_take isl_vec *vec)
419
343k
{
420
343k
  int i;
421
343k
  struct isl_vec *prod;
422
343k
423
343k
  if (!mat || !vec)
424
0
    goto error;
425
343k
426
343k
  isl_assert(mat->ctx, mat->n_col == vec->size, goto error);
427
343k
428
343k
  prod = isl_vec_alloc(mat->ctx, mat->n_row);
429
343k
  if (!prod)
430
0
    goto error;
431
343k
432
2.64M
  
for (i = 0; 343k
i < prod->size;
++i2.30M
)
433
2.30M
    isl_seq_inner_product(mat->row[i], vec->el, vec->size,
434
2.30M
          &prod->block.data[i]);
435
343k
  isl_mat_free(mat);
436
343k
  isl_vec_free(vec);
437
343k
  return prod;
438
0
error:
439
0
  isl_mat_free(mat);
440
0
  isl_vec_free(vec);
441
0
  return NULL;
442
343k
}
443
444
__isl_give isl_vec *isl_mat_vec_inverse_product(__isl_take isl_mat *mat,
445
  __isl_take isl_vec *vec)
446
160
{
447
160
  struct isl_mat *vec_mat;
448
160
  int i;
449
160
450
160
  if (!mat || !vec)
451
0
    goto error;
452
160
  vec_mat = isl_mat_alloc(vec->ctx, vec->size, 1);
453
160
  if (!vec_mat)
454
0
    goto error;
455
1.18k
  
for (i = 0; 160
i < vec->size;
++i1.02k
)
456
1.02k
    isl_int_set(vec_mat->row[i][0], vec->el[i]);
457
160
  vec_mat = isl_mat_inverse_product(mat, vec_mat);
458
160
  isl_vec_free(vec);
459
160
  if (!vec_mat)
460
0
    return NULL;
461
160
  vec = isl_vec_alloc(vec_mat->ctx, vec_mat->n_row);
462
160
  if (vec)
463
1.18k
    
for (i = 0; 160
i < vec->size;
++i1.02k
)
464
1.02k
      isl_int_set(vec->el[i], vec_mat->row[i][0]);
465
160
  isl_mat_free(vec_mat);
466
160
  return vec;
467
0
error:
468
0
  isl_mat_free(mat);
469
0
  isl_vec_free(vec);
470
0
  return NULL;
471
160
}
472
473
__isl_give isl_vec *isl_vec_mat_product(__isl_take isl_vec *vec,
474
  __isl_take isl_mat *mat)
475
18.7k
{
476
18.7k
  int i, j;
477
18.7k
  struct isl_vec *prod;
478
18.7k
479
18.7k
  if (!mat || !vec)
480
0
    goto error;
481
18.7k
482
18.7k
  isl_assert(mat->ctx, mat->n_row == vec->size, goto error);
483
18.7k
484
18.7k
  prod = isl_vec_alloc(mat->ctx, mat->n_col);
485
18.7k
  if (!prod)
486
0
    goto error;
487
18.7k
488
103k
  
for (i = 0; 18.7k
i < prod->size;
++i85.1k
) {
489
85.1k
    isl_int_set_si(prod->el[i], 0);
490
658k
    for (j = 0; j < vec->size; 
++j572k
)
491
572k
      isl_int_addmul(prod->el[i], vec->el[j], mat->row[j][i]);
492
85.1k
  }
493
18.7k
  isl_mat_free(mat);
494
18.7k
  isl_vec_free(vec);
495
18.7k
  return prod;
496
0
error:
497
0
  isl_mat_free(mat);
498
0
  isl_vec_free(vec);
499
0
  return NULL;
500
18.7k
}
501
502
__isl_give isl_mat *isl_mat_aff_direct_sum(__isl_take isl_mat *left,
503
  __isl_take isl_mat *right)
504
143k
{
505
143k
  int i;
506
143k
  struct isl_mat *sum;
507
143k
508
143k
  if (!left || !right)
509
0
    goto error;
510
143k
511
143k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
512
143k
  isl_assert(left->ctx, left->n_row >= 1, goto error);
513
143k
  isl_assert(left->ctx, left->n_col >= 1, goto error);
514
143k
  isl_assert(left->ctx, right->n_col >= 1, goto error);
515
143k
  isl_assert(left->ctx,
516
143k
      isl_seq_first_non_zero(left->row[0]+1, left->n_col-1) == -1,
517
143k
      goto error);
518
143k
  isl_assert(left->ctx,
519
143k
      isl_seq_first_non_zero(right->row[0]+1, right->n_col-1) == -1,
520
143k
      goto error);
521
143k
522
143k
  sum = isl_mat_alloc(left->ctx, left->n_row, left->n_col + right->n_col - 1);
523
143k
  if (!sum)
524
0
    goto error;
525
143k
  isl_int_lcm(sum->row[0][0], left->row[0][0], right->row[0][0]);
526
143k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
527
143k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
528
143k
529
143k
  isl_seq_clr(sum->row[0]+1, sum->n_col-1);
530
1.08M
  for (i = 1; i < sum->n_row; 
++i938k
) {
531
938k
    isl_int_mul(sum->row[i][0], left->row[0][0], left->row[i][0]);
532
938k
    isl_int_addmul(sum->row[i][0],
533
938k
        right->row[0][0], right->row[i][0]);
534
938k
    isl_seq_scale(sum->row[i]+1, left->row[i]+1, left->row[0][0],
535
938k
        left->n_col-1);
536
938k
    isl_seq_scale(sum->row[i]+left->n_col,
537
938k
        right->row[i]+1, right->row[0][0],
538
938k
        right->n_col-1);
539
938k
  }
540
143k
541
143k
  isl_int_divexact(left->row[0][0], sum->row[0][0], left->row[0][0]);
542
143k
  isl_int_divexact(right->row[0][0], sum->row[0][0], right->row[0][0]);
543
143k
  isl_mat_free(left);
544
143k
  isl_mat_free(right);
545
143k
  return sum;
546
0
error:
547
0
  isl_mat_free(left);
548
0
  isl_mat_free(right);
549
0
  return NULL;
550
143k
}
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
880k
{
555
880k
  int r;
556
6.56M
  for (r = row; r < M->n_row; 
++r5.67M
)
557
5.67M
    isl_int_swap(M->row[r][i], M->row[r][j]);
558
880k
  if (U) {
559
8.34M
    for (r = 0; r < (*U)->n_row; 
++r7.48M
)
560
7.48M
      isl_int_swap((*U)->row[r][i], (*U)->row[r][j]);
561
852k
  }
562
880k
  if (Q)
563
327k
    isl_mat_swap_rows(*Q, i, j);
564
880k
}
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
724k
{
569
724k
  int r;
570
3.02M
  for (r = row; r < M->n_row; 
++r2.29M
)
571
2.29M
    isl_int_submul(M->row[r][j], m, M->row[r][i]);
572
724k
  if (U) {
573
4.50M
    for (r = 0; r < (*U)->n_row; 
++r3.81M
)
574
3.81M
      isl_int_submul((*U)->row[r][j], m, (*U)->row[r][i]);
575
687k
  }
576
724k
  if (Q) {
577
1.32M
    for (r = 0; r < (*Q)->n_col; 
++r1.10M
)
578
1.10M
      isl_int_addmul((*Q)->row[i][r], m, (*Q)->row[j][r]);
579
215k
  }
580
724k
}
581
582
static void oppose(struct isl_mat *M, struct isl_mat **U,
583
  struct isl_mat **Q, unsigned row, unsigned col)
584
253k
{
585
253k
  int r;
586
1.08M
  for (r = row; r < M->n_row; 
++r830k
)
587
830k
    isl_int_neg(M->row[r][col], M->row[r][col]);
588
253k
  if (U) {
589
1.44M
    for (r = 0; r < (*U)->n_row; 
++r1.22M
)
590
1.22M
      isl_int_neg((*U)->row[r][col], (*U)->row[r][col]);
591
220k
  }
592
253k
  if (Q)
593
121k
    isl_seq_neg((*Q)->row[col], (*Q)->row[col], (*Q)->n_col);
594
253k
}
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
477k
{
611
477k
  isl_int c;
612
477k
  int row, col;
613
477k
614
477k
  if (U)
615
442k
    *U = NULL;
616
477k
  if (Q)
617
233k
    *Q = NULL;
618
477k
  if (!M)
619
0
    goto error;
620
477k
  M = isl_mat_cow(M);
621
477k
  if (!M)
622
0
    goto error;
623
477k
  if (U) {
624
442k
    *U = isl_mat_identity(M->ctx, M->n_col);
625
442k
    if (!*U)
626
0
      goto error;
627
477k
  }
628
477k
  if (Q) {
629
233k
    *Q = isl_mat_identity(M->ctx, M->n_col);
630
233k
    if (!*Q)
631
0
      goto error;
632
477k
  }
633
477k
634
477k
  col = 0;
635
477k
  isl_int_init(c);
636
2.22M
  for (row = 0; row < M->n_row; 
++row1.75M
) {
637
1.75M
    int first, i, off;
638
1.75M
    first = isl_seq_abs_min_non_zero(M->row[row]+col, M->n_col-col);
639
1.75M
    if (first == -1)
640
558k
      continue;
641
1.19M
    first += col;
642
1.19M
    if (first != col)
643
624k
      exchange(M, U, Q, row, first, col);
644
1.19M
    if (isl_int_is_neg(M->row[row][col]))
645
1.19M
      
oppose(M, U, Q, row, col)211k
;
646
1.19M
    first = col+1;
647
1.55M
    while ((off = isl_seq_first_non_zero(M->row[row]+first,
648
1.55M
                   M->n_col-first)) != -1) {
649
364k
      first += off;
650
364k
      isl_int_fdiv_q(c, M->row[row][first], M->row[row][col]);
651
364k
      subtract(M, U, Q, row, col, first, c);
652
364k
      if (!isl_int_is_zero(M->row[row][first]))
653
364k
        
exchange(M, U, Q, row, first, col)30.3k
;
654
333k
      else
655
333k
        ++first;
656
364k
    }
657
4.68M
    for (i = 0; i < col; 
++i3.49M
) {
658
3.49M
      if (isl_int_is_zero(M->row[row][i]))
659
3.49M
        
continue3.33M
;
660
152k
      if (neg)
661
152k
        
isl_int_cdiv_q0
(c, M->row[row][i], M->row[row][col]);
662
152k
      else
663
152k
        isl_int_fdiv_q(c, M->row[row][i], M->row[row][col]);
664
152k
      if (isl_int_is_zero(c))
665
152k
        
continue52.7k
;
666
100k
      subtract(M, U, Q, row, col, i, c);
667
100k
    }
668
1.75M
    ++col;
669
1.75M
  }
670
477k
  isl_int_clear(c);
671
477k
672
477k
  return M;
673
0
error:
674
0
  if (Q) {
675
0
    isl_mat_free(*Q);
676
0
    *Q = NULL;
677
0
  }
678
0
  if (U) {
679
0
    isl_mat_free(*U);
680
0
    *U = NULL;
681
0
  }
682
0
  isl_mat_free(M);
683
0
  return NULL;
684
477k
}
685
686
/* Use row "row" of "mat" to eliminate column "col" from all other rows.
687
 */
688
static __isl_give isl_mat *eliminate(__isl_take isl_mat *mat, int row, int col)
689
1.19k
{
690
1.19k
  int k, nr, nc;
691
1.19k
  isl_ctx *ctx;
692
1.19k
693
1.19k
  if (!mat)
694
0
    return NULL;
695
1.19k
696
1.19k
  ctx = isl_mat_get_ctx(mat);
697
1.19k
  nr = isl_mat_rows(mat);
698
1.19k
  nc = isl_mat_cols(mat);
699
1.19k
700
3.32k
  for (k = 0; k < nr; 
++k2.12k
) {
701
2.12k
    if (k == row)
702
1.19k
      continue;
703
922
    if (isl_int_is_zero(mat->row[k][col]))
704
922
      continue;
705
0
    mat = isl_mat_cow(mat);
706
0
    if (!mat)
707
0
      return NULL;
708
0
    isl_seq_elim(mat->row[k], mat->row[row], col, nc, NULL);
709
0
    isl_seq_normalize(ctx, mat->row[k], nc);
710
0
  }
711
1.19k
712
1.19k
  return mat;
713
1.19k
}
714
715
/* Perform Gaussian elimination on the rows of "mat", but start
716
 * from the final row and the final column.
717
 * Any zero rows that result from the elimination are removed.
718
 *
719
 * In particular, for each column from last to first,
720
 * look for the last row with a non-zero coefficient in that column,
721
 * move it last (but before other rows moved last in previous steps) and
722
 * use it to eliminate the column from the other rows.
723
 */
724
__isl_give isl_mat *isl_mat_reverse_gauss(__isl_take isl_mat *mat)
725
1.01k
{
726
1.01k
  int k, row, last, nr, nc;
727
1.01k
728
1.01k
  if (!mat)
729
0
    return NULL;
730
1.01k
731
1.01k
  nr = isl_mat_rows(mat);
732
1.01k
  nc = isl_mat_cols(mat);
733
1.01k
734
1.01k
  last = nc - 1;
735
2.21k
  for (row = nr - 1; row >= 0; 
--row1.19k
) {
736
1.58k
    for (; last >= 0; 
--last389
) {
737
2.05k
      for (k = row; k >= 0; 
--k469
)
738
1.66k
        if (!isl_int_is_zero(mat->row[k][last]))
739
1.66k
          
break1.19k
;
740
1.58k
      if (k >= 0)
741
1.19k
        break;
742
1.58k
    }
743
1.19k
    if (last < 0)
744
0
      break;
745
1.19k
    if (k != row)
746
0
      mat = isl_mat_swap_rows(mat, k, row);
747
1.19k
    if (!mat)
748
0
      return NULL;
749
1.19k
    if (isl_int_is_neg(mat->row[row][last]))
750
1.19k
      
mat = isl_mat_row_neg(mat, row)4
;
751
1.19k
    mat = eliminate(mat, row, last);
752
1.19k
    if (!mat)
753
0
      return NULL;
754
1.19k
  }
755
1.01k
  mat = isl_mat_drop_rows(mat, 0, row + 1);
756
1.01k
757
1.01k
  return mat;
758
1.01k
}
759
760
/* Negate the lexicographically negative rows of "mat" such that
761
 * all rows in the result are lexicographically non-negative.
762
 */
763
__isl_give isl_mat *isl_mat_lexnonneg_rows(__isl_take isl_mat *mat)
764
1.01k
{
765
1.01k
  int i, nr, nc;
766
1.01k
767
1.01k
  if (!mat)
768
0
    return NULL;
769
1.01k
770
1.01k
  nr = isl_mat_rows(mat);
771
1.01k
  nc = isl_mat_cols(mat);
772
1.01k
773
2.21k
  for (i = 0; i < nr; 
++i1.19k
) {
774
1.19k
    int pos;
775
1.19k
776
1.19k
    pos = isl_seq_first_non_zero(mat->row[i], nc);
777
1.19k
    if (pos < 0)
778
0
      continue;
779
1.19k
    if (isl_int_is_nonneg(mat->row[i][pos]))
780
1.19k
      
continue1.18k
;
781
14
    mat = isl_mat_row_neg(mat, i);
782
14
    if (!mat)
783
0
      return NULL;
784
1.19k
  }
785
1.01k
786
1.01k
  return mat;
787
1.01k
}
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 || !U)
797
0
    goto error;
798
0
799
0
  for (i = 0, rank = 0; rank < mat->n_col; ++rank) {
800
0
    while (i < mat->n_row && isl_int_is_zero(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
689k
{
820
689k
  int i;
821
689k
  struct isl_mat *mat2;
822
689k
823
689k
  if (!mat)
824
0
    return NULL;
825
689k
  mat2 = isl_mat_alloc(mat->ctx, 1+mat->n_row, 1+mat->n_col);
826
689k
  if (!mat2)
827
0
    goto error;
828
689k
  isl_int_set_si(mat2->row[0][0], 1);
829
689k
  isl_seq_clr(mat2->row[0]+1, mat->n_col);
830
4.20M
  for (i = 0; i < mat->n_row; 
++i3.51M
) {
831
3.51M
    isl_int_set_si(mat2->row[1+i][0], 0);
832
3.51M
    isl_seq_cpy(mat2->row[1+i]+1, mat->row[i], mat->n_col);
833
3.51M
  }
834
689k
  isl_mat_free(mat);
835
689k
  return mat2;
836
0
error:
837
0
  isl_mat_free(mat);
838
0
  return NULL;
839
689k
}
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
129k
{
849
129k
  int i;
850
129k
  isl_mat *mat;
851
129k
852
129k
  if (!mat1 || !mat2)
853
0
    goto error;
854
129k
855
129k
  mat = isl_mat_alloc(mat1->ctx, mat1->n_row + mat2->n_row,
856
129k
               mat1->n_col + mat2->n_col);
857
129k
  if (!mat)
858
0
    goto error;
859
259k
  
for (i = 0; 129k
i < mat1->n_row;
++i130k
) {
860
130k
    isl_seq_cpy(mat->row[i], mat1->row[i], mat1->n_col);
861
130k
    isl_seq_clr(mat->row[i] + mat1->n_col, mat2->n_col);
862
130k
  }
863
821k
  for (i = 0; i < mat2->n_row; 
++i691k
) {
864
691k
    isl_seq_clr(mat->row[mat1->n_row + i], mat1->n_col);
865
691k
    isl_seq_cpy(mat->row[mat1->n_row + i] + mat1->n_col,
866
691k
                mat2->row[i], mat2->n_col);
867
691k
  }
868
129k
  isl_mat_free(mat1);
869
129k
  isl_mat_free(mat2);
870
129k
  return mat;
871
0
error:
872
0
  isl_mat_free(mat1);
873
0
  isl_mat_free(mat2);
874
0
  return NULL;
875
129k
}
876
877
static int row_first_non_zero(isl_int **row, unsigned n_row, unsigned col)
878
949k
{
879
949k
  int i;
880
949k
881
1.60M
  for (i = 0; i < n_row; 
++i658k
)
882
1.13M
    if (!isl_int_is_zero(row[i][col]))
883
1.13M
      
return i476k
;
884
949k
  
return -1473k
;
885
949k
}
886
887
static int row_abs_min_non_zero(isl_int **row, unsigned n_row, unsigned col)
888
472k
{
889
472k
  int i, min = row_first_non_zero(row, n_row, col);
890
472k
  if (min < 0)
891
0
    return -1;
892
1.13M
  
for (i = min + 1; 472k
i < n_row;
++i659k
) {
893
659k
    if (isl_int_is_zero(row[i][col]))
894
659k
      
continue656k
;
895
3.56k
    if (isl_int_abs_lt(row[i][col], row[min][col]))
896
3.56k
      
min = i50
;
897
659k
  }
898
472k
  return min;
899
472k
}
900
901
static isl_stat inv_exchange(__isl_keep isl_mat **left,
902
  __isl_keep isl_mat **right, unsigned i, unsigned j)
903
341
{
904
341
  *left = isl_mat_swap_rows(*left, i, j);
905
341
  *right = isl_mat_swap_rows(*right, i, j);
906
341
907
341
  if (!*left || !*right)
908
0
    return isl_stat_error;
909
341
  return isl_stat_ok;
910
341
}
911
912
static void inv_oppose(
913
  struct isl_mat *left, struct isl_mat *right, unsigned row)
914
106
{
915
106
  isl_seq_neg(left->row[row]+row, left->row[row]+row, left->n_col-row);
916
106
  isl_seq_neg(right->row[row], right->row[row], right->n_col);
917
106
}
918
919
static void inv_subtract(struct isl_mat *left, struct isl_mat *right,
920
  unsigned row, unsigned i, isl_int m)
921
3.56k
{
922
3.56k
  isl_int_neg(m, m);
923
3.56k
  isl_seq_combine(left->row[i]+row,
924
3.56k
      left->ctx->one, left->row[i]+row,
925
3.56k
      m, left->row[row]+row,
926
3.56k
      left->n_col-row);
927
3.56k
  isl_seq_combine(right->row[i], right->ctx->one, right->row[i],
928
3.56k
      m, right->row[row], right->n_col);
929
3.56k
}
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
159k
{
936
159k
  int row;
937
159k
  isl_int a, b;
938
159k
939
159k
  if (!left || !right)
940
0
    goto error;
941
159k
942
159k
  isl_assert(left->ctx, left->n_row == left->n_col, goto error);
943
159k
  isl_assert(left->ctx, left->n_row == right->n_row, goto error);
944
159k
945
159k
  if (left->n_row == 0) {
946
0
    isl_mat_free(left);
947
0
    return right;
948
0
  }
949
159k
950
159k
  left = isl_mat_cow(left);
951
159k
  right = isl_mat_cow(right);
952
159k
  if (!left || !right)
953
0
    goto error;
954
159k
955
159k
  isl_int_init(a);
956
159k
  isl_int_init(b);
957
631k
  for (row = 0; row < left->n_row; 
++row472k
) {
958
472k
    int pivot, first, i, off;
959
472k
    pivot = row_abs_min_non_zero(left->row+row, left->n_row-row, row);
960
472k
    if (pivot < 0) {
961
0
      isl_int_clear(a);
962
0
      isl_int_clear(b);
963
0
      isl_assert(left->ctx, pivot >= 0, goto error);
964
0
    }
965
472k
    pivot += row;
966
472k
    if (pivot != row)
967
339
      if (inv_exchange(&left, &right, pivot, row) < 0)
968
0
        goto error;
969
472k
    if (isl_int_is_neg(left->row[row][row]))
970
472k
      
inv_oppose(left, right, row)106
;
971
472k
    first = row+1;
972
476k
    while ((off = row_first_non_zero(left->row+first,
973
476k
          left->n_row-first, row)) != -1) {
974
3.56k
      first += off;
975
3.56k
      isl_int_fdiv_q(a, left->row[first][row],
976
3.56k
          left->row[row][row]);
977
3.56k
      inv_subtract(left, right, row, first, a);
978
3.56k
      if (!isl_int_is_zero(left->row[first][row])) {
979
2
        if (inv_exchange(&left, &right, row, first) < 0)
980
0
          goto error;
981
3.56k
      } else {
982
3.56k
        ++first;
983
3.56k
      }
984
3.56k
    }
985
1.13M
    
for (i = 0; 472k
i < row;
++i660k
) {
986
660k
      if (isl_int_is_zero(left->row[i][row]))
987
660k
        
continue660k
;
988
347
      isl_int_gcd(a, left->row[row][row], left->row[i][row]);
989
347
      isl_int_divexact(b, left->row[i][row], a);
990
347
      isl_int_divexact(a, left->row[row][row], a);
991
347
      isl_int_neg(b, b);
992
660k
      isl_seq_combine(left->row[i] + i,
993
660k
          a, left->row[i] + i,
994
660k
          b, left->row[row] + i,
995
660k
          left->n_col - i);
996
660k
      isl_seq_combine(right->row[i], a, right->row[i],
997
660k
          b, right->row[row], right->n_col);
998
660k
    }
999
472k
  }
1000
159k
  isl_int_clear(b);
1001
159k
1002
159k
  isl_int_set(a, left->row[0][0]);
1003
472k
  for (row = 1; row < left->n_row; 
++row313k
)
1004
313k
    isl_int_lcm(a, a, left->row[row][row]);
1005
159k
  if (isl_int_is_zero(a)){
1006
0
    isl_int_clear(a);
1007
0
    isl_assert(left->ctx, 0, goto error);
1008
0
  }
1009
631k
  
for (row = 0; 159k
row < left->n_row;
++row472k
) {
1010
472k
    isl_int_divexact(left->row[row][row], a, left->row[row][row]);
1011
472k
    if (isl_int_is_one(left->row[row][row]))
1012
472k
      
continue465k
;
1013
6.96k
    isl_seq_scale(right->row[row], right->row[row],
1014
6.96k
        left->row[row][row], right->n_col);
1015
6.96k
  }
1016
159k
  isl_int_clear(a);
1017
159k
1018
159k
  isl_mat_free(left);
1019
159k
  return right;
1020
0
error:
1021
0
  isl_mat_free(left);
1022
0
  isl_mat_free(right);
1023
0
  return NULL;
1024
159k
}
1025
1026
void isl_mat_col_scale(struct isl_mat *mat, unsigned col, isl_int m)
1027
58.1k
{
1028
58.1k
  int i;
1029
58.1k
1030
283k
  for (i = 0; i < mat->n_row; 
++i225k
)
1031
225k
    isl_int_mul(mat->row[i][col], mat->row[i][col], m);
1032
58.1k
}
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
167k
{
1037
167k
  int i;
1038
167k
  isl_int tmp;
1039
167k
1040
167k
  isl_int_init(tmp);
1041
769k
  for (i = 0; i < mat->n_row; 
++i602k
) {
1042
602k
    isl_int_mul(tmp, m1, mat->row[i][src1]);
1043
602k
    isl_int_addmul(tmp, m2, mat->row[i][src2]);
1044
602k
    isl_int_set(mat->row[i][dst], tmp);
1045
602k
  }
1046
167k
  isl_int_clear(tmp);
1047
167k
}
1048
1049
__isl_give isl_mat *isl_mat_right_inverse(__isl_take isl_mat *mat)
1050
60.5k
{
1051
60.5k
  struct isl_mat *inv;
1052
60.5k
  int row;
1053
60.5k
  isl_int a, b;
1054
60.5k
1055
60.5k
  mat = isl_mat_cow(mat);
1056
60.5k
  if (!mat)
1057
0
    return NULL;
1058
60.5k
1059
60.5k
  inv = isl_mat_identity(mat->ctx, mat->n_col);
1060
60.5k
  inv = isl_mat_cow(inv);
1061
60.5k
  if (!inv)
1062
0
    goto error;
1063
60.5k
1064
60.5k
  isl_int_init(a);
1065
60.5k
  isl_int_init(b);
1066
245k
  for (row = 0; row < mat->n_row; 
++row185k
) {
1067
185k
    int pivot, first, i, off;
1068
185k
    pivot = isl_seq_abs_min_non_zero(mat->row[row]+row, mat->n_col-row);
1069
185k
    if (pivot < 0) {
1070
0
      isl_int_clear(a);
1071
0
      isl_int_clear(b);
1072
0
      isl_assert(mat->ctx, pivot >= 0, goto error);
1073
0
    }
1074
185k
    pivot += row;
1075
185k
    if (pivot != row)
1076
29.8k
      exchange(mat, &inv, NULL, row, pivot, row);
1077
185k
    if (isl_int_is_neg(mat->row[row][row]))
1078
185k
      
oppose(mat, &inv, NULL, row, row)41.9k
;
1079
185k
    first = row+1;
1080
445k
    while ((off = isl_seq_first_non_zero(mat->row[row]+first,
1081
445k
                mat->n_col-first)) != -1) {
1082
260k
      first += off;
1083
260k
      isl_int_fdiv_q(a, mat->row[row][first],
1084
260k
                mat->row[row][row]);
1085
260k
      subtract(mat, &inv, NULL, row, row, first, a);
1086
260k
      if (!isl_int_is_zero(mat->row[row][first]))
1087
260k
        
exchange(mat, &inv, NULL, row, row, first)196k
;
1088
64.0k
      else
1089
64.0k
        ++first;
1090
260k
    }
1091
458k
    for (i = 0; i < row; 
++i272k
) {
1092
272k
      if (isl_int_is_zero(mat->row[row][i]))
1093
272k
        
continue189k
;
1094
83.7k
      isl_int_gcd(a, mat->row[row][row], mat->row[row][i]);
1095
83.7k
      isl_int_divexact(b, mat->row[row][i], a);
1096
83.7k
      isl_int_divexact(a, mat->row[row][row], a);
1097
83.7k
      isl_int_neg(a, a);
1098
272k
      isl_mat_col_combine(mat, i, a, i, b, row);
1099
272k
      isl_mat_col_combine(inv, i, a, i, b, row);
1100
272k
    }
1101
185k
  }
1102
60.5k
  isl_int_clear(b);
1103
60.5k
1104
60.5k
  isl_int_set(a, mat->row[0][0]);
1105
185k
  for (row = 1; row < mat->n_row; 
++row124k
)
1106
124k
    isl_int_lcm(a, a, mat->row[row][row]);
1107
60.5k
  if (isl_int_is_zero(a)){
1108
0
    isl_int_clear(a);
1109
0
    goto error;
1110
0
  }
1111
245k
  
for (row = 0; 60.5k
row < mat->n_row;
++row185k
) {
1112
185k
    isl_int_divexact(mat->row[row][row], a, mat->row[row][row]);
1113
185k
    if (isl_int_is_one(mat->row[row][row]))
1114
185k
      
continue127k
;
1115
58.1k
    isl_mat_col_scale(inv, row, mat->row[row][row]);
1116
58.1k
  }
1117
60.5k
  isl_int_clear(a);
1118
60.5k
1119
60.5k
  isl_mat_free(mat);
1120
60.5k
1121
60.5k
  return inv;
1122
0
error:
1123
0
  isl_mat_free(mat);
1124
0
  isl_mat_free(inv);
1125
0
  return NULL;
1126
60.5k
}
1127
1128
__isl_give isl_mat *isl_mat_transpose(__isl_take isl_mat *mat)
1129
3.48k
{
1130
3.48k
  struct isl_mat *transpose = NULL;
1131
3.48k
  int i, j;
1132
3.48k
1133
3.48k
  if (!mat)
1134
0
    return NULL;
1135
3.48k
1136
3.48k
  if (mat->n_col == mat->n_row) {
1137
3.48k
    mat = isl_mat_cow(mat);
1138
3.48k
    if (!mat)
1139
0
      return NULL;
1140
18.7k
    
for (i = 0; 3.48k
i < mat->n_row;
++i15.2k
)
1141
54.7k
      
for (j = i + 1; 15.2k
j < mat->n_col;
++j39.5k
)
1142
39.5k
        isl_int_swap(mat->row[i][j], mat->row[j][i]);
1143
3.48k
    return mat;
1144
3.48k
  }
1145
0
  transpose = isl_mat_alloc(mat->ctx, mat->n_col, mat->n_row);
1146
0
  if (!transpose)
1147
0
    goto error;
1148
0
  for (i = 0; i < mat->n_row; ++i)
1149
0
    for (j = 0; j < mat->n_col; ++j)
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.48k
}
1157
1158
__isl_give isl_mat *isl_mat_swap_cols(__isl_take isl_mat *mat,
1159
  unsigned i, unsigned j)
1160
640k
{
1161
640k
  int r;
1162
640k
1163
640k
  mat = isl_mat_cow(mat);
1164
640k
  if (!mat)
1165
0
    return NULL;
1166
640k
  isl_assert(mat->ctx, i < mat->n_col, goto error);
1167
640k
  isl_assert(mat->ctx, j < mat->n_col, goto error);
1168
640k
1169
18.4M
  
for (r = 0; 640k
r < mat->n_row;
++r17.7M
)
1170
17.7M
    isl_int_swap(mat->row[r][i], mat->row[r][j]);
1171
640k
  return mat;
1172
0
error:
1173
0
  isl_mat_free(mat);
1174
0
  return NULL;
1175
640k
}
1176
1177
__isl_give isl_mat *isl_mat_swap_rows(__isl_take isl_mat *mat,
1178
  unsigned i, unsigned j)
1179
2.08M
{
1180
2.08M
  isl_int *t;
1181
2.08M
1182
2.08M
  if (!mat)
1183
0
    return NULL;
1184
2.08M
  mat = isl_mat_cow(mat);
1185
2.08M
  if (!mat)
1186
0
    return NULL;
1187
2.08M
  t = mat->row[i];
1188
2.08M
  mat->row[i] = mat->row[j];
1189
2.08M
  mat->row[j] = t;
1190
2.08M
  return mat;
1191
2.08M
}
1192
1193
/* Calculate the product of two matrices.
1194
 *
1195
 * This function is optimized for operand matrices that contain many zeros and
1196
 * skips multiplications where we know one of the operands is zero.
1197
 */
1198
__isl_give isl_mat *isl_mat_product(__isl_take isl_mat *left,
1199
  __isl_take isl_mat *right)
1200
2.03M
{
1201
2.03M
  int i, j, k;
1202
2.03M
  struct isl_mat *prod;
1203
2.03M
1204
2.03M
  if (!left || !right)
1205
0
    goto error;
1206
2.03M
  isl_assert(left->ctx, left->n_col == right->n_row, goto error);
1207
2.03M
  prod = isl_mat_alloc(left->ctx, left->n_row, right->n_col);
1208
2.03M
  if (!prod)
1209
0
    goto error;
1210
2.03M
  if (left->n_col == 0) {
1211
35
    for (i = 0; i < prod->n_row; 
++i8
)
1212
27
      
isl_seq_clr(prod->row[i], prod->n_col)8
;
1213
27
    isl_mat_free(left);
1214
27
    isl_mat_free(right);
1215
27
    return prod;
1216
27
  }
1217
6.84M
  
for (i = 0; 2.03M
i < prod->n_row;
++i4.80M
) {
1218
39.3M
    for (j = 0; j < prod->n_col; 
++j34.5M
)
1219
34.5M
      isl_int_mul(prod->row[i][j],
1220
4.80M
            left->row[i][0], right->row[0][j]);
1221
44.5M
    for (k = 1; k < left->n_col; 
++k39.7M
) {
1222
39.7M
      if (isl_int_is_zero(left->row[i][k]))
1223
39.7M
        
continue34.4M
;
1224
48.5M
      
for (j = 0; 5.30M
j < prod->n_col;
++j43.2M
)
1225
43.2M
        isl_int_addmul(prod->row[i][j],
1226
39.7M
              left->row[i][k], right->row[k][j]);
1227
39.7M
    }
1228
4.80M
  }
1229
2.03M
  isl_mat_free(left);
1230
2.03M
  isl_mat_free(right);
1231
2.03M
  return prod;
1232
0
error:
1233
0
  isl_mat_free(left);
1234
0
  isl_mat_free(right);
1235
0
  return NULL;
1236
2.03M
}
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.71M
{
1252
1.71M
  int i;
1253
1.71M
  struct isl_mat *t;
1254
1.71M
  int e;
1255
1.71M
1256
1.71M
  if (mat->n_col >= mat->n_row)
1257
1.07M
    e = 0;
1258
645k
  else
1259
645k
    e = mat->n_row - mat->n_col;
1260
1.71M
  if (has_div)
1261
573k
    
for (i = 0; 573k
i < n;
++i134
)
1262
573k
      
isl_int_mul134
(q[i][0], q[i][0], mat->row[0][0]);
1263
1.71M
  t = isl_mat_sub_alloc6(mat->ctx, q, 0, n, has_div, mat->n_row);
1264
1.71M
  t = isl_mat_product(t, mat);
1265
1.71M
  if (!t)
1266
0
    return -1;
1267
4.80M
  
for (i = 0; 1.71M
i < n;
++i3.08M
) {
1268
3.08M
    isl_seq_swp_or_cpy(q[i] + has_div, t->row[i], t->n_col);
1269
3.08M
    isl_seq_cpy(q[i] + has_div + t->n_col,
1270
3.08M
          q[i] + has_div + t->n_col + e, n_div);
1271
3.08M
    isl_seq_clr(q[i] + has_div + t->n_col + n_div, e);
1272
3.08M
  }
1273
1.71M
  isl_mat_free(t);
1274
1.71M
  return 0;
1275
1.71M
}
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
573k
{
1290
573k
  struct isl_ctx *ctx;
1291
573k
1292
573k
  if (!bset || !mat)
1293
0
    goto error;
1294
573k
1295
573k
  ctx = bset->ctx;
1296
573k
  bset = isl_basic_set_cow(bset);
1297
573k
  if (!bset)
1298
0
    goto error;
1299
573k
1300
573k
  isl_assert(ctx, bset->dim->nparam == 0, goto error);
1301
573k
  isl_assert(ctx, 1+bset->dim->n_out == mat->n_row, goto error);
1302
573k
  isl_assert(ctx, mat->n_col > 0, goto error);
1303
573k
1304
573k
  if (mat->n_col > mat->n_row) {
1305
103k
    bset = isl_basic_set_extend(bset, 0, mat->n_col-1, 0, 0, 0);
1306
103k
    if (!bset)
1307
0
      goto error;
1308
469k
  } else if (mat->n_col < mat->n_row) {
1309
215k
    bset->dim = isl_space_cow(bset->dim);
1310
215k
    if (!bset->dim)
1311
0
      goto error;
1312
215k
    bset->dim->n_out -= mat->n_row - mat->n_col;
1313
215k
  }
1314
573k
1315
573k
  if (preimage(ctx, bset->eq, bset->n_eq, bset->n_div, 0,
1316
573k
      isl_mat_copy(mat)) < 0)
1317
0
    goto error;
1318
573k
1319
573k
  if (preimage(ctx, bset->ineq, bset->n_ineq, bset->n_div, 0,
1320
573k
      isl_mat_copy(mat)) < 0)
1321
0
    goto error;
1322
573k
1323
573k
  if (preimage(ctx, bset->div, bset->n_div, bset->n_div, 1, mat) < 0)
1324
0
    goto error2;
1325
573k
1326
573k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_IMPLICIT);
1327
573k
  ISL_F_CLR(bset, ISL_BASIC_SET_NO_REDUNDANT);
1328
573k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED);
1329
573k
  ISL_F_CLR(bset, ISL_BASIC_SET_NORMALIZED_DIVS);
1330
573k
  ISL_F_CLR(bset, ISL_BASIC_SET_ALL_EQUALITIES);
1331
573k
1332
573k
  bset = isl_basic_set_simplify(bset);
1333
573k
  bset = isl_basic_set_finalize(bset);
1334
573k
1335
573k
  return bset;
1336
0
error:
1337
0
  isl_mat_free(mat);
1338
0
error2:
1339
0
  isl_basic_set_free(bset);
1340
0
  return NULL;
1341
573k
}
1342
1343
__isl_give isl_set *isl_set_preimage(
1344
  __isl_take isl_set *set, __isl_take isl_mat *mat)
1345
32.0k
{
1346
32.0k
  int i;
1347
32.0k
1348
32.0k
  set = isl_set_cow(set);
1349
32.0k
  if (!set)
1350
0
    goto error;
1351
32.0k
1352
91.0k
  
for (i = 0; 32.0k
i < set->n;
++i59.0k
) {
1353
59.0k
    set->p[i] = isl_basic_set_preimage(set->p[i],
1354
59.0k
                isl_mat_copy(mat));
1355
59.0k
    if (!set->p[i])
1356
0
      goto error;
1357
59.0k
  }
1358
32.0k
  if (mat->n_col != mat->n_row) {
1359
13.6k
    set->dim = isl_space_cow(set->dim);
1360
13.6k
    if (!set->dim)
1361
0
      goto error;
1362
13.6k
    set->dim->n_out += mat->n_col;
1363
13.6k
    set->dim->n_out -= mat->n_row;
1364
13.6k
  }
1365
32.0k
  isl_mat_free(mat);
1366
32.0k
  ISL_F_CLR(set, ISL_SET_NORMALIZED);
1367
32.0k
  return set;
1368
0
error:
1369
0
  isl_set_free(set);
1370
0
  isl_mat_free(mat);
1371
0
  return NULL;
1372
32.0k
}
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_row;
++i76
)
1393
144
    
isl_seq_swp_or_cpy(row[i] + first_col, t->row[i], t->n_col)76
;
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 (!mat) {
1403
0
    fprintf(out, "%*snull mat\n", indent, "");
1404
0
    return;
1405
0
  }
1406
0
1407
0
  if (mat->n_row == 0)
1408
0
    fprintf(out, "%*s[]\n", indent, "");
1409
0
1410
0
  for (i = 0; i < mat->n_row; ++i) {
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_col; ++j) {
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
87.3k
{
1435
87.3k
  int r;
1436
87.3k
1437
87.3k
  if (n == 0)
1438
2
    return mat;
1439
87.3k
1440
87.3k
  mat = isl_mat_cow(mat);
1441
87.3k
  if (!mat)
1442
0
    return NULL;
1443
87.3k
1444
87.3k
  if (col != mat->n_col-n) {
1445
66.3k
    for (r = 0; r < mat->n_row; 
++r51.8k
)
1446
51.8k
      isl_seq_cpy(mat->row[r]+col, mat->row[r]+col+n,
1447
51.8k
          mat->n_col - col - n);
1448
14.4k
  }
1449
87.3k
  mat->n_col -= n;
1450
87.3k
  return mat;
1451
87.3k
}
1452
1453
__isl_give isl_mat *isl_mat_drop_rows(__isl_take isl_mat *mat,
1454
  unsigned row, unsigned n)
1455
117k
{
1456
117k
  int r;
1457
117k
1458
117k
  mat = isl_mat_cow(mat);
1459
117k
  if (!mat)
1460
0
    return NULL;
1461
117k
1462
343k
  
for (r = row; 117k
r+n < mat->n_row;
++r225k
)
1463
225k
    mat->row[r] = mat->row[r+n];
1464
117k
1465
117k
  mat->n_row -= n;
1466
117k
  return mat;
1467
117k
}
1468
1469
__isl_give isl_mat *isl_mat_insert_cols(__isl_take isl_mat *mat,
1470
        unsigned col, unsigned n)
1471
82.8k
{
1472
82.8k
  isl_mat *ext;
1473
82.8k
1474
82.8k
  if (!mat)
1475
0
    return NULL;
1476
82.8k
  if (n == 0)
1477
435
    return mat;
1478
82.4k
1479
82.4k
  ext = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col + n);
1480
82.4k
  if (!ext)
1481
0
    goto error;
1482
82.4k
1483
82.4k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row, 0, 0, col);
1484
82.4k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, mat->n_row,
1485
82.4k
        col + n, col, mat->n_col - col);
1486
82.4k
1487
82.4k
  isl_mat_free(mat);
1488
82.4k
  return ext;
1489
0
error:
1490
0
  isl_mat_free(mat);
1491
0
  return NULL;
1492
82.8k
}
1493
1494
__isl_give isl_mat *isl_mat_insert_zero_cols(__isl_take isl_mat *mat,
1495
  unsigned first, unsigned n)
1496
82.8k
{
1497
82.8k
  int i;
1498
82.8k
1499
82.8k
  if (!mat)
1500
0
    return NULL;
1501
82.8k
  mat = isl_mat_insert_cols(mat, first, n);
1502
82.8k
  if (!mat)
1503
0
    return NULL;
1504
82.8k
1505
86.7k
  
for (i = 0; 82.8k
i < mat->n_row;
++i3.88k
)
1506
82.8k
    
isl_seq_clr(mat->row[i] + first, n)3.88k
;
1507
82.8k
1508
82.8k
  return mat;
1509
82.8k
}
1510
1511
__isl_give isl_mat *isl_mat_add_zero_cols(__isl_take isl_mat *mat, unsigned n)
1512
4.52k
{
1513
4.52k
  if (!mat)
1514
0
    return NULL;
1515
4.52k
1516
4.52k
  return isl_mat_insert_zero_cols(mat, mat->n_col, n);
1517
4.52k
}
1518
1519
__isl_give isl_mat *isl_mat_insert_rows(__isl_take isl_mat *mat,
1520
        unsigned row, unsigned n)
1521
4.97k
{
1522
4.97k
  isl_mat *ext;
1523
4.97k
1524
4.97k
  if (!mat)
1525
0
    return NULL;
1526
4.97k
  if (n == 0)
1527
433
    return mat;
1528
4.53k
1529
4.53k
  ext = isl_mat_alloc(mat->ctx, mat->n_row + n, mat->n_col);
1530
4.53k
  if (!ext)
1531
0
    goto error;
1532
4.53k
1533
4.53k
  isl_mat_sub_copy(mat->ctx, ext->row, mat->row, row, 0, 0, mat->n_col);
1534
4.53k
  isl_mat_sub_copy(mat->ctx, ext->row + row + n, mat->row + row,
1535
4.53k
        mat->n_row - row, 0, 0, mat->n_col);
1536
4.53k
1537
4.53k
  isl_mat_free(mat);
1538
4.53k
  return ext;
1539
0
error:
1540
0
  isl_mat_free(mat);
1541
0
  return NULL;
1542
4.97k
}
1543
1544
__isl_give isl_mat *isl_mat_add_rows(__isl_take isl_mat *mat, unsigned n)
1545
4.96k
{
1546
4.96k
  if (!mat)
1547
0
    return NULL;
1548
4.96k
1549
4.96k
  return isl_mat_insert_rows(mat, mat->n_row, n);
1550
4.96k
}
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 < n;
++i2
)
1562
2
    isl_seq_clr(mat->row[row + i], mat->n_col);
1563
2
1564
2
  return mat;
1565
2
}
1566
1567
__isl_give isl_mat *isl_mat_add_zero_rows(__isl_take isl_mat *mat, unsigned n)
1568
0
{
1569
0
  if (!mat)
1570
0
    return NULL;
1571
0
1572
0
  return isl_mat_insert_zero_rows(mat, mat->n_row, n);
1573
0
}
1574
1575
void isl_mat_col_submul(struct isl_mat *mat,
1576
      int dst_col, isl_int f, int src_col)
1577
5.15k
{
1578
5.15k
  int i;
1579
5.15k
1580
32.8k
  for (i = 0; i < mat->n_row; 
++i27.6k
)
1581
27.6k
    isl_int_submul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1582
5.15k
}
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_row;
++i5
)
1592
5
    isl_int_add(mat->row[i][dst_col],
1593
2
          mat->row[i][dst_col], mat->row[i][src_col]);
1594
2
}
1595
1596
void isl_mat_col_mul(struct isl_mat *mat, int dst_col, isl_int f, int src_col)
1597
10.9k
{
1598
10.9k
  int i;
1599
10.9k
1600
48.7k
  for (i = 0; i < mat->n_row; 
++i37.8k
)
1601
37.8k
    isl_int_mul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1602
10.9k
}
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 || check_col(mat, src_col) < 0)
1613
0
    return isl_mat_free(mat);
1614
6
1615
14
  
for (i = 0; 6
i < mat->n_row;
++i8
) {
1616
8
    if (isl_int_is_zero(mat->row[i][src_col]))
1617
8
      
continue2
;
1618
6
    mat = isl_mat_cow(mat);
1619
6
    if (!mat)
1620
0
      return NULL;
1621
6
    isl_int_addmul(mat->row[i][dst_col], f, mat->row[i][src_col]);
1622
6
  }
1623
6
1624
6
  return mat;
1625
6
}
1626
1627
/* Negate column "col" of "mat" and return the result.
1628
 */
1629
__isl_give isl_mat *isl_mat_col_neg(__isl_take isl_mat *mat, int col)
1630
4
{
1631
4
  int i;
1632
4
1633
4
  if (check_col(mat, col) < 0)
1634
0
    return isl_mat_free(mat);
1635
4
1636
10
  
for (i = 0; 4
i < mat->n_row;
++i6
) {
1637
6
    if (isl_int_is_zero(mat->row[i][col]))
1638
6
      
continue2
;
1639
4
    mat = isl_mat_cow(mat);
1640
4
    if (!mat)
1641
0
      return NULL;
1642
4
    isl_int_neg(mat->row[i][col], mat->row[i][col]);
1643
4
  }
1644
4
1645
4
  return mat;
1646
4
}
1647
1648
/* Negate row "row" of "mat" and return the result.
1649
 */
1650
__isl_give isl_mat *isl_mat_row_neg(__isl_take isl_mat *mat, int row)
1651
18
{
1652
18
  if (check_row(mat, row) < 0)
1653
0
    return isl_mat_free(mat);
1654
18
  if (isl_seq_first_non_zero(mat->row[row], mat->n_col) == -1)
1655
0
    return mat;
1656
18
  mat = isl_mat_cow(mat);
1657
18
  if (!mat)
1658
0
    return NULL;
1659
18
  isl_seq_neg(mat->row[row], mat->row[row], mat->n_col);
1660
18
  return mat;
1661
18
}
1662
1663
__isl_give isl_mat *isl_mat_unimodular_complete(__isl_take isl_mat *M, int row)
1664
17.8k
{
1665
17.8k
  int r;
1666
17.8k
  struct isl_mat *H = NULL, *Q = NULL;
1667
17.8k
1668
17.8k
  if (!M)
1669
0
    return NULL;
1670
17.8k
1671
17.8k
  isl_assert(M->ctx, M->n_row == M->n_col, goto error);
1672
17.8k
  M->n_row = row;
1673
17.8k
  H = isl_mat_left_hermite(isl_mat_copy(M), 0, NULL, &Q);
1674
17.8k
  M->n_row = M->n_col;
1675
17.8k
  if (!H)
1676
0
    goto error;
1677
35.6k
  
for (r = 0; 17.8k
r < row;
++r17.8k
)
1678
17.8k
    isl_assert(M->ctx, isl_int_is_one(H->row[r][r]), goto error);
1679
65.7k
  
for (r = row; 17.8k
r < M->n_row;
++r47.9k
)
1680
47.9k
    isl_seq_cpy(M->row[r], Q->row[r], M->n_col);
1681
17.8k
  isl_mat_free(H);
1682
17.8k
  isl_mat_free(Q);
1683
17.8k
  return M;
1684
0
error:
1685
0
  isl_mat_free(H);
1686
0
  isl_mat_free(Q);
1687
0
  isl_mat_free(M);
1688
0
  return NULL;
1689
17.8k
}
1690
1691
__isl_give isl_mat *isl_mat_concat(__isl_take isl_mat *top,
1692
  __isl_take isl_mat *bot)
1693
430
{
1694
430
  struct isl_mat *mat;
1695
430
1696
430
  if (!top || !bot)
1697
0
    goto error;
1698
430
1699
430
  isl_assert(top->ctx, top->n_col == bot->n_col, goto error);
1700
430
  if (top->n_row == 0) {
1701
390
    isl_mat_free(top);
1702
390
    return bot;
1703
390
  }
1704
40
  if (bot->n_row == 0) {
1705
0
    isl_mat_free(bot);
1706
0
    return top;
1707
0
  }
1708
40
1709
40
  mat = isl_mat_alloc(top->ctx, top->n_row + bot->n_row, top->n_col);
1710
40
  if (!mat)
1711
0
    goto error;
1712
40
  isl_mat_sub_copy(mat->ctx, mat->row, top->row, top->n_row,
1713
40
       0, 0, mat->n_col);
1714
40
  isl_mat_sub_copy(mat->ctx, mat->row + top->n_row, bot->row, bot->n_row,
1715
40
       0, 0, mat->n_col);
1716
40
  isl_mat_free(top);
1717
40
  isl_mat_free(bot);
1718
40
  return mat;
1719
0
error:
1720
0
  isl_mat_free(top);
1721
0
  isl_mat_free(bot);
1722
0
  return NULL;
1723
430
}
1724
1725
isl_bool isl_mat_is_equal(__isl_keep isl_mat *mat1, __isl_keep isl_mat *mat2)
1726
43.2k
{
1727
43.2k
  int i;
1728
43.2k
1729
43.2k
  if (!mat1 || !mat2)
1730
0
    return isl_bool_error;
1731
43.2k
1732
43.2k
  if (mat1->n_row != mat2->n_row)
1733
6.93k
    return isl_bool_false;
1734
36.3k
1735
36.3k
  if (mat1->n_col != mat2->n_col)
1736
0
    return isl_bool_false;
1737
36.3k
1738
70.3k
  
for (i = 0; 36.3k
i < mat1->n_row;
++i34.0k
)
1739
36.3k
    
if (34.1k
!isl_seq_eq(mat1->row[i], mat2->row[i], mat1->n_col)34.1k
)
1740
190
      return isl_bool_false;
1741
36.3k
1742
36.3k
  
return isl_bool_true36.1k
;
1743
43.2k
}
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 (row >= mat->n_row)
1773
24
    
isl_die0
(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
582
{
1793
582
  isl_mat *res;
1794
582
1795
582
  if (!mat)
1796
0
    return NULL;
1797
582
  if (n == 0 || dst_col == src_col)
1798
454
    return mat;
1799
128
1800
128
  res = isl_mat_alloc(mat->ctx, mat->n_row, mat->n_col);
1801
128
  if (!res)
1802
0
    goto error;
1803
128
1804
128
  if (dst_col < src_col) {
1805
128
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1806
128
         0, 0, dst_col);
1807
128
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1808
128
         dst_col, src_col, n);
1809
128
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1810
128
         dst_col + n, dst_col, src_col - dst_col);
1811
128
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1812
128
         src_col + n, src_col + n,
1813
128
         res->n_col - src_col - n);
1814
128
  } else {
1815
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1816
0
         0, 0, src_col);
1817
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1818
0
         src_col, src_col + n, dst_col - src_col);
1819
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1820
0
         dst_col, src_col, n);
1821
0
    isl_mat_sub_copy(res->ctx, res->row, mat->row, mat->n_row,
1822
0
         dst_col + n, dst_col + n,
1823
0
         res->n_col - dst_col - n);
1824
0
  }
1825
128
  isl_mat_free(mat);
1826
128
1827
128
  return res;
1828
0
error:
1829
0
  isl_mat_free(mat);
1830
0
  return NULL;
1831
582
}
1832
1833
/* Return the gcd of the elements in row "row" of "mat" in *gcd.
1834
 * Return isl_stat_ok on success and isl_stat_error on failure.
1835
 */
1836
isl_stat isl_mat_row_gcd(__isl_keep isl_mat *mat, int row, isl_int *gcd)
1837
7
{
1838
7
  if (check_row(mat, row) < 0)
1839
0
    return isl_stat_error;
1840
7
1841
7
  isl_seq_gcd(mat->row[row], mat->n_col, gcd);
1842
7
1843
7
  return isl_stat_ok;
1844
7
}
1845
1846
void isl_mat_gcd(__isl_keep isl_mat *mat, isl_int *gcd)
1847
771
{
1848
771
  int i;
1849
771
  isl_int g;
1850
771
1851
771
  isl_int_set_si(*gcd, 0);
1852
771
  if (!mat)
1853
0
    return;
1854
771
1855
771
  isl_int_init(g);
1856
4.72k
  for (i = 0; i < mat->n_row; 
++i3.95k
) {
1857
3.95k
    isl_seq_gcd(mat->row[i], mat->n_col, &g);
1858
3.95k
    isl_int_gcd(*gcd, *gcd, g);
1859
3.95k
  }
1860
771
  isl_int_clear(g);
1861
771
}
1862
1863
/* Return the result of scaling "mat" by a factor of "m".
1864
 */
1865
__isl_give isl_mat *isl_mat_scale(__isl_take isl_mat *mat, isl_int m)
1866
0
{
1867
0
  int i;
1868
0
1869
0
  if (isl_int_is_one(m))
1870
0
    return mat;
1871
0
1872
0
  mat = isl_mat_cow(mat);
1873
0
  if (!mat)
1874
0
    return NULL;
1875
0
1876
0
  for (i = 0; i < mat->n_row; ++i)
1877
0
    isl_seq_scale(mat->row[i], mat->row[i], m, mat->n_col);
1878
0
1879
0
  return mat;
1880
0
}
1881
1882
__isl_give isl_mat *isl_mat_scale_down(__isl_take isl_mat *mat, isl_int m)
1883
771
{
1884
771
  int i;
1885
771
1886
771
  if (isl_int_is_one(m))
1887
771
    
return mat0
;
1888
771
1889
771
  mat = isl_mat_cow(mat);
1890
771
  if (!mat)
1891
0
    return NULL;
1892
771
1893
4.72k
  
for (i = 0; 771
i < mat->n_row;
++i3.95k
)
1894
3.95k
    isl_seq_scale_down(mat->row[i], mat->row[i], m, mat->n_col);
1895
771
1896
771
  return mat;
1897
771
}
1898
1899
__isl_give isl_mat *isl_mat_scale_down_row(__isl_take isl_mat *mat, int row,
1900
  isl_int m)
1901
10
{
1902
10
  if (isl_int_is_one(m))
1903
10
    
return mat0
;
1904
10
1905
10
  mat = isl_mat_cow(mat);
1906
10
  if (!mat)
1907
0
    return NULL;
1908
10
1909
10
  isl_seq_scale_down(mat->row[row], mat->row[row], m, mat->n_col);
1910
10
1911
10
  return mat;
1912
10
}
1913
1914
__isl_give isl_mat *isl_mat_normalize(__isl_take isl_mat *mat)
1915
771
{
1916
771
  isl_int gcd;
1917
771
1918
771
  if (!mat)
1919
0
    return NULL;
1920
771
1921
771
  isl_int_init(gcd);
1922
771
  isl_mat_gcd(mat, &gcd);
1923
771
  mat = isl_mat_scale_down(mat, gcd);
1924
771
  isl_int_clear(gcd);
1925
771
1926
771
  return mat;
1927
771
}
1928
1929
__isl_give isl_mat *isl_mat_normalize_row(__isl_take isl_mat *mat, int row)
1930
1.64k
{
1931
1.64k
  mat = isl_mat_cow(mat);
1932
1.64k
  if (!mat)
1933
0
    return NULL;
1934
1.64k
1935
1.64k
  isl_seq_normalize(mat->ctx, mat->row[row], mat->n_col);
1936
1.64k
1937
1.64k
  return mat;
1938
1.64k
}
1939
1940
/* Number of initial non-zero columns.
1941
 */
1942
int isl_mat_initial_non_zero_cols(__isl_keep isl_mat *mat)
1943
1.01k
{
1944
1.01k
  int i;
1945
1.01k
1946
1.01k
  if (!mat)
1947
0
    return -1;
1948
1.01k
1949
1.47k
  
for (i = 0; 1.01k
i < mat->n_col;
++i456
)
1950
1.27k
    if (row_first_non_zero(mat->row, mat->n_row, i) < 0)
1951
818
      break;
1952
1.01k
1953
1.01k
  return i;
1954
1.01k
}