Actual source code: ex51.c
petsc-3.10.5 2019-03-28
2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for MatBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A,B,*submatA,*submatB;
9: PetscInt bs=1,m=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn,lsize;
11: PetscScalar *vals,rval;
12: IS *is1,*is2;
13: PetscRandom rdm;
14: Vec xx,s1,s2;
15: PetscReal s1norm,s2norm,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
16: PetscBool flg;
18: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
19: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
20: PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
21: PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
22: PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
23: M = m*bs;
25: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
26: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
27: MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B);
28: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
29: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
30: PetscRandomSetFromOptions(rdm);
32: PetscMalloc1(bs,&rows);
33: PetscMalloc1(bs,&cols);
34: PetscMalloc1(bs*bs,&vals);
35: PetscMalloc1(M,&idx);
37: /* Now set blocks of values */
38: for (i=0; i<20*bs; i++) {
39: PetscRandomGetValue(rdm,&rval);
40: cols[0] = bs*(int)(PetscRealPart(rval)*m);
41: PetscRandomGetValue(rdm,&rval);
42: rows[0] = bs*(int)(PetscRealPart(rval)*m);
43: for (j=1; j<bs; j++) {
44: rows[j] = rows[j-1]+1;
45: cols[j] = cols[j-1]+1;
46: }
48: for (j=0; j<bs*bs; j++) {
49: PetscRandomGetValue(rdm,&rval);
50: vals[j] = rval;
51: }
52: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
53: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
54: }
56: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
57: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
58: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
59: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
61: /* Test MatIncreaseOverlap() */
62: PetscMalloc1(nd,&is1);
63: PetscMalloc1(nd,&is2);
66: for (i=0; i<nd; i++) {
67: PetscRandomGetValue(rdm,&rval);
68: lsize = (int)(PetscRealPart(rval)*m);
69: for (j=0; j<lsize; j++) {
70: PetscRandomGetValue(rdm,&rval);
71: idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
72: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
73: }
74: ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,PETSC_COPY_VALUES,is1+i);
75: ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,PETSC_COPY_VALUES,is2+i);
76: }
77: MatIncreaseOverlap(A,nd,is1,ov);
78: MatIncreaseOverlap(B,nd,is2,ov);
80: for (i=0; i<nd; ++i) {
81: ISEqual(is1[i],is2[i],&flg);
82: if (!flg) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"i=%D, flg =%d\n",i,(int)flg);
83: }
85: for (i=0; i<nd; ++i) {
86: ISSort(is1[i]);
87: ISSort(is2[i]);
88: }
90: MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
91: MatCreateSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
93: /* Test MatMult() */
94: for (i=0; i<nd; i++) {
95: MatGetSize(submatA[i],&mm,&nn);
96: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
97: VecDuplicate(xx,&s1);
98: VecDuplicate(xx,&s2);
99: for (j=0; j<3; j++) {
100: VecSetRandom(xx,rdm);
101: MatMult(submatA[i],xx,s1);
102: MatMult(submatB[i],xx,s2);
103: VecNorm(s1,NORM_2,&s1norm);
104: VecNorm(s2,NORM_2,&s2norm);
105: rnorm = s2norm-s1norm;
106: if (rnorm<-tol || rnorm>tol) {
107: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
108: }
109: }
110: VecDestroy(&xx);
111: VecDestroy(&s1);
112: VecDestroy(&s2);
113: }
114: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
115: MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
116: MatCreateSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
118: /* Test MatMult() */
119: for (i=0; i<nd; i++) {
120: MatGetSize(submatA[i],&mm,&nn);
121: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
122: VecDuplicate(xx,&s1);
123: VecDuplicate(xx,&s2);
124: for (j=0; j<3; j++) {
125: VecSetRandom(xx,rdm);
126: MatMult(submatA[i],xx,s1);
127: MatMult(submatB[i],xx,s2);
128: VecNorm(s1,NORM_2,&s1norm);
129: VecNorm(s2,NORM_2,&s2norm);
130: rnorm = s2norm-s1norm;
131: if (rnorm<-tol || rnorm>tol) {
132: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
133: }
134: }
135: VecDestroy(&xx);
136: VecDestroy(&s1);
137: VecDestroy(&s2);
138: }
140: /* Free allocated memory */
141: for (i=0; i<nd; ++i) {
142: ISDestroy(&is1[i]);
143: ISDestroy(&is2[i]);
144: }
145: MatDestroySubMatrices(nd,&submatA);
146: MatDestroySubMatrices(nd,&submatB);
147: PetscFree(is1);
148: PetscFree(is2);
149: PetscFree(idx);
150: PetscFree(rows);
151: PetscFree(cols);
152: PetscFree(vals);
153: MatDestroy(&A);
154: MatDestroy(&B);
155: PetscRandomDestroy(&rdm);
156: PetscFinalize();
157: return ierr;
158: }
161: /*TEST
163: test:
164: args: -mat_block_size {{1 2 5 7 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}}
166: TEST*/