Skip to content

Commit

Permalink
Debug 3rd order constraint
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielDoehring
DanielDoehring committed May 30, 2023
1 parent 3fbe5e7 commit 2d385ba
Show file tree
Hide file tree
Showing 4 changed files with 120 additions and 48 deletions.
42 changes: 30 additions & 12 deletions Feasibility_Problem/src/OrderConstraints_Real.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,7 +185,9 @@ void ThirdOrder(const T* x, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand Down Expand Up @@ -224,14 +226,16 @@ void ThirdOrder(const T* x, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

g[NumEigVals+1] -= x[j] * x[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -243,7 +247,9 @@ void ThirdOrder(const T* x, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -266,7 +272,9 @@ void ThirdOrder(const std::vector<T>& x, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand Down Expand Up @@ -305,14 +313,16 @@ void ThirdOrder(const std::vector<T>& x, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

g[NumEigVals+1] -= x[j] * x[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -324,7 +334,9 @@ void ThirdOrder(const std::vector<T>& x, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -347,7 +359,9 @@ T ThirdOrder(const std::vector<T>& x, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand Down Expand Up @@ -388,14 +402,16 @@ T ThirdOrder(const std::vector<T>& x, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

g -= x[j] * x[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -407,7 +423,9 @@ T ThirdOrder(const std::vector<T>& x, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand Down
42 changes: 30 additions & 12 deletions Feasibility_Problem/src/OrderConstraints_RealImag.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -184,7 +184,9 @@ void ThirdOrder(const T* xy, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand Down Expand Up @@ -223,14 +225,16 @@ void ThirdOrder(const T* xy, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

g[NumEigVals+1] -= xy[j] * xy[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand All @@ -242,7 +246,9 @@ void ThirdOrder(const T* xy, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand All @@ -265,7 +271,9 @@ void ThirdOrder(const std::vector<T>& xy, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand Down Expand Up @@ -304,14 +312,16 @@ void ThirdOrder(const std::vector<T>& xy, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

g[NumEigVals+1] -= xy[j] * xy[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand All @@ -323,7 +333,9 @@ void ThirdOrder(const std::vector<T>& xy, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand All @@ -346,7 +358,9 @@ T ThirdOrder(const std::vector<T>& xy, const int NumRoots,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand Down Expand Up @@ -387,14 +401,16 @@ T ThirdOrder(const std::vector<T>& xy, const int NumRoots,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

g -= xy[j] * xy[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand All @@ -406,7 +422,9 @@ T ThirdOrder(const std::vector<T>& xy, const int NumRoots,
b1 = Lin_IntPol(xy[j], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[j + NumRoots];
Radius1 = xy[j]*xy[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(xy[k], RealRange, ImagRange, ImagDiff_over_RealDiff) + xy[k + NumRoots];
Radius2 = xy[k]*xy[k] + b2*b2;

Expand Down
42 changes: 30 additions & 12 deletions Optimization_Problem/src/OrderConstraints_Real.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -181,7 +181,9 @@ void ThirdOrder(const T* x, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(x[j], RealRange, ImagRange);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange);
Radius2 = x[k]*x[k] + b2*b2;

Expand Down Expand Up @@ -220,14 +222,16 @@ void ThirdOrder(const T* x, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

g[NumEigVals+1] -= x[j] * x[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -239,7 +243,9 @@ void ThirdOrder(const T* x, T* g, const int NumRoots, const int NumEigVals,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -261,7 +267,9 @@ void ThirdOrder(const std::vector<T>& x, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange);
Radius2 = x[k]*x[k] + b2*b2;

Expand Down Expand Up @@ -300,14 +308,16 @@ void ThirdOrder(const std::vector<T>& x, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

g[NumEigVals+1] -= x[j] * x[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -319,7 +329,9 @@ void ThirdOrder(const std::vector<T>& x, std::vector<T>& g, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g[NumEigVals+1] -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -341,7 +353,9 @@ T ThirdOrder(const std::vector<T>& x, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange);
Radius2 = x[k]*x[k] + b2*b2;

Expand Down Expand Up @@ -382,14 +396,16 @@ T ThirdOrder(const std::vector<T>& x, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < i_min; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < i_min; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

g -= x[j] * x[k] / (Radius1 * Radius2);
}

for(size_t k = i_min+1; j < NumRoots; j++) {
for(size_t k = i_min+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand All @@ -401,7 +417,9 @@ T ThirdOrder(const std::vector<T>& x, const int NumRoots,
b1 = Lin_IntPol(x[j], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius1 = x[j]*x[j] + b1*b1;

for(size_t k = j+1; j < NumRoots; j++) {
g -= 0.25 / Radius1;

for(size_t k = j+1; k < NumRoots; k++) {
b2 = Lin_IntPol(x[k], RealRange, ImagRange, ImagDiff_over_RealDiff);
Radius2 = x[k]*x[k] + b2*b2;

Expand Down
Loading

0 comments on commit 2d385ba

Please sign in to comment.