Loading...
Searching...
No Matches
linearAlgebra.h File Reference
`#include "kernel/structs.h"`

Go to the source code of this file.

## Functions

void | luDecomp (const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring r=currRing)
| LU-decomposition of a given (m x n)-matrix.
int | pivotScore (number n, const ring r=currRing)
| Returns a pivot score for any non-zero matrix entry.
bool | pivot (const matrix aMat, const int r1, const int r2, const int c1, const int c2, int *bestR, int *bestC, const ring r=currRing)
| Finds the best pivot element in some part of a given matrix.
bool | luInverse (const matrix aMat, matrix &iMat, const ring r=currRing)
| Computes the inverse of a given (n x n)-matrix.
bool | upperRightTriangleInverse (const matrix uMat, matrix &iMat, bool diagonalIsOne, const ring r=currRing)
| Computes the inverse of a given (n x n)-matrix in upper right triangular form.
bool | lowerLeftTriangleInverse (const matrix lMat, matrix &iMat, bool diagonalIsOne)
| Computes the inverse of a given (n x n)-matrix in lower left triangular form.
bool | luInverseFromLUDecomp (const matrix pMat, const matrix lMat, const matrix uMat, matrix &iMat, const ring r=currRing)
| Computes the inverse of an (n x n)-matrix which is given by its LU- decomposition.
int | luRank (const matrix aMat, const bool isRowEchelon, const ring r=currRing)
| Computes the rank of a given (m x n)-matrix.
bool | luSolveViaLUDecomp (const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
| Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decomposition.
bool | luSolveViaLDUDecomp (const matrix pMat, const matrix lMat, const matrix dMat, const matrix uMat, const poly l, const poly u, const poly lTimesU, const matrix bVec, matrix &xVec, matrix &H)
| Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LDU-decomposition.
bool | unitMatrix (const int n, matrix &unitMat, const ring r=currRing)
| Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.
bool | subMatrix (const matrix aMat, const int rowIndex1, const int rowIndex2, const int colIndex1, const int colIndex2, matrix &subMat)
| Creates a new matrix which is a submatrix of the first argument, and returns true in case of success.
void | swapRows (int row1, int row2, matrix &aMat)
| Swaps two rows of a given matrix and thereby modifies it.
void | swapColumns (int column1, int column2, matrix &aMat)
| Swaps two columns of a given matrix and thereby modifies it.
void | matrixBlock (const matrix aMat, const matrix bMat, matrix &block)
| Creates a new matrix which contains the first argument in the top left corner, and the second in the bottom right.
bool | charPoly (const matrix aMat, poly &charPoly)
| Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of success.
bool | realSqrt (const number n, const number tolerance, number &root)
| Computes the square root of a non-negative real number and returns it as a new number.
void | hessenberg (const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring r)
| Computes the Hessenberg form of a given square matrix.
int | quadraticSolve (const poly p, number &s1, number &s2, const number tolerance)
| Returns all solutions of a quadratic polynomial relation with real coefficients.
void | henselFactors (const int xIndex, const int yIndex, const poly h, const poly f0, const poly g0, const int d, poly &f, poly &g)
| Computes a factorization of a polynomial h(x, y) in K[[x]][y] up to a certain degree in x, whenever a factorization of h(0, y) is given.
void | lduDecomp (const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)
| LU-decomposition of a given (m x n)-matrix with performing only those divisions that yield zero remainders.
bool | qrDS (const int n, matrix *queue, int &queueL, number *eigenValues, int &eigenValuesL, const number tol1, const number tol2, const ring R)
int | similar (const number *nn, const int nnLength, const number n, const number tolerance)
| Tries to find the number n in the array nn[0..nnLength-1].

## Function Documentation

## charPoly()

bool charPoly | ( | const matrix | aMat,
|  | poly & | charPoly )

Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of success.

The method will be successful whenever the input matrix is a (2x2) matrix. In this case, the resulting polynomial will be a univariate polynomial in the first ring variable of degree 2 and it will reside in the second argument. The method assumes that the matrix entries are all numbers, i.e., elements from the ground field/ring.

Returns
  true iff the input matrix was (2x2)

Parameters

  [in] | aMat | the input matrix
  [out] | charPoly | the characteristic polynomial

Definition at line 757 of file linearAlgebra.cc.

{
if (MATROWS(aMat) != 2) return false;
if (MATCOLS(aMat) != 2) return false;
number b = nInit(0); number t;
if (MATELEM(aMat, 1, 1) != NULL)
{ t = nAdd(b, pGetCoeff(MATELEM(aMat, 1, 1))); nDelete(&b); b = t;}
if (MATELEM(aMat, 2, 2) != NULL)
{ t = nAdd(b, pGetCoeff(MATELEM(aMat, 2, 2))); nDelete(&b); b = t;}
b = nInpNeg(b);
number t1;
if ((MATELEM(aMat, 1, 1) != NULL) && (MATELEM(aMat, 2, 2) != NULL))
t1 = nMult(pGetCoeff(MATELEM(aMat, 1, 1)),
pGetCoeff(MATELEM(aMat, 2, 2)));
else t1 = nInit(0);
number t2;
if ((MATELEM(aMat, 1, 2) != NULL) && (MATELEM(aMat, 2, 1) != NULL))
t2 = nMult(pGetCoeff(MATELEM(aMat, 1, 2)),
pGetCoeff(MATELEM(aMat, 2, 1)));
else t2 = nInit(0);
number c = nSub(t1, t2); nDelete(&t1); nDelete(&t2);
poly p = pOne(); pSetExp(p, 1, 2); pSetm(p);
poly q = NULL;
if (!nIsZero(b))
{ q = pOne(); pSetExp(q, 1, 1); pSetm(q); pSetCoeff(q, b); }
poly r = NULL;
if (!nIsZero(c))
{ r = pOne(); pSetCoeff(r, c); }
p = pAdd(p, q); p = pAdd(p, r);
charPoly = p;
return true;
}
p
int p
Definition cfModGcd.cc:4086
b
CanonicalForm b
Definition cfModGcd.cc:4111
charPoly
bool charPoly(const matrix aMat, poly &charPoly)
Computes the characteristic polynomial from a quadratic (2x2) matrix and returns true in case of succ...
Definition linearAlgebra.cc:757
MATELEM
#define MATELEM(mat, i, j)
1-based access to matrix
Definition matpol.h:29
MATROWS
#define MATROWS(i)
Definition matpol.h:26
MATCOLS
#define MATCOLS(i)
Definition matpol.h:27
pGetCoeff
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition monomials.h:44
nDelete
#define nDelete(n)
Definition numbers.h:16
nInpNeg
#define nInpNeg(n)
Definition numbers.h:21
nIsZero
#define nIsZero(n)
Definition numbers.h:19
nSub
#define nSub(n1, n2)
Definition numbers.h:22
nAdd
#define nAdd(n1, n2)
Definition numbers.h:18
nInit
#define nInit(i)
Definition numbers.h:24
nMult
#define nMult(n1, n2)
Definition numbers.h:17
NULL
#define NULL
Definition omList.c:12
pAdd
#define pAdd(p, q)
Definition polys.h:204
pSetm
#define pSetm(p)
Definition polys.h:272
pSetCoeff
#define pSetCoeff(p, n)
deletes old coeff before setting the new one
Definition polys.h:32
pSetExp
#define pSetExp(p, i, v)
Definition polys.h:43
pOne
#define pOne()
Definition polys.h:316

## henselFactors()

void henselFactors | ( | const int | xIndex,
|  | const int | yIndex,
|  | const poly | h,
|  | const poly | f0,
|  | const poly | g0,
|  | const int | d,
|  | poly & | f,
|  | poly & | g )

Computes a factorization of a polynomial h(x, y) in K[[x]][y] up to a certain degree in x, whenever a factorization of h(0, y) is given.

The algorithm is based on Hensel's lemma: Let h(x, y) denote a monic polynomial in y of degree m + n with coefficients in K[[x]]. Suppose there are two monic factors f_0(y) (of degree n) and g_0(y) of degree (m) such that h(0, y) = f_0(y) * g_0(y) and <f_0, g_0> = K[y]. Fix an integer d >= 0. Then there are monic polynomials in y with coefficients in K[[x]], namely f(x, y) of degree n and g(x, y) of degree m such that h(x, y) = f(x, y) * g(x, y) modulo <x^(d+1)> (*).

This implementation expects h, f0, g0, and d as described and computes the factors f and g. Effectively, h will be given as an element of K[x, y] since all terms of h with x-degree larger than d can be ignored due to (*). The method expects the ground ring to contain at least two variables; then x is the first ring variable, specified by the input index xIndex, and y the second one, specified by yIndex.

This code was placed here since the algorithm works by successively solving d linear equation systems. It is hence an application of other methods defined in this h-file and its corresponding cc-file.

Parameters

  [in] | xIndex | index of x in {1, ..., nvars(basering)}
  [in] | yIndex | index of y in {1, ..., nvars(basering)}
  [in] | h | the polynomial h(x, y) as above
  [in] | f0 | the first univariate factor of h(0, y)
  [in] | g0 | the second univariate factor of h(0, y)
  [in] | d | the degree bound, d >= 0
  [out] | f | the first factor of h(x, y)
  [out] | g | the second factor of h(x, y)

Definition at line 1228 of file linearAlgebra.cc.

{
int n = (int)p_Deg(f0,currRing);
int m = (int)p_Deg(g0,currRing);
matrix aMat = mpNew(n + m, n + m); /* matrix A for linear system */
matrix pMat; matrix lMat; matrix uMat; /* for the decomposition of A */
f = pCopy(f0); g = pCopy(g0); /* initially: h = f*g mod <x^1> */
/* initial step: read off coefficients of f0, and g0 */
poly p = f0; poly matEntry; number c;
while (p != NULL)
{
c = nCopy(pGetCoeff(p));
matEntry = pOne(); pSetCoeff(matEntry, c);
MATELEM(aMat, pGetExp(p, yIndex) + 1, 1) = matEntry;
p = pNext(p);
}
p = g0;
while (p != NULL)
{
c = nCopy(pGetCoeff(p));
matEntry = pOne(); pSetCoeff(matEntry, c);
MATELEM(aMat, pGetExp(p, yIndex) + 1, m + 1) = matEntry;
p = pNext(p);
}
/* fill the rest of A */
for (int row = 2; row <= n + 1; row++)
for (int col = 2; col <= m; col++)
{
if (col > row) break;
MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
}
for (int row = n + 2; row <= n + m; row++)
for (int col = row - n; col <= m; col++)
MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
for (int row = 2; row <= m + 1; row++)
for (int col = m + 2; col <= m + n; col++)
{
if (col - m > row) break;
MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
}
for (int row = m + 2; row <= n + m; row++)
for (int col = row; col <= m + n; col++)
MATELEM(aMat, row, col) = pCopy(MATELEM(aMat, row - 1, col - 1));
/* constructing the LU-decomposition of A */
luDecomp(aMat, pMat, lMat, uMat);
/* Before the xExp-th loop, we know that h = f*g mod <x^xExp>.
Afterwards the algorithm ensures h = f*g mod <x^(xExp + 1)>.
Hence in the end we obtain f and g as required, i.e.,
h = f*g mod <x^(d+1)>.
The algorithm works by solving a (m+n)x(m+n) linear equation system
A*x = b with constant matrix A (as decomposed above). By theory, the
system is guaranteed to have a unique solution. */
poly fg = ppMult_qq(f, g); /* for storing the product of f and g */
for (int xExp = 1; xExp <= d; xExp++)
{
matrix bVec = mpNew(n + m, 1); /* b */
matrix xVec = mpNew(n + m, 1); /* x */
p = pCopy(fg);
p = pAdd(pCopy(h), pNeg(p)); /* p = h - f*g */
/* we collect all terms in p with x-exponent = xExp and use their
coefficients to build the vector b */
int bIsZeroVector = true;
while (p != NULL)
{
if (pGetExp(p, xIndex) == xExp)
{
number c = nCopy(pGetCoeff(p));
poly matEntry = pOne(); pSetCoeff(matEntry, c);
MATELEM(bVec, pGetExp(p, yIndex) + 1, 1) = matEntry;
if (matEntry != NULL) bIsZeroVector = false;
}
pLmDelete(&p); /* destruct leading term of p and move to next term */
}
/* solve the linear equation system */
if (!bIsZeroVector) /* otherwise x = 0 and f, g do not change */
{
matrix notUsedMat;
luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, notUsedMat);
idDelete((ideal*)&notUsedMat);
/* update f and g by newly computed terms, and update f*g */
poly fNew = NULL; poly gNew = NULL;
for (int row = 1; row <= m; row++)
{
if (MATELEM(xVec, row, 1) != NULL)
{
p = pCopy(MATELEM(xVec, row, 1)); /* p = c */
pSetExp(p, xIndex, xExp); /* p = c * x^xExp */
pSetExp(p, yIndex, row - 1); /* p = c * x^xExp * y^i */
pSetm(p);
gNew = pAdd(gNew, p);
}
}
for (int row = m + 1; row <= m + n; row++)
{
if (MATELEM(xVec, row, 1) != NULL)
{
p = pCopy(MATELEM(xVec, row, 1)); /* p = c */
pSetExp(p, xIndex, xExp); /* p = c * x^xExp */
pSetExp(p, yIndex, row - m - 1); /* p = c * x^xExp * y^i */
pSetm(p);
fNew = pAdd(fNew, p);
}
}
fg = pAdd(fg, ppMult_qq(f, gNew));
fg = pAdd(fg, ppMult_qq(g, fNew));
fg = pAdd(fg, ppMult_qq(fNew, gNew));
f = pAdd(f, fNew);
g = pAdd(g, gNew);
}
/* clean-up loop-dependent vectors */
idDelete((ideal*)&bVec); idDelete((ideal*)&xVec);
}
/* clean-up matrices A, P, L and U, and polynomial fg */
idDelete((ideal*)&aMat); idDelete((ideal*)&pMat);
idDelete((ideal*)&lMat); idDelete((ideal*)&uMat);
pDelete(&fg);
}
m
int m
Definition cfEzgcd.cc:128
g
g
Definition cfModGcd.cc:4098
f
FILE * f
Definition checklibs.c:9
idDelete
#define idDelete(H)
delete an ideal
Definition ideals.h:29
h
STATIC_VAR Poly * h
Definition janet.cc:971
luDecomp
void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat, const ring R)
LU-decomposition of a given (m x n)-matrix.
Definition linearAlgebra.cc:103
luSolveViaLUDecomp
bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, const matrix bVec, matrix &xVec, matrix &H)
Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decompositi...
Definition linearAlgebra.cc:377
mpNew
matrix mpNew(int r, int c)
create a r x c zero-matrix
Definition matpol.cc:37
matrix
ip_smatrix * matrix
Definition matpol.h:43
pNext
#define pNext(p)
Definition monomials.h:36
nCopy
#define nCopy(n)
Definition numbers.h:15
p_Deg
long p_Deg(poly a, const ring r)
Definition p_polys.cc:586
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
pDelete
#define pDelete(p_ptr)
Definition polys.h:187
pNeg
#define pNeg(p)
Definition polys.h:199
ppMult_qq
#define ppMult_qq(p, q)
Definition polys.h:209
pLmDelete
#define pLmDelete(p)
assume p != NULL, deletes Lm(p)->coef and Lm(p)
Definition polys.h:77
pGetExp
#define pGetExp(p, i)
Exponent.
Definition polys.h:42
pCopy
#define pCopy(p)
return a copy of the poly
Definition polys.h:186

