| 907 | } |
| 908 | |
| 909 | void hessenberg(const matrix aMat, matrix &pMat, matrix &hessenbergMat, |
| 910 | const number tolerance, const ring R) |
| 911 | { |
| 912 | int n = MATROWS(aMat); |
| 913 | unitMatrix(n, pMat); |
| 914 | subMatrix(aMat, 1, n, 1, n, hessenbergMat); |
| 915 | for (int c = 1; c <= n; c++) |
| 916 | { |
| 917 | /* find one or two non-zero entries in the current column */ |
| 918 | int r1 = 0; int r2 = 0; |
| 919 | for (int r = c + 1; r <= n; r++) |
| 920 | if (MATELEM(hessenbergMat, r, c) != NULL) |
| 921 | { |
| 922 | if (r1 == 0) r1 = r; |
| 923 | else if (r2 == 0) { r2 = r; break; } |
| 924 | } |
| 925 | if (r1 != 0) |
| 926 | { /* At least one entry in the current column is non-zero. */ |
| 927 | if (r1 != c + 1) |
| 928 | { /* swap rows to bring non-zero element to row with index c + 1 */ |
| 929 | swapRows(r1, c + 1, hessenbergMat); |
| 930 | /* now also swap columns to reflect action of permutation |
| 931 | from the right-hand side */ |
| 932 | swapColumns(r1, c + 1, hessenbergMat); |
| 933 | /* include action of permutation also in pMat */ |
| 934 | swapRows(r1, c + 1, pMat); |
| 935 | } |
| 936 | if (r2 != 0) |
| 937 | { /* There is at least one more non-zero element in the current |
| 938 | column. So let us perform a hessenberg step in order to make |
| 939 | all additional non-zero elements zero. */ |
| 940 | matrix v; subMatrix(hessenbergMat, c + 1, n, c, c, v); |
| 941 | matrix u; matrix pTmp; |
| 942 | number r = hessenbergStep(v, u, pTmp, tolerance); |
| 943 | idDelete((ideal*)&v); idDelete((ideal*)&u); nDelete(&r); |
| 944 | /* pTmp just acts on the lower right block of hessenbergMat; |
| 945 | i.e., it needs to be extended by a unit matrix block at the top |
| 946 | left in order to define a whole transformation matrix; |
| 947 | this code may be optimized */ |
| 948 | unitMatrix(c, u); |
| 949 | matrix pTmpFull; matrixBlock(u, pTmp, pTmpFull); |
| 950 | idDelete((ideal*)&u); idDelete((ideal*)&pTmp); |
| 951 | /* now include pTmpFull in pMat (by letting it act from the left) */ |
| 952 | pTmp = mp_Mult(pTmpFull, pMat,R); idDelete((ideal*)&pMat); |
| 953 | pMat = pTmp; |
| 954 | /* now let pTmpFull act on hessenbergMat from the left and from the |
| 955 | right (note that pTmpFull is self-inverse) */ |
| 956 | pTmp = mp_Mult(pTmpFull, hessenbergMat,R); |
| 957 | idDelete((ideal*)&hessenbergMat); |
| 958 | hessenbergMat = mp_Mult(pTmp, pTmpFull, R); |
| 959 | idDelete((ideal*)&pTmp); idDelete((ideal*)&pTmpFull); |
| 960 | /* as there may be inaccuracy, we erase those entries of hessenbergMat |
| 961 | which must have become zero by the last transformation */ |
| 962 | for (int r = c + 2; r <= n; r++) |
| 963 | pDelete(&MATELEM(hessenbergMat, r, c)); |
| 964 | } |
| 965 | } |
| 966 | } |
no test coverage detected