Changeset 19bf86 in git for kernel/linearAlgebra.cc

Timestamp:

Jun 21, 2010, 1:20:59 PM (14 years ago)

Author:

Frank Seelisch <seelisch@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

e3cc9c88085d8e217de6a163ac2bc7d891c19815

Parents:

e558c41b132228287898863201fbd5c35b4231fd

Message:

implemented and doc'ed lusolve; available thru interpreter

git-svn-id: file:///usr/local/Singular/svn/trunk@12885 2c84dea3-7e68-4137-9b89-c4e89433aadc

File:

• kernel/linearAlgebra.cc (modified) (5 diffs)

kernel/linearAlgebra.cc

@@ -103 +103 @@
   for (int r = 1; r <= rr; r++) MATELEM(lMat, r, r) = pOne();
 
-  int bestR; int bestC; int intSwap; poly pSwap;
+  int bestR; int bestC; int intSwap; poly pSwap; int cOffset = 0;
   for (int r = 1; r < rr; r++)
   {
-    if ((r <= cc) && (pivot(uMat, r, rr, r, r, &bestR, &bestC)))
+    if (r > cc) break;
+    while ((r + cOffset <= cc) &&
+           (!pivot(uMat, r, rr, r + cOffset, r + cOffset, &bestR, &bestC)))
+      cOffset++;
+    if (r + cOffset <= cc)
     {
       /* swap rows with indices r and bestR in permut */
@@ -114 +118 @@
 
       /* swap rows with indices r and bestR in uMat;
-         it is sufficient to do this for columns >= r */
-      for (int c = r; c <= cc; c++)
+         it is sufficient to do this for columns >= r + cOffset*/
+      for (int c = r + cOffset; c <= cc; c++)
       {
         pSwap = MATELEM(uMat, r, c);
@@ -134 +138 @@
          row with index r;
          we need to adjust lMat and uMat;
-         we are certain that the matrix entry at [r, r] is non-zero: */
-      number pivotElement = pGetCoeff(MATELEM(uMat, r, r));
+         we are certain that the matrix entry at [r, r + cOffset]
+         is non-zero: */
+      number pivotElement = pGetCoeff(MATELEM(uMat, r, r + cOffset));
       poly p; number n;
       for (int rGauss = r + 1; rGauss <= rr; rGauss++)
       {
-        p = MATELEM(uMat, rGauss, r);
+        p = MATELEM(uMat, rGauss, r + cOffset);
         if (p != NULL)
         {
@@ -149 +154 @@
 
           /* adjusting uMat: */
-          MATELEM(uMat, rGauss, r) = NULL; pDelete(&p);
+          MATELEM(uMat, rGauss, r + cOffset) = NULL; pDelete(&p);
           n = nNeg(n);
-          for (int cGauss = r + 1; cGauss <= cc; cGauss++)
+          for (int cGauss = r + cOffset + 1; cGauss <= cc; cGauss++)
             MATELEM(uMat, rGauss, cGauss)
               = pAdd(MATELEM(uMat, rGauss, cGauss),
@@ -316 +321 @@
   return result;
 }
+
+bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat,
+                        const matrix uMat, const matrix bVec,
+                        matrix &xVec, int* dim)
+{
+  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);
+  }
+  
+  /* 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); *dim = 0;
+    /* solve uMat * xVec = yVec and determine the dimension of the affine
+       solution space */
+    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;
+      *dim = *dim + lastNonZeroC - nonZeroC - 1;
+      poly 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)));
+      poly q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
+      MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
+      lastNonZeroC = nonZeroC;
+    }
+  }
+
+  idDelete((ideal*)&cVec);
+  idDelete((ideal*)&yVec);
+
+  return isSolvable;
+}