1 static char help[] = "Tests MatCreateLRC()\n\n"; 2 3 #include <petscmat.h> 4 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, (char *)0, 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_DEFAULT, &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/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 131 132 TEST*/ 133