Actual source code: ex54.c
petsc-3.6.4 2016-04-12
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(NULL,"-mat_block_size",&bs,NULL);
26: PetscOptionsGetInt(NULL,"-mat_size",&m,NULL);
27: PetscOptionsGetInt(NULL,"-ov",&ov,NULL);
28: PetscOptionsGetInt(NULL,"-nd",&nd,NULL);
30: /* MatCreateBAIJ(PETSC_COMM_WORLD,bs,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL,&A); */
31: MatCreate(PETSC_COMM_WORLD,&A);
32: MatSetSizes(A,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE);
33: MatSetType(A,MATBAIJ);
34: MatSeqBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL);
35: MatMPIBAIJSetPreallocation(A,bs,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
36: MatSetFromOptions(A);
37: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
39: MatCreateAIJ(PETSC_COMM_WORLD,m*bs,m*bs,PETSC_DECIDE,PETSC_DECIDE,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL,&B);
40: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
41: PetscRandomCreate(PETSC_COMM_WORLD,&rdm);
42: PetscRandomSetFromOptions(rdm);
44: MatGetOwnershipRange(A,&rstart,&rend);
45: MatGetSize(A,&M,&N);
46: Mbs = M/bs;
48: PetscMalloc1(bs,&rows);
49: PetscMalloc1(bs,&cols);
50: PetscMalloc1(bs*bs,&vals);
51: PetscMalloc1(M,&idx);
53: /* Now set blocks of values */
54: for (i=0; i<40*bs; i++) {
55: PetscRandomGetValue(rdm,&rval);
56: cols[0] = bs*(int)(PetscRealPart(rval)*Mbs);
57: PetscRandomGetValue(rdm,&rval);
58: rows[0] = rstart + bs*(int)(PetscRealPart(rval)*m);
59: for (j=1; j<bs; j++) {
60: rows[j] = rows[j-1]+1;
61: cols[j] = cols[j-1]+1;
62: }
64: for (j=0; j<bs*bs; j++) {
65: PetscRandomGetValue(rdm,&rval);
66: vals[j] = rval;
67: }
68: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
69: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
70: }
72: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
73: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
74: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
75: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
77: /* Test MatIncreaseOverlap() */
78: PetscMalloc1(nd,&is1);
79: PetscMalloc1(nd,&is2);
82: for (i=0; i<nd; i++) {
83: PetscRandomGetValue(rdm,&rval);
84: sz = (int)(PetscRealPart(rval)*m);
85: for (j=0; j<sz; j++) {
86: PetscRandomGetValue(rdm,&rval);
87: idx[j*bs] = bs*(int)(PetscRealPart(rval)*Mbs);
88: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
89: }
90: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is1+i);
91: ISCreateGeneral(PETSC_COMM_SELF,sz*bs,idx,PETSC_COPY_VALUES,is2+i);
92: }
93: MatIncreaseOverlap(A,nd,is1,ov);
94: MatIncreaseOverlap(B,nd,is2,ov);
96: for (i=0; i<nd; ++i) {
97: ISEqual(is1[i],is2[i],&flg);
99: if (!flg) {
100: 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);
101: }
102: }
104: for (i=0; i<nd; ++i) {
105: ISSort(is1[i]);
106: ISSort(is2[i]);
107: }
109: MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
110: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
113: /* Test MatMult() */
114: for (i=0; i<nd; i++) {
115: MatGetSize(submatA[i],&mm,&nn);
116: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
117: VecDuplicate(xx,&s1);
118: VecDuplicate(xx,&s2);
119: for (j=0; j<3; j++) {
120: VecSetRandom(xx,rdm);
121: MatMult(submatA[i],xx,s1);
122: MatMult(submatB[i],xx,s2);
123: VecNorm(s1,NORM_2,&s1norm);
124: VecNorm(s2,NORM_2,&s2norm);
125: rnorm = s2norm-s1norm;
126: if (rnorm<-tol || rnorm>tol) {
127: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
128: }
129: }
130: VecDestroy(&xx);
131: VecDestroy(&s1);
132: VecDestroy(&s2);
133: }
135: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
137: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
138: MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
140: /* Test MatMult() */
141: for (i=0; i<nd; i++) {
142: MatGetSize(submatA[i],&mm,&nn);
143: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
144: VecDuplicate(xx,&s1);
145: VecDuplicate(xx,&s2);
146: for (j=0; j<3; j++) {
147: VecSetRandom(xx,rdm);
148: MatMult(submatA[i],xx,s1);
149: MatMult(submatB[i],xx,s2);
150: VecNorm(s1,NORM_2,&s1norm);
151: VecNorm(s2,NORM_2,&s2norm);
152: rnorm = s2norm-s1norm;
153: if (rnorm<-tol || rnorm>tol) {
154: PetscPrintf(PETSC_COMM_SELF,"[%d]Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",rank,s1norm,s2norm);
155: }
156: }
157: VecDestroy(&xx);
158: VecDestroy(&s1);
159: VecDestroy(&s2);
160: }
162: /* Free allocated memory */
163: for (i=0; i<nd; ++i) {
164: ISDestroy(&is1[i]);
165: ISDestroy(&is2[i]);
166: MatDestroy(&submatA[i]);
167: MatDestroy(&submatB[i]);
168: }
169: PetscFree(is1);
170: PetscFree(is2);
171: PetscFree(idx);
172: PetscFree(rows);
173: PetscFree(cols);
174: PetscFree(vals);
175: MatDestroy(&A);
176: MatDestroy(&B);
177: PetscFree(submatA);
178: PetscFree(submatB);
179: PetscRandomDestroy(&rdm);
180: PetscFinalize();
181: return 0;
182: }