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