Actual source code: ex54.c
2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for parallel MatBAIJ format.\n";
4: #include petscmat.h
8: int main(int argc,char **args)
9: {
10: Mat A,B,*submatA,*submatB;
11: PetscInt bs=1,m=11,ov=1,i,j,k,*rows,*cols,nd=5,*idx,rstart,rend,sz,mm,nn,M,N,Mbs;
13: PetscMPIInt size,rank;
14: PetscScalar *vals,rval;
15: IS *is1,*is2;
16: PetscRandom rdm;
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,&rdm);
35: PetscRandomSetFromOptions(rdm);
37: MatGetOwnershipRange(A,&rstart,&rend);
38: MatGetSize(A,&M,&N);
39: Mbs = M/bs;
41: PetscMalloc(bs*sizeof(PetscInt),&rows);
42: PetscMalloc(bs*sizeof(PetscInt),&cols);
43: PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
44: PetscMalloc(M*sizeof(PetscScalar),&idx);
46: /* Now set blocks of values */
47: for (i=0; i<40*bs; i++) {
48: PetscRandomGetValue(rdm,&rval);
49: cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
50: PetscRandomGetValue(rdm,&rval);
51: rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m);
52: for (j=1; j<bs; j++) {
53: rows[j] = rows[j-1]+1;
54: cols[j] = cols[j-1]+1;
55: }
57: for (j=0; j<bs*bs; j++) {
58: PetscRandomGetValue(rdm,&rval);
59: vals[j] = rval;
60: }
61: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
62: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
63: }
65: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
67: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
68: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
70: /* Test MatIncreaseOverlap() */
71: PetscMalloc(nd*sizeof(IS **),&is1);
72: PetscMalloc(nd*sizeof(IS **),&is2);
74:
75: for (i=0; i<nd; i++) {
76: PetscRandomGetValue(rdm,&rval);
77: sz = (int)(PetscRealPart(rval)*m);
78: for (j=0; j<sz; j++) {
79: PetscRandomGetValue(rdm,&rval);
80: idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
81: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
82: }
83: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is1+i);
84: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,is2+i);
85: }
86: MatIncreaseOverlap(A,nd,is1,ov);
87: MatIncreaseOverlap(B,nd,is2,ov);
89: for (i=0; i<nd; ++i) {
90: ISEqual(is1[i],is2[i],&flg);
92: if (!flg) {
93: 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);
94: }
95: }
97: for (i=0; i<nd; ++i) {
98: ISSort(is1[i]);
99: ISSort(is2[i]);
100: }
101:
102: MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
103: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
106: /* Test MatMult() */
107: for (i=0; i<nd; i++) {
108: MatGetSize(submatA[i],&mm,&nn);
109: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
110: VecDuplicate(xx,&s1);
111: VecDuplicate(xx,&s2);
112: for (j=0; j<3; j++) {
113: VecSetRandom(xx,rdm);
114: MatMult(submatA[i],xx,s1);
115: MatMult(submatB[i],xx,s2);
116: VecNorm(s1,NORM_2,&s1norm);
117: VecNorm(s2,NORM_2,&s2norm);
118: rnorm = s2norm-s1norm;
119: if (rnorm<-tol || rnorm>tol) {
120: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
121: }
122: }
123: VecDestroy(xx);
124: VecDestroy(s1);
125: VecDestroy(s2);
126: }
128: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
129:
130: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
131: MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
133: /* Test MatMult() */
134: for (i=0; i<nd; i++) {
135: MatGetSize(submatA[i],&mm,&nn);
136: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
137: VecDuplicate(xx,&s1);
138: VecDuplicate(xx,&s2);
139: for (j=0; j<3; j++) {
140: VecSetRandom(xx,rdm);
141: MatMult(submatA[i],xx,s1);
142: MatMult(submatB[i],xx,s2);
143: VecNorm(s1,NORM_2,&s1norm);
144: VecNorm(s2,NORM_2,&s2norm);
145: rnorm = s2norm-s1norm;
146: if (rnorm<-tol || rnorm>tol) {
147: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
148: }
149: }
150: VecDestroy(xx);
151: VecDestroy(s1);
152: VecDestroy(s2);
153: }
154:
155: /* Free allocated memory */
156: for (i=0; i<nd; ++i) {
157: ISDestroy(is1[i]);
158: ISDestroy(is2[i]);
159: MatDestroy(submatA[i]);
160: MatDestroy(submatB[i]);
161: }
162: PetscFree(is1);
163: PetscFree(is2);
164: PetscFree(idx);
165: PetscFree(rows);
166: PetscFree(cols);
167: PetscFree(vals);
168: MatDestroy(A);
169: MatDestroy(B);
170: PetscFree(submatA);
171: PetscFree(submatB);
172: PetscRandomDestroy(rdm);
173: PetscFinalize();
174: return 0;
175: }