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:   a->scale *= scale;
 16:   return 0;
 17: }

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

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

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

 50:   if (scall != MAT_REUSE_MATRIX) PetscCalloc1(n, submat);
 51:   MatGetSize(B, &M, NULL);
 52:   PetscMalloc1(n, &row);
 53:   ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]);
 54:   ISSetIdentity(row[0]);
 55:   for (M = 1; M < n; ++M) row[M] = row[0];
 56:   MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba);
 57:   for (M = 0; M < n; ++M) {
 58:     MatCreateNormalHermitian(suba[M], *submat + M);
 59:     ((Mat_Normal *)(*submat)[M]->data)->scale = a->scale;
 60:   }
 61:   ISDestroy(&row[0]);
 62:   PetscFree(row);
 63:   MatDestroySubMatrices(n, &suba);
 64:   return 0;
 65: }

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

 74:   ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row);
 75:   ISSetIdentity(row);
 76:   MatPermute(Aa, row, colp, &C);
 77:   ISDestroy(&row);
 78:   MatCreateNormalHermitian(C, B);
 79:   MatDestroy(&C);
 80:   return 0;
 81: }

 83: PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B)
 84: {
 85:   Mat_Normal *a = (Mat_Normal *)A->data;
 86:   Mat         C;

 89:   MatDuplicate(a->A, op, &C);
 90:   MatCreateNormalHermitian(C, B);
 91:   MatDestroy(&C);
 92:   if (op == MAT_COPY_VALUES) ((Mat_Normal *)(*B)->data)->scale = a->scale;
 93:   return 0;
 94: }

 96: PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str)
 97: {
 98:   Mat_Normal *a = (Mat_Normal *)A->data, *b = (Mat_Normal *)B->data;

101:   MatCopy(a->A, b->A, str);
102:   b->scale = a->scale;
103:   VecDestroy(&b->left);
104:   VecDestroy(&b->right);
105:   VecDestroy(&b->leftwork);
106:   VecDestroy(&b->rightwork);
107:   return 0;
108: }

110: PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y)
111: {
112:   Mat_Normal *Na = (Mat_Normal *)N->data;
113:   Vec         in;

115:   in = x;
116:   if (Na->right) {
117:     if (!Na->rightwork) VecDuplicate(Na->right, &Na->rightwork);
118:     VecPointwiseMult(Na->rightwork, Na->right, in);
119:     in = Na->rightwork;
120:   }
121:   MatMult(Na->A, in, Na->w);
122:   MatMultHermitianTranspose(Na->A, Na->w, y);
123:   if (Na->left) VecPointwiseMult(y, Na->left, y);
124:   VecScale(y, Na->scale);
125:   return 0;
126: }

128: PetscErrorCode MatMultHermitianAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
129: {
130:   Mat_Normal *Na = (Mat_Normal *)N->data;
131:   Vec         in;

133:   in = v1;
134:   if (Na->right) {
135:     if (!Na->rightwork) VecDuplicate(Na->right, &Na->rightwork);
136:     VecPointwiseMult(Na->rightwork, Na->right, in);
137:     in = Na->rightwork;
138:   }
139:   MatMult(Na->A, in, Na->w);
140:   VecScale(Na->w, Na->scale);
141:   if (Na->left) {
142:     MatMultHermitianTranspose(Na->A, Na->w, v3);
143:     VecPointwiseMult(v3, Na->left, v3);
144:     VecAXPY(v3, 1.0, v2);
145:   } else {
146:     MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3);
147:   }
148:   return 0;
149: }

151: PetscErrorCode MatMultHermitianTranspose_Normal(Mat N, Vec x, Vec y)
152: {
153:   Mat_Normal *Na = (Mat_Normal *)N->data;
154:   Vec         in;

156:   in = x;
157:   if (Na->left) {
158:     if (!Na->leftwork) VecDuplicate(Na->left, &Na->leftwork);
159:     VecPointwiseMult(Na->leftwork, Na->left, in);
160:     in = Na->leftwork;
161:   }
162:   MatMult(Na->A, in, Na->w);
163:   MatMultHermitianTranspose(Na->A, Na->w, y);
164:   if (Na->right) VecPointwiseMult(y, Na->right, y);
165:   VecScale(y, Na->scale);
166:   return 0;
167: }

