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:   PetscObjectComposeFunction((PetscObject)N,"MatHermitianTransposeGetMat_C",NULL);
 56: #if !defined(PETSC_USE_COMPLEX)
 57:   PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
 58: #endif
 59:   PetscFree(N->data);
 60:   return(0);
 61: }

 63: PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat* m)
 64: {
 65:   Mat_HT         *Na = (Mat_HT*)N->data;

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

 78: PetscErrorCode MatCreateVecs_HT(Mat N,Vec *r, Vec *l)
 79: {
 80:   Mat_HT         *Na = (Mat_HT*)N->data;

 84:   MatCreateVecs(Na->A,l,r);
 85:   return(0);
 86: }

 88: PetscErrorCode MatAXPY_HT(Mat Y,PetscScalar a,Mat X,MatStructure str)
 89: {
 90:   Mat_HT         *Ya = (Mat_HT*)Y->data;
 91:   Mat_HT         *Xa = (Mat_HT*)X->data;
 92:   Mat              M = Ya->A;
 93:   Mat              N = Xa->A;

 97:   MatAXPY(M,a,N,str);
 98:   return(0);
 99: }

101: PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N,Mat *M)
102: {
103:   Mat_HT         *Na = (Mat_HT*)N->data;

106:   *M = Na->A;
107:   return(0);
108: }

110: /*@
111:       MatHermitianTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT

113:    Logically collective on Mat

115:    Input Parameter:
116: .   A  - the MATTRANSPOSE matrix

118:    Output Parameter:
119: .   M - the matrix object stored inside A

121:    Level: intermediate

123: .seealso: MatCreateHermitianTranspose()

125: @*/
126: PetscErrorCode MatHermitianTransposeGetMat(Mat A,Mat *M)
127: {

134:   PetscUseMethod(A,"MatHermitianTransposeGetMat_C",(Mat,Mat*),(A,M));
135:   return(0);
136: }

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

141:    Collective on Mat

143:    Input Parameter:
144: .   A  - the (possibly rectangular) matrix

146:    Output Parameter:
147: .   N - the matrix that represents A'*

149:    Level: intermediate

151:    Notes:
152:     The hermitian transpose A' is NOT actually formed! Rather the new matrix
153:           object performs the matrix-vector product by using the MatMultHermitianTranspose() on
154:           the original matrix

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

158: @*/
159: PetscErrorCode  MatCreateHermitianTranspose(Mat A,Mat *N)
160: {
162:   PetscInt       m,n;
163:   Mat_HT         *Na;

166:   MatGetLocalSize(A,&m,&n);
167:   MatCreate(PetscObjectComm((PetscObject)A),N);
168:   MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
169:   PetscLayoutSetUp((*N)->rmap);
170:   PetscLayoutSetUp((*N)->cmap);
171:   PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);

173:   PetscNewLog(*N,&Na);
174:   (*N)->data = (void*) Na;
175:   PetscObjectReference((PetscObject)A);
176:   Na->A      = A;

178:   (*N)->ops->destroy                    = MatDestroy_HT;
179:   (*N)->ops->mult                       = MatMult_HT;
180:   (*N)->ops->multadd                    = MatMultAdd_HT;
181:   (*N)->ops->multhermitiantranspose     = MatMultHermitianTranspose_HT;
182:   (*N)->ops->multhermitiantransposeadd  = MatMultHermitianTransposeAdd_HT;
183:   (*N)->ops->duplicate                  = MatDuplicate_HT;
184:   (*N)->ops->getvecs                    = MatCreateVecs_HT;
185:   (*N)->ops->axpy                       = MatAXPY_HT;
186:   (*N)->assembled                       = PETSC_TRUE;

188:   PetscObjectComposeFunction((PetscObject)(*N),"MatHermitianTransposeGetMat_C",MatHermitianTransposeGetMat_HT);
189: #if !defined(PETSC_USE_COMPLEX)
190:   PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatHermitianTransposeGetMat_HT);
191: #endif
192:   MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
193:   MatSetUp(*N);
194:   return(0);
195: }