Actual source code: htransm.c

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

  3: typedef struct {
  4:   Mat A;
  5: } Mat_HT;

  7: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_HermitianTranspose(Mat D)
  8: {
  9:   Mat            A, B, C, Ain, Bin, Cin;
 10:   PetscBool      Aistrans, Bistrans, Cistrans;
 11:   PetscInt       Atrans, Btrans, Ctrans;
 12:   MatProductType ptype;

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

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

113:   PetscFunctionBegin;
114:   PetscCall(MatMultHermitianTranspose(Na->A, x, y));
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

118: static PetscErrorCode MatMultAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
119: {
120:   Mat_HT *Na = (Mat_HT *)N->data;

122:   PetscFunctionBegin;
123:   PetscCall(MatMultHermitianTransposeAdd(Na->A, v1, v2, v3));
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

127: static PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
128: {
129:   Mat_HT *Na = (Mat_HT *)N->data;

131:   PetscFunctionBegin;
132:   PetscCall(MatMult(Na->A, x, y));
133:   PetscFunctionReturn(PETSC_SUCCESS);
134: }

136: static PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
137: {
138:   Mat_HT *Na = (Mat_HT *)N->data;

140:   PetscFunctionBegin;
141:   PetscCall(MatMultAdd(Na->A, v1, v2, v3));
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }

145: static PetscErrorCode MatDestroy_HT(Mat N)
146: {
147:   Mat_HT *Na = (Mat_HT *)N->data;

149:   PetscFunctionBegin;
150:   PetscCall(MatDestroy(&Na->A));
151:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL));
152:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL));
153: #if !defined(PETSC_USE_COMPLEX)
154:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL));
155: #endif
156:   PetscCall(PetscFree(N->data));
157:   PetscFunctionReturn(PETSC_SUCCESS);
158: }

160: static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
161: {
162:   Mat_HT *Na = (Mat_HT *)N->data;

164:   PetscFunctionBegin;
165:   if (op == MAT_COPY_VALUES) {
166:     PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, m));
167:   } else if (op == MAT_DO_NOT_COPY_VALUES) {
168:     PetscCall(MatDuplicate(Na->A, MAT_DO_NOT_COPY_VALUES, m));
169:     PetscCall(MatHermitianTranspose(*m, MAT_INPLACE_MATRIX, m));
170:   } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
171:   PetscFunctionReturn(PETSC_SUCCESS);
172: }

174: static PetscErrorCode MatCreateVecs_HT(Mat N, Vec *r, Vec *l)
175: {
176:   Mat_HT *Na = (Mat_HT *)N->data;

178:   PetscFunctionBegin;
179:   PetscCall(MatCreateVecs(Na->A, l, r));
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: static PetscErrorCode MatAXPY_HT(Mat Y, PetscScalar a, Mat X, MatStructure str)
184: {
185:   Mat_HT *Ya = (Mat_HT *)Y->data;
186:   Mat_HT *Xa = (Mat_HT *)X->data;
187:   Mat     M  = Ya->A;
188:   Mat     N  = Xa->A;

190:   PetscFunctionBegin;
191:   PetscCall(MatAXPY(M, a, N, str));
192:   PetscFunctionReturn(PETSC_SUCCESS);
193: }

195: static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
196: {
197:   Mat_HT *Na = (Mat_HT *)N->data;

199:   PetscFunctionBegin;
200:   *M = Na->A;
201:   PetscFunctionReturn(PETSC_SUCCESS);
202: }

204: /*@
205:   MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`

207:   Logically Collective

209:   Input Parameter:
210: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix

212:   Output Parameter:
213: . M - the matrix object stored inside A

215:   Level: intermediate

217: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
218: @*/
219: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
220: {
221:   PetscFunctionBegin;
224:   PetscAssertPointer(M, 2);
225:   PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
226:   PetscFunctionReturn(PETSC_SUCCESS);
227: }

229: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat);

231: static PetscErrorCode MatGetDiagonal_HT(Mat A, Vec v)
232: {
233:   Mat_HT *Na = (Mat_HT *)A->data;

235:   PetscFunctionBegin;
236:   PetscCall(MatGetDiagonal(Na->A, v));
237:   PetscCall(VecConjugate(v));
238:   PetscFunctionReturn(PETSC_SUCCESS);
239: }

241: static PetscErrorCode MatConvert_HT(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
242: {
243:   Mat_HT   *Na = (Mat_HT *)A->data;
244:   PetscBool flg;

246:   PetscFunctionBegin;
247:   PetscCall(MatHasOperation(Na->A, MATOP_HERMITIAN_TRANSPOSE, &flg));
248:   if (flg) {
249:     Mat B;

251:     PetscCall(MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, &B));
252:     if (reuse != MAT_INPLACE_MATRIX) {
253:       PetscCall(MatConvert(B, newtype, reuse, newmat));
254:       PetscCall(MatDestroy(&B));
255:     } else {
256:       PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
257:       PetscCall(MatHeaderReplace(A, &B));
258:     }
259:   } else { /* use basic converter as fallback */
260:     PetscCall(MatConvert_Basic(A, newtype, reuse, newmat));
261:   }
262:   PetscFunctionReturn(PETSC_SUCCESS);
263: }

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

268:   Level: advanced

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

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

276:   Collective

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

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

284:   Level: intermediate

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

291: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
292:           `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN`
293: @*/
294: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
295: {
296:   Mat_HT *Na;
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(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));

305:   PetscCall(PetscNew(&Na));
306:   (*N)->data = (void *)Na;
307:   PetscCall(PetscObjectReference((PetscObject)A));
308:   Na->A = A;

310:   (*N)->ops->destroy                   = MatDestroy_HT;
311:   (*N)->ops->mult                      = MatMult_HT;
312:   (*N)->ops->multadd                   = MatMultAdd_HT;
313:   (*N)->ops->multhermitiantranspose    = MatMultHermitianTranspose_HT;
314:   (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT;
315: #if !defined(PETSC_USE_COMPLEX)
316:   (*N)->ops->multtranspose    = MatMultHermitianTranspose_HT;
317:   (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_HT;
318: #endif
319:   (*N)->ops->duplicate = MatDuplicate_HT;
320:   (*N)->ops->getvecs   = MatCreateVecs_HT;
321:   (*N)->ops->axpy      = MatAXPY_HT;
322: #if !defined(PETSC_USE_COMPLEX)
323:   (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
324: #endif
325:   (*N)->ops->getdiagonal = MatGetDiagonal_HT;
326:   (*N)->ops->convert     = MatConvert_HT;
327:   (*N)->assembled        = PETSC_TRUE;

329:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
330:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HermitianTranspose));
331: #if !defined(PETSC_USE_COMPLEX)
332:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
333: #endif
334:   PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
335:   PetscCall(MatGetVecType(A, &vtype));
336:   PetscCall(MatSetVecType(*N, vtype));
337: #if defined(PETSC_HAVE_DEVICE)
338:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
339: #endif
340:   PetscCall(MatSetUp(*N));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }