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

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

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

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

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

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

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

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

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

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

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

 54:   MatDestroy(&Na->A);
 55:   PetscObjectComposeFunction((PetscObject)N,"MatTransposeGetMat_C",NULL);
 56:   PetscFree(N->data);
 57:   return(0);
 58: }

 60: PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat* m)
 61: {
 62:   Mat_Transpose  *Na = (Mat_Transpose*)N->data;

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

 75: PetscErrorCode MatCreateVecs_Transpose(Mat A,Vec *r, Vec *l)
 76: {
 77:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

 81:   MatCreateVecs(Aa->A,l,r);
 82:   return(0);
 83: }

 85: PetscErrorCode MatAXPY_Transpose(Mat Y,PetscScalar a,Mat X,MatStructure str)
 86: {
 87:   Mat_Transpose  *Ya = (Mat_Transpose*)Y->data;
 88:   Mat_Transpose  *Xa = (Mat_Transpose*)X->data;
 89:   Mat              M = Ya->A;
 90:   Mat              N = Xa->A;

 94:   MatAXPY(M,a,N,str);
 95:   return(0);
 96: }

 98: PetscErrorCode MatHasOperation_Transpose(Mat mat,MatOperation op,PetscBool *has)
 99: {
100:   Mat_Transpose  *X = (Mat_Transpose*)mat->data;

104:   *has = PETSC_FALSE;
105:   if (op == MATOP_MULT) {
106:     MatHasOperation(X->A,MATOP_MULT_TRANSPOSE,has);
107:   } else if (op == MATOP_MULT_TRANSPOSE) {
108:     MatHasOperation(X->A,MATOP_MULT,has);
109:   } else if (op == MATOP_MULT_ADD) {
110:     MatHasOperation(X->A,MATOP_MULT_TRANSPOSE_ADD,has);
111:   } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
112:     MatHasOperation(X->A,MATOP_MULT_ADD,has);
113:   } else if (((void**)mat->ops)[op]) *has = PETSC_TRUE;
114:   return(0);
115: }

117: PetscErrorCode MatTransposeGetMat_Transpose(Mat A,Mat *M)
118: {
119:   Mat_Transpose  *Aa = (Mat_Transpose*)A->data;

122:   *M = Aa->A;
123:   return(0);
124: }

126: /*@
127:       MatTransposeGetMat - Gets the Mat object stored inside a MATTRANSPOSEMAT

129:    Logically collective on Mat

131:    Input Parameter:
132: .   A  - the MATTRANSPOSE matrix

134:    Output Parameter:
135: .   M - the matrix object stored inside A

137:    Level: intermediate

139: .seealso: MatCreateTranspose()

141: @*/
142: PetscErrorCode MatTransposeGetMat(Mat A,Mat *M)
143: {

150:   PetscUseMethod(A,"MatTransposeGetMat_C",(Mat,Mat*),(A,M));
151:   return(0);
152: }

154: /*@
155:       MatCreateTranspose - Creates a new matrix object that behaves like A'

157:    Collective on Mat

159:    Input Parameter:
160: .   A  - the (possibly rectangular) matrix

162:    Output Parameter:
163: .   N - the matrix that represents A'

165:    Level: intermediate

167:    Notes:
168:     The transpose A' is NOT actually formed! Rather the new matrix
169:           object performs the matrix-vector product by using the MatMultTranspose() on
170:           the original matrix

172: .seealso: MatCreateNormal(), MatMult(), MatMultTranspose(), MatCreate()

174: @*/
175: PetscErrorCode  MatCreateTranspose(Mat A,Mat *N)
176: {
178:   PetscInt       m,n;
179:   Mat_Transpose  *Na;

182:   MatGetLocalSize(A,&m,&n);
183:   MatCreate(PetscObjectComm((PetscObject)A),N);
184:   MatSetSizes(*N,n,m,PETSC_DECIDE,PETSC_DECIDE);
185:   PetscLayoutSetUp((*N)->rmap);
186:   PetscLayoutSetUp((*N)->cmap);
187:   PetscObjectChangeTypeName((PetscObject)*N,MATTRANSPOSEMAT);

189:   PetscNewLog(*N,&Na);
190:   (*N)->data = (void*) Na;
191:   PetscObjectReference((PetscObject)A);
192:   Na->A      = A;

194:   (*N)->ops->destroy          = MatDestroy_Transpose;
195:   (*N)->ops->mult             = MatMult_Transpose;
196:   (*N)->ops->multadd          = MatMultAdd_Transpose;
197:   (*N)->ops->multtranspose    = MatMultTranspose_Transpose;
198:   (*N)->ops->multtransposeadd = MatMultTransposeAdd_Transpose;
199:   (*N)->ops->duplicate        = MatDuplicate_Transpose;
200:   (*N)->ops->getvecs          = MatCreateVecs_Transpose;
201:   (*N)->ops->axpy             = MatAXPY_Transpose;
202:   (*N)->ops->hasoperation     = MatHasOperation_Transpose;
203:   (*N)->assembled             = PETSC_TRUE;

205:   PetscObjectComposeFunction((PetscObject)(*N),"MatTransposeGetMat_C",MatTransposeGetMat_Transpose);
206:   MatSetBlockSizes(*N,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));
207:   MatSetUp(*N);
208:   return(0);
209: }