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