Actual source code: htransm.c

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

  3: PETSC_INTERN 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(MatHermitianTransposeGetMat(Ain, &Ain));
 24:     PetscCall(PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans));
 25:   }
 26:   Btrans = 0;
 27:   Bin    = B;
 28:   while (Bistrans) {
 29:     Btrans++;
 30:     PetscCall(MatHermitianTransposeGetMat(Bin, &Bin));
 31:     PetscCall(PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans));
 32:   }
 33:   Ctrans = 0;
 34:   Cin    = C;
 35:   while (Cistrans) {
 36:     Ctrans++;
 37:     PetscCall(MatHermitianTransposeGetMat(Cin, &Cin));
 38:     PetscCall(PetscObjectTypeCompare((PetscObject)Cin, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans));
 39:   }
 40:   Atrans = Atrans % 2;
 41:   Btrans = Btrans % 2;
 42:   Ctrans = Ctrans % 2;
 43:   ptype  = D->product->type; /* same product type by default */
 44:   if (Ain->symmetric == PETSC_BOOL3_TRUE) Atrans = 0;
 45:   if (Bin->symmetric == PETSC_BOOL3_TRUE) Btrans = 0;
 46:   if (Cin && Cin->symmetric == PETSC_BOOL3_TRUE) Ctrans = 0;

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

109:   PetscFunctionBegin;
110:   PetscCall(MatShellGetContext(N, &A));
111:   PetscCall(MatMultHermitianTranspose(A, x, y));
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: static PetscErrorCode MatMultHermitianTranspose_HT(Mat N, Vec x, Vec y)
116: {
117:   Mat A;

119:   PetscFunctionBegin;
120:   PetscCall(MatShellGetContext(N, &A));
121:   PetscCall(MatMult(A, x, y));
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }
124: static PetscErrorCode MatDestroy_HT(Mat N)
125: {
126:   Mat A;

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

140: static PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
141: {
142:   Mat A, C;

144:   PetscFunctionBegin;
145:   PetscCall(MatShellGetContext(N, &A));
146:   PetscCall(MatDuplicate(A, op, &C));
147:   PetscCall(MatCreateHermitianTranspose(C, m));
148:   PetscCall(MatDestroy(&C));
149:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(N, *m, SAME_NONZERO_PATTERN));
150:   PetscFunctionReturn(PETSC_SUCCESS);
151: }

153: static PetscErrorCode MatHasOperation_HT(Mat mat, MatOperation op, PetscBool *has)
154: {
155:   Mat A;

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

169: static PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
170: {
171:   PetscFunctionBegin;
172:   PetscCall(MatShellGetContext(N, M));
173:   PetscFunctionReturn(PETSC_SUCCESS);
174: }

176: /*@
177:   MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`

179:   Logically Collective

181:   Input Parameter:
182: . A - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix

184:   Output Parameter:
185: . M - the matrix object stored inside A

187:   Level: intermediate

189: .seealso: [](ch_matrices), `Mat`, `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
190: @*/
191: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
192: {
193:   PetscFunctionBegin;
196:   PetscAssertPointer(M, 2);
197:   PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
198:   PetscFunctionReturn(PETSC_SUCCESS);
199: }

201: static PetscErrorCode MatGetDiagonal_HT(Mat N, Vec v)
202: {
203:   Mat A;

205:   PetscFunctionBegin;
206:   PetscCall(MatShellGetContext(N, &A));
207:   PetscCall(MatGetDiagonal(A, v));
208:   PetscCall(VecConjugate(v));
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }

212: static PetscErrorCode MatCopy_HT(Mat A, Mat B, MatStructure str)
213: {
214:   Mat a, b;

216:   PetscFunctionBegin;
217:   PetscCall(MatShellGetContext(A, &a));
218:   PetscCall(MatShellGetContext(B, &b));
219:   PetscCall(MatCopy(a, b, str));
220:   PetscFunctionReturn(PETSC_SUCCESS);
221: }

223: static PetscErrorCode MatConvert_HT(Mat N, MatType newtype, MatReuse reuse, Mat *newmat)
224: {
225:   Mat       A;
226:   PetscBool flg;

228:   PetscFunctionBegin;
229:   PetscCall(MatShellGetContext(N, &A));
230:   PetscCall(MatHasOperation(A, MATOP_HERMITIAN_TRANSPOSE, &flg));
231:   if (flg) {
232:     Mat B;

234:     PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &B));
235:     if (reuse != MAT_INPLACE_MATRIX) {
236:       PetscCall(MatConvert(B, newtype, reuse, newmat));
237:       PetscCall(MatDestroy(&B));
238:     } else {
239:       PetscCall(MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B));
240:       PetscCall(MatHeaderReplace(N, &B));
241:     }
242:   } else { /* use basic converter as fallback */
243:     PetscCall(MatConvert_Basic(N, newtype, reuse, newmat));
244:   }
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

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

251:   Level: advanced

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

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

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

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

264:   Collective

266:   Input Parameter:
267: . A - the (possibly rectangular) matrix

269:   Output Parameter:
270: . N - the matrix that represents A'*

272:   Level: intermediate

274:   Note:
275:   The Hermitian transpose A' is NOT actually formed! Rather the new matrix
276:   object performs the matrix-vector product, `MatMult()`, by using the `MatMultHermitianTranspose()` on
277:   the original matrix

279: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
280:           `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`, `MATNORMAL`, `MATNORMALHERMITIAN`
281: @*/
282: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
283: {
284:   VecType vtype;

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

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

302:   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_HT));
303:   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_HT));
304:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
305: #if !defined(PETSC_USE_COMPLEX)
306:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMultHermitianTranspose_HT));
307: #endif
308:   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_HT));
309:   PetscCall(MatShellSetOperation(*N, MATOP_HAS_OPERATION, (void (*)(void))MatHasOperation_HT));
310:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_HT));
311:   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_HT));
312:   PetscCall(MatShellSetOperation(*N, MATOP_CONVERT, (void (*)(void))MatConvert_HT));

314:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
315: #if !defined(PETSC_USE_COMPLEX)
316:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT));
317: #endif
318:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HT));
319:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
320:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
321:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
322:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL));
323:   PetscFunctionReturn(PETSC_SUCCESS);
324: }