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 PetscCall(PetscInitialize(&argc,&args,(char*)0,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.; v[1] = 2.; v[2] = -1.; 24 for (i=2; i<n-1; i++) { 25 col[0] = i-1; col[1] = i; col[2] = i+1; 26 PetscCall(MatSetValues(A,1,&i,3,col,v,INSERT_VALUES)); 27 } 28 i = 1; col[0] = 1; col[1] = 2; 29 PetscCall(MatSetValues(A,1,&i,2,col,v+1,INSERT_VALUES)); 30 i = n - 1; col[0] = n - 2; col[1] = n - 1; 31 PetscCall(MatSetValues(A,1,&i,2,col,v,INSERT_VALUES)); 32 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 33 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 34 35 for (test = 0; test < 2; test++) { 36 PetscCall(MatProductCreate(A,A,NULL,&B)); 37 38 if (test == 0) { 39 /* Compute B = A*A; B misses 0-th diagonal */ 40 PetscCall(MatProductSetType(B,MATPRODUCT_AB)); 41 PetscCall(MatSetOptionsPrefix(B,"AB_")); 42 } else { 43 /* Compute B = A^t*A; B misses 0-th diagonal */ 44 PetscCall(MatProductSetType(B,MATPRODUCT_AtB)); 45 PetscCall(MatSetOptionsPrefix(B,"AtB_")); 46 } 47 48 /* Force allocate missing diagonal entries of B */ 49 PetscCall(MatSetOption(B,MAT_FORCE_DIAGONAL_ENTRIES,PETSC_TRUE)); 50 PetscCall(MatProductSetFromOptions(B)); 51 52 PetscCall(MatProductSymbolic(B)); 53 PetscCall(MatProductNumeric(B)); 54 55 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 56 57 /* Insert entries to diagonal of B */ 58 PetscCall(MatCreateVecs(B,NULL,&diag)); 59 PetscCall(MatGetDiagonal(B,diag)); 60 PetscCall(VecSetValue(diag,0,100.0,INSERT_VALUES)); 61 PetscCall(VecAssemblyBegin(diag)); 62 PetscCall(VecAssemblyEnd(diag)); 63 64 PetscCall(MatDiagonalSet(B,diag,INSERT_VALUES)); 65 if (test == 1) { 66 PetscCall(MatView(B,PETSC_VIEWER_STDOUT_WORLD)); 67 } 68 PetscCall(MatDestroy(&B)); 69 PetscCall(VecDestroy(&diag)); 70 } 71 72 PetscCall(MatDestroy(&A)); 73 PetscCall(PetscFinalize()); 74 return 0; 75 } 76 77 /*TEST 78 79 test: 80 output_file: output/ex81_1.out 81 82 test: 83 suffix: 2 84 args: -AtB_mat_product_algorithm at*b 85 output_file: output/ex81_1.out 86 87 test: 88 suffix: 3 89 args: -AtB_mat_product_algorithm outerproduct 90 output_file: output/ex81_1.out 91 92 test: 93 suffix: 4 94 nsize: 3 95 args: -AtB_mat_product_algorithm nonscalable 96 output_file: output/ex81_3.out 97 98 test: 99 suffix: 5 100 nsize: 3 101 args: -AtB_mat_product_algorithm scalable 102 output_file: output/ex81_3.out 103 104 test: 105 suffix: 6 106 nsize: 3 107 args: -AtB_mat_product_algorithm at*b 108 output_file: output/ex81_3.out 109 110 TEST*/ 111