diff --git a/gcc/ChangeLog b/gcc/ChangeLog
index e8d1398d3f877262a7c3b436d7f2db0d19f0e2f8..fd7b7e348b540a44ac4c1bbe1e81d99e2180166d 100644
--- a/gcc/ChangeLog
+++ b/gcc/ChangeLog
@@ -1,3 +1,20 @@
+2004-08-10 Daniel Berlin <dberlin@dberlin.org>
+
+ * lambda.h: Add matrix type, and prototypes for remainder of
+ matrix and vector functions.
+ (lambda_vector_mult_const): New function.
+ (lambda_vector_negate): Ditto.
+ (lambda_vector_add): Ditto.
+ (lambda_vector_add_mc): Ditto.
+ (lambda_vector_copy): Ditto.
+ (lambda_vector_zerop): Ditto.
+ (lambda_vector_equal): Ditto.
+ (lambda_vector_min_nz): Ditto.
+ (lambda_vector_first_nz): Ditto.
+ (lambda_vector_matrix_mult): Ditto.
+ * lambda-mat.c: New file.
+ * Makefile.in (lambda-mat.o): New.
+
2004-08-10 Andrew MacLeod <amacleod@redhat.com>
* tree-cfg.c (bsi_insert_before, bsi_insert_after): Call modify_stmt
diff --git a/gcc/Makefile.in b/gcc/Makefile.in
index 984bb9e0bdd391638565b3ebd850962dbb80b3b5..c51760a0697e6b261de33168bf913bd7fe4ec3dd 100644
--- a/gcc/Makefile.in
+++ b/gcc/Makefile.in
@@ -140,7 +140,7 @@ XCFLAGS =
TCFLAGS =
CFLAGS = -g
STAGE1_CFLAGS = -g @stage1_cflags@
-BOOT_CFLAGS = -g -O2
+BOOT_CFLAGS = -g -O2 -dU
# Flags to determine code coverage. When coverage is disabled, this will
# contain the optimization flags, as you normally want code coverage
@@ -922,7 +922,7 @@ OBJS-common = \
targhooks.o timevar.o toplev.o tracer.o tree.o tree-dump.o unroll.o \
varasm.o varray.o vec.o version.o vmsdbgout.o xcoffout.o alloc-pool.o \
et-forest.o cfghooks.o bt-load.o pretty-print.o $(GGC) web.o passes.o \
- rtl-profile.o tree-profile.o rtlhooks.o cfgexpand.o
+ rtl-profile.o tree-profile.o rtlhooks.o cfgexpand.o lambda-mat.o
OBJS-md = $(out_object_file)
OBJS-archive = $(EXTRA_OBJS) $(host_hook_obj) tree-inline.o \
@@ -2127,6 +2127,7 @@ ifcvt.o : ifcvt.c $(CONFIG_H) $(SYSTEM_H) coretypes.h $(TM_H) $(RTL_H) \
$(REGS_H) toplev.h $(FLAGS_H) insn-config.h function.h $(RECOG_H) $(TARGET_H) \
$(BASIC_BLOCK_H) $(EXPR_H) output.h except.h $(TM_P_H) real.h $(OPTABS_H) \
$(CFGLOOP_H)
+lambda-mat.o : lambda-mat.c lambda.h $(GGC_H) $(SYSTEM_H) $(CONFIG_H) $(TM_H)
params.o : params.c $(CONFIG_H) $(SYSTEM_H) coretypes.h $(TM_H) $(PARAMS_H) toplev.h
hooks.o: hooks.c $(CONFIG_H) $(SYSTEM_H) coretypes.h $(TM_H) $(HOOKS_H)
pretty-print.o: $(CONFIG_H) $(SYSTEM_H) pretty-print.c $(PRETTY_PRINT_H)
diff --git a/gcc/lambda-mat.c b/gcc/lambda-mat.c
new file mode 100644
index 0000000000000000000000000000000000000000..987fa5e1fd17c46b13c867ffc1d4c5522a6459df
--- /dev/null
+++ b/gcc/lambda-mat.c
@@ -0,0 +1,638 @@
+/* Integer matrix math routines
+ Copyright (C) 2003, 2004 Free Software Foundation, Inc.
+ Contributed by Daniel Berlin <dberlin@dberlin.org>.
+
+This file is part of GCC.
+
+GCC is free software; you can redistribute it and/or modify it under
+the terms of the GNU General Public License as published by the Free
+Software Foundation; either version 2, or (at your option) any later
+version.
+
+GCC is distributed in the hope that it will be useful, but WITHOUT ANY
+WARRANTY; without even the implied warranty of MERCHANTABILITY or
+FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+for more details.
+
+You should have received a copy of the GNU General Public License
+along with GCC; see the file COPYING. If not, write to the Free
+Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+02111-1307, USA. */
+#include "config.h"
+#include "system.h"
+#include "coretypes.h"
+#include "tm.h"
+#include "ggc.h"
+#include "varray.h"
+#include "lambda.h"
+
+static void lambda_matrix_get_column (lambda_matrix, int, int,
+ lambda_vector);
+
+/* Allocate a matrix of M rows x N cols. */
+
+lambda_matrix
+lambda_matrix_new (int m, int n)
+{
+ lambda_matrix mat;
+ int i;
+
+ mat = ggc_alloc (m * sizeof (lambda_vector));
+
+ for (i = 0; i < m; i++)
+ mat[i] = lambda_vector_new (n);
+
+ return mat;
+}
+
+/* Copy the elements of M x N matrix MAT1 to MAT2. */
+
+void
+lambda_matrix_copy (lambda_matrix mat1, lambda_matrix mat2,
+ int m, int n)
+{
+ int i;
+
+ for (i = 0; i < m; i++)
+ lambda_vector_copy (mat1[i], mat2[i], n);
+}
+
+/* Store the N x N identity matrix in MAT. */
+
+void
+lambda_matrix_id (lambda_matrix mat, int size)
+{
+ int i, j;
+
+ for (i = 0; i < size; i++)
+ for (j = 0; j < size; j++)
+ mat[i][j] = (i == j) ? 1 : 0;
+}
+
+/* Negate the elements of the M x N matrix MAT1 and store it in MAT2. */
+
+void
+lambda_matrix_negate (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
+{
+ int i;
+
+ for (i = 0; i < m; i++)
+ lambda_vector_negate (mat1[i], mat2[i], n);
+}
+
+/* Take the transpose of matrix MAT1 and store it in MAT2.
+ MAT1 is an M x N matrix, so MAT2 must be N x M. */
+
+void
+lambda_matrix_transpose (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
+{
+ int i, j;
+
+ for (i = 0; i < n; i++)
+ for (j = 0; j < m; j++)
+ mat2[i][j] = mat1[j][i];
+}
+
+
+/* Add two M x N matrices together: MAT3 = MAT1+MAT2. */
+
+void
+lambda_matrix_add (lambda_matrix mat1, lambda_matrix mat2,
+ lambda_matrix mat3, int m, int n)
+{
+ int i;
+
+ for (i = 0; i < m; i++)
+ lambda_vector_add (mat1[i], mat2[i], mat3[i], n);
+}
+
+/* MAT3 = CONST1 * MAT1 + CONST2 * MAT2. All matrices are M x N. */
+
+void
+lambda_matrix_add_mc (lambda_matrix mat1, int const1,
+ lambda_matrix mat2, int const2,
+ lambda_matrix mat3, int m, int n)
+{
+ int i;
+
+ for (i = 0; i < m; i++)
+ lambda_vector_add_mc (mat1[i], const1, mat2[i], const2, mat3[i], n);
+}
+
+/* Multiply two matrices: MAT3 = MAT1 * MAT2.
+ MAT1 is an M x R matrix, and MAT2 is R x N. The resulting MAT2
+ must therefore be M x N. */
+
+void
+lambda_matrix_mult (lambda_matrix mat1, lambda_matrix mat2,
+ lambda_matrix mat3, int m, int r, int n)
+{
+
+ int i, j, k;
+
+ for (i = 0; i < m; i++)
+ {
+ for (j = 0; j < n; j++)
+ {
+ mat3[i][j] = 0;
+ for (k = 0; k < r; k++)
+ mat3[i][j] += mat1[i][k] * mat2[k][j];
+ }
+ }
+}
+
+/* Get column COL from the matrix MAT and store it in VEC. MAT has
+ N rows, so the length of VEC must be N. */
+
+static void
+lambda_matrix_get_column (lambda_matrix mat, int n, int col,
+ lambda_vector vec)
+{
+ int i;
+
+ for (i = 0; i < n; i++)
+ vec[i] = mat[i][col];
+}
+
+/* Delete rows r1 to r2 (not including r2). */
+
+void
+lambda_matrix_delete_rows (lambda_matrix mat, int rows, int from, int to)
+{
+ int i;
+ int dist;
+ dist = to - from;
+
+ for (i = to; i < rows; i++)
+ mat[i - dist] = mat[i];
+
+ for (i = rows - dist; i < rows; i++)
+ mat[i] = NULL;
+}
+
+/* Swap rows R1 and R2 in matrix MAT. */
+
+void
+lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2)
+{
+ lambda_vector row;
+
+ row = mat[r1];
+ mat[r1] = mat[r2];
+ mat[r2] = row;
+}
+
+/* Add a multiple of row R1 of matrix MAT with N columns to row R2:
+ R2 = R2 + CONST1 * R1. */
+
+void
+lambda_matrix_row_add (lambda_matrix mat, int n, int r1, int r2, int const1)
+{
+ int i;
+
+ if (const1 == 0)
+ return;
+
+ for (i = 0; i < n; i++)
+ mat[r2][i] += const1 * mat[r1][i];
+}
+
+/* Negate row R1 of matrix MAT which has N columns. */
+
+void
+lambda_matrix_row_negate (lambda_matrix mat, int n, int r1)
+{
+ lambda_vector_negate (mat[r1], mat[r1], n);
+}
+
+/* Multiply row R1 of matrix MAT with N columns by CONST1. */
+
+void
+lambda_matrix_row_mc (lambda_matrix mat, int n, int r1, int const1)
+{
+ int i;
+
+ for (i = 0; i < n; i++)
+ mat[r1][i] *= const1;
+}
+
+/* Exchange COL1 and COL2 in matrix MAT. M is the number of rows. */
+
+void
+lambda_matrix_col_exchange (lambda_matrix mat, int m, int col1, int col2)
+{
+ int i;
+ int tmp;
+ for (i = 0; i < m; i++)
+ {
+ tmp = mat[i][col1];
+ mat[i][col1] = mat[i][col2];
+ mat[i][col2] = tmp;
+ }
+}
+
+/* Add a multiple of column C1 of matrix MAT with M rows to column C2:
+ C2 = C2 + CONST1 * C1. */
+
+void
+lambda_matrix_col_add (lambda_matrix mat, int m, int c1, int c2, int const1)
+{
+ int i;
+
+ if (const1 == 0)
+ return;
+
+ for (i = 0; i < m; i++)
+ mat[i][c2] += const1 * mat[i][c1];
+}
+
+/* Negate column C1 of matrix MAT which has M rows. */
+
+void
+lambda_matrix_col_negate (lambda_matrix mat, int m, int c1)
+{
+ int i;
+
+ for (i = 0; i < m; i++)
+ mat[i][c1] *= -1;
+}
+
+/* Multiply column C1 of matrix MAT with M rows by CONST1. */
+
+void
+lambda_matrix_col_mc (lambda_matrix mat, int m, int c1, int const1)
+{
+ int i;
+
+ for (i = 0; i < m; i++)
+ mat[i][c1] *= const1;
+}
+
+/* Compute the inverse of the N x N matrix MAT and store it in INV.
+
+ We don't _really_ compute the inverse of MAT. Instead we compute
+ det(MAT)*inv(MAT), and we return det(MAT) to the caller as the function
+ result. This is necessary to preserve accuracy, because we are dealing
+ with integer matrices here.
+
+ The algorithm used here is a column based Gauss-Jordan elimination on MAT
+ and the identity matrix in parallel. The inverse is the result of applying
+ the same operations on the identity matrix that reduce MAT to the identity
+ matrix.
+
+ When MAT is a 2 x 2 matrix, we don't go through the whole process, because
+ it is easily inverted by inspection and it is a very common case. */
+
+static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int);
+
+int
+lambda_matrix_inverse (lambda_matrix mat, lambda_matrix inv, int n)
+{
+ if (n == 2)
+ {
+ int a, b, c, d, det;
+ a = mat[0][0];
+ b = mat[1][0];
+ c = mat[0][1];
+ d = mat[1][1];
+ inv[0][0] = d;
+ inv[0][1] = -c;
+ inv[1][0] = -b;
+ inv[1][1] = a;
+ det = (a * d - b * c);
+ if (det < 0)
+ {
+ det *= -1;
+ inv[0][0] *= -1;
+ inv[1][0] *= -1;
+ inv[0][1] *= -1;
+ inv[1][1] *= -1;
+ }
+ return det;
+ }
+ else
+ return lambda_matrix_inverse_hard (mat, inv, n);
+}
+
+/* If MAT is not a special case, invert it the hard way. */
+
+static int
+lambda_matrix_inverse_hard (lambda_matrix mat, lambda_matrix inv, int n)
+{
+ lambda_vector row;
+ lambda_matrix temp;
+ int i, j;
+ int determinant;
+
+ temp = lambda_matrix_new (n, n);
+ lambda_matrix_copy (mat, temp, n, n);
+ lambda_matrix_id (inv, n);
+
+ /* Reduce TEMP to a lower triangular form, applying the same operations on
+ INV which starts as the identity matrix. N is the number of rows and
+ columns. */
+ for (j = 0; j < n; j++)
+ {
+ row = temp[j];
+
+ /* Make every element in the current row positive. */
+ for (i = j; i < n; i++)
+ if (row[i] < 0)
+ {
+ lambda_matrix_col_negate (temp, n, i);
+ lambda_matrix_col_negate (inv, n, i);
+ }
+
+ /* Sweep the upper triangle. Stop when only the diagonal element in the
+ current row is nonzero. */
+ while (lambda_vector_first_nz (row, n, j + 1) < n)
+ {
+ int min_col = lambda_vector_min_nz (row, n, j);
+ lambda_matrix_col_exchange (temp, n, j, min_col);
+ lambda_matrix_col_exchange (inv, n, j, min_col);
+
+ for (i = j + 1; i < n; i++)
+ {
+ int factor;
+
+ factor = -1 * row[i];
+ if (row[j] != 1)
+ factor /= row[j];
+
+ lambda_matrix_col_add (temp, n, j, i, factor);
+ lambda_matrix_col_add (inv, n, j, i, factor);
+ }
+ }
+ }
+
+ /* Reduce TEMP from a lower triangular to the identity matrix. Also compute
+ the determinant, which now is simply the product of the elements on the
+ diagonal of TEMP. If one of these elements is 0, the matrix has 0 as an
+ eigenvalue so it is singular and hence not invertible. */
+ determinant = 1;
+ for (j = n - 1; j >= 0; j--)
+ {
+ int diagonal;
+
+ row = temp[j];
+ diagonal = row[j];
+
+ /* If the matrix is singular, abort. */
+ if (diagonal == 0)
+ abort ();
+
+ determinant = determinant * diagonal;
+
+ /* If the diagonal is not 1, then multiply the each row by the
+ diagonal so that the middle number is now 1, rather than a
+ rational. */
+ if (diagonal != 1)
+ {
+ for (i = 0; i < j; i++)
+ lambda_matrix_col_mc (inv, n, i, diagonal);
+ for (i = j + 1; i < n; i++)
+ lambda_matrix_col_mc (inv, n, i, diagonal);
+
+ row[j] = diagonal = 1;
+ }
+
+ /* Sweep the lower triangle column wise. */
+ for (i = j - 1; i >= 0; i--)
+ {
+ if (row[i])
+ {
+ int factor = -row[i];
+ lambda_matrix_col_add (temp, n, j, i, factor);
+ lambda_matrix_col_add (inv, n, j, i, factor);
+ }
+
+ }
+ }
+
+ return determinant;
+}
+
+/* Decompose a N x N matrix MAT to a product of a lower triangular H
+ and a unimodular U matrix such that MAT = H.U. N is the size of
+ the rows of MAT. */
+
+void
+lambda_matrix_hermite (lambda_matrix mat, int n,
+ lambda_matrix H, lambda_matrix U)
+{
+ lambda_vector row;
+ int i, j, factor, minimum_col;
+
+ lambda_matrix_copy (mat, H, n, n);
+ lambda_matrix_id (U, n);
+
+ for (j = 0; j < n; j++)
+ {
+ row = H[j];
+
+ /* Make every element of H[j][j..n] positive. */
+ for (i = j; i < n; i++)
+ {
+ if (row[i] < 0)
+ {
+ lambda_matrix_col_negate (H, n, i);
+ lambda_vector_negate (U[i], U[i], n);
+ }
+ }
+
+ /* Stop when only the diagonal element is non-zero. */
+ while (lambda_vector_first_nz (row, n, j + 1) < n)
+ {
+ minimum_col = lambda_vector_min_nz (row, n, j);
+ lambda_matrix_col_exchange (H, n, j, minimum_col);
+ lambda_matrix_row_exchange (U, j, minimum_col);
+
+ for (i = j + 1; i < n; i++)
+ {
+ factor = row[i] / row[j];
+ lambda_matrix_col_add (H, n, j, i, -1 * factor);
+ lambda_matrix_row_add (U, n, i, j, factor);
+ }
+ }
+ }
+}
+
+/* Given an M x N integer matrix A, this function determines an M x
+ M unimodular matrix U, and an M x N echelon matrix S such that
+ "U.A = S". This decomposition is also known as "right Hermite".
+
+ Ref: Algorithm 2.1 page 33 in "Loop Transformations for
+ Restructuring Compilers" Utpal Banerjee. */
+
+void
+lambda_matrix_right_hermite (lambda_matrix A, int m, int n,
+ lambda_matrix S, lambda_matrix U)
+{
+ int i, j, i0 = 0;
+
+ lambda_matrix_copy (A, S, m, n);
+ lambda_matrix_id (U, m);
+
+ for (j = 0; j < n; j++)
+ {
+ if (lambda_vector_first_nz (S[j], m, i0) < m)
+ {
+ ++i0;
+ for (i = m - 1; i >= i0; i--)
+ {
+ while (S[i][j] != 0)
+ {
+ int sigma, factor, a, b;
+
+ a = S[i-1][j];
+ b = S[i][j];
+ sigma = (a * b < 0) ? -1: 1;
+ a = abs (a);
+ b = abs (b);
+ factor = sigma * (a / b);
+
+ lambda_matrix_row_add (S, n, i, i-1, -factor);
+ lambda_matrix_row_exchange (S, i, i-1);
+
+ lambda_matrix_row_add (U, m, i, i-1, -factor);
+ lambda_matrix_row_exchange (U, i, i-1);
+ }
+ }
+ }
+ }
+}
+
+/* Given an M x N integer matrix A, this function determines an M x M
+ unimodular matrix V, and an M x N echelon matrix S such that "A =
+ V.S". This decomposition is also known as "left Hermite".
+
+ Ref: Algorithm 2.2 page 36 in "Loop Transformations for
+ Restructuring Compilers" Utpal Banerjee. */
+
+void
+lambda_matrix_left_hermite (lambda_matrix A, int m, int n,
+ lambda_matrix S, lambda_matrix V)
+{
+ int i, j, i0 = 0;
+
+ lambda_matrix_copy (A, S, m, n);
+ lambda_matrix_id (V, m);
+
+ for (j = 0; j < n; j++)
+ {
+ if (lambda_vector_first_nz (S[j], m, i0) < m)
+ {
+ ++i0;
+ for (i = m - 1; i >= i0; i--)
+ {
+ while (S[i][j] != 0)
+ {
+ int sigma, factor, a, b;
+
+ a = S[i-1][j];
+ b = S[i][j];
+ sigma = (a * b < 0) ? -1: 1;
+ a = abs (a);
+ b = abs (b);
+ factor = sigma * (a / b);
+
+ lambda_matrix_row_add (S, n, i, i-1, -factor);
+ lambda_matrix_row_exchange (S, i, i-1);
+
+ lambda_matrix_col_add (V, m, i-1, i, factor);
+ lambda_matrix_col_exchange (V, m, i, i-1);
+ }
+ }
+ }
+ }
+}
+
+/* When it exists, return the first non-zero row in MAT after row
+ STARTROW. Otherwise return rowsize. */
+
+int
+lambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize,
+ int startrow)
+{
+ int j;
+ bool found = false;
+
+ for (j = startrow; (j < rowsize) && !found; j++)
+ {
+ if ((mat[j] != NULL)
+ && (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize))
+ found = true;
+ }
+
+ if (found)
+ return j - 1;
+ return rowsize;
+}
+
+/* Calculate the projection of E sub k to the null space of B. */
+
+void
+lambda_matrix_project_to_null (lambda_matrix B, int rowsize,
+ int colsize, int k, lambda_vector x)
+{
+ lambda_matrix M1, M2, M3, I;
+ int determinant;
+
+ /* compute c(I-B^T inv(B B^T) B) e sub k */
+
+ /* M1 is the transpose of B */
+ M1 = lambda_matrix_new (colsize, colsize);
+ lambda_matrix_transpose (B, M1, rowsize, colsize);
+
+ /* M2 = B * B^T */
+ M2 = lambda_matrix_new (colsize, colsize);
+ lambda_matrix_mult (B, M1, M2, rowsize, colsize, rowsize);
+
+ /* M3 = inv(M2) */
+ M3 = lambda_matrix_new (colsize, colsize);
+ determinant = lambda_matrix_inverse (M2, M3, rowsize);
+
+ /* M2 = B^T (inv(B B^T)) */
+ lambda_matrix_mult (M1, M3, M2, colsize, rowsize, rowsize);
+
+ /* M1 = B^T (inv(B B^T)) B */
+ lambda_matrix_mult (M2, B, M1, colsize, rowsize, colsize);
+ lambda_matrix_negate (M1, M1, colsize, colsize);
+
+ I = lambda_matrix_new (colsize, colsize);
+ lambda_matrix_id (I, colsize);
+
+ lambda_matrix_add_mc (I, determinant, M1, 1, M2, colsize, colsize);
+
+ lambda_matrix_get_column (M2, colsize, k - 1, x);
+
+}
+
+/* Multiply a vector VEC by a matrix MAT.
+ MAT is an M*N matrix, and VEC is a vector with length N. The result
+ is stored in DEST which must be a vector of length M. */
+
+void
+lambda_matrix_vector_mult (lambda_matrix matrix, int m, int n,
+ lambda_vector vec, lambda_vector dest)
+{
+ int i, j;
+
+ lambda_vector_clear (dest, m);
+ for (i = 0; i < m; i++)
+ for (j = 0; j < n; j++)
+ dest[i] += matrix[i][j] * vec[j];
+}
+
+/* Print out an M x N matrix MAT to OUTFILE. */
+
+void
+print_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n)
+{
+ int i;
+
+ for (i = 0; i < m; i++)
+ print_lambda_vector (outfile, matrix[i], n);
+ fprintf (outfile, "\n");
+}
+
diff --git a/gcc/lambda.h b/gcc/lambda.h
index 60ed8accbf2cbae7110fbba3417321f6362714e8..998e3f18c62c08b6c36aef96f1b024fb346c32c3 100644
--- a/gcc/lambda.h
+++ b/gcc/lambda.h
@@ -1,4 +1,4 @@
-/* Lambda matrix interface.
+/* Lambda matrix and vector interface.
Copyright (C) 2003, 2004 Free Software Foundation, Inc.
Contributed by Daniel Berlin <dberlin@dberlin.org>
@@ -18,11 +18,65 @@ You should have received a copy of the GNU General Public License
along with GCC; see the file COPYING. If not, write to the Free
Software Foundation, 59 Temple Place - Suite 330, Boston, MA
02111-1307, USA. */
-
+
#ifndef LAMBDA_H
#define LAMBDA_H
+/* An integer vector. A vector formally consists of an element of a vector
+ space. A vector space is a set that is closed under vector addition
+ and scalar multiplication. In this vector space, an element is a list of
+ integers. */
typedef int *lambda_vector;
+/* An integer matrix. A matrix consists of m vectors of length n (IE
+ all vectors are the same length). */
+typedef lambda_vector *lambda_matrix;
+
+lambda_matrix lambda_matrix_new (int, int);
+
+void lambda_matrix_id (lambda_matrix, int);
+void lambda_matrix_copy (lambda_matrix, lambda_matrix, int, int);
+void lambda_matrix_negate (lambda_matrix, lambda_matrix, int, int);
+void lambda_matrix_transpose (lambda_matrix, lambda_matrix, int, int);
+void lambda_matrix_add (lambda_matrix, lambda_matrix, lambda_matrix, int,
+ int);
+void lambda_matrix_add_mc (lambda_matrix, int, lambda_matrix, int,
+ lambda_matrix, int, int);
+void lambda_matrix_mult (lambda_matrix, lambda_matrix, lambda_matrix,
+ int, int, int);
+void lambda_matrix_delete_rows (lambda_matrix, int, int, int);
+void lambda_matrix_row_exchange (lambda_matrix, int, int);
+void lambda_matrix_row_add (lambda_matrix, int, int, int, int);
+void lambda_matrix_row_negate (lambda_matrix mat, int, int);
+void lambda_matrix_row_mc (lambda_matrix, int, int, int);
+void lambda_matrix_col_exchange (lambda_matrix, int, int, int);
+void lambda_matrix_col_add (lambda_matrix, int, int, int, int);
+void lambda_matrix_col_negate (lambda_matrix, int, int);
+void lambda_matrix_col_mc (lambda_matrix, int, int, int);
+int lambda_matrix_inverse (lambda_matrix, lambda_matrix, int);
+void lambda_matrix_hermite (lambda_matrix, int, lambda_matrix, lambda_matrix);
+void lambda_matrix_left_hermite (lambda_matrix, int, int, lambda_matrix, lambda_matrix);
+void lambda_matrix_right_hermite (lambda_matrix, int, int, lambda_matrix, lambda_matrix);
+int lambda_matrix_first_nz_vec (lambda_matrix, int, int, int);
+void lambda_matrix_project_to_null (lambda_matrix, int, int, int,
+ lambda_vector);
+void print_lambda_matrix (FILE *, lambda_matrix, int, int);
+
+void lambda_matrix_vector_mult (lambda_matrix, int, int, lambda_vector,
+ lambda_vector);
+
+static inline void lambda_vector_negate (lambda_vector, lambda_vector, int);
+static inline void lambda_vector_mult_const (lambda_vector, lambda_vector, int, int);
+static inline void lambda_vector_add (lambda_vector, lambda_vector,
+ lambda_vector, int);
+static inline void lambda_vector_add_mc (lambda_vector, int, lambda_vector, int,
+ lambda_vector, int);
+static inline void lambda_vector_copy (lambda_vector, lambda_vector, int);
+static inline bool lambda_vector_zerop (lambda_vector, int);
+static inline void lambda_vector_clear (lambda_vector, int);
+static inline bool lambda_vector_equal (lambda_vector, lambda_vector, int);
+static inline int lambda_vector_min_nz (lambda_vector, int, int);
+static inline int lambda_vector_first_nz (lambda_vector, int, int);
+static inline void print_lambda_vector (FILE *, lambda_vector, int);
/* Allocate a new vector of given SIZE. */
@@ -32,14 +86,149 @@ lambda_vector_new (int size)
return ggc_alloc_cleared (size * sizeof(int));
}
+
+
+/* Multiply vector VEC1 of length SIZE by a constant CONST1,
+ and store the result in VEC2. */
+
+static inline void
+lambda_vector_mult_const (lambda_vector vec1, lambda_vector vec2,
+ int size, int const1)
+{
+ int i;
+
+ if (const1 == 0)
+ lambda_vector_clear (vec2, size);
+ else
+ for (i = 0; i < size; i++)
+ vec2[i] = const1 * vec1[i];
+}
+
+/* Negate vector VEC1 with length SIZE and store it in VEC2. */
+
+static inline void
+lambda_vector_negate (lambda_vector vec1, lambda_vector vec2,
+ int size)
+{
+ lambda_vector_mult_const (vec1, vec2, size, -1);
+}
+
+/* VEC3 = VEC1+VEC2, where all three the vectors are of length SIZE. */
+
+static inline void
+lambda_vector_add (lambda_vector vec1, lambda_vector vec2,
+ lambda_vector vec3, int size)
+{
+ int i;
+ for (i = 0; i < size; i++)
+ vec3[i] = vec1[i] + vec2[i];
+}
+
+/* VEC3 = CONSTANT1*VEC1 + CONSTANT2*VEC2. All vectors have length SIZE. */
+
+static inline void
+lambda_vector_add_mc (lambda_vector vec1, int const1,
+ lambda_vector vec2, int const2,
+ lambda_vector vec3, int size)
+{
+ int i;
+ for (i = 0; i < size; i++)
+ vec3[i] = const1 * vec1[i] + const2 * vec2[i];
+}
+
+/* Copy the elements of vector VEC1 with length SIZE to VEC2. */
+
+static inline void
+lambda_vector_copy (lambda_vector vec1, lambda_vector vec2,
+ int size)
+{
+ memcpy (vec2, vec1, size * sizeof (*vec1));
+}
+
+/* Return true if vector VEC1 of length SIZE is the zero vector. */
+
+static inline bool
+lambda_vector_zerop (lambda_vector vec1, int size)
+{
+ int i;
+ for (i = 0; i < size; i++)
+ if (vec1[i] != 0)
+ return false;
+ return true;
+}
+
/* Clear out vector VEC1 of length SIZE. */
static inline void
lambda_vector_clear (lambda_vector vec1, int size)
{
- memset (vec1, 0, size * sizeof (int));
+ memset (vec1, 0, size * sizeof (*vec1));
}
+/* Return true if two vectors are equal. */
+
+static inline bool
+lambda_vector_equal (lambda_vector vec1, lambda_vector vec2, int size)
+{
+ int i;
+ for (i = 0; i < size; i++)
+ if (vec1[i] != vec2[i])
+ return false;
+ return true;
+}
+
+/* Return the minimum non-zero element in vector VEC1 between START and N.
+ We must have START <= N. */
+
+static inline int
+lambda_vector_min_nz (lambda_vector vec1, int n, int start)
+{
+ int j;
+ int min = -1;
+#ifdef ENABLE_CHECKING
+ if (start > n)
+ abort ();
+#endif
+ for (j = start; j < n; j++)
+ {
+ if (vec1[j])
+ if (min < 0 || vec1[j] < vec1[min])
+ min = j;
+ }
+
+ if (min < 0)
+ abort ();
+
+ return min;
+}
+
+/* Return the first nonzero element of vector VEC1 between START and N.
+ We must have START <= N. Returns N if VEC1 is the zero vector. */
+
+static inline int
+lambda_vector_first_nz (lambda_vector vec1, int n, int start)
+{
+ int j = start;
+ while (j < n && vec1[j] == 0)
+ j++;
+ return j;
+}
+
+
+/* Multiply a vector by a matrix. */
+
+static inline void
+lambda_vector_matrix_mult (lambda_vector vect, int m, lambda_matrix mat,
+ int n, lambda_vector dest)
+{
+ int i, j;
+ lambda_vector_clear (dest, n);
+ for (i = 0; i < n; i++)
+ for (j = 0; j < m; j++)
+ dest[i] += mat[j][i] * vect[j];
+}
+
+
/* Print out a vector VEC of length N to OUTFILE. */
static inline void
@@ -51,7 +240,5 @@ print_lambda_vector (FILE * outfile, lambda_vector vector, int n)
fprintf (outfile, "%3d ", vector[i]);
fprintf (outfile, "\n");
}
-
-
#endif /* LAMBDA_H */