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