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:   PetscReal   val, density = 0.3;
 18:   PetscBool   test_quad_mat = PETSC_FALSE;

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

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

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

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

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

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

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

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

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

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

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

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

100: /*TEST

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

105:   test:
106:     suffix: 0
107:     args: -tao_monitor -tao_view -tao_type nls

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

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

117: TEST*/