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,&ltog);
 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*/