Changeset 5d651d in git

Timestamp:

Sep 20, 2010, 12:55:44 PM (14 years ago)

Author:

Frank Seelisch <seelisch@…>

Branches:

(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', '38dfc5131670d387a89455159ed1e071997eec94')

Children:

b5f276e5ad58127bf983ec9cc30196b7e764d519

Parents:

6ea057a63c210b915ec61cedb5bc5fac0c8f1432

Message:

changed lusolve (incl. docu): computes now also basis of homogeneous linear system (instead of just dim as before)

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

Files:

• Singular/iparith.cc (modified) (4 diffs)
• doc/reference.doc (modified) (previous)
• kernel/linearAlgebra.cc (modified) (3 diffs)
• kernel/linearAlgebra.h (modified) (3 diffs)

Singular/iparith.cc

@@ -7742 +7742 @@
      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.
+     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 ((v == NULL) || (v->Typ() != MATRIX_CMD) ||
@@ -7760 +7761 @@
   matrix uMat = (matrix)v->next->next->Data();
   matrix bVec = (matrix)v->next->next->next->Data();
-  matrix xVec; int solvable; int dim;
+  matrix xVec; int solvable; matrix homogSolSpace;
   if (pMat->rows() != pMat->cols())
   {
@@ -7785 +7786 @@
     return TRUE;
   }
-  solvable = luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, &dim);
+  solvable = luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, homogSolSpace);
 
   /* build the return structure; a list with either one or three entries */
@@ -7794 +7795 @@
     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;
+    ll->m[2].rtyp=MATRIX_CMD; ll->m[2].data=(void *)homogSolSpace;
   }
   else

kernel/linearAlgebra.cc

@@ -363 +363 @@
 bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat,
                         const matrix uMat, const matrix bVec,
-                        matrix &xVec, int* dim)
+                        matrix &xVec, matrix &H)
 {
   int m = uMat->rows(); int n = uMat->cols();
@@ -405 +405 @@
     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 */
+    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--)
@@ -416 +421 @@
       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 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)));
+        q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
+        MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);;
+      }
+      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))));
+      q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
       MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
       pNormalize(MATELEM(xVec, nonZeroC, 1));
       lastNonZeroC = nonZeroC;
     }
+    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);
   }
 

kernel/linearAlgebra.h

@@ -230 +230 @@
  * method
  * 1) computes b' = pMat * b,
- * 2) solves the easy system L * y = b', and then
- * 3) solves the easy system U * x = y.
+ * 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
@@ -238 +238 @@
  * 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>
+ *    as well as &H 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>
+ *    H is the 1x1 matrix with zero entry.<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.
+ *    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
@@ -257 +259 @@
                                      system to be solved          */
        matrix &xVec,      /**< [out] solution of A*x = b          */
-       int* dim           /**< [out] dimension of affine solution
-                                     space                        */
+       matrix &H          /**< [out] matrix with columns spanning
+                                     homogeneous solution space   */
                           );
 
 #endif
 /* LINEAR_ALGEBRA_H */
+