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