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