Actual source code: transm.c


  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;

 12:   PetscFunctionBegin;
 13:   PetscCall(MatMultTranspose(Na->A, x, y));
 14:   PetscFunctionReturn(PETSC_SUCCESS);
 15: }

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

 21:   PetscFunctionBegin;
 22:   PetscCall(MatMultTransposeAdd(Na->A, v1, v2, v3));
 23:   PetscFunctionReturn(PETSC_SUCCESS);
 24: }

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

 30:   PetscFunctionBegin;
 31:   PetscCall(MatMult(Na->A, x, y));
 32:   PetscFunctionReturn(PETSC_SUCCESS);
 33: }

 35: PetscErrorCode MatMultTransposeAdd_Transpose(Mat N, Vec v1, Vec v2, Vec v3)
 36: {
 37:   Mat_Transpose *Na = (Mat_Transpose *)N->data;

 39:   PetscFunctionBegin;
 40:   PetscCall(MatMultAdd(Na->A, v1, v2, v3));
 41:   PetscFunctionReturn(PETSC_SUCCESS);
 42: }

 44: PetscErrorCode MatDestroy_Transpose(Mat N)
 45: {
 46:   Mat_Transpose *Na = (Mat_Transpose *)N->data;

 48:   PetscFunctionBegin;
 49:   PetscCall(MatDestroy(&Na->A));
 50:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
 51:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
 52:   PetscCall(PetscFree(N->data));
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: PetscErrorCode MatDuplicate_Transpose(Mat N, MatDuplicateOption op, Mat *m)
 57: {
 58:   Mat_Transpose *Na = (Mat_Transpose *)N->data;

 60:   PetscFunctionBegin;
 61:   if (op == MAT_COPY_VALUES) {
 62:     PetscCall(MatTranspose(Na->A, MAT_INITIAL_MATRIX, m));
 63:   } else if (op == MAT_DO_NOT_COPY_VALUES) {
 64:     PetscCall(MatDuplicate(Na->A, MAT_DO_NOT_COPY_VALUES, m));
 65:     PetscCall(MatTranspose(*m, MAT_INPLACE_MATRIX, m));
 66:   } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

 70: PetscErrorCode MatCreateVecs_Transpose(Mat A, Vec *r, Vec *l)
 71: {
 72:   Mat_Transpose *Aa = (Mat_Transpose *)A->data;

 74:   PetscFunctionBegin;
 75:   PetscCall(MatCreateVecs(Aa->A, l, r));
 76:   PetscFunctionReturn(PETSC_SUCCESS);
 77: }

 79: PetscErrorCode MatAXPY_Transpose(Mat Y, PetscScalar a, Mat X, MatStructure str)
 80: {
 81:   Mat_Transpose *Ya = (Mat_Transpose *)Y->data;
 82:   Mat_Transpose *Xa = (Mat_Transpose *)X->data;
 83:   Mat            M  = Ya->A;
 84:   Mat            N  = Xa->A;

 86:   PetscFunctionBegin;
 87:   PetscCall(MatAXPY(M, a, N, str));
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: PetscErrorCode MatHasOperation_Transpose(Mat mat, MatOperation op, PetscBool *has)
 92: {
 93:   Mat_Transpose *X = (Mat_Transpose *)mat->data;
 94:   PetscFunctionBegin;

 96:   *has = PETSC_FALSE;
 97:   if (op == MATOP_MULT) {
 98:     PetscCall(MatHasOperation(X->A, MATOP_MULT_TRANSPOSE, has));
 99:   } else if (op == MATOP_MULT_TRANSPOSE) {
100:     PetscCall(MatHasOperation(X->A, MATOP_MULT, has));
101:   } else if (op == MATOP_MULT_ADD) {
102:     PetscCall(MatHasOperation(X->A, MATOP_MULT_TRANSPOSE_ADD, has));
103:   } else if (op == MATOP_MULT_TRANSPOSE_ADD) {
104:     PetscCall(MatHasOperation(X->A, MATOP_MULT_ADD, has));
105:   } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }

109: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat D)
110: {
111:   Mat            A, B, C, Ain, Bin, Cin;
112:   PetscBool      Aistrans, Bistrans, Cistrans;
113:   PetscInt       Atrans, Btrans, Ctrans;
114:   MatProductType ptype;

116:   PetscFunctionBegin;
117:   MatCheckProduct(D, 1);
118:   A = D->product->A;
119:   B = D->product->B;
120:   C = D->product->C;
121:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATTRANSPOSEVIRTUAL, &Aistrans));
122:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATTRANSPOSEVIRTUAL, &Bistrans));
123:   PetscCall(PetscObjectTypeCompare((PetscObject)C, MATTRANSPOSEVIRTUAL, &Cistrans));
124:   PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
125:   Atrans = 0;
126:   Ain    = A;
127:   while (Aistrans) {
128:     Atrans++;
129:     PetscCall(MatTransposeGetMat(Ain, &Ain));
130:     PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATTRANSPOSEVIRTUAL, &Aistrans));
131:   }
132:   Btrans = 0;
133:   Bin    = B;
134:   while (Bistrans) {
135:     Btrans++;
136:     PetscCall(MatTransposeGetMat(Bin, &Bin));
137:     PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATTRANSPOSEVIRTUAL, &Bistrans));
138:   }
139:   Ctrans = 0;
140:   Cin    = C;
141:   while (Cistrans) {
142:     Ctrans++;
143:     PetscCall(MatTransposeGetMat(Cin, &Cin));
144:     PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATTRANSPOSEVIRTUAL, &Cistrans));
145:   }
146:   Atrans = Atrans % 2;
147:   Btrans = Btrans % 2;
148:   Ctrans = Ctrans % 2;
149:   ptype  = D->product->type; /* same product type by default */
150:   if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
151:   if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
152:   if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;

