helper for qrDoubleShift */
| 1088 | |
| 1089 | /* helper for qrDoubleShift */ |
| 1090 | bool qrDS( |
| 1091 | const int /*n*/, |
| 1092 | matrix* queue, |
| 1093 | int& queueL, |
| 1094 | number* eigenValues, |
| 1095 | int& eigenValuesL, |
| 1096 | const number tol1, |
| 1097 | const number tol2, |
| 1098 | const ring R |
| 1099 | ) |
| 1100 | { |
| 1101 | bool deflationFound = true; |
| 1102 | /* we loop until the working queue is empty, |
| 1103 | provided we always find deflation */ |
| 1104 | while (deflationFound && (queueL > 0)) |
| 1105 | { |
| 1106 | /* take out last queue entry */ |
| 1107 | matrix currentMat = queue[queueL - 1]; queueL--; |
| 1108 | int m = MATROWS(currentMat); |
| 1109 | if (m == 1) |
| 1110 | { |
| 1111 | number newEigenvalue; |
| 1112 | /* the entry at [1, 1] is the eigenvalue */ |
| 1113 | if (MATELEM(currentMat, 1, 1) == NULL) newEigenvalue = nInit(0); |
| 1114 | else newEigenvalue = nCopy(pGetCoeff(MATELEM(currentMat, 1, 1))); |
| 1115 | eigenValues[eigenValuesL++] = newEigenvalue; |
| 1116 | } |
| 1117 | else if (m == 2) |
| 1118 | { |
| 1119 | /* there are two eigenvalues which come as zeros of the characteristic |
| 1120 | polynomial */ |
| 1121 | poly p; charPoly(currentMat, p); |
| 1122 | number s1; number s2; |
| 1123 | int nSol = quadraticSolve(p, s1, s2, tol2); pDelete(&p); |
| 1124 | assume(nSol >= 2); |
| 1125 | eigenValues[eigenValuesL++] = s1; |
| 1126 | /* if nSol = 2, then s1 is a double zero, and s2 is invalid: */ |
| 1127 | if (nSol == 2) s2 = nCopy(s1); |
| 1128 | eigenValues[eigenValuesL++] = s2; |
| 1129 | } |
| 1130 | else /* m > 2 */ |
| 1131 | { |
| 1132 | /* bring currentMat into Hessenberg form to fasten computations: */ |
| 1133 | matrix mm1; matrix mm2; |
| 1134 | hessenberg(currentMat, mm1, mm2, tol2,R); |
| 1135 | idDelete((ideal*)¤tMat); idDelete((ideal*)&mm1); |
| 1136 | currentMat = mm2; |
| 1137 | int it = 1; bool doLoop = true; |
| 1138 | while (doLoop && (it <= 30 * m)) |
| 1139 | { |
| 1140 | /* search for deflation */ |
| 1141 | number w1; number w2; |
| 1142 | number test1; number test2; bool stopCriterion = false; int k; |
| 1143 | for (k = 1; k < m; k++) |
| 1144 | { |
| 1145 | test1 = absValue(MATELEM(currentMat, k + 1, k)); |
| 1146 | w1 = absValue(MATELEM(currentMat, k, k)); |
| 1147 | w2 = absValue(MATELEM(currentMat, k + 1, k + 1)); |
no test coverage detected