1 2 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for MatBAIJ format.\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=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn,lsize; 10 PetscScalar *vals,rval; 11 IS *is1,*is2; 12 PetscRandom rdm; 13 Vec xx,s1,s2; 14 PetscReal s1norm,s2norm,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON; 15 PetscBool flg; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 19 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL)); 20 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL)); 21 PetscCall(PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL)); 22 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL)); 23 M = m*bs; 24 25 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A)); 26 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 27 PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B)); 28 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 29 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rdm)); 30 PetscCall(PetscRandomSetFromOptions(rdm)); 31 32 PetscCall(PetscMalloc1(bs,&rows)); 33 PetscCall(PetscMalloc1(bs,&cols)); 34 PetscCall(PetscMalloc1(bs*bs,&vals)); 35 PetscCall(PetscMalloc1(M,&idx)); 36 37 /* Now set blocks of values */ 38 for (i=0; i<20*bs; i++) { 39 PetscCall(PetscRandomGetValue(rdm,&rval)); 40 cols[0] = bs*(int)(PetscRealPart(rval)*m); 41 PetscCall(PetscRandomGetValue(rdm,&rval)); 42 rows[0] = bs*(int)(PetscRealPart(rval)*m); 43 for (j=1; j<bs; j++) { 44 rows[j] = rows[j-1]+1; 45 cols[j] = cols[j-1]+1; 46 } 47 48 for (j=0; j<bs*bs; j++) { 49 PetscCall(PetscRandomGetValue(rdm,&rval)); 50 vals[j] = rval; 51 } 52 PetscCall(MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES)); 53 PetscCall(MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES)); 54 } 55 56 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 57 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 58 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 59 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 60 61 /* Test MatIncreaseOverlap() */ 62 PetscCall(PetscMalloc1(nd,&is1)); 63 PetscCall(PetscMalloc1(nd,&is2)); 64 65 for (i=0; i<nd; i++) { 66 PetscCall(PetscRandomGetValue(rdm,&rval)); 67 lsize = (int)(PetscRealPart(rval)*m); 68 for (j=0; j<lsize; j++) { 69 PetscCall(PetscRandomGetValue(rdm,&rval)); 70 idx[j*bs] = bs*(int)(PetscRealPart(rval)*m); 71 for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k; 72 } 73 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,PETSC_COPY_VALUES,is1+i)); 74 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,PETSC_COPY_VALUES,is2+i)); 75 } 76 PetscCall(MatIncreaseOverlap(A,nd,is1,ov)); 77 PetscCall(MatIncreaseOverlap(B,nd,is2,ov)); 78 79 for (i=0; i<nd; ++i) { 80 PetscCall(ISEqual(is1[i],is2[i],&flg)); 81 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"i=%" PetscInt_FMT ", flg =%d",i,(int)flg); 82 } 83 84 for (i=0; i<nd; ++i) { 85 PetscCall(ISSort(is1[i])); 86 PetscCall(ISSort(is2[i])); 87 } 88 89 PetscCall(MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA)); 90 PetscCall(MatCreateSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB)); 91 92 /* Test MatMult() */ 93 for (i=0; i<nd; i++) { 94 PetscCall(MatGetSize(submatA[i],&mm,&nn)); 95 PetscCall(VecCreateSeq(PETSC_COMM_SELF,mm,&xx)); 96 PetscCall(VecDuplicate(xx,&s1)); 97 PetscCall(VecDuplicate(xx,&s2)); 98 for (j=0; j<3; j++) { 99 PetscCall(VecSetRandom(xx,rdm)); 100 PetscCall(MatMult(submatA[i],xx,s1)); 101 PetscCall(MatMult(submatB[i],xx,s2)); 102 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 103 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 104 rnorm = s2norm-s1norm; 105 if (rnorm<-tol || rnorm>tol) { 106 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm)); 107 } 108 } 109 PetscCall(VecDestroy(&xx)); 110 PetscCall(VecDestroy(&s1)); 111 PetscCall(VecDestroy(&s2)); 112 } 113 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 114 PetscCall(MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA)); 115 PetscCall(MatCreateSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB)); 116 117 /* Test MatMult() */ 118 for (i=0; i<nd; i++) { 119 PetscCall(MatGetSize(submatA[i],&mm,&nn)); 120 PetscCall(VecCreateSeq(PETSC_COMM_SELF,mm,&xx)); 121 PetscCall(VecDuplicate(xx,&s1)); 122 PetscCall(VecDuplicate(xx,&s2)); 123 for (j=0; j<3; j++) { 124 PetscCall(VecSetRandom(xx,rdm)); 125 PetscCall(MatMult(submatA[i],xx,s1)); 126 PetscCall(MatMult(submatB[i],xx,s2)); 127 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 128 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 129 rnorm = s2norm-s1norm; 130 if (rnorm<-tol || rnorm>tol) { 131 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm)); 132 } 133 } 134 PetscCall(VecDestroy(&xx)); 135 PetscCall(VecDestroy(&s1)); 136 PetscCall(VecDestroy(&s2)); 137 } 138 139 /* Free allocated memory */ 140 for (i=0; i<nd; ++i) { 141 PetscCall(ISDestroy(&is1[i])); 142 PetscCall(ISDestroy(&is2[i])); 143 } 144 PetscCall(MatDestroySubMatrices(nd,&submatA)); 145 PetscCall(MatDestroySubMatrices(nd,&submatB)); 146 PetscCall(PetscFree(is1)); 147 PetscCall(PetscFree(is2)); 148 PetscCall(PetscFree(idx)); 149 PetscCall(PetscFree(rows)); 150 PetscCall(PetscFree(cols)); 151 PetscCall(PetscFree(vals)); 152 PetscCall(MatDestroy(&A)); 153 PetscCall(MatDestroy(&B)); 154 PetscCall(PetscRandomDestroy(&rdm)); 155 PetscCall(PetscFinalize()); 156 return 0; 157 } 158 159 /*TEST 160 161 test: 162 args: -mat_block_size {{1 2 5 7 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}} 163 164 TEST*/ 165