36 #ifndef UTILITIES_NEKMESH_NODEOPTI_HESSIAN
37 #define UTILITIES_NEKMESH_NODEOPTI_HESSIAN
56 template <>
int NodeOpti::IsIndefinite<2>()
69 char jobvl =
'N', jobvr =
'N';
70 int worklen = 8 * nVel, info;
78 Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
79 &wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
83 if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0)
85 if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0)
98 template <>
int NodeOpti::IsIndefinite<3>()
114 char jobvl =
'N', jobvr =
'N';
115 int worklen = 8 * nVel, info;
123 Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
124 &wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
128 if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0 || eval(2, 2) < 0.0)
130 if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0 && eval(2, 2))
161 NekDouble D = (H[0][0] - H[1][1]) * (H[0][0] - H[1][1]) + 4.0 * H[1][0] * H[1][0];
165 val = (H[0][0] + H[1][1] - Dsqrt ) / 2.0;
183 NekDouble p1 = H[0][1] * H[0][1] + H[0][2] * H[0][2] + H[1][2] * H[1][2];
187 if(H[0][0] < H[1][1])
189 if(H[0][0] < H[2][2])
200 if(H[1][1] < H[2][2])
212 NekDouble q = (H[0][0] + H[1][1] + H[2][2]) / 3.0;
213 NekDouble p2 = (H[0][0] - q)*(H[0][0] - q)
214 + (H[1][1] - q)*(H[1][1] - q)
215 + (H[2][2] - q)*(H[2][2] - q)
221 B[0][0] = pinv * (H[0][0] - q);
222 B[1][1] = pinv * (H[1][1] - q);
223 B[2][2] = pinv * (H[2][2] - q);
224 B[0][1] = pinv * H[0][1];
226 B[0][2] = pinv * H[0][2];
228 B[1][2] = pinv * H[1][2];
251 val = q + 2.0 * p * cos(phi + (2.0*M_PI/3.0));
#define ASSERTL0(condition, msg)
int IsIndefinite()
Returns 1 if Hessian matrix is indefinite and 0 otherwise.
void MinEigen(NekDouble &val)
Calculates minimum eigenvalue of Hessian matrix.
NekDouble Determinant< 3 >(NekDouble jac[][3])