Actual source code: htransm.c


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

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

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

 15:   MatCheckProduct(D, 1);
 16:   A = D->product->A;
 17:   B = D->product->B;
 18:   C = D->product->C;
 19:   PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans);
 20:   PetscObjectTypeCompare((PetscObject)B, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans);
 21:   PetscObjectTypeCompare((PetscObject)C, MATHERMITIANTRANSPOSEVIRTUAL, &Cistrans);
 23:   Atrans = 0;
 24:   Ain    = A;
 25:   while (Aistrans) {
 26:     Atrans++;
 27:     MatHermitianTransposeGetMat(Ain, &Ain);
 28:     PetscObjectTypeCompare((PetscObject)Ain, MATHERMITIANTRANSPOSEVIRTUAL, &Aistrans);
 29:   }
 30:   Btrans = 0;
 31:   Bin    = B;
 32:   while (Bistrans) {
 33:     Btrans++;
 34:     MatHermitianTransposeGetMat(Bin, &Bin);
 35:     PetscObjectTypeCompare((PetscObject)Bin, MATHERMITIANTRANSPOSEVIRTUAL, &Bistrans);
 36:   }
 37:   Ctrans = 0;
 38:   Cin    = C;
 39:   while (Cistrans) {
 40:     Ctrans++;
 41:     MatHermitianTransposeGetMat(Cin, &Cin);
 42:     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) {
 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:   MatProductReplaceMats(Ain, Bin, Cin, D);
105:   MatProductSetType(D, ptype);
106:   MatProductSetFromOptions(D);
107:   return 0;
108: }
109: PetscErrorCode MatMult_HT(Mat N, Vec x, Vec y)
110: {
111:   Mat_HT *Na = (Mat_HT *)N->data;

113:   MatMultHermitianTranspose(Na->A, x, y);
114:   return 0;
115: }

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

121:   MatMultHermitianTransposeAdd(Na->A, v1, v2, v3);
122:   return 0;
123: }

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

129:   MatMult(Na->A, x, y);
130:   return 0;
131: }

133: PetscErrorCode MatMultHermitianTransposeAdd_HT(Mat N, Vec v1, Vec v2, Vec v3)
134: {
135:   Mat_HT *Na = (Mat_HT *)N->data;

137:   MatMultAdd(Na->A, v1, v2, v3);
138:   return 0;
139: }

141: PetscErrorCode MatDestroy_HT(Mat N)
142: {
143:   Mat_HT *Na = (Mat_HT *)N->data;

145:   MatDestroy(&Na->A);
146:   PetscObjectComposeFunction((PetscObject)N, "MatHermitianTransposeGetMat_C", NULL);
147:   PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_anytype_C", NULL);
148: #if !defined(PETSC_USE_COMPLEX)
149:   PetscObjectComposeFunction((PetscObject)N, "MatTransposeGetMat_C", NULL);
150: #endif
151:   PetscFree(N->data);
152:   return 0;
153: }

155: PetscErrorCode MatDuplicate_HT(Mat N, MatDuplicateOption op, Mat *m)
156: {
157:   Mat_HT *Na = (Mat_HT *)N->data;

159:   if (op == MAT_COPY_VALUES) {
160:     MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, m);
161:   } else if (op == MAT_DO_NOT_COPY_VALUES) {
162:     MatDuplicate(Na->A, MAT_DO_NOT_COPY_VALUES, m);
163:     MatHermitianTranspose(*m, MAT_INPLACE_MATRIX, m);
164:   } else SETERRQ(PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "MAT_SHARE_NONZERO_PATTERN not supported for this matrix type");
165:   return 0;
166: }

168: PetscErrorCode MatCreateVecs_HT(Mat N, Vec *r, Vec *l)
169: {
170:   Mat_HT *Na = (Mat_HT *)N->data;

172:   MatCreateVecs(Na->A, l, r);
173:   return 0;
174: }

176: PetscErrorCode MatAXPY_HT(Mat Y, PetscScalar a, Mat X, MatStructure str)
177: {
178:   Mat_HT *Ya = (Mat_HT *)Y->data;
179:   Mat_HT *Xa = (Mat_HT *)X->data;
180:   Mat     M  = Ya->A;
181:   Mat     N  = Xa->A;

183:   MatAXPY(M, a, N, str);
184:   return 0;
185: }

187: PetscErrorCode MatHermitianTransposeGetMat_HT(Mat N, Mat *M)
188: {
189:   Mat_HT *Na = (Mat_HT *)N->data;

191:   *M = Na->A;
192:   return 0;
193: }

