Actual source code: multequal.c

      1: #define PETSCMAT_DLL

     3:  #include src/mat/matimpl.h

      7: /*@
      8:    MatMultEqual - Compares matrix-vector products of two matrices.

     10:    Collective on Mat

     12:    Input Parameters:
     13: +  A - the first matrix
     14: -  B - the second matrix
     15: -  n - number of random vectors to be tested

     17:    Output Parameter:
     18: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

     20:    Level: intermediate

     22:    Concepts: matrices^equality between
     23: @*/
     24: PetscErrorCode PETSCMAT_DLLEXPORT MatMultEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
     25: {
     27:   Vec            x,s1,s2;
     28:   PetscRandom    rctx;
     29:   PetscReal      r1,r2,tol=1.e-10;
     30:   PetscInt       am,an,bm,bn,k;
     31:   PetscScalar    none = -1.0;

     36:   MatGetLocalSize(A,&am,&an);
     37:   MatGetLocalSize(B,&bm,&bn);
     38:   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
     40:   PetscRandomCreate(A->comm,RANDOM_DEFAULT,&rctx);
     41:   VecCreate(A->comm,&x);
     42:   VecSetSizes(x,an,PETSC_DECIDE);
     43:   VecSetFromOptions(x);
     44:
     45:   VecCreate(A->comm,&s1);
     46:   VecSetSizes(s1,am,PETSC_DECIDE);
     47:   VecSetFromOptions(s1);
     48:   VecDuplicate(s1,&s2);
     49:
     50:   *flg = PETSC_TRUE;
     51:   for (k=0; k<n; k++) {
     52:     VecSetRandom(x,rctx);
     53:     MatMult(A,x,s1);
     54:     MatMult(B,x,s2);
     55:     VecNorm(s2,NORM_INFINITY,&r2);
     56:     if (r2 < tol){
     57:       VecNorm(s1,NORM_INFINITY,&r1);
     58:     } else {
     59:       VecAXPY(s2,none,s1);
     60:       VecNorm(s2,NORM_INFINITY,&r1);
     61:       r1 /= r2;
     62:     }
     63:     if (r1 > tol) {
     64:       *flg = PETSC_FALSE;
     65:       PetscLogInfo((0,"Error: %D-th MatMult() %g\n",k,r1));
     66:       break;
     67:     }
     68:   }
     69:   PetscRandomDestroy(rctx);
     70:   VecDestroy(x);
     71:   VecDestroy(s1);
     72:   VecDestroy(s2);
     73:   return(0);
     74: }

     78: /*@
     79:    MatMultAddEqual - Compares matrix-vector products of two matrices.

     81:    Collective on Mat

     83:    Input Parameters:
     84: +  A - the first matrix
     85: -  B - the second matrix
     86: -  n - number of random vectors to be tested

     88:    Output Parameter:
     89: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

     91:    Level: intermediate

     93:    Concepts: matrices^equality between
     94: @*/
     95: PetscErrorCode PETSCMAT_DLLEXPORT MatMultAddEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
     96: {
     98:   Vec            x,y,s1,s2;
     99:   PetscRandom    rctx;
    100:   PetscReal      r1,r2,tol=1.e-10;
    101:   PetscInt       am,an,bm,bn,k;
    102:   PetscScalar    none = -1.0;

    105:   MatGetLocalSize(A,&am,&an);
    106:   MatGetLocalSize(B,&bm,&bn);
    107:   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
    109:   PetscRandomCreate(A->comm,RANDOM_DEFAULT,&rctx);
    110:   VecCreate(A->comm,&x);
    111:   VecSetSizes(x,an,PETSC_DECIDE);
    112:   VecSetFromOptions(x);

    114:   VecCreate(A->comm,&s1);
    115:   VecSetSizes(s1,am,PETSC_DECIDE);
    116:   VecSetFromOptions(s1);
    117:   VecDuplicate(s1,&s2);
    118:   VecDuplicate(s1,&y);
    119:
    120:   *flg = PETSC_TRUE;
    121:   for (k=0; k<n; k++) {
    122:     VecSetRandom(x,rctx);
    123:     VecSetRandom(y,rctx);
    124:     MatMultAdd(A,x,y,s1);
    125:     MatMultAdd(B,x,y,s2);
    126:     VecNorm(s2,NORM_INFINITY,&r2);
    127:     if (r2 < tol){
    128:       VecNorm(s1,NORM_INFINITY,&r1);
    129:     } else {
    130:       VecAXPY(s2,none,s1);
    131:       VecNorm(s2,NORM_INFINITY,&r1);
    132:       r1 /= r2;
    133:     }
    134:     if (r1 > tol) {
    135:       *flg = PETSC_FALSE;
    136:       PetscLogInfo((0,"Error: %d-th MatMultAdd() %g\n",k,r1));
    137:       break;
    138:     }
    139:   }
    140:   PetscRandomDestroy(rctx);
    141:   VecDestroy(x);
    142:   VecDestroy(y);
    143:   VecDestroy(s1);
    144:   VecDestroy(s2);
    145:   return(0);
    146: }

    150: /*@
    151:    MatMultTransposeEqual - Compares matrix-vector products of two matrices.

    153:    Collective on Mat

    155:    Input Parameters:
    156: +  A - the first matrix
    157: -  B - the second matrix
    158: -  n - number of random vectors to be tested

    160:    Output Parameter:
    161: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

    163:    Level: intermediate

    165:    Concepts: matrices^equality between
    166: @*/
    167: PetscErrorCode PETSCMAT_DLLEXPORT MatMultTransposeEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
    168: {
    170:   Vec            x,s1,s2;
    171:   PetscRandom    rctx;
    172:   PetscReal      r1,r2,tol=1.e-10;
    173:   PetscInt       am,an,bm,bn,k;
    174:   PetscScalar    none = -1.0;

    177:   MatGetLocalSize(A,&am,&an);
    178:   MatGetLocalSize(B,&bm,&bn);
    179:   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
    181:   PetscRandomCreate(A->comm,RANDOM_DEFAULT,&rctx);
    182:   VecCreate(A->comm,&x);
    183:   VecSetSizes(x,am,PETSC_DECIDE);
    184:   VecSetFromOptions(x);
    185:
    186:   VecCreate(A->comm,&s1);
    187:   VecSetSizes(s1,an,PETSC_DECIDE);
    188:   VecSetFromOptions(s1);
    189:   VecDuplicate(s1,&s2);
    190:
    191:   *flg = PETSC_TRUE;
    192:   for (k=0; k<n; k++) {
    193:     VecSetRandom(x,rctx);
    194:     MatMultTranspose(A,x,s1);
    195:     MatMultTranspose(B,x,s2);
    196:     VecNorm(s2,NORM_INFINITY,&r2);
    197:     if (r2 < tol){
    198:       VecNorm(s1,NORM_INFINITY,&r1);
    199:     } else {
    200:       VecAXPY(s2,none,s1);
    201:       VecNorm(s2,NORM_INFINITY,&r1);
    202:       r1 /= r2;
    203:     }
    204:     if (r1 > tol) {
    205:       *flg = PETSC_FALSE;
    206:       PetscLogInfo((0,"Error: %d-th MatMultTranspose() %g\n",k,r1));
    207:       break;
    208:     }
    209:   }
    210:   PetscRandomDestroy(rctx);
    211:   VecDestroy(x);
    212:   VecDestroy(s1);
    213:   VecDestroy(s2);
    214:   return(0);
    215: }

    219: /*@
    220:    MatMultTransposeAddEqual - Compares matrix-vector products of two matrices.

    222:    Collective on Mat

    224:    Input Parameters:
    225: +  A - the first matrix
    226: -  B - the second matrix
    227: -  n - number of random vectors to be tested

    229:    Output Parameter:
    230: .  flg - PETSC_TRUE if the products are equal; PETSC_FALSE otherwise.

    232:    Level: intermediate

    234:    Concepts: matrices^equality between
    235: @*/
    236: PetscErrorCode PETSCMAT_DLLEXPORT MatMultTransposeAddEqual(Mat A,Mat B,PetscInt n,PetscTruth *flg)
    237: {
    239:   Vec            x,y,s1,s2;
    240:   PetscRandom    rctx;
    241:   PetscReal      r1,r2,tol=1.e-10;
    242:   PetscInt       am,an,bm,bn,k;
    243:   PetscScalar    none = -1.0;

    246:   MatGetLocalSize(A,&am,&an);
    247:   MatGetLocalSize(B,&bm,&bn);
    248:   if (am != bm || an != bn) SETERRQ4(PETSC_ERR_ARG_SIZ,"Mat A,Mat B: local dim %D %D %D %D",am,bm,an,bn);
    250:   PetscRandomCreate(A->comm,RANDOM_DEFAULT,&rctx);
    251:   VecCreate(A->comm,&x);
    252:   VecSetSizes(x,am,PETSC_DECIDE);
    253:   VecSetFromOptions(x);

    255:   VecCreate(A->comm,&s1);
    256:   VecSetSizes(s1,an,PETSC_DECIDE);
    257:   VecSetFromOptions(s1);
    258:   VecDuplicate(s1,&s2);
    259:   VecDuplicate(s1,&y);
    260:
    261:   *flg = PETSC_TRUE;
    262:   for (k=0; k<n; k++) {
    263:     VecSetRandom(x,rctx);
    264:     VecSetRandom(y,rctx);
    265:     MatMultTransposeAdd(A,x,y,s1);
    266:     MatMultTransposeAdd(B,x,y,s2);
    267:     VecNorm(s2,NORM_INFINITY,&r2);
    268:     if (r2 < tol){
    269:       VecNorm(s1,NORM_INFINITY,&r1);
    270:     } else {
    271:       VecAXPY(s2,none,s1);
    272:       VecNorm(s2,NORM_INFINITY,&r1);
    273:       r1 /= r2;
    274:     }
    275:     if (r1 > tol) {
    276:       *flg = PETSC_FALSE;
    277:       PetscLogInfo((0,"Error: %d-th MatMultTransposeAdd() %g\n",k,r1));
    278:       break;
    279:     }
    280:   }
    281:   PetscRandomDestroy(rctx);
    282:   VecDestroy(x);
    283:   VecDestroy(y);
    284:   VecDestroy(s1);
    285:   VecDestroy(s2);
    286:   return(0);
    287: }