xref: /petsc/src/mat/utils/multequal.c (revision cb5d8e9ef36bd5726d5c4e38d977490d990864c9)
1 
2 #include "src/mat/matimpl.h"  /*I   "petscmat.h"  I*/
3 
4 #undef __FUNCT__
5 #define __FUNCT__ "MatMultEqual"
6 /*@
7    MatMultEqual - Compares matrix-vector products of two matrices.
8 
9    Collective on Mat
10 
11    Input Parameters:
12 +  A - the first matrix
13 -  B - the second matrix
14 -  n - number of random vectors to be tested
15 
16    Output Parameter:
17 .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
18 
19    Level: intermediate
20 
21    Concepts: matrices^equality between
22 @*/
23 PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
24 {
25   PetscErrorCode ierr;
26   Vec            x,s1,s2;
27   PetscRandom    rctx;
28   PetscReal      r1,r2,tol=1.e-10;
29   PetscInt       am,an,bm,bn,k;
30 
31   PetscFunctionBegin;
32   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
33   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
34   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
35   PetscCheckSameComm(A,1,B,2);
36   ierr = PetscRandomCreate(A->comm,RANDOM_DEFAULT,&rctx);CHKERRQ(ierr);
37   ierr = VecCreate(A->comm,&x);CHKERRQ(ierr);
38   ierr = VecSetSizes(x,an,PETSC_DECIDE);CHKERRQ(ierr);
39   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
40 
41   ierr = VecCreate(A->comm,&s1);CHKERRQ(ierr);
42   ierr = VecSetSizes(s1,am,PETSC_DECIDE);CHKERRQ(ierr);
43   ierr = VecSetFromOptions(s1);CHKERRQ(ierr);
44   ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr);
45 
46   for (k=0; k<n; k++) {
47     ierr = VecSetRandom(rctx,x);CHKERRQ(ierr);
48     ierr = MatMult(A,x,s1);CHKERRQ(ierr);
49     ierr = MatMult(B,x,s2);CHKERRQ(ierr);
50     ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
51     ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
52     r1 -= r2;
53     if (r1<-tol || r1>tol) {
54       *flg = PETSC_FALSE;
55     } else {
56       *flg = PETSC_TRUE;
57     }
58   }
59   ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
60   ierr = VecDestroy(x);CHKERRQ(ierr);
61   ierr = VecDestroy(s1);CHKERRQ(ierr);
62   ierr = VecDestroy(s2);CHKERRQ(ierr);
63   PetscFunctionReturn(0);
64 }
65 
66 #undef __FUNCT__
67 #define __FUNCT__ "MatMultAddEqual"
68 /*@
69    MatMultAddEqual - Compares matrix-vector products of two matrices.
70 
71    Collective on Mat
72 
73    Input Parameters:
74 +  A - the first matrix
75 -  B - the second matrix
76 -  n - number of random vectors to be tested
77 
78    Output Parameter:
79 .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
80 
81    Level: intermediate
82 
83    Concepts: matrices^equality between
84 @*/
85 PetscErrorCode MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
86 {
87   PetscErrorCode ierr;
88   Vec            x,y,s1,s2;
89   PetscRandom    rctx;
90   PetscReal      r1,r2,tol=1.e-10;
91   PetscInt       am,an,bm,bn,k;
92 
93   PetscFunctionBegin;
94   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
95   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
96   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
97   PetscCheckSameComm(A,1,B,2);
98   ierr = PetscRandomCreate(A->comm,RANDOM_DEFAULT,&rctx);CHKERRQ(ierr);
99   ierr = VecCreate(A->comm,&x);CHKERRQ(ierr);
100   ierr = VecSetSizes(x,an,PETSC_DECIDE);CHKERRQ(ierr);
101   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
102   ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
103   ierr = VecDuplicate(x,&s1);CHKERRQ(ierr);
104   ierr = VecDuplicate(x,&s2);CHKERRQ(ierr);
105 
106   for (k=0; k<n; k++) {
107     ierr = VecSetRandom(rctx,x);CHKERRQ(ierr);
108     ierr = VecSetRandom(rctx,y);CHKERRQ(ierr);
109     ierr = MatMultAdd(A,x,y,s1);CHKERRQ(ierr);
110     ierr = MatMultAdd(B,x,y,s2);CHKERRQ(ierr);
111     ierr = VecNorm(s1,NORM_1,&r1);CHKERRQ(ierr);
112     ierr = VecNorm(s2,NORM_1,&r2);CHKERRQ(ierr);
113     r1 -= r2;
114     if (r1<-tol || r1>tol) {
115       *flg = PETSC_FALSE;
116     } else {
117       *flg = PETSC_TRUE;
118     }
119   }
120   ierr = PetscRandomDestroy(rctx);CHKERRQ(ierr);
121   ierr = VecDestroy(x);CHKERRQ(ierr);
122   ierr = VecDestroy(y);CHKERRQ(ierr);
123   ierr = VecDestroy(s1);CHKERRQ(ierr);
124   ierr = VecDestroy(s2);CHKERRQ(ierr);
125   PetscFunctionReturn(0);
126 }
127