Actual source code: ex54.c
2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for parallel AIJ and BAIJ formats.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A,B,*submatA,*submatB;
9: PetscInt bs=1,m=11,ov=1,i,j,k,*rows,*cols,nd=5,*idx,rstart,rend,sz,mm,nn,M,N,Mbs;
10: PetscMPIInt size,rank;
11: PetscScalar *vals,rval;
12: IS *is1,*is2;
13: PetscRandom rdm;
14: Vec xx,s1,s2;
15: PetscReal s1norm,s2norm,rnorm,tol = 100*PETSC_SMALL;
16: PetscBool flg,test_nd0=PETSC_FALSE;
18: PetscInitialize(&argc,&args,(char*)0,help);
19: MPI_Comm_size(PETSC_COMM_WORLD,&size);
20: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
22: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
23: PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
24: PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
25: PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
26: PetscOptionsGetBool(NULL,NULL,"-test_nd0",&test_nd0,NULL);
28: /* Create a AIJ matrix A */
29: MatCreate(PETSC_COMM_WORLD,&A);
30: MatSetSizes(A,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
31: MatSetType(A,MATAIJ);
32: MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,NULL);
33: MatMPIAIJSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
34: MatSetFromOptions(A);
35: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
37: /* Create a BAIJ matrix B */
38: MatCreate(PETSC_COMM_WORLD,&B);
39: MatSetSizes(B,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
40: MatSetType(B,MATBAIJ);
41: MatSeqBAIJSetPreallocation(B,bs,PETSC_DEFAULT,NULL);
42: MatMPIBAIJSetPreallocation(B,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
43: MatSetFromOptions(B);
44: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
46: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
47: PetscRandomSetFromOptions(rdm);
49: MatGetOwnershipRange(A,&rstart,&rend);
50: MatGetSize(A,&M,&N);
51: Mbs = M/bs;
53: PetscMalloc1(bs,&rows);
54: PetscMalloc1(bs,&cols);
55: PetscMalloc1(bs*bs,&vals);
56: PetscMalloc1(M,&idx);
58: /* Now set blocks of values */
59: for (i=0; i<40*bs; i++) {
60: PetscRandomGetValue(rdm,&rval);
61: cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
62: PetscRandomGetValue(rdm,&rval);
63: rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m);
64: for (j=1; j<bs; j++) {
65: rows[j] = rows[j-1]+1;
66: cols[j] = cols[j-1]+1;
67: }
69: for (j=0; j<bs*bs; j++) {
70: PetscRandomGetValue(rdm,&rval);
71: vals[j] = rval;
72: }
73: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
74: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
75: }
77: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
79: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
82: /* Test MatIncreaseOverlap() */
83: PetscMalloc1(nd,&is1);
84: PetscMalloc1(nd,&is2);
86: if (rank == 0 && test_nd0) nd = 0; /* test case */
88: for (i=0; i<nd; i++) {
89: PetscRandomGetValue(rdm,&rval);
90: sz = (int)(PetscRealPart(rval)*m);
91: for (j=0; j<sz; j++) {
92: PetscRandomGetValue(rdm,&rval);
93: idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
94: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
95: }
96: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
97: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
98: }
99: MatIncreaseOverlap(A,nd,is1,ov);
100: MatIncreaseOverlap(B,nd,is2,ov);
102: for (i=0; i<nd; ++i) {
103: ISEqual(is1[i],is2[i],&flg);
105: if (!flg) {
106: PetscPrintf(PETSC_COMM_SELF,"i=%" PetscInt_FMT ", flg=%d :bs=%" PetscInt_FMT " m=%" PetscInt_FMT " ov=%" PetscInt_FMT " nd=%" PetscInt_FMT " np=%d\n",i,flg,bs,m,ov,nd,size);
107: }
108: }
110: for (i=0; i<nd; ++i) {
111: ISSort(is1[i]);
112: ISSort(is2[i]);
113: }
115: MatCreateSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
116: MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
118: /* Test MatMult() */
119: for (i=0; i<nd; i++) {
120: MatGetSize(submatA[i],&mm,&nn);
121: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
122: VecDuplicate(xx,&s1);
123: VecDuplicate(xx,&s2);
124: for (j=0; j<3; j++) {
125: VecSetRandom(xx,rdm);
126: MatMult(submatA[i],xx,s1);
127: MatMult(submatB[i],xx,s2);
128: VecNorm(s1,NORM_2,&s1norm);
129: VecNorm(s2,NORM_2,&s2norm);
130: rnorm = s2norm-s1norm;
131: if (rnorm<-tol || rnorm>tol) {
132: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,(double)s1norm,(double)s2norm);
133: }
134: }
135: VecDestroy(&xx);
136: VecDestroy(&s1);
137: VecDestroy(&s2);
138: }
140: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
141: MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
142: MatCreateSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
144: /* Test MatMult() */
145: for (i=0; i<nd; i++) {
146: MatGetSize(submatA[i],&mm,&nn);
147: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
148: VecDuplicate(xx,&s1);
149: VecDuplicate(xx,&s2);
150: for (j=0; j<3; j++) {
151: VecSetRandom(xx,rdm);
152: MatMult(submatA[i],xx,s1);
153: MatMult(submatB[i],xx,s2);
154: VecNorm(s1,NORM_2,&s1norm);
155: VecNorm(s2,NORM_2,&s2norm);
156: rnorm = s2norm-s1norm;
157: if (rnorm<-tol || rnorm>tol) {
158: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,(double)s1norm,(double)s2norm);
159: }
160: }
161: VecDestroy(&xx);
162: VecDestroy(&s1);
163: VecDestroy(&s2);
164: }
166: /* Free allocated memory */
167: for (i=0; i<nd; ++i) {
168: ISDestroy(&is1[i]);
169: ISDestroy(&is2[i]);
170: }
171: MatDestroySubMatrices(nd,&submatA);
172: MatDestroySubMatrices(nd,&submatB);
174: PetscFree(is1);
175: PetscFree(is2);
176: PetscFree(idx);
177: PetscFree(rows);
178: PetscFree(cols);
179: PetscFree(vals);
180: MatDestroy(&A);
181: MatDestroy(&B);
182: PetscRandomDestroy(&rdm);
183: PetscFinalize();
184: return 0;
185: }
187: /*TEST
189: test:
190: nsize: {{1 3}}
191: args: -mat_block_size {{1 3 4 6 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}} ; done
192: output_file: output/ex54.out
194: test:
195: suffix: 2
196: args: -nd 2 -test_nd0
197: output_file: output/ex54.out
199: test:
200: suffix: 3
201: nsize: 3
202: args: -nd 2 -test_nd0
203: output_file: output/ex54.out
205: TEST*/