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;
15: PetscRandom rand;
16: IS is1,is2;
17: PetscReal s1norm,s2norm,rnorm,tol = 1.e-10;
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: 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,"Error:MatMult - Norm1=%16.14e Norm2=%16.14e bs = %d \n",s1norm,s2norm,bs);
127: }
128: }
129:
130: /* Test MatMultAdd() */
131: for (i=0; i<40; i++) {
132: VecSetRandom(rand,xx);
133: VecSetRandom(rand,yy);
134: MatMultAdd(A,xx,yy,s1);
135: MatMultAdd(B,xx,yy,s2);
136: VecNorm(s1,NORM_2,&s1norm);
137: VecNorm(s2,NORM_2,&s2norm);
138: rnorm = s2norm-s1norm;
139: if (rnorm<-tol || rnorm>tol) {
140: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
141: }
142: }
143:
144: /* Test MatMultTranspose() */
145: for (i=0; i<40; i++) {
146: VecSetRandom(rand,xx);
147: MatMultTranspose(A,xx,s1);
148: MatMultTranspose(B,xx,s2);
149: VecNorm(s1,NORM_2,&s1norm);
150: VecNorm(s2,NORM_2,&s2norm);
151: rnorm = s2norm-s1norm;
152: if (rnorm<-tol || rnorm>tol) {
153: PetscPrintf(PETSC_COMM_SELF,"Error:MatMultTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
154: }
155: }
156: /* Test MatMultTransposeAdd() */
157: for (i=0; i<40; i++) {
158: VecSetRandom(rand,xx);
159: VecSetRandom(rand,yy);
160: MatMultTransposeAdd(A,xx,yy,s1);
161: MatMultTransposeAdd(B,xx,yy,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:MatMultTransposeAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
167: }
168: }
169:
170:
171: /* Do LUFactor() on both the matrices */
172: PetscMalloc(M*sizeof(int),&idx);
173: for (i=0; i<M; i++) idx[i] = i;
174: ISCreateGeneral(PETSC_COMM_SELF,M,idx,&is1);
175: ISCreateGeneral(PETSC_COMM_SELF,M,idx,&is2);
176: PetscFree(idx);
177: ISSetPermutation(is1);
178: ISSetPermutation(is2);
180: info.fill = 2.0;
181: info.dtcol = 0.0;
182: info.damping = 0.0;
183: info.zeropivot = 1.e-14;
184: info.pivotinblocks = 1.0;
185: MatLUFactor(B,is1,is2,&info);
186: MatLUFactor(A,is1,is2,&info);
187:
188: /* Test MatSolveAdd() */
189: for (i=0; i<40; i++) {
190: VecSetRandom(rand,xx);
191: VecSetRandom(rand,yy);
192: MatSolveAdd(B,xx,yy,s2);
193: MatSolveAdd(A,xx,yy,s1);
194: VecNorm(s1,NORM_2,&s1norm);
195: VecNorm(s2,NORM_2,&s2norm);
196: rnorm = s2norm-s1norm;
197: if (rnorm<-tol || rnorm>tol) {
198: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
199: }
200: }
201:
202: /* Test MatSolveAdd() when x = A'b +x */
203: for (i=0; i<40; i++) {
204: VecSetRandom(rand,xx);
205: VecSetRandom(rand,s1);
206: VecCopy(s2,s1);
207: MatSolveAdd(B,xx,s2,s2);
208: MatSolveAdd(A,xx,s1,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(same) - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
214: }
215: }
216:
217: /* Test MatSolve() */
218: for (i=0; i<40; i++) {
219: VecSetRandom(rand,xx);
220: MatSolve(B,xx,s2);
221: MatSolve(A,xx,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:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
227: }
228: }
229:
230: /* Test MatSolveTranspose() */
231: if (bs < 8) {
232: for (i=0; i<40; i++) {
233: VecSetRandom(rand,xx);
234: MatSolveTranspose(B,xx,s2);
235: MatSolveTranspose(A,xx,s1);
236: VecNorm(s1,NORM_2,&s1norm);
237: VecNorm(s2,NORM_2,&s2norm);
238: rnorm = s2norm-s1norm;
239: if (rnorm<-tol || rnorm>tol) {
240: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %d\n",s1norm,s2norm,bs);
241: }
242: }
243: }
245: MatDestroy(A);
246: MatDestroy(B);
247: VecDestroy(xx);
248: VecDestroy(s1);
249: VecDestroy(s2);
250: VecDestroy(yy);
251: ISDestroy(is1);
252: ISDestroy(is2);
253: PetscRandomDestroy(rand);
254: PetscFinalize();
255: return 0;
256: }