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