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 ierr = MatSetOptionsPrefix(B,"AB_");CHKERRQ(ierr); 43 } else { 44 /* Compute B = A^t*A; B misses 0-th diagonal */ 45 ierr = MatProductSetType(B,MATPRODUCT_AtB);CHKERRQ(ierr); 46 ierr = MatSetOptionsPrefix(B,"AtB_");CHKERRQ(ierr); 47 } 48 49 /* Force allocate missing diagonal entries of B */ 50 ierr = MatSetOption(B,MAT_FORCE_DIAGONAL_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 51 ierr = MatProductSetFromOptions(B);CHKERRQ(ierr); 52 53 ierr = MatProductSymbolic(B);CHKERRQ(ierr); 54 ierr = MatProductNumeric(B);CHKERRQ(ierr); 55 56 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 57 58 /* Insert entries to diagonal of B */ 59 ierr = MatCreateVecs(B,NULL,&diag);CHKERRQ(ierr); 60 ierr = MatGetDiagonal(B,diag);CHKERRQ(ierr); 61 ierr = VecSetValue(diag,0,100.0,INSERT_VALUES);CHKERRQ(ierr); 62 ierr = VecAssemblyBegin(diag);CHKERRQ(ierr); 63 ierr = VecAssemblyEnd(diag);CHKERRQ(ierr); 64 65 ierr = MatDiagonalSet(B,diag,INSERT_VALUES);CHKERRQ(ierr); 66 if (test == 1) { 67 ierr = MatView(B,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 68 } 69 ierr = MatDestroy(&B);CHKERRQ(ierr); 70 ierr = VecDestroy(&diag);CHKERRQ(ierr); 71 } 72 73 ierr = MatDestroy(&A);CHKERRQ(ierr); 74 ierr = PetscFinalize(); 75 return ierr; 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