OpenCores
URL https://opencores.org/ocsvn/scarts/scarts/trunk

Subversion Repositories scarts

[/] [scarts/] [trunk/] [toolchain/] [scarts-gcc/] [gcc-4.1.1/] [gcc/] [lambda-mat.c] - Blame information for rev 12

Details | Compare with Previous | View Log

Line No. Rev Author Line
1 12 jlechner
/* Integer matrix math routines
2
   Copyright (C) 2003, 2004, 2005 Free Software Foundation, Inc.
3
   Contributed by Daniel Berlin <dberlin@dberlin.org>.
4
 
5
This file is part of GCC.
6
 
7
GCC is free software; you can redistribute it and/or modify it under
8
the terms of the GNU General Public License as published by the Free
9
Software Foundation; either version 2, or (at your option) any later
10
version.
11
 
12
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
13
WARRANTY; without even the implied warranty of MERCHANTABILITY or
14
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
15
for more details.
16
 
17
You should have received a copy of the GNU General Public License
18
along with GCC; see the file COPYING.  If not, write to the Free
19
Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
20
02110-1301, USA.  */
21
#include "config.h"
22
#include "system.h"
23
#include "coretypes.h"
24
#include "tm.h"
25
#include "ggc.h"
26
#include "varray.h"
27
#include "tree.h"
28
#include "lambda.h"
29
 
30
static void lambda_matrix_get_column (lambda_matrix, int, int,
31
                                      lambda_vector);
32
 
33
/* Allocate a matrix of M rows x  N cols.  */
34
 
35
lambda_matrix
36
lambda_matrix_new (int m, int n)
37
{
38
  lambda_matrix mat;
39
  int i;
40
 
41
  mat = ggc_alloc (m * sizeof (lambda_vector));
42
 
43
  for (i = 0; i < m; i++)
44
    mat[i] = lambda_vector_new (n);
45
 
46
  return mat;
47
}
48
 
49
/* Copy the elements of M x N matrix MAT1 to MAT2.  */
50
 
51
void
52
lambda_matrix_copy (lambda_matrix mat1, lambda_matrix mat2,
53
                    int m, int n)
54
{
55
  int i;
56
 
57
  for (i = 0; i < m; i++)
58
    lambda_vector_copy (mat1[i], mat2[i], n);
59
}
60
 
61
/* Store the N x N identity matrix in MAT.  */
62
 
63
void
64
lambda_matrix_id (lambda_matrix mat, int size)
65
{
66
  int i, j;
67
 
68
  for (i = 0; i < size; i++)
69
    for (j = 0; j < size; j++)
70
      mat[i][j] = (i == j) ? 1 : 0;
71
}
72
 
73
/* Return true if MAT is the identity matrix of SIZE */
74
 
75
bool
76
lambda_matrix_id_p (lambda_matrix mat, int size)
77
{
78
  int i, j;
79
  for (i = 0; i < size; i++)
80
    for (j = 0; j < size; j++)
81
      {
82
        if (i == j)
83
          {
84
            if (mat[i][j] != 1)
85
              return false;
86
          }
87
        else
88
          {
89
            if (mat[i][j] != 0)
90
              return false;
91
          }
92
      }
93
  return true;
94
}
95
 
96
/* Negate the elements of the M x N matrix MAT1 and store it in MAT2.  */
97
 
