Actual source code: ex48.c
petsc-3.4.5 2014-06-29
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;
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: 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);
167: info.fill = 2.0;
168: info.dtcol = 0.0;
169: info.zeropivot = 1.e-14;
170: info.pivotinblocks = 1.0;
172: MatLUFactor(B,is1,is2,&info);
173: MatLUFactor(A,is1,is2,&info);
175: /* Test MatSolveAdd() */
176: for (i=0; i<10; i++) {
177: VecSetRandom(xx,rdm);
178: VecSetRandom(yy,rdm);
179: MatSolveAdd(B,xx,yy,s2);
180: MatSolveAdd(A,xx,yy,s1);
181: VecNorm(s1,NORM_2,&s1norm);
182: VecNorm(s2,NORM_2,&s2norm);
183: rnorm = s2norm-s1norm;
184: if (rnorm<-tol || rnorm>tol) {
185: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
186: }
187: }
189: /* Test MatSolveAdd() when x = A'b +x */
190: for (i=0; i<10; i++) {
191: VecSetRandom(xx,rdm);
192: VecSetRandom(s1,rdm);
193: VecCopy(s2,s1);
194: MatSolveAdd(B,xx,s2,s2);
195: MatSolveAdd(A,xx,s1,s1);
196: VecNorm(s1,NORM_2,&s1norm);
197: VecNorm(s2,NORM_2,&s2norm);
198: rnorm = s2norm-s1norm;
199: if (rnorm<-tol || rnorm>tol) {
200: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
201: }
202: }
204: /* Test MatSolve() */
205: for (i=0; i<10; i++) {
206: VecSetRandom(xx,rdm);
207: MatSolve(B,xx,s2);
208: MatSolve(A,xx,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:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
214: }
215: }
217: /* Test MatSolveTranspose() */
218: if (bs < 8) {
219: for (i=0; i<10; i++) {
220: VecSetRandom(xx,rdm);
221: MatSolveTranspose(B,xx,s2);
222: MatSolveTranspose(A,xx,s1);
223: VecNorm(s1,NORM_2,&s1norm);
224: VecNorm(s2,NORM_2,&s2norm);
225: rnorm = s2norm-s1norm;
226: if (rnorm<-tol || rnorm>tol) {
227: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
228: }
229: }
230: }
232: MatDestroy(&A);
233: MatDestroy(&B);
234: VecDestroy(&xx);
235: VecDestroy(&s1);
236: VecDestroy(&s2);
237: VecDestroy(&yy);
238: ISDestroy(&is1);
239: ISDestroy(&is2);
240: PetscRandomDestroy(&rdm);
241: PetscFinalize();
242: return 0;
243: }