Actual source code: normmh.c

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

  3: typedef struct {
  4:   Mat         A;
  5:   Mat         D; /* local submatrix for diagonal part */
  6:   Vec         w, left, right, leftwork, rightwork;
  7:   PetscScalar scale;
  8: } Mat_NormalHermitian;

 10: static PetscErrorCode MatScale_NormalHermitian(Mat inA, PetscScalar scale)
 11: {
 12:   Mat_NormalHermitian *a = (Mat_NormalHermitian *)inA->data;

 14:   PetscFunctionBegin;
 15:   a->scale *= scale;
 16:   PetscFunctionReturn(PETSC_SUCCESS);
 17: }

 19: static PetscErrorCode MatDiagonalScale_NormalHermitian(Mat inA, Vec left, Vec right)
 20: {
 21:   Mat_NormalHermitian *a = (Mat_NormalHermitian *)inA->data;

 23:   PetscFunctionBegin;
 24:   if (left) {
 25:     if (!a->left) {
 26:       PetscCall(VecDuplicate(left, &a->left));
 27:       PetscCall(VecCopy(left, a->left));
 28:     } else {
 29:       PetscCall(VecPointwiseMult(a->left, left, a->left));
 30:     }
 31:   }
 32:   if (right) {
 33:     if (!a->right) {
 34:       PetscCall(VecDuplicate(right, &a->right));
 35:       PetscCall(VecCopy(right, a->right));
 36:     } else {
 37:       PetscCall(VecPointwiseMult(a->right, right, a->right));
 38:     }
 39:   }
 40:   PetscFunctionReturn(PETSC_SUCCESS);
 41: }

 43: static PetscErrorCode MatCreateSubMatrices_NormalHermitian(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
 44: {
 45:   Mat_NormalHermitian *a = (Mat_NormalHermitian *)mat->data;
 46:   Mat                  B = a->A, *suba;
 47:   IS                  *row;
 48:   PetscInt             M;

 50:   PetscFunctionBegin;
 51:   PetscCheck(!a->left && !a->right && irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented");
 52:   if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
 53:   PetscCall(MatGetSize(B, &M, NULL));
 54:   PetscCall(PetscMalloc1(n, &row));
 55:   PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]));
 56:   PetscCall(ISSetIdentity(row[0]));
 57:   for (M = 1; M < n; ++M) row[M] = row[0];
 58:   PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba));
 59:   for (M = 0; M < n; ++M) {
 60:     PetscCall(MatCreateNormalHermitian(suba[M], *submat + M));
 61:     ((Mat_NormalHermitian *)(*submat)[M]->data)->scale = a->scale;
 62:   }
 63:   PetscCall(ISDestroy(&row[0]));
 64:   PetscCall(PetscFree(row));
 65:   PetscCall(MatDestroySubMatrices(n, &suba));
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }

 69: static PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B)
 70: {
 71:   Mat_NormalHermitian *a = (Mat_NormalHermitian *)A->data;
 72:   Mat                  C, Aa = a->A;
 73:   IS                   row;

 75:   PetscFunctionBegin;
 76:   PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
 77:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row));
 78:   PetscCall(ISSetIdentity(row));
 79:   PetscCall(MatPermute(Aa, row, colp, &C));
 80:   PetscCall(ISDestroy(&row));
 81:   PetscCall(MatCreateNormalHermitian(C, B));
 82:   PetscCall(MatDestroy(&C));
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }

 86: static PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B)
 87: {
 88:   Mat_NormalHermitian *a = (Mat_NormalHermitian *)A->data;
 89:   Mat                  C;

 91:   PetscFunctionBegin;
 92:   PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented");
 93:   PetscCall(MatDuplicate(a->A, op, &C));
 94:   PetscCall(MatCreateNormalHermitian(C, B));
 95:   PetscCall(MatDestroy(&C));
 96:   if (op == MAT_COPY_VALUES) ((Mat_NormalHermitian *)(*B)->data)->scale = a->scale;
 97:   PetscFunctionReturn(PETSC_SUCCESS);
 98: }

