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;
}

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);
}

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));
      }
    }
  }
}

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);
}

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;
}

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;
}

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;
}

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;
}

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;
  }
}

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;
}

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));
}

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 )

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;
}

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;
}

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;
}

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;
}

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;
}