Actual source code: ex102.c

  1: static char help[] = "Tests MatCreateLRC()\n\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Vec       x, b, c = NULL;
  8:   Mat       A, U, V, LR, X, LRe;
  9:   PetscInt  M = 5, N = 7;
 10:   PetscBool flg;

 12:   PetscFunctionBeginUser;
 13:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 14:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &M, NULL));
 15:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));

 17:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 18:          Create the sparse matrix
 19:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 20:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 21:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
 22:   PetscCall(MatSetOptionsPrefix(A, "A_"));
 23:   PetscCall(MatSetFromOptions(A));
 24:   PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
 25:   PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
 26:   PetscCall(MatSetUp(A));
 27:   PetscCall(MatSetRandom(A, NULL));

 29:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 30:          Create the dense matrices
 31:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 32:   PetscCall(MatCreate(PETSC_COMM_WORLD, &U));
 33:   PetscCall(MatSetSizes(U, PETSC_DECIDE, PETSC_DECIDE, M, 3));
 34:   PetscCall(MatSetType(U, MATDENSE));
 35:   PetscCall(MatSetOptionsPrefix(U, "U_"));
 36:   PetscCall(MatSetFromOptions(U));
 37:   PetscCall(MatSetUp(U));
 38:   PetscCall(MatSetRandom(U, NULL));

 40:   PetscCall(MatCreate(PETSC_COMM_WORLD, &V));
 41:   PetscCall(MatSetSizes(V, PETSC_DECIDE, PETSC_DECIDE, N, 3));
 42:   PetscCall(MatSetType(V, MATDENSE));
 43:   PetscCall(MatSetOptionsPrefix(V, "V_"));
 44:   PetscCall(MatSetFromOptions(V));
 45:   PetscCall(MatSetUp(V));
 46:   PetscCall(MatSetRandom(V, NULL));

 48:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 49:          Create a vector to hold the diagonal of C
 50:          A sequential vector can be created as well on each process
 51:          It is user responsibility to ensure the data in the vector
 52:          is consistent across processors
 53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 54:   PetscCall(PetscOptionsHasName(NULL, NULL, "-use_c", &flg));
 55:   if (flg) {
 56:     PetscCall(VecCreateFromOptions(PETSC_COMM_WORLD, NULL, 1, PETSC_DECIDE, 3, &c));
 57:     PetscCall(VecSetRandom(c, NULL));
 58:   }

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:          Create low rank correction matrix
 62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 63:   PetscCall(PetscOptionsHasName(NULL, NULL, "-low_rank", &flg));
 64:   if (flg) {
 65:     /* create a low-rank matrix, with no A-matrix */
 66:     PetscCall(MatCreateLRC(NULL, U, c, V, &LR));
 67:     PetscCall(MatDestroy(&A));
 68:   } else {
 69:     PetscCall(MatCreateLRC(A, U, c, V, &LR));
 70:   }

 72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 73:          Create the low rank correction matrix explicitly to check for
 74:          correctness
 75:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 76:   PetscCall(MatHermitianTranspose(V, MAT_INITIAL_MATRIX, &X));
 77:   PetscCall(MatDiagonalScale(X, c, NULL));
 78:   PetscCall(MatMatMult(U, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &LRe));
 79:   PetscCall(MatDestroy(&X));
 80:   if (A) PetscCall(MatAYPX(LRe, 1.0, A, DIFFERENT_NONZERO_PATTERN));

 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:          Create test vectors
 84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 85:   PetscCall(MatCreateVecs(LR, &x, &b));
 86:   PetscCall(VecSetRandom(x, NULL));
 87:   PetscCall(MatMult(LR, x, b));
 88:   PetscCall(MatMultTranspose(LR, b, x));
 89:   PetscCall(VecDestroy(&x));
 90:   PetscCall(VecDestroy(&b));

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 93:          Check correctness
 94:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 95:   PetscCall(MatMultEqual(LR, LRe, 10, &flg));
 96:   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error in MatMult\n"));
 97: #if !defined(PETSC_USE_COMPLEX)
 98:   PetscCall(MatMultHermitianTransposeEqual(LR, LRe, 10, &flg));
 99:   if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error in MatMultTranspose\n"));
100: #endif

102:   PetscCall(MatDestroy(&A));
103:   PetscCall(MatDestroy(&LRe));
104:   PetscCall(MatDestroy(&U));
105:   PetscCall(MatDestroy(&V));
106:   PetscCall(VecDestroy(&c));
107:   PetscCall(MatDestroy(&LR));

109:   /*
110:      Always call PetscFinalize() before exiting a program.  This routine
111:        - finalizes the PETSc libraries as well as MPI
112:        - provides summary and diagnostic information if certain runtime
113:          options are chosen (e.g., -log_view).
114:   */
115:   PetscCall(PetscFinalize());
116:   return 0;
117: }

119: /*TEST

121:    testset:
122:       output_file: output/ex102_1.out
123:       nsize: {{1 2}}
124:       args: -low_rank {{0 1}} -use_c {{0 1}}
125:       test:
126:         suffix: standard
127:       test:
128:         suffix: cuda
129:         requires: cuda
130:         args: -A_mat_type aijcusparse -U_mat_type densecuda -V_mat_type densecuda

132: TEST*/