1 2 static char help[] = "Tests MatCreateLRC()\n\n"; 3 4 #include <petscmat.h> 5 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; 12 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)); 17 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)); 29 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)); 40 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)); 48 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 } 60 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 } 72 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)); 82 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)); 92 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 102 103 PetscCall(MatDestroy(&A)); 104 PetscCall(MatDestroy(&LRe)); 105 PetscCall(MatDestroy(&U)); 106 PetscCall(MatDestroy(&V)); 107 PetscCall(VecDestroy(&c)); 108 PetscCall(MatDestroy(&LR)); 109 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 } 119 120 /*TEST 121 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 132 133 TEST*/ 134