Actual source code: normm.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_Normal(Mat inA, PetscScalar scale)
 12: {
 13:   Mat_Normal *a = (Mat_Normal *)inA->data;

 15:   a->scale *= scale;
 16:   return 0;
 17: }

 19: PetscErrorCode MatDiagonalScale_Normal(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 MatIncreaseOverlap_Normal(Mat A, PetscInt is_max, IS is[], PetscInt ov)
 43: {
 44:   Mat_Normal *a = (Mat_Normal *)A->data;
 45:   Mat         pattern;

 48:   MatProductCreate(a->A, a->A, NULL, &pattern);
 49:   MatProductSetType(pattern, MATPRODUCT_AtB);
 50:   MatProductSetFromOptions(pattern);
 51:   MatProductSymbolic(pattern);
 52:   MatIncreaseOverlap(pattern, is_max, is, ov);
 53:   MatDestroy(&pattern);
 54:   return 0;
 55: }

 57: PetscErrorCode MatCreateSubMatrices_Normal(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
 58: {
 59:   Mat_Normal *a = (Mat_Normal *)mat->data;
 60:   Mat         B = a->A, *suba;
 61:   IS         *row;
 62:   PetscInt    M;

 65:   if (scall != MAT_REUSE_MATRIX) PetscCalloc1(n, submat);
 66:   MatGetSize(B, &M, NULL);
 67:   PetscMalloc1(n, &row);
 68:   ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]);
 69:   ISSetIdentity(row[0]);
 70:   for (M = 1; M < n; ++M) row[M] = row[0];
 71:   MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba);
 72:   for (M = 0; M < n; ++M) {
 73:     MatCreateNormal(suba[M], *submat + M);
 74:     ((Mat_Normal *)(*submat)[M]->data)->scale = a->scale;
 75:   }
 76:   ISDestroy(&row[0]);
 77:   PetscFree(row);
 78:   MatDestroySubMatrices(n, &suba);
 79:   return 0;
 80: }

 82: PetscErrorCode MatPermute_Normal(Mat A, IS rowp, IS colp, Mat *B)
 83: {
 84:   Mat_Normal *a = (Mat_Normal *)A->data;
 85:   Mat         C, Aa = a->A;
 86:   IS          row;

 89:   ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row);
 90:   ISSetIdentity(row);
 91:   MatPermute(Aa, row, colp, &C);
 92:   ISDestroy(&row);
 93:   MatCreateNormal(C, B);
 94:   MatDestroy(&C);
 95:   return 0;
 96: }

 98: PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B)
 99: {
100:   Mat_Normal *a = (Mat_Normal *)A->data;
101:   Mat         C;

104:   MatDuplicate(a->A, op, &C);
105:   MatCreateNormal(C, B);
106:   MatDestroy(&C);
107:   if (op == MAT_COPY_VALUES) ((Mat_Normal *)(*B)->data)->scale = a->scale;
108:   return 0;
109: }

111: PetscErrorCode MatCopy_Normal(Mat A, Mat B, MatStructure str)
112: {
113:   Mat_Normal *a = (Mat_Normal *)A->data, *b = (Mat_Normal *)B->data;

116:   MatCopy(a->A, b->A, str);
117:   b->scale = a->scale;
118:   VecDestroy(&b->left);
119:   VecDestroy(&b->right);
120:   VecDestroy(&b->leftwork);
121:   VecDestroy(&b->rightwork);
122:   return 0;
123: }

125: PetscErrorCode MatMult_Normal(Mat N, Vec x, Vec y)
126: {
127:   Mat_Normal *Na = (Mat_Normal *)N->data;
128:   Vec         in;

130:   in = x;
131:   if (Na->right) {
132:     if (!Na->rightwork) VecDuplicate(Na->right, &Na->rightwork);
133:     VecPointwiseMult(Na->rightwork, Na->right, in);
134:     in = Na->rightwork;
135:   }
136:   MatMult(Na->A, in, Na->w);
137:   MatMultTranspose(Na->A, Na->w, y);
138:   if (Na->left) VecPointwiseMult(y, Na->left, y);
139:   VecScale(y, Na->scale);
140:   return 0;
141: }

