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*/