98
void
99
lambda_matrix_negate (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
100
{
101
  int i;
102
 
103
  for (i = 0; i < m; i++)
104
    lambda_vector_negate (mat1[i], mat2[i], n);
105
}
106
 
107
/* Take the transpose of matrix MAT1 and store it in MAT2.
108
   MAT1 is an M x N matrix, so MAT2 must be N x M.  */
109
 
110
void
111
lambda_matrix_transpose (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
112
{
113
  int i, j;
114
 
115
  for (i = 0; i < n; i++)
116
    for (j = 0; j < m; j++)
117
      mat2[i][j] = mat1[j][i];
118
}
119
 
120
 
121
/* Add two M x N matrices together: MAT3 = MAT1+MAT2.  */
122
 
123
void
124
lambda_matrix_add (lambda_matrix mat1, lambda_matrix mat2,
125
                   lambda_matrix mat3, int m, int n)
126
{
127
  int i;
128
 
129
  for (i = 0; i < m; i++)
130
    lambda_vector_add (mat1[i], mat2[i], mat3[i], n);
131
}
132
 
133
/* MAT3 = CONST1 * MAT1 + CONST2 * MAT2.  All matrices are M x N.  */
134
 
135
void
136
lambda_matrix_add_mc (lambda_matrix mat1, int const1,
137
                      lambda_matrix mat2, int const2,
138
                      lambda_matrix mat3, int m, int n)
139
{
140
  int i;
141
 
142
  for (i = 0; i < m; i++)
143
    lambda_vector_add_mc (mat1[i], const1, mat2[i], const2, mat3[i], n);
144
}
145
 
146
/* Multiply two matrices: MAT3 = MAT1 * MAT2.
147
   MAT1 is an M x R matrix, and MAT2 is R x N.  The resulting MAT2
148
   must therefore be M x N.  */
149
 
150
void
151
lambda_matrix_mult (lambda_matrix mat1, lambda_matrix mat2,
152
                    lambda_matrix mat3, int m, int r, int n)
153
{
154
 
155
  int i, j, k;
156
 
157
  for (i = 0; i < m; i++)
158
    {
159
      for (j = 0; j < n; j++)
160
        {
161
          mat3[i][j] = 0;
162
          for (k = 0; k < r; k++)
163
            mat3[i][j] += mat1[i][k] * mat2[k][j];
164
        }
165
    }
166
}
167
 
168
/* Get column COL from the matrix MAT and store it in VEC.  MAT has
169
   N rows, so the length of VEC must be N.  */
170
 
171
static void
172
lambda_matrix_get_column (lambda_matrix mat, int n, int col,
173
                          lambda_vector vec)
174
{
175
  int i;
176
 
177
  for (i = 0; i < n; i++)
178
    vec[i] = mat[i][col];
179
}
180
 
181
/* Delete rows r1 to r2 (not including r2).  */
182
 
183
void
184
lambda_matrix_delete_rows (lambda_matrix mat, int rows, int from, int to)
185
{
186
  int i;
187
  int dist;
188
  dist = to - from;
189
 
190
  for (i = to; i < rows; i++)
191
    mat[i - dist] = mat[i];
192
 
193
  for (i = rows - dist; i < rows; i++)
194
    mat[i] = NULL;
195
}
196
 
197
/* Swap rows R1 and R2 in matrix MAT.  */
198
 
199
void
200
lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2)
201
{
202
  lambda_vector row;
203
 
204
  row = mat[r1];
205
  mat[r1] = mat[r2];
206
  mat[r2] = row;
207
}
208
 
209
/* Add a multiple of row R1 of matrix MAT with N columns to row R2:
210
   R2 = R2 + CONST1 * R1.  */
211
 
212
void
213
lambda_matrix_row_add (lambda_matrix mat, int n, int r1, int r2, int const1)
214
{
215
  int i;
216
 
217
  if (const1 == 0)
218
    return;
219
 
220
  for (i = 0; i < n; i++)
221
    mat[r2][i] += const1 * mat[r1][i];
222
}
223
 
224
/* Negate row R1 of matrix MAT which has N columns.  */
225
 
226
void
227
lambda_matrix_row_negate (lambda_matrix mat, int n, int r1)
228
{
229
  lambda_vector_negate (mat[r1], mat[r1], n);
230
}
231
 
232
/* Multiply row R1 of matrix MAT with N columns by CONST1.  */
233
 
234
void
235
lambda_matrix_row_mc (lambda_matrix mat, int n, int r1, int const1)
236
{
237
  int i;
238
 
239
  for (i = 0; i < n; i++)
240
    mat[r1][i] *= const1;
241
}
242
 
243
/* Exchange COL1 and COL2 in matrix MAT. M is the number of rows.  */
244
 
245
void
246
lambda_matrix_col_exchange (lambda_matrix mat, int m, int col1, int col2)
247
{
248
  int i;
249
  int tmp;
250
  for (i = 0; i < m; i++)
251
    {
252
      tmp = mat[i][col1];
253
      mat[i][col1] = mat[i][col2];
254
      mat[i][col2] = tmp;
255
    }
256
}
257
 
258
/* Add a multiple of column C1 of matrix MAT with M rows to column C2:
259
   C2 = C2 + CONST1 * C1.  */
260
 
261
void
262
lambda_matrix_col_add (lambda_matrix mat, int m, int c1, int c2, int const1)
263
{
264
  int i;
265
 
266
  if (const1 == 0)
267
    return;
268
 
269
  for (i = 0; i < m; i++)
270
    mat[i][c2] += const1 * mat[i][c1];
271
}
272
 
273
/* Negate column C1 of matrix MAT which has M rows.  */
274
 
275
void
276
lambda_matrix_col_negate (lambda_matrix mat, int m, int c1)
277
{
278
  int i;
279
 
280
  for (i = 0; i < m; i++)
281
    mat[i][c1] *= -1;
282
}
283
 
284
/* Multiply column C1 of matrix MAT with M rows by CONST1.  */
285
 
286
void
287
lambda_matrix_col_mc (lambda_matrix mat, int m, int c1, int const1)
288
{
289
  int i;
290
 
291
  for (i = 0; i < m; i++)
292
    mat[i][c1] *= const1;
293
}
294
 
295
/* Compute the inverse of the N x N matrix MAT and store it in INV.
296
 
297
   We don't _really_ compute the inverse of MAT.  Instead we compute
298
   det(MAT)*inv(MAT), and we return det(MAT) to the caller as the function
299
   result.  This is necessary to preserve accuracy, because we are dealing
300
   with integer matrices here.
301
 
302
   The algorithm used here is a column based Gauss-Jordan elimination on MAT
303
   and the identity matrix in parallel.  The inverse is the result of applying
304
   the same operations on the identity matrix that reduce MAT to the identity
305
   matrix.
306
 
307
   When MAT is a 2 x 2 matrix, we don't go through the whole process, because
308
   it is easily inverted by inspection and it is a very common case.  */
309
 
310
static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int);
311
 