143: PetscErrorCode MatMultAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
144: {
145:   Mat_Normal *Na = (Mat_Normal *)N->data;
146:   Vec         in;

148:   in = v1;
149:   if (Na->right) {
150:     if (!Na->rightwork) VecDuplicate(Na->right, &Na->rightwork);
151:     VecPointwiseMult(Na->rightwork, Na->right, in);
152:     in = Na->rightwork;
153:   }
154:   MatMult(Na->A, in, Na->w);
155:   VecScale(Na->w, Na->scale);
156:   if (Na->left) {
157:     MatMultTranspose(Na->A, Na->w, v3);
158:     VecPointwiseMult(v3, Na->left, v3);
159:     VecAXPY(v3, 1.0, v2);
160:   } else {
161:     MatMultTransposeAdd(Na->A, Na->w, v2, v3);
162:   }
163:   return 0;
164: }

166: PetscErrorCode MatMultTranspose_Normal(Mat N, Vec x, Vec y)
167: {
168:   Mat_Normal *Na = (Mat_Normal *)N->data;
169:   Vec         in;

171:   in = x;
172:   if (Na->left) {
173:     if (!Na->leftwork) VecDuplicate(Na->left, &Na->leftwork);
174:     VecPointwiseMult(Na->leftwork, Na->left, in);
175:     in = Na->leftwork;
176:   }
177:   MatMult(Na->A, in, Na->w);
178:   MatMultTranspose(Na->A, Na->w, y);
179:   if (Na->right) VecPointwiseMult(y, Na->right, y);
180:   VecScale(y, Na->scale);
181:   return 0;
182: }

184: PetscErrorCode MatMultTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
185: {
186:   Mat_Normal *Na = (Mat_Normal *)N->data;
187:   Vec         in;

189:   in = v1;
190:   if (Na->left) {
191:     if (!Na->leftwork) VecDuplicate(Na->left, &Na->leftwork);
192:     VecPointwiseMult(Na->leftwork, Na->left, in);
193:     in = Na->leftwork;
194:   }
195:   MatMult(Na->A, in, Na->w);
196:   VecScale(Na->w, Na->scale);
197:   if (Na->right) {
198:     MatMultTranspose(Na->A, Na->w, v3);
199:     VecPointwiseMult(v3, Na->right, v3);
200:     VecAXPY(v3, 1.0, v2);
201:   } else {
202:     MatMultTransposeAdd(Na->A, Na->w, v2, v3);
203:   }
204:   return 0;
205: }

207: PetscErrorCode MatDestroy_Normal(Mat N)
208: {
209:   Mat_Normal *Na = (Mat_Normal *)N->data;

211:   MatDestroy(&Na->A);
212:   MatDestroy(&Na->D);
213:   VecDestroy(&Na->w);
214:   VecDestroy(&Na->left);
215:   VecDestroy(&Na->right);
216:   VecDestroy(&Na->leftwork);
217:   VecDestroy(&Na->rightwork);
218:   PetscFree(N->data);
219:   PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL);
220:   PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_seqaij_C", NULL);
221:   PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_mpiaij_C", NULL);
222:   PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_seqdense_C", NULL);
223:   PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_mpidense_C", NULL);
224:   PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_dense_C", NULL);
225:   return 0;
226: }

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

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

259: PetscErrorCode MatGetDiagonalBlock_Normal(Mat N, Mat *D)
260: {
261:   Mat_Normal *Na = (Mat_Normal *)N->data;
262:   Mat         M, A = Na->A;

264:   MatGetDiagonalBlock(A, &M);
265:   MatCreateNormal(M, &Na->D);
266:   *D = Na->D;
267:   return 0;
268: }

270: PetscErrorCode MatNormalGetMat_Normal(Mat A, Mat *M)
271: {
272:   Mat_Normal *Aa = (Mat_Normal *)A->data;

274:   *M = Aa->A;
275:   return 0;
276: }

