Actual source code: ex48.c
petsc-3.14.6 2021-03-30
2: static char help[] = "Tests various routines in MatSeqBAIJ format.\n";
4: #include <petscmat.h>
6: int main(int argc,char **args)
7: {
8: Mat A,B,C,D,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: PetscArraycmp(vals1,vals2,4,&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: /* Test MatMatMult() */
154: MatCreateSeqDense(PETSC_COMM_SELF,M,40,NULL,&C);
155: MatSetRandom(C,rdm);
156: MatMatMult(A,C,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&D);
157: MatMatMultEqual(A,C,D,40,&flg);
158: if (!flg) {
159: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n");
160: }
161: MatDestroy(&D);
162: MatCreateSeqDense(PETSC_COMM_SELF,M,40,NULL,&D);
163: MatMatMult(A,C,MAT_REUSE_MATRIX,PETSC_DEFAULT,&D);
164: MatMatMultEqual(A,C,D,40,&flg);
165: if (!flg) {
166: PetscPrintf(PETSC_COMM_SELF,"Error: MatMatMult()\n");
167: }
169: /* Do LUFactor() on both the matrices */
170: PetscMalloc1(M,&idx);
171: for (i=0; i<M; i++) idx[i] = i;
172: ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is1);
173: ISCreateGeneral(PETSC_COMM_SELF,M,idx,PETSC_COPY_VALUES,&is2);
174: PetscFree(idx);
175: ISSetPermutation(is1);
176: ISSetPermutation(is2);
178: MatFactorInfoInitialize(&info);
180: info.fill = 2.0;
181: info.dtcol = 0.0;
182: info.zeropivot = 1.e-14;
183: info.pivotinblocks = 1.0;
185: if (bs < 4) {
186: MatGetFactor(A,"petsc",MAT_FACTOR_LU,&Fact);
187: MatLUFactorSymbolic(Fact,A,is1,is2,&info);
188: MatLUFactorNumeric(Fact,A,&info);
189: VecSetRandom(yy,rdm);
190: MatForwardSolve(Fact,yy,xx);
191: MatBackwardSolve(Fact,xx,s1);
192: MatDestroy(&Fact);
193: VecScale(s1,-1.0);
194: MatMultAdd(A,s1,yy,yy);
195: VecNorm(yy,NORM_2,&rnorm);
196: if (rnorm<-tol || rnorm>tol) {
197: PetscPrintf(PETSC_COMM_SELF,"Error:MatForwardSolve/MatBackwardSolve - Norm1=%16.14e bs = %D\n",rnorm,bs);
198: }
199: }
201: MatLUFactor(B,is1,is2,&info);
202: MatLUFactor(A,is1,is2,&info);
204: /* Test MatSolveAdd() */
205: for (i=0; i<10; i++) {
206: VecSetRandom(xx,rdm);
207: VecSetRandom(yy,rdm);
208: MatSolveAdd(B,xx,yy,s2);
209: MatSolveAdd(A,xx,yy,s1);
210: VecNorm(s1,NORM_2,&s1norm);
211: VecNorm(s2,NORM_2,&s2norm);
212: rnorm = s2norm-s1norm;
213: if (rnorm<-tol || rnorm>tol) {
214: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
215: }
216: }
218: /* Test MatSolveAdd() when x = A'b +x */
219: for (i=0; i<10; i++) {
220: VecSetRandom(xx,rdm);
221: VecSetRandom(s1,rdm);
222: VecCopy(s2,s1);
223: MatSolveAdd(B,xx,s2,s2);
224: MatSolveAdd(A,xx,s1,s1);
225: VecNorm(s1,NORM_2,&s1norm);
226: VecNorm(s2,NORM_2,&s2norm);
227: rnorm = s2norm-s1norm;
228: if (rnorm<-tol || rnorm>tol) {
229: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveAdd(same) - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
230: }
231: }
233: /* Test MatSolve() */
234: for (i=0; i<10; i++) {
235: VecSetRandom(xx,rdm);
236: MatSolve(B,xx,s2);
237: MatSolve(A,xx,s1);
238: VecNorm(s1,NORM_2,&s1norm);
239: VecNorm(s2,NORM_2,&s2norm);
240: rnorm = s2norm-s1norm;
241: if (rnorm<-tol || rnorm>tol) {
242: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
243: }
244: }
246: /* Test MatSolveTranspose() */
247: if (bs < 8) {
248: for (i=0; i<10; i++) {
249: VecSetRandom(xx,rdm);
250: MatSolveTranspose(B,xx,s2);
251: MatSolveTranspose(A,xx,s1);
252: VecNorm(s1,NORM_2,&s1norm);
253: VecNorm(s2,NORM_2,&s2norm);
254: rnorm = s2norm-s1norm;
255: if (rnorm<-tol || rnorm>tol) {
256: PetscPrintf(PETSC_COMM_SELF,"Error:MatSolveTranspose - Norm1=%16.14e Norm2=%16.14e bs = %D\n",s1norm,s2norm,bs);
257: }
258: }
259: }
261: MatDestroy(&A);
262: MatDestroy(&B);
263: MatDestroy(&C);
264: MatDestroy(&D);
265: VecDestroy(&xx);
266: VecDestroy(&s1);
267: VecDestroy(&s2);
268: VecDestroy(&yy);
269: ISDestroy(&is1);
270: ISDestroy(&is2);
271: PetscRandomDestroy(&rdm);
272: PetscFinalize();
273: return ierr;
274: }
277: /*TEST
279: test:
280: args: -mat_block_size {{1 2 3 4 5 6 7 8}}
282: TEST*/