Actual source code: ex102.c
2: static char help[] = "Tests MatCreateLRC()\n\n";
4: #include <petscmat.h>
6: int main(int argc, char **args)
7: {
8: Vec x, b, c = NULL;
9: Mat A, U, V, LR, X, LRe;
10: PetscInt M = 5, N = 7;
11: PetscBool flg;
13: PetscFunctionBeginUser;
14: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
15: PetscCall(PetscOptionsGetInt(NULL, NULL, "-m", &M, NULL));
16: PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &N, NULL));
18: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
19: Create the sparse matrix
20: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
21: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
22: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, M, N));
23: PetscCall(MatSetOptionsPrefix(A, "A_"));
24: PetscCall(MatSetFromOptions(A));
25: PetscCall(MatSeqAIJSetPreallocation(A, 5, NULL));
26: PetscCall(MatMPIAIJSetPreallocation(A, 5, NULL, 5, NULL));
27: PetscCall(MatSetUp(A));
28: PetscCall(MatSetRandom(A, NULL));
30: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
31: Create the dense matrices
32: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
33: PetscCall(MatCreate(PETSC_COMM_WORLD, &U));
34: PetscCall(MatSetSizes(U, PETSC_DECIDE, PETSC_DECIDE, M, 3));
35: PetscCall(MatSetType(U, MATDENSE));
36: PetscCall(MatSetOptionsPrefix(U, "U_"));
37: PetscCall(MatSetFromOptions(U));
38: PetscCall(MatSetUp(U));
39: PetscCall(MatSetRandom(U, NULL));
41: PetscCall(MatCreate(PETSC_COMM_WORLD, &V));
42: PetscCall(MatSetSizes(V, PETSC_DECIDE, PETSC_DECIDE, N, 3));
43: PetscCall(MatSetType(V, MATDENSE));
44: PetscCall(MatSetOptionsPrefix(V, "V_"));
45: PetscCall(MatSetFromOptions(V));
46: PetscCall(MatSetUp(V));
47: PetscCall(MatSetRandom(V, NULL));
49: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
50: Create a vector to hold the diagonal of C
51: A sequential vector can be created as well on each process
52: It is user responsibility to ensure the data in the vector
53: is consistent across processors
54: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55: PetscCall(PetscOptionsHasName(NULL, NULL, "-use_c", &flg));
56: if (flg) {
57: PetscCall(VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, 3, &c));
58: PetscCall(VecSetRandom(c, NULL));
59: }
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Create low rank correction matrix
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: PetscCall(PetscOptionsHasName(NULL, NULL, "-low_rank", &flg));
65: if (flg) {
66: /* create a low-rank matrix, with no A-matrix */
67: PetscCall(MatCreateLRC(NULL, U, c, V, &LR));
68: PetscCall(MatDestroy(&A));
69: } else {
70: PetscCall(MatCreateLRC(A, U, c, V, &LR));
71: }
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Create the low rank correction matrix explicitly to check for
75: correctness
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: PetscCall(MatHermitianTranspose(V, MAT_INITIAL_MATRIX, &X));
78: PetscCall(MatDiagonalScale(X, c, NULL));
79: PetscCall(MatMatMult(U, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &LRe));
80: PetscCall(MatDestroy(&X));
81: if (A) PetscCall(MatAYPX(LRe, 1.0, A, DIFFERENT_NONZERO_PATTERN));
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create test vectors
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: PetscCall(MatCreateVecs(LR, &x, &b));
87: PetscCall(VecSetRandom(x, NULL));
88: PetscCall(MatMult(LR, x, b));
89: PetscCall(MatMultTranspose(LR, b, x));
90: PetscCall(VecDestroy(&x));
91: PetscCall(VecDestroy(&b));
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Check correctness
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: PetscCall(MatMultEqual(LR, LRe, 10, &flg));
97: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error in MatMult\n"));
98: #if !defined(PETSC_USE_COMPLEX)
99: PetscCall(MatMultHermitianTransposeEqual(LR, LRe, 10, &flg));
100: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error in MatMultTranspose\n"));
101: #endif
103: PetscCall(MatDestroy(&A));
104: PetscCall(MatDestroy(&LRe));
105: PetscCall(MatDestroy(&U));
106: PetscCall(MatDestroy(&V));
107: PetscCall(VecDestroy(&c));
108: PetscCall(MatDestroy(&LR));
110: /*
111: Always call PetscFinalize() before exiting a program. This routine
112: - finalizes the PETSc libraries as well as MPI
113: - provides summary and diagnostic information if certain runtime
114: options are chosen (e.g., -log_view).
115: */
116: PetscCall(PetscFinalize());
117: return 0;
118: }
120: /*TEST
122: testset:
123: output_file: output/ex102_1.out
124: nsize: {{1 2}}
125: args: -low_rank {{0 1}} -use_c {{0 1}}
126: test:
127: suffix: standard
128: test:
129: suffix: cuda
130: requires: cuda
131: args: -A_mat_type aijcusparse -U_mat_type densecuda -V_mat_type densecuda
133: TEST*/