1 static char help[] = "Tests MatOption MAT_FORCE_DIAGONAL_ENTRIES.\n\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **args) 6 { 7 Mat A, B; 8 Vec diag; 9 PetscInt i, n = 10, col[3], test; 10 PetscScalar v[3]; 11 12 PetscFunctionBeginUser; 13 PetscCall(PetscInitialize(&argc, &args, NULL, help)); 14 PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); 15 16 /* Create A which has empty 0-th row and column */ 17 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 18 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n)); 19 PetscCall(MatSetType(A, MATAIJ)); 20 PetscCall(MatSetFromOptions(A)); 21 PetscCall(MatSetUp(A)); 22 23 v[0] = -1.; 24 v[1] = 2.; 25 v[2] = -1.; 26 for (i = 2; i < n - 1; i++) { 27 col[0] = i - 1; 28 col[1] = i; 29 col[2] = i + 1; 30 PetscCall(MatSetValues(A, 1, &i, 3, col, v, INSERT_VALUES)); 31 } 32 i = 1; 33 col[0] = 1; 34 col[1] = 2; 35 PetscCall(MatSetValues(A, 1, &i, 2, col, v + 1, INSERT_VALUES)); 36 i = n - 1; 37 col[0] = n - 2; 38 col[1] = n - 1; 39 PetscCall(MatSetValues(A, 1, &i, 2, col, v, INSERT_VALUES)); 40 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 41 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 42 43 for (test = 0; test < 2; test++) { 44 PetscCall(MatProductCreate(A, A, NULL, &B)); 45 46 if (test == 0) { 47 /* Compute B = A*A; B misses 0-th diagonal */ 48 PetscCall(MatProductSetType(B, MATPRODUCT_AB)); 49 PetscCall(MatSetOptionsPrefix(B, "AB_")); 50 } else { 51 /* Compute B = A^t*A; B misses 0-th diagonal */ 52 PetscCall(MatProductSetType(B, MATPRODUCT_AtB)); 53 PetscCall(MatSetOptionsPrefix(B, "AtB_")); 54 } 55 56 /* Force allocate missing diagonal entries of B */ 57 PetscCall(MatSetOption(B, MAT_FORCE_DIAGONAL_ENTRIES, PETSC_TRUE)); 58 PetscCall(MatProductSetFromOptions(B)); 59 60 PetscCall(MatProductSymbolic(B)); 61 PetscCall(MatProductNumeric(B)); 62 63 PetscCall(MatSetOption(B, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 64 65 /* Insert entries to diagonal of B */ 66 PetscCall(MatCreateVecs(B, NULL, &diag)); 67 PetscCall(MatGetDiagonal(B, diag)); 68 PetscCall(VecSetValue(diag, 0, 100.0, INSERT_VALUES)); 69 PetscCall(VecAssemblyBegin(diag)); 70 PetscCall(VecAssemblyEnd(diag)); 71 72 PetscCall(MatDiagonalSet(B, diag, INSERT_VALUES)); 73 if (test == 1) PetscCall(MatView(B, PETSC_VIEWER_STDOUT_WORLD)); 74 PetscCall(MatDestroy(&B)); 75 PetscCall(VecDestroy(&diag)); 76 } 77 78 PetscCall(MatDestroy(&A)); 79 PetscCall(PetscFinalize()); 80 return 0; 81 } 82 83 /*TEST 84 85 test: 86 output_file: output/ex81_1.out 87 88 test: 89 suffix: 2 90 args: -AtB_mat_product_algorithm at*b 91 output_file: output/ex81_1.out 92 93 test: 94 suffix: 3 95 args: -AtB_mat_product_algorithm outerproduct 96 output_file: output/ex81_1.out 97 98 test: 99 suffix: 4 100 nsize: 3 101 args: -AtB_mat_product_algorithm nonscalable 102 output_file: output/ex81_3.out 103 104 test: 105 suffix: 5 106 nsize: 3 107 args: -AtB_mat_product_algorithm scalable 108 output_file: output/ex81_3.out 109 110 test: 111 suffix: 6 112 nsize: 3 113 args: -AtB_mat_product_algorithm at*b 114 output_file: output/ex81_3.out 115 116 TEST*/ 117