* Performs one transformation step on the given matrix H as part of * the gouverning QR double shift algorithm. * The method will change the given matrix H side-effect-wise. The resulting * matrix H' will be in Hessenberg form. * The iteration index is needed, since for the 11th and 21st iteration, * the transformation step is different from the usual step, to avoid * convergence problems of
| 977 | * also the only caller of this method). |
| 978 | **/ |
| 979 | void mpTrafo( |
| 980 | matrix &H, /**< [in/out] the matrix to be transformed */ |
| 981 | int it, /**< [in] iteration index */ |
| 982 | const number tolerance,/**< [in] accuracy for square roots */ |
| 983 | const ring R |
| 984 | ) |
| 985 | { |
| 986 | int n = MATROWS(H); |
| 987 | number trace; number det; number tmp1; number tmp2; number tmp3; |
| 988 | |
| 989 | if ((it != 11) && (it != 21)) /* the standard case */ |
| 990 | { |
| 991 | /* in this case 'trace' will really be the trace of the lowermost |
| 992 | (2x2) block of hMat */ |
| 993 | trace = nInit(0); |
| 994 | det = nInit(0); |
| 995 | if (MATELEM(H, n - 1, n - 1) != NULL) |
| 996 | { |
| 997 | tmp1 = nAdd(trace, pGetCoeff(MATELEM(H, n - 1, n - 1))); |
| 998 | nDelete(&trace); |
| 999 | trace = tmp1; |
| 1000 | } |
| 1001 | if (MATELEM(H, n, n) != NULL) |
| 1002 | { |
| 1003 | tmp1 = nAdd(trace, pGetCoeff(MATELEM(H, n, n))); |
| 1004 | nDelete(&trace); |
| 1005 | trace = tmp1; |
| 1006 | } |
| 1007 | /* likewise 'det' will really be the determinante of the lowermost |
| 1008 | (2x2) block of hMat */ |
| 1009 | if ((MATELEM(H, n - 1, n - 1 ) != NULL) && (MATELEM(H, n, n) != NULL)) |
| 1010 | { |
| 1011 | tmp1 = nMult(pGetCoeff(MATELEM(H, n - 1, n - 1)), |
| 1012 | pGetCoeff(MATELEM(H, n, n))); |
| 1013 | tmp2 = nAdd(tmp1, det); nDelete(&tmp1); nDelete(&det); |
| 1014 | det = tmp2; |
| 1015 | } |
| 1016 | if ((MATELEM(H, n - 1, n) != NULL) && (MATELEM(H, n, n - 1) != NULL)) |
| 1017 | { |
| 1018 | tmp1 = nMult(pGetCoeff(MATELEM(H, n - 1, n)), |
| 1019 | pGetCoeff(MATELEM(H, n, n - 1))); |
| 1020 | tmp2 = nSub(det, tmp1); nDelete(&tmp1); nDelete(&det); |
| 1021 | det = tmp2; |
| 1022 | } |
| 1023 | } |
| 1024 | else |
| 1025 | { |
| 1026 | /* for it = 11 or it = 21, we use special formulae to avoid convergence |
| 1027 | problems of the gouverning QR double shift algorithm (who is the only |
| 1028 | caller of this method) */ |
| 1029 | /* trace = 3/2 * (|hMat[n, n-1]| + |hMat[n-1, n-2]|) */ |
| 1030 | tmp1 = nInit(0); |
| 1031 | if (MATELEM(H, n, n - 1) != NULL) |
| 1032 | { nDelete(&tmp1); tmp1 = nCopy(pGetCoeff(MATELEM(H, n, n - 1))); } |
| 1033 | if (!nGreaterZero(tmp1)) tmp1 = nInpNeg(tmp1); |
| 1034 | tmp2 = nInit(0); |
| 1035 | if (MATELEM(H, n - 1, n - 2) != NULL) |
| 1036 | { nDelete(&tmp2); tmp2 = nCopy(pGetCoeff(MATELEM(H, n - 1, n - 2))); } |
no test coverage detected