xref: /petsc/src/mat/tests/ex161.c (revision ccb4e88a40f0b86eaeca07ff64c64e4de2fae686)
1 static char help[] = "Test MatTransposeColoring for SeqAIJ matrices. Used for '-matmattransmult_color' on  MatMatTransposeMult \n\n";
2 
3 #include <petscmat.h>
4 #include <petsc/private/matimpl.h> /* Need struct _p_MatTransposeColoring for this test. */
5 
6 int main(int argc,char **argv)
7 {
8   Mat                  A,R,C,C_dense,C_sparse,Rt_dense,P,PtAP;
9   PetscInt             row,col,m,n;
10   PetscErrorCode       ierr;
11   MatScalar            one         =1.0,val;
12   MatColoring          mc;
13   MatTransposeColoring matcoloring = 0;
14   ISColoring           iscoloring;
15   PetscBool            equal;
16   PetscMPIInt          size;
17 
18   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
19   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
20   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!");
21 
22   /* Create seqaij A */
23   ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr);
24   ierr = MatSetSizes(A,4,4,4,4);CHKERRQ(ierr);
25   ierr = MatSetType(A,MATSEQAIJ);CHKERRQ(ierr);
26   ierr = MatSetFromOptions(A);CHKERRQ(ierr);
27   ierr = MatSetUp(A);CHKERRQ(ierr);
28   row  = 0; col=0; val=1.0; ierr = MatSetValues(A,1,&row,1,&col,&val,ADD_VALUES);CHKERRQ(ierr);
29   row  = 1; col=3; val=2.0; ierr = MatSetValues(A,1,&row,1,&col,&val,ADD_VALUES);CHKERRQ(ierr);
30   row  = 2; col=2; val=3.0; ierr = MatSetValues(A,1,&row,1,&col,&val,ADD_VALUES);CHKERRQ(ierr);
31   row  = 3; col=0; val=4.0; ierr = MatSetValues(A,1,&row,1,&col,&val,ADD_VALUES);CHKERRQ(ierr);
32   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
33   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
34   ierr = MatSetOptionsPrefix(A,"A_");CHKERRQ(ierr);
35   ierr = MatView(A,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
36   ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr);
37 
38   /* Create seqaij R */
39   ierr = MatCreate(PETSC_COMM_SELF,&R);CHKERRQ(ierr);
40   ierr = MatSetSizes(R,2,4,2,4);CHKERRQ(ierr);
41   ierr = MatSetType(R,MATSEQAIJ);CHKERRQ(ierr);
42   ierr = MatSetFromOptions(R);CHKERRQ(ierr);
43   ierr = MatSetUp(R);CHKERRQ(ierr);
44   row  = 0; col=0; ierr = MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES);CHKERRQ(ierr);
45   row  = 0; col=1; ierr = MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES);CHKERRQ(ierr);
46 
47   row  = 1; col=1; ierr = MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES);CHKERRQ(ierr);
48   row  = 1; col=2; ierr = MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES);CHKERRQ(ierr);
49   row  = 1; col=3; ierr = MatSetValues(R,1,&row,1,&col,&one,ADD_VALUES);CHKERRQ(ierr);
50   ierr = MatAssemblyBegin(R,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
51   ierr = MatAssemblyEnd(R,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
52   ierr = MatSetOptionsPrefix(R,"R_");CHKERRQ(ierr);
53   ierr = MatView(R,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
54   ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr);
55 
56   /* C = A*R^T */
57   ierr = MatMatTransposeMult(A,R,MAT_INITIAL_MATRIX,2.0,&C);CHKERRQ(ierr);
58   ierr = MatSetOptionsPrefix(C,"ARt_");CHKERRQ(ierr);
59   ierr = MatView(C,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
60   ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr);
61 
62   /* Create MatTransposeColoring from symbolic C=A*R^T */
63   ierr = MatColoringCreate(C,&mc);CHKERRQ(ierr);
64   ierr = MatColoringSetDistance(mc,2);CHKERRQ(ierr);
65   /* ierr = MatColoringSetType(mc,MATCOLORINGSL);CHKERRQ(ierr); */
66   ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr);
67   ierr = MatColoringApply(mc,&iscoloring);CHKERRQ(ierr);
68   ierr = MatColoringDestroy(&mc);CHKERRQ(ierr);
69   ierr = MatTransposeColoringCreate(C,iscoloring,&matcoloring);CHKERRQ(ierr);
70   ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
71 
72   /* Create Rt_dense */
73   ierr = MatCreate(PETSC_COMM_WORLD,&Rt_dense);CHKERRQ(ierr);
74   ierr = MatSetSizes(Rt_dense,4,matcoloring->ncolors,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
75   ierr = MatSetType(Rt_dense,MATDENSE);CHKERRQ(ierr);
76   ierr = MatSeqDenseSetPreallocation(Rt_dense,NULL);CHKERRQ(ierr);
77   ierr = MatAssemblyBegin(Rt_dense,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
78   ierr = MatAssemblyEnd(Rt_dense,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
79   ierr = MatGetLocalSize(Rt_dense,&m,&n);CHKERRQ(ierr);
80   ierr = PetscPrintf(PETSC_COMM_SELF,"Rt_dense: %D,%D\n",m,n);CHKERRQ(ierr);
81 
82   /* Get Rt_dense by Apply MatTransposeColoring to R */
83   ierr = MatTransColoringApplySpToDen(matcoloring,R,Rt_dense);CHKERRQ(ierr);
84 
85   /* C_dense = A*Rt_dense */
86   ierr = MatMatMult(A,Rt_dense,MAT_INITIAL_MATRIX,2.0,&C_dense);CHKERRQ(ierr);
87   ierr = MatSetOptionsPrefix(C_dense,"ARt_dense_");CHKERRQ(ierr);
88   /*ierr = MatView(C_dense,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
89   /*ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr); */
90 
91   /* Recover C from C_dense */
92   ierr = MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&C_sparse);CHKERRQ(ierr);
93   ierr = MatTransColoringApplyDenToSp(matcoloring,C_dense,C_sparse);CHKERRQ(ierr);
94   ierr = MatSetOptionsPrefix(C_sparse,"ARt_color_");CHKERRQ(ierr);
95   ierr = MatView(C_sparse,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
96   ierr = PetscPrintf(PETSC_COMM_SELF,"\n");CHKERRQ(ierr);
97 
98   ierr = MatDestroy(&C_dense);CHKERRQ(ierr);
99   ierr = MatDestroy(&C_sparse);CHKERRQ(ierr);
100   ierr = MatDestroy(&Rt_dense);CHKERRQ(ierr);
101   ierr = MatTransposeColoringDestroy(&matcoloring);CHKERRQ(ierr);
102   ierr = MatDestroy(&C);CHKERRQ(ierr);
103 
104   /* Test PtAP = P^T*A*P, P = R^T */
105   ierr = MatTranspose(R,MAT_INITIAL_MATRIX,&P);CHKERRQ(ierr);
106   ierr = MatPtAP(A,P,MAT_INITIAL_MATRIX,2.0,&PtAP);CHKERRQ(ierr);
107   ierr = MatSetOptionsPrefix(PtAP,"PtAP_");CHKERRQ(ierr);
108   ierr = MatView(PtAP,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
109   ierr = MatDestroy(&P);CHKERRQ(ierr);
110 
111   /* Test C = RARt */
112   ierr = MatRARt(A,R,MAT_INITIAL_MATRIX,2.0,&C);CHKERRQ(ierr);
113   ierr = MatRARt(A,R,MAT_REUSE_MATRIX,2.0,&C);CHKERRQ(ierr);
114   ierr = MatEqual(C,PtAP,&equal);CHKERRQ(ierr);
115   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: PtAP != RARt");
116 
117   /* Free spaces */
118   ierr = MatDestroy(&C);CHKERRQ(ierr);
119   ierr = MatDestroy(&A);CHKERRQ(ierr);
120   ierr = MatDestroy(&R);CHKERRQ(ierr);
121   ierr = MatDestroy(&PtAP);CHKERRQ(ierr);
122   ierr = PetscFinalize();
123   return ierr;
124 }
125 
126 /*TEST
127 
128    test:
129       output_file: output/ex161.out
130    test:
131       suffix: 2
132       args: -matmattransmult_via color
133       output_file: output/ex161.out
134 
135    test:
136       suffix: 3
137       args: -matmattransmult_via color -matden2sp_brows 3
138       output_file: output/ex161.out
139 
140    test:
141       suffix: 4
142       args: -matmattransmult_via color -matrart_via r*art
143       output_file: output/ex161.out
144 
145    test:
146       suffix: 5
147       args: -matmattransmult_via color -matrart_via coloring_rart
148       output_file: output/ex161.out
149 
150 TEST*/
151