Actual source code: transm.c

  1: #include <../src/mat/impls/shell/shell.h>

  3: static PetscErrorCode MatMult_Transpose(Mat N, Vec x, Vec y)
  4: {
  5:   Mat A;

  7:   PetscFunctionBegin;
  8:   PetscCall(MatShellGetContext(N, &A));
  9:   PetscCall(MatMultTranspose(A, x, y));
 10:   PetscFunctionReturn(PETSC_SUCCESS);
 11: }

 13: static PetscErrorCode MatMultTranspose_Transpose(Mat N, Vec x, Vec y)
 14: {
 15:   Mat A;

 17:   PetscFunctionBegin;
 18:   PetscCall(MatShellGetContext(N, &A));
 19:   PetscCall(MatMult(A, x, y));
 20:   PetscFunctionReturn(PETSC_SUCCESS);
 21: }

 23: static PetscErrorCode MatDestroy_Transpose(Mat N)
 24: {
 25:   Mat A;

 27:   PetscFunctionBegin;
 28:   PetscCall(MatShellGetContext(N, &A));
 29:   PetscCall(MatDestroy(&A));
 30:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
 31:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
 32:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
 33:   PetscFunctionReturn(PETSC_SUCCESS);
 34: }

 36: static PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m)
 37: {
 38:   Mat A, C;

 40:   PetscFunctionBegin;
 41:   PetscCall(MatShellGetContext(N, &A));
 42:   PetscCall(MatDuplicate(A, op, &C));
 43:   PetscCall(MatCreateTranspose(C, m));
 44:   PetscCall(MatDestroy(&C));
 45:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 47: }

 49: static PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has)
 50: {
 51:   Mat A;

 53:   PetscFunctionBegin;
 54:   PetscCall(MatShellGetContext(mat, &A));
 55:   *has = PETSC_FALSE;
 56:   if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
 57:     PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
 58:   } else if (op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
 59:     PetscCall(MatHasOperation(A, MATOP_MULT, has));
 60:   } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
 61:   PetscFunctionReturn(PETSC_SUCCESS);
 62: }

 64: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
 65: {
 66:   Mat            A, B, C, Ain, Bin, Cin;
 67:   PetscBool      Aistrans, Bistrans, Cistrans;
 68:   PetscInt       Atrans, Btrans, Ctrans;
 69:   MatProductType ptype;

 71:   PetscFunctionBegin;
 72:   MatCheckProduct(D, 1);
 73:   A = D->product->A;
 74:   B = D->product->B;
 75:   C = D->product->C;
 76:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
 77:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
 78:   PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
 79:   PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
 80:   Atrans = 0;
 81:   Ain    = A;
 82:   while (Aistrans) {
 83:     Atrans++;
 84:     PetscCall(MatTransposeGetMat(Ain, &Ain));
 85:     PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
 86:   }
 87:   Btrans = 0;
 88:   Bin    = B;
 89:   while (Bistrans) {
 90:     Btrans++;
 91:     PetscCall(MatTransposeGetMat(Bin, &Bin));
 92:     PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
 93:   }
 94:   Ctrans = 0;
 95:   Cin    = C;
 96:   while (Cistrans) {
 97:     Ctrans++;
 98:     PetscCall(MatTransposeGetMat(Cin, &Cin));
 99:     PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
100:   }
101:   Atrans = Atrans % 2;
102:   Btrans = Btrans % 2;
103:   Ctrans = Ctrans % 2;
104:   ptype  = D->product->type; /* same product type by default */
105:   if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
106:   if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
107:   if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;

109:   if (Atrans || Btrans || Ctrans) {
110:     ptype = MATPRODUCT_UNSPECIFIED;
111:     switch (D->product->type) {
112:     case MATPRODUCT_AB:
113:       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
114:         /* TODO custom implementation ? */
115:       } else if (Atrans) { /* At * B */
116:         ptype = MATPRODUCT_AtB;
117:       } else { /* A * Bt */
118:         ptype = MATPRODUCT_ABt;
119:       }
120:       break;
121:     case MATPRODUCT_AtB:
122:       if (Atrans && Btrans) { /* A * Bt */
123:         ptype = MATPRODUCT_ABt;
124:       } else if (Atrans) { /* A * B */
125:         ptype = MATPRODUCT_AB;
126:       } else { /* At * Bt we do not have support for this */
127:         /* TODO custom implementation ? */
128:       }
129:       break;
130:     case MATPRODUCT_ABt:
131:       if (Atrans && Btrans) { /* At * B */
132:         ptype = MATPRODUCT_AtB;
133:       } else if (Atrans) { /* At * Bt we do not have support for this */
134:         /* TODO custom implementation ? */
135:       } else { /* A * B */
136:         ptype = MATPRODUCT_AB;
137:       }
138:       break;
139:     case MATPRODUCT_PtAP:
140:       if (Atrans) { /* PtAtP */
141:         /* TODO custom implementation ? */
142:       } else { /* RARt */
143:         ptype = MATPRODUCT_RARt;
144:       }
145:       break;
146:     case MATPRODUCT_RARt:
147:       if (Atrans) { /* RAtRt */
148:         /* TODO custom implementation ? */
149:       } else { /* PtAP */
150:         ptype = MATPRODUCT_PtAP;
151:       }
152:       break;
153:     case MATPRODUCT_ABC:
154:       /* TODO custom implementation ? */
155:       break;
156:     default:
157:       SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
158:     }
159:   }
160:   PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
161:   PetscCall(MatProductSetType(D, ptype));
162:   PetscCall(MatProductSetFromOptions(D));
163:   PetscFunctionReturn(PETSC_SUCCESS);
164: }

