# Changeset 19bf86 in git

Ignore:
Timestamp:
Jun 21, 2010, 1:20:59 PM (13 years ago)
Branches:
Children:
e3cc9c88085d8e217de6a163ac2bc7d891c19815
Parents:
e558c41b132228287898863201fbd5c35b4231fd
Message:
```implemented and doc'ed lusolve; available thru interpreter

Files:
5 edited

Unmodified
Removed
• ## Singular/iparith.cc

 re558c4 { "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}, 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) { ,{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

 re558c4 LU_CMD, LUI_CMD, LUS_CMD, MEMORY_CMD, MONITOR_CMD,
• ## kernel/linearAlgebra.cc

 re558c4 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 */ /* 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); 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) { /* 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), 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

 re558c4 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)'.
* 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:
* 1) There is no solution. Then the method returns false, and &xVec as well as &dim remain unchanged.
* 2) There is a unique solution. Then the method returns true and *    the dimension of the affine solution space is zero.
* 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
Note: See TracChangeset for help on using the changeset viewer.