45 namespace NekMeshUtils
59 for (
int i = 0; i < ti.num_elements(); i++)
63 ti[i] = (xi[i] - ui[i]) / J(i, 0);
67 ti[i] = (xi[i] - li[i]) / J(i, 0);
71 ti[i] = numeric_limits<double>::max();
80 DNekMat d(xi.num_elements(), 1);
81 for (
int i = 0; i < xi.num_elements(); i++)
83 if (fabs(ti[i]) < 1E-10)
89 d(i, 0) = -1.0 * J(i, 0);
95 bool hitbounded =
false;
97 for (
int i = 0; i < xci.num_elements(); i++)
99 if (xi[i] + d(i, 0) < li[i])
103 cout <<
"hit bounded" << endl;
109 xci[i] = xi[i] + d(i, 0);
112 if (xi[i] + d(i, 0) > ui[i])
116 cout <<
"hit bounded" << endl;
122 xci[i] = xi[i] + d(i, 0);
126 DNekMat Z(xci.num_elements(), xci.num_elements(), 0.0);
129 for (
int i = 0; i < xci.num_elements(); i++)
132 if (it != Fset.end())
138 DNekMat dx(xci.num_elements(), 1, 0.0);
139 for (
int i = 0; i < xci.num_elements(); i++)
141 dx(i, 0) = xci[i] - xi[i];
149 for (it = Fset.begin(); it != Fset.end(); it++)
152 if (li[i] - xci[i] > alpha * du(i, 0))
154 alpha = min(alpha, (li[i] - xci[i]) / du(i, 0));
156 else if (ui[i] - xci[i] < alpha * du(i, 0))
158 alpha = min(alpha, (ui[i] - xci[i]) / du(i, 0));
165 for (
int i = 0; i < xci.num_elements(); i++)
170 xibar[i] = xci[i] + grad(i, 0);
194 Vmath::Vsub(xci.num_elements(), &xibar[0], 1, &xi[0], 1, &dk[0], 1);
199 for (
int i = 0; i < dk.num_elements(); i++)
201 c += 1E-4 * J(i, 0) * dk[i];
202 r += J(i, 0) * dk[i];
230 for (
int i = 0; i < xi.num_elements(); i++)
232 tst[i] = xi[i] + lam * dk[i];
240 for (
int i = 0; i < dk.num_elements(); i++)
242 l += jn(i, 0) * dk[i];
245 }
while (fn > fo + c || fabs(l) > 0.9 * fabs(r));
254 DNekMat s(dk.num_elements(), 1, 0.0);
255 for (
int i = 0; i < dk.num_elements(); i++)
257 s(i, 0) = lam * dk[i];
268 for (
int i = 0; i < dk.num_elements(); i++)
270 ynorm += y(i, 0) * y(i, 0);
273 if (d3(0, 0) > 2.2E-16 * ynorm)
275 B = B + y * yT * (1.0 / d1(0, 0)) - B * s * sT * B * (1.0 / d2(0, 0));
276 H = H + (d3(0, 0) + n1(0, 0)) / d3(0, 0) / d3(0, 0) * s * sT -
277 1.0 / d3(0, 0) * (H * y * sT + s * yT * H);
void Vsub(int n, const T *x, const int incx, const T *y, const int incy, T *z, const int incz)
Subtract vector z = x-y.
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
boost::shared_ptr< OptiObj > OptiObjSharedPtr
bool BGFSUpdate(OptiObjSharedPtr opti, DNekMat &J, DNekMat &B, DNekMat &H)