Actual source code: ex51.c
2: static char help[] = "Tests MatIncreaseOverlap(), MatGetSubMatrices() for 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=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn,lsize;
13: PetscScalar *vals,rval;
14: IS *is1,*is2;
15: PetscRandom rdm;
16: Vec xx,s1,s2;
17: PetscReal s1norm,s2norm,rnorm,tol = 1.e-10;
18: PetscTruth flg;
20: PetscInitialize(&argc,&args,(char *)0,help);
21:
23: PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
24: PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);
25: PetscOptionsGetInt(PETSC_NULL,"-ov",&ov,PETSC_NULL);
26: PetscOptionsGetInt(PETSC_NULL,"-nd",&nd,PETSC_NULL);
27: M = m*bs;
29: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,PETSC_NULL,&A);
30: MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,PETSC_NULL,&B);
31: PetscRandomCreate(PETSC_COMM_SELF,RANDOM_DEFAULT,&rdm);
33: PetscMalloc(bs*sizeof(PetscInt),&rows);
34: PetscMalloc(bs*sizeof(PetscInt),&cols);
35: PetscMalloc(bs*bs*sizeof(PetscScalar),&vals);
36: PetscMalloc(M*sizeof(PetscScalar),&idx);
37:
38: /* Now set blocks of values */
39: for (i=0; i<20*bs; i++) {
40: PetscRandomGetValue(rdm,&rval);
41: cols[0] = bs*(int)(PetscRealPart(rval)*m);
42: PetscRandomGetValue(rdm,&rval);
43: rows[0] = bs*(int)(PetscRealPart(rval)*m);
44: for (j=1; j<bs; j++) {
45: rows[j] = rows[j-1]+1;
46: cols[j] = cols[j-1]+1;
47: }
49: for (j=0; j<bs*bs; j++) {
50: PetscRandomGetValue(rdm,&rval);
51: vals[j] = rval;
52: }
53: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
54: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
55: }
57: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
58: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
59: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
60: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
62: /* Test MatIncreaseOverlap() */
63: PetscMalloc(nd*sizeof(IS **),&is1);
64: PetscMalloc(nd*sizeof(IS **),&is2);
66:
67: for (i=0; i<nd; i++) {
68: PetscRandomGetValue(rdm,&rval);
69: lsize = (int)(PetscRealPart(rval)*m);
70: for (j=0; j<lsize; j++) {
71: PetscRandomGetValue(rdm,&rval);
72: idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
73: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
74: }
75: ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,is1+i);
76: ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,is2+i);
77: }
78: MatIncreaseOverlap(A,nd,is1,ov);
79: MatIncreaseOverlap(B,nd,is2,ov);
81: for (i=0; i<nd; ++i) {
82: ISEqual(is1[i],is2[i],&flg);
83: PetscPrintf(PETSC_COMM_SELF,"i=%D, flg =%d\n",i,(int)flg);
84: }
86: for (i=0; i<nd; ++i) {
87: ISSort(is1[i]);
88: ISSort(is2[i]);
89: }
90:
91: MatGetSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
92: MatGetSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
94: /* Test MatMult() */
95: for (i=0; i<nd; i++) {
96: MatGetSize(submatA[i],&mm,&nn);
97: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
98: VecDuplicate(xx,&s1);
99: VecDuplicate(xx,&s2);
100: for (j=0; j<3; j++) {
101: VecSetRandom(xx,rdm);
102: MatMult(submatA[i],xx,s1);
103: MatMult(submatB[i],xx,s2);
104: VecNorm(s1,NORM_2,&s1norm);
105: VecNorm(s2,NORM_2,&s2norm);
106: rnorm = s2norm-s1norm;
107: if (rnorm<-tol || rnorm>tol) {
108: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
109: }
110: }
111: VecDestroy(xx);
112: VecDestroy(s1);
113: VecDestroy(s2);
114: }
115: /* Now test MatGetSubmatrices with MAT_REUSE_MATRIX option */
116: MatGetSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
117: MatGetSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
118:
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,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",s1norm,s2norm);
134: }
135: }
136: VecDestroy(xx);
137: VecDestroy(s1);
138: VecDestroy(s2);
139: }
140:
141: /* Free allocated memory */
142: for (i=0; i<nd; ++i) {
143: ISDestroy(is1[i]);
144: ISDestroy(is2[i]);
145: MatDestroy(submatA[i]);
146: MatDestroy(submatB[i]);
147: }
148: PetscFree(is1);
149: PetscFree(is2);
150: PetscFree(idx);
151: PetscFree(rows);
152: PetscFree(cols);
153: PetscFree(vals);
154: MatDestroy(A);
155: MatDestroy(B);
156: PetscFree(submatA);
157: PetscFree(submatB);
158: PetscRandomDestroy(rdm);
159: PetscFinalize();
160: return 0;
161: }