From: Daniel Berlin Date: Tue, 10 Aug 2004 20:43:05 +0000 (+0000) Subject: lambda.h: Add matrix type, and prototypes for remainder of matrix and vector functions. X-Git-Url: https://git.libre-soc.org/?a=commitdiff_plain;h=98975653223837f9835c7b626dde25c3647191a5;p=gcc.git lambda.h: Add matrix type, and prototypes for remainder of matrix and vector functions. 2004-08-10 Daniel Berlin * 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. From-SVN: r85767 --- diff --git a/gcc/ChangeLog b/gcc/ChangeLog index e8d1398d3f8..fd7b7e348b5 100644 --- a/gcc/ChangeLog +++ b/gcc/ChangeLog @@ -1,3 +1,20 @@ +2004-08-10 Daniel Berlin + + * 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 * tree-cfg.c (bsi_insert_before, bsi_insert_after): Call modify_stmt diff --git a/gcc/Makefile.in b/gcc/Makefile.in index 984bb9e0bdd..c51760a0697 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 00000000000..987fa5e1fd1 --- /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 . + +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 60ed8accbf2..998e3f18c62 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 @@ -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 */