Index: Singular/extra.cc
===================================================================
--- Singular/extra.cc	(revision 2fd30cb3d4e41a93238b92df65a30fb88c8276de)
+++ Singular/extra.cc	(revision 9e269a0a45fc3fa26da7b381bc6593a5b4b7df4d)
@@ -667,4 +667,137 @@
         }
       }
+  /*==================== lduDecomp ======================*/
+      if(strcmp(sys_cmd, "lduDecomp")==0)
+      {
+        if ((h != NULL) && (h->Typ() == MATRIX_CMD) && (h->next == NULL))
+        {
+          matrix aMat = (matrix)h->Data();
+          matrix pMat; matrix lMat; matrix dMat; matrix uMat;
+          poly l; poly u; poly prodLU;
+          lduDecomp(aMat, pMat, lMat, dMat, uMat, l, u, prodLU);
+          lists L = (lists)omAllocBin(slists_bin);
+          L->Init(7);
+          L->m[0].rtyp = MATRIX_CMD; L->m[0].data=(void*)pMat;
+          L->m[1].rtyp = MATRIX_CMD; L->m[1].data=(void*)lMat;
+          L->m[2].rtyp = MATRIX_CMD; L->m[2].data=(void*)dMat;
+          L->m[3].rtyp = MATRIX_CMD; L->m[3].data=(void*)uMat;
+          L->m[4].rtyp = POLY_CMD; L->m[4].data=(void*)l;
+          L->m[5].rtyp = POLY_CMD; L->m[5].data=(void*)u;
+          L->m[6].rtyp = POLY_CMD; L->m[6].data=(void*)prodLU;
+          res->rtyp = LIST_CMD;
+          res->data = (char *)L;
+          return FALSE;
+        }
+        else
+        {
+          Werror( "expected argument list (int, int, poly, poly, poly, int)");
+          return TRUE;
+        }
+      }
+  /*==================== lduSolve ======================*/
+      if(strcmp(sys_cmd, "lduSolve")==0)
+      {
+        /* for solving a linear equation system A * x = b, via the
+           given LDU-decomposition of the matrix A;
+           There is one valid parametrisation:
+           1) exactly eight arguments P, L, D, U, l, u, lTimesU, b;
+              P, L, D, and U realise the LDU-decomposition of A, that is,
+              P * A = L * D^(-1) * U, and P, L, D, and U satisfy the
+              properties decribed in method 'luSolveViaLDUDecomp' in
+              linearAlgebra.h; see there;
+              l, u, and lTimesU are as described in the same location;
+              b is the right-hand side vector of the linear equation system;
+           The method will return a list of either 1 entry or three entries:
+           1) [0] if there is no solution to the system;
+           2) [1, x, H] if there is at least one solution;
+              x is any solution of the given linear system,
+              H is the matrix with column vectors spanning the homogeneous
+              solution space.
+           The method produces an error if matrix and vector sizes do not
+           fit. */
+        if ((h == NULL) || (h->Typ() != MATRIX_CMD) ||
+            (h->next == NULL) || (h->next->Typ() != MATRIX_CMD) ||
+            (h->next->next == NULL) || (h->next->next->Typ() != MATRIX_CMD) ||
+            (h->next->next->next == NULL) ||
+              (h->next->next->next->Typ() != MATRIX_CMD) ||
+            (h->next->next->next->next == NULL) ||
+              (h->next->next->next->next->Typ() != POLY_CMD) ||
+            (h->next->next->next->next->next == NULL) ||
+              (h->next->next->next->next->next->Typ() != POLY_CMD) ||
+            (h->next->next->next->next->next->next == NULL) ||
+              (h->next->next->next->next->next->next->Typ() != POLY_CMD) ||
+            (h->next->next->next->next->next->next->next == NULL) ||
+              (h->next->next->next->next->next->next->next->Typ()
+                != MATRIX_CMD) ||
+            (h->next->next->next->next->next->next->next->next != NULL))
+        {
+          Werror("expected input (matrix, matrix, matrix, matrix, %s",
+                                 "poly, poly, poly, matrix)");
+          return TRUE;
+        }
+        matrix pMat  = (matrix)h->Data();
+        matrix lMat  = (matrix)h->next->Data();
+        matrix dMat  = (matrix)h->next->next->Data();
+        matrix uMat  = (matrix)h->next->next->next->Data();
+        poly l       = (poly)  h->next->next->next->next->Data();
+        poly u       = (poly)  h->next->next->next->next->next->Data();
+        poly lTimesU = (poly)  h->next->next->next->next->next->next
+                                                              ->Data();
+        matrix bVec  = (matrix)h->next->next->next->next->next->next
+                                                        ->next->Data();
+        matrix xVec; int solvable; matrix homogSolSpace;
+        if (pMat->rows() != pMat->cols())
+        {
+          Werror("first matrix (%d x %d) is not quadratic",
+                 pMat->rows(), pMat->cols());
+          return TRUE;
+        }
+        if (lMat->rows() != lMat->cols())
+        {
+          Werror("second matrix (%d x %d) is not quadratic",
+                 lMat->rows(), lMat->cols());
+          return TRUE;
+        }
+        if (dMat->rows() != dMat->cols())
+        {
+          Werror("third matrix (%d x %d) is not quadratic",
+                 dMat->rows(), dMat->cols());
+          return TRUE;
+        }
+        if (dMat->cols() != uMat->rows())
+        {
+          Werror("third matrix (%d x %d) and fourth matrix (%d x %d) %s",
+                 "do not t",
+                 dMat->rows(), dMat->cols(), uMat->rows(), uMat->cols());
+          return TRUE;
+        }
+        if (uMat->rows() != bVec->rows())
+        {
+          Werror("fourth matrix (%d x %d) and vector (%d x 1) do not fit",
+                 uMat->rows(), uMat->cols(), bVec->rows());
+          return TRUE;
+        }
+        solvable = luSolveViaLDUDecomp(pMat, lMat, dMat, uMat, l, u, lTimesU,
+                                       bVec, xVec, homogSolSpace);
+
+        /* build the return structure; a list with either one or
+           three entries */
+        lists ll = (lists)omAllocBin(slists_bin);
+        if (solvable)
+        {
+          ll->Init(3);
+          ll->m[0].rtyp=INT_CMD;    ll->m[0].data=(void *)solvable;
+          ll->m[1].rtyp=MATRIX_CMD; ll->m[1].data=(void *)xVec;
+          ll->m[2].rtyp=MATRIX_CMD; ll->m[2].data=(void *)homogSolSpace;
+        }
+        else
+        {
+          ll->Init(1);
+          ll->m[0].rtyp=INT_CMD;    ll->m[0].data=(void *)solvable;
+        }
+        res->rtyp = LIST_CMD;
+        res->data=(char*)ll;
+        return FALSE;
+      }
   /*==================== forking experiments ======================*/
       if(strcmp(sys_cmd, "waitforssilinks")==0)
