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