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