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