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>

 11: int main(int argc,char **argv)
 12: {
 13:   Mat            A,B,C,D;
 14:   PetscInt       i,M=10,N=5,j,nrows,ncols;
 16:   PetscRandom    r;
 17:   PetscBool      equal,iselemental;
 18:   PetscReal      fill = 1.0;
 19:   IS             isrows,iscols;
 20:   const PetscInt *rows,*cols;
 21:   PetscScalar    *v,rval;

 23:   PetscInitialize(&argc,&argv,(char*)0,help);
 24:   PetscOptionsGetInt(NULL,"-M",&M,NULL);
 25:   PetscOptionsGetInt(NULL,"-N",&N,NULL);
 26:   MatCreate(PETSC_COMM_WORLD,&A);
 27:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N);
 28:   MatSetType(A,MATDENSE);
 29:   MatSetFromOptions(A);
 30:   MatSetUp(A);
 31:   PetscRandomCreate(PETSC_COMM_WORLD,&r);
 32:   PetscRandomSetFromOptions(r);

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

 57:   /* Test MatMatMult() */
 58:   MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
 59:   MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */

 61:   MatMatMultSymbolic(B,A,fill,&D); /* D = B*A = A^T*A */
 62:   for (i=0; i<2; i++) {
 63:     /* Repeat the numeric product to test reuse of the previous symbolic product */
 64:     MatMatMultNumeric(B,A,D);
 65:   }
 66:   MatMultEqual(C,D,10,&equal);
 67:   if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"C != D");
 68:   MatDestroy(&D);

 70:   /* Test MatTransposeMatMult() */
 71:   PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&iselemental);
 72:   if (!iselemental) {
 73:     MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A^T*A */
 74:     MatMultEqual(C,D,10,&equal);
 75:     if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"C != D");
 76:     MatDestroy(&D);
 77:   }

 79:   MatDestroy(&C);
 80:   MatDestroy(&B);
 81:   MatDestroy(&A);
 82:   PetscFinalize();
 83:   return(0);
 84: }