1 2 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel MatSBAIJ format.\n"; 3 /* Example of usage: 4 mpiexec -n 2 ./ex92 -nd 2 -ov 3 -mat_block_size 2 -view_id 0 -test_overlap -test_submat 5 */ 6 #include <petscmat.h> 7 8 int main(int argc,char **args) 9 { 10 Mat A,Atrans,sA,*submatA,*submatsA; 11 PetscMPIInt size,rank; 12 PetscInt bs=1,mbs=10,ov=1,i,j,k,*rows,*cols,nd=2,*idx,rstart,rend,sz,M,N,Mbs; 13 PetscScalar *vals,rval,one=1.0; 14 IS *is1,*is2; 15 PetscRandom rand; 16 PetscBool flg,TestOverlap,TestSubMat,TestAllcols,test_sorted=PETSC_FALSE; 17 PetscInt vid = -1; 18 #if defined(PETSC_USE_LOG) 19 PetscLogStage stages[2]; 20 #endif 21 22 PetscFunctionBeginUser; 23 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 24 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 25 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 26 27 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL)); 28 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mat_mbs",&mbs,NULL)); 29 PetscCall(PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL)); 30 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL)); 31 PetscCall(PetscOptionsGetInt(NULL,NULL,"-view_id",&vid,NULL)); 32 PetscCall(PetscOptionsHasName(NULL,NULL, "-test_overlap", &TestOverlap)); 33 PetscCall(PetscOptionsHasName(NULL,NULL, "-test_submat", &TestSubMat)); 34 PetscCall(PetscOptionsHasName(NULL,NULL, "-test_allcols", &TestAllcols)); 35 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_sorted",&test_sorted,NULL)); 36 37 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 38 PetscCall(MatSetSizes(A,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE)); 39 PetscCall(MatSetType(A,MATBAIJ)); 40 PetscCall(MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL)); 41 PetscCall(MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL)); 42 43 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 44 PetscCall(PetscRandomSetFromOptions(rand)); 45 46 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 47 PetscCall(MatGetSize(A,&M,&N)); 48 Mbs = M/bs; 49 50 PetscCall(PetscMalloc1(bs,&rows)); 51 PetscCall(PetscMalloc1(bs,&cols)); 52 PetscCall(PetscMalloc1(bs*bs,&vals)); 53 PetscCall(PetscMalloc1(M,&idx)); 54 55 /* Now set blocks of values */ 56 for (j=0; j<bs*bs; j++) vals[j] = 0.0; 57 for (i=0; i<Mbs; i++) { 58 cols[0] = i*bs; rows[0] = i*bs; 59 for (j=1; j<bs; j++) { 60 rows[j] = rows[j-1]+1; 61 cols[j] = cols[j-1]+1; 62 } 63 PetscCall(MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES)); 64 } 65 /* second, add random blocks */ 66 for (i=0; i<20*bs; i++) { 67 PetscCall(PetscRandomGetValue(rand,&rval)); 68 cols[0] = bs*(PetscInt)(PetscRealPart(rval)*Mbs); 69 PetscCall(PetscRandomGetValue(rand,&rval)); 70 rows[0] = rstart + bs*(PetscInt)(PetscRealPart(rval)*mbs); 71 for (j=1; j<bs; j++) { 72 rows[j] = rows[j-1]+1; 73 cols[j] = cols[j-1]+1; 74 } 75 76 for (j=0; j<bs*bs; j++) { 77 PetscCall(PetscRandomGetValue(rand,&rval)); 78 vals[j] = rval; 79 } 80 PetscCall(MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES)); 81 } 82 83 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 84 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 85 86 /* make A a symmetric matrix: A <- A^T + A */ 87 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans)); 88 PetscCall(MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN)); 89 PetscCall(MatDestroy(&Atrans)); 90 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans)); 91 PetscCall(MatEqual(A, Atrans, &flg)); 92 if (flg) { 93 PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 94 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A+A^T is non-symmetric"); 95 PetscCall(MatDestroy(&Atrans)); 96 97 /* create a SeqSBAIJ matrix sA (= A) */ 98 PetscCall(MatConvert(A,MATSBAIJ,MAT_INITIAL_MATRIX,&sA)); 99 if (vid >= 0 && vid < size) { 100 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"A:\n")); 101 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_WORLD)); 102 PetscCall(PetscViewerASCIIPrintf(PETSC_VIEWER_STDOUT_WORLD,"sA:\n")); 103 PetscCall(MatView(sA,PETSC_VIEWER_STDOUT_WORLD)); 104 } 105 106 /* Test sA==A through MatMult() */ 107 PetscCall(MatMultEqual(A,sA,10,&flg)); 108 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Error in MatConvert(): A != sA"); 109 110 /* Test MatIncreaseOverlap() */ 111 PetscCall(PetscMalloc1(nd,&is1)); 112 PetscCall(PetscMalloc1(nd,&is2)); 113 114 for (i=0; i<nd; i++) { 115 if (!TestAllcols) { 116 PetscCall(PetscRandomGetValue(rand,&rval)); 117 sz = (PetscInt)((0.5+0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */ 118 119 for (j=0; j<sz; j++) { 120 PetscCall(PetscRandomGetValue(rand,&rval)); 121 idx[j*bs] = bs*(PetscInt)(PetscRealPart(rval)*Mbs); 122 for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k; 123 } 124 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i)); 125 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i)); 126 if (rank == vid) { 127 PetscCall(PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%" PetscInt_FMT "]: %" PetscInt_FMT "\n",rank,i,sz)); 128 PetscCall(ISView(is2[i],PETSC_VIEWER_STDOUT_SELF)); 129 } 130 } else { /* Test all rows and columns */ 131 sz = M; 132 PetscCall(ISCreateStride(PETSC_COMM_SELF,sz,0,1,is1+i)); 133 PetscCall(ISCreateStride(PETSC_COMM_SELF,sz,0,1,is2+i)); 134 135 if (rank == vid) { 136 PetscBool colflag; 137 PetscCall(ISIdentity(is2[i],&colflag)); 138 PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d] is2[%" PetscInt_FMT "], colflag %d\n",rank,i,colflag)); 139 PetscCall(ISView(is2[i],PETSC_VIEWER_STDOUT_SELF)); 140 } 141 } 142 } 143 144 PetscCall(PetscLogStageRegister("MatOv_SBAIJ",&stages[0])); 145 PetscCall(PetscLogStageRegister("MatOv_BAIJ",&stages[1])); 146 147 /* Test MatIncreaseOverlap */ 148 if (TestOverlap) { 149 PetscCall(PetscLogStagePush(stages[0])); 150 PetscCall(MatIncreaseOverlap(sA,nd,is2,ov)); 151 PetscCall(PetscLogStagePop()); 152 153 PetscCall(PetscLogStagePush(stages[1])); 154 PetscCall(MatIncreaseOverlap(A,nd,is1,ov)); 155 PetscCall(PetscLogStagePop()); 156 157 if (rank == vid) { 158 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n[%d] IS from BAIJ:\n",rank)); 159 PetscCall(ISView(is1[0],PETSC_VIEWER_STDOUT_SELF)); 160 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n[%d] IS from SBAIJ:\n",rank)); 161 PetscCall(ISView(is2[0],PETSC_VIEWER_STDOUT_SELF)); 162 } 163 164 for (i=0; i<nd; ++i) { 165 PetscCall(ISEqual(is1[i],is2[i],&flg)); 166 if (!flg) { 167 if (rank == 0) { 168 PetscCall(ISSort(is1[i])); 169 PetscCall(ISSort(is2[i])); 170 } 171 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"i=%" PetscInt_FMT ", is1 != is2",i); 172 } 173 } 174 } 175 176 /* Test MatCreateSubmatrices */ 177 if (TestSubMat) { 178 if (test_sorted) { 179 for (i = 0; i < nd; ++i) { 180 PetscCall(ISSort(is1[i])); 181 } 182 } 183 PetscCall(MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA)); 184 PetscCall(MatCreateSubMatrices(sA,nd,is1,is1,MAT_INITIAL_MATRIX,&submatsA)); 185 186 PetscCall(MatMultEqual(A,sA,10,&flg)); 187 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"A != sA"); 188 189 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 190 PetscCall(MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA)); 191 PetscCall(MatCreateSubMatrices(sA,nd,is1,is1,MAT_REUSE_MATRIX,&submatsA)); 192 PetscCall(MatMultEqual(A,sA,10,&flg)); 193 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"MatCreateSubmatrices(): A != sA"); 194 195 PetscCall(MatDestroySubMatrices(nd,&submatA)); 196 PetscCall(MatDestroySubMatrices(nd,&submatsA)); 197 } 198 199 /* Free allocated memory */ 200 for (i=0; i<nd; ++i) { 201 PetscCall(ISDestroy(&is1[i])); 202 PetscCall(ISDestroy(&is2[i])); 203 } 204 PetscCall(PetscFree(is1)); 205 PetscCall(PetscFree(is2)); 206 PetscCall(PetscFree(idx)); 207 PetscCall(PetscFree(rows)); 208 PetscCall(PetscFree(cols)); 209 PetscCall(PetscFree(vals)); 210 PetscCall(MatDestroy(&A)); 211 PetscCall(MatDestroy(&sA)); 212 PetscCall(PetscRandomDestroy(&rand)); 213 PetscCall(PetscFinalize()); 214 return 0; 215 } 216 217 /*TEST 218 219 test: 220 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat 221 output_file: output/ex92_1.out 222 223 test: 224 suffix: 2 225 nsize: {{3 4}} 226 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_submat 227 output_file: output/ex92_1.out 228 229 test: 230 suffix: 3 231 nsize: {{3 4}} 232 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols 233 output_file: output/ex92_1.out 234 235 test: 236 suffix: 3_sorted 237 nsize: {{3 4}} 238 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_overlap -test_allcols -test_sorted 239 output_file: output/ex92_1.out 240 241 test: 242 suffix: 4 243 nsize: {{3 4}} 244 args: -ov {{1 3}} -mat_block_size {{2 8}} -test_submat -test_allcols 245 output_file: output/ex92_1.out 246 247 TEST*/ 248