1: static char help[] = "Tests TAO utilities with one-sided variable bounds.\n";

  3: #include <petsctao.h>

  5: static PetscErrorCode SetSolution(Vec x)
  6: {
  7:   PetscScalar *values;
  8:   PetscInt     i, rstart, rend;

 10:   PetscFunctionBeginUser;
 11:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
 12:   PetscCall(VecGetArray(x, &values));
 13:   for (i = rstart; i < rend; ++i) values[i - rstart] = i == 0 ? -2.0 : (i == 3 ? 2.0 : 0.0);
 14:   PetscCall(VecRestoreArray(x, &values));
 15:   PetscFunctionReturn(PETSC_SUCCESS);
 16: }

 18: static PetscErrorCode CheckSolution(Vec x, PetscBool lower)
 19: {
 20:   const PetscScalar *values;
 21:   PetscInt           i, rstart, rend;

 23:   PetscFunctionBeginUser;
 24:   PetscCall(VecGetOwnershipRange(x, &rstart, &rend));
 25:   PetscCall(VecGetArrayRead(x, &values));
 26:   for (i = rstart; i < rend; ++i) {
 27:     PetscScalar expected = i == 0 ? (lower ? -1.0 : -2.0) : (i == 3 ? (lower ? 2.0 : 1.0) : 0.0);

 29:     PetscCheck(values[i - rstart] == expected, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Unexpected value at index %" PetscInt_FMT, i);
 30:   }
 31:   PetscCall(VecRestoreArrayRead(x, &values));
 32:   PetscFunctionReturn(PETSC_SUCCESS);
 33: }

 35: static PetscErrorCode TestActiveBounds(Vec x, Vec bound, Vec g, Vec s, Vec work, PetscBool lower)
 36: {
 37:   IS        active_lower = NULL, active_upper = NULL, active_fixed = NULL, active = NULL, inactive = NULL;
 38:   IS        active_bound;
 39:   PetscReal bound_tol = 0.0;
 40:   PetscInt  n;

 42:   PetscFunctionBeginUser;
 43:   PetscCall(SetSolution(x));
 44:   PetscCall(VecSet(g, lower ? 1.0 : -1.0));
 45:   PetscCall(VecSet(s, 0.0));
 46:   PetscCall(TaoEstimateActiveBounds(x, lower ? bound : NULL, lower ? NULL : bound, g, s, work, 1.0, &bound_tol, &active_lower, &active_upper, &active_fixed, &active, &inactive));
 47:   active_bound = lower ? active_lower : active_upper;
 48:   PetscCheck(active_bound, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected an active bound");
 49:   PetscCall(ISGetSize(active_bound, &n));
 50:   PetscCheck(n == 1, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected one active bound, got %" PetscInt_FMT, n);
 51:   PetscCheck(!(lower ? active_upper : active_lower), PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Unexpected active bound on the unbounded side");
 52:   PetscCall(ISDestroy(&active_lower));
 53:   PetscCall(ISDestroy(&active_upper));
 54:   PetscCall(ISDestroy(&active_fixed));
 55:   PetscCall(ISDestroy(&active));
 56:   PetscCall(ISDestroy(&inactive));
 57:   PetscFunctionReturn(PETSC_SUCCESS);
 58: }

 60: int main(int argc, char **argv)
 61: {
 62:   Vec      x, xl, xu, g, s, work;
 63:   PetscInt nDiff;

 65:   PetscFunctionBeginUser;
 66:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 67:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 68:   PetscCall(VecSetSizes(x, PETSC_DECIDE, 4));
 69:   PetscCall(VecSetFromOptions(x));
 70:   PetscCall(VecDuplicate(x, &xl));
 71:   PetscCall(VecDuplicate(x, &xu));
 72:   PetscCall(VecDuplicate(x, &g));
 73:   PetscCall(VecDuplicate(x, &s));
 74:   PetscCall(VecDuplicate(x, &work));
 75:   PetscCall(VecSet(xl, -1.0));
 76:   PetscCall(VecSet(xu, 1.0));

 78:   PetscCall(SetSolution(x));
 79:   PetscCall(TaoBoundSolution(x, xl, NULL, 0.0, &nDiff, x));
 80:   PetscCheck(nDiff == 1, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected one lower-bound correction, got %" PetscInt_FMT, nDiff);
 81:   PetscCall(CheckSolution(x, PETSC_TRUE));

 83:   PetscCall(SetSolution(x));
 84:   PetscCall(TaoBoundSolution(x, NULL, xu, 0.0, &nDiff, x));
 85:   PetscCheck(nDiff == 1, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Expected one upper-bound correction, got %" PetscInt_FMT, nDiff);
 86:   PetscCall(CheckSolution(x, PETSC_FALSE));

 88:   PetscCall(TestActiveBounds(x, xl, g, s, work, PETSC_TRUE));
 89:   PetscCall(TestActiveBounds(x, xu, g, s, work, PETSC_FALSE));

 91:   PetscCall(VecDestroy(&x));
 92:   PetscCall(VecDestroy(&xl));
 93:   PetscCall(VecDestroy(&xu));
 94:   PetscCall(VecDestroy(&g));
 95:   PetscCall(VecDestroy(&s));
 96:   PetscCall(VecDestroy(&work));
 97:   PetscCall(PetscFinalize());
 98:   return 0;
 99: }

101: /*TEST

103:   build:
104:     requires: !complex

106:   test:
107:     output_file: output/empty.out

109:   test:
110:     suffix: 2
111:     nsize: 2
112:     output_file: output/empty.out

114: TEST*/