1 2 static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat A,Atrans,sA,*submatA,*submatsA; 9 PetscInt bs=1,m=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn; 10 PetscMPIInt size; 11 PetscScalar *vals,rval,one=1.0; 12 IS *is1,*is2; 13 PetscRandom rand; 14 Vec xx,s1,s2; 15 PetscReal s1norm,s2norm,rnorm,tol = 10*PETSC_SMALL; 16 PetscBool flg; 17 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 24 /* create a SeqBAIJ matrix A */ 25 M = m*bs; 26 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A)); 27 PetscCall(MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE)); 28 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rand)); 29 PetscCall(PetscRandomSetFromOptions(rand)); 30 31 PetscCall(PetscMalloc1(bs,&rows)); 32 PetscCall(PetscMalloc1(bs,&cols)); 33 PetscCall(PetscMalloc1(bs*bs,&vals)); 34 PetscCall(PetscMalloc1(M,&idx)); 35 36 /* Now set blocks of random values */ 37 /* first, set diagonal blocks as zero */ 38 for (j=0; j<bs*bs; j++) vals[j] = 0.0; 39 for (i=0; i<m; i++) { 40 cols[0] = i*bs; rows[0] = i*bs; 41 for (j=1; j<bs; j++) { 42 rows[j] = rows[j-1]+1; 43 cols[j] = cols[j-1]+1; 44 } 45 PetscCall(MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES)); 46 } 47 /* second, add random blocks */ 48 for (i=0; i<20*bs; i++) { 49 PetscCall(PetscRandomGetValue(rand,&rval)); 50 cols[0] = bs*(int)(PetscRealPart(rval)*m); 51 PetscCall(PetscRandomGetValue(rand,&rval)); 52 rows[0] = bs*(int)(PetscRealPart(rval)*m); 53 for (j=1; j<bs; j++) { 54 rows[j] = rows[j-1]+1; 55 cols[j] = cols[j-1]+1; 56 } 57 58 for (j=0; j<bs*bs; j++) { 59 PetscCall(PetscRandomGetValue(rand,&rval)); 60 vals[j] = rval; 61 } 62 PetscCall(MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES)); 63 } 64 65 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 66 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 67 68 /* make A a symmetric matrix: A <- A^T + A */ 69 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans)); 70 PetscCall(MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN)); 71 PetscCall(MatDestroy(&Atrans)); 72 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans)); 73 PetscCall(MatEqual(A, Atrans, &flg)); 74 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"A+A^T is non-symmetric"); 75 PetscCall(MatDestroy(&Atrans)); 76 77 /* create a SeqSBAIJ matrix sA (= A) */ 78 PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 79 PetscCall(MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA)); 80 81 /* Test sA==A through MatMult() */ 82 for (i=0; i<nd; i++) { 83 PetscCall(MatGetSize(A,&mm,&nn)); 84 PetscCall(VecCreateSeq(PETSC_COMM_SELF,mm,&xx)); 85 PetscCall(VecDuplicate(xx,&s1)); 86 PetscCall(VecDuplicate(xx,&s2)); 87 for (j=0; j<3; j++) { 88 PetscCall(VecSetRandom(xx,rand)); 89 PetscCall(MatMult(A,xx,s1)); 90 PetscCall(MatMult(sA,xx,s2)); 91 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 92 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 93 rnorm = s2norm-s1norm; 94 if (rnorm<-tol || rnorm>tol) { 95 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm)); 96 } 97 } 98 PetscCall(VecDestroy(&xx)); 99 PetscCall(VecDestroy(&s1)); 100 PetscCall(VecDestroy(&s2)); 101 } 102 103 /* Test MatIncreaseOverlap() */ 104 PetscCall(PetscMalloc1(nd,&is1)); 105 PetscCall(PetscMalloc1(nd,&is2)); 106 107 for (i=0; i<nd; i++) { 108 PetscCall(PetscRandomGetValue(rand,&rval)); 109 size = (int)(PetscRealPart(rval)*m); 110 for (j=0; j<size; j++) { 111 PetscCall(PetscRandomGetValue(rand,&rval)); 112 idx[j*bs] = bs*(int)(PetscRealPart(rval)*m); 113 for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k; 114 } 115 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is1+i)); 116 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is2+i)); 117 } 118 /* for debugging */ 119 /* 120 PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF)); 121 PetscCall(MatView(sA,PETSC_VIEWER_STDOUT_SELF)); 122 */ 123 124 PetscCall(MatIncreaseOverlap(A,nd,is1,ov)); 125 PetscCall(MatIncreaseOverlap(sA,nd,is2,ov)); 126 127 for (i=0; i<nd; ++i) { 128 PetscCall(ISSort(is1[i])); 129 PetscCall(ISSort(is2[i])); 130 } 131 132 for (i=0; i<nd; ++i) { 133 PetscCall(ISEqual(is1[i],is2[i],&flg)); 134 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_PLIB,"i=%" PetscInt_FMT ", is1 != is2",i); 135 } 136 137 PetscCall(MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA)); 138 PetscCall(MatCreateSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA)); 139 140 /* Test MatMult() */ 141 for (i=0; i<nd; i++) { 142 PetscCall(MatGetSize(submatA[i],&mm,&nn)); 143 PetscCall(VecCreateSeq(PETSC_COMM_SELF,mm,&xx)); 144 PetscCall(VecDuplicate(xx,&s1)); 145 PetscCall(VecDuplicate(xx,&s2)); 146 for (j=0; j<3; j++) { 147 PetscCall(VecSetRandom(xx,rand)); 148 PetscCall(MatMult(submatA[i],xx,s1)); 149 PetscCall(MatMult(submatsA[i],xx,s2)); 150 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 151 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 152 rnorm = s2norm-s1norm; 153 if (rnorm<-tol || rnorm>tol) { 154 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm)); 155 } 156 } 157 PetscCall(VecDestroy(&xx)); 158 PetscCall(VecDestroy(&s1)); 159 PetscCall(VecDestroy(&s2)); 160 } 161 162 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 163 PetscCall(MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA)); 164 PetscCall(MatCreateSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA)); 165 166 /* Test MatMult() */ 167 for (i=0; i<nd; i++) { 168 PetscCall(MatGetSize(submatA[i],&mm,&nn)); 169 PetscCall(VecCreateSeq(PETSC_COMM_SELF,mm,&xx)); 170 PetscCall(VecDuplicate(xx,&s1)); 171 PetscCall(VecDuplicate(xx,&s2)); 172 for (j=0; j<3; j++) { 173 PetscCall(VecSetRandom(xx,rand)); 174 PetscCall(MatMult(submatA[i],xx,s1)); 175 PetscCall(MatMult(submatsA[i],xx,s2)); 176 PetscCall(VecNorm(s1,NORM_2,&s1norm)); 177 PetscCall(VecNorm(s2,NORM_2,&s2norm)); 178 rnorm = s2norm-s1norm; 179 if (rnorm<-tol || rnorm>tol) { 180 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm)); 181 } 182 } 183 PetscCall(VecDestroy(&xx)); 184 PetscCall(VecDestroy(&s1)); 185 PetscCall(VecDestroy(&s2)); 186 } 187 188 /* Free allocated memory */ 189 for (i=0; i<nd; ++i) { 190 PetscCall(ISDestroy(&is1[i])); 191 PetscCall(ISDestroy(&is2[i])); 192 } 193 PetscCall(MatDestroySubMatrices(nd,&submatA)); 194 PetscCall(MatDestroySubMatrices(nd,&submatsA)); 195 196 PetscCall(PetscFree(is1)); 197 PetscCall(PetscFree(is2)); 198 PetscCall(PetscFree(idx)); 199 PetscCall(PetscFree(rows)); 200 PetscCall(PetscFree(cols)); 201 PetscCall(PetscFree(vals)); 202 PetscCall(MatDestroy(&A)); 203 PetscCall(MatDestroy(&sA)); 204 PetscCall(PetscRandomDestroy(&rand)); 205 PetscCall(PetscFinalize()); 206 return 0; 207 } 208 209 /*TEST 210 211 test: 212 args: -ov 2 213 214 TEST*/ 215