## hessenberg()

void hessenberg | ( | const matrix | aMat,
|  | matrix & | pMat,
|  | matrix & | hessenbergMat,
|  | const number | tolerance,
|  | const ring | r )

Computes the Hessenberg form of a given square matrix.

The method assumes that all matrix entries are numbers coming from some subfield of the reals.. Afterwards, the following conditions will hold: 1) hessenbergMat = pMat * aMat * pMat^{-1}, 2) hessenbergMat is in Hessenberg form. During the algorithm, pMat will be constructed as the product of self- inverse matrices. The algorithm relies on computing square roots of floating point numbers. These will be combuted by Heron's iteration formula, with iteration stopping when two successive approximations of the square root differ by no more than the given tolerance, which is assumed to be a positive real number.

Parameters

  [in] | aMat | the square input matrix
  [out] | pMat | the transformation matrix
  [out] | hessenbergMat | the Hessenberg form of aMat
  [in] | tolerance | accuracy

Definition at line 918 of file linearAlgebra.cc.

{
int n = MATROWS(aMat);
unitMatrix(n, pMat);
subMatrix(aMat, 1, n, 1, n, hessenbergMat);
for (int c = 1; c <= n; c++)
{
/* find one or two non-zero entries in the current column */
int r1 = 0; int r2 = 0;
for (int r = c + 1; r <= n; r++)
if (MATELEM(hessenbergMat, r, c) != NULL)
{
if (r1 == 0) r1 = r;
else if (r2 == 0) { r2 = r; break; }
}
if (r1 != 0)
{ /* At least one entry in the current column is non-zero. */
if (r1 != c + 1)
{ /* swap rows to bring non-zero element to row with index c + 1 */
swapRows(r1, c + 1, hessenbergMat);
/* now also swap columns to reflect action of permutation
from the right-hand side */
swapColumns(r1, c + 1, hessenbergMat);
/* include action of permutation also in pMat */
swapRows(r1, c + 1, pMat);
}
if (r2 != 0)
{ /* There is at least one more non-zero element in the current
column. So let us perform a hessenberg step in order to make
all additional non-zero elements zero. */
matrix v; subMatrix(hessenbergMat, c + 1, n, c, c, v);
matrix u; matrix pTmp;
number r = hessenbergStep(v, u, pTmp, tolerance);
idDelete((ideal*)&v); idDelete((ideal*)&u); nDelete(&r);
/* pTmp just acts on the lower right block of hessenbergMat;
i.e., it needs to be extended by a unit matrix block at the top
left in order to define a whole transformation matrix;
this code may be optimized */
unitMatrix(c, u);
matrix pTmpFull; matrixBlock(u, pTmp, pTmpFull);
idDelete((ideal*)&u); idDelete((ideal*)&pTmp);
/* now include pTmpFull in pMat (by letting it act from the left) */
pTmp = mp_Mult(pTmpFull, pMat,R); idDelete((ideal*)&pMat);
pMat = pTmp;
/* now let pTmpFull act on hessenbergMat from the left and from the
right (note that pTmpFull is self-inverse) */
pTmp = mp_Mult(pTmpFull, hessenbergMat,R);
idDelete((ideal*)&hessenbergMat);
hessenbergMat = mp_Mult(pTmp, pTmpFull, R);
idDelete((ideal*)&pTmp); idDelete((ideal*)&pTmpFull);
/* as there may be inaccuracy, we erase those entries of hessenbergMat
which must have become zero by the last transformation */
for (int r = c + 2; r <= n; r++)
pDelete(&MATELEM(hessenbergMat, r, c));
}
}
}
}
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
hessenbergStep
number hessenbergStep(const matrix vVec, matrix &uVec, matrix &pMat, const number tolerance)
Computes information related to one householder transformation step for constructing the Hessenberg f...
Definition linearAlgebra.cc:846
swapColumns
void swapColumns(int column1, int column2, matrix &aMat)
Swaps two columns of a given matrix and thereby modifies it.
Definition linearAlgebra.cc:802
matrixBlock
void matrixBlock(const matrix aMat, const matrix bMat, matrix &block)
Creates a new matrix which contains the first argument in the top left corner, and the second in the ...
Definition linearAlgebra.cc:814
unitMatrix
bool unitMatrix(const int n, matrix &unitMat, const ring R)
Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.
Definition linearAlgebra.cc:95
swapRows
void swapRows(int row1, int row2, matrix &aMat)
Swaps two rows of a given matrix and thereby modifies it.
Definition linearAlgebra.cc:790
subMatrix
bool subMatrix(const matrix aMat, const int rowIndex1, const int rowIndex2, const int colIndex1, const int colIndex2, matrix &subMat)
Creates a new matrix which is a submatrix of the first argument, and returns true in case of success.
Definition linearAlgebra.cc:742
mp_Mult
matrix mp_Mult(matrix a, matrix b, const ring R)
Definition matpol.cc:206
R
#define R
Definition sirandom.c:27