154:   if (Atrans || Btrans || Ctrans) {
155:     ptype = MATPRODUCT_UNSPECIFIED;
156:     switch (D->product->type) {
157:     case MATPRODUCT_AB:
158:       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
159:         /* TODO custom implementation ? */
160:       } else if (Atrans) { /* At * B */
161:         ptype = MATPRODUCT_AtB;
162:       } else { /* A * Bt */
163:         ptype = MATPRODUCT_ABt;
164:       }
165:       break;
166:     case MATPRODUCT_AtB:
167:       if (Atrans && Btrans) { /* A * Bt */
168:         ptype = MATPRODUCT_ABt;
169:       } else if (Atrans) { /* A * B */
170:         ptype = MATPRODUCT_AB;
171:       } else { /* At * Bt we do not have support for this */
172:         /* TODO custom implementation ? */
173:       }
174:       break;
175:     case MATPRODUCT_ABt:
176:       if (Atrans && Btrans) { /* At * B */
177:         ptype = MATPRODUCT_AtB;
178:       } else if (Atrans) { /* At * Bt we do not have support for this */
179:         /* TODO custom implementation ? */
180:       } else { /* A * B */
181:         ptype = MATPRODUCT_AB;
182:       }
183:       break;
184:     case MATPRODUCT_PtAP:
185:       if (Atrans) { /* PtAtP */
186:         /* TODO custom implementation ? */
187:       } else { /* RARt */
188:         ptype = MATPRODUCT_RARt;
189:       }
190:       break;
191:     case MATPRODUCT_RARt:
192:       if (Atrans) { /* RAtRt */
193:         /* TODO custom implementation ? */
194:       } else { /* PtAP */
195:         ptype = MATPRODUCT_PtAP;
196:       }
197:       break;
198:     case MATPRODUCT_ABC:
199:       /* TODO custom implementation ? */
200:       break;
201:     default:
202:       SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
203:     }
204:   }
205:   PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
206:   PetscCall(MatProductSetType(D, ptype));
207:   PetscCall(MatProductSetFromOptions(D));
208:   PetscFunctionReturn(PETSC_SUCCESS);
209: }

