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