195: /*@
196:       MatHermitianTransposeGetMat - Gets the `Mat` object stored inside a `MATHERMITIANTRANSPOSEVIRTUAL`

198:    Logically collective

200:    Input Parameter:
201: .   A  - the `MATHERMITIANTRANSPOSEVIRTUAL` matrix

203:    Output Parameter:
204: .   M - the matrix object stored inside A

206:    Level: intermediate

208: .seealso: `MATHERMITIANTRANSPOSEVIRTUAL`, `MatCreateHermitianTranspose()`
209: @*/
210: PetscErrorCode MatHermitianTransposeGetMat(Mat A, Mat *M)
211: {
215:   PetscUseMethod(A, "MatHermitianTransposeGetMat_C", (Mat, Mat *), (A, M));
216:   return 0;
217: }

219: PETSC_INTERN PetscErrorCode MatProductSetFromOptions_Transpose(Mat);

221: PetscErrorCode MatGetDiagonal_HT(Mat A, Vec v)
222: {
223:   Mat_HT *Na = (Mat_HT *)A->data;

225:   MatGetDiagonal(Na->A, v);
226:   VecConjugate(v);
227:   return 0;
228: }

230: PetscErrorCode MatConvert_HT(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
231: {
232:   Mat_HT   *Na = (Mat_HT *)A->data;
233:   PetscBool flg;

235:   MatHasOperation(Na->A, MATOP_HERMITIAN_TRANSPOSE, &flg);
236:   if (flg) {
237:     Mat B;

239:     MatHermitianTranspose(Na->A, MAT_INITIAL_MATRIX, &B);
240:     if (reuse != MAT_INPLACE_MATRIX) {
241:       MatConvert(B, newtype, reuse, newmat);
242:       MatDestroy(&B);
243:     } else {
244:       MatConvert(B, newtype, MAT_INPLACE_MATRIX, &B);
245:       MatHeaderReplace(A, &B);
246:     }
247:   } else { /* use basic converter as fallback */
248:     MatConvert_Basic(A, newtype, reuse, newmat);
249:   }
250:   return 0;
251: }

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

256:   Level: advanced

258: .seealso: `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: `MatCreateNormal()`, `MatMult()`, `MatMultHermitianTranspose()`, `MatCreate()`,
280:           `MATTRANSPOSEVIRTUAL`, `MatCreateTranspose()`, `MatHermitianTransposeGetMat()`
281: @*/
282: PetscErrorCode MatCreateHermitianTranspose(Mat A, Mat *N)
283: {
284:   PetscInt m, n;
285:   Mat_HT  *Na;
286:   VecType  vtype;

288:   MatGetLocalSize(A, &m, &n);
289:   MatCreate(PetscObjectComm((PetscObject)A), N);
290:   MatSetSizes(*N, n, m, PETSC_DECIDE, PETSC_DECIDE);
291:   PetscLayoutSetUp((*N)->rmap);
292:   PetscLayoutSetUp((*N)->cmap);
293:   PetscObjectChangeTypeName((PetscObject)*N, MATHERMITIANTRANSPOSEVIRTUAL);

295:   PetscNew(&Na);
296:   (*N)->data = (void *)Na;
297:   PetscObjectReference((PetscObject)A);
298:   Na->A = A;

300:   (*N)->ops->destroy                   = MatDestroy_HT;
301:   (*N)->ops->mult                      = MatMult_HT;
302:   (*N)->ops->multadd                   = MatMultAdd_HT;
303:   (*N)->ops->multhermitiantranspose    = MatMultHermitianTranspose_HT;
304:   (*N)->ops->multhermitiantransposeadd = MatMultHermitianTransposeAdd_HT;
305: #if !defined(PETSC_USE_COMPLEX)
306:   (*N)->ops->multtranspose    = MatMultHermitianTranspose_HT;
307:   (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_HT;
308: #endif
309:   (*N)->ops->duplicate = MatDuplicate_HT;
310:   (*N)->ops->getvecs   = MatCreateVecs_HT;
311:   (*N)->ops->axpy      = MatAXPY_HT;
312: #if !defined(PETSC_USE_COMPLEX)
313:   (*N)->ops->productsetfromoptions = MatProductSetFromOptions_Transpose;
314: #endif
315:   (*N)->ops->getdiagonal = MatGetDiagonal_HT;
316:   (*N)->ops->convert     = MatConvert_HT;
317:   (*N)->assembled        = PETSC_TRUE;

319:   PetscObjectComposeFunction((PetscObject)(*N), "MatHermitianTransposeGetMat_C", MatHermitianTransposeGetMat_HT);
320:   PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_anytype_C", MatProductSetFromOptions_HermitianTranspose);
321: #if !defined(PETSC_USE_COMPLEX)
322:   PetscObjectComposeFunction((PetscObject)(*N), "MatTransposeGetMat_C", MatHermitianTransposeGetMat_HT);
323: #endif
324:   MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs));
325:   MatGetVecType(A, &vtype);
326:   MatSetVecType(*N, vtype);
327: #if defined(PETSC_HAVE_DEVICE)
328:   MatBindToCPU(*N, A->boundtocpu);
329: #endif
330:   MatSetUp(*N);
331:   return 0;
332: }