## lduDecomp()

void lduDecomp | ( | const matrix | aMat,
|  | matrix & | pMat,
|  | matrix & | lMat,
|  | matrix & | dMat,
|  | matrix & | uMat,
|  | poly & | l,
|  | poly & | u,
|  | poly & | lTimesU )

LU-decomposition of a given (m x n)-matrix with performing only those divisions that yield zero remainders.

Given an (m x n) matrix A, the method computes its LDU-decomposition, that is, it computes matrices P, L, D, and U such that

- P * A = L * D^(-1) * U,
- P is an (m x m) permutation matrix, i.e., its row/columns form the standard basis of K^m,
- L is an (m x m) matrix in lower triangular form of full rank,
- D is an (m x m) diagonal matrix with full rank, and
- U is an (m x n) matrix in upper row echelon form.
  From these conditions, it follows immediately that also A = P * L * D^(-1) * U, since P is self-inverse.

In contrast to luDecomp, this method only performs those divisions which yield zero remainders. Hence, when e.g. computing over a rational function field and starting with polynomial entries only (or over Q and starting with integer entries only), then any newly computed matrix entry will again be polynomial (or integer).

The method will modify all argument matrices except aMat, so that they will finally contain the matrices P, L, D, and U as specified above. Moreover, in order to further speed up subsequent calls of luSolveViaLDUDecomp, the two denominators and their product will also be returned.

Parameters

  [in] | aMat | the initial matrix A,
  [out] | pMat | the row permutation matrix P
  [out] | lMat | the lower triangular matrix L
  [out] | dMat | the diagonal matrix D
  [out] | uMat | the upper row echelon matrix U
  [out] | l | the product of pivots of L
  [out] | u | the product of pivots of U
  [out] | lTimesU | the product of l and u

Definition at line 1352 of file linearAlgebra.cc.

{
int rr = aMat->rows();
int cc = aMat->cols();
/* we use an int array to store all row permutations;
note that we only make use of the entries [1..rr] */
int* permut = new int[rr + 1];
for (int i = 1; i <= rr; i++) permut[i] = i;
/* fill lMat and dMat with the (rr x rr) unit matrix */
unitMatrix(rr, lMat);
unitMatrix(rr, dMat);
uMat = mpNew(rr, cc);
/* copy aMat into uMat: */
for (int r = 1; r <= rr; r++)
for (int c = 1; c <= cc; c++)
MATELEM(uMat, r, c) = pCopy(MATELEM(aMat, r, c));
u = pOne(); l = pOne();
int col = 1; int row = 1;
while ((col <= cc) & (row < rr))
{
int pivotR; int pivotC; bool pivotValid = false;
while (col <= cc)
{
pivotValid = pivot(uMat, row, rr, col, col, &pivotR, &pivotC);
if (pivotValid) break;
col++;
}
if (pivotValid)
{
if (pivotR != row)
{
swapRows(row, pivotR, uMat);
poly p = MATELEM(dMat, row, row);
MATELEM(dMat, row, row) = MATELEM(dMat, pivotR, pivotR);
MATELEM(dMat, pivotR, pivotR) = p;
swapColumns(row, pivotR, lMat);
swapRows(row, pivotR, lMat);
int temp = permut[row];
permut[row] = permut[pivotR]; permut[pivotR] = temp;
}
/* in gg, we compute the gcd of all non-zero elements in
uMat[row..rr, col];
the next number is the pivot and thus guaranteed to be different
from zero: */
number gg = nCopy(pGetCoeff(MATELEM(uMat, row, col))); number t;
for (int r = row + 1; r <= rr; r++)
{
if (MATELEM(uMat, r, col) != NULL)
{
t = gg;
gg = n_Gcd(t, pGetCoeff(MATELEM(uMat, r, col)),currRing->cf);
nDelete(&t);
}
}
t = nDiv(pGetCoeff(MATELEM(uMat, row, col)), gg);
nNormalize(t); /* this division works without remainder */
if (!nIsOne(t))
{
for (int r = row; r <= rr; r++)
MATELEM(dMat, r, r)=__p_Mult_nn(MATELEM(dMat, r, r), t,currRing);
MATELEM(lMat, row, row)=__p_Mult_nn(MATELEM(lMat, row, row), t,currRing);
}
l = pMult(l, pCopy(MATELEM(lMat, row, row)));
u = pMult(u, pCopy(MATELEM(uMat, row, col)));
for (int r = row + 1; r <= rr; r++)
{
if (MATELEM(uMat, r, col) != NULL)
{
number g = n_Gcd(pGetCoeff(MATELEM(uMat, row, col)),
pGetCoeff(MATELEM(uMat, r, col)),
currRing->cf);
number f1 = nDiv(pGetCoeff(MATELEM(uMat, r, col)), g);
nNormalize(f1); /* this division works without remainder */
number f2 = nDiv(pGetCoeff(MATELEM(uMat, row, col)), g);
nNormalize(f2); /* this division works without remainder */
pDelete(&MATELEM(uMat, r, col)); MATELEM(uMat, r, col) = NULL;
for (int c = col + 1; c <= cc; c++)
{
poly p = MATELEM(uMat, r, c);
p=__p_Mult_nn(p, f2,currRing);
poly q = pCopy(MATELEM(uMat, row, c));
q=__p_Mult_nn(q, f1,currRing); q = pNeg(q);
MATELEM(uMat, r, c) = pAdd(p, q);
}
number tt = nDiv(g, gg);
nNormalize(tt); /* this division works without remainder */
MATELEM(lMat, r, r)=__p_Mult_nn(MATELEM(lMat, r, r), tt, currRing);
nDelete(&tt);
MATELEM(lMat, r, row) = pCopy(MATELEM(lMat, r, r));
MATELEM(lMat, r, row)=__p_Mult_nn(MATELEM(lMat, r, row), f1,currRing);
nDelete(&f1); nDelete(&f2); nDelete(&g);
}
else
MATELEM(lMat, r, r)=__p_Mult_nn(MATELEM(lMat, r, r), t, currRing);
}
nDelete(&t); nDelete(&gg);
col++; row++;
}
}
/* one factor in the product u might be missing: */
if (row == rr)
{
while ((col <= cc) && (MATELEM(uMat, row, col) == NULL)) col++;
if (col <= cc) u = pMult(u, pCopy(MATELEM(uMat, row, col)));
}
/* building the permutation matrix from 'permut' and computing l */
pMat = mpNew(rr, rr);
for (int r = 1; r <= rr; r++)
MATELEM(pMat, r, permut[r]) = pOne();
delete[] permut;
lTimesU = ppMult_qq(l, u);
}
l
int l
Definition cfEzgcd.cc:100
i
int i
Definition cfEzgcd.cc:132
ip_smatrix::cols
int & cols()
Definition matpol.h:24
ip_smatrix::rows
int & rows()
Definition matpol.h:23
n_Gcd
static FORCE_INLINE number n_Gcd(number a, number b, const coeffs r)
in Z: return the gcd of 'a' and 'b' in Z/nZ, Z/2^kZ: computed as in the case Z in Z/pZ,...
Definition coeffs.h:667
pivot
bool pivot(const matrix aMat, const int r1, const int r2, const int c1, const int c2, int *bestR, int *bestC, const ring R)
This code computes a score for each non-zero matrix entry in aMat[r1..r2, c1..c2].
Definition linearAlgebra.cc:68
nDiv
#define nDiv(a, b)
Definition numbers.h:32
nIsOne
#define nIsOne(n)
Definition numbers.h:25
nNormalize
#define nNormalize(n)
Definition numbers.h:30
__p_Mult_nn
#define __p_Mult_nn(p, n, r)
Definition p_polys.h:973
pMult
#define pMult(p, q)
Definition polys.h:208

## lowerLeftTriangleInverse()

bool lowerLeftTriangleInverse | ( | const matrix | lMat,
|  | matrix & | iMat,
|  | bool | diagonalIsOne )

Computes the inverse of a given (n x n)-matrix in lower left triangular form.

This method expects an (n x n)-matrix, that is, it must have as many rows as columns. Moreover, lMat[i,j] = 0, at least for all j > i, that ism lMat is in lower left triangular form.
If the argument diagonalIsOne is true, then we know additionally, that lMat[i, i] = 1, for all i. In this case lMat is invertible. Contrariwise, if diagonalIsOne is false, we do not know anything about the diagonal entries. (Note that they may still all be 1.)
In general, there are two cases:
1) The matrix is not invertible. Then the method returns false, and &iMat remains unchanged.
2) The matrix is invertible. Then the method returns true, and the content of iMat is the inverse of lMat.

Returns
  true iff lMat is invertible, false otherwise

Parameters

  [in] | lMat | (n x n)-matrix in lower left triangular form
  [out] | iMat | inverse of lMat if invertible
  [in] | diagonalIsOne | if true, then all diagonal entries of lMat are 1

Definition at line 300 of file linearAlgebra.cc.

