Changeset 5d651d in git


Ignore:
Timestamp:
Sep 20, 2010, 12:55:44 PM (14 years ago)
Author:
Frank Seelisch <seelisch@…>
Branches:
(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b52fc4b2495505785981d640dcf7eb3e456778ef')
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:
4 edited

Legend:

Unmodified
Added
Removed
  • Singular/iparith.cc

    r6ea057 r5d651d  
    77427742     The method will return a list of either 1 entry or three entries:
    77437743     1) [0] if there is no solution to the system;
    7744      2) [1, x, dim] if there is at least one solution;
    7745         x is any solution, dim is the dimension of the affine solution
    7746         space.
     7744     2) [1, x, H] if there is at least one solution;
     7745        x is any solution of the given linear system,
     7746        H is the matrix with column vectors spanning the homogeneous
     7747        solution space.
    77477748     The method produces an error if matrix and vector sizes do not fit. */
    77487749  if ((v == NULL) || (v->Typ() != MATRIX_CMD) ||
     
    77607761  matrix uMat = (matrix)v->next->next->Data();
    77617762  matrix bVec = (matrix)v->next->next->next->Data();
    7762   matrix xVec; int solvable; int dim;
     7763  matrix xVec; int solvable; matrix homogSolSpace;
    77637764  if (pMat->rows() != pMat->cols())
    77647765  {
     
    77857786    return TRUE;
    77867787  }
    7787   solvable = luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, &dim);
     7788  solvable = luSolveViaLUDecomp(pMat, lMat, uMat, bVec, xVec, homogSolSpace);
    77887789
    77897790  /* build the return structure; a list with either one or three entries */
     
    77947795    ll->m[0].rtyp=INT_CMD;    ll->m[0].data=(void *)solvable;
    77957796    ll->m[1].rtyp=MATRIX_CMD; ll->m[1].data=(void *)xVec;
    7796     ll->m[2].rtyp=INT_CMD;    ll->m[2].data=(void *)dim;
     7797    ll->m[2].rtyp=MATRIX_CMD; ll->m[2].data=(void *)homogSolSpace;
    77977798  }
    77987799  else
  • kernel/linearAlgebra.cc

    r6ea057 r5d651d  
    363363bool luSolveViaLUDecomp(const matrix pMat, const matrix lMat,
    364364                        const matrix uMat, const matrix bVec,
    365                         matrix &xVec, int* dim)
     365                        matrix &xVec, matrix &H)
    366366{
    367367  int m = uMat->rows(); int n = uMat->cols();
     
    405405    else { nonZeroRowIndex = r; break; }
    406406  }
    407  
     407
    408408  if (isSolvable)
    409409  {
    410     xVec = mpNew(n, 1); *dim = 0;
    411     /* solve uMat * xVec = yVec and determine the dimension of the affine
    412        solution space */
     410    xVec = mpNew(n, 1);
     411    matrix N = mpNew(n, n); int dim = 0;
     412    poly p; poly q;
     413    /* solve uMat * xVec = yVec and determine a basis of the solution
     414       space of the homogeneous system uMat * xVec = 0;
     415       We do not know in advance what the dimension (dim) of the latter
     416       solution space will be. Thus, we start with the possibly too wide
     417       matrix N and later copy the relevant columns of N into H. */
    413418    int nonZeroC; int lastNonZeroC = n + 1;
    414419    for (int r = nonZeroRowIndex; r >= 1; r--)
     
    416421      for (nonZeroC = 1; nonZeroC <= n; nonZeroC++)
    417422        if (MATELEM(uMat, r, nonZeroC) != NULL) break;
    418       *dim = *dim + lastNonZeroC - nonZeroC - 1;
    419       poly p = pNeg(pCopy(MATELEM(yVec, r, 1)));
     423      for (int w = lastNonZeroC - 1; w >= nonZeroC + 1; w--)
     424      {
     425        /* this loop will only be done when the given linear system has
     426           more than one, i.e., infinitely many solutions */
     427        dim++;
     428        /* now we fill two entries of the dim-th column of N */
     429        MATELEM(N, w, dim) = pNeg(pCopy(MATELEM(uMat, r, nonZeroC)));
     430        MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w));
     431      }
     432      for (int d = 1; d <= dim; d++)
     433      {
     434        /* here we fill the entry of N at [r, d], for each valid vector
     435           that we already have in N;
     436           again, this will only be done when the given linear system has
     437           more than one, i.e., infinitely many solutions */
     438        p = NULL;
     439        for (int c = nonZeroC + 1; c <= n; c++)
     440          if (MATELEM(N, c, d) != NULL)
     441            p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d)));
     442        q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
     443        MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);;
     444      }
     445      p = pNeg(pCopy(MATELEM(yVec, r, 1)));
    420446      for (int c = nonZeroC + 1; c <= n; c++)
    421447        if (MATELEM(xVec, c, 1) != NULL)
    422448          p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1)));
    423       poly q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
     449      q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC))));
    424450      MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q);
    425451      pNormalize(MATELEM(xVec, nonZeroC, 1));
    426452      lastNonZeroC = nonZeroC;
    427453    }
     454    if (dim == 0)
     455    {
     456      /* that means the given linear system has exactly one solution;
     457         we return as H the 1x1 matrix with entry zero */
     458      H = mpNew(1, 1);
     459    }
     460    else
     461    {
     462      /* copy the first 'dim' columns of N into H */
     463      H = mpNew(n, dim);
     464      for (int r = 1; r <= n; r++)
     465        for (int c = 1; c <= dim; c++)
     466          MATELEM(H, r, c) = pCopy(MATELEM(N, r, c));
     467    }
     468    /* clean up N */
     469    idDelete((ideal*)&N);
    428470  }
    429471
  • kernel/linearAlgebra.h

    r6ea057 r5d651d  
    230230 * method
    231231 * 1) computes b' = pMat * b,
    232  * 2) solves the easy system L * y = b', and then
    233  * 3) solves the easy system U * x = y.
     232 * 2) solves the simple system L * y = b', and then
     233 * 3) solves the simple system U * x = y.
    234234 * Note that steps 1) and 2) will always work, as L is in any case
    235235 * invertible. Moreover, y is uniquely determined. Step 3) will only
     
    238238 * Thus, there are three cases:<br>
    239239 * 1) There is no solution. Then the method returns false, and &xVec
    240       as well as &dim remain unchanged.<br>
     240 *    as well as &H remain unchanged.<br>
    241241 * 2) There is a unique solution. Then the method returns true and
    242  *    the dimension of the affine solution space is zero.<br>
     242 *    H is the 1x1 matrix with zero entry.<br>
    243243 * 3) There are infinitely many solutions. Then the method returns
    244  *    true and some solution. Furthermore, the dimension of the
    245       affine solution space is set accordingly.
     244 *    true and some solution of the given original linear system.
     245 *    Furthermore, the columns of H span the solution space of the
     246 *    homogeneous linear system. The dimension of the solution
     247 *    space is then the number of columns of H.
    246248 *
    247249 * @return true if there is at least one solution, false otherwise
     
    257259                                     system to be solved          */
    258260       matrix &xVec,      /**< [out] solution of A*x = b          */
    259        int* dim           /**< [out] dimension of affine solution
    260                                      space                        */
     261       matrix &H          /**< [out] matrix with columns spanning
     262                                     homogeneous solution space   */
    261263                          );
    262264
    263265#endif
    264266/* LINEAR_ALGEBRA_H */
     267
Note: See TracChangeset for help on using the changeset viewer.