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: PetscBool 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,PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL,&A); */
31: MatCreate(PETSC_COMM_WORLD,&A);
32: MatSetSizes(A,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
33: MatSetType(A,MATBAIJ);
34: MatSetBlockSize(A,bs);
35: MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,PETSC_NULL);
36: MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL);
37: MatSetFromOptions(A);
39: MatCreateMPIAIJ(PETSC_COMM_WORLD,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,
40: PETSC_DEFAULT,PETSC_NULL,PETSC_DEFAULT,PETSC_NULL,&B);
42: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
43: PetscRandomSetFromOptions(rdm);
45: MatGetOwnershipRange(A,&rstart,&rend);
46: MatGetSize(A,&M,&N);
47: Mbs = M/bs;
49: PetscMalloc(bs*sizeof(PetscInt),&rows);
50: PetscMalloc(bs*sizeof(PetscInt),&cols);
51: PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
52: PetscMalloc(M*sizeof(PetscScalar),&idx);
54: /* Now set blocks of values */
55: for (i=0; i<40*bs; i++) {
56: PetscRandomGetValue(rdm,&rval);
57: cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
58: PetscRandomGetValue(rdm,&rval);
59: rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m);
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(rdm,&rval);
67: vals[j] = rval;
68: }
69: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
70: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
71: }
73: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
74: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
75: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
78: /* Test MatIncreaseOverlap() */
79: PetscMalloc(nd*sizeof(IS **),&is1);
80: PetscMalloc(nd*sizeof(IS **),&is2);
82:
83: for (i=0; i<nd; i++) {
84: PetscRandomGetValue(rdm,&rval);
85: sz = (int)(PetscRealPart(rval)*m);
86: for (j=0; j<sz; j++) {
87: PetscRandomGetValue(rdm,&rval);
88: idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
89: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
90: }
91: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
92: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
93: }
94: MatIncreaseOverlap(A,nd,is1,ov);
95: MatIncreaseOverlap(B,nd,is2,ov);
97: for (i=0; i<nd; ++i) {
98: ISEqual(is1[i],is2[i],&flg);
100: if (!flg) {
101: 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);
102: }
103: }
105: for (i=0; i<nd; ++i) {
106: ISSort(is1[i]);
107: ISSort(is2[i]);
108: }
109:
110: MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
111: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
114: /* Test MatMult() */
115: for (i=0; i<nd; i++) {
116: MatGetSize(submatA[i],&mm,&nn);
117: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
118: VecDuplicate(xx,&s1);
119: VecDuplicate(xx,&s2);
120: for (j=0; j<3; j++) {
121: VecSetRandom(xx,rdm);
122: MatMult(submatA[i],xx,s1);
123: MatMult(submatB[i],xx,s2);
124: VecNorm(s1,NORM_2,&s1norm);
125: VecNorm(s2,NORM_2,&s2norm);
126: rnorm = s2norm-s1norm;
127: if (rnorm<-tol || rnorm>tol) {
128: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
129: }
130: }
131: VecDestroy(&xx);
132: VecDestroy(&s1);
133: VecDestroy(&s2);
134: }
136: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
137:
138: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
139: MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
141: /* Test MatMult() */
142: for (i=0; i<nd; i++) {
143: MatGetSize(submatA[i],&mm,&nn);
144: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
145: VecDuplicate(xx,&s1);
146: VecDuplicate(xx,&s2);
147: for (j=0; j<3; j++) {
148: VecSetRandom(xx,rdm);
149: MatMult(submatA[i],xx,s1);
150: MatMult(submatB[i],xx,s2);
151: VecNorm(s1,NORM_2,&s1norm);
152: VecNorm(s2,NORM_2,&s2norm);
153: rnorm = s2norm-s1norm;
154: if (rnorm<-tol || rnorm>tol) {
155: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
156: }
157: }
158: VecDestroy(&xx);
159: VecDestroy(&s1);
160: VecDestroy(&s2);
161: }
162:
163: /* Free allocated memory */
164: for (i=0; i<nd; ++i) {
165: ISDestroy(&is1[i]);
166: ISDestroy(&is2[i]);
167: MatDestroy(&submatA[i]);
168: MatDestroy(&submatB[i]);
169: }
170: PetscFree(is1);
171: PetscFree(is2);
172: PetscFree(idx);
173: PetscFree(rows);
174: PetscFree(cols);
175: PetscFree(vals);
176: MatDestroy(&A);
177: MatDestroy(&B);
178: PetscFree(submatA);
179: PetscFree(submatB);
180: PetscRandomDestroy(&rdm);
181: PetscFinalize();
182: return 0;
183: }