Index: Singular/iparith.cc
===================================================================
--- Singular/iparith.cc	(revision 2fd30cb3d4e41a93238b92df65a30fb88c8276de)
+++ Singular/iparith.cc	(revision 9e269a0a45fc3fa26da7b381bc6593a5b4b7df4d)
@@ -7462,5 +7462,5 @@
       (v->next->next->next->next != NULL))
   {
-    Werror("expected exactly three matrices and one vector as input");
+    WerrorS("expected exactly three matrices and one vector as input");
     return TRUE;
   }
Index: kernel/linearAlgebra.cc
===================================================================
--- kernel/linearAlgebra.cc	(revision 2fd30cb3d4e41a93238b92df65a30fb88c8276de)
+++ kernel/linearAlgebra.cc	(revision 9e269a0a45fc3fa26da7b381bc6593a5b4b7df4d)
@@ -98,5 +98,4 @@
   int cc = aMat->cols();
   pMat = mpNew(rr, rr);
-  lMat = mpNew(rr, rr);
   uMat = mpNew(rr, cc);
   /* copy aMat into uMat: */
@@ -440,4 +439,5 @@
         MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w));
       }
+      number z;
       for (int d = 1; d <= dim; d++)
       {
@@ -450,6 +450,12 @@
           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);;
+        MATELEM(N, nonZeroC, d) = pNeg(p);
+        if (MATELEM(N, nonZeroC, d) != NULL)
+        {
+          z = nDiv(pGetCoeff(MATELEM(N, nonZeroC, d)),
+                   pGetCoeff(MATELEM(uMat, r, nonZeroC)));
+          nNormalize(z); pDelete(&MATELEM(N, nonZeroC, d));
+          MATELEM(N, nonZeroC, d) = pNSet(z);
+        }
       }
       p = pNeg(pCopy(MATELEM(yVec, r, 1)));
