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