Actual source code: ex48.c
petsc-3.6.1 2015-08-06
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,Fact;
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;
21: PetscInitialize(&argc,&args,(char*)0,help);
23: /* Test MatSetValues() and MatGetValues() */
24: PetscOptionsGetInt(NULL,"-mat_block_size",&bs,NULL);
25: PetscOptionsGetInt(NULL,"-mat_size",&m,NULL);
26: M = m*bs;
27: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
28: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
29: MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,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);
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: }
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: }
66: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
67: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
68: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
69: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
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: }
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: }
138: /* Test MatMultAdd() */
139: MatMultAddEqual(A,B,10,&flg);
140: if (!flg) {
141: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n");
142: }
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: PetscMalloc1(M,&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);
167: info.fill = 2.0;
168: info.dtcol = 0.0;
169: info.zeropivot = 1.e-14;
170: info.pivotinblocks = 1.0;
172: if (bs < 4) {
173: MatGetFactor(A,"petsc",MAT_FACTOR_LU,&Fact);
174: MatLUFactorSymbolic(Fact,A,is1,is2,&info);
175: MatLUFactorNumeric(Fact,A,&info);
176: VecSetRandom(yy,rdm);
177: MatForwardSolve(Fact,yy,xx);
178: MatBackwardSolve(Fact,xx,s1);
179: MatDestroy(&Fact);
180: VecScale(s1,-1.0);
181: MatMultAdd(A,s1,yy,yy);
182: VecNorm(yy,NORM_2,&rnorm);
183: if (rnorm<-tol || rnorm>tol) {
184: PetscPrintf(PETSC_COMM_SELF,"Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %D\n",rnorm,bs);
185: }
186: }
188: MatLUFactor(B,is1,is2,&info);
189: MatLUFactor(A,is1,is2,&info);
191: /* Test MatSolveAdd() */
192: for (i=0; i<10; i++) {
193: VecSetRandom(xx,rdm);
194: VecSetRandom(yy,rdm);
195: MatSolveAdd(B,xx,yy,s2);
196: MatSolveAdd(A,xx,yy,s1);
197: VecNorm(s1,NORM_2,&s1norm);
198: VecNorm(s2,NORM_2,&s2norm);
199: rnorm = s2norm-s1norm;
200: if (rnorm<-tol || rnorm>tol) {
201: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
202: }
203: }
205: /* Test MatSolveAdd() when x = A'b +x */
206: for (i=0; i<10; i++) {
207: VecSetRandom(xx,rdm);
208: VecSetRandom(s1,rdm);
209: VecCopy(s2,s1);
210: MatSolveAdd(B,xx,s2,s2);
211: MatSolveAdd(A,xx,s1,s1);
212: VecNorm(s1,NORM_2,&s1norm);
213: VecNorm(s2,NORM_2,&s2norm);
214: rnorm = s2norm-s1norm;
215: if (rnorm<-tol || rnorm>tol) {
216: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
217: }
218: }
220: /* Test MatSolve() */
221: for (i=0; i<10; i++) {
222: VecSetRandom(xx,rdm);
223: MatSolve(B,xx,s2);
224: MatSolve(A,xx,s1);
225: VecNorm(s1,NORM_2,&s1norm);
226: VecNorm(s2,NORM_2,&s2norm);
227: rnorm = s2norm-s1norm;
228: if (rnorm<-tol || rnorm>tol) {
229: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
230: }
231: }
233: /* Test MatSolveTranspose() */
234: if (bs < 8) {
235: for (i=0; i<10; i++) {
236: VecSetRandom(xx,rdm);
237: MatSolveTranspose(B,xx,s2);
238: MatSolveTranspose(A,xx,s1);
239: VecNorm(s1,NORM_2,&s1norm);
240: VecNorm(s2,NORM_2,&s2norm);
241: rnorm = s2norm-s1norm;
242: if (rnorm<-tol || rnorm>tol) {
243: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
244: }
245: }
246: }
248: MatDestroy(&A);
249: MatDestroy(&B);
250: VecDestroy(&xx);
251: VecDestroy(&s1);
252: VecDestroy(&s2);
253: VecDestroy(&yy);
254: ISDestroy(&is1);
255: ISDestroy(&is2);
256: PetscRandomDestroy(&rdm);
257: PetscFinalize();
258: return 0;
259: }