Actual source code: ex91.c
petsc-3.7.3 2016-08-01
2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n";
4: #include <petscmat.h>
8: int main(int argc,char **args)
9: {
10: Mat A,Atrans,sA,*submatA,*submatsA;
11: PetscInt bs=1,m=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn;
13: PetscMPIInt size;
14: PetscScalar *vals,rval,one=1.0;
15: IS *is1,*is2;
16: PetscRandom rand;
17: Vec xx,s1,s2;
18: PetscReal s1norm,s2norm,rnorm,tol = 1.e-10;
19: PetscBool flg;
21: PetscInitialize(&argc,&args,(char*)0,help);
24: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
25: PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
26: PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
27: PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
29: /* create a SeqBAIJ matrix A */
30: M = m*bs;
31: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
32: PetscRandomCreate(PETSC_COMM_SELF,&rand);
33: PetscRandomSetFromOptions(rand);
35: PetscMalloc1(bs,&rows);
36: PetscMalloc1(bs,&cols);
37: PetscMalloc1(bs*bs,&vals);
38: PetscMalloc1(M,&idx);
40: /* Now set blocks of random values */
41: /* first, set diagonal blocks as zero */
42: for (j=0; j<bs*bs; j++) vals[j] = 0.0;
43: for (i=0; i<m; i++) {
44: cols[0] = i*bs; rows[0] = i*bs;
45: for (j=1; j<bs; j++) {
46: rows[j] = rows[j-1]+1;
47: cols[j] = cols[j-1]+1;
48: }
49: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
50: }
51: /* second, add random blocks */
52: for (i=0; i<20*bs; i++) {
53: PetscRandomGetValue(rand,&rval);
54: cols[0] = bs*(int)(PetscRealPart(rval)*m);
55: PetscRandomGetValue(rand,&rval);
56: rows[0] = bs*(int)(PetscRealPart(rval)*m);
57: for (j=1; j<bs; j++) {
58: rows[j] = rows[j-1]+1;
59: cols[j] = cols[j-1]+1;
60: }
62: for (j=0; j<bs*bs; j++) {
63: PetscRandomGetValue(rand,&rval);
64: vals[j] = rval;
65: }
66: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
67: }
69: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
70: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
72: /* make A a symmetric matrix: A <- A^T + A */
73: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
74: MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
75: MatDestroy(&Atrans);
76: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
77: MatEqual(A, Atrans, &flg);
78: if (!flg) SETERRQ(PETSC_COMM_SELF,1,"A+A^T is non-symmetric");
79: MatDestroy(&Atrans);
81: /* create a SeqSBAIJ matrix sA (= A) */
82: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
84: /* Test sA==A through MatMult() */
85: for (i=0; i<nd; i++) {
86: MatGetSize(A,&mm,&nn);
87: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
88: VecDuplicate(xx,&s1);
89: VecDuplicate(xx,&s2);
90: for (j=0; j<3; j++) {
91: VecSetRandom(xx,rand);
92: MatMult(A,xx,s1);
93: MatMult(sA,xx,s2);
94: VecNorm(s1,NORM_2,&s1norm);
95: VecNorm(s2,NORM_2,&s2norm);
96: rnorm = s2norm-s1norm;
97: if (rnorm<-tol || rnorm>tol) {
98: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
99: }
100: }
101: VecDestroy(&xx);
102: VecDestroy(&s1);
103: VecDestroy(&s2);
104: }
106: /* Test MatIncreaseOverlap() */
107: PetscMalloc1(nd,&is1);
108: PetscMalloc1(nd,&is2);
111: for (i=0; i<nd; i++) {
112: PetscRandomGetValue(rand,&rval);
113: size = (int)(PetscRealPart(rval)*m);
114: for (j=0; j<size; j++) {
115: PetscRandomGetValue(rand,&rval);
116: idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
117: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
118: }
119: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is1+i);
120: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is2+i);
121: }
122: /* for debugging */
123: /*
124: MatView(A,PETSC_VIEWER_STDOUT_SELF);
125: MatView(sA,PETSC_VIEWER_STDOUT_SELF);
126: */
128: MatIncreaseOverlap(A,nd,is1,ov);
129: MatIncreaseOverlap(sA,nd,is2,ov);
131: for (i=0; i<nd; ++i) {
132: ISSort(is1[i]);
133: ISSort(is2[i]);
134: }
136: for (i=0; i<nd; ++i) {
137: ISEqual(is1[i],is2[i],&flg);
138: if (!flg) SETERRQ1(PETSC_COMM_SELF,1,"i=%d, is1 != is2",i);
139: }
141: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
142: MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);
144: /* Test MatMult() */
145: for (i=0; i<nd; i++) {
146: MatGetSize(submatA[i],&mm,&nn);
147: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
148: VecDuplicate(xx,&s1);
149: VecDuplicate(xx,&s2);
150: for (j=0; j<3; j++) {
151: VecSetRandom(xx,rand);
152: MatMult(submatA[i],xx,s1);
153: MatMult(submatsA[i],xx,s2);
154: VecNorm(s1,NORM_2,&s1norm);
155: VecNorm(s2,NORM_2,&s2norm);
156: rnorm = s2norm-s1norm;
157: if (rnorm<-tol || rnorm>tol) {
158: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
159: }
160: }
161: VecDestroy(&xx);
162: VecDestroy(&s1);
163: VecDestroy(&s2);
164: }
166: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
167: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
168: MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
170: /* Test MatMult() */
171: for (i=0; i<nd; i++) {
172: MatGetSize(submatA[i],&mm,&nn);
173: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
174: VecDuplicate(xx,&s1);
175: VecDuplicate(xx,&s2);
176: for (j=0; j<3; j++) {
177: VecSetRandom(xx,rand);
178: MatMult(submatA[i],xx,s1);
179: MatMult(submatsA[i],xx,s2);
180: VecNorm(s1,NORM_2,&s1norm);
181: VecNorm(s2,NORM_2,&s2norm);
182: rnorm = s2norm-s1norm;
183: if (rnorm<-tol || rnorm>tol) {
184: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
185: }
186: }
187: VecDestroy(&xx);
188: VecDestroy(&s1);
189: VecDestroy(&s2);
190: }
192: /* Free allocated memory */
193: for (i=0; i<nd; ++i) {
194: ISDestroy(&is1[i]);
195: ISDestroy(&is2[i]);
197: MatDestroy(&submatA[i]);
198: MatDestroy(&submatsA[i]);
200: }
202: PetscFree(submatA);
203: PetscFree(submatsA);
205: PetscFree(is1);
206: PetscFree(is2);
207: PetscFree(idx);
208: PetscFree(rows);
209: PetscFree(cols);
210: PetscFree(vals);
211: MatDestroy(&A);
212: MatDestroy(&sA);
213: PetscRandomDestroy(&rand);
214: PetscFinalize();
215: return 0;
216: }