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

  4: typedef struct {
  5:   Mat A;
  6: } Mat_HT;

  8: PetscErrorCode MatMult_HT(Mat N,Vec x,Vec y)
  9: {
 10:   Mat_HT         *Na = (Mat_HT*)N->data;

 14:   MatMultHermitianTranspose(Na->A,x,y);
 15:   return(0);
 16: }

 18: PetscErrorCode MatMultAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)
 19: {
 20:   Mat_HT         *Na = (Mat_HT*)N->data;

 24:   MatMultHermitianTransposeAdd(Na->A,v1,v2,v3);
 25:   return(0);
 26: }

 28: PetscErrorCode MatMultHermitianTranspose_HT(Mat N,Vec x,Vec y)
 29: {
 30:   Mat_HT         *Na = (Mat_HT*)N->data;

 34:   MatMult(Na->A,x,y);
 35:   return(0);
 36: }

 38: PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N,Vec v1,Vec v2,Vec v3)
 39: {
 40:   Mat_HT         *Na = (Mat_HT*)N->data;

 44:   MatMultAdd(Na->A,v1,v2,v3);
 45:   return(0);
 46: }

 48: PetscErrorCode MatDestroy_HT(Mat N)
 49: {
 50:   Mat_HT         *Na = (Mat_HT*)N->data;

 54:   MatDestroy(&Na->A);
 55:   PetscFree(N->data);
 56:   return(0);
 57: }

 59: PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat* m)
 60: {
 61:   Mat_HT         *Na = (Mat_HT*)N->data;

 65:   if (op == MAT_COPY_VALUES) {
 66:     MatHermitianTranspose(Na->A,MAT_INITIAL_MATRIX,m);
 67:   } else if (op == MAT_DO_NOT_COPY_VALUES) {
 68:     MatDuplicate(Na->A,MAT_DO_NOT_COPY_VALUES,m);
 69:     MatHermitianTranspose(*m,MAT_INPLACE_MATRIX,m);
 70:   } else SETERRQ(PetscObjectComm((PetscObject)N),PETSC_ERR_SUP,"MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
 71:   return(0);
 72: }

 74: /*@
 75:       MatCreateHermitianTranspose - Creates a new matrix object that behaves like A'*

 77:    Collective on Mat

 79:    Input Parameter:
 80: .   A  - the (possibly rectangular) matrix

 82:    Output Parameter:
 83: .   N - the matrix that represents A'*

 85:    Level: intermediate

 87:    Notes:
 88:     The hermitian transpose A' is NOT actually formed! Rather the new matrix
 89:           object performs the matrix-vector product by using the MatMultHermitianTranspose() on
 90:           the original matrix

 92: .seealso: MatCreateNormal(), MatMult(), MatMultHermitianTranspose(), MatCreate()

 94: @*/
 95: PetscErrorCode  MatCreateHermitianTranspose(Mat A,Mat *N)
 96: {
 98:   PetscInt       m,n;
 99:   Mat_HT         *Na;

102:   MatGetLocalSize(A,&m,&n);
103:   MatCreate(PetscObjectComm((PetscObject)A),N);
104:   MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
105:   PetscLayoutSetUp((*N)->rmap);
106:   PetscLayoutSetUp((*N)->cmap);
107:   PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);

109:   PetscNewLog(*N,&Na);
110:   (*N)->data = (void*) Na;
111:   PetscObjectReference((PetscObject)A);
112:   Na->A      = A;

114:   (*N)->ops->destroy                    = MatDestroy_HT;
115:   (*N)->ops->mult                       = MatMult_HT;
116:   (*N)->ops->multadd                    = MatMultAdd_HT;
117:   (*N)->ops->multhermitiantranspose     = MatMultHermitianTranspose_HT;
118:   (*N)->ops->multhermitiantransposeadd  = MatMultHermitianTransposeAdd_HT;
119:   (*N)->ops->duplicate                  = MatDuplicate_HT;
120:   (*N)->assembled                       = PETSC_TRUE;

122:   MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
123:   MatSetUp(*N);
124:   return(0);
125: }