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,am,an,rstart,rend;
 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;
 22: #if defined(PETSC_HAVE_ELEMENTAL)
 23:   PetscBool      Test_MatMatMult=PETSC_TRUE;
 24: #else
 25:   PetscBool      Test_MatMatMult=PETSC_FALSE;
 26: #endif
 27:   PetscMPIInt    size;

 29:   PetscInitialize(&argc,&argv,(char*)0,help);
 30:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

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

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

 64:   /* Test MatMatMult() */
 65:   if (Test_MatMatMult) {
 66: #if !defined(PETSC_HAVE_ELEMENTAL)
 67:     if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This test requires ELEMENTAL");
 68: #endif
 69:     MatTranspose(A,MAT_INITIAL_MATRIX,&B); /* B = A^T */
 70:     MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C); /* C = B*A = A^T*A */
 71:     MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C);

 73:     /* Test B*A*x = C*x for n random vector x */
 74:     MatMatMultEqual(B,A,C,10,&equal);
 75:     if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
 76:     MatDestroy(&C);

 78:     MatMatMultSymbolic(B,A,fill,&C);
 79:     for (i=0; i<2; i++) {
 80:       /* Repeat the numeric product to test reuse of the previous symbolic product */
 81:       MatMatMultNumeric(B,A,C);
 82: 
 83:       MatMatMultEqual(B,A,C,10,&equal);
 84:       if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x");
 85:     }
 86:     MatDestroy(&C);
 87:     MatDestroy(&B);
 88:   }

 90:   /* Test MatTransposeMatMult() */
 91:   PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&iselemental);
 92:   if (!iselemental) {
 93:     MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&D); /* D = A^T*A */
 94:     MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&D);
 95:     /* MatView(D,PETSC_VIEWER_STDOUT_WORLD); */
 96:     MatTransposeMatMultEqual(A,A,D,10,&equal);
 97:     if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x");
 98:     MatDestroy(&D);

100:     /* Test D*x = A^T*C*A*x, where C is in AIJ format */
101:     MatGetLocalSize(A,&am,&an);
102:     MatCreate(PETSC_COMM_WORLD,&C);
103:     if (size == 1) {
104:       MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,am,am);
105:     } else {
106:       MatSetSizes(C,am,am,PETSC_DECIDE,PETSC_DECIDE);
107:     }
108:     MatSetFromOptions(C);
109:     MatSetUp(C);
110:     MatGetOwnershipRange(C,&rstart,&rend);
111:     v[0] = 1.0;
112:     for (i=rstart; i<rend; i++) {
113:       MatSetValues(C,1,&i,1,&i,v,INSERT_VALUES);
114:     }
115:     MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
116:     MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);

118:     /* B = C*A, D = A^T*B */
119:     MatMatMult(C,A,MAT_INITIAL_MATRIX,1.0,&B);
120:     MatTransposeMatMult(A,B,MAT_INITIAL_MATRIX,fill,&D);
121:     MatTransposeMatMultEqual(A,B,D,10,&equal);
122:     if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*B*x");

124:     MatDestroy(&D);
125:     MatDestroy(&C);
126:     MatDestroy(&B);
127:   }

129:   MatDestroy(&A);
130:   PetscFree(v);
131:   PetscFinalize();
132:   return(0);
133: }