{
int d = lMat->rows(); poly p; poly q;
/* check whether lMat is invertible */
bool invertible = diagonalIsOne;
if (!invertible)
{
invertible = true;
for (int r = 1; r <= d; r++)
{
if (MATELEM(lMat, r, r) == NULL)
{
invertible = false;
break;
}
}
}
if (invertible)
{
iMat = mpNew(d, d);
for (int c = d; c >= 1; c--)
{
if (diagonalIsOne)
MATELEM(iMat, c, c) = pOne();
else
MATELEM(iMat, c, c) = pNSet(nInvers(pGetCoeff(MATELEM(lMat, c, c))));
for (int r = c + 1; r <= d; r++)
{
p = NULL;
for (int k = c; k <= r - 1; k++)
{
q = ppMult_qq(MATELEM(lMat, r, k), MATELEM(iMat, k, c));
p = pAdd(p, q);
}
p = pNeg(p);
p = pMult(p, pCopy(MATELEM(iMat, c, c)));
pNormalize(p);
MATELEM(iMat, r, c) = p;
}
}
}
return invertible;
}
k
int k
Definition cfEzgcd.cc:99
nInvers
#define nInvers(a)
Definition numbers.h:33
pNSet
#define pNSet(n)
Definition polys.h:314
pNormalize
#define pNormalize(p)
Definition polys.h:318

## luDecomp()

void luDecomp | ( | const matrix | aMat,
|  | matrix & | pMat,
|  | matrix & | lMat,
|  | matrix & | uMat,
|  | const ring | r = currRing )

LU-decomposition of a given (m x n)-matrix.

Given an (m x n) matrix A, the method computes its LU-decomposition, that is, it computes matrices P, L, and U such that

- P * A = L * U,
- P is an (m x m) permutation matrix, i.e., its row/columns form the standard basis of K^m,
- L is an (m x m) matrix in lower triangular form with all diagonal entries equal to 1,
- U is an (m x n) matrix in upper row echelon form.
  From these conditions, it follows immediately that also A = P * L * U, since P is self-inverse.

The method will modify all argument matrices except aMat, so that they will finally contain the matrices P, L, and U as specified above.

Parameters

  [in] | aMat | the initial matrix A,
  [out] | pMat | the row permutation matrix P
  [out] | lMat | the lower triangular matrix L
  [out] | uMat | the upper row echelon matrix U
  [in] | r | current ring

Definition at line 103 of file linearAlgebra.cc.

{
int rr = aMat->rows();
int cc = aMat->cols();
pMat = mpNew(rr, rr);
uMat = mp_Copy(aMat,R); /* copy aMat into uMat: */
/* we use an int array to store all row permutations;
note that we only make use of the entries [1..rr] */
int* permut = (int*)omAlloc((rr + 1)*sizeof(int));
for (int i = 1; i <= rr; i++) permut[i] = i;
/* fill lMat with the (rr x rr) unit matrix */
unitMatrix(rr, lMat,R);
int bestR; int bestC; int intSwap; poly pSwap; int cOffset = 0;
for (int r = 1; r < rr; r++)
{
if (r > cc) break;
while ((r + cOffset <= cc) &&
(!pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC, R)))
cOffset++;
if (r + cOffset <= cc)
{
/* swap rows with indices r and bestR in permut */
intSwap = permut[r];
permut[r] = permut[bestR];
permut[bestR] = intSwap;
/* swap rows with indices r and bestR in uMat;
it is sufficient to do this for columns >= r + cOffset*/
for (int c = r + cOffset; c <= cc; c++)
{
pSwap = MATELEM(uMat, r, c);
MATELEM(uMat, r, c) = MATELEM(uMat, bestR, c);
MATELEM(uMat, bestR, c) = pSwap;
}
/* swap rows with indices r and bestR in lMat;
we must do this only for columns < r */
for (int c = 1; c < r; c++)
{
pSwap = MATELEM(lMat, r, c);
MATELEM(lMat, r, c) = MATELEM(lMat, bestR, c);
MATELEM(lMat, bestR, c) = pSwap;
}
/* perform next Gauss elimination step, i.e., below the
row with index r;
we need to adjust lMat and uMat;
we are certain that the matrix entry at [r, r + cOffset]
is non-zero: */
number pivotElement = pGetCoeff(MATELEM(uMat, r, r + cOffset));
poly p;
for (int rGauss = r + 1; rGauss <= rr; rGauss++)
{
p = MATELEM(uMat, rGauss, r + cOffset);
if (p != NULL)
{
number n = n_Div(pGetCoeff(p), pivotElement, R->cf);
n_Normalize(n,R->cf);
/* filling lMat;
old entry was zero, so no need to call pDelete(.) */
MATELEM(lMat, rGauss, r) = p_NSet(n_Copy(n,R->cf),R);
/* adjusting uMat: */
MATELEM(uMat, rGauss, r + cOffset) = NULL; p_Delete(&p,R);
n = n_InpNeg(n,R->cf);
for (int cGauss = r + cOffset + 1; cGauss <= cc; cGauss++)
{
MATELEM(uMat, rGauss, cGauss)
= p_Add_q(MATELEM(uMat, rGauss, cGauss),
pp_Mult_nn(MATELEM(uMat, r, cGauss), n, R), R);
p_Normalize(MATELEM(uMat, rGauss, cGauss),R);
}
n_Delete(&n,R->cf); /* clean up n */
}
}
}
}
/* building the permutation matrix from 'permut' */
for (int r = 1; r <= rr; r++)
MATELEM(pMat, r, permut[r]) = p_One(R);
omFreeSize(permut,(rr + 1)*sizeof(int));
return;
}
n_Copy
static FORCE_INLINE number n_Copy(number n, const coeffs r)
return a copy of 'n'
Definition coeffs.h:457
n_InpNeg
static FORCE_INLINE number n_InpNeg(number n, const coeffs r)
in-place negation of n MUST BE USED: n = n_InpNeg(n) (no copy is returned)
Definition coeffs.h:560
n_Div
static FORCE_INLINE number n_Div(number a, number b, const coeffs r)
return the quotient of 'a' and 'b', i.e., a/b; raises an error if 'b' is not invertible in r exceptio...
Definition coeffs.h:618
n_Delete
static FORCE_INLINE void n_Delete(number *p, const coeffs r)
delete 'p'
Definition coeffs.h:461
n_Normalize
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition coeffs.h:581
mp_Copy
matrix mp_Copy(matrix a, const ring r)
copies matrix a (from ring r to r)
Definition matpol.cc:57
omFreeSize
#define omFreeSize(addr, size)
Definition omAllocDecl.h:260
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
p_Normalize
void p_Normalize(poly p, const ring r)
Definition p_polys.cc:3939
p_One
poly p_One(const ring r)
Definition p_polys.cc:1314
p_NSet
poly p_NSet(number n, const ring r)
returns the poly representing the number n, destroys n
Definition p_polys.cc:1474
p_Add_q
static poly p_Add_q(poly p, poly q, const ring r)
Definition p_polys.h:938
pp_Mult_nn
static poly pp_Mult_nn(poly p, number n, const ring r)
Definition p_polys.h:994
p_Delete
static void p_Delete(poly *p, const ring r)
Definition p_polys.h:903

## luInverse()

bool luInverse | ( | const matrix | aMat,
|  | matrix & | iMat,
|  | const ring | R )

Computes the inverse of a given (n x n)-matrix.

This method expects an (n x n)-matrix, that is, it must have as many rows as columns. Inversion of the first argument is attempted via the LU-decomposition. There are two cases:
1) The matrix is not invertible. Then the method returns false, and &iMat remains unchanged.
2) The matrix is invertible. Then the method returns true, and the content of iMat is the inverse of aMat.

Returns
  true iff aMat is invertible, false otherwise

See also
  luInverseFromLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, matrix &iMat)

Computes the inverse of a given (n x n)-matrix.

Parameters

  [in] | aMat | matrix to be inverted
  [out] | iMat | inverse of aMat if invertible
  [in] | R | current ring

Definition at line 200 of file linearAlgebra.cc.

{ /* aMat is guaranteed to be an (n x n)-matrix */
matrix pMat;
matrix lMat;
matrix uMat;
luDecomp(aMat, pMat, lMat, uMat, R);
bool result = luInverseFromLUDecomp(pMat, lMat, uMat, iMat, R);
/* clean-up */
id_Delete((ideal*)&pMat,R);
id_Delete((ideal*)&lMat,R);
id_Delete((ideal*)&uMat,R);
return result;
}
result
return result
Definition facAbsBiFact.cc:76
luInverseFromLUDecomp
bool luInverseFromLUDecomp(const matrix pMat, const matrix lMat, const matrix uMat, matrix &iMat, const ring R)
This code computes the inverse by inverting lMat and uMat, and then performing two matrix multiplicat...
Definition linearAlgebra.cc:352
id_Delete
void id_Delete(ideal *h, ring r)
deletes an ideal/module/matrix
Definition simpleideals.cc:123

## luInverseFromLUDecomp()

bool luInverseFromLUDecomp | ( | const matrix | pMat,
|  | const matrix | lMat,
|  | const matrix | uMat,
|  | matrix & | iMat,
|  | const ring | R )

Computes the inverse of an (n x n)-matrix which is given by its LU- decomposition.

With A denoting the matrix to be inverted, the method expects the LU-decomposition of A, that is, pMat * A = lMat * uMat, where the argument matrices have the appropriate properties; see method 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)'.
Furthermore, uMat is expected to be an (n x n)-matrix. Then A^(-1) is computed according to A^(-1) = uMat^(-1) * lMat^(-1) * pMat, since pMat is self-inverse. This will work if and only if uMat is invertible, because lMat and pMat are in any case invertible.
In general, there are two cases:
1) uMat and hence A is not invertible. Then the method returns false, and &iMat remains unchanged.
2) uMat and hence A is invertible. Then the method returns true, and the content of iMat is the inverse of A.

Returns
  true if A is invertible, false otherwise

See also
  luInverse(const matrix aMat, matrix &iMat)

Computes the inverse of an (n x n)-matrix which is given by its LU- decomposition.

