Actual source code: multequal.c

      1: #include <petsc/private/matimpl.h>

      3: /*
      4:   n; try the MatMult variant n times
      5:   flg: return the boolean result, equal or not
      6:   t: 0 => no transpose; 1 => transpose; 2 => Hermitian transpose
      7:   add:  0 => no add (e.g., y = Ax);  1 => add third vector (e.g., z = Ax + y); 2 => add update (e.g., y = Ax + y)
      8: */
      9: static PetscErrorCode MatMultEqual_Private(Mat A, Mat B, PetscInt n, PetscBool *flg, PetscInt t, PetscInt add)
     10: {
     11:   Vec         Ax = NULL, Bx = NULL, s1 = NULL, s2 = NULL, Ay = NULL, By = NULL;
     12:   PetscRandom rctx;
     13:   PetscReal   r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
     14:   PetscInt    am, an, bm, bn, k;
     15:   PetscScalar none = -1.0;
     16: #if defined(PETSC_USE_INFO)
     17:   const char *sops[] = {"MatMult", "MatMultAdd", "MatMultAdd (update)", "MatMultTranspose", "MatMultTransposeAdd", "MatMultTransposeAdd (update)", "MatMultHermitianTranspose", "MatMultHermitianTransposeAdd", "MatMultHermitianTransposeAdd (update)"};
     18:   const char *sop;
     19: #endif

     21:   PetscFunctionBegin;
     24:   PetscCheckSameComm(A, 1, B, 2);
     26:   PetscAssertPointer(flg, 4);
     29:   PetscCall(MatGetLocalSize(A, &am, &an));
     30:   PetscCall(MatGetLocalSize(B, &bm, &bn));
     31:   PetscCheck(am == bm && an == bn, 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);
     32: #if defined(PETSC_USE_INFO)
     33:   sop = sops[add + 3 * t];
     34: #endif
     35:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)A), &rctx));
     36:   PetscCall(PetscRandomSetFromOptions(rctx));
     37:   if (t) {
     38:     PetscCall(MatCreateVecs(A, &s1, &Ax));
     39:     PetscCall(MatCreateVecs(B, &s2, &Bx));
     40:   } else {
     41:     PetscCall(MatCreateVecs(A, &Ax, &s1));
     42:     PetscCall(MatCreateVecs(B, &Bx, &s2));
     43:   }
     44:   if (add) {
     45:     PetscCall(VecDuplicate(s1, &Ay));
     46:     PetscCall(VecDuplicate(s2, &By));
     47:   }

     49:   *flg = PETSC_TRUE;
     50:   for (k = 0; k < n; k++) {
     51:     Vec Aadd = NULL, Badd = NULL;

     53:     PetscCall(VecSetRandom(Ax, rctx));
     54:     PetscCall(VecCopy(Ax, Bx));
     55:     if (add) {
     56:       PetscCall(VecSetRandom(Ay, rctx));
     57:       PetscCall(VecCopy(Ay, By));
     58:       Aadd = Ay;
     59:       Badd = By;
     60:       if (add == 2) {
     61:         PetscCall(VecCopy(Ay, s1));
     62:         PetscCall(VecCopy(By, s2));
     63:         Aadd = s1;
     64:         Badd = s2;
     65:       }
     66:     }
     67:     if (t == 1) {
     68:       if (add) {
     69:         PetscCall(MatMultTransposeAdd(A, Ax, Aadd, s1));
     70:         PetscCall(MatMultTransposeAdd(B, Bx, Badd, s2));
     71:       } else {
     72:         PetscCall(MatMultTranspose(A, Ax, s1));
     73:         PetscCall(MatMultTranspose(B, Bx, s2));
     74:       }
     75:     } else if (t == 2) {
     76:       if (add) {
     77:         PetscCall(MatMultHermitianTransposeAdd(A, Ax, Aadd, s1));
     78:         PetscCall(MatMultHermitianTransposeAdd(B, Bx, Badd, s2));
     79:       } else {
     80:         PetscCall(MatMultHermitianTranspose(A, Ax, s1));
     81:         PetscCall(MatMultHermitianTranspose(B, Bx, s2));
     82:       }
     83:     } else {
     84:       if (add) {
     85:         PetscCall(MatMultAdd(A, Ax, Aadd, s1));
     86:         PetscCall(MatMultAdd(B, Bx, Badd, s2));
     87:       } else {
     88:         PetscCall(MatMult(A, Ax, s1));
     89:         PetscCall(MatMult(B, Bx, s2));
     90:       }
     91:     }
     92:     PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
     93:     if (r2 < tol) {
     94:       PetscCall(VecNorm(s1, NORM_INFINITY, &r1));
     95:     } else {
     96:       PetscCall(VecAXPY(s2, none, s1));
     97:       PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
     98:       r1 /= r2;
     99:     }
    100:     if (r1 > tol) {
    101:       *flg = PETSC_FALSE;
    102:       PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s() %g\n", k, sop, (double)r1));
    103:       break;
    104:     }
    105:   }
    106:   PetscCall(PetscRandomDestroy(&rctx));
    107:   PetscCall(VecDestroy(&Ax));
    108:   PetscCall(VecDestroy(&Bx));
    109:   PetscCall(VecDestroy(&Ay));
    110:   PetscCall(VecDestroy(&By));
    111:   PetscCall(VecDestroy(&s1));
    112:   PetscCall(VecDestroy(&s2));
    113:   PetscFunctionReturn(PETSC_SUCCESS);
    114: }

    116: static PetscErrorCode MatMatMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg, PetscBool At, PetscBool Bt)
    117: {
    118:   Vec         Ax, Bx, Cx, s1, s2, s3;
    119:   PetscRandom rctx;
    120:   PetscReal   r1, r2, tol = PETSC_SQRT_MACHINE_EPSILON;
    121:   PetscInt    am, an, bm, bn, cm, cn, k;
    122:   PetscScalar none = -1.0;
    123: #if defined(PETSC_USE_INFO)
    124:   const char *sops[] = {"MatMatMult", "MatTransposeMatMult", "MatMatTransposeMult", "MatTransposeMatTransposeMult"};
    125:   const char *sop;
    126: #endif

    128:   PetscFunctionBegin;
    131:   PetscCheckSameComm(A, 1, B, 2);
    133:   PetscCheckSameComm(A, 1, C, 3);
    135:   PetscAssertPointer(flg, 5);
    138:   PetscCall(MatGetLocalSize(A, &am, &an));
    139:   PetscCall(MatGetLocalSize(B, &bm, &bn));
    140:   PetscCall(MatGetLocalSize(C, &cm, &cn));
    141:   if (At) {
    142:     PetscInt tt = an;
    143:     an          = am;
    144:     am          = tt;
    145:   }
    146:   if (Bt) {
    147:     PetscInt tt = bn;
    148:     bn          = bm;
    149:     bm          = tt;
    150:   }
    151:   PetscCheck(an == bm && am == cm && bn == cn, 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);

    153: #if defined(PETSC_USE_INFO)
    154:   sop = sops[(At ? 1 : 0) + 2 * (Bt ? 1 : 0)];
    155: #endif
    156:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)C), &rctx));
    157:   PetscCall(PetscRandomSetFromOptions(rctx));
    158:   if (Bt) {
    159:     PetscCall(MatCreateVecs(B, &s1, &Bx));
    160:   } else {
    161:     PetscCall(MatCreateVecs(B, &Bx, &s1));
    162:   }
    163:   if (At) {
    164:     PetscCall(MatCreateVecs(A, &s2, &Ax));
    165:   } else {
    166:     PetscCall(MatCreateVecs(A, &Ax, &s2));
    167:   }
    168:   PetscCall(MatCreateVecs(C, &Cx, &s3));

    170:   *flg = PETSC_TRUE;
    171:   for (k = 0; k < n; k++) {
    172:     PetscCall(VecSetRandom(Bx, rctx));
    173:     if (Bt) {
    174:       PetscCall(MatMultTranspose(B, Bx, s1));
    175:     } else {
    176:       PetscCall(MatMult(B, Bx, s1));
    177:     }
    178:     PetscCall(VecCopy(s1, Ax));
    179:     if (At) {
    180:       PetscCall(MatMultTranspose(A, Ax, s2));
    181:     } else {
    182:       PetscCall(MatMult(A, Ax, s2));
    183:     }
    184:     PetscCall(VecCopy(Bx, Cx));
    185:     PetscCall(MatMult(C, Cx, s3));

    187:     PetscCall(VecNorm(s2, NORM_INFINITY, &r2));
    188:     if (r2 < tol) {
    189:       PetscCall(VecNorm(s3, NORM_INFINITY, &r1));
    190:     } else {
    191:       PetscCall(VecAXPY(s2, none, s3));
    192:       PetscCall(VecNorm(s2, NORM_INFINITY, &r1));
    193:       r1 /= r2;
    194:     }
    195:     if (r1 > tol) {
    196:       *flg = PETSC_FALSE;
    197:       PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th %s %g\n", k, sop, (double)r1));
    198:       break;
    199:     }
    200:   }
    201:   PetscCall(PetscRandomDestroy(&rctx));
    202:   PetscCall(VecDestroy(&Ax));
    203:   PetscCall(VecDestroy(&Bx));
    204:   PetscCall(VecDestroy(&Cx));
    205:   PetscCall(VecDestroy(&s1));
    206:   PetscCall(VecDestroy(&s2));
    207:   PetscCall(VecDestroy(&s3));
    208:   PetscFunctionReturn(PETSC_SUCCESS);
    209: }

    211: /*@
    212:   MatMultEqual - Compares matrix-vector products of two matrices using `n` random vectors

    214:   Collective

    216:   Input Parameters:
    217: + A - the first matrix
    218: . B - the second matrix
    219: - n - number of random vectors to be tested

    221:   Output Parameter:
    222: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    224:   Level: intermediate

    226:   Note:
    227:   The tolerance for equality is a generous `PETSC_SQRT_MACHINE_EPSILON` in the norm of the difference of the two computed vectors to
    228:   allow for differences in the numerical computations. Hence this routine may indicate equality even if there is a small systematic difference
    229:   between the two matrices.

    231: .seealso: `Mat`, `MatMultAddEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`, `MatIsLinear()`, `MatEqual()`
    232: @*/
    233: PetscErrorCode MatMultEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
    234: {
    235:   PetscFunctionBegin;
    236:   PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 0));
    237:   PetscFunctionReturn(PETSC_SUCCESS);
    238: }

    240: /*@
    241:   MatMultAddEqual - Compares matrix-vector product plus vector add of two matrices.

    243:   Collective

    245:   Input Parameters:
    246: + A - the first matrix
    247: . B - the second matrix
    248: - n - number of random vectors to be tested

    250:   Output Parameter:
    251: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    253:   Level: intermediate

    255: .seealso: `Mat`, `MatMultEqual()`, `MatMultTransposeEqual()`, `MatMultTransposeAddEqual()`
    256: @*/
    257: PetscErrorCode MatMultAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
    258: {
    259:   PetscFunctionBegin;
    260:   PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 1));
    261:   PetscCall(MatMultEqual_Private(A, B, n, flg, 0, 2));
    262:   PetscFunctionReturn(PETSC_SUCCESS);
    263: }

    265: /*@
    266:   MatMultTransposeEqual - Compares matrix-vector products of two matrices.

    268:   Collective

    270:   Input Parameters:
    271: + A - the first matrix
    272: . B - the second matrix
    273: - n - number of random vectors to be tested

    275:   Output Parameter:
    276: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    278:   Level: intermediate

    280: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeAddEqual()`
    281: @*/
    282: PetscErrorCode MatMultTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
    283: {
    284:   PetscFunctionBegin;
    285:   PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 0));
    286:   PetscFunctionReturn(PETSC_SUCCESS);
    287: }

    289: /*@
    290:   MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

    292:   Collective

    294:   Input Parameters:
    295: + A - the first matrix
    296: . B - the second matrix
    297: - n - number of random vectors to be tested

    299:   Output Parameter:
    300: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    302:   Level: intermediate

    304: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    305: @*/
    306: PetscErrorCode MatMultTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
    307: {
    308:   PetscFunctionBegin;
    309:   PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 1));
    310:   PetscCall(MatMultEqual_Private(A, B, n, flg, 1, 2));
    311:   PetscFunctionReturn(PETSC_SUCCESS);
    312: }

    314: /*@
    315:   MatMultHermitianTransposeEqual - Compares matrix-vector products of two matrices.

    317:   Collective

    319:   Input Parameters:
    320: + A - the first matrix
    321: . B - the second matrix
    322: - n - number of random vectors to be tested

    324:   Output Parameter:
    325: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    327:   Level: intermediate

    329: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    330: @*/
    331: PetscErrorCode MatMultHermitianTransposeEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
    332: {
    333:   PetscFunctionBegin;
    334:   PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 0));
    335:   PetscFunctionReturn(PETSC_SUCCESS);
    336: }

    338: /*@
    339:   MatMultHermitianTransposeAddEqual - Compares matrix-vector products of two matrices.

    341:   Collective

    343:   Input Parameters:
    344: + A - the first matrix
    345: . B - the second matrix
    346: - n - number of random vectors to be tested

    348:   Output Parameter:
    349: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    351:   Level: intermediate

    353: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    354: @*/
    355: PetscErrorCode MatMultHermitianTransposeAddEqual(Mat A, Mat B, PetscInt n, PetscBool *flg)
    356: {
    357:   PetscFunctionBegin;
    358:   PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 1));
    359:   PetscCall(MatMultEqual_Private(A, B, n, flg, 2, 2));
    360:   PetscFunctionReturn(PETSC_SUCCESS);
    361: }

    363: /*@
    364:   MatMatMultEqual - Test A*B*x = C*x for n random vector x

    366:   Collective

    368:   Input Parameters:
    369: + A - the first matrix
    370: . B - the second matrix
    371: . C - the third matrix
    372: - n - number of random vectors to be tested

    374:   Output Parameter:
    375: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    377:   Level: intermediate

    379: .seealso: `Mat`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    380: @*/
    381: PetscErrorCode MatMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
    382: {
    383:   PetscFunctionBegin;
    384:   PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_FALSE));
    385:   PetscFunctionReturn(PETSC_SUCCESS);
    386: }

    388: /*@
    389:   MatTransposeMatMultEqual - Test A^T*B*x = C*x for n random vector x

    391:   Collective

    393:   Input Parameters:
    394: + A - the first matrix
    395: . B - the second matrix
    396: . C - the third matrix
    397: - n - number of random vectors to be tested

    399:   Output Parameter:
    400: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    402:   Level: intermediate

    404: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    405: @*/
    406: PetscErrorCode MatTransposeMatMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
    407: {
    408:   PetscFunctionBegin;
    409:   PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_TRUE, PETSC_FALSE));
    410:   PetscFunctionReturn(PETSC_SUCCESS);
    411: }

    413: /*@
    414:   MatMatTransposeMultEqual - Test A*B^T*x = C*x for n random vector x

    416:   Collective

    418:   Input Parameters:
    419: + A - the first matrix
    420: . B - the second matrix
    421: . C - the third matrix
    422: - n - number of random vectors to be tested

    424:   Output Parameter:
    425: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    427:   Level: intermediate

    429: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    430: @*/
    431: PetscErrorCode MatMatTransposeMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
    432: {
    433:   PetscFunctionBegin;
    434:   PetscCall(MatMatMultEqual_Private(A, B, C, n, flg, PETSC_FALSE, PETSC_TRUE));
    435:   PetscFunctionReturn(PETSC_SUCCESS);
    436: }

    438: static PetscErrorCode MatProjMultEqual_Private(Mat A, Mat B, Mat C, PetscInt n, PetscBool rart, PetscBool *flg)
    439: {
    440:   Vec         x, v1, v2, v3, v4, Cx, Bx;
    441:   PetscReal   norm_abs, norm_rel, tol = PETSC_SQRT_MACHINE_EPSILON;
    442:   PetscInt    i, am, an, bm, bn, cm, cn;
    443:   PetscRandom rdm;
    444:   PetscScalar none = -1.0;

    446:   PetscFunctionBegin;
    447:   PetscCall(MatGetLocalSize(A, &am, &an));
    448:   PetscCall(MatGetLocalSize(B, &bm, &bn));
    449:   if (rart) {
    450:     PetscInt t = bm;
    451:     bm         = bn;
    452:     bn         = t;
    453:   }
    454:   PetscCall(MatGetLocalSize(C, &cm, &cn));
    455:   PetscCheck(an == bm && bn == cm && bn == cn, 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);

    457:   /* Create left vector of A: v2 */
    458:   PetscCall(MatCreateVecs(A, &Bx, &v2));

    460:   /* Create right vectors of B: x, v3, v4 */
    461:   if (rart) {
    462:     PetscCall(MatCreateVecs(B, &v1, &x));
    463:   } else {
    464:     PetscCall(MatCreateVecs(B, &x, &v1));
    465:   }
    466:   PetscCall(VecDuplicate(x, &v3));

    468:   PetscCall(MatCreateVecs(C, &Cx, &v4));
    469:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rdm));
    470:   PetscCall(PetscRandomSetFromOptions(rdm));

    472:   *flg = PETSC_TRUE;
    473:   for (i = 0; i < n; i++) {
    474:     PetscCall(VecSetRandom(x, rdm));
    475:     PetscCall(VecCopy(x, Cx));
    476:     PetscCall(MatMult(C, Cx, v4)); /* v4 = C*x   */
    477:     if (rart) {
    478:       PetscCall(MatMultTranspose(B, x, v1));
    479:     } else {
    480:       PetscCall(MatMult(B, x, v1));
    481:     }
    482:     PetscCall(VecCopy(v1, Bx));
    483:     PetscCall(MatMult(A, Bx, v2)); /* v2 = A*B*x */
    484:     PetscCall(VecCopy(v2, v1));
    485:     if (rart) {
    486:       PetscCall(MatMult(B, v1, v3)); /* v3 = R*A*R^t*x */
    487:     } else {
    488:       PetscCall(MatMultTranspose(B, v1, v3)); /* v3 = Bt*A*B*x */
    489:     }
    490:     PetscCall(VecNorm(v4, NORM_2, &norm_abs));
    491:     PetscCall(VecAXPY(v4, none, v3));
    492:     PetscCall(VecNorm(v4, NORM_2, &norm_rel));

    494:     if (norm_abs > tol) norm_rel /= norm_abs;
    495:     if (norm_rel > tol) {
    496:       *flg = PETSC_FALSE;
    497:       PetscCall(PetscInfo(A, "Error: %" PetscInt_FMT "-th Mat%sMult() %g\n", i, rart ? "RARt" : "PtAP", (double)norm_rel));
    498:       break;
    499:     }
    500:   }

    502:   PetscCall(PetscRandomDestroy(&rdm));
    503:   PetscCall(VecDestroy(&x));
    504:   PetscCall(VecDestroy(&Bx));
    505:   PetscCall(VecDestroy(&Cx));
    506:   PetscCall(VecDestroy(&v1));
    507:   PetscCall(VecDestroy(&v2));
    508:   PetscCall(VecDestroy(&v3));
    509:   PetscCall(VecDestroy(&v4));
    510:   PetscFunctionReturn(PETSC_SUCCESS);
    511: }

    513: /*@
    514:   MatPtAPMultEqual - Compares matrix-vector products of C = Bt*A*B

    516:   Collective

    518:   Input Parameters:
    519: + A - the first matrix
    520: . B - the second matrix
    521: . C - the third matrix
    522: - n - number of random vectors to be tested

    524:   Output Parameter:
    525: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    527:   Level: intermediate

    529: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    530: @*/
    531: PetscErrorCode MatPtAPMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
    532: {
    533:   PetscFunctionBegin;
    534:   PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_FALSE, flg));
    535:   PetscFunctionReturn(PETSC_SUCCESS);
    536: }

    538: /*@
    539:   MatRARtMultEqual - Compares matrix-vector products of C = B*A*B^t

    541:   Collective

    543:   Input Parameters:
    544: + A - the first matrix
    545: . B - the second matrix
    546: . C - the third matrix
    547: - n - number of random vectors to be tested

    549:   Output Parameter:
    550: . flg - `PETSC_TRUE` if the products are equal; `PETSC_FALSE` otherwise.

    552:   Level: intermediate

    554: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    555: @*/
    556: PetscErrorCode MatRARtMultEqual(Mat A, Mat B, Mat C, PetscInt n, PetscBool *flg)
    557: {
    558:   PetscFunctionBegin;
    559:   PetscCall(MatProjMultEqual_Private(A, B, C, n, PETSC_TRUE, flg));
    560:   PetscFunctionReturn(PETSC_SUCCESS);
    561: }

    563: /*@
    564:   MatIsLinear - Check if a shell matrix `A` is a linear operator.

    566:   Collective

    568:   Input Parameters:
    569: + A - the shell matrix
    570: - n - number of random vectors to be tested

    572:   Output Parameter:
    573: . flg - `PETSC_TRUE` if the shell matrix is linear; `PETSC_FALSE` otherwise.

    575:   Level: intermediate

    577: .seealso: `Mat`, `MatMatMultEqual()`, `MatMultEqual()`, `MatMultAddEqual()`, `MatMultTransposeEqual()`
    578: @*/
    579: PetscErrorCode MatIsLinear(Mat A, PetscInt n, PetscBool *flg)
    580: {
    581:   Vec         x, y, s1, s2;
    582:   PetscRandom rctx;
    583:   PetscScalar a;
    584:   PetscInt    k;
    585:   PetscReal   norm, normA;
    586:   MPI_Comm    comm;
    587:   PetscMPIInt rank;

    589:   PetscFunctionBegin;
    591:   PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
    592:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

    594:   PetscCall(PetscRandomCreate(comm, &rctx));
    595:   PetscCall(PetscRandomSetFromOptions(rctx));
    596:   PetscCall(MatCreateVecs(A, &x, &s1));
    597:   PetscCall(VecDuplicate(x, &y));
    598:   PetscCall(VecDuplicate(s1, &s2));

    600:   *flg = PETSC_TRUE;
    601:   for (k = 0; k < n; k++) {
    602:     PetscCall(VecSetRandom(x, rctx));
    603:     PetscCall(VecSetRandom(y, rctx));
    604:     if (rank == 0) PetscCall(PetscRandomGetValue(rctx, &a));
    605:     PetscCallMPI(MPI_Bcast(&a, 1, MPIU_SCALAR, 0, comm));

    607:     /* s2 = a*A*x + A*y */
    608:     PetscCall(MatMult(A, y, s2));  /* s2 = A*y */
    609:     PetscCall(MatMult(A, x, s1));  /* s1 = A*x */
    610:     PetscCall(VecAXPY(s2, a, s1)); /* s2 = a s1 + s2 */

    612:     /* s1 = A * (a x + y) */
    613:     PetscCall(VecAXPY(y, a, x)); /* y = a x + y */
    614:     PetscCall(MatMult(A, y, s1));
    615:     PetscCall(VecNorm(s1, NORM_INFINITY, &normA));

    617:     PetscCall(VecAXPY(s2, -1.0, s1)); /* s2 = - s1 + s2 */
    618:     PetscCall(VecNorm(s2, NORM_INFINITY, &norm));
    619:     if (norm / normA > 100. * PETSC_MACHINE_EPSILON) {
    620:       *flg = PETSC_FALSE;
    621:       PetscCall(PetscInfo(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)));
    622:       break;
    623:     }
    624:   }
    625:   PetscCall(PetscRandomDestroy(&rctx));
    626:   PetscCall(VecDestroy(&x));
    627:   PetscCall(VecDestroy(&y));
    628:   PetscCall(VecDestroy(&s1));
    629:   PetscCall(VecDestroy(&s2));
    630:   PetscFunctionReturn(PETSC_SUCCESS);
    631: }