312
int
313
lambda_matrix_inverse (lambda_matrix mat, lambda_matrix inv, int n)
314
{
315
  if (n == 2)
316
    {
317
      int a, b, c, d, det;
318
      a = mat[0][0];
319
      b = mat[1][0];
320
      c = mat[0][1];
321
      d = mat[1][1];
322
      inv[0][0] =  d;
323
      inv[0][1] = -c;
324
      inv[1][0] = -b;
325
      inv[1][1] =  a;
326
      det = (a * d - b * c);
327
      if (det < 0)
328
        {
329
          det *= -1;
330
          inv[0][0] *= -1;
331
          inv[1][0] *= -1;
332
          inv[0][1] *= -1;
333
          inv[1][1] *= -1;
334
        }
335
      return det;
336
    }
337
  else
338
    return lambda_matrix_inverse_hard (mat, inv, n);
339
}
340
 
341
/* If MAT is not a special case, invert it the hard way.  */
342
 
343
static int
344
lambda_matrix_inverse_hard (lambda_matrix mat, lambda_matrix inv, int n)
345
{
346
  lambda_vector row;
347
  lambda_matrix temp;
348
  int i, j;
349
  int determinant;
350
 
351
  temp = lambda_matrix_new (n, n);
352
  lambda_matrix_copy (mat, temp, n, n);
353
  lambda_matrix_id (inv, n);
354
 
355
  /* Reduce TEMP to a lower triangular form, applying the same operations on
356
     INV which starts as the identity matrix.  N is the number of rows and
357
     columns.  */
358
  for (j = 0; j < n; j++)
359
    {
360
      row = temp[j];
361
 
362
      /* Make every element in the current row positive.  */
363
      for (i = j; i < n; i++)
364
        if (row[i] < 0)
365
          {
366
            lambda_matrix_col_negate (temp, n, i);
367
            lambda_matrix_col_negate (inv, n, i);
368
          }
369
 
370
      /* Sweep the upper triangle.  Stop when only the diagonal element in the
371
         current row is nonzero.  */
372
      while (lambda_vector_first_nz (row, n, j + 1) < n)
373
        {
374
          int min_col = lambda_vector_min_nz (row, n, j);
375
          lambda_matrix_col_exchange (temp, n, j, min_col);
376
          lambda_matrix_col_exchange (inv, n, j, min_col);
377
 
378
          for (i = j + 1; i < n; i++)
379
            {
380
              int factor;
381
 
382
              factor = -1 * row[i];
383
              if (row[j] != 1)
384
                factor /= row[j];
385
 
386
              lambda_matrix_col_add (temp, n, j, i, factor);
387
              lambda_matrix_col_add (inv, n, j, i, factor);
388
            }
389
        }
390
    }
391
 
392
  /* Reduce TEMP from a lower triangular to the identity matrix.  Also compute
393
     the determinant, which now is simply the product of the elements on the
394
     diagonal of TEMP.  If one of these elements is 0, the matrix has 0 as an
395
     eigenvalue so it is singular and hence not invertible.  */
396
  determinant = 1;
397
  for (j = n - 1; j >= 0; j--)
398
    {
399
      int diagonal;
400
 
401
      row = temp[j];
402
      diagonal = row[j];
403
 
404
      /* The matrix must not be singular.  */
405
      gcc_assert (diagonal);
406
 
407
      determinant = determinant * diagonal;
408
 
409
      /* If the diagonal is not 1, then multiply the each row by the
410
         diagonal so that the middle number is now 1, rather than a
411
         rational.  */
412
      if (diagonal != 1)
413
        {
414
          for (i = 0; i < j; i++)
415
            lambda_matrix_col_mc (inv, n, i, diagonal);
416
          for (i = j + 1; i < n; i++)
417
            lambda_matrix_col_mc (inv, n, i, diagonal);
418
 
419
          row[j] = diagonal = 1;
420
        }
421
 
422
      /* Sweep the lower triangle column wise.  */
423
      for (i = j - 1; i >= 0; i--)
424
        {
425
          if (row[i])
426
            {
427
              int factor = -row[i];
428
              lambda_matrix_col_add (temp, n, j, i, factor);
429
              lambda_matrix_col_add (inv, n, j, i, factor);
430
            }
431
 
432
        }
433
    }
434
 
435
  return determinant;
436
}
437
 
