qp() in Octave 3.0.0 returns result egregiously violating input constraints
Joshua Redstone
redstone at gmail.com
Sun Apr 6 16:34:53 CDT 2008
HI Ben,
So I added bunch of printf's around that 'if' block as follows:
RowVector tmp_row = Ain.row (i);
if (iter > 180 && iter < 200) {
fprintf(stderr, "[%d %d] tmp_row,p: ", iter, i);
for (int ii = 0 ; ii < tmp_row.length() ; ++ii) {
fprintf(stderr, "%d:%g:%g ", ii, tmp_row(ii),
p(ii));
}
fprintf(stderr, "\n");
}
double tmp = tmp_row * p;
double res = tmp_row * x;
if (tmp < 0.0)
{
fprintf(stderr, "Line %d: bin(i) = %g, res=%g,
tmp=%g i=%d\n", __LINE__, bin(i), res, tmp, i);
double alpha_tmp = (bin(i) - res) / tmp;
if (alpha_tmp < alpha)
{
alpha = alpha_tmp;
fprintf(stderr, "Line %d: alpha=%g i=%d\n",
__LINE__, alpha, i);
is_block = i;
}
}
and got the output (I added '*****' to mark the value of interest):
[188 61] tmp_row,p: 0:0:-0.00438204 1:0:-2.59478e-16 2:0:-0.000523908
3:0:0.00490594 4:0:-0.101927 5:0:3.05229e-16 6:0:-1.16656e-16
7:0:1.64152e-16 8:0:4.29474e-17 9:0:-3.20386e-16 10:0:-1.46512e-16
11:0:-0.105331 12:0:-5.54741e-17 13:0:2.56617e-16 14:0:-0 15:0:-6.36029e-16
16:0:-2.84945e-16 17:0:2.05942e-17 18:0:1.89786e-16 19:0:-1.45598e-16
20:0:-0.13203 21:0:-0.1688 22:0:-0 23:0:4.4292e-17 24:0:0.0719175
25:0:-0.00535049 26:0:-3.05951e-17 27:0:0.917857 28:0:-0 29:0:-8.73269e-18
30:0:1.86084e-17 31:0:1.0408e-16 32:0:-0 33:0:-9.67336e-17 34:0:-0
35:0:6.90424e-17 36:0:1.80412e-17 37:0:4.77894e-17 38:0:-0 39:0:-1.621e-16
40:0:-0 41:0:-1.43261e-16 42:0:-1.42703e-16 43:0:-3.51633e-17 44:0:-0
45:0:-5.01889e-17 46:0:-9.70589e-17 47:0:-0 48:0:-2.81307e-17
49:0:6.05628e-17 50:0:5.36531e-17 51:0:-2.07794e-16 52:0:-5.66085e-17
53:0:-3.47238e-17 54:0:-0 55:0:-5.36905e-17 56:0:-0 57:0:-0
58:0:-7.28251e-18 59:0:4.16334e-17 60:0:-0 *******61:1:-6.93889e-17*********
62:0:-4.4353e-17 63:0:-4.66874e-17 64:0:-0 65:0:-1.91163e-18 66:0:-0
67:0:-2.01221e-17 68:0:-0 69:0:-0 70:0:-0 71:0:1.72961e-17 72:0:-0
73:0:-1.99217e-16 74:0:-5.36531e-17 75:0:6.27992e-17 76:0:-1.43757e-16
77:0:1.2659e-17 78:0:-3.03284e-17 79:0:-0.0243573 80:0:-8.51108e-17 81:0:-0
82:0:-0 83:0:-2.60964e-17 84:0:-7.46959e-17 85:0:-0 86:0:-0 87:0:5.29591e-17
88:0:3.27027e-17 89:0:-1.5377e-17 90:0:-7.93097e-17 91:0:-0
92:0:-2.96691e-17 93:0:-3.48337e-17 94:0:9.33816e-17 95:0:-6.52868e-17
96:0:-0 97:0:-7.77222e-17 98:0:-0.164806 99:0:-8.66702e-17 100:0:-0
101:0:-1.62188e-16 102:0:-0.167668 103:0:3.71916e-17 104:0:-0.143863
105:0:0.0243573 106:0:-3.59903e-17 107:0:-9.30389e-17 108:0:4.39729e-16
109:0:1.38307e-16 110:0:7.62439e-17 111:0:2.42705e-16 112:0:-0.0621635
113:0:-0 114:0:0.0621635 115:0:-6.34118e-16 116:0:-5.53839e-16
117:0:2.15295e-17 118:0:-0 119:0:-2.80813e-17 120:0:-3.26276e-16
121:0:1.19071e-18 122:0:-1.23195e-16 123:0:-3.7967e-16 124:0:2.31629e-18
125:0:-7.36494e-16 126:0:6.07431e-16 127:0:-3.16624e-16 128:0:-3.70403e-17
129:0:1.11591e-16 130:0:4.45572e-17 131:0:-1.12e-16 132:0:1.58505e-15
133:0:1.0766e-15 134:0:5.22372e-17 135:0:2.28611e-16 136:0:4.59645e-17
137:0:1.75177e-16 138:0:4.84766e-17 139:0:-3.57924e-17 140:0:-2.81022e-17
141:0:-9.8548e-17 142:0:-0 143:0:1.10398e-16 144:0:-7.18162e-17
145:0:2.42834e-17 146:0:-0 147:0:-9.46025e-17 148:0:-1.39964e-16
149:0:5.97696e-17 150:0:-2.59568e-16 151:0:4.55345e-17 152:0:-2.63416e-17
153:0:1.631e-17 154:0:1.42538e-17 155:0:2.90691e-17 156:0:9.70695e-17
157:0:-2.8143e-17 158:0:-0
Line 409: bin(i) = 0, res=-0.731546, tmp=-6.93889e-17 i=61
Line 415: alpha=-1.05427e+16 i=61
so the dot product that tmp is computing is producing a -1e-17, which is
producing the huge alpha, which steps the solution outside of feasibility.
Could someone who understands what the code is computing (which I don't)
perhaps comment on if the small value of 'tmp' is to be expected or where to
look further?
Josh
On Sun, Apr 6, 2008 at 1:32 PM, Ben Abbott <bpabbott at mac.com> wrote:
>
> On Apr 6, 2008, at 3:50 PM, Joshua Redstone wrote:
>
> > I finally got Octave to compile on my mac!
> > So I instrumented __qp__.cc and have narrowed it down to the following
> > so far:
> > There a blurb of code around line 400:
> > double alpha_tmp = (bin(i) - res) / tmp;
> > if (alpha_tmp < alpha)
> > {
> > alpha = alpha_tmp;
> > is_block = i;
> > }
> >
> > I put a printf around the computation of alpha_tmp, and around algorithm
> > iteration 168 and at loop iteration i=61:
> > bin(i) = 0.000000, res=-0.731546, tmp=-0.000000
> > and alpha = alpha_tmp sets alpha to:
> > alpha=-10542688472441812.000000
> > This is used in the update of x, which results in an "infeasible" value
> > for x. Things go downhill from there.
> > so it looks like 'tmp' is getting set to a value that's too small.
> >
> > Next step is to put more printfs and see why tmp is set weird. I gotta
> > sign off till later today/tmrw, but perhaps this helps?
> > Josh
> >
>
> hmmm ... from the release-3-0-x archive, I have __qp__.cc:414-423 as
> being
>
> if (tmp < 0.0)
> {
> double alpha_tmp = (bin(i) - res) / tmp;
>
> if (alpha_tmp < alpha)
> {
> alpha = alpha_tmp;
> is_block = i;
> }
> }
>
> Perhaps a better choice would be "tmp < -eps"?
>
> Ben
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://www.cae.wisc.edu/pipermail/bug-octave/attachments/20080406/756edd7b/attachment-0001.html
More information about the Bug-octave
mailing list