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