Actual source code: ex54.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 MatBAIJ format.\n";
5: #include petscmat.h
9: int main(int argc,char **args)
10: {
11: Mat A,B,*submatA,*submatB;
12: int bs=1,m=11,ov=1,i,j,k,*rows,*cols,ierr,nd=5,*idx,size;
13: int rank,rstart,rend,sz,mm,nn,M,N,Mbs;
14: PetscScalar *vals,rval;
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: MPI_Comm_size(PETSC_COMM_WORLD,&size);
23: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
25: PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
26: PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);
27: PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
28: PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);
30: MatCreateMPIBAIJ(PETSC_COMM_WORLD,bs,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,
31: PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL,&A);
32: MatCreateMPIAIJ(PETSC_COMM_WORLD,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,
33: PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL,&B);
34: PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rand);
36: MatGetOwnershipRange(A,&rstart,&rend);
37: MatGetSize(A,&M,&N);
38: Mbs = M/bs;
40: PetscMalloc(bs*sizeof(int),&rows);
41: PetscMalloc(bs*sizeof(int),&cols);
42: PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
43: PetscMalloc(M*sizeof(PetscScalar),&idx);
45: /* Now set blocks of values */
46: for (i=0; i<40*bs; i++) {
47: PetscRandomGetValue(rand,&rval);
48: cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
49: PetscRandomGetValue(rand,&rval);
50: rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m);
51: for (j=1; j<bs; j++) {
52: rows[j] = rows[j-1]+1;
53: cols[j] = cols[j-1]+1;
54: }
56: for (j=0; j<bs*bs; j++) {
57: PetscRandomGetValue(rand,&rval);
58: vals[j] = rval;
59: }
60: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
61: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
62: }
64: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
65: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
66: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
69: /* Test MatIncreaseOverlap() */
70: PetscMalloc(nd*sizeof(IS **),&is1);
71: PetscMalloc(nd*sizeof(IS **),&is2);
73:
74: for (i=0; i<nd; i++) {
75: PetscRandomGetValue(rand,&rval);
76: sz = (int)(PetscRealPart(rval)*m);
77: for (j=0; j<sz; j++) {
78: PetscRandomGetValue(rand,&rval);
79: idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
80: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
81: }
82: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is1+i);
83: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is2+i);
84: }
85: MatIncreaseOverlap(A,nd,is1,ov);
86: MatIncreaseOverlap(B,nd,is2,ov);
88: for (i=0; i<nd; ++i) {
89: ISEqual(is1[i],is2[i],&flg);
91: if (!flg) {
92: PetscPrintf(PETSC_COMM_SELF,"i=%d, flg=%d :bs=%d m=%d ov=%d nd=%d np=%d\n",i,flg,bs,m,ov,nd,size);
93: }
94: }
96: for (i=0; i<nd; ++i) {
97: ISSort(is1[i]);
98: ISSort(is2[i]);
99: }
100:
101: MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
102: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
105: /* Test MatMult() */
106: for (i=0; i<nd; i++) {
107: MatGetSize(submatA[i],&mm,&nn);
108: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
109: VecDuplicate(xx,&s1);
110: VecDuplicate(xx,&s2);
111: for (j=0; j<3; j++) {
112: VecSetRandom(rand,xx);
113: MatMult(submatA[i],xx,s1);
114: MatMult(submatB[i],xx,s2);
115: VecNorm(s1,NORM_2,&s1norm);
116: VecNorm(s2,NORM_2,&s2norm);
117: rnorm = s2norm-s1norm;
118: if (rnorm<-tol || rnorm>tol) {
119: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
120: }
121: }
122: VecDestroy(xx);
123: VecDestroy(s1);
124: VecDestroy(s2);
125: }
127: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
128:
129: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
130: MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
132: /* Test MatMult() */
133: for (i=0; i<nd; i++) {
134: MatGetSize(submatA[i],&mm,&nn);
135: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
136: VecDuplicate(xx,&s1);
137: VecDuplicate(xx,&s2);
138: for (j=0; j<3; j++) {
139: VecSetRandom(rand,xx);
140: MatMult(submatA[i],xx,s1);
141: MatMult(submatB[i],xx,s2);
142: VecNorm(s1,NORM_2,&s1norm);
143: VecNorm(s2,NORM_2,&s2norm);
144: rnorm = s2norm-s1norm;
145: if (rnorm<-tol || rnorm>tol) {
146: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
147: }
148: }
149: VecDestroy(xx);
150: VecDestroy(s1);
151: VecDestroy(s2);
152: }
153:
154: /* Free allocated memory */
155: for (i=0; i<nd; ++i) {
156: ISDestroy(is1[i]);
157: ISDestroy(is2[i]);
158: MatDestroy(submatA[i]);
159: MatDestroy(submatB[i]);
160: }
161: PetscFree(is1);
162: PetscFree(is2);
163: PetscFree(idx);
164: PetscFree(rows);
165: PetscFree(cols);
166: PetscFree(vals);
167: MatDestroy(A);
168: MatDestroy(B);
169: PetscFree(submatA);
170: PetscFree(submatB);
171: PetscRandomDestroy(rand);
172: PetscFinalize();
173: return 0;
174: }