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