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