Changeset 9e269a0 in git
- Timestamp:
- Dec 22, 2010, 1:54:29 PM (12 years ago)
- Branches:
- (u'spielwiese', 'd1ba061a762c62d3a25159d8da8b6e17332291fa')
- Children:
- e2fd5d453ccf5a7982b273cf039c5c7f41f41080
- Parents:
- 2fd30cb3d4e41a93238b92df65a30fb88c8276de
- Files:
-
- 4 edited
-
Singular/extra.cc (modified) (1 diff)
-
Singular/iparith.cc (modified) (1 diff)
-
kernel/linearAlgebra.cc (modified) (6 diffs)
-
kernel/linearAlgebra.h (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
Singular/extra.cc
r2fd30cb r9e269a0 667 667 } 668 668 } 669 /*==================== lduDecomp ======================*/ 670 if(strcmp(sys_cmd, "lduDecomp")==0) 671 { 672 if ((h != NULL) && (h->Typ() == MATRIX_CMD) && (h->next == NULL)) 673 { 674 matrix aMat = (matrix)h->Data(); 675 matrix pMat; matrix lMat; matrix dMat; matrix uMat; 676 poly l; poly u; poly prodLU; 677 lduDecomp(aMat, pMat, lMat, dMat, uMat, l, u, prodLU); 678 lists L = (lists)omAllocBin(slists_bin); 679 L->Init(7); 680 L->m[0].rtyp = MATRIX_CMD; L->m[0].data=(void*)pMat; 681 L->m[1].rtyp = MATRIX_CMD; L->m[1].data=(void*)lMat; 682 L->m[2].rtyp = MATRIX_CMD; L->m[2].data=(void*)dMat; 683 L->m[3].rtyp = MATRIX_CMD; L->m[3].data=(void*)uMat; 684 L->m[4].rtyp = POLY_CMD; L->m[4].data=(void*)l; 685 L->m[5].rtyp = POLY_CMD; L->m[5].data=(void*)u; 686 L->m[6].rtyp = POLY_CMD; L->m[6].data=(void*)prodLU; 687 res->rtyp = LIST_CMD; 688 res->data = (char *)L; 689 return FALSE; 690 } 691 else 692 { 693 Werror( "expected argument list (int, int, poly, poly, poly, int)"); 694 return TRUE; 695 } 696 } 697 /*==================== lduSolve ======================*/ 698 if(strcmp(sys_cmd, "lduSolve")==0) 699 { 700 /* for solving a linear equation system A * x = b, via the 701 given LDU-decomposition of the matrix A; 702 There is one valid parametrisation: 703 1) exactly eight arguments P, L, D, U, l, u, lTimesU, b; 704 P, L, D, and U realise the LDU-decomposition of A, that is, 705 P * A = L * D^(-1) * U, and P, L, D, and U satisfy the 706 properties decribed in method 'luSolveViaLDUDecomp' in 707 linearAlgebra.h; see there; 708 l, u, and lTimesU are as described in the same location; 709 b is the right-hand side vector of the linear equation system; 710 The method will return a list of either 1 entry or three entries: 711 1) [0] if there is no solution to the system; 712 2) [1, x, H] if there is at least one solution; 713 x is any solution of the given linear system, 714 H is the matrix with column vectors spanning the homogeneous 715 solution space. 716 The method produces an error if matrix and vector sizes do not 717 fit. */ 718 if ((h == NULL) || (h->Typ() != MATRIX_CMD) || 719 (h->next == NULL) || (h->next->Typ() != MATRIX_CMD) || 720 (h->next->next == NULL) || (h->next->next->Typ() != MATRIX_CMD) || 721 (h->next->next->next == NULL) || 722 (h->next->next->next->Typ() != MATRIX_CMD) || 723 (h->next->next->next->next == NULL) || 724 (h->next->next->next->next->Typ() != POLY_CMD) || 725 (h->next->next->next->next->next == NULL) || 726 (h->next->next->next->next->next->Typ() != POLY_CMD) || 727 (h->next->next->next->next->next->next == NULL) || 728 (h->next->next->next->next->next->next->Typ() != POLY_CMD) || 729 (h->next->next->next->next->next->next->next == NULL) || 730 (h->next->next->next->next->next->next->next->Typ() 731 != MATRIX_CMD) || 732 (h->next->next->next->next->next->next->next->next != NULL)) 733 { 734 Werror("expected input (matrix, matrix, matrix, matrix, %s", 735 "poly, poly, poly, matrix)"); 736 return TRUE; 737 } 738 matrix pMat = (matrix)h->Data(); 739 matrix lMat = (matrix)h->next->Data(); 740 matrix dMat = (matrix)h->next->next->Data(); 741 matrix uMat = (matrix)h->next->next->next->Data(); 742 poly l = (poly) h->next->next->next->next->Data(); 743 poly u = (poly) h->next->next->next->next->next->Data(); 744 poly lTimesU = (poly) h->next->next->next->next->next->next 745 ->Data(); 746 matrix bVec = (matrix)h->next->next->next->next->next->next 747 ->next->Data(); 748 matrix xVec; int solvable; matrix homogSolSpace; 749 if (pMat->rows() != pMat->cols()) 750 { 751 Werror("first matrix (%d x %d) is not quadratic", 752 pMat->rows(), pMat->cols()); 753 return TRUE; 754 } 755 if (lMat->rows() != lMat->cols()) 756 { 757 Werror("second matrix (%d x %d) is not quadratic", 758 lMat->rows(), lMat->cols()); 759 return TRUE; 760 } 761 if (dMat->rows() != dMat->cols()) 762 { 763 Werror("third matrix (%d x %d) is not quadratic", 764 dMat->rows(), dMat->cols()); 765 return TRUE; 766 } 767 if (dMat->cols() != uMat->rows()) 768 { 769 Werror("third matrix (%d x %d) and fourth matrix (%d x %d) %s", 770 "do not t", 771 dMat->rows(), dMat->cols(), uMat->rows(), uMat->cols()); 772 return TRUE; 773 } 774 if (uMat->rows() != bVec->rows()) 775 { 776 Werror("fourth matrix (%d x %d) and vector (%d x 1) do not fit", 777 uMat->rows(), uMat->cols(), bVec->rows()); 778 return TRUE; 779 } 780 solvable = luSolveViaLDUDecomp(pMat, lMat, dMat, uMat, l, u, lTimesU, 781 bVec, xVec, homogSolSpace); 782 783 /* build the return structure; a list with either one or 784 three entries */ 785 lists ll = (lists)omAllocBin(slists_bin); 786 if (solvable) 787 { 788 ll->Init(3); 789 ll->m[0].rtyp=INT_CMD; ll->m[0].data=(void *)solvable; 790 ll->m[1].rtyp=MATRIX_CMD; ll->m[1].data=(void *)xVec; 791 ll->m[2].rtyp=MATRIX_CMD; ll->m[2].data=(void *)homogSolSpace; 792 } 793 else 794 { 795 ll->Init(1); 796 ll->m[0].rtyp=INT_CMD; ll->m[0].data=(void *)solvable; 797 } 798 res->rtyp = LIST_CMD; 799 res->data=(char*)ll; 800 return FALSE; 801 } 669 802 /*==================== forking experiments ======================*/ 670 803 if(strcmp(sys_cmd, "waitforssilinks")==0) -
Singular/iparith.cc
r2fd30cb r9e269a0 7462 7462 (v->next->next->next->next != NULL)) 7463 7463 { 7464 Werror ("expected exactly three matrices and one vector as input");7464 WerrorS("expected exactly three matrices and one vector as input"); 7465 7465 return TRUE; 7466 7466 } -
kernel/linearAlgebra.cc
r2fd30cb r9e269a0 98 98 int cc = aMat->cols(); 99 99 pMat = mpNew(rr, rr); 100 lMat = mpNew(rr, rr);101 100 uMat = mpNew(rr, cc); 102 101 /* copy aMat into uMat: */ … … 440 439 MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w)); 441 440 } 441 number z; 442 442 for (int d = 1; d <= dim; d++) 443 443 { … … 450 450 if (MATELEM(N, c, d) != NULL) 451 451 p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d))); 452 q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC)))); 453 MATELEM(N, nonZeroC, d) = pMult(pNeg(p), q);; 452 MATELEM(N, nonZeroC, d) = pNeg(p); 453 if (MATELEM(N, nonZeroC, d) != NULL) 454 { 455 z = nDiv(pGetCoeff(MATELEM(N, nonZeroC, d)), 456 pGetCoeff(MATELEM(uMat, r, nonZeroC))); 457 nNormalize(z); pDelete(&MATELEM(N, nonZeroC, d)); 458 MATELEM(N, nonZeroC, d) = pNSet(z); 459 } 454 460 } 455 461 p = pNeg(pCopy(MATELEM(yVec, r, 1))); … … 457 463 if (MATELEM(xVec, c, 1) != NULL) 458 464 p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1))); 459 q = pNSet(nInvers(pGetCoeff(MATELEM(uMat, r, nonZeroC)))); 460 MATELEM(xVec, nonZeroC, 1) = pMult(pNeg(p), q); 461 pNormalize(MATELEM(xVec, nonZeroC, 1)); 465 MATELEM(xVec, nonZeroC, 1) = pNeg(p); 466 if (MATELEM(xVec, nonZeroC, 1) != NULL) 467 { 468 z = nDiv(pGetCoeff(MATELEM(xVec, nonZeroC, 1)), 469 pGetCoeff(MATELEM(uMat, r, nonZeroC))); 470 nNormalize(z); pDelete(&MATELEM(xVec, nonZeroC, 1)); 471 MATELEM(xVec, nonZeroC, 1) = pNSet(z); 472 } 462 473 lastNonZeroC = nonZeroC; 463 474 } … … 491 502 void printNumber(const number z) 492 503 { 493 494 495 504 if (nIsZero(z)) printf("number = 0\n"); 505 else 506 { 496 507 poly p = pOne(); 497 508 pSetCoeff(p, nCopy(z)); … … 1389 1400 } 1390 1401 1402 void lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat, 1403 matrix &uMat, poly &l, poly &u, poly &lTimesU) 1404 { 1405 int rr = aMat->rows(); 1406 int cc = aMat->cols(); 1407 /* we use an int array to store all row permutations; 1408 note that we only make use of the entries [1..rr] */ 1409 int* permut = new int[rr + 1]; 1410 for (int i = 1; i <= rr; i++) permut[i] = i; 1411 /* fill lMat and dMat with the (rr x rr) unit matrix */ 1412 unitMatrix(rr, lMat); 1413 unitMatrix(rr, dMat); 1414 uMat = mpNew(rr, cc); 1415 /* copy aMat into uMat: */ 1416 for (int r = 1; r <= rr; r++) 1417 for (int c = 1; c <= cc; c++) 1418 MATELEM(uMat, r, c) = pCopy(MATELEM(aMat, r, c)); 1419 u = pOne(); l = pOne(); 1420 1421 int col = 1; int row = 1; 1422 while ((col <= cc) & (row < rr)) 1423 { 1424 int pivotR; int pivotC; bool pivotValid = false; 1425 while (col <= cc) 1426 { 1427 pivotValid = pivot(uMat, row, rr, col, col, &pivotR, &pivotC); 1428 if (pivotValid) break; 1429 col++; 1430 } 1431 if (pivotValid) 1432 { 1433 if (pivotR != row) 1434 { 1435 swapRows(row, pivotR, uMat); 1436 poly p = MATELEM(dMat, row, row); 1437 MATELEM(dMat, row, row) = MATELEM(dMat, pivotR, pivotR); 1438 MATELEM(dMat, pivotR, pivotR) = p; 1439 swapColumns(row, pivotR, lMat); 1440 swapRows(row, pivotR, lMat); 1441 int temp = permut[row]; 1442 permut[row] = permut[pivotR]; permut[pivotR] = temp; 1443 } 1444 /* in gg, we compute the gcd of all non-zero elements in 1445 uMat[row..rr, col]; 1446 the next number is the pivot and thus guaranteed to be different 1447 from zero: */ 1448 number gg = nCopy(pGetCoeff(MATELEM(uMat, row, col))); number t; 1449 for (int r = row + 1; r <= rr; r++) 1450 { 1451 if (MATELEM(uMat, r, col) != NULL) 1452 { 1453 t = gg; 1454 gg = nGcd(t, pGetCoeff(MATELEM(uMat, r, col)), currRing); 1455 nDelete(&t); 1456 } 1457 } 1458 t = nDiv(pGetCoeff(MATELEM(uMat, row, col)), gg); 1459 nNormalize(t); /* this division works without remainder */ 1460 if (!nIsOne(t)) 1461 { 1462 for (int r = row; r <= rr; r++) 1463 pMult_nn(MATELEM(dMat, r, r), t); 1464 pMult_nn(MATELEM(lMat, row, row), t); 1465 } 1466 l = pMult(l, pCopy(MATELEM(lMat, row, row))); 1467 u = pMult(u, pCopy(MATELEM(uMat, row, col))); 1468 1469 for (int r = row + 1; r <= rr; r++) 1470 { 1471 if (MATELEM(uMat, r, col) != NULL) 1472 { 1473 number g = nGcd(pGetCoeff(MATELEM(uMat, row, col)), 1474 pGetCoeff(MATELEM(uMat, r, col)), currRing); 1475 number f1 = nDiv(pGetCoeff(MATELEM(uMat, r, col)), g); 1476 nNormalize(f1); /* this division works without remainder */ 1477 number f2 = nDiv(pGetCoeff(MATELEM(uMat, row, col)), g); 1478 nNormalize(f2); /* this division works without remainder */ 1479 pDelete(&MATELEM(uMat, r, col)); MATELEM(uMat, r, col) = NULL; 1480 for (int c = col + 1; c <= cc; c++) 1481 { 1482 poly p = MATELEM(uMat, r, c); 1483 pMult_nn(p, f2); 1484 poly q = pCopy(MATELEM(uMat, row, c)); 1485 pMult_nn(q, f1); q = pNeg(q); 1486 MATELEM(uMat, r, c) = pAdd(p, q); 1487 } 1488 number tt = nDiv(g, gg); 1489 nNormalize(tt); /* this division works without remainder */ 1490 pMult_nn(MATELEM(lMat, r, r), tt); nDelete(&tt); 1491 MATELEM(lMat, r, row) = pCopy(MATELEM(lMat, r, r)); 1492 pMult_nn(MATELEM(lMat, r, row), f1); 1493 nDelete(&f1); nDelete(&f2); nDelete(&g); 1494 } 1495 else pMult_nn(MATELEM(lMat, r, r), t); 1496 } 1497 nDelete(&t); nDelete(&gg); 1498 col++; row++; 1499 } 1500 } 1501 /* one factor in the product u might be missing: */ 1502 if (row == rr) 1503 { 1504 while ((col <= cc) && (MATELEM(uMat, row, col) == NULL)) col++; 1505 if (col <= cc) u = pMult(u, pCopy(MATELEM(uMat, row, col))); 1506 } 1507 1508 /* building the permutation matrix from 'permut' and computing l */ 1509 pMat = mpNew(rr, rr); 1510 for (int r = 1; r <= rr; r++) 1511 MATELEM(pMat, r, permut[r]) = pOne(); 1512 delete[] permut; 1513 1514 lTimesU = ppMult_qq(l, u); 1515 } 1516 1517 bool luSolveViaLDUDecomp(const matrix pMat, const matrix lMat, 1518 const matrix dMat, const matrix uMat, 1519 const poly l, const poly u, const poly lTimesU, 1520 const matrix bVec, matrix &xVec, matrix &H) 1521 { 1522 int m = uMat->rows(); int n = uMat->cols(); 1523 matrix cVec = mpNew(m, 1); /* for storing l * pMat * bVec */ 1524 matrix yVec = mpNew(m, 1); /* for storing the unique solution of 1525 lMat * yVec = cVec */ 1526 1527 /* compute cVec = l * pMat * bVec but without actual matrix mult. */ 1528 for (int r = 1; r <= m; r++) 1529 { 1530 for (int c = 1; c <= m; c++) 1531 { 1532 if (MATELEM(pMat, r, c) != NULL) 1533 { MATELEM(cVec, r, 1) = ppMult_qq(l, MATELEM(bVec, c, 1)); break; } 1534 } 1535 } 1536 1537 /* solve lMat * yVec = cVec; this will always work as lMat is invertible; 1538 moreover, all divisions are guaranteed to be without remainder */ 1539 number z; 1540 for (int r = 1; r <= m; r++) 1541 { 1542 poly p = pNeg(pCopy(MATELEM(cVec, r, 1))); 1543 for (int c = 1; c < r; c++) 1544 p = pAdd(p, ppMult_qq(MATELEM(lMat, r, c), MATELEM(yVec, c, 1) )); 1545 MATELEM(yVec, r, 1) = pNeg(p); 1546 if (MATELEM(yVec, r, 1) != NULL) 1547 { 1548 z = nDiv(pGetCoeff(MATELEM(yVec, r, 1)), pGetCoeff(MATELEM(lMat, r, r))); 1549 nNormalize(z); /* division is without remainder */ 1550 pDelete(&MATELEM(yVec, r, 1)); 1551 MATELEM(yVec, r, 1) = pNSet(z); 1552 } 1553 } 1554 1555 /* compute u * dMat * yVec and put result into yVec */ 1556 poly p; 1557 for (int r = 1; r <= m; r++) 1558 { 1559 p = ppMult_qq(u, MATELEM(dMat, r, r)); 1560 MATELEM(yVec, r, 1) = pMult(p, MATELEM(yVec, r, 1)); 1561 } 1562 1563 /* determine whether uMat * xVec = yVec is solvable */ 1564 bool isSolvable = true; 1565 bool isZeroRow; int nonZeroRowIndex; 1566 for (int r = m; r >= 1; r--) 1567 { 1568 isZeroRow = true; 1569 for (int c = 1; c <= n; c++) 1570 if (MATELEM(uMat, r, c) != NULL) { isZeroRow = false; break; } 1571 if (isZeroRow) 1572 { 1573 if (MATELEM(yVec, r, 1) != NULL) { isSolvable = false; break; } 1574 } 1575 else { nonZeroRowIndex = r; break; } 1576 } 1577 1578 if (isSolvable) 1579 { 1580 xVec = mpNew(n, 1); 1581 matrix N = mpNew(n, n); int dim = 0; 1582 poly p; poly q; 1583 /* solve uMat * xVec = yVec and determine a basis of the solution 1584 space of the homogeneous system uMat * xVec = 0; 1585 We do not know in advance what the dimension (dim) of the latter 1586 solution space will be. Thus, we start with the possibly too wide 1587 matrix N and later copy the relevant columns of N into H. */ 1588 int nonZeroC; int lastNonZeroC = n + 1; 1589 for (int r = nonZeroRowIndex; r >= 1; r--) 1590 { 1591 for (nonZeroC = 1; nonZeroC <= n; nonZeroC++) 1592 if (MATELEM(uMat, r, nonZeroC) != NULL) break; 1593 for (int w = lastNonZeroC - 1; w >= nonZeroC + 1; w--) 1594 { 1595 /* this loop will only be done when the given linear system has 1596 more than one, i.e., infinitely many solutions */ 1597 dim++; 1598 /* now we fill two entries of the dim-th column of N */ 1599 MATELEM(N, w, dim) = pNeg(pCopy(MATELEM(uMat, r, nonZeroC))); 1600 MATELEM(N, nonZeroC, dim) = pCopy(MATELEM(uMat, r, w)); 1601 } 1602 for (int d = 1; d <= dim; d++) 1603 { 1604 /* here we fill the entry of N at [r, d], for each valid vector 1605 that we already have in N; 1606 again, this will only be done when the given linear system has 1607 more than one, i.e., infinitely many solutions */ 1608 p = NULL; 1609 for (int c = nonZeroC + 1; c <= n; c++) 1610 if (MATELEM(N, c, d) != NULL) 1611 p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(N, c, d))); 1612 MATELEM(N, nonZeroC, d) = pNeg(p); 1613 if (MATELEM(N, nonZeroC, d) != NULL) 1614 { 1615 z = nDiv(pGetCoeff(MATELEM(N, nonZeroC, d)), 1616 pGetCoeff(MATELEM(uMat, r, nonZeroC))); 1617 nNormalize(z); /* division may be with remainder but only takes 1618 place for dim > 0 */ 1619 pDelete(&MATELEM(N, nonZeroC, d)); 1620 MATELEM(N, nonZeroC, d) = pNSet(z); 1621 } 1622 } 1623 p = pNeg(pCopy(MATELEM(yVec, r, 1))); 1624 for (int c = nonZeroC + 1; c <= n; c++) 1625 if (MATELEM(xVec, c, 1) != NULL) 1626 p = pAdd(p, ppMult_qq(MATELEM(uMat, r, c), MATELEM(xVec, c, 1))); 1627 MATELEM(xVec, nonZeroC, 1) = pNeg(p); 1628 if (MATELEM(xVec, nonZeroC, 1) != NULL) 1629 { 1630 z = nDiv(pGetCoeff(MATELEM(xVec, nonZeroC, 1)), 1631 pGetCoeff(MATELEM(uMat, r, nonZeroC))); 1632 nNormalize(z); /* division is without remainder */ 1633 pDelete(&MATELEM(xVec, nonZeroC, 1)); 1634 MATELEM(xVec, nonZeroC, 1) = pNSet(z); 1635 } 1636 lastNonZeroC = nonZeroC; 1637 } 1638 1639 /* divide xVec by l*u = lTimesU and put result in xVec */ 1640 number zz = pGetCoeff(lTimesU); 1641 for (int c = 1; c <= n; c++) 1642 { 1643 if (MATELEM(xVec, c, 1) != NULL) 1644 { 1645 z = nDiv(pGetCoeff(MATELEM(xVec, c, 1)), zz); 1646 nNormalize(z); 1647 pDelete(&MATELEM(xVec, c, 1)); 1648 MATELEM(xVec, c, 1) = pNSet(z); 1649 } 1650 } 1651 1652 if (dim == 0) 1653 { 1654 /* that means the given linear system has exactly one solution; 1655 we return as H the 1x1 matrix with entry zero */ 1656 H = mpNew(1, 1); 1657 } 1658 else 1659 { 1660 /* copy the first 'dim' columns of N into H */ 1661 H = mpNew(n, dim); 1662 for (int r = 1; r <= n; r++) 1663 for (int c = 1; c <= dim; c++) 1664 MATELEM(H, r, c) = pCopy(MATELEM(N, r, c)); 1665 } 1666 /* clean up N */ 1667 idDelete((ideal*)&N); 1668 } 1669 1670 idDelete((ideal*)&cVec); 1671 idDelete((ideal*)&yVec); 1672 1673 return isSolvable; 1674 } 1675 -
kernel/linearAlgebra.h
r2fd30cb r9e269a0 39 39 * entries equal to 1,<br> 40 40 * - U is an (m x n) matrix in upper row echelon form.<br> 41 * From these conditions, it easily followsthat also A = P * L * U,41 * From these conditions, it follows immediately that also A = P * L * U, 42 42 * since P is self-inverse.<br> 43 * 43 44 * The method will modify all argument matrices except aMat, so that 44 45 * they will finally contain the matrices P, L, and U as specified … … 264 265 265 266 /** 267 * Solves the linear system A * x = b, where A is an (m x n)-matrix 268 * which is given by its LDU-decomposition. 269 * 270 * The method expects the LDU-decomposition of A, that is, 271 * pMat * A = lMat * dMat^(-1) * uMat, where the argument matrices have 272 * the appropriate proteries; see method 273 * 'lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, 274 * matrix &dMat, matrix &uMat, poly &l, poly &u, poly &lTimesU)'.<br> 275 * Instead of trying to invert A and return x = A^(-1)*b, this 276 * method 277 * 1) computes b' = l * pMat * b, 278 * 2) solves the simple system L * y = b', 279 * 3) computes y' = u * dMat * y, 280 * 4) solves the simple system U * y'' = y', 281 * 5) computes x = 1/(lTimesU) * y''. 282 * Note that steps 1), 2) and 3) will always work, as L is in any case 283 * invertible. Moreover, y and thus y' are uniquely determined. 284 * Step 4) will only work if and only if y' is in the column span of U. 285 * In that case, there may be infinitely many solutions. 286 * In contrast to luSolveViaLUDecomp, this algorithm guarantees that 287 * all divisions which have to be performed in steps 2) and 4) will 288 * work without remainder. This is due to properties of the given LDU- 289 * decomposition. Only the final step 5) performs a division of a vector 290 * by a member of the ground field. (Suppose the ground field is Q or 291 * some rational function field. Then, if A contains only integers or 292 * polynomials, respectively, then steps 1) - 4) will keep all data 293 * integer or polynomial, respectively. This may speed up computations 294 * considerably.) 295 * For the outcome, there are three cases:<br> 296 * 1) There is no solution. Then the method returns false, and &xVec 297 * as well as &H remain unchanged.<br> 298 * 2) There is a unique solution. Then the method returns true and 299 * H is the 1x1 matrix with zero entry.<br> 300 * 3) There are infinitely many solutions. Then the method returns 301 * true and some solution of the given original linear system. 302 * Furthermore, the columns of H span the solution space of the 303 * homogeneous linear system. The dimension of the solution 304 * space is then the number of columns of H. 305 * 306 * @return true if there is at least one solution, false otherwise 307 **/ 308 bool luSolveViaLDUDecomp( 309 const matrix pMat, /**< [in] permutation matrix of an LDU- 310 decomposition */ 311 const matrix lMat, /**< [in] lower left matrix of an LDU- 312 decomposition */ 313 const matrix dMat, /**< [in] diagonal matrix of an LDU- 314 decomposition */ 315 const matrix uMat, /**< [in] upper right matrix of an LDU- 316 decomposition */ 317 const poly l, /**< [in] pivot product l of an LDU decomp. */ 318 const poly u, /**< [in] pivot product u of an LDU decomp. */ 319 const poly lTimesU, /**< [in] product of l and u */ 320 const matrix bVec, /**< [in] right-hand side of the linear 321 system to be solved */ 322 matrix &xVec, /**< [out] solution of A*x = b */ 323 matrix &H /**< [out] matrix with columns spanning 324 homogeneous solution space */ 325 ); 326 327 /** 266 328 * Creates a new matrix which is the (nxn) unit matrix, and returns true 267 329 * in case of success. … … 489 551 ); 490 552 553 /** 554 * LU-decomposition of a given (m x n)-matrix with performing only those 555 * divisions that yield zero remainders. 556 * 557 * Given an (m x n) matrix A, the method computes its LDU-decomposition, 558 * that is, it computes matrices P, L, D, and U such that<br> 559 * - P * A = L * D^(-1) * U,<br> 560 * - P is an (m x m) permutation matrix, i.e., its row/columns form the 561 * standard basis of K^m,<br> 562 * - L is an (m x m) matrix in lower triangular form of full rank,<br> 563 * - D is an (m x m) diagonal matrix with full rank, and<br> 564 * - U is an (m x n) matrix in upper row echelon form.<br> 565 * From these conditions, it follows immediately that also 566 * A = P * L * D^(-1) * U, since P is self-inverse.<br> 567 * 568 * In contrast to luDecomp, this method only performs those divisions which 569 * yield zero remainders. Hence, when e.g. computing over a rational function 570 * field and starting with polynomial entries only (or over Q and starting 571 * with integer entries only), then any newly computed matrix entry will again 572 * be polynomial (or integer). 573 * 574 * The method will modify all argument matrices except aMat, so that 575 * they will finally contain the matrices P, L, D, and U as specified 576 * above. Moreover, in order to further speed up subsequent calls of 577 * luSolveViaLDUDecomp, the two denominators and their product will also be 578 * returned. 579 **/ 580 void lduDecomp( 581 const matrix aMat, /**< [in] the initial matrix A, */ 582 matrix &pMat, /**< [out] the row permutation matrix P */ 583 matrix &lMat, /**< [out] the lower triangular matrix L */ 584 matrix &dMat, /**< [out] the diagonal matrix D */ 585 matrix &uMat, /**< [out] the upper row echelon matrix U */ 586 poly &l, /**< [out] the product of pivots of L */ 587 poly &u, /**< [out] the product of pivots of U */ 588 poly &lTimesU /**< [out] the product of l and u */ 589 ); 590 491 591 #endif 492 592 /* LINEAR_ALGEBRA_H */
Note: See TracChangeset
for help on using the changeset viewer.