211: PetscErrorCode MatGetDiagonal_Transpose(Mat A, Vec v)
212: {
213:   Mat_Transpose *Aa = (Mat_Transpose *)A->data;

215:   PetscFunctionBegin;
216:   PetscCall(MatGetDiagonal(Aa->A, v));
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: PetscErrorCode MatConvert_Transpose(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
221: {
222:   Mat_Transpose *Aa = (Mat_Transpose *)A->data;
223:   PetscBool      flg;

225:   PetscFunctionBegin;
226:   PetscCall(MatHasOperation(Aa->A, MATOP_TRANSPOSE, &flg));
227:   if (flg) {
228:     Mat B;

230:     PetscCall(MatTranspose(Aa->A, MAT_INITIAL_MATRIX, &B));
231:     if (reuse != MAT_INPLACE_MATRIX) {
232:       PetscCall(MatConvert(B, newtype, reuse, newmat));
233:       PetscCall(MatDestroy(&B));
234:     } else {
235:       PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
236:       PetscCall(MatHeaderReplace(A, &B));
237:     }
238:   } else { /* use basic converter as fallback */
239:     PetscCall(MatConvert_Basic(A, newtype, reuse, newmat));
240:   }
241:   PetscFunctionReturn(PETSC_SUCCESS);
242: }

244: PetscErrorCode MatTransposeGetMat_Transpose(Mat A, Mat *M)
245: {
246:   Mat_Transpose *Aa = (Mat_Transpose *)A->data;

248:   PetscFunctionBegin;
249:   *M = Aa->A;
250:   PetscFunctionReturn(PETSC_SUCCESS);
251: }

253: /*@
254:       MatTransposeGetMat - Gets the `Mat` object stored inside a `MATTRANSPOSEVIRTUAL`

256:    Logically Collective

258:    Input Parameter:
259: .   A  - the `MATTRANSPOSEVIRTUAL` matrix

261:    Output Parameter:
262: .   M - the matrix object stored inside `A`

264:    Level: intermediate

266: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`
267: @*/
268: PetscErrorCode MatTransposeGetMat(Mat A, Mat *M)
269: {
270:   PetscFunctionBegin;
274:   PetscUseMethod(A, "MatTransposeGetMat_C", (Mat, Mat *), (A, M));
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

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

281:   Level: advanced

283: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`,
284:           `MATNORMALHERMITIAN`, `MATNORMAL`
285: M*/

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

290:    Collective

292:    Input Parameter:
293: .   A  - the (possibly rectangular) matrix

295:    Output Parameter:
296: .   N - the matrix that represents A'

298:    Level: intermediate

300:    Note:
301:     The transpose A' is NOT actually formed! Rather the new matrix
302:           object performs the matrix-vector product by using the `MatMultTranspose()` on
303:           the original matrix

305: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `MatCreateNormal()`, `MatMult()`, `MatMultTranspose()`, `MatCreate()`,
306:           `MATNORMALHERMITIAN`
307: @*/
308: PetscErrorCode MatCreateTranspose(Mat A, Mat *N)
309: {
310:   PetscInt       m, n;
311:   Mat_Transpose *Na;
312:   VecType        vtype;

314:   PetscFunctionBegin;
315:   PetscCall(MatGetLocalSize(A, &m, &n));
316:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
317:   PetscCall(MatSetSizes(*N, n, m, PETSC_DECIDE, PETSC_DECIDE));
318:   PetscCall(PetscLayoutSetUp((*N)->rmap));
319:   PetscCall(PetscLayoutSetUp((*N)->cmap));
320:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATTRANSPOSEVIRTUAL));

322:   PetscCall(PetscNew(&Na));
323:   (*N)->data = (void *)Na;
324:   PetscCall(PetscObjectReference((PetscObject)A));
325:   Na->A = A;

327:   (*N)->ops->destroy               = MatDestroy_Transpose;
328:   (*N)->ops->mult                  = MatMult_Transpose;
329:   (*N)->ops->multadd               = MatMultAdd_Transpose;
330:   (*N)->ops->multtranspose         = MatMultTranspose_Transpose;
331:   (*N)->ops->multtransposeadd      = MatMultTransposeAdd_Transpose;
332:   (*N)->ops->duplicate             = MatDuplicate_Transpose;
333:   (*N)->ops->getvecs               = MatCreateVecs_Transpose;
334:   (*N)->ops->axpy                  = MatAXPY_Transpose;
335:   (*N)->ops->hasoperation          = MatHasOperation_Transpose;
336:   (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
337:   (*N)->ops->getdiagonal           = MatGetDiagonal_Transpose;
338:   (*N)->ops->convert               = MatConvert_Transpose;
339:   (*N)->assembled                  = PETSC_TRUE;

341:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatTransposeGetMat_Transpose));
342:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_Transpose));
343:   PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
344:   PetscCall(MatGetVecType(A, &vtype));
345:   PetscCall(MatSetVecType(*N, vtype));
346: #if defined(PETSC_HAVE_DEVICE)
347:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
348: #endif
349:   PetscCall(MatSetUp(*N));
350:   PetscFunctionReturn(PETSC_SUCCESS);
351: }