100: static PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str)
101: {
102:   Mat_NormalHermitian *a = (Mat_NormalHermitian *)A->data, *b = (Mat_NormalHermitian *)B->data;

104:   PetscFunctionBegin;
105:   PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented");
106:   PetscCall(MatCopy(a->A, b->A, str));
107:   b->scale = a->scale;
108:   PetscCall(VecDestroy(&b->left));
109:   PetscCall(VecDestroy(&b->right));
110:   PetscCall(VecDestroy(&b->leftwork));
111:   PetscCall(VecDestroy(&b->rightwork));
112:   PetscFunctionReturn(PETSC_SUCCESS);
113: }

115: static PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y)
116: {
117:   Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
118:   Vec                  in;

120:   PetscFunctionBegin;
121:   in = x;
122:   if (Na->right) {
123:     if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork));
124:     PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in));
125:     in = Na->rightwork;
126:   }
127:   PetscCall(MatMult(Na->A, in, Na->w));
128:   PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y));
129:   if (Na->left) PetscCall(VecPointwiseMult(y, Na->left, y));
130:   PetscCall(VecScale(y, Na->scale));
131:   PetscFunctionReturn(PETSC_SUCCESS);
132: }

134: static PetscErrorCode MatMultHermitianAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
135: {
136:   Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
137:   Vec                  in;

139:   PetscFunctionBegin;
140:   in = v1;
141:   if (Na->right) {
142:     if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork));
143:     PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in));
144:     in = Na->rightwork;
145:   }
146:   PetscCall(MatMult(Na->A, in, Na->w));
147:   PetscCall(VecScale(Na->w, Na->scale));
148:   if (Na->left) {
149:     PetscCall(MatMultHermitianTranspose(Na->A, Na->w, v3));
150:     PetscCall(VecPointwiseMult(v3, Na->left, v3));
151:     PetscCall(VecAXPY(v3, 1.0, v2));
152:   } else {
153:     PetscCall(MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3));
154:   }
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }

158: static PetscErrorCode MatMultHermitianTranspose_Normal(Mat N, Vec x, Vec y)
159: {
160:   Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
161:   Vec                  in;

163:   PetscFunctionBegin;
164:   in = x;
165:   if (Na->left) {
166:     if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork));
167:     PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in));
168:     in = Na->leftwork;
169:   }
170:   PetscCall(MatMult(Na->A, in, Na->w));
171:   PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y));
172:   if (Na->right) PetscCall(VecPointwiseMult(y, Na->right, y));
173:   PetscCall(VecScale(y, Na->scale));
174:   PetscFunctionReturn(PETSC_SUCCESS);
175: }

177: static PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
178: {
179:   Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
180:   Vec                  in;

182:   PetscFunctionBegin;
183:   in = v1;
184:   if (Na->left) {
185:     if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork));
186:     PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in));
187:     in = Na->leftwork;
188:   }
189:   PetscCall(MatMult(Na->A, in, Na->w));
190:   PetscCall(VecScale(Na->w, Na->scale));
191:   if (Na->right) {
192:     PetscCall(MatMultHermitianTranspose(Na->A, Na->w, v3));
193:     PetscCall(VecPointwiseMult(v3, Na->right, v3));
194:     PetscCall(VecAXPY(v3, 1.0, v2));
195:   } else {
196:     PetscCall(MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3));
197:   }
198:   PetscFunctionReturn(PETSC_SUCCESS);
199: }

