Actual source code: ex48.c
petsc-3.3-p6 2013-02-11
2: static char help[] = "Tests the vatious routines in MatSeqBAIJ format.\n";
4: #include <petscmat.h>
8: int main(int argc,char **args)
9: {
10: Mat A,B;
11: Vec xx,s1,s2,yy;
13: PetscInt m=45,rows[2],cols[2],bs=1,i,row,col,*idx,M;
14: PetscScalar rval,vals1[4],vals2[4];
15: PetscRandom rdm;
16: IS is1,is2;
17: PetscReal s1norm,s2norm,rnorm,tol = 1.e-4;
18: PetscBool flg;
19: MatFactorInfo info;
20:
21: PetscInitialize(&argc,&args,(char *)0,help);
22:
23: /* Test MatSetValues() and MatGetValues() */
24: PetscOptionsGetInt(PETSC_NULL,"-mat_block_size",&bs,PETSC_NULL);
25: PetscOptionsGetInt(PETSC_NULL,"-mat_size",&m,PETSC_NULL);
26: M = m*bs;
27: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,PETSC_NULL,&A);
28: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
29: MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,PETSC_NULL,&B);
30: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
31: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
32: PetscRandomSetFromOptions(rdm);
33: VecCreateSeq(PETSC_COMM_SELF,M,&xx);
34: VecDuplicate(xx,&s1);
35: VecDuplicate(xx,&s2);
36: VecDuplicate(xx,&yy);
37:
38: /* For each row add atleast 15 elements */
39: for (row=0; row<M; row++) {
40: for (i=0; i<25*bs; i++) {
41: PetscRandomGetValue(rdm,&rval);
42: col = (PetscInt)(PetscRealPart(rval)*M);
43: MatSetValues(A,1,&row,1,&col,&rval,INSERT_VALUES);
44: MatSetValues(B,1,&row,1,&col,&rval,INSERT_VALUES);
45: }
46: }
47:
48: /* Now set blocks of values */
49: for (i=0; i<20*bs; i++) {
50: PetscRandomGetValue(rdm,&rval);
51: cols[0] = (PetscInt)(PetscRealPart(rval)*M);
52: vals1[0] = rval;
53: PetscRandomGetValue(rdm,&rval);
54: cols[1] = (PetscInt)(PetscRealPart(rval)*M);
55: vals1[1] = rval;
56: PetscRandomGetValue(rdm,&rval);
57: rows[0] = (PetscInt)(PetscRealPart(rval)*M);
58: vals1[2] = rval;
59: PetscRandomGetValue(rdm,&rval);
60: rows[1] = (PetscInt)(PetscRealPart(rval)*M);
61: vals1[3] = rval;
62: MatSetValues(A,2,rows,2,cols,vals1,INSERT_VALUES);
63: MatSetValues(B,2,rows,2,cols,vals1,INSERT_VALUES);
64: }
65:
66: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
68: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
70:
71: /* Test MatNorm() */
72: MatNorm(A,NORM_FROBENIUS,&s1norm);
73: MatNorm(B,NORM_FROBENIUS,&s2norm);
74: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
75: if ( rnorm>tol ) {
76: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);
77: }
78: MatNorm(A,NORM_INFINITY,&s1norm);
79: MatNorm(B,NORM_INFINITY,&s2norm);
80: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
81: if ( rnorm>tol ) {
82: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);
83: }
84: MatNorm(A,NORM_1,&s1norm);
85: MatNorm(B,NORM_1,&s2norm);
86: rnorm = PetscAbsScalar(s2norm-s1norm)/s2norm;
87: if ( rnorm>tol ) {
88: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);
89: }
91: /* MatShift() */
92: rval = 10*s1norm;
93: MatShift(A,rval);
94: MatShift(B,rval);
96: /* Test MatTranspose() */
97: MatTranspose(A,MAT_REUSE_MATRIX,&A);
98: MatTranspose(B,MAT_REUSE_MATRIX,&B);
100: /* Now do MatGetValues() */
101: for (i=0; i<30; i++) {
102: PetscRandomGetValue(rdm,&rval);
103: cols[0] = (PetscInt)(PetscRealPart(rval)*M);
104: PetscRandomGetValue(rdm,&rval);
105: cols[1] = (PetscInt)(PetscRealPart(rval)*M);
106: PetscRandomGetValue(rdm,&rval);
107: rows[0] = (PetscInt)(PetscRealPart(rval)*M);
108: PetscRandomGetValue(rdm,&rval);
109: rows[1] = (PetscInt)(PetscRealPart(rval)*M);
110: MatGetValues(A,2,rows,2,cols,vals1);
111: MatGetValues(B,2,rows,2,cols,vals2);
112: PetscMemcmp(vals1,vals2,4*sizeof(PetscScalar),&flg);
113: if (!flg) {
114: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %D\n",bs);
115: }
116: }
117:
118: /* Test MatMult(), MatMultAdd() */
119: for (i=0; i<40; i++) {
120: VecSetRandom(xx,rdm);
121: VecSet(s2,0.0);
122: MatMult(A,xx,s1);
123: MatMultAdd(A,xx,s2,s2);
124: VecNorm(s1,NORM_2,&s1norm);
125: VecNorm(s2,NORM_2,&s2norm);
126: rnorm = s2norm-s1norm;
127: if (rnorm<-tol || rnorm>tol) {
128: PetscPrintf(PETSC_COMM_SELF,"MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %D\n",s1norm,s2norm,bs);
129: }
130: }
132: /* Test MatMult() */
133: MatMultEqual(A,B,10,&flg);
134: if (!flg){
135: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult()\n");
136: }
137:
138: /* Test MatMultAdd() */
139: MatMultAddEqual(A,B,10,&flg);
140: if (!flg){
141: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n");
142: }
143:
144: /* Test MatMultTranspose() */
145: MatMultTransposeEqual(A,B,10,&flg);
146: if (!flg){
147: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTranspose()\n");
148: }
150: /* Test MatMultTransposeAdd() */
151: MatMultTransposeAddEqual(A,B,10,&flg);
152: if (!flg){
153: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTransposeAdd()\n");
154: }
156: /* Do LUFactor() on both the matrices */
157: PetscMalloc(M*sizeof(PetscInt),&idx);
158: for (i=0; i<M; i++) idx[i] = i;
159: ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is1);
160: ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is2);
161: PetscFree(idx);
162: ISSetPermutation(is1);
163: ISSetPermutation(is2);
165: MatFactorInfoInitialize(&info);
166: info.fill = 2.0;
167: info.dtcol = 0.0;
168: info.zeropivot = 1.e-14;
169: info.pivotinblocks = 1.0;
170: MatLUFactor(B,is1,is2,&info);
171: MatLUFactor(A,is1,is2,&info);
172:
173: /* Test MatSolveAdd() */
174: for (i=0; i<10; i++) {
175: VecSetRandom(xx,rdm);
176: VecSetRandom(yy,rdm);
177: MatSolveAdd(B,xx,yy,s2);
178: MatSolveAdd(A,xx,yy,s1);
179: VecNorm(s1,NORM_2,&s1norm);
180: VecNorm(s2,NORM_2,&s2norm);
181: rnorm = s2norm-s1norm;
182: if (rnorm<-tol || rnorm>tol) {
183: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
184: }
185: }
186:
187: /* Test MatSolveAdd() when x = A'b +x */
188: for (i=0; i<10; i++) {
189: VecSetRandom(xx,rdm);
190: VecSetRandom(s1,rdm);
191: VecCopy(s2,s1);
192: MatSolveAdd(B,xx,s2,s2);
193: MatSolveAdd(A,xx,s1,s1);
194: VecNorm(s1,NORM_2,&s1norm);
195: VecNorm(s2,NORM_2,&s2norm);
196: rnorm = s2norm-s1norm;
197: if (rnorm<-tol || rnorm>tol) {
198: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
199: }
200: }
201:
202: /* Test MatSolve() */
203: for (i=0; i<10; i++) {
204: VecSetRandom(xx,rdm);
205: MatSolve(B,xx,s2);
206: MatSolve(A,xx,s1);
207: VecNorm(s1,NORM_2,&s1norm);
208: VecNorm(s2,NORM_2,&s2norm);
209: rnorm = s2norm-s1norm;
210: if (rnorm<-tol || rnorm>tol) {
211: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
212: }
213: }
214:
215: /* Test MatSolveTranspose() */
216: if (bs < 8) {
217: for (i=0; i<10; i++) {
218: VecSetRandom(xx,rdm);
219: MatSolveTranspose(B,xx,s2);
220: MatSolveTranspose(A,xx,s1);
221: VecNorm(s1,NORM_2,&s1norm);
222: VecNorm(s2,NORM_2,&s2norm);
223: rnorm = s2norm-s1norm;
224: if (rnorm<-tol || rnorm>tol) {
225: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
226: }
227: }
228: }
230: MatDestroy(&A);
231: MatDestroy(&B);
232: VecDestroy(&xx);
233: VecDestroy(&s1);
234: VecDestroy(&s2);
235: VecDestroy(&yy);
236: ISDestroy(&is1);
237: ISDestroy(&is2);
238: PetscRandomDestroy(&rdm);
239: PetscFinalize();
240: return 0;
241: }