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