201: static PetscErrorCode MatDestroy_NormalHermitian(Mat N)
202: {
203:   Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;

205:   PetscFunctionBegin;
206:   PetscCall(MatDestroy(&Na->A));
207:   PetscCall(MatDestroy(&Na->D));
208:   PetscCall(VecDestroy(&Na->w));
209:   PetscCall(VecDestroy(&Na->left));
210:   PetscCall(VecDestroy(&Na->right));
211:   PetscCall(VecDestroy(&Na->leftwork));
212:   PetscCall(VecDestroy(&Na->rightwork));
213:   PetscCall(PetscFree(N->data));
214:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalHermitianGetMat_C", NULL));
215:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL));
216:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL));
217: #if defined(PETSC_HAVE_HYPRE)
218:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_hypre_C", NULL));
219: #endif
220:   PetscFunctionReturn(PETSC_SUCCESS);
221: }

223: /*
224:       Slow, nonscalable version
225: */
226: static PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
227: {
228:   Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
229:   Mat                  A  = Na->A;
230:   PetscInt             i, j, rstart, rend, nnz;
231:   const PetscInt      *cols;
232:   PetscScalar         *diag, *work, *values;
233:   const PetscScalar   *mvalues;

235:   PetscFunctionBegin;
236:   PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
237:   PetscCall(PetscArrayzero(work, A->cmap->N));
238:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
239:   for (i = rstart; i < rend; i++) {
240:     PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
241:     for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
242:     PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
243:   }
244:   PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
245:   rstart = N->cmap->rstart;
246:   rend   = N->cmap->rend;
247:   PetscCall(VecGetArray(v, &values));
248:   PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
249:   PetscCall(VecRestoreArray(v, &values));
250:   PetscCall(PetscFree2(diag, work));
251:   PetscCall(VecScale(v, Na->scale));
252:   PetscFunctionReturn(PETSC_SUCCESS);
253: }

255: static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
256: {
257:   Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
258:   Mat                  M, A = Na->A;

260:   PetscFunctionBegin;
261:   PetscCall(MatGetDiagonalBlock(A, &M));
262:   PetscCall(MatCreateNormalHermitian(M, &Na->D));
263:   *D = Na->D;
264:   PetscFunctionReturn(PETSC_SUCCESS);
265: }

267: static PetscErrorCode MatNormalHermitianGetMat_NormalHermitian(Mat A, Mat *M)
268: {
269:   Mat_NormalHermitian *Aa = (Mat_NormalHermitian *)A->data;

271:   PetscFunctionBegin;
272:   *M = Aa->A;
273:   PetscFunctionReturn(PETSC_SUCCESS);
274: }

276: /*@
277:   MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`

279:   Logically Collective

281:   Input Parameter:
282: . A - the `MATNORMALHERMITIAN` matrix

284:   Output Parameter:
285: . M - the matrix object stored inside A

287:   Level: intermediate

289: .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
290: @*/
291: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
292: {
293:   PetscFunctionBegin;
296:   PetscAssertPointer(M, 2);
297:   PetscUseMethod(A, "MatNormalHermitianGetMat_C", (Mat, Mat *), (A, M));
298:   PetscFunctionReturn(PETSC_SUCCESS);
299: }

301: static PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
302: {
303:   Mat_NormalHermitian *Aa = (Mat_NormalHermitian *)A->data;
304:   Mat                  B, conjugate;
305:   PetscInt             m, n, M, N;

307:   PetscFunctionBegin;
308:   PetscCall(MatGetSize(A, &M, &N));
309:   PetscCall(MatGetLocalSize(A, &m, &n));
310:   if (reuse == MAT_REUSE_MATRIX) {
311:     B = *newmat;
312:     PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
313:   } else {
314:     PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
315:     PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
316:     PetscCall(MatProductSetFromOptions(B));
317:     PetscCall(MatProductSymbolic(B));
318:     PetscCall(MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE));
319:   }
320:   if (PetscDefined(USE_COMPLEX)) {
321:     PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate));
322:     PetscCall(MatConjugate(conjugate));
323:     PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B));
324:   }
325:   PetscCall(MatProductNumeric(B));
326:   if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate));
327:   if (reuse == MAT_INPLACE_MATRIX) {
328:     PetscCall(MatHeaderReplace(A, &B));
329:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
330:   PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
331:   PetscFunctionReturn(PETSC_SUCCESS);
332: }

