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