Actual source code: ex48.c
1: /*$Id: ex48.c,v 1.25 2001/08/07 21:30:08 bsmith Exp $*/
3: static char help[] = "Tests the vatious routines in MatBAIJ format.\n";
5: #include petscmat.h
9: int main(int argc,char **args)
10: {
11: Mat A,B;
12: Vec xx,s1,s2,yy;
13: int m=45,ierr,rows[2],cols[2],bs=1,i,row,col,*idx,M;
14: PetscScalar rval,vals1[4],vals2[4],zero=0.0,neg_one=-1.0;
15: PetscRandom rand;
16: IS is1,is2;
17: PetscReal s1norm,s2norm,rnorm,tol = 1.e-8;
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,RANDOM_DEFAULT,&rand);
30: VecCreateSeq(PETSC_COMM_SELF,M,&xx);
31: VecDuplicate(xx,&s1);
32: VecDuplicate(xx,&s2);
33: VecDuplicate(xx,&yy);
34:
35: /* For each row add atleast 15 elements */
36: for (row=0; row<M; row++) {
37: for (i=0; i<25*bs; i++) {
38: PetscRandomGetValue(rand,&rval);
39: col = (int)(PetscRealPart(rval)*M);
40: MatSetValues(A,1,&row,1,&col,&rval,ADD_VALUES);
41: MatSetValues(B,1,&row,1,&col,&rval,ADD_VALUES);
42: }
43: }
44:
45: /* Now set blocks of values */
46: for (i=0; i<20*bs; i++) {
47: PetscRandomGetValue(rand,&rval);
48: cols[0] = (int)(PetscRealPart(rval)*M);
49: vals1[0] = rval;
50: PetscRandomGetValue(rand,&rval);
51: cols[1] = (int)(PetscRealPart(rval)*M);
52: vals1[1] = rval;
53: PetscRandomGetValue(rand,&rval);
54: rows[0] = (int)(PetscRealPart(rval)*M);
55: vals1[2] = rval;
56: PetscRandomGetValue(rand,&rval);
57: rows[1] = (int)(PetscRealPart(rval)*M);
58: vals1[3] = rval;
59: MatSetValues(A,2,rows,2,cols,vals1,ADD_VALUES);
60: MatSetValues(B,2,rows,2,cols,vals1,ADD_VALUES);
61: }
62:
63: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
64: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
65: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
67:
68: /* Test MatNorm() */
69: MatNorm(A,NORM_FROBENIUS,&s1norm);
70: MatNorm(B,NORM_FROBENIUS,&s2norm);
71: rnorm = s2norm-s1norm;
72: if (rnorm<-tol || rnorm>tol) {
73: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm()- Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
74: }
75: /* MatScale() */
76: rval = 10*s1norm;
77: MatShift(&rval,A);
78: MatShift(&rval,B);
79:
80: /* Test MatTranspose() */
81: MatTranspose(A,PETSC_NULL);
82: MatTranspose(B,PETSC_NULL);
83:
84:
85: /* Now do MatGetValues() */
86: for (i=0; i<30; i++) {
87: PetscRandomGetValue(rand,&rval);
88: cols[0] = (int)(PetscRealPart(rval)*M);
89: PetscRandomGetValue(rand,&rval);
90: cols[1] = (int)(PetscRealPart(rval)*M);
91: PetscRandomGetValue(rand,&rval);
92: rows[0] = (int)(PetscRealPart(rval)*M);
93: PetscRandomGetValue(rand,&rval);
94: rows[1] = (int)(PetscRealPart(rval)*M);
95: MatGetValues(A,2,rows,2,cols,vals1);
96: MatGetValues(B,2,rows,2,cols,vals2);
97: PetscMemcmp(vals1,vals2,4*sizeof(PetscScalar),&flg);
98: if (!flg) {
99: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %d\n",bs);
100: }
101: }
102:
103: /* Test MatMult(), MatMultAdd() */
104: for (i=0; i<40; i++) {
105: VecSetRandom(rand,xx);
106: VecSet(&zero,s2);
107: MatMult(A,xx,s1);
108: MatMultAdd(A,xx,s2,s2);
109: VecNorm(s1,NORM_2,&s1norm);
110: VecNorm(s2,NORM_2,&s2norm);
111: rnorm = s2norm-s1norm;
112: if (rnorm<-tol || rnorm>tol) {
113: PetscPrintf(PETSC_COMM_SELF,"MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %d\n",s1norm,s2norm,bs);
114: }
115: }
117: /* Test MatMult() */
118: for (i=0; i<40; i++) {
119: VecSetRandom(rand,xx);
120: MatMult(A,xx,s1);
121: MatMult(B,xx,s2);
122: VecAXPY(&neg_one,s1,s2);
123: VecNorm(s2,NORM_1,&s2norm);
124: if (s2norm >tol) {
125: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), ||s1-s2||=%g\n",s2norm);
126: }
127: /*
128: VecNorm(s1,NORM_2,&s1norm);
129: VecNorm(s2,NORM_2,&s2norm);
130: rnorm = s2norm-s1norm;
131: if (rnorm<-tol || rnorm>tol) {
132: PetscPrintf(PETSC_COMM_SELF,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d \n",s1norm,s2norm,bs);
133: */
134: }
135:
136: /* Test MatMultAdd() */
137: for (i=0; i<40; i++) {
138: VecSetRandom(rand,xx);
139: VecSetRandom(rand,yy);
140: MatMultAdd(A,xx,yy,s1);
141: MatMultAdd(B,xx,yy,s2);
142: VecAXPY(&neg_one,s1,s2);
143: VecNorm(s2,NORM_1,&s2norm);
144: if (s2norm >tol) {
145: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd(), ||s1-s2||=%g\n",s2norm);
146: }
147: /*
148: VecNorm(s1,NORM_2,&s1norm);
149: VecNorm(s2,NORM_2,&s2norm);
150: rnorm = s2norm-s1norm;
151: if (rnorm<-tol || rnorm>tol) {
152: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
153: }
154: */
155: }
156:
157: /* Test MatMultTranspose() */
158: for (i=0; i<40; i++) {
159: VecSetRandom(rand,xx);
160: MatMultTranspose(A,xx,s1);
161: MatMultTranspose(B,xx,s2);
162: VecNorm(s1,NORM_2,&s1norm);
163: VecNorm(s2,NORM_2,&s2norm);
164: rnorm = s2norm-s1norm;
165: if (rnorm<-tol || rnorm>tol) {
166: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
167: }
168: }
169: /* Test MatMultTransposeAdd() */
170: for (i=0; i<40; i++) {
171: VecSetRandom(rand,xx);
172: VecSetRandom(rand,yy);
173: MatMultTransposeAdd(A,xx,yy,s1);
174: MatMultTransposeAdd(B,xx,yy,s2);
175: VecNorm(s1,NORM_2,&s1norm);
176: VecNorm(s2,NORM_2,&s2norm);
177: rnorm = s2norm-s1norm;
178: if (rnorm<-tol || rnorm>tol) {
179: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
180: }
181: }
182:
183:
184: /* Do LUFactor() on both the matrices */
185: PetscMalloc(M*sizeof(int),&idx);
186: for (i=0; i<M; i++) idx[i] = i;
187: ISCreateGeneral(PETSC_COMM_SELF,M,idx,&is1);
188: ISCreateGeneral(PETSC_COMM_SELF,M,idx,&is2);
189: PetscFree(idx);
190: ISSetPermutation(is1);
191: ISSetPermutation(is2);
193: info.fill = 2.0;
194: info.dtcol = 0.0;
195: info.damping = 0.0;
196: info.zeropivot = 1.e-14;
197: info.pivotinblocks = 1.0;
198: MatLUFactor(B,is1,is2,&info);
199: MatLUFactor(A,is1,is2,&info);
200:
201: /* Test MatSolveAdd() */
202: for (i=0; i<40; i++) {
203: VecSetRandom(rand,xx);
204: VecSetRandom(rand,yy);
205: MatSolveAdd(B,xx,yy,s2);
206: MatSolveAdd(A,xx,yy,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:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
212: }
213: }
214:
215: /* Test MatSolveAdd() when x = A'b +x */
216: for (i=0; i<40; i++) {
217: VecSetRandom(rand,xx);
218: VecSetRandom(rand,s1);
219: VecCopy(s2,s1);
220: MatSolveAdd(B,xx,s2,s2);
221: MatSolveAdd(A,xx,s1,s1);
222: VecNorm(s1,NORM_2,&s1norm);
223: VecNorm(s2,NORM_2,&s2norm);
224: rnorm = s2norm-s1norm;
225: if (rnorm<-tol || rnorm>tol) {
226: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
227: }
228: }
229:
230: /* Test MatSolve() */
231: for (i=0; i<40; i++) {
232: VecSetRandom(rand,xx);
233: MatSolve(B,xx,s2);
234: MatSolve(A,xx,s1);
235: VecNorm(s1,NORM_2,&s1norm);
236: VecNorm(s2,NORM_2,&s2norm);
237: rnorm = s2norm-s1norm;
238: if (rnorm<-tol || rnorm>tol) {
239: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
240: }
241: }
242:
243: /* Test MatSolveTranspose() */
244: if (bs < 8) {
245: for (i=0; i<40; i++) {
246: VecSetRandom(rand,xx);
247: MatSolveTranspose(B,xx,s2);
248: MatSolveTranspose(A,xx,s1);
249: VecNorm(s1,NORM_2,&s1norm);
250: VecNorm(s2,NORM_2,&s2norm);
251: rnorm = s2norm-s1norm;
252: if (rnorm<-tol || rnorm>tol) {
253: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
254: }
255: }
256: }
258: MatDestroy(A);
259: MatDestroy(B);
260: VecDestroy(xx);
261: VecDestroy(s1);
262: VecDestroy(s2);
263: VecDestroy(yy);
264: ISDestroy(is1);
265: ISDestroy(is2);
266: PetscRandomDestroy(rand);
267: PetscFinalize();
268: return 0;
269: }