Actual source code: ex38.c
petsc-3.14.6 2021-03-30
2: static char help[] = "Test interface of Elemental. \n\n";
4: #include <petscmat.h>
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;
13: PetscBool flg,Test_MatMatMult=PETSC_FALSE,mats_view=PETSC_FALSE;
14: PetscScalar *v;
15: PetscMPIInt rank,size;
17: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
18: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
19: MPI_Comm_size(PETSC_COMM_WORLD,&size);
20: PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);
22: /* Get local block or element size*/
23: PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
24: n = m;
25: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
27: MatCreate(PETSC_COMM_WORLD,&C);
28: MatSetSizes(C,m,n,PETSC_DECIDE,PETSC_DECIDE);
29: MatSetType(C,MATELEMENTAL);
30: MatSetFromOptions(C);
31: MatSetUp(C);
33: PetscOptionsHasName(NULL,NULL,"-row_oriented",&flg);
34: if (flg) {MatSetOption(C,MAT_ROW_ORIENTED,PETSC_TRUE);}
35: MatGetOwnershipIS(C,&isrows,&iscols);
36: PetscOptionsHasName(NULL,NULL,"-Cexp_view_ownership",&flg);
37: if (flg) { /* View ownership of explicit C */
38: IS tmp;
39: PetscPrintf(PETSC_COMM_WORLD,"Ownership of explicit C:\n");
40: PetscPrintf(PETSC_COMM_WORLD,"Row index set:\n");
41: ISOnComm(isrows,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp);
42: ISView(tmp,PETSC_VIEWER_STDOUT_WORLD);
43: ISDestroy(&tmp);
44: PetscPrintf(PETSC_COMM_WORLD,"Column index set:\n");
45: ISOnComm(iscols,PETSC_COMM_WORLD,PETSC_USE_POINTER,&tmp);
46: ISView(tmp,PETSC_VIEWER_STDOUT_WORLD);
47: ISDestroy(&tmp);
48: }
50: /* Set local matrix entries */
51: ISGetLocalSize(isrows,&nrows);
52: ISGetIndices(isrows,&rows);
53: ISGetLocalSize(iscols,&ncols);
54: ISGetIndices(iscols,&cols);
55: 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: MatSetValues(C,nrows,rows,ncols,cols,v,INSERT_VALUES);
63: ISRestoreIndices(isrows,&rows);
64: ISRestoreIndices(iscols,&cols);
65: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
68: /* Test MatView() */
69: if (mats_view) {
70: MatView(C,PETSC_VIEWER_STDOUT_WORLD);
71: }
73: /* Test MatMissingDiagonal() */
74: MatMissingDiagonal(C,&flg,NULL);
75: if (flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"MatMissingDiagonal() did not return false!\n");
77: /* Set unowned matrix entries - add subdiagonals and diagonals from proc[0] */
78: if (!rank) {
79: PetscInt M,N,cols[2];
80: 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: MatSetValues(C,1,&i,2,cols,v,ADD_VALUES);
86: } else {
87: MatSetValues(C,1,&i,1,&i,v,ADD_VALUES);
88: }
89: }
90: }
91: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
92: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
94: /* Test MatMult() */
95: MatComputeOperator(C,MATAIJ,&Caij);
96: MatMultEqual(C,Caij,5,&flg);
97: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultEqual() fails");
98: MatMultTransposeEqual(C,Caij,5,&flg);
99: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeEqual() fails");
101: /* Test MatMultAdd() and MatMultTransposeAddEqual() */
102: MatMultAddEqual(C,Caij,5,&flg);
103: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultAddEqual() fails");
104: MatMultTransposeAddEqual(C,Caij,5,&flg);
105: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"C != Caij. MatMultTransposeAddEqual() fails");
107: /* Test MatMatMult() */
108: PetscOptionsHasName(NULL,NULL,"-test_matmatmult",&Test_MatMatMult);
109: if (Test_MatMatMult) {
110: Mat CCelem,CCaij;
111: MatMatMult(C,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCelem);
112: MatMatMult(Caij,Caij,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&CCaij);
113: MatMultEqual(CCelem,CCaij,5,&flg);
114: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"CCelem != CCaij. MatMatMult() fails");
115: MatDestroy(&CCaij);
116: MatDestroy(&CCelem);
117: }
119: PetscFree(v);
120: ISDestroy(&isrows);
121: ISDestroy(&iscols);
122: MatDestroy(&C);
123: MatDestroy(&Caij);
124: PetscFinalize();
125: return ierr;
126: }
131: /*TEST
133: test:
134: nsize: 2
135: requires: elemental
136: args: -mat_type elemental -m 2 -n 3
138: test:
139: suffix: 2
140: nsize: 6
141: requires: elemental
142: args: -mat_type elemental -m 2 -n 2
144: test:
145: suffix: 3
146: nsize: 6
147: requires: elemental
148: args: -mat_type elemental -m 2 -n 2 -test_matmatmult
149: output_file: output/ex38_2.out
151: TEST*/