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) PetscCall(MatFDColoringSetBlockSize(fdcoloring,PETSC_DEFAULT,nc)); 77 PetscCall(PetscOptionsHasName(NULL,NULL,"-two_block",&two)); 78 if (two) PetscCall(MatFDColoringSetBlockSize(fdcoloring,PETSC_DEFAULT,2)); 79 PetscCall(MatFDColoringSetFromOptions(fdcoloring)); 80 PetscCall(MatFDColoringSetUp(A,iscoloring,fdcoloring)); 81 82 PetscCall(MatFDColoringSetValues(A,fdcoloring,cm)); 83 PetscCall(MatView(A,NULL)); 84 85 /* check the values were put in the correct locations */ 86 PetscCall(MatAXPY(A,-1.0,B,SAME_NONZERO_PATTERN)); 87 PetscCall(MatView(A,NULL)); 88 PetscCall(MatNorm(A,NORM_FROBENIUS,&norm)); 89 if (norm > PETSC_MACHINE_EPSILON) { 90 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Matrix is not identical, problem with MatFDColoringSetValues()\n")); 91 } 92 PetscCall(PetscFree(colors)); 93 PetscCall(PetscFree(cm)); 94 PetscCall(ISColoringDestroy(&iscoloring)); 95 PetscCall(MatFDColoringDestroy(&fdcoloring)); 96 PetscCall(MatDestroy(&A)); 97 PetscCall(MatDestroy(&B)); 98 PetscCall(DMDestroy(&da)); 99 PetscCall(PetscFinalize()); 100 return 0; 101 } 102 103 /*TEST 104 105 test: 106 nsize: 2 107 requires: !complex 108 109 test: 110 suffix: single 111 requires: !complex 112 nsize: 2 113 args: -single_block 114 output_file: output/ex240_1.out 115 116 test: 117 suffix: two 118 requires: !complex 119 nsize: 2 120 args: -two_block 121 output_file: output/ex240_1.out 122 123 test: 124 suffix: 2 125 requires: !complex 126 nsize: 5 127 128 TEST*/ 129