278: /*@
279:       MatNormalGetMat - Gets the `Mat` object stored inside a `MATNORMAL`

281:    Logically collective on A

283:    Input Parameter:
284: .   A  - the `MATNORMAL` matrix

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

289:    Level: intermediate

291: .seealso: `MATNORMAL`, `MATNORMALHERMITIAN`, `MatCreateNormal()`
292: @*/
293: PetscErrorCode MatNormalGetMat(Mat A, Mat *M)
294: {
298:   PetscUseMethod(A, "MatNormalGetMat_C", (Mat, Mat *), (A, M));
299:   return 0;
300: }

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

308:   MatGetSize(A, &M, &N);
309:   MatGetLocalSize(A, &m, &n);
310:   if (reuse == MAT_REUSE_MATRIX) {
311:     B = *newmat;
312:     MatProductReplaceMats(Aa->A, Aa->A, NULL, B);
313:   } else {
314:     MatProductCreate(Aa->A, Aa->A, NULL, &B);
315:     MatProductSetType(B, MATPRODUCT_AtB);
316:     MatProductSetFromOptions(B);
317:     MatProductSymbolic(B);
318:     MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE);
319:   }
320:   MatProductNumeric(B);
321:   if (reuse == MAT_INPLACE_MATRIX) {
322:     MatHeaderReplace(A, &B);
323:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
324:   MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat);
325:   return 0;
326: }

328: typedef struct {
329:   Mat work[2];
330: } Normal_Dense;

332: PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
333: {
334:   Mat           A, B;
335:   Normal_Dense *contents;
336:   Mat_Normal   *a;
337:   PetscScalar  *array;

339:   MatCheckProduct(C, 1);
340:   A        = C->product->A;
341:   a        = (Mat_Normal *)A->data;
342:   B        = C->product->B;
343:   contents = (Normal_Dense *)C->product->data;
345:   if (a->right) {
346:     MatCopy(B, C, SAME_NONZERO_PATTERN);
347:     MatDiagonalScale(C, a->right, NULL);
348:   }
349:   MatProductNumeric(contents->work[0]);
350:   MatDenseGetArrayWrite(C, &array);
351:   MatDensePlaceArray(contents->work[1], array);
352:   MatProductNumeric(contents->work[1]);
353:   MatDenseRestoreArrayWrite(C, &array);
354:   MatDenseResetArray(contents->work[1]);
355:   MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
356:   MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
357:   MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
358:   MatScale(C, a->scale);
359:   return 0;
360: }

362: PetscErrorCode MatNormal_DenseDestroy(void *ctx)
363: {
364:   Normal_Dense *contents = (Normal_Dense *)ctx;

366:   MatDestroy(contents->work);
367:   MatDestroy(contents->work + 1);
368:   PetscFree(contents);
369:   return 0;
370: }

372: PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
373: {
374:   Mat           A, B;
375:   Normal_Dense *contents = NULL;
376:   Mat_Normal   *a;
377:   PetscScalar  *array;
378:   PetscInt      n, N, m, M;

380:   MatCheckProduct(C, 1);
382:   A = C->product->A;
383:   a = (Mat_Normal *)A->data;
385:   B = C->product->B;
386:   MatGetLocalSize(C, &m, &n);
387:   MatGetSize(C, &M, &N);
388:   if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
389:     MatGetLocalSize(B, NULL, &n);
390:     MatGetSize(B, NULL, &N);
391:     MatGetLocalSize(A, &m, NULL);
392:     MatGetSize(A, &M, NULL);
393:     MatSetSizes(C, m, n, M, N);
394:   }
395:   MatSetType(C, ((PetscObject)B)->type_name);
396:   MatSetUp(C);
397:   PetscNew(&contents);
398:   C->product->data    = contents;
399:   C->product->destroy = MatNormal_DenseDestroy;
400:   if (a->right) {
401:     MatProductCreate(a->A, C, NULL, contents->work);
402:   } else {
403:     MatProductCreate(a->A, B, NULL, contents->work);
404:   }
405:   MatProductSetType(contents->work[0], MATPRODUCT_AB);
406:   MatProductSetFromOptions(contents->work[0]);
407:   MatProductSymbolic(contents->work[0]);
408:   MatProductCreate(a->A, contents->work[0], NULL, contents->work + 1);
409:   MatProductSetType(contents->work[1], MATPRODUCT_AtB);
410:   MatProductSetFromOptions(contents->work[1]);
411:   MatProductSymbolic(contents->work[1]);
412:   MatDenseGetArrayWrite(C, &array);
413:   MatSeqDenseSetPreallocation(contents->work[1], array);
414:   MatMPIDenseSetPreallocation(contents->work[1], array);
415:   MatDenseRestoreArrayWrite(C, &array);
416:   C->ops->productnumeric = MatProductNumeric_Normal_Dense;
417:   return 0;
418: }

