1: static char help[] = "test least-squares problem created from a mapped taoterm quadratic";

  3: #include <petsctao.h>

  5: int main(int argc, char **argv)
  6: {
  7:   MPI_Comm    comm;
  8:   Mat         A;       // data matrix
  9:   Mat         W;       // weight matrix
 10:   Vec         w;       // observation vector
 11:   Vec         b;       // observation vector
 12:   PetscInt    m = 100; // data size
 13:   PetscInt    n = 20;  // model size
 14:   TaoTerm     data_term;
 15:   PetscRandom rand;
 16:   Tao         tao;
 17:   PetscInt    i, j;
 18:   PetscReal   val, density = 0.3;
 19:   PetscBool   test_quad_mat = PETSC_FALSE;

 21:   PetscFunctionBeginUser;
 22:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 23:   comm = PETSC_COMM_WORLD;

 25:   PetscOptionsBegin(comm, "", help, "none");
 26:   PetscCall(PetscOptionsBoundedInt("-m", "data size", "", m, &m, NULL, 0));
 27:   PetscCall(PetscOptionsBoundedInt("-n", "model size", "", n, &n, NULL, 0));
 28:   PetscCall(PetscOptionsBool("-test_quad_mat", "Test if quadratic term matrix matches W matrix", "", test_quad_mat, &test_quad_mat, NULL));
 29:   PetscOptionsEnd();

 31:   PetscCall(TaoCreate(comm, &tao));

 33:   PetscCall(PetscRandomCreate(comm, &rand));
 34:   PetscCall(PetscRandomSetInterval(rand, -1.0, 1.0));
 35:   PetscCall(PetscRandomSetFromOptions(rand));

 37:   // create the model data, A, W and b
 38:   PetscCall(MatCreate(comm, &A));
 39:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m, n));
 40:   PetscCall(MatSetType(A, MATAIJ));
 41:   PetscCall(MatSetFromOptions(A));
 42:   PetscCall(MatSetUp(A));
 43:   for (i = 0; i < m; i++) {
 44:     for (j = 0; j < n; j++) {
 45:       PetscCall(PetscRandomGetValueReal(rand, &val));
 46:       // Optionally make it sparse: only insert some entries
 47:       if (val < density) PetscCall(MatSetValue(A, i, j, val, INSERT_VALUES));
 48:     }
 49:   }
 50:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 51:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 52:   PetscCall(VecCreateMPI(comm, PETSC_DECIDE, m, &b));
 53:   PetscCall(VecSetRandom(b, rand));
 54:   PetscCall(VecDuplicate(b, &w));
 55:   PetscCall(VecSetRandom(w, rand));
 56:   PetscCall(VecAbs(w));
 57:   PetscCall(VecShift(w, 1.0));
 58:   PetscCall(MatCreateDiagonal(w, &W));
 59:   PetscCall(VecDestroy(&w));

 61:   // the model term,  (1/2) || Ax - b ||_W^2
 62:   PetscCall(TaoTermCreateQuadratic(W, &data_term));
 63:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)data_term, "data_"));
 64:   PetscCall(TaoAddTerm(tao, "data_", 3.0, data_term, b, A));
 65:   PetscCall(TaoTermDestroy(&data_term));

 67:   PetscCall(TaoSetFromOptions(tao));
 68:   PetscCall(TaoSolve(tao));

 70:   if (test_quad_mat) {
 71:     PetscReal   scale;
 72:     TaoTerm     term;
 73:     Vec         params;
 74:     Mat         map;
 75:     Mat         quad_mat;
 76:     TaoTermType term_type;
 77:     PetscBool   is_quad, mat_equal;

 79:     PetscCall(TaoGetTerm(tao, &scale, &term, &params, &map));
 80:     PetscCall(TaoTermGetType(term, &term_type));
 81:     PetscCall(PetscStrcmp(term_type, TAOTERMQUADRATIC, &is_quad));
 82:     PetscCheck(is_quad, comm, PETSC_ERR_ARG_WRONG, "Term from TaoGetTerm is not a quadratic term");

 84:     PetscCall(TaoTermQuadraticGetMat(term, &quad_mat));
 85:     PetscCheck(quad_mat != NULL, comm, PETSC_ERR_ARG_NULL, "Quadratic term matrix is NULL");

 87:     PetscCall(MatEqual(W, quad_mat, &mat_equal));
 88:     PetscCheck(mat_equal, comm, PETSC_ERR_PLIB, "Quadratic term matrix does not match W matrix");
 89:     PetscCall(PetscPrintf(comm, "Test passed: Quadratic term matrix matches W matrix\n"));
 90:   }

 92:   PetscCall(VecDestroy(&b));
 93:   PetscCall(MatDestroy(&W));
 94:   PetscCall(MatDestroy(&A));
 95:   PetscCall(PetscRandomDestroy(&rand));
 96:   PetscCall(TaoDestroy(&tao));
 97:   PetscCall(PetscFinalize());
 98:   return 0;
 99: }

101: /*TEST

103:   build:
104:     requires: !complex !single !quad !defined(PETSC_USE_64BIT_INDICES) !__float128

106:   test:
107:     suffix: 0
108:     args: -tao_monitor_short -tao_view -tao_type nls

110:   test:
111:     suffix: 1
112:     args: -tao_view ::ascii_info_detail -tao_type nls

114:   test:
115:     suffix: test_quad_mat
116:     args: -test_quad_mat 1

118: TEST*/