Parameters

  [in] | pMat | permutation matrix of an LU- decomposition
  [in] | lMat | lower left matrix of an LU- decomposition
  [in] | uMat | upper right matrix of an LU- decomposition
  [out] | iMat | inverse of A if invertible

Definition at line 352 of file linearAlgebra.cc.

{ /* uMat is guaranteed to be quadratic */
//int d = uMat->rows();
matrix lMatInverse; /* for storing the inverse of lMat;
this inversion will always work */
matrix uMatInverse; /* for storing the inverse of uMat, if invertible */
bool result = upperRightTriangleInverse(uMat, uMatInverse, false);
if (result)
{
/* next will always work, since lMat is known to have all diagonal
entries equal to 1 */
lowerLeftTriangleInverse(lMat, lMatInverse, true);
iMat = mp_Mult(mp_Mult(uMatInverse, lMatInverse,R), pMat,R);
/* clean-up */
idDelete((ideal*)&lMatInverse);
idDelete((ideal*)&uMatInverse);
}
return result;
}
upperRightTriangleInverse
bool upperRightTriangleInverse(const matrix uMat, matrix &iMat, bool diagonalIsOne, const ring R)
Computes the inverse of a given (n x n)-matrix in upper right triangular form.
Definition linearAlgebra.cc:251
lowerLeftTriangleInverse
bool lowerLeftTriangleInverse(const matrix lMat, matrix &iMat, bool diagonalIsOne)
Computes the inverse of a given (n x n)-matrix in lower left triangular form.
Definition linearAlgebra.cc:300

## luRank()

int luRank | ( | const matrix | aMat,
|  | const bool | isRowEchelon,
|  | const ring | r = currRing )

Computes the rank of a given (m x n)-matrix.

The matrix may already be given in row echelon form, e.g., when the user has before called luDecomp and passes the upper triangular matrix U to luRank. In this case, the second argument can be set to true in order to pass this piece of information to the method. Otherwise, luDecomp will be called first to compute the matrix U. The rank is then read off the matrix U.

Returns
  the rank as an int

See also
  luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)

Parameters

  [in] | aMat | input matrix
  [in] | isRowEchelon | if true then aMat is assumed to be already in row echelon form

Definition at line 230 of file linearAlgebra.cc.

{
if (isRowEchelon) return rankFromRowEchelonForm(aMat);
else
{ /* compute the LU-decomposition and read off the rank from
the upper triangular matrix of that decomposition */
matrix pMat;
matrix lMat;
matrix uMat;
luDecomp(aMat, pMat, lMat, uMat,R);
int result = rankFromRowEchelonForm(uMat);
/* clean-up */
id_Delete((ideal*)&pMat,R);
id_Delete((ideal*)&lMat,R);
id_Delete((ideal*)&uMat,R);
return result;
}
}
rankFromRowEchelonForm
int rankFromRowEchelonForm(const matrix aMat)
Definition linearAlgebra.cc:217

## luSolveViaLDUDecomp()

bool luSolveViaLDUDecomp | ( | const matrix | pMat,
|  | const matrix | lMat,
|  | const matrix | dMat,
|  | const matrix | uMat,
|  | const poly | l,
|  | const poly | u,
|  | const poly | lTimesU,
|  | const matrix | bVec,
|  | matrix & | xVec,
|  | matrix & | H )

Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LDU-decomposition.

The method expects the LDU-decomposition of A, that is, pMat * A = lMat * dMat^(-1) * uMat, where the argument matrices have the appropriate properties; see method 'lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)'.
Instead of trying to invert A and return x = A^(-1)*b, this method 1) computes b' = l * pMat * b, 2) solves the simple system L * y = b', 3) computes y' = u * dMat * y, 4) solves the simple system U * y'' = y', 5) computes x = 1/(lTimesU) * y''. Note that steps 1), 2) and 3) will always work, as L is in any case invertible. Moreover, y and thus y' are uniquely determined. Step 4) will only work if and only if y' is in the column span of U. In that case, there may be infinitely many solutions. In contrast to luSolveViaLUDecomp, this algorithm guarantees that all divisions which have to be performed in steps 2) and 4) will work without remainder. This is due to properties of the given LDU- decomposition. Only the final step 5) performs a division of a vector by a member of the ground field. (Suppose the ground field is Q or some rational function field. Then, if A contains only integers or polynomials, respectively, then steps 1) - 4) will keep all data integer or polynomial, respectively. This may speed up computations considerably.) For the outcome, there are three cases:
1) There is no solution. Then the method returns false, and &xVec as well as &H remain unchanged.
2) There is a unique solution. Then the method returns true and H is the 1x1 matrix with zero entry.
3) There are infinitely many solutions. Then the method returns true and some solution of the given original linear system. Furthermore, the columns of H span the solution space of the homogeneous linear system. The dimension of the solution space is then the number of columns of H.

Returns
  true if there is at least one solution, false otherwise

Parameters

  [in] | pMat | permutation matrix of an LDU- decomposition
  [in] | lMat | lower left matrix of an LDU- decomposition
  [in] | dMat | diagonal matrix of an LDU- decomposition
  [in] | uMat | upper right matrix of an LDU- decomposition
  [in] | l | pivot product l of an LDU decomp.
  [in] | u | pivot product u of an LDU decomp.
  [in] | lTimesU | product of l and u
  [in] | bVec | right-hand side of the linear system to be solved
  [out] | xVec | solution of A*x = b
  [out] | H | matrix with columns spanning homogeneous solution space

Definition at line 1470 of file linearAlgebra.cc.

{
int m = uMat->rows(); int n = uMat->cols();
matrix cVec = mpNew(m, 1); /* for storing l * pMat * bVec */
matrix yVec = mpNew(m, 1); /* for storing the unique solution of
lMat * yVec = cVec */
/* compute cVec = l * pMat * bVec but without actual matrix mult. */
for (int r = 1; r <= m; r++)
{
for (int c = 1; c <= m; c++)
{
if (MATELEM(pMat, r, c) != NULL)
{ MATELEM(cVec, r, 1) = ppMult_qq(l, MATELEM(bVec, c, 1)); break; }
}
}
/* solve lMat * yVec = cVec; this will always work as lMat is invertible;
moreover, all divisions are guaranteed to be without remainder */
poly p; poly q;
for (int r = 1; r <= m; r++)
{
p = pNeg(pCopy(MATELEM(cVec, r, 1)));
for (int c = 1; c < r; c++)
p = pAdd(p, ppMult_qq(MATELEM(lMat, r, c), MATELEM(yVec, c, 1) ));
/* The following division is without remainder. */
q = pNSet(nInvers(pGetCoeff(MATELEM(lMat, r, r))));
MATELEM(yVec, r, 1) = pMult(pNeg(p), q);
pNormalize(MATELEM(yVec, r, 1));
}
/* compute u * dMat * yVec and put result into yVec */
for (int r = 1; r <= m; r++)
{
p = ppMult_qq(u, MATELEM(dMat, r, r));
MATELEM(yVec, r, 1) = pMult(p, MATELEM(yVec, r, 1));
}
/* determine whether uMat * xVec = yVec is solvable */
bool isSolvable = true;
bool isZeroRow; int nonZeroRowIndex;
for (int r = m; r >= 1; r--)
{
isZeroRow = true;
for (int c = 1; c <= n; c++)
if (MATELEM(uMat, r, c) != NULL) { isZeroRow = false; break; }
if (isZeroRow)
{
if (MATELEM(yVec, r, 1) != NULL) { isSolvable = false; break; }
}
else { nonZeroRowIndex = r; break; }
}
if (isSolvable)
{
xVec = mpNew(n, 1);
matrix N = mpNew(n, n); int dim = 0;
/* solve uMat * xVec = yVec and determine a basis of the solution
space of the homogeneous system uMat * xVec = 0;
We do not know in advance what the dimension (dim) of the latter
solution space will be. Thus, we start with the possibly too wide
matrix N and later copy the relevant columns of N into H. */
int nonZeroC; int lastNonZeroC = n + 1;
for (int r = nonZeroRowIndex; r >= 1; r--)
{
for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
if (MATELEM(uMat, r, nonZeroC) != NULL) break;
for (int w = lastNonZeroC - 1; w >= nonZeroC + 1; w--)
{
/* this loop will only be done when the given linear system has
more than one, i.e., infinitely many solutions */
dim++;
/* now we fill two entries of the dim-th column of N */
MATELEM(N, w, dim) = pNeg(pCopy(MATELEM(uMat, r, nonZeroC)));
MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w));
}
for (int d = 1; d <= dim; d++)
{
/* here we fill the entry of N at [r, d], for each valid vector
that we already have in N;
again, this will only be done when the given linear system has
more than one, i.e., infinitely many solutions */
p = NULL;
for (int c = nonZeroC + 1; c <= n; c++)
if (MATELEM(N, c, d) != NULL)
p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d)));
/* The following division may be with remainder but only takes place
when dim > 0. */
q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);
pNormalize(MATELEM(N, nonZeroC, d));
}
p = pNeg(pCopy(MATELEM(yVec, r, 1)));
for (int c = nonZeroC + 1; c <= n; c++)
if (MATELEM(xVec, c, 1) != NULL)
p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1)));
/* The following division is without remainder. */
q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
pNormalize(MATELEM(xVec, nonZeroC, 1));
lastNonZeroC = nonZeroC;
}
/* divide xVec by l*u = lTimesU and put result in xVec */
number z = nInvers(pGetCoeff(lTimesU));
for (int c = 1; c <= n; c++)
{
MATELEM(xVec, c, 1)=__p_Mult_nn(MATELEM(xVec, c, 1), z,currRing);
pNormalize(MATELEM(xVec, c, 1));
}
nDelete(&z);
if (dim == 0)
{
/* that means the given linear system has exactly one solution;
we return as H the 1x1 matrix with entry zero */
H = mpNew(1, 1);
}
else
{
/* copy the first 'dim' columns of N into H */
H = mpNew(n, dim);
for (int r = 1; r <= n; r++)
for (int c = 1; c <= dim; c++)
MATELEM(H, r, c) = pCopy(MATELEM(N, r, c));
}
/* clean up N */
idDelete((ideal*)&N);
}
idDelete((ideal*)&cVec);
idDelete((ideal*)&yVec);
return isSolvable;
}
N
const CanonicalForm CFMap CFMap & N
Definition cfEzgcd.cc:56
w
const CanonicalForm & w
Definition facAbsFact.cc:51
H
CanonicalForm H
Definition facAbsFact.cc:60
dim
int dim(ideal I, ring r)
Definition tropicalStrategy.cc:23

