Actual source code: ex240.c
1: static char help[] ="Tests MatFDColoringSetValues()\n\n";
3: #include <petscdm.h>
4: #include <petscdmda.h>
6: int main(int argc,char **argv)
7: {
8: DM da;
9: PetscInt N, mx = 5,my = 4,i,j,nc,nrow,n,ncols,rstart,*colors,*map;
10: const PetscInt *cols;
11: const PetscScalar *vals;
12: Mat A,B;
13: PetscReal norm;
14: MatFDColoring fdcoloring;
15: ISColoring iscoloring;
16: PetscScalar *cm;
17: const ISColoringValue *icolors;
18: PetscMPIInt rank;
19: ISLocalToGlobalMapping ltog;
20: PetscBool single,two;
22: PetscInitialize(&argc,&argv,NULL,help);
23: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
24: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,mx,my,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
25: DMSetUp(da);
26: DMCreateMatrix(da,&A);
28: /* as a test copy the matrices from the standard format to the compressed format; this is not scalable but is ok because just for testing */
29: /* first put the coloring in the global ordering */
30: DMCreateColoring(da,IS_COLORING_LOCAL,&iscoloring);
31: ISColoringGetColors(iscoloring,&n,&nc,&icolors);
32: DMGetLocalToGlobalMapping(da,<og);
33: PetscMalloc1(n,&map);
34: for (i=0; i<n; i++) map[i] = i;
35: ISLocalToGlobalMappingApply(ltog,n,map,map);
36: MatGetSize(A,&N,NULL);
37: PetscMalloc1(N,&colors);
38: for (i=0; i<N; i++) colors[i] = -1;
39: for (i=0; i<n; i++) colors[map[i]]= icolors[i];
40: PetscFree(map);
41: PetscSynchronizedPrintf(MPI_COMM_WORLD,"[%d]Global colors \n",rank);
42: for (i=0; i<N; i++) PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT " %" PetscInt_FMT "\n",i,colors[i]);
43: PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
45: /* second, compress the A matrix */
46: MatSetRandom(A,NULL);
47: MatView(A,NULL);
48: MatGetLocalSize(A,&nrow,NULL);
49: MatGetOwnershipRange(A,&rstart,NULL);
50: PetscCalloc1(nrow*nc,&cm);
51: for (i=0; i<nrow; i++) {
52: MatGetRow(A,rstart+i,&ncols,&cols,&vals);
53: for (j=0; j<ncols; j++) {
55: cm[i + nrow*colors[cols[j]]] = vals[j];
56: }
57: MatRestoreRow(A,rstart+i,&ncols,&cols,&vals);
58: }
60: /* print compressed matrix */
61: PetscSynchronizedPrintf(MPI_COMM_WORLD,"[%d] Compressed matrix \n",rank);
62: for (i=0; i<nrow; i++) {
63: for (j=0; j<nc; j++) {
64: PetscSynchronizedPrintf(MPI_COMM_WORLD,"%12.4e ",(double)PetscAbsScalar(cm[i+nrow*j]));
65: }
66: PetscSynchronizedPrintf(MPI_COMM_WORLD,"\n");
67: }
68: PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
70: /* put the compressed matrix into the standard matrix */
71: MatDuplicate(A,MAT_COPY_VALUES,&B);
72: MatZeroEntries(A);
73: MatView(B,0);
74: MatFDColoringCreate(A,iscoloring,&fdcoloring);
75: PetscOptionsHasName(NULL,NULL,"-single_block",&single);
76: if (single) {
77: MatFDColoringSetBlockSize(fdcoloring,PETSC_DEFAULT,nc);
78: }
79: PetscOptionsHasName(NULL,NULL,"-two_block",&two);
80: if (two) {
81: MatFDColoringSetBlockSize(fdcoloring,PETSC_DEFAULT,2);
82: }
83: MatFDColoringSetFromOptions(fdcoloring);
84: MatFDColoringSetUp(A,iscoloring,fdcoloring);
86: MatFDColoringSetValues(A,fdcoloring,cm);
87: MatView(A,NULL);
89: /* check the values were put in the correct locations */
90: MatAXPY(A,-1.0,B,SAME_NONZERO_PATTERN);
91: MatView(A,NULL);
92: MatNorm(A,NORM_FROBENIUS,&norm);
93: if (norm > PETSC_MACHINE_EPSILON) {
94: PetscPrintf(PETSC_COMM_WORLD,"Matrix is not identical, problem with MatFDColoringSetValues()\n");
95: }
96: PetscFree(colors);
97: PetscFree(cm);
98: ISColoringDestroy(&iscoloring);
99: MatFDColoringDestroy(&fdcoloring);
100: MatDestroy(&A);
101: MatDestroy(&B);
102: DMDestroy(&da);
103: PetscFinalize();
104: return 0;
105: }
107: /*TEST
109: test:
110: nsize: 2
111: requires: !complex
113: test:
114: suffix: single
115: requires: !complex
116: nsize: 2
117: args: -single_block
118: output_file: output/ex240_1.out
120: test:
121: suffix: two
122: requires: !complex
123: nsize: 2
124: args: -two_block
125: output_file: output/ex240_1.out
127: test:
128: suffix: 2
129: requires: !complex
130: nsize: 5
132: TEST*/