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