I'm working in standard C language.

My problem is to solve a equation of order 10. And "zrhqr" in NR gives the way to do it. In "zrhqr", it calls the subfunction "hqr".
I have the problem like "Too many iterations in hqr ...now exiting to system..." for some input parameters.
The bad coefficients of the 10-order equation is :
n[0]=-47.12339401,n[1]= 6.64312363,n[2]= 85.31824493,n[3]=-44.96127701,n[4]= 29.37084198,n[5]= 4.72641611,n[6]=-83.75795746,n[7]= 61.78692245,n[8]=-12.52025509,n[9]= -0.03707334,n[10]= 0.07153397,
with n[i] represents the i-order coefficient.

In addition, for "svdcmp" I have also the same problem for some matrixes.


The topic "Possible convergence problems in svdcmp, jacobi, tqli, hqr"(http://www.nr.com/forum/showthread.php?t=292) gives some suggestion to correct it, but I don't understand the keywords "DP" and "volatile" in "DP volatile temp... ". Maybe it is C++. But I'm working in C.
And I just tried this solution(declare a volatile float variable "temp", and do what the topic suggests), but the convergence problem still exists.

But for me, ALL THE THINGS work well in MATLAB. I don't understand what causes the difference.

If everyone has some suggestions, I would be appreaciated.



Best regards,
tangfrch@gmail.com

...
I have the problem like "Too many iterations in hqr ...now exiting to system..." for some input parameters.

I used a test program like zxrhqr.c, but with your polynomial:

#define M 10 /* degree of polynomial */
#define MP1 (M+1) /* no. of polynomial coefficients */

int main(void)
{
static float a[MP1] = {
-47.12339401, 6.64312363, 85.31824493, -44.96127701,
29.37084198, 4.72641611,-83.75795746, 61.78692245,
-12.52025509, -0.03707334, 0.07153397
};
.
.
.
printf("\nRoots of the polynomial\n");
printf("\n%15s %15s\n","Real","Complex");
zrhqr(a,M,rtr,rti);
for (i=1;i<=M;i++) {
printf("%5d %12.6f %12.6f\n",i,rtr[i],rti[i]);
}
.
.
.


With this program, the function quits with the message that you indicated. (Numerical recipes C version 2.11 and GNU gcc compiler version 4.1.2)

Indeed, I observe that, in the file hqr.c, there are the lines:


if (its == 30) nrerror("Too many iterations in hqr");
if (its == 10 || its == 20) {
.
. // Other stuff
.
}



If I just change the first line that I showed to

if (its == 100) nrerror("Too many iterations in hqr");




Then the output is:

Roots of the polynomial

Real Complex
1 -15.143620 0.000000
2 -0.778530 -0.039171
3 -0.778530 0.039171
4 -0.040623 -1.046191
5 -0.040623 1.046191
6 1.049256 -0.031325
7 1.049256 0.031325
8 2.132582 0.000000
9 2.673366 0.000000
10 10.395741 0.000000



When I used the same polynomial in octave, I got:

vv =

Columns 1 through 6:

0.071534 -0.037073 -12.520255 61.786922 -83.757957 4.726416

Columns 7 through 11:

29.370842 -44.961277 85.318245 6.643124 -47.123394

octave:6> roots(vv)
ans =

-15.14363 + 0.00000i
10.39574 + 0.00000i
2.67337 + 0.00000i
2.13258 + 0.00000i
1.04926 + 0.03133i
1.04926 - 0.03133i
-0.04062 + 1.04619i
-0.04062 - 1.04619i
-0.77853 + 0.03917i
-0.77853 - 0.03917i


Which looks to me to be spot on. (I assume that the Matlab output would be similar to octave's.)

I notice that even with the latest, greatest, NR3 release, in the Unsymmeig::hqr() function in eigen_unsym.cpp, there is exactly the same thing with maximum iteration value hard-coded in:

if (its == 30) throw("Too many iterations in hqr");
if (its == 10 || its == 20) {
.
. // other stuff
.
}


And the results are the same as with the NR2 C program with your polynomial: The program quits before getting there, but when I bump up the maximum allowable iteration count, it gets through.

Maybe the maximum iteration count should be made to be a parameter or a #defined constant (or maybe a struct member for the NR3 Unsymmeig class) that could be changed by user programs???

Regards,

Dave

Footnote: I also notice that the maximum iteration count has a hard-coded value of 30 in NR2 svdcmp.c, and may cause it to bail out prematurely for some "perfectly good" cases. Same for the SVD::decompose() function in svd.h of the NR3 distribution.

But there is a little thing different: in "svdcmp.c", the variable "its" is already coded in "for (its=1;its<=30;its++)". If I increased "30" to something like "3000". Correspondingly, in "if (its == 30) nrerror("no convergence in 30 svdcmp iterations")" , I should also change the "if (its == 30)" to "if (its == 3000)". I DON'T KNOW WHETHER THE CHANGE FOR "svdcmp.c" IS GOOD OR NOT. Maybe for some really "bad" matrixes, the number of iteration is not enough, and perhaps there is the convergence problem in the topic "Possible convergence problems in svdcmp, jacobi, tqli, hqr"(http://www.nr.com/forum/showthread.php?t=292) ". If I meet the problem, I will send you the "bad" matrix and could you test it?

And I also think I should convert all the code in "double" precision.

In addition, I have another bad case, if you could test:
the problem is to use "gaussj.c" to solve the linear equation A1*x = A2
where [A1 | A2] is a 10x20 matrix as shown in following, and A1 is the first 10x10 submatrix of the following 10x20 matrix and A2 is the right 10x10 submatrix.


(row1) 0.029585 0.064582 -0.136313 0.016689 0.281869 0.066044 0.108775 -0.229152 0.311887 0.496570 -0.150138 0.509618 -0.078060 -0.041084 0.121468 0.029311 -0.139500 0.271525 0.060817 0.459730 ;
(row2) -0.260452 -0.042350 -0.483982 -0.250033 -0.147693 -0.053897 -0.123454 -0.139404 -0.214106 0.223113 0.080756 -0.269116 -0.388333 0.047137 -0.259193 -0.221410 0.045450 0.404866 -0.176126 -0.008634 ;
(row3) -0.085590 0.051184 0.147889 0.173514 0.095090 -0.067482 -0.089320 -0.000316 -0.065823 -0.126384 0.066882 -0.391447 0.067935 -0.044878 0.378131 0.020737 -0.020286 0.053951 -0.000414 -0.169133 ;
(row4) 0.036925 0.034172 -0.080981 0.044695 -0.118007 0.294856 -0.121116 0.060961 -0.015042 0.414192 -0.409640 -0.137326 0.104641 -0.274831 -0.134284 -0.158203 0.290374 -0.091750 0.351658 0.050031 ;
(row5) -0.202182 -0.077243 -0.203941 -0.264501 0.216076 -0.106958 0.066678 0.059296 0.148005 0.128631 -0.233994 -0.541190 0.183339 -0.132256 -0.294542 0.069509 0.234084 -0.266996 -0.269032 0.028121 ;
(row6) 0.087603 -0.078205 -0.022944 -0.025538 0.153599 -0.109277 -0.168392 -0.038037 -0.139858 -0.105195 -0.092441 0.197748 0.104181 0.243110 0.028158 -0.217813 -0.053039 0.103551 -0.211494 0.009501 ;
(row7) 0.071637 -0.014438 -0.003638 -0.125248 -0.132962 -0.226811 0.115727 0.259975 0.179338 0.285443 0.113412 -0.101548 -0.063871 -0.195493 0.256049 0.107518 -0.011001 0.006913 0.012162 0.002873 ;
(row8) 0.132616 -0.108027 -0.151645 -0.135052 0.057015 0.140889 -0.139902 -0.138292 -0.023229 -0.117421 -0.124299 0.036397 -0.167339 0.180499 -0.066320 0.303815 0.009978 0.001544 0.016294 0.008975 ;
(row9) 0.007155 0.058741 -0.110605 0.126532 -0.027478 0.075405 -0.050014 -0.081092 0.086079 -0.099965 -0.013881 0.017345 0.051630 0.027404 -0.054167 -0.085966 0.002964 0.001239 -0.004421 0.000924 ;
(row10) -0.049370 -0.063147 0.120073 -0.025486 0.018326 0.104872 -0.001737 0.317693 -0.025024 -0.364270 0.029100 -0.100099 0.017466 -0.053948 0.170527 -0.027289 -0.003587 0.004069 -0.000674 -0.000734 ;

I've used gaussj(A1,10,A2,10) in float precision and C language, the result will be save in A2. In my program, it gives:

11759.421875 8463.134766 -17863.427734 1851.753174 21722.306641 6886.758301 -11748.499023 20895.472656 3591.740967 5145.439453;
74676.085938 53719.257812 -113438.507812 11752.653320 137945.890625 43740.746094 -74601.531250 132687.953125 22810.121094 32668.277344;
3008.955811 2161.075439 -4571.005371 472.677063 5559.666992 1765.385498 -3005.328857 5345.070312 920.094910 1314.889160;
-34228.566406 -24620.236328 51996.683594 -5386.615723 -63229.417969 -20051.882812 34193.941406 -60818.117188 -10455.627930 -14973.028320;
8322.989258 5989.784180 -12643.424805 1310.557129 15375.372070 4873.874023 -8315.611328 14790.283203 2541.737305 3642.416748;
1709.468018 1230.516235 -2597.799072 268.898499 3158.409180 1001.792236 -1708.074097 3038.358398 523.031128 748.485840;
-10607.082031 -7625.562012 16113.992188 -1668.553345 -19596.982422 -6216.067871 10595.964844 -18847.105469 -3240.905518 -4638.572266;
4698.194336 3381.125732 -7136.770996 739.330261 8679.241211 2750.301270 -4693.820801 8348.922852 1434.877563 2056.073975;
-5011.217285 -3608.489746 7612.143555 -789.661133 -9255.598633 -2932.342285 5006.794922 -8905.336914 -1530.046265 -2193.302734;
-5749.338379 -4135.557129 8733.211914 -905.145752 -10619.764648 -3367.686768 5743.425781 -10214.930664 -1755.746216 -2514.785156;

But in matlab, I used the command "rref( [A1 | A2] )" to do gauss elimination of THE WHOLE [A1 | A2], AND A1 will become a 10x10 identity, and A2 should be the same result as the result of the "gaussj(A1,10,A2,10)". The result in Matlab is:

Columns 1 through 6

1.168391000540243e+04 8.408786710523322e+03 -1.774872177941630e+04 1.839854398595625e+03 2.158283323918593e+04 6.842577714213660e+03
7.419670817625764e+04 5.337423289321869e+04 -1.127103101068602e+05 1.167711496975762e+04 1.370604581125594e+05 4.346027061114516e+04
2.989694378662352e+03 2.147212208448494e+03 -4.541746323059907e+03 4.696419108705258e+02 5.524090322471344e+03 1.754116055544115e+03
-3.400891311320989e+04 -2.446214258740538e+04 5.166301781983009e+04 -5.352003370809163e+03 -6.282370560200904e+04 -1.992336740935288e+04
8.269563750437208e+03 5.951332462972821e+03 -1.256226914048220e+04 1.302138575845526e+03 1.527669322652071e+04 4.842616164801992e+03
1.698497300035017e+03 1.222620199134296e+03 -2.581133926291135e+03 2.671697795114558e+02 3.138145679089312e+03 9.953734229316907e+02
-1.053896413291427e+04 -7.576535023316529e+03 1.601051671203197e+04 -1.657819589679238e+03 -1.947116489798937e+04 -6.176212953595595e+03
4.667997136171566e+03 3.359391801002050e+03 -7.090899904749986e+03 7.345719176428800e+02 8.623466212120629e+03 2.732633475866378e+03
-4.979057559398943e+03 -3.585343432930499e+03 7.563291551304102e+03 -7.845935859564956e+02 -9.196198022854753e+03 -2.913526256358901e+03
-5.712412931106455e+03 -4.108980681081712e+03 8.677119943995052e+03 -8.993271981883232e+02 -1.055156176591107e+04 -3.346082334708111e+03

Columns 7 through 10

-1.167306927780407e+04 2.076126573981268e+04 3.568691316095103e+03 5.112377545362924e+03
-7.412267427846456e+04 1.318359509407706e+05 2.266379266411605e+04 3.245838590227705e+04
-2.986088326285114e+03 5.310836557188489e+03 9.142154463065374e+02 1.306455742424449e+03
3.397452411536969e+04 -6.042772360203556e+04 -1.038857886482953e+04 -1.487685388520792e+04
-8.262244307799081e+03 1.469533072128265e+04 2.525429535585281e+03 3.619025129417119e+03
-1.697115216042915e+03 3.018860054123719e+03 5.196823126218866e+02 7.436823865129286e+02
1.052792053210346e+04 -1.872603874182132e+04 -3.220112643235071e+03 -4.608747004648517e+03
-4.663656357547317e+03 8.295253696661181e+03 1.425659950533758e+03 2.042852505964923e+03
4.974670067850078e+03 -8.848178961943418e+03 -1.520229706055314e+03 -2.179221856405220e+03
5.706540313085756e+03 -1.014930208113666e+04 -1.744474769779243e+03 -2.498617606664762e+03

It is different from the result of "gaussj"

Could you explain the difference for me? I think it is from the precision(Matlab works by default in Double). Maybe I should use "gaussj" in double precision. But I'm not sure.

Hope for your help.

...I DON'T KNOW WHETHER THE CHANGE FOR "svdcmp.c" IS GOOD OR NOT.I don't either. I would never recommend that someone should change the Numerical Recipes routines. The book and the code are there for learning purposes. I can't get into the middle of your research and tell you whether something is OK or not. I suggested something that you can investigate in trying to get your NR program to run to completion. What happens next can't depend on me or my personal opinions.

perhaps there is the convergence problem in the topic "Possible convergence problems in svdcmp,
I have no idea, and no way of testing all of those things. Maybe someone else can give some insight...


And I also think I should convert all the code in "double" precision.That might be a start.


In addition, I have another bad case...Matlab...is different from the result of "gaussj"...the difference...?

If you perform SIngular Value Decomposition for the matrix of coefficients that you posted I think that you will find that the condition number (ratio of largest to smallest singular values) is something like 3.9e5, which makes solutions of systems of equations with this coefficient matrix verrrrrry susceptible to roundoff error. (Especially when dealing with single-precision floating point numbers which are probably good for six or seven significant decimal digits at best.)

Going to double precision may help. It may or may not give acceptable results. Note that, even if Matlab uses double precision and if you change gaussj to use double precision, the results may not be bit-for-bit exactly the same, depending on several implementation-dependent (compiler and library) conditions. Harsh reality.

You might try calculating the residuals of both cases

If xhat is the (approximate) solution you obtained from the system A*x = b, then calculate the residual matrix as res = (A * xhat) - b. Of course, you hope to see "very small" values (whatever "small" means in your context).

However, roundoff error comes into play when you calculate the residuals also.

The result is that, given a large condition number, it is possible to obtain two xhat approximate solutions (from different programs) that look quite different but have the result that both residual matrices have comparably small values. So, even with small residuals you still might have an approximate solution that isn't necessarily "close" to the actual solution.

The key here is not the values of the numbers the program gives you (or whether two programs give similar results), but knowing how to interpret the results.


Regards,

Dave