Actual source code: ex51.c
2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for MatBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A,B,*submatA,*submatB;
9: PetscInt bs=1,m=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn,lsize;
10: PetscScalar *vals,rval;
11: IS *is1,*is2;
12: PetscRandom rdm;
13: Vec xx,s1,s2;
14: PetscReal s1norm,s2norm,rnorm,tol = PETSC_SQRT_MACHINE_EPSILON;
15: PetscBool flg;
17: PetscInitialize(&argc,&args,(char*)0,help);
18: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
19: PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
20: PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
21: PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
22: M = m*bs;
24: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
25: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
26: MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B);
27: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
28: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
29: PetscRandomSetFromOptions(rdm);
31: PetscMalloc1(bs,&rows);
32: PetscMalloc1(bs,&cols);
33: PetscMalloc1(bs*bs,&vals);
34: PetscMalloc1(M,&idx);
36: /* Now set blocks of values */
37: for (i=0; i<20*bs; i++) {
38: PetscRandomGetValue(rdm,&rval);
39: cols[0] = bs*(int)(PetscRealPart(rval)*m);
40: PetscRandomGetValue(rdm,&rval);
41: rows[0] = bs*(int)(PetscRealPart(rval)*m);
42: for (j=1; j<bs; j++) {
43: rows[j] = rows[j-1]+1;
44: cols[j] = cols[j-1]+1;
45: }
47: for (j=0; j<bs*bs; j++) {
48: PetscRandomGetValue(rdm,&rval);
49: vals[j] = rval;
50: }
51: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
52: MatSetValues(B,bs,rows,bs,cols,vals,ADD_VALUES);
53: }
55: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
56: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
57: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
58: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
60: /* Test MatIncreaseOverlap() */
61: PetscMalloc1(nd,&is1);
62: PetscMalloc1(nd,&is2);
64: for (i=0; i<nd; i++) {
65: PetscRandomGetValue(rdm,&rval);
66: lsize = (int)(PetscRealPart(rval)*m);
67: for (j=0; j<lsize; j++) {
68: PetscRandomGetValue(rdm,&rval);
69: idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
70: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
71: }
72: ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,PETSC_COPY_VALUES,is1+i);
73: ISCreateGeneral(PETSC_COMM_SELF,lsize*bs,idx,PETSC_COPY_VALUES,is2+i);
74: }
75: MatIncreaseOverlap(A,nd,is1,ov);
76: MatIncreaseOverlap(B,nd,is2,ov);
78: for (i=0; i<nd; ++i) {
79: ISEqual(is1[i],is2[i],&flg);
81: }
83: for (i=0; i<nd; ++i) {
84: ISSort(is1[i]);
85: ISSort(is2[i]);
86: }
88: MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
89: MatCreateSubMatrices(B,nd,is2,is2,MAT_INITIAL_MATRIX,&submatB);
91: /* Test MatMult() */
92: for (i=0; i<nd; i++) {
93: MatGetSize(submatA[i],&mm,&nn);
94: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
95: VecDuplicate(xx,&s1);
96: VecDuplicate(xx,&s2);
97: for (j=0; j<3; j++) {
98: VecSetRandom(xx,rdm);
99: MatMult(submatA[i],xx,s1);
100: MatMult(submatB[i],xx,s2);
101: VecNorm(s1,NORM_2,&s1norm);
102: VecNorm(s2,NORM_2,&s2norm);
103: rnorm = s2norm-s1norm;
104: if (rnorm<-tol || rnorm>tol) {
105: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm);
106: }
107: }
108: VecDestroy(&xx);
109: VecDestroy(&s1);
110: VecDestroy(&s2);
111: }
112: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
113: MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
114: MatCreateSubMatrices(B,nd,is2,is2,MAT_REUSE_MATRIX,&submatB);
116: /* Test MatMult() */
117: for (i=0; i<nd; i++) {
118: MatGetSize(submatA[i],&mm,&nn);
119: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
120: VecDuplicate(xx,&s1);
121: VecDuplicate(xx,&s2);
122: for (j=0; j<3; j++) {
123: VecSetRandom(xx,rdm);
124: MatMult(submatA[i],xx,s1);
125: MatMult(submatB[i],xx,s2);
126: VecNorm(s1,NORM_2,&s1norm);
127: VecNorm(s2,NORM_2,&s2norm);
128: rnorm = s2norm-s1norm;
129: if (rnorm<-tol || rnorm>tol) {
130: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm);
131: }
132: }
133: VecDestroy(&xx);
134: VecDestroy(&s1);
135: VecDestroy(&s2);
136: }
138: /* Free allocated memory */
139: for (i=0; i<nd; ++i) {
140: ISDestroy(&is1[i]);
141: ISDestroy(&is2[i]);
142: }
143: MatDestroySubMatrices(nd,&submatA);
144: MatDestroySubMatrices(nd,&submatB);
145: PetscFree(is1);
146: PetscFree(is2);
147: PetscFree(idx);
148: PetscFree(rows);
149: PetscFree(cols);
150: PetscFree(vals);
151: MatDestroy(&A);
152: MatDestroy(&B);
153: PetscRandomDestroy(&rdm);
154: PetscFinalize();
155: return 0;
156: }
158: /*TEST
160: test:
161: args: -mat_block_size {{1 2 5 7 8}} -ov {{1 3}} -mat_size {{11 13}} -nd {{7}}
163: TEST*/