Changeset 19bf86 in git

Timestamp:

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

Author:

Frank Seelisch <seelisch@…>

Branches:

(u'jengelh-datetime', 'ceac47cbc86fe4a15902392bdbb9bd2ae0ea02c6')(u'spielwiese', 'f875bbaccd0831e36aaed09ff6adeb3eb45aeb94')

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

Files:

• Singular/iparith.cc (modified) (3 diffs)
• Singular/tok.h (modified) (1 diff)
• doc/reference.doc (modified) (previous)
• kernel/linearAlgebra.cc (modified) (5 diffs)
• kernel/linearAlgebra.h (modified) (1 diff)

Singular/iparith.cc

@@ -317 +317 @@
   { "ludecomp",    0, LU_CMD ,            CMD_1},
   { "luinverse",   0, LUI_CMD ,           CMD_M},
+  { "lusolve",     0, LUS_CMD ,           CMD_M},
   { "map",         0, MAP_CMD ,           RING_DECL},
   { "matrix",      0, MATRIX_CMD ,        MATRIX_CMD},
@@ -7109 +7110 @@
   return FALSE;
 }
+static BOOLEAN jjLU_SOLVE(leftv res, leftv v)
+{
+  /* for solving a linear equation system A * x = b, via the
+     given LU-decomposition of the matrix A;
+     There is one valid parametrisation:
+     1) exactly four arguments P, L, U, b;
+        P, L, and U realise the L-U-decomposition of A, that is,
+        P * A = L * U, and P, L, and U satisfy the
+        properties decribed in method 'jjLU_DECOMP';
+        see there;
+        b is the right-hand side vector of the 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, dim] if there is at least one solution;
+        x is any solution, dim is the dimension of the affine solution
+        space.
+     The method produces an error if matrix and vector sizes do not fit. */
+  if ((v == NULL) || (v->Typ() != MATRIX_CMD) ||
+      (v->next == NULL) || (v->next->Typ() != MATRIX_CMD) ||
+      (v->next->next == NULL) || (v->next->next->Typ() != MATRIX_CMD) ||
+      (v->next->next->next == NULL) ||
+      (v->next->next->next->Typ() != MATRIX_CMD) ||
+      (v->next->next->next->next != NULL))
+  {
+    Werror("expected exactly three matrices and one vector as input");
+    return TRUE;
+  }
+  matrix pMat = (matrix)v->Data();
+  matrix lMat = (matrix)v->next->Data();
+  matrix uMat = (matrix)v->next->next->Data();
+  matrix bVec = (matrix)v->next->next->next->Data();
+  matrix xVec; int solvable; int dim;
+  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 (lMat->rows() != uMat->rows())
+  {
+    Werror("second matrix (%d x %d) and third matrix (%d x %d) do not fit",
+           lMat->rows(), lMat->cols(), uMat->rows(), uMat->cols());
+    return TRUE;
+  }
+  if (uMat->rows() != bVec->rows())
+  {
+    Werror("third matrix (%d x %d) and vector (%d x 1) do not fit",
+           uMat->rows(), uMat->cols(), bVec->rows());
+    return TRUE;
+  }
+  solvable = luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, &dim);
+
+  /* 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=INT_CMD;    ll->m[2].data=(void *)dim;
+  }
+  else
+  {
+    ll->Init(1);
+    ll->m[0].rtyp=INT_CMD;    ll->m[0].data=(void *)solvable;
+  }
+
+  res->data=(char*)ll;
+  return FALSE;
+}
 static BOOLEAN jjINTVEC_PL(leftv res, leftv v)
 {
@@ -7592 +7669 @@
 ,{jjLIST_PL,   LIST_CMD,        LIST_CMD,           -1      , ALLOW_PLURAL |ALLOW_RING}
 ,{jjLU_INVERSE,LUI_CMD,         LIST_CMD,           -2      , NO_PLURAL |NO_RING}
+,{jjLU_SOLVE,  LUS_CMD,         LIST_CMD,           -2      , NO_PLURAL |NO_RING}
 ,{jjWRONG,     MINOR_CMD,       NONE,               1       , ALLOW_PLURAL |ALLOW_RING}
 ,{jjMINOR_M,   MINOR_CMD,       IDEAL_CMD,          -2      , ALLOW_PLURAL |ALLOW_RING}

Singular/tok.h

@@ -97 +97 @@
   LU_CMD,
   LUI_CMD,
+  LUS_CMD,
   MEMORY_CMD,
   MONITOR_CMD,

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

kernel/linearAlgebra.h

@@ -198 +198 @@
        matrix &iMat       /**< [out] inverse of A if invertible   */
                           );
+                          
+/**
+ * Solves the linear system A*x = b, where A is an (n 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 proteries; see method
+ * 'luDecomp(const matrix aMat, matrix &pMat, matrix &lMat,
+ * matrix &uMat)'.<br>
+ * Instead of trying to invert A and return x = A^(-1)*b, this
+ * method
+ * 1) computes b' = pMat * b,
+ * 2) solves the easy system L * y = b', and then
+ * 3) solves the easy 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:<br>
+ * 1) There is no solution. Then the method returns false, and &xVec
+      as well as &dim remain unchanged.<br>
+ * 2) There is a unique solution. Then the method returns true and
+ *    the dimension of the affine solution space is zero.<br>
+ * 3) There are infinitely many solutions. Then the method returns
+ *    true and some solution. Furthermore, the dimension of the
+      affine solution space is set accordingly.
+ *
+ * @return true if there is at least one solution, false otherwise
+ **/
+bool luSolveViaLUDecomp(
+       const matrix pMat, /**< [in]  permutation matrix of an LU-
+                                     decomposition                */
+       const matrix lMat, /**< [in]  lower left matrix of an LU-
+                                     decomposition                */
+       const matrix uMat, /**< [in]  upper right matrix of an LU-
+                                     decomposition                */
+       const matrix bVec, /**< [in]  right-hand side of the linear
+                                     system to be solved          */
+       matrix &xVec,      /**< [out] solution of A*x = b          */
+       int* dim           /**< [out] dimension of affine solution
+                                     space                        */
+                          );
 
 #endif