Actual source code: ex48.c
2: static char help[] = "Tests various routines in MatSeqBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A,B,C,D,Fact;
9: Vec xx,s1,s2,yy;
10: PetscInt m=45,rows[2],cols[2],bs=1,i,row,col,*idx,M;
11: PetscScalar rval,vals1[4],vals2[4];
12: PetscRandom rdm;
13: IS is1,is2;
14: PetscReal s1norm,s2norm,rnorm,tol = 1.e-4;
15: PetscBool flg;
16: MatFactorInfo info;
18: PetscInitialize(&argc,&args,(char*)0,help);
19: /* Test MatSetValues() and MatGetValues() */
20: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
21: PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
22: M = m*bs;
23: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
24: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
25: MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B);
26: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
27: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
28: PetscRandomSetFromOptions(rdm);
29: VecCreateSeq(PETSC_COMM_SELF,M,&xx);
30: VecDuplicate(xx,&s1);
31: VecDuplicate(xx,&s2);
32: VecDuplicate(xx,&yy);
34: /* For each row add at least 15 elements */
35: for (row=0; row<M; row++) {
36: for (i=0; i<25*bs; i++) {
37: PetscRandomGetValue(rdm,&rval);
38: col = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
39: MatSetValues(A,1,&row,1,&col,&rval,INSERT_VALUES);
40: MatSetValues(B,1,&row,1,&col,&rval,INSERT_VALUES);
41: }
42: }
44: /* Now set blocks of values */
45: for (i=0; i<20*bs; i++) {
46: PetscRandomGetValue(rdm,&rval);
47: cols[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
48: vals1[0] = rval;
49: PetscRandomGetValue(rdm,&rval);
50: cols[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
51: vals1[1] = rval;
52: PetscRandomGetValue(rdm,&rval);
53: rows[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
54: vals1[2] = rval;
55: PetscRandomGetValue(rdm,&rval);
56: rows[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
57: vals1[3] = rval;
58: MatSetValues(A,2,rows,2,cols,vals1,INSERT_VALUES);
59: MatSetValues(B,2,rows,2,cols,vals1,INSERT_VALUES);
60: }
62: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
63: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
64: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
65: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
67: /* Test MatNorm() */
68: MatNorm(A,NORM_FROBENIUS,&s1norm);
69: MatNorm(B,NORM_FROBENIUS,&s2norm);
70: rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
71: if (rnorm>tol) {
72: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
73: }
74: MatNorm(A,NORM_INFINITY,&s1norm);
75: MatNorm(B,NORM_INFINITY,&s2norm);
76: rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
77: if (rnorm>tol) {
78: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
79: }
80: MatNorm(A,NORM_1,&s1norm);
81: MatNorm(B,NORM_1,&s2norm);
82: rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
83: if (rnorm>tol) {
84: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
85: }
87: /* MatShift() */
88: rval = 10*s1norm;
89: MatShift(A,rval);
90: MatShift(B,rval);
92: /* Test MatTranspose() */
93: MatTranspose(A,MAT_INPLACE_MATRIX,&A);
94: MatTranspose(B,MAT_INPLACE_MATRIX,&B);
96: /* Now do MatGetValues() */
97: for (i=0; i<30; i++) {
98: PetscRandomGetValue(rdm,&rval);
99: cols[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
100: PetscRandomGetValue(rdm,&rval);
101: cols[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
102: PetscRandomGetValue(rdm,&rval);
103: rows[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
104: PetscRandomGetValue(rdm,&rval);
105: rows[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
106: MatGetValues(A,2,rows,2,cols,vals1);
107: MatGetValues(B,2,rows,2,cols,vals2);
108: PetscArraycmp(vals1,vals2,4,&flg);
109: if (!flg) {
110: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %" PetscInt_FMT "\n",bs);
111: }
112: }
114: /* Test MatMult(), MatMultAdd() */
115: for (i=0; i<40; i++) {
116: VecSetRandom(xx,rdm);
117: VecSet(s2,0.0);
118: MatMult(A,xx,s1);
119: MatMultAdd(A,xx,s2,s2);
120: VecNorm(s1,NORM_2,&s1norm);
121: VecNorm(s2,NORM_2,&s2norm);
122: rnorm = s2norm-s1norm;
123: if (rnorm<-tol || rnorm>tol) {
124: PetscPrintf(PETSC_COMM_SELF,"MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
125: }
126: }
128: /* Test MatMult() */
129: MatMultEqual(A,B,10,&flg);
130: if (!flg) {
131: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult()\n");
132: }
134: /* Test MatMultAdd() */
135: MatMultAddEqual(A,B,10,&flg);
136: if (!flg) {
137: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n");
138: }
140: /* Test MatMultTranspose() */
141: MatMultTransposeEqual(A,B,10,&flg);
142: if (!flg) {
143: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTranspose()\n");
144: }
146: /* Test MatMultTransposeAdd() */
147: MatMultTransposeAddEqual(A,B,10,&flg);
148: if (!flg) {
149: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTransposeAdd()\n");
150: }
152: /* Test MatMatMult() */
153: MatCreateSeqDense(PETSC_COMM_SELF,M,40,NULL,&C);
154: MatSetRandom(C,rdm);
155: MatMatMult(A,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
156: MatMatMultEqual(A,C,D,40,&flg);
157: if (!flg) {
158: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n");
159: }
160: MatDestroy(&D);
161: MatCreateSeqDense(PETSC_COMM_SELF,M,40,NULL,&D);
162: MatMatMult(A,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&D);
163: MatMatMultEqual(A,C,D,40,&flg);
164: if (!flg) {
165: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n");
166: }
168: /* Do LUFactor() on both the matrices */
169: PetscMalloc1(M,&idx);
170: for (i=0; i<M; i++) idx[i] = i;
171: ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is1);
172: ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is2);
173: PetscFree(idx);
174: ISSetPermutation(is1);
175: ISSetPermutation(is2);
177: MatFactorInfoInitialize(&info);
179: info.fill = 2.0;
180: info.dtcol = 0.0;
181: info.zeropivot = 1.e-14;
182: info.pivotinblocks = 1.0;
184: if (bs < 4) {
185: MatGetFactor(A,"petsc",MAT_FACTOR_LU,&Fact);
186: MatLUFactorSymbolic(Fact,A,is1,is2,&info);
187: MatLUFactorNumeric(Fact,A,&info);
188: VecSetRandom(yy,rdm);
189: MatForwardSolve(Fact,yy,xx);
190: MatBackwardSolve(Fact,xx,s1);
191: MatDestroy(&Fact);
192: VecScale(s1,-1.0);
193: MatMultAdd(A,s1,yy,yy);
194: VecNorm(yy,NORM_2,&rnorm);
195: if (rnorm<-tol || rnorm>tol) {
196: PetscPrintf(PETSC_COMM_SELF,"Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %" PetscInt_FMT "\n",(double)rnorm,bs);
197: }
198: }
200: MatLUFactor(B,is1,is2,&info);
201: MatLUFactor(A,is1,is2,&info);
203: /* Test MatSolveAdd() */
204: for (i=0; i<10; i++) {
205: VecSetRandom(xx,rdm);
206: VecSetRandom(yy,rdm);
207: MatSolveAdd(B,xx,yy,s2);
208: MatSolveAdd(A,xx,yy,s1);
209: VecNorm(s1,NORM_2,&s1norm);
210: VecNorm(s2,NORM_2,&s2norm);
211: rnorm = s2norm-s1norm;
212: if (rnorm<-tol || rnorm>tol) {
213: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
214: }
215: }
217: /* Test MatSolveAdd() when x = A'b +x */
218: for (i=0; i<10; i++) {
219: VecSetRandom(xx,rdm);
220: VecSetRandom(s1,rdm);
221: VecCopy(s2,s1);
222: MatSolveAdd(B,xx,s2,s2);
223: MatSolveAdd(A,xx,s1,s1);
224: VecNorm(s1,NORM_2,&s1norm);
225: VecNorm(s2,NORM_2,&s2norm);
226: rnorm = s2norm-s1norm;
227: if (rnorm<-tol || rnorm>tol) {
228: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
229: }
230: }
232: /* Test MatSolve() */
233: for (i=0; i<10; i++) {
234: VecSetRandom(xx,rdm);
235: MatSolve(B,xx,s2);
236: MatSolve(A,xx,s1);
237: VecNorm(s1,NORM_2,&s1norm);
238: VecNorm(s2,NORM_2,&s2norm);
239: rnorm = s2norm-s1norm;
240: if (rnorm<-tol || rnorm>tol) {
241: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
242: }
243: }
245: /* Test MatSolveTranspose() */
246: if (bs < 8) {
247: for (i=0; i<10; i++) {
248: VecSetRandom(xx,rdm);
249: MatSolveTranspose(B,xx,s2);
250: MatSolveTranspose(A,xx,s1);
251: VecNorm(s1,NORM_2,&s1norm);
252: VecNorm(s2,NORM_2,&s2norm);
253: rnorm = s2norm-s1norm;
254: if (rnorm<-tol || rnorm>tol) {
255: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %" PetscInt_FMT "\n",(double)s1norm,(double)s2norm,bs);
256: }
257: }
258: }
260: MatDestroy(&A);
261: MatDestroy(&B);
262: MatDestroy(&C);
263: MatDestroy(&D);
264: VecDestroy(&xx);
265: VecDestroy(&s1);
266: VecDestroy(&s2);
267: VecDestroy(&yy);
268: ISDestroy(&is1);
269: ISDestroy(&is2);
270: PetscRandomDestroy(&rdm);
271: PetscFinalize();
272: return 0;
273: }
275: /*TEST
277: test:
278: args: -mat_block_size {{1 2 3 4 5 6 7 8}}
280: TEST*/