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: PetscErrorCode ierr;
10: PetscInt N, mx = 5,my = 4,i,j,nc,nrow,n,ncols,rstart,*colors,*map;
11: const PetscInt *cols;
12: const PetscScalar *vals;
13: Mat A,B;
14: PetscReal norm;
15: MatFDColoring fdcoloring;
16: ISColoring iscoloring;
17: PetscScalar *cm;
18: const ISColoringValue *icolors;
19: PetscMPIInt rank;
20: ISLocalToGlobalMapping ltog;
21: PetscBool single,two;
23: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
24: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
25: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,mx,my,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
26: DMSetUp(da);
27: DMCreateMatrix(da,&A);
29: /* 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 */
30: /* first put the coloring in the global ordering */
31: DMCreateColoring(da,IS_COLORING_LOCAL,&iscoloring);
32: ISColoringGetColors(iscoloring,&n,&nc,&icolors);
33: DMGetLocalToGlobalMapping(da,<og);
34: PetscMalloc1(n,&map);
35: for (i=0; i<n; i++) map[i] = i;
36: ISLocalToGlobalMappingApply(ltog,n,map,map);
37: MatGetSize(A,&N,NULL);
38: PetscMalloc1(N,&colors);
39: for (i=0; i<N; i++) colors[i] = -1;
40: for (i=0; i<n; i++) colors[map[i]]= icolors[i];
41: PetscFree(map);
42: PetscSynchronizedPrintf(MPI_COMM_WORLD,"[%d]Global colors \n",rank);
43: for (i=0; i<N; i++) {PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%D %D\n",i,colors[i]);}
44: PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
46: /* second, compress the A matrix */
47: MatSetRandom(A,NULL);
48: MatView(A,NULL);
49: MatGetLocalSize(A,&nrow,NULL);
50: MatGetOwnershipRange(A,&rstart,NULL);
51: PetscCalloc1(nrow*nc,&cm);
52: for (i=0; i<nrow; i++) {
53: MatGetRow(A,rstart+i,&ncols,&cols,&vals);
54: for (j=0; j<ncols; j++) {
55: if (colors[cols[j]] < 0) SETERRQ1(PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Global column %D had no color",cols[j]);
56: cm[i + nrow*colors[cols[j]]] = vals[j];
57: }
58: MatRestoreRow(A,rstart+i,&ncols,&cols,&vals);
59: }
61: /* print compressed matrix */
62: PetscSynchronizedPrintf(MPI_COMM_WORLD,"[%d] Compressed matrix \n",rank);
63: for (i=0; i<nrow; i++) {
64: for (j=0; j<nc; j++) {
65: PetscSynchronizedPrintf(MPI_COMM_WORLD,"%12.4e ",cm[i+nrow*j]);
66: }
67: PetscSynchronizedPrintf(MPI_COMM_WORLD,"\n");
68: }
69: PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout);
71: /* put the compressed matrix into the standard matrix */
72: MatDuplicate(A,MAT_COPY_VALUES,&B);
73: MatZeroEntries(A);
74: MatView(B,0);
75: MatFDColoringCreate(A,iscoloring,&fdcoloring);
76: PetscOptionsHasName(NULL,NULL,"-single_block",&single);
77: if (single) {
78: MatFDColoringSetBlockSize(fdcoloring,PETSC_DEFAULT,nc);
79: }
80: PetscOptionsHasName(NULL,NULL,"-two_block",&two);
81: if (two) {
82: MatFDColoringSetBlockSize(fdcoloring,PETSC_DEFAULT,2);
83: }
84: MatFDColoringSetFromOptions(fdcoloring);
85: MatFDColoringSetUp(A,iscoloring,fdcoloring);
87: MatFDColoringSetValues(A,fdcoloring,cm);
88: MatView(A,NULL);
90: /* check the values were put in the correct locations */
91: MatAXPY(A,-1.0,B,SAME_NONZERO_PATTERN);
92: MatView(A,NULL);
93: MatNorm(A,NORM_FROBENIUS,&norm);
94: if (norm > PETSC_MACHINE_EPSILON) {
95: PetscPrintf(PETSC_COMM_WORLD,"Matrix is not identical, problem with MatFDColoringSetValues()\n");
96: }
97: PetscFree(colors);
98: PetscFree(cm);
99: ISColoringDestroy(&iscoloring);
100: MatFDColoringDestroy(&fdcoloring);
101: MatDestroy(&A);
102: MatDestroy(&B);
103: DMDestroy(&da);
104: PetscFinalize();
105: return ierr;
106: }
108: /*TEST
110: test:
111: nsize: 2
112: requires: !complex
114: test:
115: suffix: single
116: requires: !complex
117: nsize: 2
118: args: -single_block
119: output_file: output/ex240_1.out
121: test:
122: suffix: two
123: requires: !complex
124: nsize: 2
125: args: -two_block
126: output_file: output/ex240_1.out
128: test:
129: suffix: 2
130: requires: !complex
131: nsize: 5
133: TEST*/