xref: /petsc/src/mat/utils/multequal.c (revision c0aa6a63a7860d309a895cb4aa0f9e11e7859f3a)
186a22c91SHong Zhang 
2af0996ceSBarry Smith #include <petsc/private/matimpl.h>  /*I   "petscmat.h"  I*/
386a22c91SHong Zhang 
4447fed29SStefano Zampini static PetscErrorCode MatMultEqual_Private(Mat A,Mat B,PetscInt n,PetscBool *flg,PetscBool t,PetscBool add)
5447fed29SStefano Zampini {
6447fed29SStefano Zampini   PetscErrorCode ierr;
7b84f494bSStefano Zampini   Vec            Ax = NULL,Bx = NULL,s1 = NULL,s2 = NULL,Ay = NULL, By = NULL;
8447fed29SStefano Zampini   PetscRandom    rctx;
9447fed29SStefano Zampini   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
10447fed29SStefano Zampini   PetscInt       am,an,bm,bn,k;
11447fed29SStefano Zampini   PetscScalar    none = -1.0;
12447fed29SStefano Zampini   const char*    sops[] = {"MatMult","MatMultAdd","MatMultTranspose","MatMultTranposeAdd"};
13447fed29SStefano Zampini   const char*    sop;
14447fed29SStefano Zampini 
15447fed29SStefano Zampini   PetscFunctionBegin;
16447fed29SStefano Zampini   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
17447fed29SStefano Zampini   PetscValidHeaderSpecific(B,MAT_CLASSID,2);
18447fed29SStefano Zampini   PetscCheckSameComm(A,1,B,2);
19447fed29SStefano Zampini   PetscValidLogicalCollectiveInt(A,n,3);
20447fed29SStefano Zampini   PetscValidPointer(flg,4);
21447fed29SStefano Zampini   PetscValidLogicalCollectiveBool(A,t,5);
22447fed29SStefano Zampini   PetscValidLogicalCollectiveBool(A,add,6);
23447fed29SStefano Zampini   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
24447fed29SStefano Zampini   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
25*c0aa6a63SJacob Faibussowitsch   if (PetscUnlikely(am != bm || an != bn)) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT,am,bm,an,bn);
26447fed29SStefano Zampini   sop  = sops[(add ? 1 : 0) + 2 * (t ? 1 : 0)];
27447fed29SStefano Zampini   ierr = PetscRandomCreate(PetscObjectComm((PetscObject)A),&rctx);CHKERRQ(ierr);
28447fed29SStefano Zampini   ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
29447fed29SStefano Zampini   if (t) {
30447fed29SStefano Zampini     ierr = MatCreateVecs(A,&s1,&Ax);CHKERRQ(ierr);
31447fed29SStefano Zampini     ierr = MatCreateVecs(B,&s2,&Bx);CHKERRQ(ierr);
32447fed29SStefano Zampini   } else {
33447fed29SStefano Zampini     ierr = MatCreateVecs(A,&Ax,&s1);CHKERRQ(ierr);
34447fed29SStefano Zampini     ierr = MatCreateVecs(B,&Bx,&s2);CHKERRQ(ierr);
35447fed29SStefano Zampini   }
36447fed29SStefano Zampini   if (add) {
37447fed29SStefano Zampini     ierr = VecDuplicate(s1,&Ay);CHKERRQ(ierr);
38447fed29SStefano Zampini     ierr = VecDuplicate(s2,&By);CHKERRQ(ierr);
39447fed29SStefano Zampini   }
40447fed29SStefano Zampini 
41447fed29SStefano Zampini   *flg = PETSC_TRUE;
42447fed29SStefano Zampini   for (k=0; k<n; k++) {
43447fed29SStefano Zampini     ierr = VecSetRandom(Ax,rctx);CHKERRQ(ierr);
44447fed29SStefano Zampini     ierr = VecCopy(Ax,Bx);CHKERRQ(ierr);
45447fed29SStefano Zampini     if (add) {
46447fed29SStefano Zampini       ierr = VecSetRandom(Ay,rctx);CHKERRQ(ierr);
47447fed29SStefano Zampini       ierr = VecCopy(Ay,By);CHKERRQ(ierr);
48447fed29SStefano Zampini     }
49447fed29SStefano Zampini     if (t) {
50447fed29SStefano Zampini       if (add) {
51447fed29SStefano Zampini         ierr = MatMultTransposeAdd(A,Ax,Ay,s1);CHKERRQ(ierr);
52447fed29SStefano Zampini         ierr = MatMultTransposeAdd(B,Bx,By,s2);CHKERRQ(ierr);
53447fed29SStefano Zampini       } else {
54447fed29SStefano Zampini         ierr = MatMultTranspose(A,Ax,s1);CHKERRQ(ierr);
55447fed29SStefano Zampini         ierr = MatMultTranspose(B,Bx,s2);CHKERRQ(ierr);
56447fed29SStefano Zampini       }
57447fed29SStefano Zampini     } else {
58447fed29SStefano Zampini       if (add) {
59447fed29SStefano Zampini         ierr = MatMultAdd(A,Ax,Ay,s1);CHKERRQ(ierr);
60447fed29SStefano Zampini         ierr = MatMultAdd(B,Bx,By,s2);CHKERRQ(ierr);
61447fed29SStefano Zampini       } else {
62447fed29SStefano Zampini         ierr = MatMult(A,Ax,s1);CHKERRQ(ierr);
63447fed29SStefano Zampini         ierr = MatMult(B,Bx,s2);CHKERRQ(ierr);
64447fed29SStefano Zampini       }
65447fed29SStefano Zampini     }
66447fed29SStefano Zampini     ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr);
67447fed29SStefano Zampini     if (r2 < tol) {
68447fed29SStefano Zampini       ierr = VecNorm(s1,NORM_INFINITY,&r1);CHKERRQ(ierr);
69447fed29SStefano Zampini     } else {
70447fed29SStefano Zampini       ierr = VecAXPY(s2,none,s1);CHKERRQ(ierr);
71447fed29SStefano Zampini       ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr);
72447fed29SStefano Zampini       r1  /= r2;
73447fed29SStefano Zampini     }
74447fed29SStefano Zampini     if (r1 > tol) {
75447fed29SStefano Zampini       *flg = PETSC_FALSE;
76*c0aa6a63SJacob Faibussowitsch       ierr = PetscInfo3(A,"Error: %" PetscInt_FMT "-th %s() %g\n",k,sop,(double)r1);CHKERRQ(ierr);
77447fed29SStefano Zampini       break;
78447fed29SStefano Zampini     }
79447fed29SStefano Zampini   }
80447fed29SStefano Zampini   ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
81447fed29SStefano Zampini   ierr = VecDestroy(&Ax);CHKERRQ(ierr);
82447fed29SStefano Zampini   ierr = VecDestroy(&Bx);CHKERRQ(ierr);
83447fed29SStefano Zampini   ierr = VecDestroy(&Ay);CHKERRQ(ierr);
84447fed29SStefano Zampini   ierr = VecDestroy(&By);CHKERRQ(ierr);
85447fed29SStefano Zampini   ierr = VecDestroy(&s1);CHKERRQ(ierr);
86447fed29SStefano Zampini   ierr = VecDestroy(&s2);CHKERRQ(ierr);
87447fed29SStefano Zampini   PetscFunctionReturn(0);
88447fed29SStefano Zampini }
89447fed29SStefano Zampini 
90447fed29SStefano Zampini static PetscErrorCode MatMatMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg,PetscBool At,PetscBool Bt)
91447fed29SStefano Zampini {
92447fed29SStefano Zampini   PetscErrorCode ierr;
93447fed29SStefano Zampini   Vec            Ax,Bx,Cx,s1,s2,s3;
94447fed29SStefano Zampini   PetscRandom    rctx;
95447fed29SStefano Zampini   PetscReal      r1,r2,tol=PETSC_SQRT_MACHINE_EPSILON;
96447fed29SStefano Zampini   PetscInt       am,an,bm,bn,cm,cn,k;
97447fed29SStefano Zampini   PetscScalar    none = -1.0;
98447fed29SStefano Zampini   const char*    sops[] = {"MatMatMult","MatTransposeMatMult","MatMatTransposeMult","MatTransposeMatTranposeMult"};
99447fed29SStefano Zampini   const char*    sop;
100447fed29SStefano Zampini 
101447fed29SStefano Zampini   PetscFunctionBegin;
102447fed29SStefano Zampini   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
103447fed29SStefano Zampini   PetscValidHeaderSpecific(B,MAT_CLASSID,2);
104447fed29SStefano Zampini   PetscCheckSameComm(A,1,B,2);
105447fed29SStefano Zampini   PetscValidHeaderSpecific(C,MAT_CLASSID,3);
106447fed29SStefano Zampini   PetscCheckSameComm(A,1,C,3);
107447fed29SStefano Zampini   PetscValidLogicalCollectiveInt(A,n,4);
108447fed29SStefano Zampini   PetscValidPointer(flg,5);
109447fed29SStefano Zampini   PetscValidLogicalCollectiveBool(A,At,6);
110447fed29SStefano Zampini   PetscValidLogicalCollectiveBool(B,Bt,7);
111447fed29SStefano Zampini   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
112447fed29SStefano Zampini   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
113447fed29SStefano Zampini   ierr = MatGetLocalSize(C,&cm,&cn);CHKERRQ(ierr);
114447fed29SStefano Zampini   if (At) { PetscInt tt = an; an = am; am = tt; };
115447fed29SStefano Zampini   if (Bt) { PetscInt tt = bn; bn = bm; bm = tt; };
116*c0aa6a63SJacob Faibussowitsch   if (PetscUnlikely(an != bm || am != cm || bn != cn)) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A, B, C local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT,am,an,bm,bn,cm,cn);
117447fed29SStefano Zampini 
118447fed29SStefano Zampini   sop  = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
119447fed29SStefano Zampini   ierr = PetscRandomCreate(PetscObjectComm((PetscObject)C),&rctx);CHKERRQ(ierr);
120447fed29SStefano Zampini   ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
121447fed29SStefano Zampini   if (Bt) {
122447fed29SStefano Zampini     ierr = MatCreateVecs(B,&s1,&Bx);CHKERRQ(ierr);
123447fed29SStefano Zampini   } else {
124447fed29SStefano Zampini     ierr = MatCreateVecs(B,&Bx,&s1);CHKERRQ(ierr);
125447fed29SStefano Zampini   }
126447fed29SStefano Zampini   if (At) {
127447fed29SStefano Zampini     ierr = MatCreateVecs(A,&s2,&Ax);CHKERRQ(ierr);
128447fed29SStefano Zampini   } else {
129447fed29SStefano Zampini     ierr = MatCreateVecs(A,&Ax,&s2);CHKERRQ(ierr);
130447fed29SStefano Zampini   }
131447fed29SStefano Zampini   ierr = MatCreateVecs(C,&Cx,&s3);CHKERRQ(ierr);
132447fed29SStefano Zampini 
133447fed29SStefano Zampini   *flg = PETSC_TRUE;
134447fed29SStefano Zampini   for (k=0; k<n; k++) {
135447fed29SStefano Zampini     ierr = VecSetRandom(Bx,rctx);CHKERRQ(ierr);
136447fed29SStefano Zampini     if (Bt) {
137447fed29SStefano Zampini       ierr = MatMultTranspose(B,Bx,s1);CHKERRQ(ierr);
138447fed29SStefano Zampini     } else {
139447fed29SStefano Zampini       ierr = MatMult(B,Bx,s1);CHKERRQ(ierr);
140447fed29SStefano Zampini     }
141447fed29SStefano Zampini     ierr = VecCopy(s1,Ax);CHKERRQ(ierr);
142447fed29SStefano Zampini     if (At) {
143447fed29SStefano Zampini       ierr = MatMultTranspose(A,Ax,s2);CHKERRQ(ierr);
144447fed29SStefano Zampini     } else {
145447fed29SStefano Zampini       ierr = MatMult(A,Ax,s2);CHKERRQ(ierr);
146447fed29SStefano Zampini     }
147447fed29SStefano Zampini     ierr = VecCopy(Bx,Cx);CHKERRQ(ierr);
148447fed29SStefano Zampini     ierr = MatMult(C,Cx,s3);CHKERRQ(ierr);
149447fed29SStefano Zampini 
150447fed29SStefano Zampini     ierr = VecNorm(s2,NORM_INFINITY,&r2);CHKERRQ(ierr);
151447fed29SStefano Zampini     if (r2 < tol) {
152447fed29SStefano Zampini       ierr = VecNorm(s3,NORM_INFINITY,&r1);CHKERRQ(ierr);
153447fed29SStefano Zampini     } else {
154447fed29SStefano Zampini       ierr = VecAXPY(s2,none,s3);CHKERRQ(ierr);
155447fed29SStefano Zampini       ierr = VecNorm(s2,NORM_INFINITY,&r1);CHKERRQ(ierr);
156447fed29SStefano Zampini       r1  /= r2;
157447fed29SStefano Zampini     }
158447fed29SStefano Zampini     if (r1 > tol) {
159447fed29SStefano Zampini       *flg = PETSC_FALSE;
160*c0aa6a63SJacob Faibussowitsch       ierr = PetscInfo3(A,"Error: %" PetscInt_FMT "-th %s %g\n",k,sop,(double)r1);CHKERRQ(ierr);
161447fed29SStefano Zampini       break;
162447fed29SStefano Zampini     }
163447fed29SStefano Zampini   }
164447fed29SStefano Zampini   ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
165447fed29SStefano Zampini   ierr = VecDestroy(&Ax);CHKERRQ(ierr);
166447fed29SStefano Zampini   ierr = VecDestroy(&Bx);CHKERRQ(ierr);
167447fed29SStefano Zampini   ierr = VecDestroy(&Cx);CHKERRQ(ierr);
168447fed29SStefano Zampini   ierr = VecDestroy(&s1);CHKERRQ(ierr);
169447fed29SStefano Zampini   ierr = VecDestroy(&s2);CHKERRQ(ierr);
170447fed29SStefano Zampini   ierr = VecDestroy(&s3);CHKERRQ(ierr);
171447fed29SStefano Zampini   PetscFunctionReturn(0);
172447fed29SStefano Zampini }
173447fed29SStefano Zampini 
17486a22c91SHong Zhang /*@
17586a22c91SHong Zhang    MatMultEqual - Compares matrix-vector products of two matrices.
17686a22c91SHong Zhang 
17786a22c91SHong Zhang    Collective on Mat
17886a22c91SHong Zhang 
17986a22c91SHong Zhang    Input Parameters:
18086a22c91SHong Zhang +  A - the first matrix
1814222ddf1SHong Zhang .  B - the second matrix
18286a22c91SHong Zhang -  n - number of random vectors to be tested
18386a22c91SHong Zhang 
18486a22c91SHong Zhang    Output Parameter:
18586a22c91SHong Zhang .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
18686a22c91SHong Zhang 
18786a22c91SHong Zhang    Level: intermediate
18886a22c91SHong Zhang 
18986a22c91SHong Zhang @*/
1907087cfbeSBarry Smith PetscErrorCode MatMultEqual(Mat A,Mat B,PetscInt n,PetscBool *flg)
19186a22c91SHong Zhang {
19286a22c91SHong Zhang   PetscErrorCode ierr;
19386a22c91SHong Zhang 
19486a22c91SHong Zhang   PetscFunctionBegin;
195447fed29SStefano Zampini   ierr = MatMultEqual_Private(A,B,n,flg,PETSC_FALSE,PETSC_FALSE);CHKERRQ(ierr);
19686a22c91SHong Zhang   PetscFunctionReturn(0);
19786a22c91SHong Zhang }
19886a22c91SHong Zhang 
19986a22c91SHong Zhang /*@
20086a22c91SHong Zhang    MatMultAddEqual - Compares matrix-vector products of two matrices.
20186a22c91SHong Zhang 
20286a22c91SHong Zhang    Collective on Mat
20386a22c91SHong Zhang 
20486a22c91SHong Zhang    Input Parameters:
20586a22c91SHong Zhang +  A - the first matrix
2064222ddf1SHong Zhang .  B - the second matrix
20786a22c91SHong Zhang -  n - number of random vectors to be tested
20886a22c91SHong Zhang 
20986a22c91SHong Zhang    Output Parameter:
21086a22c91SHong Zhang .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
21186a22c91SHong Zhang 
21286a22c91SHong Zhang    Level: intermediate
21386a22c91SHong Zhang 
21486a22c91SHong Zhang @*/
2157087cfbeSBarry Smith PetscErrorCode  MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
21686a22c91SHong Zhang {
21786a22c91SHong Zhang   PetscErrorCode ierr;
21886a22c91SHong Zhang 
21986a22c91SHong Zhang   PetscFunctionBegin;
220447fed29SStefano Zampini   ierr = MatMultEqual_Private(A,B,n,flg,PETSC_FALSE,PETSC_TRUE);CHKERRQ(ierr);
22163879571SHong Zhang   PetscFunctionReturn(0);
22263879571SHong Zhang }
22363879571SHong Zhang 
22463879571SHong Zhang /*@
22563879571SHong Zhang    MatMultTransposeEqual - Compares matrix-vector products of two matrices.
22663879571SHong Zhang 
22763879571SHong Zhang    Collective on Mat
22863879571SHong Zhang 
22963879571SHong Zhang    Input Parameters:
23063879571SHong Zhang +  A - the first matrix
2314222ddf1SHong Zhang .  B - the second matrix
23263879571SHong Zhang -  n - number of random vectors to be tested
23363879571SHong Zhang 
23463879571SHong Zhang    Output Parameter:
23563879571SHong Zhang .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
23663879571SHong Zhang 
23763879571SHong Zhang    Level: intermediate
23863879571SHong Zhang 
23963879571SHong Zhang @*/
2407087cfbeSBarry Smith PetscErrorCode  MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
24163879571SHong Zhang {
24263879571SHong Zhang   PetscErrorCode ierr;
24363879571SHong Zhang 
24463879571SHong Zhang   PetscFunctionBegin;
245447fed29SStefano Zampini   ierr = MatMultEqual_Private(A,B,n,flg,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr);
24663879571SHong Zhang   PetscFunctionReturn(0);
24763879571SHong Zhang }
24863879571SHong Zhang 
24963879571SHong Zhang /*@
25063879571SHong Zhang    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.
25163879571SHong Zhang 
25263879571SHong Zhang    Collective on Mat
25363879571SHong Zhang 
25463879571SHong Zhang    Input Parameters:
25563879571SHong Zhang +  A - the first matrix
2564222ddf1SHong Zhang .  B - the second matrix
25763879571SHong Zhang -  n - number of random vectors to be tested
25863879571SHong Zhang 
25963879571SHong Zhang    Output Parameter:
26063879571SHong Zhang .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
26163879571SHong Zhang 
26263879571SHong Zhang    Level: intermediate
26363879571SHong Zhang 
26463879571SHong Zhang @*/
2657087cfbeSBarry Smith PetscErrorCode  MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscBool  *flg)
26663879571SHong Zhang {
26763879571SHong Zhang   PetscErrorCode ierr;
26863879571SHong Zhang 
26963879571SHong Zhang   PetscFunctionBegin;
270447fed29SStefano Zampini   ierr = MatMultEqual_Private(A,B,n,flg,PETSC_TRUE,PETSC_TRUE);CHKERRQ(ierr);
27186a22c91SHong Zhang   PetscFunctionReturn(0);
27286a22c91SHong Zhang }
273a52676ecSHong Zhang 
274a52676ecSHong Zhang /*@
275a52676ecSHong Zhang    MatMatMultEqual - Test A*B*x = C*x for n random vector x
276a52676ecSHong Zhang 
277a52676ecSHong Zhang    Collective on Mat
278a52676ecSHong Zhang 
279a52676ecSHong Zhang    Input Parameters:
280a52676ecSHong Zhang +  A - the first matrix
281c2916339SPierre Jolivet .  B - the second matrix
282c2916339SPierre Jolivet .  C - the third matrix
283a52676ecSHong Zhang -  n - number of random vectors to be tested
284a52676ecSHong Zhang 
285a52676ecSHong Zhang    Output Parameter:
286f0fc11ceSJed Brown .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
287a52676ecSHong Zhang 
288a52676ecSHong Zhang    Level: intermediate
289a52676ecSHong Zhang 
290a52676ecSHong Zhang @*/
291a52676ecSHong Zhang PetscErrorCode MatMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
292a52676ecSHong Zhang {
293a52676ecSHong Zhang   PetscErrorCode ierr;
294a52676ecSHong Zhang 
295a52676ecSHong Zhang   PetscFunctionBegin;
296447fed29SStefano Zampini   ierr = MatMatMultEqual_Private(A,B,C,n,flg,PETSC_FALSE,PETSC_FALSE);CHKERRQ(ierr);
297a52676ecSHong Zhang   PetscFunctionReturn(0);
298a52676ecSHong Zhang }
299a52676ecSHong Zhang 
300a52676ecSHong Zhang /*@
301a52676ecSHong Zhang    MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x
302a52676ecSHong Zhang 
303a52676ecSHong Zhang    Collective on Mat
304a52676ecSHong Zhang 
305a52676ecSHong Zhang    Input Parameters:
306a52676ecSHong Zhang +  A - the first matrix
3074222ddf1SHong Zhang .  B - the second matrix
3084222ddf1SHong Zhang .  C - the third matrix
309a52676ecSHong Zhang -  n - number of random vectors to be tested
310a52676ecSHong Zhang 
311a52676ecSHong Zhang    Output Parameter:
312a52676ecSHong Zhang .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
313a52676ecSHong Zhang 
314a52676ecSHong Zhang    Level: intermediate
315a52676ecSHong Zhang 
316a52676ecSHong Zhang @*/
317a52676ecSHong Zhang PetscErrorCode MatTransposeMatMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
318a52676ecSHong Zhang {
319a52676ecSHong Zhang   PetscErrorCode ierr;
320a52676ecSHong Zhang 
321a52676ecSHong Zhang   PetscFunctionBegin;
322447fed29SStefano Zampini   ierr = MatMatMultEqual_Private(A,B,C,n,flg,PETSC_TRUE,PETSC_FALSE);CHKERRQ(ierr);
323a52676ecSHong Zhang   PetscFunctionReturn(0);
324a52676ecSHong Zhang }
32586919fd6SHong Zhang 
32626546c1bSToby Isaac /*@
32726546c1bSToby Isaac    MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x
32826546c1bSToby Isaac 
32926546c1bSToby Isaac    Collective on Mat
33026546c1bSToby Isaac 
33126546c1bSToby Isaac    Input Parameters:
33226546c1bSToby Isaac +  A - the first matrix
3334222ddf1SHong Zhang .  B - the second matrix
3344222ddf1SHong Zhang .  C - the third matrix
33526546c1bSToby Isaac -  n - number of random vectors to be tested
33626546c1bSToby Isaac 
33726546c1bSToby Isaac    Output Parameter:
33826546c1bSToby Isaac .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
33926546c1bSToby Isaac 
34026546c1bSToby Isaac    Level: intermediate
34126546c1bSToby Isaac 
34226546c1bSToby Isaac @*/
343cc48ffa7SToby Isaac PetscErrorCode MatMatTransposeMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
344cc48ffa7SToby Isaac {
345cc48ffa7SToby Isaac   PetscErrorCode ierr;
346447fed29SStefano Zampini 
347447fed29SStefano Zampini   PetscFunctionBegin;
348447fed29SStefano Zampini   ierr = MatMatMultEqual_Private(A,B,C,n,flg,PETSC_FALSE,PETSC_TRUE);CHKERRQ(ierr);
349447fed29SStefano Zampini   PetscFunctionReturn(0);
350447fed29SStefano Zampini }
351447fed29SStefano Zampini 
352447fed29SStefano Zampini static PetscErrorCode MatProjMultEqual_Private(Mat A,Mat B,Mat C,PetscInt n,PetscBool rart,PetscBool *flg)
353447fed29SStefano Zampini {
354447fed29SStefano Zampini   PetscErrorCode ierr;
355447fed29SStefano Zampini   Vec            x,v1,v2,v3,v4,Cx,Bx;
356447fed29SStefano Zampini   PetscReal      norm_abs,norm_rel,tol=PETSC_SQRT_MACHINE_EPSILON;
357447fed29SStefano Zampini   PetscInt       i,am,an,bm,bn,cm,cn;
358447fed29SStefano Zampini   PetscRandom    rdm;
359cc48ffa7SToby Isaac   PetscScalar    none = -1.0;
360cc48ffa7SToby Isaac 
361cc48ffa7SToby Isaac   PetscFunctionBegin;
362cc48ffa7SToby Isaac   ierr = MatGetLocalSize(A,&am,&an);CHKERRQ(ierr);
363cc48ffa7SToby Isaac   ierr = MatGetLocalSize(B,&bm,&bn);CHKERRQ(ierr);
364447fed29SStefano Zampini   if (rart) { PetscInt t = bm; bm = bn; bn = t; }
365cc48ffa7SToby Isaac   ierr = MatGetLocalSize(C,&cm,&cn);CHKERRQ(ierr);
366*c0aa6a63SJacob Faibussowitsch   if (PetscUnlikely(an != bm || bn != cm || bn != cn)) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat A, B, C local dim %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT,am,an,bm,bn,cm,cn);
367cc48ffa7SToby Isaac 
368447fed29SStefano Zampini   /* Create left vector of A: v2 */
369447fed29SStefano Zampini   ierr = MatCreateVecs(A,&Bx,&v2);CHKERRQ(ierr);
370447fed29SStefano Zampini 
371447fed29SStefano Zampini   /* Create right vectors of B: x, v3, v4 */
372447fed29SStefano Zampini   if (rart) {
373447fed29SStefano Zampini     ierr = MatCreateVecs(B,&v1,&x);CHKERRQ(ierr);
374447fed29SStefano Zampini   } else {
375447fed29SStefano Zampini     ierr = MatCreateVecs(B,&x,&v1);CHKERRQ(ierr);
376447fed29SStefano Zampini   }
377447fed29SStefano Zampini   ierr = VecDuplicate(x,&v3);CHKERRQ(ierr);
378447fed29SStefano Zampini 
379447fed29SStefano Zampini   ierr = MatCreateVecs(C,&Cx,&v4);CHKERRQ(ierr);
380447fed29SStefano Zampini   ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rdm);CHKERRQ(ierr);
381447fed29SStefano Zampini   ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr);
382cc48ffa7SToby Isaac 
383cc48ffa7SToby Isaac   *flg = PETSC_TRUE;
384447fed29SStefano Zampini   for (i=0; i<n; i++) {
385447fed29SStefano Zampini     ierr = VecSetRandom(x,rdm);CHKERRQ(ierr);
386447fed29SStefano Zampini     ierr = VecCopy(x,Cx);CHKERRQ(ierr);
387447fed29SStefano Zampini     ierr = MatMult(C,Cx,v4);CHKERRQ(ierr);           /* v4 = C*x   */
388447fed29SStefano Zampini     if (rart) {
389447fed29SStefano Zampini       ierr = MatMultTranspose(B,x,v1);CHKERRQ(ierr);
390cc48ffa7SToby Isaac     } else {
391447fed29SStefano Zampini       ierr = MatMult(B,x,v1);CHKERRQ(ierr);
392cc48ffa7SToby Isaac     }
393447fed29SStefano Zampini     ierr = VecCopy(v1,Bx);CHKERRQ(ierr);
394447fed29SStefano Zampini     ierr = MatMult(A,Bx,v2);CHKERRQ(ierr);          /* v2 = A*B*x */
395447fed29SStefano Zampini     ierr = VecCopy(v2,v1);CHKERRQ(ierr);
396447fed29SStefano Zampini     if (rart) {
397447fed29SStefano Zampini       ierr = MatMult(B,v1,v3);CHKERRQ(ierr); /* v3 = R*A*R^t*x */
398447fed29SStefano Zampini     } else {
399447fed29SStefano Zampini       ierr = MatMultTranspose(B,v1,v3);CHKERRQ(ierr); /* v3 = Bt*A*B*x */
400447fed29SStefano Zampini     }
401447fed29SStefano Zampini     ierr = VecNorm(v4,NORM_2,&norm_abs);CHKERRQ(ierr);
402447fed29SStefano Zampini     ierr = VecAXPY(v4,none,v3);CHKERRQ(ierr);
403447fed29SStefano Zampini     ierr = VecNorm(v4,NORM_2,&norm_rel);CHKERRQ(ierr);
404447fed29SStefano Zampini 
405447fed29SStefano Zampini     if (norm_abs > tol) norm_rel /= norm_abs;
406447fed29SStefano Zampini     if (norm_rel > tol) {
407cc48ffa7SToby Isaac       *flg = PETSC_FALSE;
408*c0aa6a63SJacob Faibussowitsch       ierr = PetscInfo3(A,"Error: %" PetscInt_FMT "-th Mat%sMult() %g\n",i,rart ? "RARt" : "PtAP",(double)norm_rel);CHKERRQ(ierr);
409cc48ffa7SToby Isaac       break;
410cc48ffa7SToby Isaac     }
411cc48ffa7SToby Isaac   }
412447fed29SStefano Zampini 
413447fed29SStefano Zampini   ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr);
414cc48ffa7SToby Isaac   ierr = VecDestroy(&x);CHKERRQ(ierr);
415447fed29SStefano Zampini   ierr = VecDestroy(&Bx);CHKERRQ(ierr);
416447fed29SStefano Zampini   ierr = VecDestroy(&Cx);CHKERRQ(ierr);
417447fed29SStefano Zampini   ierr = VecDestroy(&v1);CHKERRQ(ierr);
418447fed29SStefano Zampini   ierr = VecDestroy(&v2);CHKERRQ(ierr);
419447fed29SStefano Zampini   ierr = VecDestroy(&v3);CHKERRQ(ierr);
420447fed29SStefano Zampini   ierr = VecDestroy(&v4);CHKERRQ(ierr);
421cc48ffa7SToby Isaac   PetscFunctionReturn(0);
422cc48ffa7SToby Isaac }
423cc48ffa7SToby Isaac 
42486919fd6SHong Zhang /*@
4254222ddf1SHong Zhang    MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B
4264222ddf1SHong Zhang 
4274222ddf1SHong Zhang    Collective on Mat
4284222ddf1SHong Zhang 
4294222ddf1SHong Zhang    Input Parameters:
4304222ddf1SHong Zhang +  A - the first matrix
4314222ddf1SHong Zhang .  B - the second matrix
4324222ddf1SHong Zhang .  C - the third matrix
4334222ddf1SHong Zhang -  n - number of random vectors to be tested
4344222ddf1SHong Zhang 
4354222ddf1SHong Zhang    Output Parameter:
4364222ddf1SHong Zhang .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
4374222ddf1SHong Zhang 
4384222ddf1SHong Zhang    Level: intermediate
4394222ddf1SHong Zhang 
4404222ddf1SHong Zhang @*/
4414222ddf1SHong Zhang PetscErrorCode MatPtAPMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
4424222ddf1SHong Zhang {
4434222ddf1SHong Zhang   PetscErrorCode ierr;
4444222ddf1SHong Zhang 
4454222ddf1SHong Zhang   PetscFunctionBegin;
446447fed29SStefano Zampini   ierr = MatProjMultEqual_Private(A,B,C,n,PETSC_FALSE,flg);CHKERRQ(ierr);
447447fed29SStefano Zampini   PetscFunctionReturn(0);
4484222ddf1SHong Zhang }
4494222ddf1SHong Zhang 
450447fed29SStefano Zampini /*@
451447fed29SStefano Zampini    MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t
452447fed29SStefano Zampini 
453447fed29SStefano Zampini    Collective on Mat
454447fed29SStefano Zampini 
455447fed29SStefano Zampini    Input Parameters:
456447fed29SStefano Zampini +  A - the first matrix
457447fed29SStefano Zampini .  B - the second matrix
458447fed29SStefano Zampini .  C - the third matrix
459447fed29SStefano Zampini -  n - number of random vectors to be tested
460447fed29SStefano Zampini 
461447fed29SStefano Zampini    Output Parameter:
462447fed29SStefano Zampini .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.
463447fed29SStefano Zampini 
464447fed29SStefano Zampini    Level: intermediate
465447fed29SStefano Zampini 
466447fed29SStefano Zampini @*/
467447fed29SStefano Zampini PetscErrorCode MatRARtMultEqual(Mat A,Mat B,Mat C,PetscInt n,PetscBool *flg)
468447fed29SStefano Zampini {
469447fed29SStefano Zampini   PetscErrorCode ierr;
470447fed29SStefano Zampini 
471447fed29SStefano Zampini   PetscFunctionBegin;
472447fed29SStefano Zampini   ierr = MatProjMultEqual_Private(A,B,C,n,PETSC_TRUE,flg);CHKERRQ(ierr);
4734222ddf1SHong Zhang   PetscFunctionReturn(0);
4744222ddf1SHong Zhang }
4754222ddf1SHong Zhang 
4764222ddf1SHong Zhang /*@
47786919fd6SHong Zhang    MatIsLinear - Check if a shell matrix A is a linear operator.
47886919fd6SHong Zhang 
47986919fd6SHong Zhang    Collective on Mat
48086919fd6SHong Zhang 
48186919fd6SHong Zhang    Input Parameters:
48286919fd6SHong Zhang +  A - the shell matrix
48386919fd6SHong Zhang -  n - number of random vectors to be tested
48486919fd6SHong Zhang 
48586919fd6SHong Zhang    Output Parameter:
48686919fd6SHong Zhang .  flg - PETSC_TRUE if the shell matrix is linear; PETSC_FALSE otherwise.
48786919fd6SHong Zhang 
48886919fd6SHong Zhang    Level: intermediate
48986919fd6SHong Zhang @*/
49086919fd6SHong Zhang PetscErrorCode MatIsLinear(Mat A,PetscInt n,PetscBool  *flg)
49186919fd6SHong Zhang {
49286919fd6SHong Zhang   PetscErrorCode ierr;
49386919fd6SHong Zhang   Vec            x,y,s1,s2;
49486919fd6SHong Zhang   PetscRandom    rctx;
49586919fd6SHong Zhang   PetscScalar    a;
49686919fd6SHong Zhang   PetscInt       k;
49786919fd6SHong Zhang   PetscReal      norm,normA;
49886919fd6SHong Zhang   MPI_Comm       comm;
49986919fd6SHong Zhang   PetscMPIInt    rank;
50086919fd6SHong Zhang 
50186919fd6SHong Zhang   PetscFunctionBegin;
50286919fd6SHong Zhang   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
50386919fd6SHong Zhang   ierr = PetscObjectGetComm((PetscObject)A,&comm);CHKERRQ(ierr);
504ffc4695bSBarry Smith   ierr = MPI_Comm_rank(comm,&rank);CHKERRMPI(ierr);
50586919fd6SHong Zhang 
50686919fd6SHong Zhang   ierr = PetscRandomCreate(comm,&rctx);CHKERRQ(ierr);
50786919fd6SHong Zhang   ierr = PetscRandomSetFromOptions(rctx);CHKERRQ(ierr);
50886919fd6SHong Zhang   ierr = MatCreateVecs(A,&x,&s1);CHKERRQ(ierr);
50986919fd6SHong Zhang   ierr = VecDuplicate(x,&y);CHKERRQ(ierr);
51086919fd6SHong Zhang   ierr = VecDuplicate(s1,&s2);CHKERRQ(ierr);
51186919fd6SHong Zhang 
51286919fd6SHong Zhang   *flg = PETSC_TRUE;
51386919fd6SHong Zhang   for (k=0; k<n; k++) {
51486919fd6SHong Zhang     ierr = VecSetRandom(x,rctx);CHKERRQ(ierr);
5156d5db48cSHong Zhang     ierr = VecSetRandom(y,rctx);CHKERRQ(ierr);
516dd400576SPatrick Sanan     if (rank == 0) {
51786919fd6SHong Zhang       ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr);
51886919fd6SHong Zhang     }
519ffc4695bSBarry Smith     ierr = MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm);CHKERRMPI(ierr);
52086919fd6SHong Zhang 
52186919fd6SHong Zhang     /* s2 = a*A*x + A*y */
52286919fd6SHong Zhang     ierr = MatMult(A,y,s2);CHKERRQ(ierr); /* s2 = A*y */
52386919fd6SHong Zhang     ierr = MatMult(A,x,s1);CHKERRQ(ierr); /* s1 = A*x */
52486919fd6SHong Zhang     ierr = VecAXPY(s2,a,s1);CHKERRQ(ierr); /* s2 = a s1 + s2 */
52586919fd6SHong Zhang 
52686919fd6SHong Zhang     /* s1 = A * (a x + y) */
52786919fd6SHong Zhang     ierr = VecAXPY(y,a,x);CHKERRQ(ierr); /* y = a x + y */
52886919fd6SHong Zhang     ierr = MatMult(A,y,s1);CHKERRQ(ierr);
52986919fd6SHong Zhang     ierr = VecNorm(s1,NORM_INFINITY,&normA);CHKERRQ(ierr);
53086919fd6SHong Zhang 
53186919fd6SHong Zhang     ierr = VecAXPY(s2,-1.0,s1);CHKERRQ(ierr); /* s2 = - s1 + s2 */
53286919fd6SHong Zhang     ierr = VecNorm(s2,NORM_INFINITY,&norm);CHKERRQ(ierr);
5331b8dac88SHong Zhang     if (norm/normA > 100.*PETSC_MACHINE_EPSILON) {
53486919fd6SHong Zhang       *flg = PETSC_FALSE;
535*c0aa6a63SJacob Faibussowitsch       ierr = PetscInfo3(A,"Error: %" PetscInt_FMT "-th |A*(ax+y) - (a*A*x+A*y)|/|A(ax+y)| %g > tol %g\n",k,(double)(norm/normA),(double)(100.*PETSC_MACHINE_EPSILON));CHKERRQ(ierr);
53686919fd6SHong Zhang       break;
53786919fd6SHong Zhang     }
53886919fd6SHong Zhang   }
53986919fd6SHong Zhang   ierr = PetscRandomDestroy(&rctx);CHKERRQ(ierr);
54086919fd6SHong Zhang   ierr = VecDestroy(&x);CHKERRQ(ierr);
54186919fd6SHong Zhang   ierr = VecDestroy(&y);CHKERRQ(ierr);
54286919fd6SHong Zhang   ierr = VecDestroy(&s1);CHKERRQ(ierr);
54386919fd6SHong Zhang   ierr = VecDestroy(&s2);CHKERRQ(ierr);
54486919fd6SHong Zhang   PetscFunctionReturn(0);
54586919fd6SHong Zhang }
546