xref: /petsc/src/mat/tests/ex38.c (revision 28b400f66ebc7ae0049166a2294dfcd3df27e64b)
1 
2 static char help[] = "Test interface of Elemental. \n\n";
3 
4 #include <petscmat.h>
5 
6 int main(int argc,char **args)
7 {
8   Mat            C,Caij;
9   PetscInt       i,j,m = 5,n,nrows,ncols;
10   const PetscInt *rows,*cols;
11   IS             isrows,iscols;
12   PetscErrorCode ierr;
13   PetscBool      flg,Test_MatMatMult=PETSC_FALSE,mats_view=PETSC_FALSE;
14   PetscScalar    *v;
15   PetscMPIInt    rank,size;
16 
17   ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
18   CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
19   CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
20   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view));
21 
22   /* Get local block or element size*/
23   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL));
24   n    = m;
25   CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
26 
27   CHKERRQ(MatCreate(PETSC_COMM_WORLD,&C));
28   CHKERRQ(MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE));
29   CHKERRQ(MatSetType(C,MATELEMENTAL));
30   CHKERRQ(MatSetFromOptions(C));
31   CHKERRQ(MatSetUp(C));
32 
33   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-row_oriented",&flg));
34   if (flg) CHKERRQ(MatSetOption(C,MAT_ROW_ORIENTED,PETSC_TRUE));
35   CHKERRQ(MatGetOwnershipIS(C,&isrows,&iscols));
36   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-Cexp_view_ownership",&flg));
37   if (flg) { /* View ownership of explicit C */
38     IS tmp;
39     CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Ownership of explicit C:\n"));
40     CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Row index set:\n"));
41     CHKERRQ(ISOnComm(isrows,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp));
42     CHKERRQ(ISView(tmp,PETSC_VIEWER_STDOUT_WORLD));
43     CHKERRQ(ISDestroy(&tmp));
44     CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Column index set:\n"));
45     CHKERRQ(ISOnComm(iscols,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp));
46     CHKERRQ(ISView(tmp,PETSC_VIEWER_STDOUT_WORLD));
47     CHKERRQ(ISDestroy(&tmp));
48   }
49 
50   /* Set local matrix entries */
51   CHKERRQ(ISGetLocalSize(isrows,&nrows));
52   CHKERRQ(ISGetIndices(isrows,&rows));
53   CHKERRQ(ISGetLocalSize(iscols,&ncols));
54   CHKERRQ(ISGetIndices(iscols,&cols));
55   CHKERRQ(PetscMalloc1(nrows*ncols,&v));
56   for (i=0; i<nrows; i++) {
57     for (j=0; j<ncols; j++) {
58       /*v[i*ncols+j] = (PetscReal)(rank);*/
59       v[i*ncols+j] = (PetscReal)(rank*10000+100*rows[i]+cols[j]);
60     }
61   }
62   CHKERRQ(MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES));
63   CHKERRQ(ISRestoreIndices(isrows,&rows));
64   CHKERRQ(ISRestoreIndices(iscols,&cols));
65   CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
66   CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
67 
68   /* Test MatView() */
69   if (mats_view) {
70     CHKERRQ(MatView(C,PETSC_VIEWER_STDOUT_WORLD));
71   }
72 
73   /* Test MatMissingDiagonal() */
74   CHKERRQ(MatMissingDiagonal(C,&flg,NULL));
75   PetscCheck(!flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"MatMissingDiagonal() did not return false!");
76 
77   /* Set unowned matrix entries - add subdiagonals and diagonals from proc[0] */
78   if (rank == 0) {
79     PetscInt M,N,cols[2];
80     CHKERRQ(MatGetSize(C,&M,&N));
81     for (i=0; i<M; i++) {
82       cols[0] = i;   v[0] = i + 0.5;
83       cols[1] = i-1; v[1] = 0.5;
84       if (i) {
85         CHKERRQ(MatSetValues(C,1,&i,2,cols,v,ADD_VALUES));
86       } else {
87         CHKERRQ(MatSetValues(C,1,&i,1,&i,v,ADD_VALUES));
88       }
89     }
90   }
91   CHKERRQ(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
92   CHKERRQ(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
93 
94   /* Test MatMult() */
95   CHKERRQ(MatComputeOperator(C,MATAIJ,&Caij));
96   CHKERRQ(MatMultEqual(C,Caij,5,&flg));
97   PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultEqual() fails");
98   CHKERRQ(MatMultTransposeEqual(C,Caij,5,&flg));
99   PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeEqual() fails");
100 
101   /* Test MatMultAdd() and MatMultTransposeAddEqual() */
102   CHKERRQ(MatMultAddEqual(C,Caij,5,&flg));
103   PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultAddEqual() fails");
104   CHKERRQ(MatMultTransposeAddEqual(C,Caij,5,&flg));
105   PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeAddEqual() fails");
106 
107   /* Test MatMatMult() */
108   CHKERRQ(PetscOptionsHasName(NULL,NULL,"-test_matmatmult",&Test_MatMatMult));
109   if (Test_MatMatMult) {
110     Mat CCelem,CCaij;
111     CHKERRQ(MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCelem));
112     CHKERRQ(MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij));
113     CHKERRQ(MatMultEqual(CCelem,CCaij,5,&flg));
114     PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"CCelem != CCaij. MatMatMult() fails");
115     CHKERRQ(MatDestroy(&CCaij));
116     CHKERRQ(MatDestroy(&CCelem));
117   }
118 
119   CHKERRQ(PetscFree(v));
120   CHKERRQ(ISDestroy(&isrows));
121   CHKERRQ(ISDestroy(&iscols));
122   CHKERRQ(MatDestroy(&C));
123   CHKERRQ(MatDestroy(&Caij));
124   ierr = PetscFinalize();
125   return ierr;
126 }
127 
128 /*TEST
129 
130    test:
131       nsize: 2
132       requires: elemental
133       args: -mat_type elemental -m 2 -n 3
134 
135    test:
136       suffix: 2
137       nsize: 6
138       requires: elemental
139       args: -mat_type elemental -m 2 -n 2
140 
141    test:
142       suffix: 3
143       nsize: 6
144       requires: elemental
145       args: -mat_type elemental -m 2 -n 2 -test_matmatmult
146       output_file: output/ex38_2.out
147 
148 TEST*/
149