1 static char help[] = "Tests MatCreateLRC()\n\n";
2
3 #include <petscmat.h>
4
main(int argc,char ** args)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;
11
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));
16
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));
28
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));
39
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));
47
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 }
59
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 }
71
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));
81
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));
91
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
101
102 PetscCall(MatDestroy(&A));
103 PetscCall(MatDestroy(&LRe));
104 PetscCall(MatDestroy(&U));
105 PetscCall(MatDestroy(&V));
106 PetscCall(VecDestroy(&c));
107 PetscCall(MatDestroy(&LR));
108
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 }
118
119 /*TEST
120
121 testset:
122 output_file: output/empty.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
131
132 TEST*/
133