Rosetta
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
dixon.cc
Go to the documentation of this file.
1 // -*- mode:c++;tab-width:2;indent-tabs-mode:t;show-trailing-whitespace:t;rm-trailing-spaces:t -*-
2 // vi: set ts=2 noet:
3 //
4 // (c) Copyright Rosetta Commons Member Institutions.
5 // (c) This file is part of the Rosetta software suite and is made available under license.
6 // (c) The Rosetta software is developed by the contributing members of the Rosetta Commons.
7 // (c) For more information, see http://www.rosettacommons.org. Questions about this can be
8 // (c) addressed to University of Washington UW TechTransfer, email: license@u.washington.edu.
9 
10 /// @file dixon.cc
11 /// @brief computes the dixon resultant
12 /// @author Evangelos A. Coutsias
13 /// @author Daniel J. Mandell
14 
15 // C++ headers
16 // Utility headers
17 //#include <basic/Tracer.hh> // tracer output
18 
19 // Unit headers
20 #include <numeric/kinematic_closure/dixon.hh>
21 
22 // Numeric headers
23 #include <numeric/types.hh>
24 #include <numeric/kinematic_closure/sturm.hh>
25 #include <numeric/kinematic_closure/kinematic_closure_helpers.hh>
26 #include <numeric/linear_algebra/GeneralizedEigenSolver.hh>
27 
28 // External headers
29 #include <Eigen/Dense>
30 
31 // C++ headers
32 #include <cmath>
33 #include <iostream>
34 
35 // Constants
36 #define PP4x2_VECSIZE 7
37 #define DIXON_SIZE 8
38 #define DIXON_RESULTANT_SIZE 3
39 
40 using numeric::Real;
41 using utility::vector1;
42 
43 namespace numeric {
44 namespace kinematic_closure {
45 
46 using namespace std;
47 Eigen::IOFormat scipy(8, 0, ", ", ",\n", "[", "]", "[", "]");
48 string skipl = "\n\n";
49 
50 typedef Eigen::Matrix<Real, 8, 8> Matrix8;
51 typedef Eigen::Matrix<Real, 16, 16> Matrix16;
52 typedef Eigen::Matrix<Real, Eigen::Dynamic, 1, 0, 16, 1> Vector16;
53 typedef linear_algebra::GeneralizedEigenSolver<Matrix16> SolverType;
54 
55 // If I can rant for a little bit, using STL vectors as a matrix type was a
56 // fucking horrible idea. For one thing, the API presented by the std::vector
57 // class is totally inappropriate for doing linear algebra. More importantly,
58 // you run the risk of getting really terrible caching performance because
59 // these matrices aren't necessarily stored contiguously in memory.
60 
61 typedef vector1<Real> PseudoVector;
62 typedef vector1<PseudoVector> PseudoMatrix;
63 
64 /// dixon functions ///
65 
66 // point_value2 {{{1
67 /* C becomes the point value 2 between A and t */
68 void point_value2(const utility::vector1<Real>& A, const utility::vector1<Real>& t, utility::vector1<Real>& C) {
69  C.resize(t.size());
70  for ( unsigned i=1; i<=t.size(); i++ ) {
71  C[i] = A[1]+t[i] * (A[2]+t[i]*A[3]);
72  }
73 }
74 
75 // point_value4 {{{1
76 /* C becomes the point value 4 between A and t */
77 void point_value4(const utility::vector1<Real>& A, const utility::vector1<Real>& t, utility::vector1<Real>& C) {
78  C.resize(t.size());
79  for ( unsigned i=1; i<=t.size(); i++ ) {
80  C[i] = A[1]+t[i] * (A[2] + t[i] * (A[3] + t[i] * (A[4] + t[i]*A[5])));
81  }
82 }
83 
84 // point_value6 {{{1
85 /* C becomes the point value 6 between B and t */
86 void point_value6(const utility::vector1<Real>& B, const utility::vector1<Real>& t, utility::vector1<Real>& C) {
87  C.resize(t.size());
88  for ( unsigned i=1; i<=t.size(); i++ ) {
89  C[i] = B[1]+t[i]*(B[2]+t[i]*(B[3]+t[i]*(B[4]+t[i]*(B[5]+t[i]*(B[6]+t[i]*B[7])))));
90  }
91 }
92 
93 // point_value8 {{{1
94 /* C becomes the point value 8 between A and t */
95 void point_value8(const utility::vector1<Real>& A, const utility::vector1<Real>& t, utility::vector1<Real>& C) {
96  C.resize(t.size());
97  for ( unsigned i=1; i<=t.size(); i++ ) {
98  C[i] = A[1]+t[i]*(A[2]+t[i]*(A[3]+t[i]*(A[4]+t[i]*(A[5]+t[i]*(A[6]+t[i]*(A[7]+t[i]*(A[8]+t[i]*A[9])))))));
99  }
100 }
101 
102 // point_value16 {{{1
103 /* C becomes the point value 16 between A and t */
104 void point_value16(const utility::vector1<Real>& A, const utility::vector1<Real>& t, utility::vector1<Real>& C) {
105  C.resize(t.size());
106  for ( unsigned i=1; i<=t.size(); i++ ) {
107  C[i] = A[1]+t[i]*(A[2]+t[i]*(A[3]+t[i]*(A[4]+t[i]*(A[5]+t[i]*(A[6]+t[i]*(A[7]+t[i]*(A[8]+t[i]*(A[9]+t[i]*(A[10]+t[i]*(A[11]+t[i]*(A[12]+t[i]*(A[13]+t[i]*(A[14]+t[i]*(A[15]+t[i]*(A[16]+t[i]*A[17])))))))))))))));
108  }
109 }
110 // }}}1
111 
112 // polyProduct2x2 {{{1
113 /* C becomes the polyProduct2x2 of A and B */
114 void polyProduct2x2(const utility::vector1<Real>& A, const utility::vector1<Real>& B, utility::vector1<Real>& C) {
115  C.resize(PP4x2_VECSIZE);
116  C[5] = A[3]*B[3];
117  C[4] = A[3]*B[2]+A[2]*B[3];
118  C[3] = A[3]*B[1]+A[2]*B[2]+A[1]*B[3];
119  C[2] = A[2]*B[1]+A[1]*B[2];
120  C[1] = A[1]*B[1];
121  return;
122 }
123 
124 // polyProduct4x2 {{{1
125 /* C becomes the polyProduct4x2 of A and B */
126 void polyProduct4x2(const utility::vector1<Real>& A, const utility::vector1<Real>& B, utility::vector1<Real>& C) {
127  C.resize(7);
128  C[7] = A[5]*B[3];
129  C[6] = A[5]*B[2]+A[4]*B[3];
130  C[5] = A[5]*B[1]+A[4]*B[2]+A[3]*B[3];
131  C[4] = A[4]*B[1]+A[3]*B[2]+A[2]*B[3];
132  C[3] = A[3]*B[1]+A[2]*B[2]+A[1]*B[3];
133  C[2] = A[2]*B[1]+A[1]*B[2];
134  C[1] = A[1]*B[1];
135  return;
136 }
137 
138 // polyProduct4x4 {{{1
139 /* C becomes the polyProduct4x4 of A and B */
140 void polyProduct4x4(const utility::vector1<Real>& A, const utility::vector1<Real>& B, utility::vector1<Real>& C) {
141  C.resize(9);
142  C[9] = A[5]*B[5];
143  C[8] = A[5]*B[4]+A[4]*B[5];
144  C[7] = A[5]*B[3]+A[4]*B[4]+A[3]*B[5];
145  C[6] = A[5]*B[2]+A[4]*B[3]+A[3]*B[4]+A[2]*B[5];
146  C[5] = A[5]*B[1]+A[4]*B[2]+A[3]*B[3]+A[2]*B[4]+A[1]*B[5];
147  C[4] = A[4]*B[1]+A[3]*B[2]+A[2]*B[3]+A[1]*B[4];
148  C[3] = A[3]*B[1]+A[2]*B[2]+A[1]*B[3];
149  C[2] = A[2]*B[1]+A[1]*B[2];
150  C[1] = A[1]*B[1];
151  return;
152 }
153 
154 // polyProduct4sq {{{1
155 /* C becomes the polyProduct4sq of A*A */
156 void polyProduct4sq(const utility::vector1<Real>& A, utility::vector1<Real>& C) {
157  C.resize(9);
158  C[9] = A[5]*A[5];
159  C[8] = 2* A[5]*A[4];
160  C[7] = A[4]*A[4] + 2* A[5]*A[3];
161  C[6] = 2*(A[5]*A[2]+A[4]*A[3]);
162  C[5] = A[3]*A[3] + 2*(A[5]*A[1]+A[4]*A[2]);
163  C[4] = 2*(A[4]*A[1]+A[3]*A[2]);
164  C[3] = A[2]*A[2] + 2* A[3]*A[1];
165  C[2] = 2* A[2]*A[1];
166  C[1] = A[1]*A[1];
167  return;
168 }
169 
170 // polyProduct6x6 {{{1
171 /* C becomes the polyProduct6x6 of A and B */
172 void polyProduct6x6(const utility::vector1<Real>& A, const utility::vector1<Real>& B, utility::vector1<Real>& C) {
173  C.resize(13);
174  C[13] = A[7]*B[7];
175  C[12] = A[7]*B[6]+A[6]*B[7];
176  C[11] = A[7]*B[5]+A[6]*B[6]+A[5]*B[7];
177  C[10] = A[7]*B[4]+A[6]*B[5]+A[5]*B[6]+A[4]*B[7];
178  C[9] = A[7]*B[3]+A[6]*B[4]+A[5]*B[5]+A[4]*B[6]+A[3]*B[7];
179  C[8] = A[7]*B[2]+A[6]*B[3]+A[5]*B[4]+A[4]*B[5]+A[3]*B[6]+A[2]*B[7];
180  C[7] = A[7]*B[1]+A[6]*B[2]+A[5]*B[3]+A[4]*B[4]+A[3]*B[5]+A[2]*B[6]+A[1]*B[7];
181  C[6] = A[6]*B[1]+A[5]*B[2]+A[4]*B[3]+A[3]*B[4]+A[2]*B[5]+A[1]*B[6];
182  C[5] = A[5]*B[1]+A[4]*B[2]+A[3]*B[3]+A[2]*B[4]+A[1]*B[5];
183  C[4] = A[4]*B[1]+A[3]*B[2]+A[2]*B[3]+A[1]*B[4];
184  C[3] = A[3]*B[1]+A[2]*B[2]+A[1]*B[3];
185  C[2] = A[2]*B[1]+A[1]*B[2];
186  C[1] = A[1]*B[1];
187  return;
188 }
189 
190 // polyProduct12x4 {{{1
191 /* C becomes the polyProduct12x4 of A and B */
192 void polyProduct12x4(const utility::vector1<Real>& A, const utility::vector1<Real>& B, utility::vector1<Real>& C) {
193  C.resize(17);
194  C[17] = A[13]*B[5];
195  C[16] = A[13]*B[4]+A[12]*B[5];
196  C[15] = A[13]*B[3]+A[12]*B[4]+A[11]*B[5];
197  C[14] = A[13]*B[2]+A[12]*B[3]+A[11]*B[4]+A[10]*B[5];
198  C[13] = A[13]*B[1]+A[12]*B[2]+A[11]*B[3]+A[10]*B[4]+A[9]*B[5];
199  C[12] = A[12]*B[1]+A[11]*B[2]+A[10]*B[3]+A[ 9]*B[4]+A[8]*B[5];
200  C[11] = A[11]*B[1]+A[10]*B[2]+A[ 9]*B[3]+A[ 8]*B[4]+A[7]*B[5];
201  C[10] = A[10]*B[1]+A[ 9]*B[2]+A[ 8]*B[3]+A[ 7]*B[4]+A[6]*B[5];
202  C[ 9] = A[ 9]*B[1]+A[ 8]*B[2]+A[ 7]*B[3]+A[ 6]*B[4]+A[5]*B[5];
203  C[ 8] = A[ 8]*B[1]+A[ 7]*B[2]+A[ 6]*B[3]+A[ 5]*B[4]+A[4]*B[5];
204  C[ 7] = A[ 7]*B[1]+A[ 6]*B[2]+A[ 5]*B[3]+A[ 4]*B[4]+A[3]*B[5];
205  C[ 6] = A[ 6]*B[1]+A[ 5]*B[2]+A[ 4]*B[3]+A[ 3]*B[4]+A[2]*B[5];
206  C[ 5] = A[ 5]*B[1]+A[ 4]*B[2]+A[ 3]*B[3]+A[ 2]*B[4]+A[1]*B[5];
207  C[ 4] = A[ 4]*B[1]+A[ 3]*B[2]+A[ 2]*B[3]+A[ 1]*B[4];
208  C[ 3] = A[ 3]*B[1]+A[ 2]*B[2]+A[ 1]*B[3];
209  C[ 2] = A[ 2]*B[1]+A[ 1]*B[2];
210  C[ 1] = A[ 1]*B[1];
211  return;
212 }
213 
214 // polyProduct8x8 {{{1
215 /* C becomes the polyProduct8x8 of A and B */
216 void polyProduct8x8(const utility::vector1<Real>& A, const utility::vector1<Real>& B, utility::vector1<Real>& C) {
217  C.resize(17);
218  C[17] = A[9]*B[9];
219  C[16] = A[9]*B[8]+A[8]*B[9];
220  C[15] = A[9]*B[7]+A[8]*B[8]+A[7]*B[9];
221  C[14] = A[9]*B[6]+A[8]*B[7]+A[7]*B[8]+A[6]*B[9];
222  C[13] = A[9]*B[5]+A[8]*B[6]+A[7]*B[7]+A[6]*B[8]+A[5]*B[9];
223  C[12] = A[9]*B[4]+A[8]*B[5]+A[7]*B[6]+A[6]*B[7]+A[5]*B[8]+A[4]*B[9];
224  C[11] = A[9]*B[3]+A[8]*B[4]+A[7]*B[5]+A[6]*B[6]+A[5]*B[7]+A[4]*B[8]+A[3]*B[9];
225  C[10] = A[9]*B[2]+A[8]*B[3]+A[7]*B[4]+A[6]*B[5]+A[5]*B[6]+A[4]*B[7]+A[3]*B[8]+A[2]*B[9];
226  C[ 9] = A[9]*B[1]+A[8]*B[2]+A[7]*B[3]+A[6]*B[4]+A[5]*B[5]+A[4]*B[6]+A[3]*B[7]+A[2]*B[8]+A[1]*B[9];
227  C[ 8] = A[8]*B[1]+A[7]*B[2]+A[6]*B[3]+A[5]*B[4]+A[4]*B[5]+A[3]*B[6]+A[2]*B[7]+A[1]*B[8];
228  C[ 7] = A[7]*B[1]+A[6]*B[2]+A[5]*B[3]+A[4]*B[4]+A[3]*B[5]+A[2]*B[6]+A[1]*B[7];
229  C[ 6] = A[6]*B[1]+A[5]*B[2]+A[4]*B[3]+A[3]*B[4]+A[2]*B[5]+A[1]*B[6];
230  C[ 5] = A[5]*B[1]+A[4]*B[2]+A[3]*B[3]+A[2]*B[4]+A[1]*B[5];
231  C[ 4] = A[4]*B[1]+A[3]*B[2]+A[2]*B[3]+A[ 1]*B[4];
232  C[ 3] = A[3]*B[1]+A[2]*B[2]+A[1]*B[3];
233  C[ 2] = A[2]*B[1]+A[1]*B[2];
234  C[ 1] = A[1]*B[1];
235  return;
236 }
237 
238 // polyProduct8sq {{{1
239 /* C becomes the polyProduct8sq of A*A */
240 void polyProduct8sq(const utility::vector1<Real>& A, utility::vector1<Real>& C) {
241  C.resize(17);
242  C[17] = A[9]*A[9];
243  C[16] = 2* A[9]*A[8];
244  C[15] = A[8]*A[8] + 2* A[9]*A[7];
245  C[14] = 2*(A[9]*A[6]+A[8]*A[7]);
246  C[13] = A[7]*A[7] + 2*(A[9]*A[5]+A[8]*A[6]);
247  C[12] = 2*(A[9]*A[4]+A[8]*A[5]+A[7]*A[6]);
248  C[11] = A[6]*A[6] + 2*(A[9]*A[3]+A[8]*A[4]+A[7]*A[5]);
249  C[10] = 2*(A[9]*A[2]+A[8]*A[3]+A[7]*A[4]+A[6]*A[5]);
250  C[ 9] = A[5]*A[5] + 2*(A[9]*A[1]+A[8]*A[2]+A[7]*A[3]+A[6]*A[4]);
251  C[ 8] = 2*(A[8]*A[1]+A[7]*A[2]+A[6]*A[3]+A[5]*A[4]);
252  C[ 7] = A[4]*A[4] + 2*(A[7]*A[1]+A[6]*A[2]+A[5]*A[3]);
253  C[ 6] = 2*(A[6]*A[1]+A[5]*A[2]+A[4]*A[3]);
254  C[ 5] = A[3]*A[3] + 2*(A[5]*A[1]+A[4]*A[2]);
255  C[ 4] = 2*(A[4]*A[1]+A[3]*A[2]);
256  C[ 3] = A[2]*A[2] + 2* A[3]*A[1];
257  C[ 2] = 2* A[2]*A[1];
258  C[ 1] = A[1]*A[1];
259  return;
260 }
261 // }}}1
262 
263 // vectorDiff {{{1
264 /* C becomes the difference of A - B */
265 /* Assumes A and B are the same size */
266 void vectorDiff(const utility::vector1<Real>& A, const utility::vector1<Real>& B, utility::vector1<Real>& C) {
267  C.resize(A.size());
268  for ( unsigned i=1; i<=A.size(); i++ ) {
269  C[i]=A[i]-B[i];
270  }
271  return;
272 }
273 
274 
275 // dixonResultant {{{1
276 /* R becomes the Dixon resultant of the determinant of the matrix polynomial
277 D(t3) := R{1} + R{2} * t3 + R{3} * t3^2
278 */
279 void dixonResultant(const utility::vector1<utility::vector1<Real> >& A,
280  const utility::vector1<utility::vector1<Real> >& B,
281  const utility::vector1<utility::vector1<Real> >& C,
282  const utility::vector1<utility::vector1<Real> >& D,
283  utility::vector1<utility::vector1<utility::vector1<Real> > >& R) {
284  R.resize(DIXON_RESULTANT_SIZE);
285  for ( int i=1; i<=3; i++ ) {
286  R[i].resize(DIXON_SIZE);
287  //first row
288  R[i][1].resize(DIXON_SIZE);
289  R[i][1][1]=0.0;
290  R[i][1][2]=A[1][i];
291  R[i][1][3]=A[2][i];
292  R[i][1][4]=A[3][i];
293  R[i][1][5]=0.0;
294  R[i][1][6]=B[1][i];
295  R[i][1][7]=B[2][i];
296  R[i][1][8]=B[3][i];
297  // second row
298  R[i][2].resize(DIXON_SIZE);
299  R[i][2][1]=A[1][i];
300  R[i][2][2]=A[2][i];
301  R[i][2][3]=A[3][i];
302  R[i][2][4]=0.0;
303  R[i][2][5]=B[1][i];
304  R[i][2][6]=B[2][i];
305  R[i][2][7]=B[3][i];
306  R[i][2][8]=0.0;
307  // third row
308  R[i][3].resize(DIXON_SIZE);
309  R[i][3][1]=0.0;
310  R[i][3][2]=B[1][i];
311  R[i][3][3]=B[2][i];
312  R[i][3][4]=B[3][i];
313  R[i][3][5]=0.0;
314  R[i][3][6]=C[1][i];
315  R[i][3][7]=C[2][i];
316  R[i][3][8]=C[3][i];
317  // forth row
318  R[i][4].resize(DIXON_SIZE);
319  R[i][4][1]=B[1][i];
320  R[i][4][2]=B[2][i];
321  R[i][4][3]=B[3][i];
322  R[i][4][4]=0.0;
323  R[i][4][5]=C[1][i];
324  R[i][4][6]=C[2][i];
325  R[i][4][7]=C[3][i];
326  R[i][4][8]=0.0;
327  // fifth row
328  R[i][5].resize(DIXON_SIZE);
329  R[i][5][1]=0.0;
330  R[i][5][2]=0.0;
331  R[i][5][3]=0.0;
332  R[i][5][4]=0.0;
333  R[i][5][5]=0.0;
334  R[i][5][6]=D[1][i];
335  R[i][5][7]=D[2][i];
336  R[i][5][8]=D[3][i];
337  // sixth row
338  R[i][6].resize(DIXON_SIZE);
339  R[i][6][1]=0.0;
340  R[i][6][2]=0.0;
341  R[i][6][3]=0.0;
342  R[i][6][4]=0.0;
343  R[i][6][5]=D[1][i];
344  R[i][6][6]=D[2][i];
345  R[i][6][7]=D[3][i];
346  R[i][6][8]=0.0;
347  // seventh row
348  R[i][7].resize(DIXON_SIZE);
349  R[i][7][1]=0.0;
350  R[i][7][2]=D[1][i];
351  R[i][7][3]=D[2][i];
352  R[i][7][4]=D[3][i];
353  R[i][7][5]=0.0;
354  R[i][7][6]=0.0;
355  R[i][7][7]=0.0;
356  R[i][7][8]=0.0;
357  // eighth row
358  R[i][8].resize(DIXON_SIZE);
359  R[i][8][1]=D[1][i];
360  R[i][8][2]=D[2][i];
361  R[i][8][3]=D[3][i];
362  R[i][8][4]=0.0;
363  R[i][8][5]=0.0;
364  R[i][8][6]=0.0;
365  R[i][8][7]=0.0;
366  R[i][8][8]=0.0;
367  }
368 }
369 // }}}1
370 
371 // build_dixon_matrices {{{1
372 void build_dixon_matrices(
373  PseudoMatrix const & A, PseudoMatrix const & B,
374  PseudoMatrix const & C, PseudoMatrix const & D,
375  Matrix8 & R0, Matrix8 & R1, Matrix8 & R2) {
376 
377  // The variable names used here mostly correspond to those used by:
378  // Coutsias, Seok, Wester, Dill; Int. J. Quant. Chem. 106:176-189 (2006).
379 
380  Matrix8 * R[3] = {&R0, &R1, &R2};
381 
382  for ( int i = 1; i <= 3; i++ ) {
383  Matrix8 & Ri = *R[i-1];
384  Ri.fill(0);
385 
386  Ri(0, 1) = A[1][i]; Ri(0, 2) = A[2][i]; Ri(0, 3) = A[3][i];
387  Ri(1, 0) = A[1][i]; Ri(1, 1) = A[2][i]; Ri(1, 2) = A[3][i];
388 
389  Ri(0, 5) = B[1][i]; Ri(0, 6) = B[2][i]; Ri(0, 7) = B[3][i];
390  Ri(1, 4) = B[1][i]; Ri(1, 5) = B[2][i]; Ri(1, 6) = B[3][i];
391  Ri(2, 1) = B[1][i]; Ri(2, 2) = B[2][i]; Ri(2, 3) = B[3][i];
392  Ri(3, 0) = B[1][i]; Ri(3, 1) = B[2][i]; Ri(3, 2) = B[3][i];
393 
394  Ri(2, 5) = C[1][i]; Ri(2, 6) = C[2][i]; Ri(2, 7) = C[3][i];
395  Ri(3, 4) = C[1][i]; Ri(3, 5) = C[2][i]; Ri(3, 6) = C[3][i];
396 
397  Ri(4, 5) = D[1][i]; Ri(4, 6) = D[2][i]; Ri(4, 7) = D[3][i];
398  Ri(5, 4) = D[1][i]; Ri(5, 5) = D[2][i]; Ri(5, 6) = D[3][i];
399  Ri(6, 1) = D[1][i]; Ri(6, 2) = D[2][i]; Ri(6, 3) = D[3][i];
400  Ri(7, 0) = D[1][i]; Ri(7, 1) = D[2][i]; Ri(7, 2) = D[3][i];
401  }
402 }
403 
404 // build_sin_and_cos {{{1
405 void build_sin_and_cos(
406  PseudoMatrix const & u, PseudoMatrix & sin, PseudoMatrix & cos) {
407 
408  // The variable names used here mostly correspond to those used by:
409  // Coutsias, Seok, Wester, Dill; Int. J. Quant. Chem. 106:176-189 (2006).
410 
411  int num_solutions = u.size();
412  double u_squared, u_squared_plus_1;
413 
414  cos.resize(num_solutions);
415  sin.resize(num_solutions);
416 
417  for ( int i = 1; i <= num_solutions; i++ ) {
418  cos[i].resize(3);
419  sin[i].resize(3);
420 
421  for ( int j = 1; j <= 3; j++ ) {
422  u_squared = u[i][j] * u[i][j];
423  u_squared_plus_1 = u_squared + 1;
424 
425  cos[i][j] = (1 - u_squared) / u_squared_plus_1;
426  sin[i][j] = 2 * u[i][j] / u_squared_plus_1;
427  }
428  }
429 }
430 
431 void dixon_eig( // {{{1
432  PseudoMatrix const & A,
433  PseudoMatrix const & B,
434  PseudoMatrix const & C,
435  PseudoMatrix const & D,
436  vector1<int> const & /*order*/,
437  PseudoMatrix & cos,
438  PseudoMatrix & sin,
439  PseudoMatrix & u,
440  int & num_solutions) {
441 
442  // The variable names used here mostly correspond to those used by:
443  // Coutsias, Seok, Wester, Dill; Int. J. Quant. Chem. 106:176-189 (2006).
444  //
445  // An exception is made for the script A and B variables, which represent the
446  // 16x16 matrices used to setup to generalized eigenvalue problem, because
447  // the non-script A and B variables are already in use. Instead, script A
448  // and B will be named U and V, respectively.
449 
450  // Fill in the Dixon matrices.
451  Matrix8 R0, R1, R2;
452  build_dixon_matrices(A, B, C, D, R0, R1, R2);
453 
454  // Setup the eigenvalue problem.
455  Matrix16 U = Matrix16::Zero();
456  Matrix16 V = Matrix16::Zero();
457  Matrix8 I = Matrix8::Identity();
458 
459  U.topRightCorner(8, 8) = I;
460  U.bottomLeftCorner(8, 8) = -R0;
461  U.bottomRightCorner(8, 8) = -R1;
462 
463  V.topLeftCorner(8, 8) = I;
464  V.bottomRightCorner(8, 8) = R2;
465 
466  // Solve the eigenvalue problem.
467  SolverType solver(U, V);
468  SolverType::RealEigenvalueType eigenvalues = solver.real_eigenvalues();
469  SolverType::RealEigenvectorType eigenvectors = solver.real_eigenvectors();
470 
471  // Extract the closure variables.
472  num_solutions = solver.num_real_solutions();
473  u.resize(num_solutions);
474 
475  for ( int i = 0; i < num_solutions; i++ ) {
476  u[i+1].resize(3);
477 
478  if ( std::abs(eigenvectors(0, i)) > 1e-8 ) {
479  u[i+1][1] = eigenvectors(1, i) / eigenvectors(0, i);
480  u[i+1][2] = eigenvectors(4, i) / eigenvectors(0, i);
481  } else {
482  u[i+1][1] = eigenvectors(3, i) / eigenvectors(2, i);
483  u[i+1][2] = eigenvectors(6, i) / eigenvectors(2, i);
484  }
485 
486  u[i+1][3] = eigenvalues(i);
487  }
488 
489  build_sin_and_cos(u, sin, cos);
490 }
491 
492 void dixon_sturm( // {{{1
493  const utility::vector1<utility::vector1<Real> >& A,
494  const utility::vector1<utility::vector1<Real> >& B,
495  const utility::vector1<utility::vector1<Real> >& C,
496  const utility::vector1<utility::vector1<Real> >& D,
497  const utility::vector1<int>& order,
498  utility::vector1<utility::vector1<Real> >& cos,
499  utility::vector1<utility::vector1<Real> >& sin,
500  utility::vector1<utility::vector1<Real> >& tau, int& nsol) {
501 
502  using utility::vector1;
503  using namespace numeric::kinematic_closure;
504 
505  utility::vector1<Real> ddet (17); // Dixon determinant
506  utility::vector1<Real> A0D1, A1D0, A1D2, A2D1, A0D2, A2D0, B0D1, B0D2, B1D0, B1D2, B2D0, B2D1, C0D1, C0D2, C1D0, C1D2, C2D0, C2D1, D0D2; // 1. Products of 2 quadratics
507  utility::vector1<Real> Det11, Det12, Det13, Det21, Det22, Det23, Det31, Det32, Det33; // 2. Sums of quartics;
508  utility::vector1<Real> D1131, D1132, D1231, D1232, D1233, D1332, D1333, D2122, D2223; // 3. Products of 2 quartics
509  utility::vector1<Real> D2121, D2222, D2323; // 4. Squares of quartics
510  utility::vector1<Real> Det1256, Det3478, Det1257, Det3468, Det1356, Det2478, Det1357; // 5. Sum of octics
511  utility::vector1<Real> C0Det23, C1Det22, C2Det21, A0Det23, A1Det22, A2Det21; // 6. Products of a quadratic with a quartic
512  utility::vector1<Real> Det678, Det123; // 7. Sums of sextics
513  utility::vector1<Real> Det123_678; // 8. Product of sextics
514  utility::vector1<Real> D1235_4678; // 9. Product of a quartic with a degree-12 poly
515  utility::vector1<Real> D1256_3478, D1257_3468, D1356_2478; // 10. Products of octics
516  utility::vector1<Real> D1357_2468; // 11. square of octic
517  utility::vector1<Real> pd1, pd2, pd3; // point_value2 results
518  utility::vector1<Real> pdet21, pdet22, pdet23, pdet31, pdet32, pdet33, pd0d2; // point_value4 results
519  utility::vector1<Real> cram12, crn112, crn262; // point_value6 results
520  utility::vector1<Real> cram22, cram32, cram42, cram52, cram122, cram132, cram222, cram232, cram242, cram252; // point_value8 results
521  utility::vector1<Real> cram11, cram21, cram31, cram41, cram51, crden; // sums and differences of products of point_value results used to find denominator crden
522  utility::vector1<Real> crn111, crn121, crn131, crn122, crn132; // sums and differences of products of point_value results used to find Cramer numerator of tau[1]
523  utility::vector1<Real> crn221, crn222, crn231, crn232, crn241, crn242, crn251, crn252, crn261; // products of point_value results used to find Cramer numerator of tau[2]
524  utility::vector1<Real> cram_num1, cram_num2; // Cramer numerators for tau[1] and tau[2]
525  utility::vector1<Real> z1, z2; // tau[1] and tau[2] numerator divided element-wise by denominator crden
526  //utility::vector1<utility::vector1<Real> > tau (16); // half tangents
527  utility::vector1<utility::vector1<Real> > tsq (16); // tau.^2
528  utility::vector1<utility::vector1<Real> > tsq1 (16); // tsq+1
529  // sturm initialization constants
530  Real tol_secant=1.0e-15;
531  int max_iter_sturm=100, max_iter_secant=20;
532  int p_order=16; // order of polynomial for sturm
533  utility::vector1<Real> Q (17); // polynomial coefficients for sturm solver
534  utility::vector1<Real> roots; // roots found by sturm solver
535 
536  ////----------------------------------------------------------------------------
537  // A: Characteristic polynomial via Lagrange 4x4-expansion of Dixon determinant.
538  // Roots of characteristic polynomial will give tau[3]. In B, we'll use
539  // Cramer's rule to derive tau[1] and tau[2] from tau[3].
540  ////----------------------------------------------------------------------------
541 
542  //-----------------------------------
543  //A1. Products of 2 quadratics
544  // (19 * 13ops = 247ops]
545  //-----------------------------------
546  polyProduct2x2(A[1],D[2],A0D1); // A0*D1
547  polyProduct2x2(A[2],D[1],A1D0); // A1*D0:
548  polyProduct2x2(A[2],D[3],A1D2); // A1*D2:
549  polyProduct2x2(A[3],D[2],A2D1); // A2*D1:
550  polyProduct2x2(A[1],D[3],A0D2); // A0*D2:
551  polyProduct2x2(A[3],D[1],A2D0); // A2*D0:
552  polyProduct2x2(B[1],D[2],B0D1); // B0*D1:
553  polyProduct2x2(B[1],D[3],B0D2); // B0*D2:
554  polyProduct2x2(B[2],D[1],B1D0); // B1*D0:
555  polyProduct2x2(B[2],D[3],B1D2); // B1*D2:
556  polyProduct2x2(B[3],D[1],B2D0); // B2*D0:
557  polyProduct2x2(B[3],D[2],B2D1); // B2*D1:
558  polyProduct2x2(C[1],D[2],C0D1); // C0*D1:
559  polyProduct2x2(C[1],D[3],C0D2); // C0*D2:
560  polyProduct2x2(C[2],D[1],C1D0); // C1*D0:
561  polyProduct2x2(C[2],D[3],C1D2); // C1*D2:
562  polyProduct2x2(C[3],D[1],C2D0); // C2*D0:
563  polyProduct2x2(C[3],D[2],C2D1); // C2*D1:
564  polyProduct2x2(D[1],D[3],D0D2); // D0*D2:
565 
566  //-----------------------------------
567  //A2. Sums of quartics
568  // (9 * 5ops = 45ops)
569  //-----------------------------------
570 
571  vectorDiff(A0D1, A1D0, Det11);
572  vectorDiff(A0D2, A2D0, Det12);
573  vectorDiff(A1D2, A2D1, Det13);
574  vectorDiff(B0D1, B1D0, Det21);
575  vectorDiff(B0D2, B2D0, Det22);
576  vectorDiff(B1D2, B2D1, Det23);
577  vectorDiff(C0D1, C1D0, Det31);
578  vectorDiff(C0D2, C2D0, Det32);
579  vectorDiff(C1D2, C2D1, Det33);
580 
581  //-----------------------------------
582  //A3. Products of 2 quartics
583  // (9 * 41ops = 369ops)
584  //-----------------------------------
585 
586  polyProduct4x4(Det11,Det31,D1131); // Det11*Det31:
587  polyProduct4x4(Det11,Det32,D1132); // Det11*Det32:
588  polyProduct4x4(Det12,Det31,D1231); // Det12*Det31:
589  polyProduct4x4(Det12,Det32,D1232); // Det12*Det32:
590  polyProduct4x4(Det12,Det33,D1233); // Det12*Det33:
591  polyProduct4x4(Det13,Det32,D1332); // Det13*Det32:
592  polyProduct4x4(Det13,Det33,D1333); // Det13*Det33:
593  polyProduct4x4(Det21,Det22,D2122); // Det21*Det22:
594  polyProduct4x4(Det22,Det23,D2223); // Det22*Det23:
595 
596  //-----------------------------------
597  //A4. Squares of quartics
598  // (3 * 28ops = 84ops)
599  //-----------------------------------
600 
601  polyProduct4sq(Det21, D2121); // Det21^2:
602  polyProduct4sq(Det22, D2222); // Det22^2:
603  polyProduct4sq(Det23, D2323); // Det23^2:
604 
605  //-----------------------------------
606  //A5. sums of octics
607  // (7 * 9ops = 63ops)
608  //-----------------------------------
609  // second factor of 5th product identical with first
610  // not used (shown for clarity)
611  // Det2468 = Det1357;
612 
613  vectorDiff(D2121,D1131,Det1256);
614  vectorDiff(D2323,D1333,Det3478);
615  vectorDiff(D2122,D1132,Det1257);
616  vectorDiff(D2223,D1332,Det3468);
617  vectorDiff(D2122,D1231,Det1356);
618  vectorDiff(D2223,D1233,Det2478);
619  vectorDiff(D2222,D1232,Det1357);
620 
621  //-----------------------------------
622  //A6. products of a quadratic with a quartic
623  // (6 * 23 = 138ops)
624  //-----------------------------------
625 
626  polyProduct4x2(Det23,C[1],C0Det23); // C0*Det23:
627  polyProduct4x2(Det22,C[2],C1Det22); // C1*Det22:
628  polyProduct4x2(Det21,C[3],C2Det21); // C2*Det21:
629  polyProduct4x2(Det23,A[1],A0Det23); // A0*Det23:
630  polyProduct4x2(Det22,A[2],A1Det22); // A1*Det22:
631  polyProduct4x2(Det21,A[3],A2Det21); // A2*Det21:
632 
633  //-----------------------------------
634  //A7. sums of sextics
635  // (4 * 7ops = 28ops)
636  //-----------------------------------
637 
638  Det678.resize(PP4x2_VECSIZE); // same size as polyProduct4x2 output
639  Det123.resize(PP4x2_VECSIZE); // same size as polyProduct4x2 output
640  // compute Det678 and Det123
641  for ( int i=1; i<=PP4x2_VECSIZE; i++ ) {
642  Det678[i]=(-C0Det23[i] + C1Det22[i] - C2Det21[i]);
643  Det123[i]=( A0Det23[i] - A1Det22[i] + A2Det21[i]);
644  }
645 
646  //-----------------------------------
647  //A8. product of sextics
648  // ( 85ops)
649  //-----------------------------------
650 
651  polyProduct6x6(Det123,Det678,Det123_678); // Det123*Det678;
652 
653  //-----------------------------------
654  //A9. product of a quartic with a degree-12 poly
655  // ( 113ops)
656  //-----------------------------------
657 
658  polyProduct12x4(Det123_678,D0D2,D1235_4678); // D0D2*Det123_678:
659 
660  //-----------------------------------
661  //A10. products of octics
662  // ( 3 * 145ops = 435ops)
663  //-----------------------------------
664 
665  polyProduct8x8(Det1256,Det3478,D1256_3478); // Det1256*Det3478:
666  polyProduct8x8(Det1257,Det3468,D1257_3468); // Det1257*Det3468:
667  polyProduct8x8(Det1356,Det2478,D1356_2478); // Det1356*Det2478:
668 
669  //-----------------------------------
670  //A11. square of octic -- since Det2468 = Det1357
671  // ( 88ops)
672  //-----------------------------------
673  // 6th product is identical to the first
674  // Not used (shown for clarity)
675  // D1567_2348 = D1235_4678;
676 
677  polyProduct8sq(Det1357,D1357_2468); // Det1357*Det2468:
678 
679  //-----------------------------------
680  //A12. Sums of 16th degree polynomials
681  // ( 5 * 17ops = 85ops)
682  //-----------------------------------
683 
684  for ( unsigned i=1; i<=ddet.size(); i++ ) {
685  ddet[i] = (-2*D1235_4678[i] + D1256_3478[i] - D1257_3468[i] - D1356_2478[i] + D1357_2468[i]);
686  }
687 
688  //-----------------------------------
689  //A13. Solve for the roots
690  //-----------------------------------
691 
692  // Get the coefficients
693  for ( unsigned i=1; i<=ddet.size(); i++ ) {
694  Q[i]=ddet[i] / ddet[17];
695  }
696 
697  // Get the roots and put them into tau[3]
698  initialize_sturm(&tol_secant, &max_iter_sturm, &max_iter_secant);
699  solve_sturm(p_order, nsol, Q, roots);
700  if ( nsol==0 ) return; // no solns found
701 
702  // remove negative roots // DJM: 7-10-2007: KEEP ALL ROOTS!
703  //utility::vector1<Real>::iterator rootIterator=roots.end();
704  //while(rootIterator != roots.begin()) {
705  // rootIterator--;
706  // if (*rootIterator < 0) {
707  // roots.erase(rootIterator);
708  // }
709  //}
710 
711  nsol=roots.size();
712  tau.resize(nsol);
713  for ( int i=1; i<=nsol; i++ ) {
714  tau[i].resize(3);
715  tau[i][order[3]]=roots[i]; // restoring column-major indexing!
716  }
717 
718  ////---------------------------------------------------------------------
719  // B: Determination of tau[1] and tau[2] from tau[3] using Cramers' rule
720  ////---------------------------------------------------------------------
721 
722  //-----------------------------------
723  //B1. Compute the denominator
724  //-----------------------------------
725 
726  point_value4(Det23,roots,pdet23);
727  point_value4(D0D2,roots,pd0d2);
728  cram11.resize(nsol);
729  for ( int i=1; i<=nsol; i++ ) {
730  cram11[i] = pdet23[i]*pd0d2[i];
731  }
732  point_value6(Det678, roots, cram12);
733  point_value4(Det31, roots, pdet31);
734  point_value2(D[2], roots, pd2);
735  cram21.resize(nsol);
736  for ( int i=1; i<=nsol; i++ ) {
737  cram21[i] = -pdet31[i]*pd2[i];
738  }
739  point_value8(Det3478, roots, cram22);
740  point_value4(Det32, roots, pdet32);
741  cram31.resize(nsol);
742  for ( int i=1; i<=nsol; i++ ) {
743  cram31[i] = -pdet32[i]*pd2[i];
744  }
745  point_value8(Det3468, roots, cram32);
746  point_value2(D[3], roots, pd3);
747  cram41.resize(nsol);
748  for ( int i=1; i<=nsol; i++ ) {
749  cram41[i] = -pdet31[i]*pd3[i];
750  }
751  point_value8(Det2478, roots, cram42);
752  cram51.resize(nsol);
753  for ( int i=1; i<=nsol; i++ ) {
754  cram51[i] = -pdet32[i]*pd3[i];
755  }
756  point_value8(Det1357, roots, cram52);
757  crden.resize(nsol);
758  for ( int i=1; i<=nsol; i++ ) {
759  crden[i] = -cram11[i]*cram12[i] + cram21[i]*cram22[i] - cram31[i]*cram32[i] - cram41[i]*cram42[i] + cram51[i]*cram52[i];
760  }
761 
762  //-------------------------------------------------------------
763  //B2. Compute the numerator for tau[1] and store tau[1] values
764  //-------------------------------------------------------------
765 
766  point_value4(Det22, roots, pdet22);
767  crn111.resize(nsol);
768  for ( int i=1; i<=nsol; i++ ) {
769  crn111[i] = pdet22[i]*pd0d2[i];
770  }
771  point_value6(Det678, roots, crn112);
772  point_value2(D[1], roots, pd1);
773  crn121.resize(nsol);
774  for ( int i=1; i<=nsol; i++ ) {
775  crn121[i] = -pdet31[i]*pd1[i];
776  }
777  point_value8(Det3478, roots, crn122);
778  crn131.resize(nsol);
779  for ( int i=1; i<=nsol; i++ ) {
780  crn131[i] = -pdet32[i]*pd1[i];
781  }
782  point_value8(Det3468, roots, crn132);
783  cram_num1.resize(nsol);
784  for ( int i=1; i<=nsol; i++ ) {
785  cram_num1[i] = -crn111[i]*crn112[i]+crn121[i]*crn122[i]-crn131[i]*crn132[i];
786  tau[i][order[1]] = -cram_num1[i] / crden[i]; // restoring column-major indexing!
787  }
788 
789  //-------------------------------------------------------------
790  //B3. Compute the numerator for tau[2] and store tau[2] values
791  //-------------------------------------------------------------
792 
793  point_value4(Det21, roots, pdet21);
794  crn221.resize(nsol);
795  for ( int i=1; i<=nsol; i++ ) {
796  crn221[i] = -pdet21[i]*pd2[i];
797  }
798  point_value8(Det3478, roots, crn222);
799  crn231.resize(nsol);
800  for ( int i=1; i<=nsol; i++ ) {
801  crn231[i] = -pdet21[i]*pd3[i];
802  }
803  point_value8(Det3468, roots, crn232);
804  crn241.resize(nsol);
805  for ( int i=1; i<=nsol; i++ ) {
806  crn241[i] = -pdet22[i]*pd2[i];
807  }
808  point_value8(Det2478, roots, crn242);
809  crn251.resize(nsol);
810  for ( int i=1; i<=nsol; i++ ) {
811  crn251[i] = -pdet22[i]*pd3[i];
812  }
813  point_value8(Det1357, roots, crn252);
814  point_value4(Det33, roots, pdet33);
815  crn261.resize(nsol);
816  for ( int i=1; i<=nsol; i++ ) {
817  crn261[i] = -pdet33[i]*pd0d2[i];
818  }
819  point_value6(Det123, roots, crn262);
820  cram_num2.resize(nsol);
821  for ( int i=1; i<=nsol; i++ ) {
822  cram_num2[i] = crn221[i]*crn222[i] - crn231[i]*crn232[i] - crn241[i]*crn242[i] + crn251[i]*crn252[i] - crn261[i]*crn262[i];
823  tau[i][order[2]] = -cram_num2[i] / crden[i]; // DJM: division inside for-loop!
824  } // and restoring column-major indexing!
825 
826  //----------------------------------------
827  //B4. Compute the output sines and cosines
828  //----------------------------------------
829  tsq.resize(nsol);
830  tsq1.resize(nsol);
831  cos.resize(nsol);
832  sin.resize(nsol);
833  for ( int i=1; i<=nsol; i++ ) {
834  tsq[i].resize(3);
835  tsq1[i].resize(3);
836  cos[i].resize(3);
837  sin[i].resize(3);
838  for ( int j=1; j<=3; j++ ) {
839  tsq[i][j]=(tau[i][j])*(tau[i][j]);
840  tsq1[i][j]=tsq[i][j]+1;
841  cos[i][j]=(1-tsq[i][j]) / tsq1[i][j]; // DJM: division inside
842  sin[i][j]=2 * tau[i][j] / tsq1[i][j]; // of for-loops!
843  }
844  }
845 
846  return;
847 }
848 // }}}1
849 
850 // test_point_value2 {{{1
851 void test_point_value2() {
852  utility::vector1<Real> A(3), t(16), C;
853  A[1]=0.0;
854  A[2]=-60;
855  A[3]=0.0;
856  t[1]=7.4248;
857  t[2]=5.3100;
858  t[3]=0.0;
859  t[4]=0.0;
860  t[5]=0.0;
861  t[6]=0.0;
862  t[7]=0.0;
863  t[8]=0.0;
864  t[9]=0.0;
865  t[10]=0.0;
866  t[11]=0.0;
867  t[12]=0.0;
868  t[13]=0.0;
869  t[14]=0.0;
870  t[15]=0.0;
871  t[16]=0.0;
872  point_value2(A,t,C);
873  printVector(C);
874  // output should be { -445.4891, -318.6002, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,0.0,0.0, 0.0, 0.0}
875 }
876 
877 // test_polyProduct6x6 {{{1
878 void test_polyProduct6x6() {
879  utility::vector1<Real> A (7), B (7), C;
880  A[1]=366130000.0;
881  A[2]=0.0;
882  A[3]=316600000.0;
883  A[4]=0.0;
884  A[5]=-47950000.0;
885  A[6]=0.0;
886  A[7]=920000.0;
887  B[1]=1790200000.0;
888  B[2]=0.0;
889  B[3]=2113200000.0;
890  B[4]=0.0;
891  B[5]=317100000.0;
892  B[6]=0.0;
893  B[7]=-8300000.0;
894  polyProduct6x6(A,B,C);
895  printVector(C);
896  // output should be 1.0e+18*{0.6554, 0.0, 1.3405, 0.0, 0.6993, 0.0, -0.0023, 0.0, -0.0159, 0.0}, 6.89717e+14, 0.0, -7.636e+12
897  return;
898 }
899 
900 // test_polyProduct4sq {{{1
901 void test_polyProduct4sq() {
902  utility::vector1<Real> A (5), C;
903  A[1]=126600.0;
904  A[2]=0.0;
905  A[3]=-28200.0;
906  A[4]=0.0;
907  A[5]=-2890.0;
908  polyProduct4sq(A,C);
909  printVector(C);
910  // output should be 1.0e+10*{1.6027, 0.0, -0.7140, 0.0, 0.0063, 0.0, 0.0163, 0.0, 0.0008}
911  return;
912 }
913 
914 // test_polyProduct4x4 {{{1
915 void test_polyProduct4x4() {
916  utility::vector1<Real> A (5), B (5), C;
917  A[1]=126940.0;
918  A[2]=0.0;
919  A[3]=-28110.0;
920  A[4]=0.0;
921  A[5]=-2890.0;
922  B[1]=0.0;
923  B[2]=172730.0;
924  B[3]=0.0;
925  B[4]=20800.0;
926  B[5]=0.0;
927  polyProduct4x4(A,B,C);
928  printVector(C);
929  // output should be 1.0e+10*{0.0, 2.1927, 0.0, -0.2216, 0.0, -0.1085, 0.0, -0.0060, 0.0}
930  return;
931 }
932 
933 // test_polyProduct4x2 {{{1
934 void test_polyProduct4x2() {
935  utility::vector1<Real> A (5), B (3), C;
936  A[1]=0.0;
937  A[2]=-108950.0;
938  A[3]=0.0;
939  A[4]=43180.0;
940  A[5]=0.0;
941  B[1]=0.0;
942  B[2]=-1799.10;
943  B[3]=0.0;
944  polyProduct4x2(A,B,C);
945  printVector(C);
946  // output should be 1.0e+08 * {0.0, 0.0, 1.9601, 0.0, -0.7769, 0.0, 0.0}
947  return;
948 }
949 
950 // test_polyProduct2x2 {{{1
951 void test_polyProduct2x2() {
952  utility::vector1<Real> A (3), B (3), C;
953  A[1]=15.999410431075338;
954  A[2]=0.0;
955  A[3]=-8.999904923464781;
956  B[1]=72.015965384650187;
957  B[2]=0.0;
958  B[3]=7.000641402062060;
959  polyProduct2x2(A,B,C);
960  printVector(C);
961  // output should be 1.03e+03 * {1.152212987779133, 0.0, -0.536130706361013, 0.0, -0.063005107021830}
962  return;
963 }
964 
965 // test_dixon() {{{1
966 void test_dixon() {
967  Real Avals[]={-0.061472096577788, -2.424098126176366, -0.016332287075674, 2.421954432062721, -0.090242315482933, 0.643479843225995, 0.061472096577788, -0.646358402699405, 0.016332287075674};
968  Real Bvals[]={-0.000625561301665, 0.0, 1.070810990722253, 0.054344479784912, 0.0, 0.054308028868740, -1.071191403371629, 0.0, 0.000963218011333};
969  Real Cvals[]={0.063394336809916, 2.498922707477288, -0.025585174231195, -2.497689253369826, 0.177963037883907, 1.008036647097698, -0.063394336809916, -1.006882302775246, 0.025585174231195};
970  Real Dvals[]={0.673457026339820, -0.015821144672541, 1.005285612756159, -0.078123031130501, 3.277064399971811, 0.078123031130501, 1.004894243686905, 0.075272539656448, -0.573852487168651};
971  //utility::vector1<utility::vector1<utility::vector1<Real> > > P(3); // dixon input
972  utility::vector1<utility::vector1<Real> > A (3);
973  utility::vector1<utility::vector1<Real> > B (3);
974  utility::vector1<utility::vector1<Real> > C (3);
975  utility::vector1<utility::vector1<Real> > D (3);
976  utility::vector1<int> order (3); // dixon input
977  utility::vector1<utility::vector1<Real> > tau; // dixon output
978  utility::vector1<utility::vector1<Real> > cos; // dixon output
979  utility::vector1<utility::vector1<Real> > sin; // dixon output
980  int nsol = 0;
981 
982  // Allocate A,B,C,D
983  for ( int i=1; i<=3; i++ ) {
984  A[i].resize(3);
985  B[i].resize(3);
986  C[i].resize(3);
987  D[i].resize(3);
988  }
989 
990  // Fill A,B,C,D
991  int n=0;
992  for ( int i=1; i<=3; i++ ) {
993  for ( int j=1; j<=3; j++ ) {
994  A[i][j]=Avals[n];
995  B[i][j]=Bvals[n];
996  C[i][j]=Cvals[n];
997  D[i][j]=Dvals[n];
998  n++;
999  }
1000  }
1001 
1002  // initialize order
1003  order[1]=1;
1004  order[2]=2;
1005  order[3]=3;
1006 
1007  // compute Dixon resultant
1008  //dixon(A, B, C, D, order, cos, sin, tau, nsol);
1009 
1010  cout << "number of solutions" << nsol << std::endl;
1011  printMatrix(cos);
1012  printMatrix(sin);
1013  return;
1014 }
1015 // }}}1
1016 
1017 // main (commented out) {{{1
1018 // int main(int argc, char** argv)
1019 // {
1020 // //test_polyProduct2x2();
1021 // //test_polyProduct4x2();
1022 // //test_polyProduct4x4();
1023 // //test_polyProduct4sq();
1024 // //test_polyProduct6x6();
1025 // //test_point_value2();
1026 // // for (int i=1; i<10000; i++) {
1027 // test_dixon();
1028 // //}
1029 // }
1030 // }}}1
1031 
1032 } // end namespace kinematic_closure
1033 } // end namespace numeric
numeric::kinematic_closure::initialize_sturm
void initialize_sturm(double *tol_secant, int *max_iter_sturm, int *max_iter_secant)
Definition: sturm.hh:54
numeric::kinematic_closure::PseudoMatrix
vector1< PseudoVector > PseudoMatrix
Definition: dixon.cc:62
numeric::kinematic_closure::test_polyProduct4x2
void test_polyProduct4x2()
Definition: dixon.cc:934
numeric::kinematic_closure::point_value8
void point_value8(const utility::vector1< Real > &A, const utility::vector1< Real > &t, utility::vector1< Real > &C)
Definition: dixon.cc:95
numeric::kinematic_closure::polyProduct4sq
void polyProduct4sq(const utility::vector1< Real > &A, utility::vector1< Real > &C)
Definition: dixon.cc:156
numeric::kinematic_closure::polyProduct2x2
void polyProduct2x2(const utility::vector1< Real > &A, const utility::vector1< Real > &B, utility::vector1< Real > &C)
Definition: dixon.cc:114
numeric::kinematic_closure::Matrix8
Eigen::Matrix< Real, 8, 8 > Matrix8
Definition: dixon.cc:50
numeric::kinematic_closure::test_point_value2
void test_point_value2()
Definition: dixon.cc:851
numeric::kinematic_closure::solve_sturm
void solve_sturm(const int &p_order, int &n_root, const utility::vector1< double > &poly_coeffs, utility::vector1< double > &roots)
Definition: sturm.hh:491
numeric::kinematic_closure::polyProduct6x6
void polyProduct6x6(const utility::vector1< Real > &A, const utility::vector1< Real > &B, utility::vector1< Real > &C)
Definition: dixon.cc:172
GeneralizedEigenSolver.hh
Header file for GeneralizedEigenSolver.
numeric::kinematic_closure::Vector16
Eigen::Matrix< Real, Eigen::Dynamic, 1, 0, 16, 1 > Vector16
Definition: dixon.cc:52
numeric::linear_algebra::GeneralizedEigenSolver::num_real_solutions
int num_real_solutions() const
Return the number of real solutions that were calculated.
Definition: GeneralizedEigenSolver.hh:219
sturm.hh
implements sturm chain solver
numeric::kinematic_closure::polyProduct8x8
void polyProduct8x8(const utility::vector1< Real > &A, const utility::vector1< Real > &B, utility::vector1< Real > &C)
Definition: dixon.cc:216
numeric::kinematic_closure::test_dixon
void test_dixon()
Definition: dixon.cc:966
numeric::kinematic_closure::test_polyProduct4x4
void test_polyProduct4x4()
Definition: dixon.cc:915
numeric::kinematic_closure::PseudoVector
vector1< Real > PseudoVector
Definition: dixon.cc:61
numeric::kinematic_closure::point_value4
void point_value4(const utility::vector1< Real > &A, const utility::vector1< Real > &t, utility::vector1< Real > &C)
Definition: dixon.cc:77
numeric::linear_algebra::GeneralizedEigenSolver::RealEigenvalueType
Eigen::Matrix< RealScalar, Eigen::Dynamic, 1, Options, MaxColsAtCompileTime, 1 > RealEigenvalueType
Vector type used to represent the real subset of calculated eigenvalues.
Definition: GeneralizedEigenSolver.hh:103
DIXON_SIZE
#define DIXON_SIZE
Definition: dixon.cc:37
numeric::kinematic_closure::polyProduct4x2
void polyProduct4x2(const utility::vector1< Real > &A, const utility::vector1< Real > &B, utility::vector1< Real > &C)
Definition: dixon.cc:126
numeric::kinematic_closure::dixonResultant
void dixonResultant(const utility::vector1< utility::vector1< Real > > &A, const utility::vector1< utility::vector1< Real > > &B, const utility::vector1< utility::vector1< Real > > &C, const utility::vector1< utility::vector1< Real > > &D, utility::vector1< utility::vector1< utility::vector1< Real > > > &R)
Definition: dixon.cc:279
numeric::kinematic_closure::point_value6
void point_value6(const utility::vector1< Real > &B, const utility::vector1< Real > &t, utility::vector1< Real > &C)
Definition: dixon.cc:86
numeric::linear_algebra::GeneralizedEigenSolver::real_eigenvalues
RealEigenvalueType real_eigenvalues() const
Return any real eigenvalues that were calculated.
Definition: GeneralizedEigenSolver.hh:424
numeric::kinematic_closure::vectorDiff
void vectorDiff(const utility::vector1< Real > &A, const utility::vector1< Real > &B, utility::vector1< Real > &C)
Definition: dixon.cc:266
numeric::kinematic_closure::printMatrix
void printMatrix(const utility::vector1< utility::vector1< Real > > &M)
prints the matrix
Definition: kinematic_closure_helpers.cc:42
C
#define C(a, b)
Definition: functions.cc:27
ObjexxFCL::abs
T abs(T const &x)
std::abs( x ) == | x |
Definition: Fmath.hh:295
utility::vector1
std::vector with 1-based indexing
Definition: vector1.fwd.hh:44
dixon.hh
Header file for dixon code for ALC.
DIXON_RESULTANT_SIZE
#define DIXON_RESULTANT_SIZE
Definition: dixon.cc:38
numeric::linear_algebra::GeneralizedEigenSolver
Solves generalized eigenvalue problems.
Definition: GeneralizedEigenSolver.fwd.hh:21
numeric::kinematic_closure::Matrix16
Eigen::Matrix< Real, 16, 16 > Matrix16
Definition: dixon.cc:51
types.hh
rosetta project type declarations. Should be kept updated with core/types.hh. This exists because num...
numeric::kinematic_closure::skipl
string skipl
Definition: dixon.cc:48
numeric::kinematic_closure::dixon_eig
void dixon_eig(PseudoMatrix const &A, PseudoMatrix const &B, PseudoMatrix const &C, PseudoMatrix const &D, vector1< int > const &, PseudoMatrix &cos, PseudoMatrix &sin, PseudoMatrix &u, int &num_solutions)
Definition: dixon.cc:431
numeric::Real
double Real
Definition: types.hh:39
numeric::kinematic_closure::polyProduct4x4
void polyProduct4x4(const utility::vector1< Real > &A, const utility::vector1< Real > &B, utility::vector1< Real > &C)
Definition: dixon.cc:140
numeric::kinematic_closure::point_value2
void point_value2(const utility::vector1< Real > &A, const utility::vector1< Real > &t, utility::vector1< Real > &C)
dixon functions ///
Definition: dixon.cc:68
numeric::kinematic_closure::test_polyProduct6x6
void test_polyProduct6x6()
Definition: dixon.cc:878
utility::io::oc::cout
ocstream cout(std::cout)
Wrapper around std::cout.
Definition: ocstream.hh:287
numeric::kinematic_closure::point_value16
void point_value16(const utility::vector1< Real > &A, const utility::vector1< Real > &t, utility::vector1< Real > &C)
Definition: dixon.cc:104
numeric::kinematic_closure::printVector
void printVector(const utility::vector1< Real > &V)
helper functions ///
Definition: kinematic_closure_helpers.cc:30
numeric::kinematic_closure::scipy
Eigen::IOFormat scipy(8, 0,", ",",\n","[","]","[","]")
numeric::linear_algebra::GeneralizedEigenSolver::real_eigenvectors
RealEigenvectorType real_eigenvectors() const
Return any real eigenvectors that were calculated.
Definition: GeneralizedEigenSolver.hh:450
numeric::kinematic_closure::polyProduct8sq
void polyProduct8sq(const utility::vector1< Real > &A, utility::vector1< Real > &C)
Definition: dixon.cc:240
numeric::kinematic_closure::build_dixon_matrices
void build_dixon_matrices(PseudoMatrix const &A, PseudoMatrix const &B, PseudoMatrix const &C, PseudoMatrix const &D, Matrix8 &R0, Matrix8 &R1, Matrix8 &R2)
Definition: dixon.cc:372
ObjexxFCL::format::I
std::string I(int const w, T const &t)
Integer Format.
Definition: format.hh:758
numeric::kinematic_closure::test_polyProduct2x2
void test_polyProduct2x2()
Definition: dixon.cc:951
numeric::kinematic_closure::SolverType
linear_algebra::GeneralizedEigenSolver< Matrix16 > SolverType
Definition: dixon.cc:53
numeric::linear_algebra::GeneralizedEigenSolver::RealEigenvectorType
Eigen::Matrix< RealScalar, RowsAtCompileTime, Eigen::Dynamic, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime > RealEigenvectorType
Matrix type used to represent the real subset of calculated eigenvectors.
Definition: GeneralizedEigenSolver.hh:108
kinematic_closure_helpers.hh
Header file for kinematic_closure_helpers.cc.
ObjexxFCL::format::A
std::string A(int const w, char const c)
char Format
Definition: format.cc:175
numeric::kinematic_closure::polyProduct12x4
void polyProduct12x4(const utility::vector1< Real > &A, const utility::vector1< Real > &B, utility::vector1< Real > &C)
Definition: dixon.cc:192
numeric::kinematic_closure::test_polyProduct4sq
void test_polyProduct4sq()
Definition: dixon.cc:901
numeric::kinematic_closure::build_sin_and_cos
void build_sin_and_cos(PseudoMatrix const &u, PseudoMatrix &sin, PseudoMatrix &cos)
Definition: dixon.cc:405
PP4x2_VECSIZE
#define PP4x2_VECSIZE
Definition: dixon.cc:36
numeric::kinematic_closure::dixon_sturm
void dixon_sturm(const utility::vector1< utility::vector1< Real > > &A, const utility::vector1< utility::vector1< Real > > &B, const utility::vector1< utility::vector1< Real > > &C, const utility::vector1< utility::vector1< Real > > &D, const utility::vector1< int > &order, utility::vector1< utility::vector1< Real > > &cos, utility::vector1< utility::vector1< Real > > &sin, utility::vector1< utility::vector1< Real > > &tau, int &nsol)
Definition: dixon.cc:492
Generated on Sun Jan 17 2016 15:21:19 for Rosetta by   doxygen 1.8.7