Actual source code: ex92.c
1: /*$Id: ex54.c,v 1.22 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for parallel MatSBAIJ format.\n";
5: #include petscmat.h
9: int main(int argc,char **args)
10: {
11: Mat A,Atrans,sA,*submatA,*submatsA;
12: int bs=1,mbs=11,ov=1,i,j,k,*rows,*cols,ierr,nd=5,*idx,size,
13: rank,rstart,rend,sz,mm,M,N,Mbs;
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;
20: int stages[2];
22: PetscInitialize(&argc,&args,(char *)0,help);
23: MPI_Comm_size(PETSC_COMM_WORLD,&size);
24: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
26: PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
27: PetscOptionsGetInt(PETSC_NULL,"-mat_size",&mbs,PETSC_NULL);
28: PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
29: PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);
31: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,mbs*bs,mbs*bs,PETSC_DECIDE,PETSC_DECIDE,
32: PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL,&A);
33: PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rand);
35: MatGetOwnershipRange(A,&rstart,&rend);
36: MatGetSize(A,&M,&N);
37: Mbs = M/bs;
39: PetscMalloc(bs*sizeof(int),&rows);
40: PetscMalloc(bs*sizeof(int),&cols);
41: PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
42: PetscMalloc(M*sizeof(PetscScalar),&idx);
44: /* Now set blocks of values */
45: for (j=0; j<bs*bs; j++) vals[j] = 0.0;
46: for (i=0; i<Mbs; i++){
47: cols[0] = i*bs; rows[0] = i*bs;
48: for (j=1; j<bs; j++) {
49: rows[j] = rows[j-1]+1;
50: cols[j] = cols[j-1]+1;
51: }
52: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
53: }
54: /* second, add random blocks */
55: for (i=0; i<20*bs; i++) {
56: PetscRandomGetValue(rand,&rval);
57: cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
58: PetscRandomGetValue(rand,&rval);
59: rows[0] = rstart + bs*(int)(PetscRealPart(rval)*mbs);
60: for (j=1; j<bs; j++) {
61: rows[j] = rows[j-1]+1;
62: cols[j] = cols[j-1]+1;
63: }
65: for (j=0; j<bs*bs; j++) {
66: PetscRandomGetValue(rand,&rval);
67: vals[j] = rval;
68: }
69: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
70: }
72: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
74:
75: /* make A a symmetric matrix: A <- A^T + A */
76: MatTranspose(A, &Atrans);
77: MatAXPY(&one,Atrans,A,DIFFERENT_NONZERO_PATTERN);
78: MatDestroy(Atrans);
79: MatTranspose(A, &Atrans);
80: MatEqual(A, Atrans, &flg);
81: if (!flg) {
82: SETERRQ(1,"A+A^T is non-symmetric");
83: }
84: MatDestroy(Atrans);
85: /* MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
87: /* create a SeqSBAIJ matrix sA (= A) */
88: MatConvert(A,MATMPISBAIJ,&sA);
89: /* MatView(sA,PETSC_VIEWER_STDOUT_WORLD); */
91: /* Test sA==A through MatMult() */
92: for (i=0; i<nd; i++) {
93: MatGetLocalSize(A, &mm,PETSC_NULL);
94: VecCreate(PETSC_COMM_WORLD,&xx);
95: VecSetSizes(xx,mm,PETSC_DECIDE);
96: VecSetFromOptions(xx);
97: VecDuplicate(xx,&s1);
98: VecDuplicate(xx,&s2);
99: for (j=0; j<3; j++) {
100: VecSetRandom(rand,xx);
101: MatMult(A,xx,s1);
102: MatMult(sA,xx,s2);
103: VecNorm(s1,NORM_2,&s1norm);
104: VecNorm(s2,NORM_2,&s2norm);
105: rnorm = s2norm-s1norm;
106: if (rnorm<-tol || rnorm>tol) {
107: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
108: }
109: }
110: VecDestroy(xx);
111: VecDestroy(s1);
112: VecDestroy(s2);
113: }
114:
115: /* Test MatIncreaseOverlap() */
116: PetscMalloc(nd*sizeof(IS **),&is1);
117: PetscMalloc(nd*sizeof(IS **),&is2);
119: for (i=0; i<nd; i++) {
120: PetscRandomGetValue(rand,&rval);
121: sz = (int)((0.5 + 0.2*PetscRealPart(rval))*mbs); /* 0.5*mbs < sz < 0.7*mbs */
122: sz /= (size*nd*10);
123: /*
124: if (!rank){
125: PetscPrintf(PETSC_COMM_SELF," [%d] IS sz[%d]: %d\n",rank,i,sz);
126: }
127: */
128: for (j=0; j<sz; j++) {
129: PetscRandomGetValue(rand,&rval);
130: idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
131: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
132: }
133: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is1+i);
134: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is2+i);
135: }
137: PetscLogStageRegister(&stages[0],"MatOv_SBAIJ");
138: PetscLogStageRegister(&stages[1],"MatOv_ BAIJ");
140: PetscLogStagePush(stages[0]);
141: MatIncreaseOverlap(sA,nd,is2,ov);
142: PetscLogStagePop();
144: PetscLogStagePush(stages[1]);
145: MatIncreaseOverlap(A,nd,is1,ov);
146: PetscLogStagePop();
148: for (i=0; i<nd; ++i) {
149: ISEqual(is1[i],is2[i],&flg);
150: if (!flg ){
151: int prid = size;
152: if (rank == prid){
153: ISSort(is1[i]);
154: ISView(is1[i],PETSC_VIEWER_STDOUT_SELF);
155: ISSort(is2[i]);
156: ISView(is2[i],PETSC_VIEWER_STDOUT_SELF);
157: }
158: SETERRQ1(1,"i=%d, is1 != is2",i);
159: }
160: }
161:
162: for (i=0; i<nd; ++i) {
163: ISSort(is1[i]);
164: ISSort(is2[i]);
165: }
166: MatGetSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);
167: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
169: /* Test MatMult() */
170: for (i=0; i<nd; i++) {
171: MatGetSize(submatA[i],&mm,PETSC_NULL);
172: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
173: VecDuplicate(xx,&s1);
174: VecDuplicate(xx,&s2);
175: for (j=0; j<3; j++) {
176: VecSetRandom(rand,xx);
177: MatMult(submatA[i],xx,s1);
178: MatMult(submatsA[i],xx,s2);
179: VecNorm(s1,NORM_2,&s1norm);
180: VecNorm(s2,NORM_2,&s2norm);
181: rnorm = s2norm-s1norm;
182: if (rnorm<-tol || rnorm>tol) {
183: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
184: }
185: }
186: VecDestroy(xx);
187: VecDestroy(s1);
188: VecDestroy(s2);
189: }
191: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
192: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
193: MatGetSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
195: /* Test MatMult() */
196: for (i=0; i<nd; i++) {
197: MatGetSize(submatA[i],&mm,PETSC_NULL);
198: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
199: VecDuplicate(xx,&s1);
200: VecDuplicate(xx,&s2);
201: for (j=0; j<3; j++) {
202: VecSetRandom(rand,xx);
203: MatMult(submatA[i],xx,s1);
204: MatMult(submatsA[i],xx,s2);
205: VecNorm(s1,NORM_2,&s1norm);
206: VecNorm(s2,NORM_2,&s2norm);
207: rnorm = s2norm-s1norm;
208: if (rnorm<-tol || rnorm>tol) {
209: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
210: }
211: }
212: VecDestroy(xx);
213: VecDestroy(s1);
214: VecDestroy(s2);
215: }
217: /* Free allocated memory */
218: for (i=0; i<nd; ++i) {
219: ISDestroy(is1[i]);
220: ISDestroy(is2[i]);
221: MatDestroy(submatA[i]);
222: MatDestroy(submatsA[i]);
223: }
224: PetscFree(submatA);
225: PetscFree(submatsA);
226: PetscFree(is1);
227: PetscFree(is2);
228: PetscFree(idx);
229: PetscFree(rows);
230: PetscFree(cols);
231: PetscFree(vals);
232: MatDestroy(A);
233: MatDestroy(sA);
234: PetscRandomDestroy(rand);
235: PetscFinalize();
236: return 0;
237: }