Actual source code: ex91.c
2: static char help[] = "Tests MatIncreaseOverlap(), MatCreateSubMatrices() for sequential MatSBAIJ format. Derived from ex51.c\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A,Atrans,sA,*submatA,*submatsA;
9: PetscInt bs=1,m=43,ov=1,i,j,k,*rows,*cols,M,nd=5,*idx,mm,nn;
10: PetscMPIInt size;
11: PetscScalar *vals,rval,one=1.0;
12: IS *is1,*is2;
13: PetscRandom rand;
14: Vec xx,s1,s2;
15: PetscReal s1norm,s2norm,rnorm,tol = 10*PETSC_SMALL;
16: PetscBool flg;
18: PetscInitialize(&argc,&args,(char*)0,help);
19: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
20: PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
21: PetscOptionsGetInt(NULL,NULL,"-ov",&ov,NULL);
22: PetscOptionsGetInt(NULL,NULL,"-nd",&nd,NULL);
24: /* create a SeqBAIJ matrix A */
25: M = m*bs;
26: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
27: MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
28: PetscRandomCreate(PETSC_COMM_SELF,&rand);
29: PetscRandomSetFromOptions(rand);
31: PetscMalloc1(bs,&rows);
32: PetscMalloc1(bs,&cols);
33: PetscMalloc1(bs*bs,&vals);
34: PetscMalloc1(M,&idx);
36: /* Now set blocks of random values */
37: /* first, set diagonal blocks as zero */
38: for (j=0; j<bs*bs; j++) vals[j] = 0.0;
39: for (i=0; i<m; i++) {
40: cols[0] = i*bs; rows[0] = i*bs;
41: for (j=1; j<bs; j++) {
42: rows[j] = rows[j-1]+1;
43: cols[j] = cols[j-1]+1;
44: }
45: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
46: }
47: /* second, add random blocks */
48: for (i=0; i<20*bs; i++) {
49: PetscRandomGetValue(rand,&rval);
50: cols[0] = bs*(int)(PetscRealPart(rval)*m);
51: PetscRandomGetValue(rand,&rval);
52: rows[0] = bs*(int)(PetscRealPart(rval)*m);
53: for (j=1; j<bs; j++) {
54: rows[j] = rows[j-1]+1;
55: cols[j] = cols[j-1]+1;
56: }
58: for (j=0; j<bs*bs; j++) {
59: PetscRandomGetValue(rand,&rval);
60: vals[j] = rval;
61: }
62: MatSetValues(A,bs,rows,bs,cols,vals,ADD_VALUES);
63: }
65: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
68: /* make A a symmetric matrix: A <- A^T + A */
69: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
70: MatAXPY(A,one,Atrans,DIFFERENT_NONZERO_PATTERN);
71: MatDestroy(&Atrans);
72: MatTranspose(A,MAT_INITIAL_MATRIX, &Atrans);
73: MatEqual(A, Atrans, &flg);
75: MatDestroy(&Atrans);
77: /* create a SeqSBAIJ matrix sA (= A) */
78: MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
79: MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);
81: /* Test sA==A through MatMult() */
82: for (i=0; i<nd; i++) {
83: MatGetSize(A,&mm,&nn);
84: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
85: VecDuplicate(xx,&s1);
86: VecDuplicate(xx,&s2);
87: for (j=0; j<3; j++) {
88: VecSetRandom(xx,rand);
89: MatMult(A,xx,s1);
90: MatMult(sA,xx,s2);
91: VecNorm(s1,NORM_2,&s1norm);
92: VecNorm(s2,NORM_2,&s2norm);
93: rnorm = s2norm-s1norm;
94: if (rnorm<-tol || rnorm>tol) {
95: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm);
96: }
97: }
98: VecDestroy(&xx);
99: VecDestroy(&s1);
100: VecDestroy(&s2);
101: }
103: /* Test MatIncreaseOverlap() */
104: PetscMalloc1(nd,&is1);
105: PetscMalloc1(nd,&is2);
107: for (i=0; i<nd; i++) {
108: PetscRandomGetValue(rand,&rval);
109: size = (int)(PetscRealPart(rval)*m);
110: for (j=0; j<size; j++) {
111: PetscRandomGetValue(rand,&rval);
112: idx[j*bs] = bs*(int)(PetscRealPart(rval)*m);
113: for (k=1; k<bs; k++) idx[j*bs+k] = idx[j*bs]+k;
114: }
115: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is1+i);
116: ISCreateGeneral(PETSC_COMM_SELF,size*bs,idx,PETSC_COPY_VALUES,is2+i);
117: }
118: /* for debugging */
119: /*
120: MatView(A,PETSC_VIEWER_STDOUT_SELF);
121: MatView(sA,PETSC_VIEWER_STDOUT_SELF);
122: */
124: MatIncreaseOverlap(A,nd,is1,ov);
125: MatIncreaseOverlap(sA,nd,is2,ov);
127: for (i=0; i<nd; ++i) {
128: ISSort(is1[i]);
129: ISSort(is2[i]);
130: }
132: for (i=0; i<nd; ++i) {
133: ISEqual(is1[i],is2[i],&flg);
135: }
137: MatCreateSubMatrices(A,nd,is1,is1,MAT_INITIAL_MATRIX,&submatA);
138: MatCreateSubMatrices(sA,nd,is2,is2,MAT_INITIAL_MATRIX,&submatsA);
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,rand);
148: MatMult(submatA[i],xx,s1);
149: MatMult(submatsA[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,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm);
155: }
156: }
157: VecDestroy(&xx);
158: VecDestroy(&s1);
159: VecDestroy(&s2);
160: }
162: /* Now test MatCreateSubmatrices with MAT_REUSE_MATRIX option */
163: MatCreateSubMatrices(A,nd,is1,is1,MAT_REUSE_MATRIX,&submatA);
164: MatCreateSubMatrices(sA,nd,is2,is2,MAT_REUSE_MATRIX,&submatsA);
166: /* Test MatMult() */
167: for (i=0; i<nd; i++) {
168: MatGetSize(submatA[i],&mm,&nn);
169: VecCreateSeq(PETSC_COMM_SELF,mm,&xx);
170: VecDuplicate(xx,&s1);
171: VecDuplicate(xx,&s2);
172: for (j=0; j<3; j++) {
173: VecSetRandom(xx,rand);
174: MatMult(submatA[i],xx,s1);
175: MatMult(submatsA[i],xx,s2);
176: VecNorm(s1,NORM_2,&s1norm);
177: VecNorm(s2,NORM_2,&s2norm);
178: rnorm = s2norm-s1norm;
179: if (rnorm<-tol || rnorm>tol) {
180: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e\n",(double)s1norm,(double)s2norm);
181: }
182: }
183: VecDestroy(&xx);
184: VecDestroy(&s1);
185: VecDestroy(&s2);
186: }
188: /* Free allocated memory */
189: for (i=0; i<nd; ++i) {
190: ISDestroy(&is1[i]);
191: ISDestroy(&is2[i]);
192: }
193: MatDestroySubMatrices(nd,&submatA);
194: MatDestroySubMatrices(nd,&submatsA);
196: PetscFree(is1);
197: PetscFree(is2);
198: PetscFree(idx);
199: PetscFree(rows);
200: PetscFree(cols);
201: PetscFree(vals);
202: MatDestroy(&A);
203: MatDestroy(&sA);
204: PetscRandomDestroy(&rand);
205: PetscFinalize();
206: return 0;
207: }
209: /*TEST
211: test:
212: args: -ov 2
214: TEST*/