Actual source code: ex104.c
petsc-3.12.4 2020-02-04
Report Typos and Errors

      1: static char help[] = "Test MatMatMult(), MatTranspose(), MatTransposeMatMult() for Dense and Elemental matrices.\n\n";
      2: /*
      3:  Example:
      4:    mpiexec -n <np> ./ex104 -mat_type elemental
      5: */

      7:  #include <petscmat.h>

      9: int main(int argc,char **argv)
     10: {
     11:   Mat            A,B,C,D;
     12:   PetscInt       i,M=10,N=5,j,nrows,ncols,am,an,rstart,rend;
     14:   PetscRandom    r;
     15:   PetscBool      equal,iselemental;
     16:   PetscReal      fill = 1.0;
     17:   IS             isrows,iscols;
     18:   const PetscInt *rows,*cols;
     19:   PetscScalar    *v,rval;
     20: #if defined(PETSC_HAVE_ELEMENTAL)
     21:   PetscBool      Test_MatMatMult=PETSC_TRUE;
     22: #else
     23:   PetscBool      Test_MatMatMult=PETSC_FALSE;
     24: #endif
     25:   PetscMPIInt    size;

     27:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     28:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

     30:   PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL);
     31:   PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL);
     32:   MatCreate(PETSC_COMM_WORLD,&A);
     33:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
     34:   MatSetType(A,MATDENSE);
     35:   MatSetFromOptions(A);
     36:   MatSetUp(A);
     37:   PetscRandomCreate(PETSC_COMM_WORLD,&r);
     38:   PetscRandomSetFromOptions(r);

     40:   /* Set local matrix entries */
     41:   MatGetOwnershipIS(A,&isrows,&iscols);
     42:   ISGetLocalSize(isrows,&nrows);
     43:   ISGetIndices(isrows,&rows);
     44:   ISGetLocalSize(iscols,&ncols);
     45:   ISGetIndices(iscols,&cols);
     46:   PetscMalloc1(nrows*ncols,&v);
     47:   for (i=0; i<nrows; i++) {
     48:     for (j=0; j<ncols; j++) {
     49:       PetscRandomGetValue(r,&rval);
     50:       v[i*ncols+j] = rval;
     51:     }
     52:   }
     53:   MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES);
     54:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
     55:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
     56:   ISRestoreIndices(isrows,&rows);
     57:   ISRestoreIndices(iscols,&cols);
     58:   ISDestroy(&isrows);
     59:   ISDestroy(&iscols);
     60:   PetscRandomDestroy(&r);

     62:   /* Test MatTranspose() */
     63:   MatCreateTranspose(A,&C);
     64:   MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
     65:   MatMultEqual(C,B,10,&equal);
     66:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
     67:   MatTranspose(A,MAT_REUSE_MATRIX,&B); /* B = A^T */
     68:   MatMultEqual(C,B,10,&equal);
     69:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
     70:   MatDestroy(&B);
     71:   MatDuplicate(A,MAT_COPY_VALUES,&B);
     72:   MatTranspose(B,MAT_INPLACE_MATRIX,&B);
     73:   MatMultEqual(C,B,10,&equal);
     74:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T");
     75:   MatDestroy(&B);
     76:   MatDestroy(&C);

     78:   /* Test MatMatMult() */
     79:   if (Test_MatMatMult) {
     80: #if !defined(PETSC_HAVE_ELEMENTAL)
     81:     if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This test requires ELEMENTAL");
     82: #endif
     83:     MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
     84:     MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */
     85:     MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C);

     87:     /* Test MatDuplicate for matrix product */
     88:     MatDuplicate(C,MAT_COPY_VALUES,&D);
     89:     MatDestroy(&D);

     91:     /* Test B*A*x = C*x for n random vector x */
     92:     MatMatMultEqual(B,A,C,10,&equal);
     93:     if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
     94:     MatDestroy(&C);

     96:     MatMatMultSymbolic(B,A,fill,&C);
     97:     for (i=0; i<2; i++) {
     98:       /* Repeat the numeric product to test reuse of the previous symbolic product */
     99:       MatMatMultNumeric(B,A,C);

    101:       MatMatMultEqual(B,A,C,10,&equal);
    102:       if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
    103:     }
    104:     MatDestroy(&C);
    105:     MatDestroy(&B);
    106:   }

    108:   /* Test MatTransposeMatMult() */
    109:   PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&iselemental);
    110:   if (!iselemental) {
    111:     MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A^T*A */
    112:     MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&D);

    114:     /* Test MatDuplicate for matrix product */
    115:     MatDuplicate(D,MAT_COPY_VALUES,&C);
    116:     MatDestroy(&C);

    118:     MatTransposeMatMultEqual(A,A,D,10,&equal);
    119:     if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x");
    120:     MatDestroy(&D);

    122:     /* Test D*x = A^T*C*A*x, where C is in AIJ format */
    123:     MatGetLocalSize(A,&am,&an);
    124:     MatCreate(PETSC_COMM_WORLD,&C);
    125:     if (size == 1) {
    126:       MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,am,am);
    127:     } else {
    128:       MatSetSizes(C,am,am,PETSC_DECIDE,PETSC_DECIDE);
    129:     }
    130:     MatSetFromOptions(C);
    131:     MatSetUp(C);
    132:     MatGetOwnershipRange(C,&rstart,&rend);
    133:     v[0] = 1.0;
    134:     for (i=rstart; i<rend; i++) {
    135:       MatSetValues(C,1,&i,1,&i,v,INSERT_VALUES);
    136:     }
    137:     MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
    138:     MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

    140:     /* B = C*A, D = A^T*B */
    141:     MatMatMult(C,A,MAT_INITIAL_MATRIX,1.0,&B);
    142:     MatTransposeMatMult(A,B,MAT_INITIAL_MATRIX,fill,&D);
    143:     MatTransposeMatMultEqual(A,B,D,10,&equal);
    144:     if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*B*x");

    146:     MatDestroy(&D);
    147:     MatDestroy(&C);
    148:     MatDestroy(&B);
    149:   }

    151:   /* Test MatMatTransposeMult() */
    152:   if (!iselemental) {
    153:     PetscReal diff, scale;
    154:     PetscInt  am, an, aM, aN;

    156:     MatGetLocalSize(A, &am, &an);
    157:     MatGetSize(A, &aM, &aN);
    158:     MatCreateDense(PetscObjectComm((PetscObject)A),PETSC_DECIDE, an, aM + 10, aN, NULL, &B);
    159:     MatSetRandom(B, NULL);
    160:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
    161:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
    162:     MatMatTransposeMult(A,B,MAT_INITIAL_MATRIX,fill,&D); /* D = A*A^T */

    164:     /* Test MatDuplicate for matrix product */
    165:     MatDuplicate(D,MAT_COPY_VALUES,&C);

    167:     MatMatTransposeMult(A,B,MAT_REUSE_MATRIX,fill,&D);
    168:     MatAXPY(C, -1., D, SAME_NONZERO_PATTERN);

    170:     MatNorm(C, NORM_FROBENIUS, &diff);
    171:     MatNorm(D, NORM_FROBENIUS, &scale);
    172:     if (diff > PETSC_SMALL * scale) SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatMatTransposeMult() differs between MAT_INITIAL_MATRIX and MAT_REUSE_MATRIX");
    173:     MatDestroy(&C);

    175:     MatMatTransposeMultEqual(A,B,D,10,&equal);
    176:     if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x");
    177:     MatDestroy(&D);
    178:     MatDestroy(&B);

    180:   }

    182:   MatDestroy(&A);
    183:   PetscFree(v);
    184:   PetscFinalize();
    185:   return ierr;
    186: }

    188: /*TEST

    190:     test:
    191:       output_file: output/ex104.out

    193:     test:
    194:       suffix: 2
    195:       nsize: 2
    196:       output_file: output/ex104.out

    198:     test:
    199:       suffix: 3
    200:       nsize: 4
    201:       output_file: output/ex104.out
    202:       args: -M 23 -N 31

    204:     test:
    205:       suffix: 4
    206:       nsize: 4
    207:       output_file: output/ex104.out
    208:       args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via cyclic

    210:     test:
    211:       suffix: 5
    212:       nsize: 4
    213:       output_file: output/ex104.out
    214:       args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via allgatherv

    216: TEST*/