* Computes information related to one householder transformation step for * constructing the Hessenberg form of a given non-derogative matrix. * * The method assumes that all matrix entries are numbers coming from some * subfield of the reals. And that v has a non-zero first entry v_1 and a * second non-zero entry somewhere else. * Given such a vector v, it computes a number r (which will be
| 835 | * @return the number r |
| 836 | **/ |
| 837 | number hessenbergStep( |
| 838 | const matrix vVec, /**< [in] the input vector v */ |
| 839 | matrix &uVec, /**< [out] the output vector u */ |
| 840 | matrix &pMat, /**< [out] the output matrix P */ |
| 841 | const number tolerance /**< [in] accuracy for square roots */ |
| 842 | ) |
| 843 | { |
| 844 | int rr = MATROWS(vVec); |
| 845 | number vNormSquared = euclideanNormSquared(vVec); |
| 846 | number vNorm; realSqrt(vNormSquared, tolerance, vNorm); |
| 847 | /* v1 is guaranteed to be non-zero: */ |
| 848 | number v1 = pGetCoeff(MATELEM(vVec, 1, 1)); |
| 849 | bool v1Sign = true; if (nGreaterZero(v1)) v1Sign = false; |
| 850 | |
| 851 | number v1Abs = nCopy(v1); if (v1Sign) v1Abs = nInpNeg(v1Abs); |
| 852 | number t1 = nDiv(v1Abs, vNorm); |
| 853 | number one = nInit(1); |
| 854 | number t2 = nAdd(t1, one); nDelete(&t1); |
| 855 | number denominator; realSqrt(t2, tolerance, denominator); nDelete(&t2); |
| 856 | uVec = mpNew(rr, 1); |
| 857 | t1 = nDiv(v1Abs, vNorm); |
| 858 | t2 = nAdd(t1, one); nDelete(&t1); |
| 859 | t1 = nDiv(t2, denominator); nDelete(&t2); |
| 860 | MATELEM(uVec, 1, 1) = pOne(); |
| 861 | pSetCoeff(MATELEM(uVec, 1, 1), t1); /* we know that t1 != 0 */ |
| 862 | for (int r = 2; r <= rr; r++) |
| 863 | { |
| 864 | if (MATELEM(vVec, r, 1) != NULL) |
| 865 | t1 = nCopy(pGetCoeff(MATELEM(vVec, r, 1))); |
| 866 | else t1 = nInit(0); |
| 867 | if (v1Sign) t1 = nInpNeg(t1); |
| 868 | t2 = nDiv(t1, vNorm); nDelete(&t1); |
| 869 | t1 = nDiv(t2, denominator); nDelete(&t2); |
| 870 | if (!nIsZero(t1)) |
| 871 | { |
| 872 | MATELEM(uVec, r, 1) = pOne(); |
| 873 | pSetCoeff(MATELEM(uVec, r, 1), t1); |
| 874 | } |
| 875 | else nDelete(&t1); |
| 876 | } |
| 877 | nDelete(&denominator); |
| 878 | |
| 879 | /* finished building vector u; now turn to pMat */ |
| 880 | pMat = mpNew(rr, rr); |
| 881 | /* we set P := E - u * u^T, as desired */ |
| 882 | for (int r = 1; r <= rr; r++) |
| 883 | for (int c = 1; c <= rr; c++) |
| 884 | { |
| 885 | if ((MATELEM(uVec, r, 1) != NULL) && (MATELEM(uVec, c, 1) != NULL)) |
| 886 | t1 = nMult(pGetCoeff(MATELEM(uVec, r, 1)), |
| 887 | pGetCoeff(MATELEM(uVec, c, 1))); |
| 888 | else t1 = nInit(0); |
| 889 | if (r == c) { t2 = nSub(one, t1); nDelete(&t1); } |
| 890 | else t2 = nInpNeg(t1); |
| 891 | if (!nIsZero(t2)) |
| 892 | { |
| 893 | MATELEM(pMat, r, c) = pOne(); |
| 894 | pSetCoeff(MATELEM(pMat, r, c), t2); |
no test coverage detected