## luSolveViaLUDecomp()

bool luSolveViaLUDecomp | ( | const matrix | pMat,
|  | const matrix | lMat,
|  | const matrix | uMat,
|  | const matrix | bVec,
|  | matrix & | xVec,
|  | matrix & | H )

Solves the linear system A * x = b, where A is an (m x n)-matrix which is given by its LU-decomposition.

The method expects the LU-decomposition of A, that is, pMat * A = lMat * uMat, where the argument matrices have the appropriate properties; see method 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat)'.
Instead of trying to invert A and return x = A^(-1)*b, this method 1) computes b' = pMat * b, 2) solves the simple system L * y = b', and then 3) solves the simple system U * x = y. Note that steps 1) and 2) will always work, as L is in any case invertible. Moreover, y is uniquely determined. Step 3) will only work if and only if y is in the column span of U. In that case, there may be infinitely many solutions. Thus, there are three cases:
1) There is no solution. Then the method returns false, and &xVec as well as &H remain unchanged.
2) There is a unique solution. Then the method returns true and H is the 1x1 matrix with zero entry.
3) There are infinitely many solutions. Then the method returns true and some solution of the given original linear system. Furthermore, the columns of H span the solution space of the homogeneous linear system. The dimension of the solution space is then the number of columns of H.

Returns
  true if there is at least one solution, false otherwise

Parameters

  [in] | pMat | permutation matrix of an LU- decomposition
  [in] | lMat | lower left matrix of an LU- decomposition
  [in] | uMat | upper right matrix of an LU- decomposition
  [in] | bVec | right-hand side of the linear system to be solved
  [out] | xVec | solution of A*x = b
  [out] | H | matrix with columns spanning homogeneous solution space

Definition at line 377 of file linearAlgebra.cc.

{
int m = uMat->rows(); int n = uMat->cols();
matrix cVec = mpNew(m, 1); /* for storing pMat * bVec */
matrix yVec = mpNew(m, 1); /* for storing the unique solution of
lMat * yVec = cVec */
/* compute cVec = pMat * bVec but without actual multiplications */
for (int r = 1; r <= m; r++)
{
for (int c = 1; c <= m; c++)
{
if (MATELEM(pMat, r, c) != NULL)
{ MATELEM(cVec, r, 1) = pCopy(MATELEM(bVec, c, 1)); break; }
}
}
/* solve lMat * yVec = cVec; this will always work as lMat is invertible;
moreover, no divisions are needed, since lMat[i, i] = 1, for all i */
for (int r = 1; r <= m; r++)
{
poly p = pNeg(pCopy(MATELEM(cVec, r, 1)));
for (int c = 1; c < r; c++)
p = pAdd(p, ppMult_qq(MATELEM(lMat, r, c), MATELEM(yVec, c, 1) ));
MATELEM(yVec, r, 1) = pNeg(p);
pNormalize(MATELEM(yVec, r, 1));
}
/* determine whether uMat * xVec = yVec is solvable */
bool isSolvable = true;
bool isZeroRow;
int nonZeroRowIndex = 0 ; // handle case that the matrix is zero
for (int r = m; r >= 1; r--)
{
isZeroRow = true;
for (int c = 1; c <= n; c++)
if (MATELEM(uMat, r, c) != NULL) { isZeroRow = false; break; }
if (isZeroRow)
{
if (MATELEM(yVec, r, 1) != NULL) { isSolvable = false; break; }
}
else { nonZeroRowIndex = r; break; }
}
if (isSolvable)
{
xVec = mpNew(n, 1);
matrix N = mpNew(n, n); int dim = 0;
poly p; poly q;
/* solve uMat * xVec = yVec and determine a basis of the solution
space of the homogeneous system uMat * xVec = 0;
We do not know in advance what the dimension (dim) of the latter
solution space will be. Thus, we start with the possibly too wide
matrix N and later copy the relevant columns of N into H. */
int nonZeroC = 0 ;
int lastNonZeroC = n + 1;
for (int r = nonZeroRowIndex; r >= 1; r--)
{
for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
if (MATELEM(uMat, r, nonZeroC) != NULL) break;
for (int w = lastNonZeroC - 1; w > nonZeroC; w--)
{
/* this loop will only be done when the given linear system has
more than one, i.e., infinitely many solutions */
dim++;
/* now we fill two entries of the dim-th column of N */
MATELEM(N, w, dim) = pNeg(pCopy(MATELEM(uMat, r, nonZeroC)));
MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w));
}
for (int d = 1; d <= dim; d++)
{
/* here we fill the entry of N at [r, d], for each valid vector
that we already have in N;
again, this will only be done when the given linear system has
more than one, i.e., infinitely many solutions */
p = NULL;
for (int c = nonZeroC + 1; c <= n; c++)
if (MATELEM(N, c, d) != NULL)
p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d)));
q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);
pNormalize(MATELEM(N, nonZeroC, d));
}
p = pNeg(pCopy(MATELEM(yVec, r, 1)));
for (int c = nonZeroC + 1; c <= n; c++)
if (MATELEM(xVec, c, 1) != NULL)
p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1)));
q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
pNormalize(MATELEM(xVec, nonZeroC, 1));
lastNonZeroC = nonZeroC;
}
for (int w = lastNonZeroC - 1; w >= 1; w--)
{
// remaining variables are free
dim++;
MATELEM(N, w, dim) = pOne();
}
if (dim == 0)
{
/* that means the given linear system has exactly one solution;
we return as H the 1x1 matrix with entry zero */
H = mpNew(1, 1);
}
else
{
/* copy the first 'dim' columns of N into H */
H = mpNew(n, dim);
for (int r = 1; r <= n; r++)
for (int c = 1; c <= dim; c++)
MATELEM(H, r, c) = pCopy(MATELEM(N, r, c));
}
/* clean up N */
idDelete((ideal*)&N);
}
idDelete((ideal*)&cVec);
idDelete((ideal*)&yVec);
return isSolvable;
}

## matrixBlock()

void matrixBlock | ( | const matrix | aMat,
|  | const matrix | bMat,
|  | matrix & | block )

Creates a new matrix which contains the first argument in the top left corner, and the second in the bottom right.

All other entries of the resulting matrix which will be created in the third argument, are zero.

Parameters

  [in] | aMat | the top left input matrix
  [in] | bMat | the bottom right input matrix
  [out] | block | the new block matrix

Definition at line 814 of file linearAlgebra.cc.

{
int rowsA = MATROWS(aMat);
int rowsB = MATROWS(bMat);
int n = rowsA + rowsB;
block = mpNew(n, n);
for (int i = 1; i <= rowsA; i++)
for (int j = 1; j <= rowsA; j++)
MATELEM(block, i, j) = pCopy(MATELEM(aMat, i, j));
for (int i = 1; i <= rowsB; i++)
for (int j = 1; j <= rowsB; j++)
MATELEM(block, i + rowsA, j + rowsA) = pCopy(MATELEM(bMat, i, j));
}
j
int j
Definition facHensel.cc:110
block
#define block
Definition scanner.cc:646

## pivot()

bool pivot | ( | const matrix | aMat,
|  | const int | r1,
|  | const int | r2,
|  | const int | c1,
|  | const int | c2,
|  | int * | bestR,
|  | int * | bestC,
|  | const ring | R )

Finds the best pivot element in some part of a given matrix.

Given any matrix A with valid row indices r1..r2 and valid column indices c1..c2, this method finds the best pivot element for subsequent Gauss elimination steps in A[r1..r2, c1..c2]. "Best" here means best with respect to the implementation of the method 'pivotScore(number n)'.
In the case that all elements in A[r1..r2, c1..c2] are zero, the method returns false, otherwise true.

Returns
  false if all relevant matrix entries are zero, true otherwise

See also
  pivotScore(number n)

Finds the best pivot element in some part of a given matrix.

If all entries are zero, false is returned, otherwise true. In the latter case, the minimum of all scores is sought, and the row and column indices of the corresponding matrix entry are stored in bestR and bestC.

Parameters

  [in] | aMat | any matrix with number entries
  [in] | r1 | lower row index
  [in] | r2 | upper row index
  [in] | c1 | lower column index
  [in] | c2 | upper column index
  [out] | bestR | address of row index of best pivot element
  [out] | bestC | address of column index of best pivot element
  [in] | R | current ring

Definition at line 68 of file linearAlgebra.cc.

