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