Actual source code: ex91.c
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: PetscTruth flg;
21: PetscInitialize(&argc,&args,(char *)0,help);
22:
24: PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
25: PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);
26: PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
27: PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);
29: /* create a SeqBAIJ matrix A */
30: M = m*bs;
31: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,PETSC_NULL,&A);
32: PetscRandomCreate(PETSC_COMM_SELF,&rand);
33: PetscRandomSetFromOptions(rand);
35: PetscMalloc(bs*sizeof(PetscInt),&rows);
36: PetscMalloc(bs*sizeof(PetscInt),&cols);
37: PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
38: PetscMalloc(M*sizeof(PetscScalar),&idx);
39:
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, &Atrans);
74: MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
75: MatDestroy(Atrans);
76: MatTranspose(A, &Atrans);
77: MatEqual(A, Atrans, &flg);
78: if (!flg) {
79: SETERRQ(1,"A+A^T is non-symmetric");
80: }
81: MatDestroy(Atrans);
83: /* create a SeqSBAIJ matrix sA (= A) */
84: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
85:
86: /* Test sA==A through MatMult() */
87: for (i=0; i<nd; i++) {
88: MatGetSize(A,&mm,&nn);
89: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
90: VecDuplicate(xx,&s1);
91: VecDuplicate(xx,&s2);
92: for (j=0; j<3; j++) {
93: VecSetRandom(xx,rand);
94: MatMult(A,xx,s1);
95: MatMult(sA,xx,s2);
96: VecNorm(s1,NORM_2,&s1norm);
97: VecNorm(s2,NORM_2,&s2norm);
98: rnorm = s2norm-s1norm;
99: if (rnorm<-tol || rnorm>tol) {
100: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
101: }
102: }
103: VecDestroy(xx);
104: VecDestroy(s1);
105: VecDestroy(s2);
106: }
108: /* Test MatIncreaseOverlap() */
109: PetscMalloc(nd*sizeof(IS **),&is1);
110: PetscMalloc(nd*sizeof(IS **),&is2);
112:
113: for (i=0; i<nd; i++) {
114: PetscRandomGetValue(rand,&rval);
115: size = (int)(PetscRealPart(rval)*m);
116: for (j=0; j<size; j++) {
117: PetscRandomGetValue(rand,&rval);
118: idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
119: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
120: }
121: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,is1+i);
122: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,is2+i);
123: }
124: /* for debugging */
125: /*
126: MatView(A,PETSC_VIEWER_STDOUT_SELF);
127: MatView(sA,PETSC_VIEWER_STDOUT_SELF);
128: */
130: MatIncreaseOverlap(A,nd,is1,ov);
131: MatIncreaseOverlap(sA,nd,is2,ov);
133: for (i=0; i<nd; ++i) {
134: ISSort(is1[i]);
135: ISSort(is2[i]);
136: }
138: for (i=0; i<nd; ++i) {
139: ISEqual(is1[i],is2[i],&flg);
140: if (!flg){
141: /* ISView(is1[i],PETSC_VIEWER_STDOUT_SELF);
142: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF); */
143: SETERRQ1(1,"i=%d, is1 != is2",i);
144: }
145: }
146:
147: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
148: MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);
150: /* Test MatMult() */
151: for (i=0; i<nd; i++) {
152: MatGetSize(submatA[i],&mm,&nn);
153: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
154: VecDuplicate(xx,&s1);
155: VecDuplicate(xx,&s2);
156: for (j=0; j<3; j++) {
157: VecSetRandom(xx,rand);
158: MatMult(submatA[i],xx,s1);
159: MatMult(submatsA[i],xx,s2);
160: VecNorm(s1,NORM_2,&s1norm);
161: VecNorm(s2,NORM_2,&s2norm);
162: rnorm = s2norm-s1norm;
163: if (rnorm<-tol || rnorm>tol) {
164: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
165: }
166: }
167: VecDestroy(xx);
168: VecDestroy(s1);
169: VecDestroy(s2);
170: }
172: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
173: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
174: MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
175:
176: /* Test MatMult() */
177: for (i=0; i<nd; i++) {
178: MatGetSize(submatA[i],&mm,&nn);
179: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
180: VecDuplicate(xx,&s1);
181: VecDuplicate(xx,&s2);
182: for (j=0; j<3; j++) {
183: VecSetRandom(xx,rand);
184: MatMult(submatA[i],xx,s1);
185: MatMult(submatsA[i],xx,s2);
186: VecNorm(s1,NORM_2,&s1norm);
187: VecNorm(s2,NORM_2,&s2norm);
188: rnorm = s2norm-s1norm;
189: if (rnorm<-tol || rnorm>tol) {
190: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
191: }
192: }
193: VecDestroy(xx);
194: VecDestroy(s1);
195: VecDestroy(s2);
196: }
197:
198: /* Free allocated memory */
199: for (i=0; i<nd; ++i) {
200: ISDestroy(is1[i]);
201: ISDestroy(is2[i]);
202:
203: MatDestroy(submatA[i]);
204: MatDestroy(submatsA[i]);
205:
206: }
207:
208: PetscFree(submatA);
209: PetscFree(submatsA);
210:
211: PetscFree(is1);
212: PetscFree(is2);
213: PetscFree(idx);
214: PetscFree(rows);
215: PetscFree(cols);
216: PetscFree(vals);
217: MatDestroy(A);
218: MatDestroy(sA);
219: PetscRandomDestroy(rand);
220: PetscFinalize();
221: return 0;
222: }