{
int bestScore; int score; bool foundBestScore = false; poly matEntry;
for (int c = c1; c <= c2; c++)
{
for (int r = r1; r <= r2; r++)
{
matEntry = MATELEM(aMat, r, c);
if (matEntry != NULL)
{
score = pivotScore(pGetCoeff(matEntry), R);
if ((!foundBestScore) || (score < bestScore))
{
bestScore = score;
*bestR = r;
*bestC = c;
}
foundBestScore = true;
}
}
}
return foundBestScore;
}
pivotScore
int pivotScore(number n, const ring r)
The returned score is based on the implementation of 'nSize' for numbers (, see numbers....
Definition linearAlgebra.cc:50

## pivotScore()

int pivotScore | ( | number | n,
|  | const ring | r )

Returns a pivot score for any non-zero matrix entry.

The smaller the score the better will n serve as a pivot element in subsequent Gauss elimination steps.

Returns
  the pivot score of n

Returns a pivot score for any non-zero matrix entry.

Thus, less complex pivot elements will be preferred, and get therefore a smaller pivot score. Consequently, we simply return the value of nSize. An exception to this rule are the ground fields R, long R, and long C: In these, the result of nSize relates to |n|. And since a larger modulus of the pivot element ensures a numerically more stable Gauss elimination, we would like to have a smaller score for larger values of |n|. Therefore, in R, long R, and long C, the negative of nSize will be returned.

Parameters

  [in] | n | a non-zero matrix entry
  [in] | r | current ring

Definition at line 50 of file linearAlgebra.cc.

{
int s = n_Size(n, r->cf);
if (rField_is_long_C(r) ||
rField_is_long_R(r) ||
rField_is_R(r))
return -s;
else
return s;
}
n_Size
static FORCE_INLINE int n_Size(number n, const coeffs r)
return a non-negative measure for the complexity of n; return 0 only when n represents zero; (used fo...
Definition coeffs.h:573
s
const CanonicalForm int s
Definition facAbsFact.cc:51
rField_is_R
static BOOLEAN rField_is_R(const ring r)
Definition ring.h:524
rField_is_long_C
static BOOLEAN rField_is_long_C(const ring r)
Definition ring.h:551
rField_is_long_R
static BOOLEAN rField_is_long_R(const ring r)
Definition ring.h:548

## qrDS()

bool qrDS | ( | const int | n,
|  | matrix * | queue,
|  | int & | queueL,
|  | number * | eigenValues,
|  | int & | eigenValuesL,
|  | const number | tol1,
|  | const number | tol2,
|  | const ring | R )

Definition at line 1099 of file linearAlgebra.cc.

{
bool deflationFound = true;
/* we loop until the working queue is empty,
provided we always find deflation */
while (deflationFound && (queueL > 0))
{
/* take out last queue entry */
matrix currentMat = queue[queueL - 1]; queueL--;
int m = MATROWS(currentMat);
if (m == 1)
{
number newEigenvalue;
/* the entry at [1, 1] is the eigenvalue */
if (MATELEM(currentMat, 1, 1) == NULL) newEigenvalue = nInit(0);
else newEigenvalue = nCopy(pGetCoeff(MATELEM(currentMat, 1, 1)));
eigenValues[eigenValuesL++] = newEigenvalue;
}
else if (m == 2)
{
/* there are two eigenvalues which come as zeros of the characteristic
polynomial */
poly p; charPoly(currentMat, p);
number s1; number s2;
int nSol = quadraticSolve(p, s1, s2, tol2); pDelete(&p);
assume(nSol >= 2);
eigenValues[eigenValuesL++] = s1;
/* if nSol = 2, then s1 is a double zero, and s2 is invalid: */
if (nSol == 2) s2 = nCopy(s1);
eigenValues[eigenValuesL++] = s2;
}
else /* m > 2 */
{
/* bring currentMat into Hessenberg form to fasten computations: */
matrix mm1; matrix mm2;
hessenberg(currentMat, mm1, mm2, tol2,R);
idDelete((ideal*)&currentMat); idDelete((ideal*)&mm1);
currentMat = mm2;
int it = 1; bool doLoop = true;
while (doLoop && (it <= 30 * m))
{
/* search for deflation */
number w1; number w2;
number test1; number test2; bool stopCriterion = false; int k;
for (k = 1; k < m; k++)
{
test1 = absValue(MATELEM(currentMat, k + 1, k));
w1 = absValue(MATELEM(currentMat, k, k));
w2 = absValue(MATELEM(currentMat, k + 1, k + 1));
test2 = nMult(tol1, nAdd(w1, w2));
nDelete(&w1); nDelete(&w2);
if (!nGreater(test1, test2)) stopCriterion = true;
nDelete(&test1); nDelete(&test2);
if (stopCriterion) break;
}
if (k < m) /* found deflation at position (k + 1, k) */
{
pDelete(&MATELEM(currentMat, k + 1, k)); /* make this entry zero */
subMatrix(currentMat, 1, k, 1, k, queue[queueL++]);
subMatrix(currentMat, k + 1, m, k + 1, m, queue[queueL++]);
doLoop = false;
}
else /* no deflation found yet */
{
mpTrafo(currentMat, it, tol2,R);
it++; /* try again */
}
}
if (doLoop) /* could not find deflation for currentMat */
{
deflationFound = false;
}
idDelete((ideal*)&currentMat);
}
}
return deflationFound;
}
absValue
number absValue(poly p)
Definition linearAlgebra.cc:734
hessenberg
void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat, const number tolerance, const ring R)
Computes the Hessenberg form of a given square matrix.
Definition linearAlgebra.cc:918
mpTrafo
void mpTrafo(matrix &H, int it, const number tolerance, const ring R)
Performs one transformation step on the given matrix H as part of the gouverning QR double shift algo...
Definition linearAlgebra.cc:988
quadraticSolve
int quadraticSolve(const poly p, number &s1, number &s2, const number tolerance)
Returns all solutions of a quadratic polynomial relation with real coefficients.
Definition linearAlgebra.cc:635
assume
#define assume(x)
Definition mod2.h:389
nGreater
#define nGreater(a, b)
Definition numbers.h:28

## quadraticSolve()

int quadraticSolve | ( | const poly | p,
|  | number & | s1,
|  | number & | s2,
|  | const number | tolerance )

Returns all solutions of a quadratic polynomial relation with real coefficients.

The method assumes that the polynomial is univariate in the first ring variable, and that the current ground field is the complex numbers. The polynomial has degree <= 2. Thus, there may be a) infinitely many zeros, when p == 0; then -1 is returned; b) no zero, when p == constant <> 0; then 0 is returned; c) one zero, when p is linear; then 1 is returned; d) one double zero; then 2 is returned; e) two distinct zeros; then 3 is returned; In the case e), s1 and s2 will contain the two distinct solutions. In cases c) and d) s1 will contain the single/double solution.

Returns
  the number of distinct zeros

Parameters

  [in] | p | the polynomial
  [out] | s1 | first zero, if any
  [out] | s2 | second zero, if any
  [in] | tolerance | accuracy

Definition at line 635 of file linearAlgebra.cc.

{
poly q = pCopy(p);
int result;
if (q == NULL) result = -1;
else
{
int degree = pGetExp(q, 1);
if (degree == 0) result = 0; /* constant polynomial <> 0 */
else
{
number c2 = nInit(0); /* coefficient of var(1)^2 */
number c1 = nInit(0); /* coefficient of var(1)^1 */
number c0 = nInit(0); /* coefficient of var(1)^0 */
if (pGetExp(q, 1) == 2)
{ nDelete(&c2); c2 = nCopy(pGetCoeff(q)); q = q->next; }
if ((q != NULL) && (pGetExp(q, 1) == 1))
{ nDelete(&c1); c1 = nCopy(pGetCoeff(q)); q = q->next; }
if ((q != NULL) && (pGetExp(q, 1) == 0))
{ nDelete(&c0); c0 = nCopy(pGetCoeff(q)); q = q->next; }
if (degree == 1)
{
c0 = nInpNeg(c0);
s1 = nDiv(c0, c1);
result = 1;
}
else
{
number tmp = nMult(c0, c2);
number tmp2 = nAdd(tmp, tmp); nDelete(&tmp);
number tmp4 = nAdd(tmp2, tmp2); nDelete(&tmp2);
number discr = nSub(nMult(c1, c1), tmp4); nDelete(&tmp4);
if (nIsZero(discr))
{
tmp = nAdd(c2, c2);
s1 = nDiv(c1, tmp); nDelete(&tmp);
s1 = nInpNeg(s1);
result = 2;
}
else if (nGreaterZero(discr))
{
realSqrt(discr, tolerance, tmp); /* sqrt of the discriminant */
tmp2 = nSub(tmp, c1);
tmp4 = nAdd(c2, c2);
s1 = nDiv(tmp2, tmp4); nDelete(&tmp2);
tmp = nInpNeg(tmp);
tmp2 = nSub(tmp, c1); nDelete(&tmp);
s2 = nDiv(tmp2, tmp4); nDelete(&tmp2); nDelete(&tmp4);
result = 3;
}
else
{
discr = nInpNeg(discr);
realSqrt(discr, tolerance, tmp); /* sqrt of |discriminant| */
tmp2 = nAdd(c2, c2);
tmp4 = nDiv(tmp, tmp2); nDelete(&tmp);
tmp = nDiv(c1, tmp2); nDelete(&tmp2);
tmp = nInpNeg(tmp);
s1 = (number)new gmp_complex(((gmp_complex*)tmp)->real(),
((gmp_complex*)tmp4)->real());
tmp4 = nInpNeg(tmp4);
s2 = (number)new gmp_complex(((gmp_complex*)tmp)->real(),
((gmp_complex*)tmp4)->real());
nDelete(&tmp); nDelete(&tmp4);
result = 3;
}
nDelete(&discr);
}
nDelete(&c0); nDelete(&c1); nDelete(&c2);
}
}
pDelete(&q);
return result;
}
degree
int degree(const CanonicalForm &f)
Definition canonicalform.h:316
gmp_complex
gmp_complex numbers based on
Definition mpr_complex.h:181
tmp2
CFList tmp2
Definition facFqBivar.cc:75
realSqrt
bool realSqrt(const number n, const number tolerance, number &root)
Computes the square root of a non-negative real number and returns it as a new number.
Definition linearAlgebra.cc:607
nGreaterZero
#define nGreaterZero(n)
Definition numbers.h:27

## realSqrt()

bool realSqrt | ( | const number | n,
|  | const number | tolerance,
|  | number & | root )

Computes the square root of a non-negative real number and returns it as a new number.

The method assumes that the current ground field is the complex numbers, and that the argument has imaginary part zero. If the real part is negative, then false is returned, otherwise true. The square root will be computed in the last argument by Heron's iteration formula with the first argument as the starting value. The iteration will stop as soon as two successive values (in the resulting sequence) differ by no more than the given tolerance, which is assumed to be a positive real number.

