Actual source code: htransm.c

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

  3: static PetscErrorCode MatProductSetFromOptions_HT(Mat D)
  4: {
  5:   Mat            A, B, C, Ain, Bin, Cin;
  6:   PetscBool      Aistrans, Bistrans, Cistrans;
  7:   PetscInt       Atrans, Btrans, Ctrans;
  8:   MatProductType ptype;

 10:   PetscFunctionBegin;
 11:   MatCheckProduct(D, 1);
 12:   A = D->product->A;
 13:   B = D->product->B;
 14:   C = D->product->C;
 15:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
 16:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
 17:   PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
 18:   PetscCheck(Aistrans || Bistrans || Cistrans, PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "This should not happen");
 19:   Atrans = 0;
 20:   Ain    = A;
 21:   while (Aistrans) {
 22:     Atrans++;
 23:     PetscCall(MatShellGetScalingShifts(Ain, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
 24:     PetscCall(MatHermitianTransposeGetMat(Ain, &Ain));
 25:     PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
 26:   }
 27:   Btrans = 0;
 28:   Bin    = B;
 29:   while (Bistrans) {
 30:     Btrans++;
 31:     PetscCall(MatShellGetScalingShifts(Bin, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
 32:     PetscCall(MatHermitianTransposeGetMat(Bin, &Bin));
 33:     PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
 34:   }
 35:   Ctrans = 0;
 36:   Cin    = C;
 37:   while (Cistrans) {
 38:     Ctrans++;
 39:     PetscCall(MatShellGetScalingShifts(Cin, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
 40:     PetscCall(MatHermitianTransposeGetMat(Cin, &Cin));
 41:     PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
 42:   }
 43:   Atrans = Atrans % 2;
 44:   Btrans = Btrans % 2;
 45:   Ctrans = Ctrans % 2;
 46:   ptype  = D->product->type; /* same product type by default */
 47:   if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
 48:   if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
 49:   if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;

 51:   if (Atrans || Btrans || Ctrans) {
 52:     PetscCheck(!PetscDefined(USE_COMPLEX), PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for complex Hermitian transpose matrices");
 53:     ptype = MATPRODUCT_UNSPECIFIED;
 54:     switch (D->product->type) {
 55:     case MATPRODUCT_AB:
 56:       if (Atrans && Btrans) { /* At * Bt we do not have support for this */
 57:         /* TODO custom implementation ? */
 58:       } else if (Atrans) { /* At * B */
 59:         ptype = MATPRODUCT_AtB;
 60:       } else { /* A * Bt */
 61:         ptype = MATPRODUCT_ABt;
 62:       }
 63:       break;
 64:     case MATPRODUCT_AtB:
 65:       if (Atrans && Btrans) { /* A * Bt */
 66:         ptype = MATPRODUCT_ABt;
 67:       } else if (Atrans) { /* A * B */
 68:         ptype = MATPRODUCT_AB;
 69:       } else { /* At * Bt we do not have support for this */
 70:         /* TODO custom implementation ? */
 71:       }
 72:       break;
 73:     case MATPRODUCT_ABt:
 74:       if (Atrans && Btrans) { /* At * B */
 75:         ptype = MATPRODUCT_AtB;
 76:       } else if (Atrans) { /* At * Bt we do not have support for this */
 77:         /* TODO custom implementation ? */
 78:       } else { /* A * B */
 79:         ptype = MATPRODUCT_AB;
 80:       }
 81:       break;
 82:     case MATPRODUCT_PtAP:
 83:       if (Atrans) { /* PtAtP */
 84:         /* TODO custom implementation ? */
 85:       } else { /* RARt */
 86:         ptype = MATPRODUCT_RARt;
 87:       }
 88:       break;
 89:     case MATPRODUCT_RARt:
 90:       if (Atrans) { /* RAtRt */
 91:         /* TODO custom implementation ? */
 92:       } else { /* PtAP */
 93:         ptype = MATPRODUCT_PtAP;
 94:       }
 95:       break;
 96:     case MATPRODUCT_ABC:
 97:       /* TODO custom implementation ? */
 98:       break;
 99:     default:
100:       SETERRQ(PetscObjectComm((PetscObject)D), PETSC_ERR_SUP, "ProductType %s is not supported", MatProductTypes[D->product->type]);
101:     }
102:   }
103:   PetscCall(MatProductReplaceMats(Ain, Bin, Cin, D));
104:   PetscCall(MatProductSetType(D, ptype));
105:   PetscCall(MatProductSetFromOptions(D));
106:   PetscFunctionReturn(PETSC_SUCCESS);
107: }
108: static PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
109: {
110:   Mat A;

112:   PetscFunctionBegin;
113:   PetscCall(MatShellGetContext(N, &A));
114:   PetscCall(MatMultHermitianTranspose(A, x, y));
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

118: static PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
119: {
120:   Mat A;

122:   PetscFunctionBegin;
123:   PetscCall(MatShellGetContext(N, &A));
124:   PetscCall(MatMult(A, x, y));
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

128: static PetscErrorCode MatDestroy_HT(Mat N)
129: {
130:   Mat A;

132:   PetscFunctionBegin;
133:   PetscCall(MatShellGetContext(N, &A));
134:   PetscCall(MatDestroy(&A));
135:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL));
136: #if !defined(PETSC_USE_COMPLEX)
137:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
138: #endif
139:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
140:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
141:   PetscFunctionReturn(PETSC_SUCCESS);
142: }

144: static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
145: {
146:   Mat A, C;

148:   PetscFunctionBegin;
149:   PetscCall(MatShellGetContext(N, &A));
150:   PetscCall(MatDuplicate(A, op, &C));
151:   PetscCall(MatCreateHermitianTranspose(C, m));
152:   PetscCall(MatDestroy(&C));
153:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
154:   PetscFunctionReturn(PETSC_SUCCESS);
155: }

157: static PetscErrorCode MatHasOperation_HT(Mat mat, MatOperation op, PetscBool *has)
158: {
159:   Mat A;

161:   PetscFunctionBegin;
162:   PetscCall(MatShellGetContext(mat, &A));
163:   *has = PETSC_FALSE;
164:   if (op == MATOP_MULT || op == MATOP_MULT_ADD) {
165:     PetscCall(MatHasOperation(A, MATOP_MULT_HERMITIAN_TRANSPOSE, has));
166:     if (!*has) PetscCall(MatHasOperation(A, MATOP_MULT_TRANSPOSE, has));
167:   } else if (op == MATOP_MULT_HERMITIAN_TRANSPOSE || op == MATOP_MULT_HERMITIAN_TRANS_ADD || op == MATOP_MULT_TRANSPOSE || op == MATOP_MULT_TRANSPOSE_ADD) {
168:     PetscCall(MatHasOperation(A, MATOP_MULT, has));
169:   } else if (((void **)mat->ops)[op]) *has = PETSC_TRUE;
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }

173: static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
174: {
175:   PetscFunctionBegin;
176:   PetscCall(MatShellGetContext(N, M));
177:   PetscFunctionReturn(PETSC_SUCCESS);
178: }

180: /*@
181:   MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`

183:   Logically Collective

185:   Input Parameter:
186: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix

188:   Output Parameter:
189: . M - the matrix object stored inside A

191:   Level: intermediate

193: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
194: @*/
195: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
196: {
197:   PetscFunctionBegin;
200:   PetscAssertPointer(M, 2);
201:   PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
202:   PetscFunctionReturn(PETSC_SUCCESS);
203: }

205: static PetscErrorCode MatGetDiagonal_HT(Mat N, Vec v)
206: {
207:   Mat A;

209:   PetscFunctionBegin;
210:   PetscCall(MatShellGetContext(N, &A));
211:   PetscCall(MatGetDiagonal(A, v));
212:   PetscCall(VecConjugate(v));
213:   PetscFunctionReturn(PETSC_SUCCESS);
214: }

216: static PetscErrorCode MatCopy_HT(Mat A, Mat B, MatStructure str)
217: {
218:   Mat a, b;

220:   PetscFunctionBegin;
221:   PetscCall(MatShellGetContext(A, &a));
222:   PetscCall(MatShellGetContext(B, &b));
223:   PetscCall(MatCopy(a, b, str));
224:   PetscFunctionReturn(PETSC_SUCCESS);
225: }

227: static PetscErrorCode MatConvert_HT(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
228: {
229:   Mat         A;
230:   PetscScalar vscale = 1.0, vshift = 0.0;
231:   PetscBool   flg;

233:   PetscFunctionBegin;
234:   PetscCall(MatShellGetContext(N, &A));
235:   PetscCall(MatHasOperation(A, MATOP_HERMITIAN_TRANSPOSE, &flg));
236:   if (flg || N->ops->getrow) { /* if this condition is false, MatConvert_Shell() will be called in MatConvert_Basic(), so the following checks are not needed */
237:     PetscCall(MatShellGetScalingShifts(N, &vshift, &vscale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
238:   }
239:   if (flg) {
240:     Mat B;

242:     PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &B));
243:     if (reuse != MAT_INPLACE_MATRIX) {
244:       PetscCall(MatConvert(B, newtype, reuse, newmat));
245:       PetscCall(MatDestroy(&B));
246:     } else {
247:       PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
248:       PetscCall(MatHeaderReplace(N, &B));
249:     }
250:   } else { /* use basic converter as fallback */
251:     flg = (PetscBool)(N->ops->getrow != NULL);
252:     PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
253:   }
254:   if (flg) {
255:     PetscCall(MatScale(*newmat, vscale));
256:     PetscCall(MatShift(*newmat, vshift));
257:   }
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: /*MC
262:    MATHERMITIANTRANSPOSEVIRTUAL - "hermitiantranspose" - A matrix type that represents a virtual transpose of a matrix

264:   Level: advanced

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

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

271: .seealso: [](ch_matrices), `Mat`, `MATTRANSPOSEVIRTUAL`, `Mat`, `MatCreateHermitianTranspose()`, `MatCreateTranspose()`
272: M*/

274: /*@
275:   MatCreateHermitianTranspose - Creates a new matrix object of `MatType` `MATHERMITIANTRANSPOSEVIRTUAL` that behaves like A'*

277:   Collective

279:   Input Parameter:
280: . A - the (possibly rectangular) matrix

282:   Output Parameter:
283: . N - the matrix that represents A'*

285:   Level: intermediate

287:   Note:
288:   The Hermitian transpose A' is NOT actually formed! Rather the new matrix
289:   object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on
290:   the original matrix

292: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
293:           `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN`
294: @*/
295: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
296: {
297:   VecType vtype;

299:   PetscFunctionBegin;
300:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
301:   PetscCall(PetscLayoutReference(A->rmap, &((*N)->cmap)));
302:   PetscCall(PetscLayoutReference(A->cmap, &((*N)->rmap)));
303:   PetscCall(MatSetType(*N, MATSHELL));
304:   PetscCall(MatShellSetContext(*N, A));
305:   PetscCall(PetscObjectReference((PetscObject)A));

307:   PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
308:   PetscCall(MatGetVecType(A, &vtype));
309:   PetscCall(MatSetVecType(*N, vtype));
310: #if defined(PETSC_HAVE_DEVICE)
311:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
312: #endif
313:   PetscCall(MatSetUp(*N));

315:   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_HT));
316:   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_HT));
317:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
318: #if !defined(PETSC_USE_COMPLEX)
319:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
320: #endif
321:   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_HT));
322:   PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_HT));
323:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_HT));
324:   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_HT));
325:   PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_HT));

327:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
328: #if !defined(PETSC_USE_COMPLEX)
329:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
330: #endif
331:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HT));
332:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
333:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
334:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
335:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }