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