@@ -457,7 +463,12 @@
         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));
+      MATELEM(xVec, nonZeroC, 1) = pNeg(p);
+      if (MATELEM(xVec, nonZeroC, 1) != NULL)
+      {
+        z = nDiv(pGetCoeff(MATELEM(xVec, nonZeroC, 1)),
+                 pGetCoeff(MATELEM(uMat, r, nonZeroC)));
+        nNormalize(z); pDelete(&MATELEM(xVec, nonZeroC, 1));
+        MATELEM(xVec, nonZeroC, 1) = pNSet(z);       
+      }
       lastNonZeroC = nonZeroC;
     }
@@ -491,7 +502,7 @@
 void printNumber(const number z)
 {
-	if (nIsZero(z)) printf("number = 0\n");
-	else
-	{
+  if (nIsZero(z)) printf("number = 0\n");
+  else
+  {
     poly p = pOne();
     pSetCoeff(p, nCopy(z));
@@ -1389,2 +1400,276 @@
 }
 
+void lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat,
+               matrix &uMat, poly &l, poly &u, poly &lTimesU)
+{
+  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 = nGcd(t, pGetCoeff(MATELEM(uMat, r, col)), currRing);
+          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++)
+          pMult_nn(MATELEM(dMat, r, r), t);
+        pMult_nn(MATELEM(lMat, row, row), t);
+      }
+      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 = nGcd(pGetCoeff(MATELEM(uMat, row, col)),
+                          pGetCoeff(MATELEM(uMat, r, col)), currRing);
+          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);
+            pMult_nn(p, f2);
+            poly q = pCopy(MATELEM(uMat, row, c));
+            pMult_nn(q, f1); q = pNeg(q);
+            MATELEM(uMat, r, c) = pAdd(p, q);
+          }
+          number tt = nDiv(g, gg);
+          nNormalize(tt);   /* this division works without remainder */
+          pMult_nn(MATELEM(lMat, r, r), tt); nDelete(&tt);
+          MATELEM(lMat, r, row) = pCopy(MATELEM(lMat, r, r));
+          pMult_nn(MATELEM(lMat, r, row), f1);
+          nDelete(&f1); nDelete(&f2); nDelete(&g);
+        }
+        else pMult_nn(MATELEM(lMat, r, r), t);
+      }
+      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);
+}
+
+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)
+{
+  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 */
+  number z;
+  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);
+    if (MATELEM(yVec, r, 1) != NULL)
+    {
+      z = nDiv(pGetCoeff(MATELEM(yVec, r, 1)), pGetCoeff(MATELEM(lMat, r, r)));
+      nNormalize(z);   /* division is without remainder */
+      pDelete(&MATELEM(yVec, r, 1));
+      MATELEM(yVec, r, 1) = pNSet(z);
+    }
+  }
+  
+  /* compute u * dMat * yVec and put result into yVec */
+  poly p;
+  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;
+    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; 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)));
+        MATELEM(N, nonZeroC, d) = pNeg(p);
+        if (MATELEM(N, nonZeroC, d) != NULL)
+        {
+          z = nDiv(pGetCoeff(MATELEM(N, nonZeroC, d)),
+                   pGetCoeff(MATELEM(uMat, r, nonZeroC)));
+          nNormalize(z);   /* division may be with remainder but only takes
+                              place for dim > 0 */
+          pDelete(&MATELEM(N, nonZeroC, d));
+          MATELEM(N, nonZeroC, d) = pNSet(z);
+        }
+      }
+      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)));
+      MATELEM(xVec, nonZeroC, 1) = pNeg(p);
+      if (MATELEM(xVec, nonZeroC, 1) != NULL)
+      {
+        z = nDiv(pGetCoeff(MATELEM(xVec, nonZeroC, 1)),
+                 pGetCoeff(MATELEM(uMat, r, nonZeroC)));
+        nNormalize(z);   /* division is without remainder */
+        pDelete(&MATELEM(xVec, nonZeroC, 1));
+        MATELEM(xVec, nonZeroC, 1) = pNSet(z);
+      }
+      lastNonZeroC = nonZeroC;
+    }
+    
+    /* divide xVec by l*u = lTimesU and put result in xVec */
+    number zz = pGetCoeff(lTimesU);
+    for (int c = 1; c <= n; c++)
+    {
+      if (MATELEM(xVec, c, 1) != NULL)
+      {
+        z = nDiv(pGetCoeff(MATELEM(xVec, c, 1)), zz);
+        nNormalize(z);
+        pDelete(&MATELEM(xVec, c, 1));
+        MATELEM(xVec, c, 1) = pNSet(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;
+}
+
Index: kernel/linearAlgebra.h
===================================================================
--- kernel/linearAlgebra.h	(revision 2fd30cb3d4e41a93238b92df65a30fb88c8276de)
+++ kernel/linearAlgebra.h	(revision 9e269a0a45fc3fa26da7b381bc6593a5b4b7df4d)
@@ -39,6 +39,7 @@
  *   entries equal to 1,<br>
  * - U is an (m x n) matrix in upper row echelon form.<br>
- * From these conditions, it easily follows that also A = P * L * U,
+ * From these conditions, it follows immediately that also A = P * L * U,
  * since P is self-inverse.<br>
+ *
  * The method will modify all argument matrices except aMat, so that
  * they will finally contain the matrices P, L, and U as specified
@@ -264,4 +265,65 @@
 
 /**
+ * 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 proteries; see method
+ * 'lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
+ * matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)'.<br>
+ * 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:<br>
+ * 1) There is no solution. Then the method returns false, and &xVec
+ *    as well as &H remain unchanged.<br>
+ * 2) There is a unique solution. Then the method returns true and
+ *    H is the 1x1 matrix with zero entry.<br>
+ * 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.
+ *
+ * @return true if there is at least one solution, false otherwise
+ **/
+bool luSolveViaLDUDecomp(
+       const matrix pMat,  /**< [in]  permutation matrix of an LDU-
+                                      decomposition                     */
+       const matrix lMat,  /**< [in]  lower left matrix of an LDU-
+                                      decomposition                     */
+       const matrix dMat,  /**< [in]  diagonal matrix of an LDU-
+                                      decomposition                     */
+       const matrix uMat,  /**< [in]  upper right matrix of an LDU-
+                                      decomposition                     */
+       const poly l,       /**< [in]  pivot product l of an LDU decomp. */
+       const poly u,       /**< [in]  pivot product u of an LDU decomp. */
+       const poly lTimesU, /**< [in]  product of l and u                */
+       const matrix bVec,  /**< [in]  right-hand side of the linear
+                                      system to be solved               */
+       matrix &xVec,       /**< [out] solution of A*x = b               */
+       matrix &H           /**< [out] matrix with columns spanning
+                                      homogeneous solution space        */
+                          );
+
+/**
  * Creates a new matrix which is the (nxn) unit matrix, and returns true
  * in case of success.
@@ -489,4 +551,42 @@
                               );
 
+/**
+ * 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<br>
+ * - P * A = L * D^(-1) * U,<br>
+ * - P is an (m x m) permutation matrix, i.e., its row/columns form the
+ *   standard basis of K^m,<br>
+ * - L is an (m x m) matrix in lower triangular form of full rank,<br>
+ * - D is an (m x m) diagonal matrix with full rank, and<br>
+ * - U is an (m x n) matrix in upper row echelon form.<br>
+ * From these conditions, it follows immediately that also
+ * A = P * L * D^(-1) * U, since P is self-inverse.<br>
+ *
+ * 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.
+ **/
+void lduDecomp(
+       const matrix aMat, /**< [in]  the initial matrix A,          */
+       matrix &pMat,      /**< [out] the row permutation matrix P   */
+       matrix &lMat,      /**< [out] the lower triangular matrix L  */
+       matrix &dMat,      /**< [out] the diagonal matrix D          */
+       matrix &uMat,      /**< [out] the upper row echelon matrix U */
+       poly &l,           /**< [out] the product of pivots of L     */
+       poly &u,           /**< [out] the product of pivots of U     */
+       poly &lTimesU      /**< [out] the product of l and u         */
+             );
+
 #endif
 /* LINEAR_ALGEBRA_H */
