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 for (i=0; i<nd; i++) { 109 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr); 110 size = (int)(PetscRealPart(rval)*m); 111 for (j=0; j<size; j++) { 112 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr); 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 ierr = ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is1+i);CHKERRQ(ierr); 117 ierr = ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is2+i);CHKERRQ(ierr); 118 } 119 /* for debugging */ 120 /* 121 ierr = MatView(A,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 122 ierr = MatView(sA,PETSC_VIEWER_STDOUT_SELF);CHKERRQ(ierr); 123 */ 124 125 ierr = MatIncreaseOverlap(A,nd,is1,ov);CHKERRQ(ierr); 126 ierr = MatIncreaseOverlap(sA,nd,is2,ov);CHKERRQ(ierr); 127 128 for (i=0; i<nd; ++i) { 129 ierr = ISSort(is1[i]);CHKERRQ(ierr); 130 ierr = ISSort(is2[i]);CHKERRQ(ierr); 131 } 132 133 for (i=0; i<nd; ++i) { 134 ierr = ISEqual(is1[i],is2[i],&flg);CHKERRQ(ierr); 135 if (!flg) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"i=%d, is1 != is2",i); 136 } 137 138 ierr = MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);CHKERRQ(ierr); 139 ierr = MatCreateSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);CHKERRQ(ierr); 140 141 /* Test MatMult() */ 142 for (i=0; i<nd; i++) { 143 ierr = MatGetSize(submatA[i],&mm,&nn);CHKERRQ(ierr); 144 ierr = VecCreateSeq(PETSC_COMM_SELF,mm,&xx);CHKERRQ(ierr); 145 ierr = VecDuplicate(xx,&s1);CHKERRQ(ierr); 146 ierr = VecDuplicate(xx,&s2);CHKERRQ(ierr); 147 for (j=0; j<3; j++) { 148 ierr = VecSetRandom(xx,rand);CHKERRQ(ierr); 149 ierr = MatMult(submatA[i],xx,s1);CHKERRQ(ierr); 150 ierr = MatMult(submatsA[i],xx,s2);CHKERRQ(ierr); 151 ierr = VecNorm(s1,NORM_2,&s1norm);CHKERRQ(ierr); 152 ierr = VecNorm(s2,NORM_2,&s2norm);CHKERRQ(ierr); 153 rnorm = s2norm-s1norm; 154 if (rnorm<-tol || rnorm>tol) { 155 ierr = PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);CHKERRQ(ierr); 156 } 157 } 158 ierr = VecDestroy(&xx);CHKERRQ(ierr); 159 ierr = VecDestroy(&s1);CHKERRQ(ierr); 160 ierr = VecDestroy(&s2);CHKERRQ(ierr); 161 } 162 163 /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */ 164 ierr = MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);CHKERRQ(ierr); 165 ierr = MatCreateSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);CHKERRQ(ierr); 166 167 /* Test MatMult() */ 168 for (i=0; i<nd; i++) { 169 ierr = MatGetSize(submatA[i],&mm,&nn);CHKERRQ(ierr); 170 ierr = VecCreateSeq(PETSC_COMM_SELF,mm,&xx);CHKERRQ(ierr); 171 ierr = VecDuplicate(xx,&s1);CHKERRQ(ierr); 172 ierr = VecDuplicate(xx,&s2);CHKERRQ(ierr); 173 for (j=0; j<3; j++) { 174 ierr = VecSetRandom(xx,rand);CHKERRQ(ierr); 175 ierr = MatMult(submatA[i],xx,s1);CHKERRQ(ierr); 176 ierr = MatMult(submatsA[i],xx,s2);CHKERRQ(ierr); 177 ierr = VecNorm(s1,NORM_2,&s1norm);CHKERRQ(ierr); 178 ierr = VecNorm(s2,NORM_2,&s2norm);CHKERRQ(ierr); 179 rnorm = s2norm-s1norm; 180 if (rnorm<-tol || rnorm>tol) { 181 ierr = PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);CHKERRQ(ierr); 182 } 183 } 184 ierr = VecDestroy(&xx);CHKERRQ(ierr); 185 ierr = VecDestroy(&s1);CHKERRQ(ierr); 186 ierr = VecDestroy(&s2);CHKERRQ(ierr); 187 } 188 189 /* Free allocated memory */ 190 for (i=0; i<nd; ++i) { 191 ierr = ISDestroy(&is1[i]);CHKERRQ(ierr); 192 ierr = ISDestroy(&is2[i]);CHKERRQ(ierr); 193 } 194 ierr = MatDestroySubMatrices(nd,&submatA);CHKERRQ(ierr); 195 ierr = MatDestroySubMatrices(nd,&submatsA);CHKERRQ(ierr); 196 197 ierr = PetscFree(is1);CHKERRQ(ierr); 198 ierr = PetscFree(is2);CHKERRQ(ierr); 199 ierr = PetscFree(idx);CHKERRQ(ierr); 200 ierr = PetscFree(rows);CHKERRQ(ierr); 201 ierr = PetscFree(cols);CHKERRQ(ierr); 202 ierr = PetscFree(vals);CHKERRQ(ierr); 203 ierr = MatDestroy(&A);CHKERRQ(ierr); 204 ierr = MatDestroy(&sA);CHKERRQ(ierr); 205 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 206 ierr = PetscFinalize(); 207 return ierr; 208 } 209 210 /*TEST 211 212 test: 213 args: -ov 2 214 215 TEST*/ 216