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 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 25 CHKERRQ(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 CHKERRQ(DMSetUp(da)); 27 CHKERRQ(DMCreateMatrix(da,&A)); 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 CHKERRQ(DMCreateColoring(da,IS_COLORING_LOCAL,&iscoloring)); 32 CHKERRQ(ISColoringGetColors(iscoloring,&n,&nc,&icolors)); 33 CHKERRQ(DMGetLocalToGlobalMapping(da,<og)); 34 CHKERRQ(PetscMalloc1(n,&map)); 35 for (i=0; i<n; i++) map[i] = i; 36 CHKERRQ(ISLocalToGlobalMappingApply(ltog,n,map,map)); 37 CHKERRQ(MatGetSize(A,&N,NULL)); 38 CHKERRQ(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 CHKERRQ(PetscFree(map)); 42 CHKERRQ(PetscSynchronizedPrintf(MPI_COMM_WORLD,"[%d]Global colors \n",rank)); 43 for (i=0; i<N; i++) CHKERRQ(PetscSynchronizedPrintf(PETSC_COMM_WORLD,"%" PetscInt_FMT " %" PetscInt_FMT "\n",i,colors[i])); 44 CHKERRQ(PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout)); 45 46 /* second, compress the A matrix */ 47 CHKERRQ(MatSetRandom(A,NULL)); 48 CHKERRQ(MatView(A,NULL)); 49 CHKERRQ(MatGetLocalSize(A,&nrow,NULL)); 50 CHKERRQ(MatGetOwnershipRange(A,&rstart,NULL)); 51 CHKERRQ(PetscCalloc1(nrow*nc,&cm)); 52 for (i=0; i<nrow; i++) { 53 CHKERRQ(MatGetRow(A,rstart+i,&ncols,&cols,&vals)); 54 for (j=0; j<ncols; j++) { 55 PetscCheckFalse(colors[cols[j]] < 0,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Global column %" PetscInt_FMT " had no color",cols[j]); 56 cm[i + nrow*colors[cols[j]]] = vals[j]; 57 } 58 CHKERRQ(MatRestoreRow(A,rstart+i,&ncols,&cols,&vals)); 59 } 60 61 /* print compressed matrix */ 62 CHKERRQ(PetscSynchronizedPrintf(MPI_COMM_WORLD,"[%d] Compressed matrix \n",rank)); 63 for (i=0; i<nrow; i++) { 64 for (j=0; j<nc; j++) { 65 CHKERRQ(PetscSynchronizedPrintf(MPI_COMM_WORLD,"%12.4e ",(double)PetscAbsScalar(cm[i+nrow*j]))); 66 } 67 CHKERRQ(PetscSynchronizedPrintf(MPI_COMM_WORLD,"\n")); 68 } 69 CHKERRQ(PetscSynchronizedFlush(PETSC_COMM_WORLD,stdout)); 70 71 /* put the compressed matrix into the standard matrix */ 72 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&B)); 73 CHKERRQ(MatZeroEntries(A)); 74 CHKERRQ(MatView(B,0)); 75 CHKERRQ(MatFDColoringCreate(A,iscoloring,&fdcoloring)); 76 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-single_block",&single)); 77 if (single) { 78 CHKERRQ(MatFDColoringSetBlockSize(fdcoloring,PETSC_DEFAULT,nc)); 79 } 80 CHKERRQ(PetscOptionsHasName(NULL,NULL,"-two_block",&two)); 81 if (two) { 82 CHKERRQ(MatFDColoringSetBlockSize(fdcoloring,PETSC_DEFAULT,2)); 83 } 84 CHKERRQ(MatFDColoringSetFromOptions(fdcoloring)); 85 CHKERRQ(MatFDColoringSetUp(A,iscoloring,fdcoloring)); 86 87 CHKERRQ(MatFDColoringSetValues(A,fdcoloring,cm)); 88 CHKERRQ(MatView(A,NULL)); 89 90 /* check the values were put in the correct locations */ 91 CHKERRQ(MatAXPY(A,-1.0,B,SAME_NONZERO_PATTERN)); 92 CHKERRQ(MatView(A,NULL)); 93 CHKERRQ(MatNorm(A,NORM_FROBENIUS,&norm)); 94 if (norm > PETSC_MACHINE_EPSILON) { 95 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Matrix is not identical, problem with MatFDColoringSetValues()\n")); 96 } 97 CHKERRQ(PetscFree(colors)); 98 CHKERRQ(PetscFree(cm)); 99 CHKERRQ(ISColoringDestroy(&iscoloring)); 100 CHKERRQ(MatFDColoringDestroy(&fdcoloring)); 101 CHKERRQ(MatDestroy(&A)); 102 CHKERRQ(MatDestroy(&B)); 103 CHKERRQ(DMDestroy(&da)); 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