Actual source code: ex91.c
1: /*$Id: ex91.c,v 1.22 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n";
5: #include petscmat.h
9: int main(int argc,char **args)
10: {
11: Mat A,Atrans,sA,*submatA,*submatsA;
12: int bs=1,m=43,ov=1,i,j,k,*rows,*cols,ierr,M,nd=5,*idx,size,mm,nn;
13: PetscScalar *vals,rval,one=1.0;
14: IS *is1,*is2;
15: PetscRandom rand;
16: Vec xx,s1,s2;
17: PetscReal s1norm,s2norm,rnorm,tol = 1.e-10;
18: PetscTruth flg;
20: PetscInitialize(&argc,&args,(char *)0,help);
21:
23: PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
24: PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);
25: PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
26: PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);
28: /* create a SeqBAIJ matrix A */
29: M = m*bs;
30: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,PETSC_NULL,&A);
31: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rand);
33: PetscMalloc(bs*sizeof(int),&rows);
34: PetscMalloc(bs*sizeof(int),&cols);
35: PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
36: PetscMalloc(M*sizeof(PetscScalar),&idx);
37:
38: /* Now set blocks of random values */
39: /* first, set diagonal blocks as zero */
40: for (j=0; j<bs*bs; j++) vals[j] = 0.0;
41: for (i=0; i<m; i++){
42: cols[0] = i*bs; rows[0] = i*bs;
43: for (j=1; j<bs; j++) {
44: rows[j] = rows[j-1]+1;
45: cols[j] = cols[j-1]+1;
46: }
47: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
48: }
49: /* second, add random blocks */
50: for (i=0; i<20*bs; i++) {
51: PetscRandomGetValue(rand,&rval);
52: cols[0] = bs*(int)(PetscRealPart(rval)*m);
53: PetscRandomGetValue(rand,&rval);
54: rows[0] = bs*(int)(PetscRealPart(rval)*m);
55: for (j=1; j<bs; j++) {
56: rows[j] = rows[j-1]+1;
57: cols[j] = cols[j-1]+1;
58: }
60: for (j=0; j<bs*bs; j++) {
61: PetscRandomGetValue(rand,&rval);
62: vals[j] = rval;
63: }
64: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
65: }
67: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
70: /* make A a symmetric matrix: A <- A^T + A */
71: MatTranspose(A, &Atrans);
72: MatAXPY(&one,Atrans,A,DIFFERENT_NONZERO_PATTERN);
73: MatDestroy(Atrans);
74: MatTranspose(A, &Atrans);
75: MatEqual(A, Atrans, &flg);
76: if (!flg) {
77: SETERRQ(1,"A+A^T is non-symmetric");
78: }
79: MatDestroy(Atrans);
81: /* create a SeqSBAIJ matrix sA (= A) */
82: MatConvert(A,MATSEQSBAIJ,&sA);
83:
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(rand,xx);
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: PetscMalloc(nd*sizeof(IS **),&is1);
108: PetscMalloc(nd*sizeof(IS **),&is2);
110:
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,is1+i);
120: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,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){
139: /* ISView(is1[i],PETSC_VIEWER_STDOUT_SELF);
140: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF); */
141: SETERRQ1(1,"i=%d, is1 != is2",i);
142: }
143: }
144:
145: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
146: MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);
148: /* Test MatMult() */
149: for (i=0; i<nd; i++) {
150: MatGetSize(submatA[i],&mm,&nn);
151: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
152: VecDuplicate(xx,&s1);
153: VecDuplicate(xx,&s2);
154: for (j=0; j<3; j++) {
155: VecSetRandom(rand,xx);
156: MatMult(submatA[i],xx,s1);
157: MatMult(submatsA[i],xx,s2);
158: VecNorm(s1,NORM_2,&s1norm);
159: VecNorm(s2,NORM_2,&s2norm);
160: rnorm = s2norm-s1norm;
161: if (rnorm<-tol || rnorm>tol) {
162: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
163: }
164: }
165: VecDestroy(xx);
166: VecDestroy(s1);
167: VecDestroy(s2);
168: }
170: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
171: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
172: MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
173:
174: /* Test MatMult() */
175: for (i=0; i<nd; i++) {
176: MatGetSize(submatA[i],&mm,&nn);
177: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
178: VecDuplicate(xx,&s1);
179: VecDuplicate(xx,&s2);
180: for (j=0; j<3; j++) {
181: VecSetRandom(rand,xx);
182: MatMult(submatA[i],xx,s1);
183: MatMult(submatsA[i],xx,s2);
184: VecNorm(s1,NORM_2,&s1norm);
185: VecNorm(s2,NORM_2,&s2norm);
186: rnorm = s2norm-s1norm;
187: if (rnorm<-tol || rnorm>tol) {
188: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
189: }
190: }
191: VecDestroy(xx);
192: VecDestroy(s1);
193: VecDestroy(s2);
194: }
195:
196: /* Free allocated memory */
197: for (i=0; i<nd; ++i) {
198: ISDestroy(is1[i]);
199: ISDestroy(is2[i]);
200:
201: MatDestroy(submatA[i]);
202: MatDestroy(submatsA[i]);
203:
204: }
205:
206: PetscFree(submatA);
207: PetscFree(submatsA);
208:
209: PetscFree(is1);
210: PetscFree(is2);
211: PetscFree(idx);
212: PetscFree(rows);
213: PetscFree(cols);
214: PetscFree(vals);
215: MatDestroy(A);
216: MatDestroy(sA);
217: PetscRandomDestroy(rand);
218: PetscFinalize();
219: return 0;
220: }