334: #if defined(PETSC_HAVE_HYPRE)
335: static PetscErrorCode MatConvert_NormalHermitian_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
336: {
337:   PetscFunctionBegin;
338:   if (reuse == MAT_INITIAL_MATRIX) {
339:     PetscCall(MatConvert(A, MATAIJ, reuse, B));
340:     PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
341:   } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
342:   PetscFunctionReturn(PETSC_SUCCESS);
343: }
344: #endif

346: /*MC
347:   MATNORMALHERMITIAN - a matrix that behaves like (A*)'*A for `MatMult()` while only containing A

349:   Level: intermediate

351: .seealso: [](ch_matrices), `Mat`, `MatCreateNormalHermitian()`, `MatMult()`, `MatNormalHermitianGetMat()`, `MATNORMAL`, `MatCreateNormal()`
352: M*/

354: /*@
355:   MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like (A*)'*A.

357:   Collective

359:   Input Parameter:
360: . A - the (possibly rectangular complex) matrix

362:   Output Parameter:
363: . N - the matrix that represents (A*)'*A

365:   Level: intermediate

367:   Note:
368:   The product (A*)'*A is NOT actually formed! Rather the new matrix
369:   object performs the matrix-vector product, `MatMult()`, by first multiplying by
370:   A and then (A*)'

372: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
373: @*/
374: PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N)
375: {
376:   PetscInt             m, n;
377:   Mat_NormalHermitian *Na;
378:   VecType              vtype;

380:   PetscFunctionBegin;
381:   PetscCall(MatGetLocalSize(A, &m, &n));
382:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
383:   PetscCall(MatSetSizes(*N, n, n, PETSC_DECIDE, PETSC_DECIDE));
384:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN));
385:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
386:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));

388:   PetscCall(PetscNew(&Na));
389:   (*N)->data = (void *)Na;
390:   PetscCall(PetscObjectReference((PetscObject)A));
391:   Na->A     = A;
392:   Na->scale = 1.0;

394:   PetscCall(MatCreateVecs(A, NULL, &Na->w));

396:   (*N)->ops->destroy           = MatDestroy_NormalHermitian;
397:   (*N)->ops->mult              = MatMult_NormalHermitian;
398:   (*N)->ops->multtranspose     = MatMultHermitianTranspose_Normal;
399:   (*N)->ops->multtransposeadd  = MatMultHermitianTransposeAdd_Normal;
400:   (*N)->ops->multadd           = MatMultHermitianAdd_Normal;
401:   (*N)->ops->getdiagonal       = MatGetDiagonal_NormalHermitian;
402:   (*N)->ops->getdiagonalblock  = MatGetDiagonalBlock_NormalHermitian;
403:   (*N)->ops->scale             = MatScale_NormalHermitian;
404:   (*N)->ops->diagonalscale     = MatDiagonalScale_NormalHermitian;
405:   (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian;
406:   (*N)->ops->permute           = MatPermute_NormalHermitian;
407:   (*N)->ops->duplicate         = MatDuplicate_NormalHermitian;
408:   (*N)->ops->copy              = MatCopy_NormalHermitian;
409:   (*N)->assembled              = PETSC_TRUE;
410:   (*N)->preallocated           = PETSC_TRUE;

412:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatNormalHermitianGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
413:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ));
414:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ));
415: #if defined(PETSC_HAVE_HYPRE)
416:   PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_hypre_C", MatConvert_NormalHermitian_HYPRE));
417: #endif
418:   PetscCall(MatSetOption(*N, MAT_HERMITIAN, PETSC_TRUE));
419:   PetscCall(MatGetVecType(A, &vtype));
420:   PetscCall(MatSetVecType(*N, vtype));
421: #if defined(PETSC_HAVE_DEVICE)
422:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
423: #endif
424:   PetscFunctionReturn(PETSC_SUCCESS);
425: }