Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Hessian.hxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: Hessian.hxx
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: Utility functions for Hessian matrices
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
36 #ifndef UTILITIES_NEKMESH_NODEOPTI_HESSIAN
37 #define UTILITIES_NEKMESH_NODEOPTI_HESSIAN
38 
39 
40 namespace Nektar
41 {
42 namespace Utilities
43 {
44 
45 /**
46  * @brief Returns 1 if Hessian matrix is indefinite and 0 otherwise.
47  *
48  * Specialised versions of this function exist only for 2x2 and 3x3 matrices.
49  */
50 template <int DIM> int NodeOpti::IsIndefinite()
51 {
52  ASSERTL0(false, "DIM error");
53  return 0;
54 }
55 
56 template <> int NodeOpti::IsIndefinite<2>()
57 {
58  Array<OneD, NekDouble> eigR(2);
59  Array<OneD, NekDouble> eigI(2);
60  NekMatrix<NekDouble> H(2, 2);
61  H(0, 0) = m_grad[2];
62  H(1, 0) = m_grad[3];
63  H(0, 1) = H(1, 0);
64  H(1, 1) = m_grad[4];
65 
66  // cout << H << endl << endl;
67 
68  int nVel = 2;
69  char jobvl = 'N', jobvr = 'N';
70  int worklen = 8 * nVel, info;
71  NekDouble dum;
72 
73  DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
74  Array<OneD, NekDouble> vl(nVel * nVel);
75  Array<OneD, NekDouble> work(worklen);
76  Array<OneD, NekDouble> wi(nVel);
77 
78  Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
79  &wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
80 
81  ASSERTL0(!info, "dgeev failed");
82 
83  if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0)
84  {
85  if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0)
86  {
87  return 2;
88  }
89  else
90  {
91  return 1;
92  }
93  }
94 
95  return 0;
96 }
97 
98 template <> int NodeOpti::IsIndefinite<3>()
99 {
100  Array<OneD, NekDouble> eigR(3);
101  Array<OneD, NekDouble> eigI(3);
102  NekMatrix<NekDouble> H(3, 3);
103  H(0, 0) = m_grad[3];
104  H(1, 0) = m_grad[4];
105  H(0, 1) = H(1, 0);
106  H(2, 0) = m_grad[5];
107  H(0, 2) = H(2, 0);
108  H(1, 1) = m_grad[6];
109  H(2, 1) = m_grad[7];
110  H(1, 2) = H(2, 1);
111  H(2, 2) = m_grad[8];
112 
113  int nVel = 3;
114  char jobvl = 'N', jobvr = 'N';
115  int worklen = 8 * nVel, info;
116  NekDouble dum;
117 
118  DNekMat eval(nVel, nVel, 0.0, eDIAGONAL);
119  Array<OneD, NekDouble> vl(nVel * nVel);
120  Array<OneD, NekDouble> work(worklen);
121  Array<OneD, NekDouble> wi(nVel);
122 
123  Lapack::Dgeev(jobvl, jobvr, nVel, H.GetRawPtr(), nVel, &(eval.GetPtr())[0],
124  &wi[0], &vl[0], nVel, &dum, nVel, &work[0], worklen, info);
125 
126  ASSERTL0(!info, "dgeev failed");
127 
128  if (eval(0, 0) < 0.0 || eval(1, 1) < 0.0 || eval(2, 2) < 0.0)
129  {
130  if (eval(0, 0) < 0.0 && eval(1, 1) < 0.0 && eval(2, 2))
131  {
132  return 2;
133  }
134  else
135  {
136  return 1;
137  }
138  }
139 
140  return 0;
141 }
142 
143 /**
144  * @brief Calculates minimum eigenvalue of Hessian matrix.
145  *
146  * Specialised versions of this function exist only for 2x2 and 3x3 matrices.
147  */
148 template <int DIM> void NodeOpti::MinEigen(NekDouble &val)
149 {
150  ASSERTL0(false, "DIM error");
151 }
152 
153 template <> void NodeOpti::MinEigen<2>(NekDouble &val)
154 {
155  NekDouble H[2][2];
156  H[0][0] = m_grad[2];
157  H[1][0] = m_grad[3];
158  //H[0][1] = H[1][0];
159  H[1][1] = m_grad[4];
160 
161  NekDouble D = (H[0][0] - H[1][1]) * (H[0][0] - H[1][1]) + 4.0 * H[1][0] * H[1][0];
162  NekDouble Dsqrt = sqrt(D);
163 
164  //eval[0] = (H[0][0] + H[1][1] + Dsqrt ) / 2.0;
165  val = (H[0][0] + H[1][1] - Dsqrt ) / 2.0; // the minimum Eigenvalue
166 }
167 
168 template <> void NodeOpti::MinEigen<3>(NekDouble &val)
169 {
170  NekDouble H[3][3];
171  H[0][0] = m_grad[3];
172  H[1][0] = m_grad[4];
173  H[0][1] = H[1][0];
174  H[2][0] = m_grad[5];
175  H[0][2] = H[2][0];
176  H[1][1] = m_grad[6];
177  H[2][1] = m_grad[7];
178  H[1][2] = H[2][1];
179  H[2][2] = m_grad[8];
180 
181  //double eval[3]; // the eigenvalues
182 
183  NekDouble p1 = H[0][1] * H[0][1] + H[0][2] * H[0][2] + H[1][2] * H[1][2];
184  if (p1 == 0.0) // H is diagonal
185  {
186  // find the minimum Eigenvalue
187  if(H[0][0] < H[1][1])
188  {
189  if(H[0][0] < H[2][2])
190  {
191  val = H[0][0];
192  }
193  else
194  {
195  val = H[2][2];
196  }
197  }
198  else
199  {
200  if(H[1][1] < H[2][2])
201  {
202  val = H[1][1];
203  }
204  else
205  {
206  val = H[2][2];
207  }
208  }
209  }
210  else
211  {
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)
216  + 2.0 * p1;
217  NekDouble p = sqrt(p2 / 6.0);
218 
219  NekDouble B[3][3]; // B = (1.0 / p) * (H - q * I) with I being the identity matrix
220  NekDouble pinv = 1.0 / p;
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];
225  B[1][0] = B[0][1];
226  B[0][2] = pinv * H[0][2];
227  B[2][0] = B[0][2];
228  B[1][2] = pinv * H[1][2];
229  B[2][1] = B[1][2];
230 
231  NekDouble r = Determinant<3>(B) / 2.0;
232 
233  // In exact arithmetic for h symmetric matrix -1 <= r <= 1
234  // but computation error can leave it slightly outside this range.
235  NekDouble phi;
236  if (r <= -1)
237  {
238  phi = M_PI / 3.0;
239  }
240  else if (r >= 1)
241  {
242  phi = 0.0;
243  }
244  else
245  {
246  phi = acos(r) / 3.0;
247  }
248 
249  // the eigenvalues satisfy eval[2] <= eval[1] <= eval[0]
250  //eval[0] = q + 2.0 * p * cos(phi);
251  val = q + 2.0 * p * cos(phi + (2.0*M_PI/3.0));
252  //eval[1] = 3.0 * q - eval[0] - eval[2]; // since trace(H) = eval[0] + eval[1] + eval[2]
253  }
254 }
255 
256 }
257 }
258 
259 #endif
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
int IsIndefinite()
Returns 1 if Hessian matrix is indefinite and 0 otherwise.
Definition: Hessian.hxx:50
void MinEigen(NekDouble &val)
Calculates minimum eigenvalue of Hessian matrix.
Definition: Hessian.hxx:148
double NekDouble
NekDouble Determinant< 3 >(NekDouble jac[][3])
Definition: Evaluator.hxx:65