Returns
  the square root of the non-negative argument number

Parameters

  [in] | n | the input number
  [in] | tolerance | accuracy of iteration
  [out] | root | the root of n

Definition at line 607 of file linearAlgebra.cc.

{
if (!nGreaterZero(n)) return false;
if (nIsZero(n)) return nInit(0);
number oneHalf = complexNumber(0.5, 0.0);
number nHalf = nMult(n, oneHalf);
root = nCopy(n);
number nOld = complexNumber(10.0, 0.0);
number nDiff = nCopy(nOld);
while (nGreater(nDiff, tolerance))
{
nDelete(&nOld);
nOld = root;
number t1=nMult(oneHalf, nOld);
number t2=nDiv(nHalf, nOld);
root = nAdd(t1, t2);
nDelete(&t1);nDelete(&t2);
nDelete(&nDiff);
nDiff = nSub(nOld, root);
if (!nGreaterZero(nDiff)) nDiff = nInpNeg(nDiff);
}
nDelete(&nOld); nDelete(&nDiff); nDelete(&oneHalf); nDelete(&nHalf);
return true;
}
complexNumber
static number complexNumber(const double r, const double i)
Creates a new complex number from real and imaginary parts given by doubles.
Definition linearAlgebra.cc:545

## similar()

int similar | ( | const number * | nn,
|  | const int | nnLength,
|  | const number | n,
|  | const number | tolerance )

Tries to find the number n in the array nn[0..nnLength-1].

The method assumes that the ground field is the complex numbers. n and an entry of nn will be regarded equal when the absolute value of their difference is not larger than the given tolerance. In this case, the index of the respective entry of nn is returned, otherwise -1.

Returns
  the index of n in nn (up to tolerance) or -1

Parameters

  [in] | nn | array of numbers to look-up
  [in] | nnLength | length of nn
  [in] | n | number to loop-up in nn
  [in] | tolerance | tolerance for comparison

Definition at line 1197 of file linearAlgebra.cc.

{
int result = -1;
number tt = nMult(tolerance, tolerance);
number nr = (number)new gmp_complex(((gmp_complex*)n)->real());
number ni = (number)new gmp_complex(((gmp_complex*)n)->imag());
number rr; number ii;
number w1; number w2; number w3; number w4; number w5;
for (int i = 0; i < nnLength; i++)
{
rr = (number)new gmp_complex(((gmp_complex*)nn[i])->real());
ii = (number)new gmp_complex(((gmp_complex*)nn[i])->imag());
w1 = nSub(nr, rr); w2 = nMult(w1, w1);
w3 = nSub(ni, ii); w4 = nMult(w3, w3);
w5 = nAdd(w2, w4);
if (!nGreater(w5, tt)) result = i;
nDelete(&w1); nDelete(&w2); nDelete(&w3); nDelete(&w4);
nDelete(&w5); nDelete(&rr); nDelete(&ii);
if (result != -1) break;
}
nDelete(&tt); nDelete(&nr); nDelete(&ni);
return result;
}

## subMatrix()

bool subMatrix | ( | const matrix | aMat,
|  | const int | rowIndex1,
|  | const int | rowIndex2,
|  | const int | colIndex1,
|  | const int | colIndex2,
|  | matrix & | subMat )

Creates a new matrix which is a submatrix of the first argument, and returns true in case of success.

The method will be successful whenever rowIndex1 <= rowIndex2 and colIndex1 <= colIndex2. In this case and only then true will be returned and the last argument will afterwards contain a copy of the respective submatrix of the first argument. Note that this method may also be used to copy an entire matrix.

Returns
  true iff the submatrix could be build

Parameters

  [in] | aMat | the input matrix
  [in] | rowIndex1 | lower row index
  [in] | rowIndex2 | higher row index
  [in] | colIndex1 | lower column index
  [in] | colIndex2 | higher column index
  [out] | subMat | the new submatrix

Definition at line 742 of file linearAlgebra.cc.

{
if (rowIndex1 > rowIndex2) return false;
if (colIndex1 > colIndex2) return false;
int rr = rowIndex2 - rowIndex1 + 1;
int cc = colIndex2 - colIndex1 + 1;
subMat = mpNew(rr, cc);
for (int r = 1; r <= rr; r++)
for (int c = 1; c <= cc; c++)
MATELEM(subMat, r, c) =
pCopy(MATELEM(aMat, rowIndex1 + r - 1, colIndex1 + c - 1));
return true;
}

## swapColumns()

void swapColumns | ( | int | column1,
|  | int | column2,
|  | matrix & | aMat )

Swaps two columns of a given matrix and thereby modifies it.

The method expects two valid, distinct column indices, i.e. in [1..n], where n is the number of columns in aMat.

Parameters

  [in] | column1 | index of first column to swap
  [in] | column2 | index of second column to swap
  | aMat | [in/out] matrix subject to swapping

Definition at line 802 of file linearAlgebra.cc.

{
poly p;
int rr = MATROWS(aMat);
for (int r = 1; r <= rr; r++)
{
p = MATELEM(aMat, r, column1);
MATELEM(aMat, r, column1) = MATELEM(aMat, r, column2);
MATELEM(aMat, r, column2) = p;
}
}

## swapRows()

void swapRows | ( | int | row1,
|  | int | row2,
|  | matrix & | aMat )

Swaps two rows of a given matrix and thereby modifies it.

The method expects two valid, distinct row indices, i.e. in [1..n], where n is the number of rows in aMat.

Parameters

  [in] | row1 | index of first row to swap
  [in] | row2 | index of second row to swap
  | aMat | [in/out] matrix subject to swapping

Definition at line 790 of file linearAlgebra.cc.

{
poly p;
int cc = MATCOLS(aMat);
for (int c = 1; c <= cc; c++)
{
p = MATELEM(aMat, row1, c);
MATELEM(aMat, row1, c) = MATELEM(aMat, row2, c);
MATELEM(aMat, row2, c) = p;
}
}

## unitMatrix()

bool unitMatrix | ( | const int | n,
|  | matrix & | unitMat,
|  | const ring | r = currRing )

Creates a new matrix which is the (nxn) unit matrix, and returns true in case of success.

The method will be successful whenever n >= 1. In this case and only then true will be returned and the new (nxn) unit matrix will reside inside the second argument.

Returns
  true iff the (nxn) unit matrix could be build

Parameters

  [in] | n | size of the matrix
  [out] | unitMat | the new (nxn) unit matrix
  [in] | r | current ring

Definition at line 95 of file linearAlgebra.cc.

{
if (n < 1) return false;
unitMat = mpNew(n, n);
for (int r = 1; r <= n; r++) MATELEM(unitMat, r, r) = p_One(R);
return true;
}

## upperRightTriangleInverse()

bool upperRightTriangleInverse | ( | const matrix | uMat,
|  | matrix & | iMat,
|  | bool | diagonalIsOne,
|  | const ring | r = currRing )

Computes the inverse of a given (n x n)-matrix in upper right triangular form.

This method expects a quadratic matrix, that is, it must have as many rows as columns. Moreover, uMat[i, j] = 0, at least for all i > j, that is, u is in upper right triangular form.
If the argument diagonalIsOne is true, then we know additionally, that uMat[i, i] = 1, for all i. In this case uMat is invertible. Contrariwise, if diagonalIsOne is false, we do not know anything about the diagonal entries. (Note that they may still all be 1.)
In general, there are two cases:
1) The matrix is not invertible. Then the method returns false, and &iMat remains unchanged.
2) The matrix is invertible. Then the method returns true, and the content of iMat is the inverse of uMat.

Returns
  true iff uMat is invertible, false otherwise

Parameters

  [in] | uMat | (n x n)-matrix in upper right triangular form
  [out] | iMat | inverse of uMat if invertible
  [in] | diagonalIsOne | if true, then all diagonal entries of uMat are 1
  [in] | r | current ring

Definition at line 251 of file linearAlgebra.cc.

{
int d = uMat->rows();
poly p; poly q;
/* check whether uMat is invertible */
bool invertible = diagonalIsOne;
if (!invertible)
{
invertible = true;
for (int r = 1; r <= d; r++)
{
if (MATELEM(uMat, r, r) == NULL)
{
invertible = false;
break;
}
}
}
if (invertible)
{
iMat = mpNew(d, d);
for (int r = d; r >= 1; r--)
{
if (diagonalIsOne)
MATELEM(iMat, r, r) = p_One(R);
else
MATELEM(iMat, r, r) = p_NSet(n_Invers(p_GetCoeff(MATELEM(uMat, r, r),R),R->cf),R);
for (int c = r + 1; c <= d; c++)
{
p = NULL;
for (int k = r + 1; k <= c; k++)
{
q = pp_Mult_qq(MATELEM(uMat, r, k), MATELEM(iMat, k, c),R);
p = p_Add_q(p, q,R);
}
p = p_Neg(p,R);
p = p_Mult_q(p, p_Copy(MATELEM(iMat, r, r),R),R);
p_Normalize(p,R);
MATELEM(iMat, r, c) = p;
}
}
}
return invertible;
}
n_Invers
static FORCE_INLINE number n_Invers(number a, const coeffs r)
return the multiplicative inverse of 'a'; raise an error if 'a' is not invertible
Definition coeffs.h:567
p_GetCoeff
#define p_GetCoeff(p, r)
Definition monomials.h:50
p_Neg
static poly p_Neg(poly p, const ring r)
Definition p_polys.h:1114
p_Mult_q
static poly p_Mult_q(poly p, poly q, const ring r)
Definition p_polys.h:1125
pp_Mult_qq
static poly pp_Mult_qq(poly p, poly q, const ring r)
Definition p_polys.h:1167
p_Copy
static poly p_Copy(poly p, const ring r)
returns a copy of p
Definition p_polys.h:848