169: PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
170: {
171:   Mat_Normal *Na = (Mat_Normal *)N->data;
172:   Vec         in;

174:   in = v1;
175:   if (Na->left) {
176:     if (!Na->leftwork) VecDuplicate(Na->left, &Na->leftwork);
177:     VecPointwiseMult(Na->leftwork, Na->left, in);
178:     in = Na->leftwork;
179:   }
180:   MatMult(Na->A, in, Na->w);
181:   VecScale(Na->w, Na->scale);
182:   if (Na->right) {
183:     MatMultHermitianTranspose(Na->A, Na->w, v3);
184:     VecPointwiseMult(v3, Na->right, v3);
185:     VecAXPY(v3, 1.0, v2);
186:   } else {
187:     MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3);
188:   }
189:   return 0;
190: }

192: PetscErrorCode MatDestroy_NormalHermitian(Mat N)
193: {
194:   Mat_Normal *Na = (Mat_Normal *)N->data;

196:   MatDestroy(&Na->A);
197:   MatDestroy(&Na->D);
198:   VecDestroy(&Na->w);
199:   VecDestroy(&Na->left);
200:   VecDestroy(&Na->right);
201:   VecDestroy(&Na->leftwork);
202:   VecDestroy(&Na->rightwork);
203:   PetscFree(N->data);
204:   PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMatHermitian_C", NULL);
205:   PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL);
206:   PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL);
207:   return 0;
208: }

210: /*
211:       Slow, nonscalable version
212: */
213: PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
214: {
215:   Mat_Normal        *Na = (Mat_Normal *)N->data;
216:   Mat                A  = Na->A;
217:   PetscInt           i, j, rstart, rend, nnz;
218:   const PetscInt    *cols;
219:   PetscScalar       *diag, *work, *values;
220:   const PetscScalar *mvalues;

222:   PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work);
223:   PetscArrayzero(work, A->cmap->N);
224:   MatGetOwnershipRange(A, &rstart, &rend);
225:   for (i = rstart; i < rend; i++) {
226:     MatGetRow(A, i, &nnz, &cols, &mvalues);
227:     for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
228:     MatRestoreRow(A, i, &nnz, &cols, &mvalues);
229:   }
230:   MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N));
231:   rstart = N->cmap->rstart;
232:   rend   = N->cmap->rend;
233:   VecGetArray(v, &values);
234:   PetscArraycpy(values, diag + rstart, rend - rstart);
235:   VecRestoreArray(v, &values);
236:   PetscFree2(diag, work);
237:   VecScale(v, Na->scale);
238:   return 0;
239: }

241: PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
242: {
243:   Mat_Normal *Na = (Mat_Normal *)N->data;
244:   Mat         M, A = Na->A;

246:   MatGetDiagonalBlock(A, &M);
247:   MatCreateNormalHermitian(M, &Na->D);
248:   *D = Na->D;
249:   return 0;
250: }

252: PetscErrorCode MatNormalGetMat_NormalHermitian(Mat A, Mat *M)
253: {
254:   Mat_Normal *Aa = (Mat_Normal *)A->data;

256:   *M = Aa->A;
257:   return 0;
258: }

260: /*@
261:       MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`

263:    Logically collective

265:    Input Parameter:
266: .   A  - the `MATNORMALHERMITIAN` matrix

268:    Output Parameter:
269: .   M - the matrix object stored inside A

271:    Level: intermediate

273: .seealso: `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
274: @*/
275: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
276: {
280:   PetscUseMethod(A, "MatNormalGetMatHermitian_C", (Mat, Mat *), (A, M));
281:   return 0;
282: }