438
/* Decompose a N x N matrix MAT to a product of a lower triangular H
439
   and a unimodular U matrix such that MAT = H.U.  N is the size of
440
   the rows of MAT.  */
441
 
442
void
443
lambda_matrix_hermite (lambda_matrix mat, int n,
444
                       lambda_matrix H, lambda_matrix U)
445
{
446
  lambda_vector row;
447
  int i, j, factor, minimum_col;
448
 
449
  lambda_matrix_copy (mat, H, n, n);
450
  lambda_matrix_id (U, n);
451
 
452
  for (j = 0; j < n; j++)
453
    {
454
      row = H[j];
455
 
456
      /* Make every element of H[j][j..n] positive.  */
457
      for (i = j; i < n; i++)
458
        {
459
          if (row[i] < 0)
460
            {
461
              lambda_matrix_col_negate (H, n, i);
462
              lambda_vector_negate (U[i], U[i], n);
463
            }
464
        }
465
 
466
      /* Stop when only the diagonal element is nonzero.  */
467
      while (lambda_vector_first_nz (row, n, j + 1) < n)
468
        {
469
          minimum_col = lambda_vector_min_nz (row, n, j);
470
          lambda_matrix_col_exchange (H, n, j, minimum_col);
471
          lambda_matrix_row_exchange (U, j, minimum_col);
472
 
473
          for (i = j + 1; i < n; i++)
474
            {
475
              factor = row[i] / row[j];
476
              lambda_matrix_col_add (H, n, j, i, -1 * factor);
477
              lambda_matrix_row_add (U, n, i, j, factor);
478
            }
479
        }
480
    }
481
}
482
 
483
/* Given an M x N integer matrix A, this function determines an M x
484
   M unimodular matrix U, and an M x N echelon matrix S such that
485
   "U.A = S".  This decomposition is also known as "right Hermite".
486
 
487
   Ref: Algorithm 2.1 page 33 in "Loop Transformations for
488
   Restructuring Compilers" Utpal Banerjee.  */
489
 
490
void
491
lambda_matrix_right_hermite (lambda_matrix A, int m, int n,
492
                             lambda_matrix S, lambda_matrix U)
493
{
494
  int i, j, i0 = 0;
495
 
496
  lambda_matrix_copy (A, S, m, n);
497
  lambda_matrix_id (U, m);
498
 
499
  for (j = 0; j < n; j++)
500
    {
501
      if (lambda_vector_first_nz (S[j], m, i0) < m)
502
        {
503
          ++i0;
504
          for (i = m - 1; i >= i0; i--)
505
            {
506
              while (S[i][j] != 0)
507
                {
508
                  int sigma, factor, a, b;
509
 
510
                  a = S[i-1][j];
511
                  b = S[i][j];
512
                  sigma = (a * b < 0) ? -1: 1;
513
                  a = abs (a);
514
                  b = abs (b);
515
                  factor = sigma * (a / b);
516
 
517
                  lambda_matrix_row_add (S, n, i, i-1, -factor);
518
                  lambda_matrix_row_exchange (S, i, i-1);
519
 
520
                  lambda_matrix_row_add (U, m, i, i-1, -factor);
521
                  lambda_matrix_row_exchange (U, i, i-1);
522
                }
523
            }
524
        }
525
    }
526
}
527
 
528
/* Given an M x N integer matrix A, this function determines an M x M
529
   unimodular matrix V, and an M x N echelon matrix S such that "A =
530
   V.S".  This decomposition is also known as "left Hermite".
531
 
532
   Ref: Algorithm 2.2 page 36 in "Loop Transformations for
533
   Restructuring Compilers" Utpal Banerjee.  */
