1 2 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel AIJ and BAIJ formats.\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat A,B,*submatA,*submatB; 9 PetscInt bs=1,m=11,ov=1,i,j,k,*rows,*cols,nd=5,*idx,rstart,rend,sz,mm,nn,M,N,Mbs; 10 PetscMPIInt size,rank; 11 PetscScalar *vals,rval; 12 IS *is1,*is2; 13 PetscRandom rdm; 14 Vec xx,s1,s2; 15 PetscReal s1norm,s2norm,rnorm,tol = 100*PETSC_SMALL; 16 PetscBool flg,test_nd0=PETSC_FALSE; 17 18 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 20 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 21 22 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL)); 23 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL)); 24 PetscCall(PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL)); 25 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL)); 26 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_nd0",&test_nd0,NULL)); 27 28 /* Create a AIJ matrix A */ 29 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 30 PetscCall(MatSetSizes(A,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE)); 31 PetscCall(MatSetType(A,MATAIJ)); 32 PetscCall(MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,NULL)); 33 PetscCall(MatMPIAIJSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL)); 34 PetscCall(MatSetFromOptions(A)); 35 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 36 37 /* Create a BAIJ matrix B */ 38 PetscCall(MatCreate(PETSC_COMM_WORLD,&B)); 39 PetscCall(MatSetSizes(B,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE)); 40 PetscCall(MatSetType(B,MATBAIJ)); 41 PetscCall(MatSeqBAIJSetPreallocation(B,bs,PETSC_DEFAULT,NULL)); 42 PetscCall(MatMPIBAIJSetPreallocation(B,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL)); 43 PetscCall(MatSetFromOptions(B)); 44 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 45 46 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rdm)); 47 PetscCall(PetscRandomSetFromOptions(rdm)); 48 49 PetscCall(MatGetOwnershipRange(A,&rstart,&rend)); 50 PetscCall(MatGetSize(A,&M,&N)); 51 Mbs = M/bs; 52 53 PetscCall(PetscMalloc1(bs,&rows)); 54 PetscCall(PetscMalloc1(bs,&cols)); 55 PetscCall(PetscMalloc1(bs*bs,&vals)); 56 PetscCall(PetscMalloc1(M,&idx)); 57 58 /* Now set blocks of values */ 59 for (i=0; i<40*bs; i++) { 60 PetscCall(PetscRandomGetValue(rdm,&rval)); 61 cols[0] = bs*(int)(PetscRealPart(rval)*Mbs); 62 PetscCall(PetscRandomGetValue(rdm,&rval)); 63 rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m); 64 for (j=1; j<bs; j++) { 65 rows[j] = rows[j-1]+1; 66 cols[j] = cols[j-1]+1; 67 } 68 69 for (j=0; j<bs*bs; j++) { 70 PetscCall(PetscRandomGetValue(rdm,&rval)); 71 vals[j] = rval; 72 } 73 PetscCall(MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES)); 74 PetscCall(MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES)); 75 } 76 77 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 78 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 79 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 80 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 81 82 /* Test MatIncreaseOverlap() */ 83 PetscCall(PetscMalloc1(nd,&is1)); 84 PetscCall(PetscMalloc1(nd,&is2)); 85 86 if (rank == 0 && test_nd0) nd = 0; /* test case */ 87 88 for (i=0; i<nd; i++) { 89 PetscCall(PetscRandomGetValue(rdm,&rval)); 90 sz = (int)(PetscRealPart(rval)*m); 91 for (j=0; j<sz; j++) { 92 PetscCall(PetscRandomGetValue(rdm,&rval)); 93 idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs); 94 for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k; 95 } 96 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i)); 97 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i)); 98 } 99 PetscCall(MatIncreaseOverlap(A,nd,is1,ov)); 100 PetscCall(MatIncreaseOverlap(B,nd,is2,ov)); 101 102 for (i=0; i<nd; ++i) { 103 PetscCall(ISEqual(is1[i],is2[i],&flg)); 104 105 if (!flg) { 106 PetscCall(PetscPrintf(PETSC_COMM_SELF,"i=%" PetscInt_FMT ", flg=%d :bs=%" PetscInt_FMT " m=%" PetscInt_FMT " ov=%" PetscInt_FMT " nd=%" PetscInt_FMT " np=%d\n",i,flg,bs,m,ov,nd,size)); 107 } 108 } 109 110 for (i=0; i<nd; ++i) { 111 PetscCall(ISSort(is1[i])); 112 PetscCall(ISSort(is2[i])); 113 } 114 115 PetscCall(MatCreateSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB)); 116 PetscCall(MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA)); 117 118 /* Test MatMult() */ 119 for (i=0; i<nd; i++) { 120 PetscCall(MatGetSize(submatA[i],&mm,&nn)); 121 PetscCall(VecCreateSeq(PETSC_COMM_SELF,mm,&xx)); 122 PetscCall(VecDuplicate(xx,&s1)); 123 PetscCall(VecDuplicate(xx,&s2)); 124 for (j=0; j<3; j++) { 125 PetscCall(VecSetRandom(xx,rdm)); 126 PetscCall(MatMult(submatA[i],xx,s1)); 127 PetscCall(MatMult(submatB[i],xx,s2)); 128 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 129 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 130 rnorm = s2norm-s1norm; 131 if (rnorm<-tol || rnorm>tol) { 132 PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,(double)s1norm,(double)s2norm)); 133 } 134 } 135 PetscCall(VecDestroy(&xx)); 136 PetscCall(VecDestroy(&s1)); 137 PetscCall(VecDestroy(&s2)); 138 } 139 140 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 141 PetscCall(MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA)); 142 PetscCall(MatCreateSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB)); 143 144 /* Test MatMult() */ 145 for (i=0; i<nd; i++) { 146 PetscCall(MatGetSize(submatA[i],&mm,&nn)); 147 PetscCall(VecCreateSeq(PETSC_COMM_SELF,mm,&xx)); 148 PetscCall(VecDuplicate(xx,&s1)); 149 PetscCall(VecDuplicate(xx,&s2)); 150 for (j=0; j<3; j++) { 151 PetscCall(VecSetRandom(xx,rdm)); 152 PetscCall(MatMult(submatA[i],xx,s1)); 153 PetscCall(MatMult(submatB[i],xx,s2)); 154 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 155 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 156 rnorm = s2norm-s1norm; 157 if (rnorm<-tol || rnorm>tol) { 158 PetscCall(PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,(double)s1norm,(double)s2norm)); 159 } 160 } 161 PetscCall(VecDestroy(&xx)); 162 PetscCall(VecDestroy(&s1)); 163 PetscCall(VecDestroy(&s2)); 164 } 165 166 /* Free allocated memory */ 167 for (i=0; i<nd; ++i) { 168 PetscCall(ISDestroy(&is1[i])); 169 PetscCall(ISDestroy(&is2[i])); 170 } 171 PetscCall(MatDestroySubMatrices(nd,&submatA)); 172 PetscCall(MatDestroySubMatrices(nd,&submatB)); 173 174 PetscCall(PetscFree(is1)); 175 PetscCall(PetscFree(is2)); 176 PetscCall(PetscFree(idx)); 177 PetscCall(PetscFree(rows)); 178 PetscCall(PetscFree(cols)); 179 PetscCall(PetscFree(vals)); 180 PetscCall(MatDestroy(&A)); 181 PetscCall(MatDestroy(&B)); 182 PetscCall(PetscRandomDestroy(&rdm)); 183 PetscCall(PetscFinalize()); 184 return 0; 185 } 186 187 /*TEST 188 189 test: 190 nsize: {{1 3}} 191 args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}} ; done 192 output_file: output/ex54.out 193 194 test: 195 suffix: 2 196 args: -nd 2 -test_nd0 197 output_file: output/ex54.out 198 199 test: 200 suffix: 3 201 nsize: 3 202 args: -nd 2 -test_nd0 203 output_file: output/ex54.out 204 205 TEST*/ 206