284: PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
285: {
286:   Mat_Normal *Aa = (Mat_Normal *)A->data;
287:   Mat         B, conjugate;
288:   PetscInt    m, n, M, N;

290:   MatGetSize(A, &M, &N);
291:   MatGetLocalSize(A, &m, &n);
292:   if (reuse == MAT_REUSE_MATRIX) {
293:     B = *newmat;
294:     MatProductReplaceMats(Aa->A, Aa->A, NULL, B);
295:   } else {
296:     MatProductCreate(Aa->A, Aa->A, NULL, &B);
297:     MatProductSetType(B, MATPRODUCT_AtB);
298:     MatProductSetFromOptions(B);
299:     MatProductSymbolic(B);
300:     MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE);
301:   }
302:   if (PetscDefined(USE_COMPLEX)) {
303:     MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate);
304:     MatConjugate(conjugate);
305:     MatProductReplaceMats(conjugate, Aa->A, NULL, B);
306:   }
307:   MatProductNumeric(B);
308:   if (PetscDefined(USE_COMPLEX)) MatDestroy(&conjugate);
309:   if (reuse == MAT_INPLACE_MATRIX) {
310:     MatHeaderReplace(A, &B);
311:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
312:   MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat);
313:   return 0;
314: }

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

319:    Collective

321:    Input Parameter:
322: .   A  - the (possibly rectangular complex) matrix

324:    Output Parameter:
325: .   N - the matrix that represents (A*)'*A

327:    Level: intermediate

329:    Note:
330:     The product (A*)'*A is NOT actually formed! Rather the new matrix
331:           object performs the matrix-vector product, `MatMult()`, by first multiplying by
332:           A and then (A*)'

334: .seealso: `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
335: @*/
336: PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N)
337: {
338:   PetscInt    m, n;
339:   Mat_Normal *Na;
340:   VecType     vtype;

342:   MatGetLocalSize(A, &m, &n);
343:   MatCreate(PetscObjectComm((PetscObject)A), N);
344:   MatSetSizes(*N, n, n, PETSC_DECIDE, PETSC_DECIDE);
345:   PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN);
346:   PetscLayoutReference(A->cmap, &(*N)->rmap);
347:   PetscLayoutReference(A->cmap, &(*N)->cmap);

349:   PetscNew(&Na);
350:   (*N)->data = (void *)Na;
351:   PetscObjectReference((PetscObject)A);
352:   Na->A     = A;
353:   Na->scale = 1.0;

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

357:   (*N)->ops->destroy           = MatDestroy_NormalHermitian;
358:   (*N)->ops->mult              = MatMult_NormalHermitian;
359:   (*N)->ops->multtranspose     = MatMultHermitianTranspose_Normal;
360:   (*N)->ops->multtransposeadd  = MatMultHermitianTransposeAdd_Normal;
361:   (*N)->ops->multadd           = MatMultHermitianAdd_Normal;
362:   (*N)->ops->getdiagonal       = MatGetDiagonal_NormalHermitian;
363:   (*N)->ops->getdiagonalblock  = MatGetDiagonalBlock_NormalHermitian;
364:   (*N)->ops->scale             = MatScale_NormalHermitian;
365:   (*N)->ops->diagonalscale     = MatDiagonalScale_NormalHermitian;
366:   (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian;
367:   (*N)->ops->permute           = MatPermute_NormalHermitian;
368:   (*N)->ops->duplicate         = MatDuplicate_NormalHermitian;
369:   (*N)->ops->copy              = MatCopy_NormalHermitian;
370:   (*N)->assembled              = PETSC_TRUE;
371:   (*N)->preallocated           = PETSC_TRUE;

373:   PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMatHermitian_C", MatNormalGetMat_NormalHermitian);
374:   PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ);
375:   PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ);
376:   MatSetOption(*N, MAT_HERMITIAN, PETSC_TRUE);
377:   MatGetVecType(A, &vtype);
378:   MatSetVecType(*N, vtype);
379: #if defined(PETSC_HAVE_DEVICE)
380:   MatBindToCPU(*N, A->boundtocpu);
381: #endif
382:   return 0;
383: }