420: PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C)
421: {
422:   C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
423:   return 0;
424: }

426: PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
427: {
428:   Mat_Product *product = C->product;

430:   if (product->type == MATPRODUCT_AB) MatProductSetFromOptions_Normal_Dense_AB(C);
431:   return 0;
432: }

434: /*@
435:       MatCreateNormal - Creates a new `MATNORMAL` matrix object that behaves like A'*A.

437:    Collective

439:    Input Parameter:
440: .   A  - the (possibly rectangular) matrix

442:    Output Parameter:
443: .   N - the matrix that represents A'*A

445:    Level: intermediate

447:    Notes:
448:     The product A'*A is NOT actually formed! Rather the new matrix
449:           object performs the matrix-vector product, `MatMult()`, by first multiplying by
450:           A and then A'

452: .seealso: `MATNORMAL`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`,
453: @*/
454: PetscErrorCode MatCreateNormal(Mat A, Mat *N)
455: {
456:   PetscInt    n, nn;
457:   Mat_Normal *Na;
458:   VecType     vtype;

460:   MatGetSize(A, NULL, &nn);
461:   MatGetLocalSize(A, NULL, &n);
462:   MatCreate(PetscObjectComm((PetscObject)A), N);
463:   MatSetSizes(*N, n, n, nn, nn);
464:   PetscObjectChangeTypeName((PetscObject)*N, MATNORMAL);
465:   PetscLayoutReference(A->cmap, &(*N)->rmap);
466:   PetscLayoutReference(A->cmap, &(*N)->cmap);

468:   PetscNew(&Na);
469:   (*N)->data = (void *)Na;
470:   PetscObjectReference((PetscObject)A);
471:   Na->A     = A;
472:   Na->scale = 1.0;

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

476:   (*N)->ops->destroy           = MatDestroy_Normal;
477:   (*N)->ops->mult              = MatMult_Normal;
478:   (*N)->ops->multtranspose     = MatMultTranspose_Normal;
479:   (*N)->ops->multtransposeadd  = MatMultTransposeAdd_Normal;
480:   (*N)->ops->multadd           = MatMultAdd_Normal;
481:   (*N)->ops->getdiagonal       = MatGetDiagonal_Normal;
482:   (*N)->ops->getdiagonalblock  = MatGetDiagonalBlock_Normal;
483:   (*N)->ops->scale             = MatScale_Normal;
484:   (*N)->ops->diagonalscale     = MatDiagonalScale_Normal;
485:   (*N)->ops->increaseoverlap   = MatIncreaseOverlap_Normal;
486:   (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
487:   (*N)->ops->permute           = MatPermute_Normal;
488:   (*N)->ops->duplicate         = MatDuplicate_Normal;
489:   (*N)->ops->copy              = MatCopy_Normal;
490:   (*N)->assembled              = PETSC_TRUE;
491:   (*N)->preallocated           = PETSC_TRUE;

493:   PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMat_C", MatNormalGetMat_Normal);
494:   PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_seqaij_C", MatConvert_Normal_AIJ);
495:   PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_mpiaij_C", MatConvert_Normal_AIJ);
496:   PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_seqdense_C", MatProductSetFromOptions_Normal_Dense);
497:   PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_mpidense_C", MatProductSetFromOptions_Normal_Dense);
498:   PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_dense_C", MatProductSetFromOptions_Normal_Dense);
499:   MatSetOption(*N, MAT_SYMMETRIC, PETSC_TRUE);
500:   MatGetVecType(A, &vtype);
501:   MatSetVecType(*N, vtype);
502: #if defined(PETSC_HAVE_DEVICE)
503:   MatBindToCPU(*N, A->boundtocpu);
504: #endif
505:   return 0;
506: }