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 = MatProductNumeric(B);CHKERRQ(ierr); 53 54 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 55 56 /* Insert entries to diagonal of B */ 57 ierr = MatCreateVecs(B,NULL,&diag);CHKERRQ(ierr); 58 ierr = MatGetDiagonal(B,diag);CHKERRQ(ierr); 59 ierr = VecSetValue(diag,0,100.0,INSERT_VALUES);CHKERRQ(ierr); 60 ierr = VecAssemblyBegin(diag);CHKERRQ(ierr); 61 ierr = VecAssemblyEnd(diag);CHKERRQ(ierr); 62 63 ierr = MatDiagonalSet(B,diag,INSERT_VALUES);CHKERRQ(ierr); 64 if (test == 1) { 65 ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 66 } 67 ierr = MatDestroy(&B);CHKERRQ(ierr); 68 ierr = VecDestroy(&diag);CHKERRQ(ierr); 69 } 70 71 ierr = MatDestroy(&A);CHKERRQ(ierr); 72 ierr = PetscFinalize(); 73 return ierr; 74 } 75 76 /*TEST 77 78 test: 79 output_file: output/ex81_1.out 80 81 test: 82 suffix: 2 83 args: -matproduct_atb_via at*b 84 output_file: output/ex81_1.out 85 86 test: 87 suffix: 3 88 args: -matproduct_atb_via outerproduct 89 output_file: output/ex81_1.out 90 91 test: 92 suffix: 4 93 nsize: 3 94 args: -matproduct_atb_via nonscalable 95 output_file: output/ex81_3.out 96 97 test: 98 suffix: 5 99 nsize: 3 100 args: -matproduct_atb_via scalable 101 output_file: output/ex81_3.out 102 103 test: 104 suffix: 6 105 nsize: 3 106 args: -matproduct_atb_via at*b 107 output_file: output/ex81_3.out 108 109 TEST*/ 110