Bugfix at order 3 (hopefully was not affecting final results, thanks to Pablo Winant for spotting this):

- in the 3rd derivatives matrix, among symmetric elements, the first one had
the right values but the following ones were set to zero
- moreover, the k-order DLL was trying to add all the symmetric elements in the
folded tensor, instead of only keeping one value among all the symmetric ones
- hopefully, Ondra's tensor library was (silently) refusing to add symmetric
elements after the first (and right) value had been added
- so the final result was correct
time-shift
Sébastien Villemot 2010-12-06 17:00:46 +01:00
parent 0d1fd700c7
commit a9b4331a95
2 changed files with 4 additions and 24 deletions

View File

@ -211,27 +211,9 @@ KordpDynare::populateDerivativesContainer(const TwoDMatrix &g, int ord, const ve
s[2] = revOrder[s2];
else
s[2] = s2;
// reordering symmetry indices in increasing order
if (s[0] > s[1])
{
int temp = s[0];
s[0] = s[1];
s[1] = temp;
}
if (s[1] > s[2])
{
int temp = s[1];
s[1] = s[2];
s[2] = temp;
}
if (s[0] > s[1])
{
int temp = s[0];
s[0] = s[1];
s[1] = temp;
}
double x = g.get(i, 2);
mdTi->insert(s, j, x);
double x = g.get(i, 2);
if (s.isSorted())
mdTi->insert(s, j, x);
}
}

View File

@ -2137,8 +2137,6 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
d3->writeOutput(third_derivatives_output, output_type, temporary_terms);
third_derivatives_output << ";" << endl;
k++;
// Compute the column numbers for the 5 other permutations of (id1,id2,id3) and store them in a set (to avoid duplicates if two indexes are equal)
set<int> cols;
cols.insert(id1 * hessianColsNbr + id3 * dynJacobianColsNbr + id2);
@ -2147,7 +2145,7 @@ DynamicModel::writeDynamicModel(ostream &DynamicOutput, bool use_dll) const
cols.insert(id3 * hessianColsNbr + id1 * dynJacobianColsNbr + id2);
cols.insert(id3 * hessianColsNbr + id2 * dynJacobianColsNbr + id1);
int k2 = 0; // Keeps the offset of the permutation relative to k
int k2 = 1; // Keeps the offset of the permutation relative to k
for (set<int>::iterator it2 = cols.begin(); it2 != cols.end(); it2++)
if (*it2 != ref_col)
{