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