166: static PetscErrorCode MatGetDiagonal_Transpose(Mat N, Vec v)
167: {
168:   Mat A;

170:   PetscFunctionBegin;
171:   PetscCall(MatShellGetContext(N, &A));
172:   PetscCall(MatGetDiagonal(A, v));
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

176: static PetscErrorCode MatCopy_Transpose(Mat A, Mat B, MatStructure str)
177: {
178:   Mat a, b;

180:   PetscFunctionBegin;
181:   PetscCall(MatShellGetContext(A, &a));
182:   PetscCall(MatShellGetContext(B, &b));
183:   PetscCall(MatCopy(a, b, str));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: static PetscErrorCode MatConvert_Transpose(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
188: {
189:   Mat       A;
190:   PetscBool flg;

192:   PetscFunctionBegin;
193:   PetscCall(MatShellGetContext(N, &A));
194:   PetscCall(MatHasOperation(A, MATOP_TRANSPOSE, &flg));
195:   if (flg) {
196:     Mat B;

198:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &B));
199:     if (reuse != MAT_INPLACE_MATRIX) {
200:       PetscCall(MatConvert(B, newtype, reuse, newmat));
201:       PetscCall(MatDestroy(&B));
202:     } else {
203:       PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
204:       PetscCall(MatHeaderReplace(N, &B));
205:     }
206:   } else { /* use basic converter as fallback */
207:     PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
208:   }
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }

212: static PetscErrorCode MatTransposeGetMat_Transpose(Mat N, Mat *M)
213: {
214:   PetscFunctionBegin;
215:   PetscCall(MatShellGetContext(N, M));
216:   PetscFunctionReturn(PETSC_SUCCESS);
217: }

219: /*@
220:   MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`

222:   Logically Collective

224:   Input Parameter:
225: . A - the `MATTRANSPOSEVIRTUAL` matrix

227:   Output Parameter:
228: . M - the matrix object stored inside `A`

230:   Level: intermediate

232: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
233: @*/
234: PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
235: {
236:   PetscFunctionBegin;
239:   PetscAssertPointer(M, 2);
240:   PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
241:   PetscFunctionReturn(PETSC_SUCCESS);
242: }

244: /*MC
245:    MATTRANSPOSEVIRTUAL - "transpose" - A matrix type that represents a virtual transpose of a matrix

247:   Level: advanced

249:   Developer Notes:
250:   This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code

252:   Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage

254: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
255:           `MATNORMALHERMITIAN`, `MATNORMAL`
256: M*/

258: /*@
259:   MatCreateTranspose - Creates a new matrix `MATTRANSPOSEVIRTUAL` object that behaves like A'

261:   Collective

263:   Input Parameter:
264: . A - the (possibly rectangular) matrix

266:   Output Parameter:
267: . N - the matrix that represents A'

269:   Level: intermediate

271:   Note:
272:   The transpose A' is NOT actually formed! Rather the new matrix
273:   object performs the matrix-vector product by using the `MatMultTranspose()` on
274:   the original matrix

276: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
277:           `MATNORMALHERMITIAN`
278: @*/
279: PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
280: {
281:   VecType vtype;

283:   PetscFunctionBegin;
284:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
285:   PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
286:   PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
287:   PetscCall(MatSetType(*N, MATSHELL));
288:   PetscCall(MatShellSetContext(*N, A));
289:   PetscCall(PetscObjectReference((PetscObject)A));

291:   PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
292:   PetscCall(MatGetVecType(A, &vtype));
293:   PetscCall(MatSetVecType(*N, vtype));
294: #if defined(PETSC_HAVE_DEVICE)
295:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
296: #endif
297:   PetscCall(MatSetUp(*N));

299:   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Transpose));
300:   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Transpose));
301:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultTranspose_Transpose));
302:   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Transpose));
303:   PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_Transpose));
304:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Transpose));
305:   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Transpose));
306:   PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_Transpose));

308:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
309:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
310:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
311:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
312:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
313:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL));
314:   PetscFunctionReturn(PETSC_SUCCESS);
315: }