amirvahid
10-21-2009, 08:05 PM
Hi,
I am using the code zrhqr.h in order to obtain all of the roots for a 4-th order polynomial. After some MD simulations, it gives me the following error:
ERROR: Too many iterations in hqr
in file Ewald.cpp at line 201
terminate called after throwing an instance of 'int'
Aborted
Here is my code:
void PredictEvent (Int na, Int nb)//cf. Rapaport chap14
{
VecR dr, dv, rs, shift, tm, w, da;
VecI m1v, m2v;
Doub b, d, t, tInt, vv, aa, ar, av, minrr = 0;
Int dir, evCode, m1x, m1y, m1z, n;
const int pol = 4;
const int poll = pol + 1;
VecComplex rt(pol);
VecDoub Coeffs(poll);
Int indxmin = -1;
VCopy (w, mol[na].inCell);
if (mol[na].rv.x > 0.) ++ w.x;
if (mol[na].rv.y > 0.) ++ w.y;
if (mol[na].rv.z > 0.) ++ w.z;
VMul (w, w, region);
VDiv (w, w, cells);
VSAdd (rs, mol[na].r, 0.5, region);
VVSub (w, rs);
WhenCross (x);
WhenCross (y);
WhenCross (z);
if (tm.y < tm.z) dir = (tm.x < tm.y) ? 0 : 1;
else dir = (tm.x < tm.z) ? 0 : 2;
evCode = 100 + dir;
ScheduleEvent (na, MOL_LIMIT + evCode, timeNow + VComp (tm, dir));
for (m1z = cellRange[0].z; m1z <= cellRange[1].z; m1z ++) {
m1v.z = m1z;
VWrapEvC (z);
for (m1y = cellRange[0].y; m1y <= cellRange[1].y; m1y ++) {
m1v.y = m1y;
VWrapEvC (y);
for (m1x = cellRange[0].x; m1x <= cellRange[1].x; m1x ++) {
m1v.x = m1x;
VWrapEvC (x);
n = VLinear (m2v, cells) + nMol;
for (n = cellList[n]; n >= 0; n = cellList[n]) {
if (n != na && n != nb && (nb >= -1 || n < na)) {
tInt = timeNow - mol[n].time;
VSub (dr, mol[na].r, mol[n].r);
VVSAdd (dr, - tInt, mol[n].rv);
VVSAdd (dr, -0.5 * Sqr (tInt), mol[n].ra);
VVSub (dr, shift);
VSub (dv, mol[na].rv, mol[n].rv);
VVSAdd (dv, - tInt, mol[n].ra);
VSub (da, mol[na].ra, mol[n].ra);
b = VDot (dr, dv);
vv = VLenSq (dv);
aa = VLenSq (da);
ar = VDot (da, dr);
av = VDot (da, dv);
Coeffs[0] =VLenSq(dr)-1;
Coeffs[1] = 2*b;
Coeffs[2] = vv + ar;
Coeffs[3] = av;
Coeffs[4] = 0.25 * aa;
if (b < 0. && abs(Coeffs[4]) < 1e-8) {
vv = VLenSq (dv);
d = Sqr (b) - vv * (VLenSq (dr) - 1.);
if (d >= 0.) {
t = - (sqrt (d) + b) / vv;
ScheduleEvent (na, n, timeNow + t);
}
}
if (b < 0. && abs(Coeffs[4])>1e-8) {
zrhqr(Coeffs,rt);
for (int m1 = 0; m1 < pol; m1++){
if (imag(rt[m1])== 0 && real(rt[m1]) > 0) {
minrr = real (rt[m1]);
indxmin = m1;
//break;
}
if (imag(rt[m1]) == 0 && real(rt[m1]) < minrr && real(rt[m1]) > 0){
minrr = real(rt[m1]);
indxmin = m1;
}
}
if (indxmin != -1){
t = minrr;
ScheduleEvent (na, n, timeNow + t);
}
}
}
}
}
}
}
}
I marked it with "red".
Here is a general code that one can use for finding polynomial root using numerical recipes
#include "nr3.h"
#include "eigen_unsym.h"
#include "zrhqr.h"
int main()
{
const int n = 6; // degree of polynomial
Doub coeff[] = { 10.0, -18.0, 25.0, -24.0, 16.0, -6.0, 1.0};
const int N = n+1;
VecDoub a(N, coeff);
//cout << a[2] << endl;
VecComplex rt(n);
cout << endl << "Roots of polynomial x^6-6x^5+16x^4-24x^3+25x^2-18x+10" << endl;
zrhqr(a,rt);
cout << fixed << setprecision(6);
for (int i = 0; i < n; i++){
cout << setw(5) << noshowpos << i;
cout << setw(25) << showpos << rt[i] << endl;
}
return 0;
}
I'll get the following results:
Roots of polynomial x^6-6x^5+16x^4-24x^3+25x^2-18x+10
0 (+2.000000,+1.000000)
1 (+2.000000,-1.000000)
2 (+1.000000,+1.000000)
3 (+1.000000,-1.000000)
4 (+0.000000,+1.000000)
5 (+0.000000,-1.000000)
I wonder if you could help me.
I am using the code zrhqr.h in order to obtain all of the roots for a 4-th order polynomial. After some MD simulations, it gives me the following error:
ERROR: Too many iterations in hqr
in file Ewald.cpp at line 201
terminate called after throwing an instance of 'int'
Aborted
Here is my code:
void PredictEvent (Int na, Int nb)//cf. Rapaport chap14
{
VecR dr, dv, rs, shift, tm, w, da;
VecI m1v, m2v;
Doub b, d, t, tInt, vv, aa, ar, av, minrr = 0;
Int dir, evCode, m1x, m1y, m1z, n;
const int pol = 4;
const int poll = pol + 1;
VecComplex rt(pol);
VecDoub Coeffs(poll);
Int indxmin = -1;
VCopy (w, mol[na].inCell);
if (mol[na].rv.x > 0.) ++ w.x;
if (mol[na].rv.y > 0.) ++ w.y;
if (mol[na].rv.z > 0.) ++ w.z;
VMul (w, w, region);
VDiv (w, w, cells);
VSAdd (rs, mol[na].r, 0.5, region);
VVSub (w, rs);
WhenCross (x);
WhenCross (y);
WhenCross (z);
if (tm.y < tm.z) dir = (tm.x < tm.y) ? 0 : 1;
else dir = (tm.x < tm.z) ? 0 : 2;
evCode = 100 + dir;
ScheduleEvent (na, MOL_LIMIT + evCode, timeNow + VComp (tm, dir));
for (m1z = cellRange[0].z; m1z <= cellRange[1].z; m1z ++) {
m1v.z = m1z;
VWrapEvC (z);
for (m1y = cellRange[0].y; m1y <= cellRange[1].y; m1y ++) {
m1v.y = m1y;
VWrapEvC (y);
for (m1x = cellRange[0].x; m1x <= cellRange[1].x; m1x ++) {
m1v.x = m1x;
VWrapEvC (x);
n = VLinear (m2v, cells) + nMol;
for (n = cellList[n]; n >= 0; n = cellList[n]) {
if (n != na && n != nb && (nb >= -1 || n < na)) {
tInt = timeNow - mol[n].time;
VSub (dr, mol[na].r, mol[n].r);
VVSAdd (dr, - tInt, mol[n].rv);
VVSAdd (dr, -0.5 * Sqr (tInt), mol[n].ra);
VVSub (dr, shift);
VSub (dv, mol[na].rv, mol[n].rv);
VVSAdd (dv, - tInt, mol[n].ra);
VSub (da, mol[na].ra, mol[n].ra);
b = VDot (dr, dv);
vv = VLenSq (dv);
aa = VLenSq (da);
ar = VDot (da, dr);
av = VDot (da, dv);
Coeffs[0] =VLenSq(dr)-1;
Coeffs[1] = 2*b;
Coeffs[2] = vv + ar;
Coeffs[3] = av;
Coeffs[4] = 0.25 * aa;
if (b < 0. && abs(Coeffs[4]) < 1e-8) {
vv = VLenSq (dv);
d = Sqr (b) - vv * (VLenSq (dr) - 1.);
if (d >= 0.) {
t = - (sqrt (d) + b) / vv;
ScheduleEvent (na, n, timeNow + t);
}
}
if (b < 0. && abs(Coeffs[4])>1e-8) {
zrhqr(Coeffs,rt);
for (int m1 = 0; m1 < pol; m1++){
if (imag(rt[m1])== 0 && real(rt[m1]) > 0) {
minrr = real (rt[m1]);
indxmin = m1;
//break;
}
if (imag(rt[m1]) == 0 && real(rt[m1]) < minrr && real(rt[m1]) > 0){
minrr = real(rt[m1]);
indxmin = m1;
}
}
if (indxmin != -1){
t = minrr;
ScheduleEvent (na, n, timeNow + t);
}
}
}
}
}
}
}
}
I marked it with "red".
Here is a general code that one can use for finding polynomial root using numerical recipes
#include "nr3.h"
#include "eigen_unsym.h"
#include "zrhqr.h"
int main()
{
const int n = 6; // degree of polynomial
Doub coeff[] = { 10.0, -18.0, 25.0, -24.0, 16.0, -6.0, 1.0};
const int N = n+1;
VecDoub a(N, coeff);
//cout << a[2] << endl;
VecComplex rt(n);
cout << endl << "Roots of polynomial x^6-6x^5+16x^4-24x^3+25x^2-18x+10" << endl;
zrhqr(a,rt);
cout << fixed << setprecision(6);
for (int i = 0; i < n; i++){
cout << setw(5) << noshowpos << i;
cout << setw(25) << showpos << rt[i] << endl;
}
return 0;
}
I'll get the following results:
Roots of polynomial x^6-6x^5+16x^4-24x^3+25x^2-18x+10
0 (+2.000000,+1.000000)
1 (+2.000000,-1.000000)
2 (+1.000000,+1.000000)
3 (+1.000000,-1.000000)
4 (+0.000000,+1.000000)
5 (+0.000000,-1.000000)
I wonder if you could help me.