534
 
535
void
536
lambda_matrix_left_hermite (lambda_matrix A, int m, int n,
537
                             lambda_matrix S, lambda_matrix V)
538
{
539
  int i, j, i0 = 0;
540
 
541
  lambda_matrix_copy (A, S, m, n);
542
  lambda_matrix_id (V, m);
543
 
544
  for (j = 0; j < n; j++)
545
    {
546
      if (lambda_vector_first_nz (S[j], m, i0) < m)
547
        {
548
          ++i0;
549
          for (i = m - 1; i >= i0; i--)
550
            {
551
              while (S[i][j] != 0)
552
                {
553
                  int sigma, factor, a, b;
554
 
555
                  a = S[i-1][j];
556
                  b = S[i][j];
557
                  sigma = (a * b < 0) ? -1: 1;
558
                  a = abs (a);
559
      b = abs (b);
560
                  factor = sigma * (a / b);
561
 
562
                  lambda_matrix_row_add (S, n, i, i-1, -factor);
563
                  lambda_matrix_row_exchange (S, i, i-1);
564
 
565
                  lambda_matrix_col_add (V, m, i-1, i, factor);
566
                  lambda_matrix_col_exchange (V, m, i, i-1);
567
                }
568
            }
569
        }
570
    }
571
}
572
 
573
/* When it exists, return the first nonzero row in MAT after row
574
   STARTROW.  Otherwise return rowsize.  */
575
 
576
int
577
lambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize,
578
                            int startrow)
579
{
580
  int j;
581
  bool found = false;
582
 
583
  for (j = startrow; (j < rowsize) && !found; j++)
584
    {
585
      if ((mat[j] != NULL)
586
          && (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize))
587
        found = true;
588
    }
589
 
590
  if (found)
591
    return j - 1;
592
  return rowsize;
593
}
594
 
595
/* Calculate the projection of E sub k to the null space of B.  */
596
 
597
void
598
lambda_matrix_project_to_null (lambda_matrix B, int rowsize,
599
                               int colsize, int k, lambda_vector x)
600
{
601
  lambda_matrix M1, M2, M3, I;
602
  int determinant;
603
 
604
  /* Compute c(I-B^T inv(B B^T) B) e sub k.  */
605
 
606
  /* M1 is the transpose of B.  */
607
  M1 = lambda_matrix_new (colsize, colsize);
608
  lambda_matrix_transpose (B, M1, rowsize, colsize);
609
 
610
  /* M2 = B * B^T */
611
  M2 = lambda_matrix_new (colsize, colsize);
612
  lambda_matrix_mult (B, M1, M2, rowsize, colsize, rowsize);
613
 
614
  /* M3 = inv(M2) */
615
  M3 = lambda_matrix_new (colsize, colsize);
616
  determinant = lambda_matrix_inverse (M2, M3, rowsize);
617
 
618
  /* M2 = B^T (inv(B B^T)) */
619
  lambda_matrix_mult (M1, M3, M2, colsize, rowsize, rowsize);
620
 
621
  /* M1 = B^T (inv(B B^T)) B */
622
  lambda_matrix_mult (M2, B, M1, colsize, rowsize, colsize);
623
  lambda_matrix_negate (M1, M1, colsize, colsize);
624
 
625
  I = lambda_matrix_new (colsize, colsize);
626
  lambda_matrix_id (I, colsize);
627
 
628
  lambda_matrix_add_mc (I, determinant, M1, 1, M2, colsize, colsize);
629
 
630
  lambda_matrix_get_column (M2, colsize, k - 1, x);
631
 
632
}
633
 
634
/* Multiply a vector VEC by a matrix MAT.
635
   MAT is an M*N matrix, and VEC is a vector with length N.  The result
636
   is stored in DEST which must be a vector of length M.  */
637
 
638
void
639
lambda_matrix_vector_mult (lambda_matrix matrix, int m, int n,
640
                           lambda_vector vec, lambda_vector dest)
641
{
642
  int i, j;
643
 
644
  lambda_vector_clear (dest, m);
645
  for (i = 0; i < m; i++)
646
    for (j = 0; j < n; j++)
647
      dest[i] += matrix[i][j] * vec[j];
648
}
649
 
650
/* Print out an M x N matrix MAT to OUTFILE.  */
651
 
652
void
653
print_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n)
654
{
655
  int i;
656
 
657
  for (i = 0; i < m; i++)
658
    print_lambda_vector (outfile, matrix[i], n);
659
  fprintf (outfile, "\n");
660
}
661
 

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.