Changeset 9e269a0 in git for kernel/linearAlgebra.cc
 Timestamp:
 Dec 22, 2010, 1:54:29 PM (13 years ago)
 Branches:
 (u'spielwiese', 'd0474371d8c5d8068ab70bfb42719c97936b18a6')
 Children:
 e2fd5d453ccf5a7982b273cf039c5c7f41f41080
 Parents:
 2fd30cb3d4e41a93238b92df65a30fb88c8276de
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

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 nonzero 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 dimth 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
Note: See TracChangeset
for help on using the changeset viewer.