| 1341 | } |
| 1342 | |
| 1343 | void lduDecomp(const matrix aMat, matrix &pMat, matrix &lMat, matrix &dMat, |
| 1344 | matrix &uMat, poly &l, poly &u, poly &lTimesU) |
| 1345 | { |
| 1346 | int rr = aMat->rows(); |
| 1347 | int cc = aMat->cols(); |
| 1348 | /* we use an int array to store all row permutations; |
| 1349 | note that we only make use of the entries [1..rr] */ |
| 1350 | int* permut = new int[rr + 1]; |
| 1351 | for (int i = 1; i <= rr; i++) permut[i] = i; |
| 1352 | /* fill lMat and dMat with the (rr x rr) unit matrix */ |
| 1353 | unitMatrix(rr, lMat); |
| 1354 | unitMatrix(rr, dMat); |
| 1355 | uMat = mpNew(rr, cc); |
| 1356 | /* copy aMat into uMat: */ |
| 1357 | for (int r = 1; r <= rr; r++) |
| 1358 | for (int c = 1; c <= cc; c++) |
| 1359 | MATELEM(uMat, r, c) = pCopy(MATELEM(aMat, r, c)); |
| 1360 | u = pOne(); l = pOne(); |
| 1361 | |
| 1362 | int col = 1; int row = 1; |
| 1363 | while ((col <= cc) & (row < rr)) |
| 1364 | { |
| 1365 | int pivotR; int pivotC; bool pivotValid = false; |
| 1366 | while (col <= cc) |
| 1367 | { |
| 1368 | pivotValid = pivot(uMat, row, rr, col, col, &pivotR, &pivotC); |
| 1369 | if (pivotValid) break; |
| 1370 | col++; |
| 1371 | } |
| 1372 | if (pivotValid) |
| 1373 | { |
| 1374 | if (pivotR != row) |
| 1375 | { |
| 1376 | swapRows(row, pivotR, uMat); |
| 1377 | poly p = MATELEM(dMat, row, row); |
| 1378 | MATELEM(dMat, row, row) = MATELEM(dMat, pivotR, pivotR); |
| 1379 | MATELEM(dMat, pivotR, pivotR) = p; |
| 1380 | swapColumns(row, pivotR, lMat); |
| 1381 | swapRows(row, pivotR, lMat); |
| 1382 | int temp = permut[row]; |
| 1383 | permut[row] = permut[pivotR]; permut[pivotR] = temp; |
| 1384 | } |
| 1385 | /* in gg, we compute the gcd of all non-zero elements in |
| 1386 | uMat[row..rr, col]; |
| 1387 | the next number is the pivot and thus guaranteed to be different |
| 1388 | from zero: */ |
| 1389 | number gg = nCopy(pGetCoeff(MATELEM(uMat, row, col))); number t; |
| 1390 | for (int r = row + 1; r <= rr; r++) |
| 1391 | { |
| 1392 | if (MATELEM(uMat, r, col) != NULL) |
| 1393 | { |
| 1394 | t = gg; |
| 1395 | gg = n_Gcd(t, pGetCoeff(MATELEM(uMat, r, col)),currRing->cf); |
| 1396 | nDelete(&t); |
| 1397 | } |
| 1398 | } |
| 1399 | t = nDiv(pGetCoeff(MATELEM(uMat, row, col)), gg); |
| 1400 | nNormalize(t); /* this division works without remainder */ |
no test coverage detected