Actual source code: ex54.c
petsc-3.8.4 2018-03-24
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;
11: PetscMPIInt size,rank;
12: PetscScalar *vals,rval;
13: IS *is1,*is2;
14: PetscRandom rdm;
15: Vec xx,s1,s2;
16: PetscReal s1norm,s2norm,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
17: PetscBool flg,test_nd0=PETSC_FALSE;
19: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20: MPI_Comm_size(PETSC_COMM_WORLD,&size);
21: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
23: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
24: PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
25: PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
26: PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
27: PetscOptionsGetBool(NULL,NULL,"-test_nd0",&test_nd0,NULL);
29: /* Create a AIJ matrix A */
30: MatCreate(PETSC_COMM_WORLD,&A);
31: MatSetSizes(A,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
32: MatSetType(A,MATAIJ);
33: MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,NULL);
34: MatMPIAIJSetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
35: MatSetFromOptions(A);
36: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
38: /* Create a BAIJ matrix B */
39: MatCreate(PETSC_COMM_WORLD,&B);
40: MatSetSizes(B,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
41: MatSetType(B,MATBAIJ);
42: MatSeqBAIJSetPreallocation(B,bs,PETSC_DEFAULT,NULL);
43: MatMPIBAIJSetPreallocation(B,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
44: MatSetFromOptions(B);
45: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
47: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
48: PetscRandomSetFromOptions(rdm);
50: MatGetOwnershipRange(A,&rstart,&rend);
51: MatGetSize(A,&M,&N);
52: Mbs = M/bs;
54: PetscMalloc1(bs,&rows);
55: PetscMalloc1(bs,&cols);
56: PetscMalloc1(bs*bs,&vals);
57: PetscMalloc1(M,&idx);
59: /* Now set blocks of values */
60: for (i=0; i<40*bs; i++) {
61: PetscRandomGetValue(rdm,&rval);
62: cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
63: PetscRandomGetValue(rdm,&rval);
64: rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m);
65: for (j=1; j<bs; j++) {
66: rows[j] = rows[j-1]+1;
67: cols[j] = cols[j-1]+1;
68: }
70: for (j=0; j<bs*bs; j++) {
71: PetscRandomGetValue(rdm,&rval);
72: vals[j] = rval;
73: }
74: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
75: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
76: }
78: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
79: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
80: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
83: /* Test MatIncreaseOverlap() */
84: PetscMalloc1(nd,&is1);
85: PetscMalloc1(nd,&is2);
87: if (!rank && test_nd0) nd = 0; /* test case */
89: for (i=0; i<nd; i++) {
90: PetscRandomGetValue(rdm,&rval);
91: sz = (int)(PetscRealPart(rval)*m);
92: for (j=0; j<sz; j++) {
93: PetscRandomGetValue(rdm,&rval);
94: idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
95: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
96: }
97: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
98: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
99: }
100: MatIncreaseOverlap(A,nd,is1,ov);
101: MatIncreaseOverlap(B,nd,is2,ov);
103: for (i=0; i<nd; ++i) {
104: ISEqual(is1[i],is2[i],&flg);
106: if (!flg) {
107: 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);
108: }
109: }
111: for (i=0; i<nd; ++i) {
112: ISSort(is1[i]);
113: ISSort(is2[i]);
114: }
116: MatCreateSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
117: MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
119: /* Test MatMult() */
120: for (i=0; i<nd; i++) {
121: MatGetSize(submatA[i],&mm,&nn);
122: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
123: VecDuplicate(xx,&s1);
124: VecDuplicate(xx,&s2);
125: for (j=0; j<3; j++) {
126: VecSetRandom(xx,rdm);
127: MatMult(submatA[i],xx,s1);
128: MatMult(submatB[i],xx,s2);
129: VecNorm(s1,NORM_2,&s1norm);
130: VecNorm(s2,NORM_2,&s2norm);
131: rnorm = s2norm-s1norm;
132: if (rnorm<-tol || rnorm>tol) {
133: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
134: }
135: }
136: VecDestroy(&xx);
137: VecDestroy(&s1);
138: VecDestroy(&s2);
139: }
141: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
142: MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
143: MatCreateSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
145: /* Test MatMult() */
146: for (i=0; i<nd; i++) {
147: MatGetSize(submatA[i],&mm,&nn);
148: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
149: VecDuplicate(xx,&s1);
150: VecDuplicate(xx,&s2);
151: for (j=0; j<3; j++) {
152: VecSetRandom(xx,rdm);
153: MatMult(submatA[i],xx,s1);
154: MatMult(submatB[i],xx,s2);
155: VecNorm(s1,NORM_2,&s1norm);
156: VecNorm(s2,NORM_2,&s2norm);
157: rnorm = s2norm-s1norm;
158: if (rnorm<-tol || rnorm>tol) {
159: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
160: }
161: }
162: VecDestroy(&xx);
163: VecDestroy(&s1);
164: VecDestroy(&s2);
165: }
167: /* Free allocated memory */
168: for (i=0; i<nd; ++i) {
169: ISDestroy(&is1[i]);
170: ISDestroy(&is2[i]);
171: }
172: MatDestroySubMatrices(nd,&submatA);
173: MatDestroySubMatrices(nd,&submatB);
175: PetscFree(is1);
176: PetscFree(is2);
177: PetscFree(idx);
178: PetscFree(rows);
179: PetscFree(cols);
180: PetscFree(vals);
181: MatDestroy(&A);
182: MatDestroy(&B);
183: PetscRandomDestroy(&rdm);
184: PetscFinalize();
185: return ierr;
186: }