Function luDecomp

kernel/linear_algebra/linearAlgebra.cc:103–193

}

void luDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &uMat,
              const ring R)
{
  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 = new int[rr + 1];
  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);