// \$Date: 2001/07/07 15:53:59 \$ \$Author: Manuel Kauers \$ \$Revision: 1.2 \$ /** This file contains the algorithmical parts of the parametric Gaussian Elimination package. It provides a function elim that determines the solution set of a parametric linear system over all the partitions of its parameter space.

The algorithm is customizable by setting particular flags. Those flags are listed below.

• COLUMN_SWAP_ALLOWED
TRUE to allow column swapping during the elimination process.
• USE_MARKOWITZ
TRUE to use an adapted version of the Markowitz criterion for pivot selection.
• USE_COLUMN_SIMPLIFICATION
TRUE to simplify a column if no constant pivot is found. The idea is to try to reduce the degree of the entries in a particular column in order to recieve an entry with zero degree which then can serve as a pivot element.
• SIMPLIFY_AFTER_BRANCH
TRUE to specify that affecting the constraint store will lead to a simplification of all the polynomial entries of the remaining submatrix. Note that this option may lead to a fill-in of new variables.
• USE_FURTHER_IMPROVEMENTS
TRUE to allow a number of non-documented heuristics which slightly improve the algorithm's efficiency.
By default, all the above flaged are set to TRUE.

It is easy to analyze the behaviour of this implementation by use of the methods defined in Test.mu. Those methods are based on some portions of code that were inserted into the methods of this file. In particular, the implementation here calls methods of Test.mu to pass several statistical results to the testing architecture. Each line of code in this file that exists for the purpose just described ends with the character sequence "// ###". So if you want to slightly speed up the implementation, you can feel free to delete all the lines ending with "// ###" -- they are not a part of the algorithm. (Design pattern "Template Method") */ /*++ elim -- normalizes a matrix w.r.t. a set of constraints by means of Gaussian Elimination. elim(A [, C]) A - matrix with entries e that are liftable to DOM_POLY via poly(e) C - (optional) set of constraints given as a ConstraintStore object If no C is given, C is assumed to be empty. elim is the user interface to invoke the Gaussian Elimination on a matrix A. The method returns a list of pairs [A', C'] where C' are specialisations of C and A' is the Gaussian normal form of A with respect to C'. ++*/ elim := proc(A : Dom::Matrix(), C = null()) local L, vars, v, i, j; save zero, width, height; begin // Once for all, what are the dimensions of the matrix? // Those values are stored in global variables to avoid frequent recalculations width := linalg::ncols(A); height := linalg::nrows(A); // This is a function that determines the set of all identifiers within an expression, // e.g. vars( f(x, g(y)) ) = {x, y}. vars := x -> if nops(x) > 1 then if type(x) = DOM_POLY then _union(vars(op(x, 1))); else _union(vars(op(x, i)) \$ i=1..nops(x)); end_if; elif op(x, 1) <> x then vars(op(x, 1)); elif type(x) = DOM_IDENT then {x}; else {}; end_if; // Use the above defined function to determine all variables inside the given matrix. v := _union(vars(A[i, j]) \$ i=1..height \$ j=1..width); // add private variables for the constraint store. v := v union {cs_var.i \$ i = 1..height * width + 1}; // add the variables of the given constraint store (if it exists) if C = null() then // No constraints were given, initialize an empty store. C := ConstraintStore(poly(0, [op(v)])); v := [op(v)]; else // The constraint store's variable list has to contain all the variables of v. //print(v); //print({op(varlist(C))}); if nops(v minus {op(varlist(C))}) > 0 then error("The given constraint store does not have enough variables."); end_if; v := varlist(C); end_if; // v is a list now. zero := poly(0, v); // this is the zero element of our polynomial ring (global definition) // lift the matrix' elements to polynomials over R[v]; for i from 1 to height do for j from 1 to width do A[i, j] := poly(A[i, j], v); if A[i, j] = FAIL then // FIXME error("Up to the moment, only matrices with polynomial entries are supported.\n". "Did not succeed to convert A[".expr2text(i).", ".expr2text(j)."] to a ". "polynomial."); end_if; end_for; end_for; // NYI: Elimination of denominators and filling them into the constraint store. elimRec(A, C, 1, 1); // Start the recursion. end_proc: /*-- elimRec -- an auxiliary method used by elim. elimRec(matrix, constraints, results, row, col, d) matrix - the whole polynomial matrix that is to be normalized at the moment constraints - ConstraintStore object representing the set of currently active constraints row - row index of the top left element of the submatrix currently active col - column index of the top left element of the submatrix currently active elimRec returns a list of pairs [A', C'] where C' is a specialization of the given constraints and A' is the normal form of the given matrix with respect to C'. This method is for internal use only -- no checking is done. Though the method is declared to use for recursive use, most of the work will be done iteratively. Recursive calls are only needed if no proper pivot element is found, i.e. the selected pivot is neither zero nor constant w.r.t. the current constraint store. In this case, one of the two branches is done by a recursive call of elimRec, while the calling procedure will continue to process the other branch itself. --*/ elimRec := proc(matrix : Dom::Matrix(), constraints : ConstraintStore, row : DOM_INT, col : DOM_INT) local r, c, pivot, tmp, i, old, // ### result, b, d; save mat; begin //print(constraints, matrix); mat := matrix; // define the matrix global. result := []; // initial list of solutions. // Simplify the matrix's entries with respect to the current constraints simplifySubmatrix(row, col, constraints); // As long as possible, do the elimination job iteratively while col <= width and row <= height do // select a pivot element pivot := choosePivot(constraints, row, col); r := op(pivot, 1); c := op(pivot, 2); pivot := mat[r, c]; // swap rows and columns to get pivot at [row, col] reimplementing // mat := linalg::swapRow(mat, row, r, 1..width); // mat := linalg::swapCol(mat, col, c, 1..height); // to avoid matrix copy if row <> r then for i from 1 to width do tmp := mat[row, i]; mat[row, i] := mat[r, i]; mat[r, i] := tmp; end_for; end_if; if col <> c then for i from 1 to height do tmp := mat[i, col]; mat[i, col] := mat[i, c]; mat[i, c] := tmp; end_for; end_if; // What do we know about the selected pivot? if isConstant(constraints, pivot) then // 1. pivot is constant, thus we can use it for an elimination step. beginElimination(); // ### // Using fraction free elimination for r from row + 1 to height do // for all rows... tmp := mat[r, col]; if not iszero(tmp) then for c from col + 1 to width do // ...eliminate row old := mat[r, c]; // ### mat[r, c] := pivot * mat[r, c] - mat[row, c] * tmp; eliminationStep(old, mat[r, c]); // ### end_for; // divide row into its gccd d := gccd(constraints, mat[r, c] \$ c=col+1..width); if degree(d) > 0 then for c from col + 1 to width do mat[r, c] := mat[r, c] / d; end_for; end_if; mat[r, col] := zero; end_if; end_for; endElimination(); // ### row := row + 1; // next row elif not isZero(constraints, pivot) then // 2. pivot is neither constant nor zero, thus we have to branch on those // two possibilities. // We will pass the 'isConstant'-case to a recursive call and do the 'isZero' // case ourself. (Due to technical reasons, this is not easily exchangable!) result := result . elimRec(mat, addConstant(constraints, pivot), row, col); constraints := addZero(constraints, pivot); // 'our case' // print(constraints); //print(constraints, matrix); // this is not necessarily covered by simplifySubmatrix mat[row, col] := zero; // Simplify the matrix's entries with respect to the current constraints simplifySubmatrix(row, col, constraints); // There may be other non-constant, non-zero entries in the current column, // thus, we will visit the current column again. col := col - 1; // incremented immediately // else // 3. pivot is zero, thus nothing has to be done. end_if; col := col + 1; // go to next column end_while; // matches "while col <= width and row <= height do ..." // Leaf of this search path is reached, add the solution to the list of results // At this point, you may insert code that does further calculations with the matrix, // e.g. determining its rank, determinant, upper echolon form, etc., and put that // value into the result's list. // This is the only leaf of the seach tree in the whole algorithm. // NYI: Verify results append(result, [mat, constraints]); end_proc: /*-- choosePivot -- manages the choose of a privot element choosePivot(constraints, row, col) constraints - ConstriantStore object with the set of current constriants. row - the row of the top left element of the submatrix of interest col - the column of the top left element of the submatrix of interest choosePivot returns the index of a proper pivot element. The pivot index is returned as a sequence r, c containing the row index and the column index. The selection strategie is as follows:

1. Select p := row, col as temporary pivot.
2. If the Markowitz criterion is allowed, update p to the Markowitz' suggestion. If p is constant w.r.t. the ConstraintStore, return p as final selection.
3. Otherwise, try to find a strict constant entry in p's column. If there is one, return this item as final selection.
4. Otherwise, if column simplifcation is allowed, simplify p's column. After that, try again to find a non-zero entry of degree 0 in the current column. If there is one, return this item as final selection.
5. Otherwise, ask the constraint store whether one of the polynomials in the current row represents a constant. If so, return such an entry as final selection.
6. Otherwise, if the markowitz suggestion was not zero w.r.t the constraint store, return that value.
7. Otherwise, try to find the simplest polynomial which is not zero with respect to the constraint store. If you find one, return this item as final selection.
8. Otherwise, the whole column does only contain zero entries. In this case, return an arbitrary entry.