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