Actual source code: normmh.c


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

202: PetscErrorCode MatDestroy_NormalHermitian(Mat N)
203: {
204:   Mat_Normal *Na = (Mat_Normal *)N->data;

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

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

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

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

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

268: PetscErrorCode MatNormalGetMat_NormalHermitian(Mat A, Mat *M)
269: {
270:   Mat_Normal *Aa = (Mat_Normal *)A->data;

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

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

280:    Logically Collective

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

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

288:    Level: intermediate

290: .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
291: @*/
292: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
293: {
294:   PetscFunctionBegin;
298:   PetscUseMethod(A, "MatNormalGetMatHermitian_C", (Mat, Mat *), (A, M));
299:   PetscFunctionReturn(PETSC_SUCCESS);
300: }

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

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

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

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

350:   Level: intermediate

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

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

358:    Collective

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

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

366:    Level: intermediate

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

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

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

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

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

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

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