Actual source code: ex48.c
petsc-3.8.4 2018-03-24
2: static char help[] = "Tests the vatious routines in MatSeqBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A,B,Fact;
9: Vec xx,s1,s2,yy;
11: PetscInt m=45,rows[2],cols[2],bs=1,i,row,col,*idx,M;
12: PetscScalar rval,vals1[4],vals2[4];
13: PetscRandom rdm;
14: IS is1,is2;
15: PetscReal s1norm,s2norm,rnorm,tol = 1.e-4;
16: PetscBool flg;
17: MatFactorInfo info;
19: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
20: /* Test MatSetValues() and MatGetValues() */
21: PetscOptionsGetInt(NULL,NULL,"-mat_block_size",&bs,NULL);
22: PetscOptionsGetInt(NULL,NULL,"-mat_size",&m,NULL);
23: M = m*bs;
24: MatCreateSeqBAIJ(PETSC_COMM_SELF,bs,M,M,1,NULL,&A);
25: MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
26: MatCreateSeqAIJ(PETSC_COMM_SELF,M,M,15,NULL,&B);
27: MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
28: PetscRandomCreate(PETSC_COMM_SELF,&rdm);
29: PetscRandomSetFromOptions(rdm);
30: VecCreateSeq(PETSC_COMM_SELF,M,&xx);
31: VecDuplicate(xx,&s1);
32: VecDuplicate(xx,&s2);
33: VecDuplicate(xx,&yy);
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(rdm,&rval);
39: col = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
40: MatSetValues(A,1,&row,1,&col,&rval,INSERT_VALUES);
41: MatSetValues(B,1,&row,1,&col,&rval,INSERT_VALUES);
42: }
43: }
45: /* Now set blocks of values */
46: for (i=0; i<20*bs; i++) {
47: PetscRandomGetValue(rdm,&rval);
48: cols[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
49: vals1[0] = rval;
50: PetscRandomGetValue(rdm,&rval);
51: cols[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
52: vals1[1] = rval;
53: PetscRandomGetValue(rdm,&rval);
54: rows[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
55: vals1[2] = rval;
56: PetscRandomGetValue(rdm,&rval);
57: rows[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
58: vals1[3] = rval;
59: MatSetValues(A,2,rows,2,cols,vals1,INSERT_VALUES);
60: MatSetValues(B,2,rows,2,cols,vals1,INSERT_VALUES);
61: }
63: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
64: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
65: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
66: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
68: /* Test MatNorm() */
69: MatNorm(A,NORM_FROBENIUS,&s1norm);
70: MatNorm(B,NORM_FROBENIUS,&s2norm);
71: rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
72: if (rnorm>tol) {
73: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);
74: }
75: MatNorm(A,NORM_INFINITY,&s1norm);
76: MatNorm(B,NORM_INFINITY,&s2norm);
77: rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
78: if (rnorm>tol) {
79: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);
80: }
81: MatNorm(A,NORM_1,&s1norm);
82: MatNorm(B,NORM_1,&s2norm);
83: rnorm = PetscAbsReal(s2norm-s1norm)/s2norm;
84: if (rnorm>tol) {
85: PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_NORM_1()- NormA=%16.14e NormB=%16.14e bs = %D\n",s1norm,s2norm,bs);
86: }
88: /* MatShift() */
89: rval = 10*s1norm;
90: MatShift(A,rval);
91: MatShift(B,rval);
93: /* Test MatTranspose() */
94: MatTranspose(A,MAT_INPLACE_MATRIX,&A);
95: MatTranspose(B,MAT_INPLACE_MATRIX,&B);
97: /* Now do MatGetValues() */
98: for (i=0; i<30; i++) {
99: PetscRandomGetValue(rdm,&rval);
100: cols[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
101: PetscRandomGetValue(rdm,&rval);
102: cols[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
103: PetscRandomGetValue(rdm,&rval);
104: rows[0] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
105: PetscRandomGetValue(rdm,&rval);
106: rows[1] = PetscMin(M-1,(PetscInt)(PetscRealPart(rval)*M));
107: MatGetValues(A,2,rows,2,cols,vals1);
108: MatGetValues(B,2,rows,2,cols,vals2);
109: PetscMemcmp(vals1,vals2,4*sizeof(PetscScalar),&flg);
110: if (!flg) {
111: PetscPrintf(PETSC_COMM_SELF,"Error: MatGetValues bs = %D\n",bs);
112: }
113: }
115: /* Test MatMult(), MatMultAdd() */
116: for (i=0; i<40; i++) {
117: VecSetRandom(xx,rdm);
118: VecSet(s2,0.0);
119: MatMult(A,xx,s1);
120: MatMultAdd(A,xx,s2,s2);
121: VecNorm(s1,NORM_2,&s1norm);
122: VecNorm(s2,NORM_2,&s2norm);
123: rnorm = s2norm-s1norm;
124: if (rnorm<-tol || rnorm>tol) {
125: PetscPrintf(PETSC_COMM_SELF,"MatMult not equal to MatMultAdd Norm1=%e Norm2=%e bs = %D\n",s1norm,s2norm,bs);
126: }
127: }
129: /* Test MatMult() */
130: MatMultEqual(A,B,10,&flg);
131: if (!flg) {
132: PetscPrintf(PETSC_COMM_SELF,"Error: MatMult()\n");
133: }
135: /* Test MatMultAdd() */
136: MatMultAddEqual(A,B,10,&flg);
137: if (!flg) {
138: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultAdd()\n");
139: }
141: /* Test MatMultTranspose() */
142: MatMultTransposeEqual(A,B,10,&flg);
143: if (!flg) {
144: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTranspose()\n");
145: }
147: /* Test MatMultTransposeAdd() */
148: MatMultTransposeAddEqual(A,B,10,&flg);
149: if (!flg) {
150: PetscPrintf(PETSC_COMM_SELF,"Error: MatMultTransposeAdd()\n");
151: }
153: /* Do LUFactor() on both the matrices */
154: PetscMalloc1(M,&idx);
155: for (i=0; i<M; i++) idx[i] = i;
156: ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is1);
157: ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is2);
158: PetscFree(idx);
159: ISSetPermutation(is1);
160: ISSetPermutation(is2);
162: MatFactorInfoInitialize(&info);
164: info.fill = 2.0;
165: info.dtcol = 0.0;
166: info.zeropivot = 1.e-14;
167: info.pivotinblocks = 1.0;
169: if (bs < 4) {
170: MatGetFactor(A,"petsc",MAT_FACTOR_LU,&Fact);
171: MatLUFactorSymbolic(Fact,A,is1,is2,&info);
172: MatLUFactorNumeric(Fact,A,&info);
173: VecSetRandom(yy,rdm);
174: MatForwardSolve(Fact,yy,xx);
175: MatBackwardSolve(Fact,xx,s1);
176: MatDestroy(&Fact);
177: VecScale(s1,-1.0);
178: MatMultAdd(A,s1,yy,yy);
179: VecNorm(yy,NORM_2,&rnorm);
180: if (rnorm<-tol || rnorm>tol) {
181: PetscPrintf(PETSC_COMM_SELF,"Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %D\n",rnorm,bs);
182: }
183: }
185: MatLUFactor(B,is1,is2,&info);
186: MatLUFactor(A,is1,is2,&info);
188: /* Test MatSolveAdd() */
189: for (i=0; i<10; i++) {
190: VecSetRandom(xx,rdm);
191: VecSetRandom(yy,rdm);
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: }
202: /* Test MatSolveAdd() when x = A'b +x */
203: for (i=0; i<10; i++) {
204: VecSetRandom(xx,rdm);
205: VecSetRandom(s1,rdm);
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: }
217: /* Test MatSolve() */
218: for (i=0; i<10; i++) {
219: VecSetRandom(xx,rdm);
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: }
230: /* Test MatSolveTranspose() */
231: if (bs < 8) {
232: for (i=0; i<10; i++) {
233: VecSetRandom(xx,rdm);
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(&rdm);
254: PetscFinalize();
255: return ierr;
256: }