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