Actual source code: ex91.c
petsc-3.8.4 2018-03-24
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 = 1.e-10;
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: PetscRandomCreate(PETSC_COMM_SELF,&rand);
29: PetscRandomSetFromOptions(rand);
31: PetscMalloc1(bs,&rows);
32: PetscMalloc1(bs,&cols);
33: PetscMalloc1(bs*bs,&vals);
34: PetscMalloc1(M,&idx);
36: /* Now set blocks of random values */
37: /* first, set diagonal blocks as zero */
38: for (j=0; j<bs*bs; j++) vals[j] = 0.0;
39: for (i=0; i<m; i++) {
40: cols[0] = i*bs; rows[0] = i*bs;
41: for (j=1; j<bs; j++) {
42: rows[j] = rows[j-1]+1;
43: cols[j] = cols[j-1]+1;
44: }
45: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
46: }
47: /* second, add random blocks */
48: for (i=0; i<20*bs; i++) {
49: PetscRandomGetValue(rand,&rval);
50: cols[0] = bs*(int)(PetscRealPart(rval)*m);
51: PetscRandomGetValue(rand,&rval);
52: rows[0] = bs*(int)(PetscRealPart(rval)*m);
53: for (j=1; j<bs; j++) {
54: rows[j] = rows[j-1]+1;
55: cols[j] = cols[j-1]+1;
56: }
58: for (j=0; j<bs*bs; j++) {
59: PetscRandomGetValue(rand,&rval);
60: vals[j] = rval;
61: }
62: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
63: }
65: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
68: /* make A a symmetric matrix: A <- A^T + A */
69: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
70: MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
71: MatDestroy(&Atrans);
72: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
73: MatEqual(A, Atrans, &flg);
74: if (!flg) SETERRQ(PETSC_COMM_SELF,1,"A+A^T is non-symmetric");
75: MatDestroy(&Atrans);
77: /* create a SeqSBAIJ matrix sA (= A) */
78: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
80: /* Test sA==A through MatMult() */
81: for (i=0; i<nd; i++) {
82: MatGetSize(A,&mm,&nn);
83: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
84: VecDuplicate(xx,&s1);
85: VecDuplicate(xx,&s2);
86: for (j=0; j<3; j++) {
87: VecSetRandom(xx,rand);
88: MatMult(A,xx,s1);
89: MatMult(sA,xx,s2);
90: VecNorm(s1,NORM_2,&s1norm);
91: VecNorm(s2,NORM_2,&s2norm);
92: rnorm = s2norm-s1norm;
93: if (rnorm<-tol || rnorm>tol) {
94: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
95: }
96: }
97: VecDestroy(&xx);
98: VecDestroy(&s1);
99: VecDestroy(&s2);
100: }
102: /* Test MatIncreaseOverlap() */
103: PetscMalloc1(nd,&is1);
104: PetscMalloc1(nd,&is2);
107: for (i=0; i<nd; i++) {
108: PetscRandomGetValue(rand,&rval);
109: size = (int)(PetscRealPart(rval)*m);
110: for (j=0; j<size; j++) {
111: PetscRandomGetValue(rand,&rval);
112: idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
113: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
114: }
115: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is1+i);
116: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is2+i);
117: }
118: /* for debugging */
119: /*
120: MatView(A,PETSC_VIEWER_STDOUT_SELF);
121: MatView(sA,PETSC_VIEWER_STDOUT_SELF);
122: */
124: MatIncreaseOverlap(A,nd,is1,ov);
125: MatIncreaseOverlap(sA,nd,is2,ov);
127: for (i=0; i<nd; ++i) {
128: ISSort(is1[i]);
129: ISSort(is2[i]);
130: }
132: for (i=0; i<nd; ++i) {
133: ISEqual(is1[i],is2[i],&flg);
134: if (!flg) SETERRQ1(PETSC_COMM_SELF,1,"i=%d, is1 != is2",i);
135: }
137: MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
138: MatCreateSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);
140: /* Test MatMult() */
141: for (i=0; i<nd; i++) {
142: MatGetSize(submatA[i],&mm,&nn);
143: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
144: VecDuplicate(xx,&s1);
145: VecDuplicate(xx,&s2);
146: for (j=0; j<3; j++) {
147: VecSetRandom(xx,rand);
148: MatMult(submatA[i],xx,s1);
149: MatMult(submatsA[i],xx,s2);
150: VecNorm(s1,NORM_2,&s1norm);
151: VecNorm(s2,NORM_2,&s2norm);
152: rnorm = s2norm-s1norm;
153: if (rnorm<-tol || rnorm>tol) {
154: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
155: }
156: }
157: VecDestroy(&xx);
158: VecDestroy(&s1);
159: VecDestroy(&s2);
160: }
162: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
163: MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
164: MatCreateSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
166: /* Test MatMult() */
167: for (i=0; i<nd; i++) {
168: MatGetSize(submatA[i],&mm,&nn);
169: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
170: VecDuplicate(xx,&s1);
171: VecDuplicate(xx,&s2);
172: for (j=0; j<3; j++) {
173: VecSetRandom(xx,rand);
174: MatMult(submatA[i],xx,s1);
175: MatMult(submatsA[i],xx,s2);
176: VecNorm(s1,NORM_2,&s1norm);
177: VecNorm(s2,NORM_2,&s2norm);
178: rnorm = s2norm-s1norm;
179: if (rnorm<-tol || rnorm>tol) {
180: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
181: }
182: }
183: VecDestroy(&xx);
184: VecDestroy(&s1);
185: VecDestroy(&s2);
186: }
188: /* Free allocated memory */
189: for (i=0; i<nd; ++i) {
190: ISDestroy(&is1[i]);
191: ISDestroy(&is2[i]);
192: }
193: MatDestroySubMatrices(nd,&submatA);
194: MatDestroySubMatrices(nd,&submatsA);
196: PetscFree(is1);
197: PetscFree(is2);
198: PetscFree(idx);
199: PetscFree(rows);
200: PetscFree(cols);
201: PetscFree(vals);
202: MatDestroy(&A);
203: MatDestroy(&sA);
204: PetscRandomDestroy(&rand);
205: PetscFinalize();
206: return ierr;
207: }