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