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