Actual source code: svd.c

      1: #include <petsc/private/pcimpl.h>
      2: #include <petscblaslapack.h>

      4: /*
      5:    Private context (data structure) for the SVD preconditioner.
      6: */
      7: typedef struct {
      8:   Vec               diag, work;
      9:   Mat               A, U, Vt;
     10:   PetscInt          nzero;
     11:   PetscReal         zerosing; /* measure of smallest singular value treated as nonzero */
     12:   PetscInt          essrank;  /* essential rank of operator */
     13:   VecScatter        left2red, right2red;
     14:   Vec               leftred, rightred;
     15:   PetscViewer       monitor;
     16:   PetscViewerFormat monitorformat;
     17: } PC_SVD;

     19: typedef enum {
     20:   READ       = 1,
     21:   WRITE      = 2,
     22:   READ_WRITE = 3
     23: } AccessMode;

     25: /*
     26:    PCSetUp_SVD - Prepares for the use of the SVD preconditioner
     27:                     by setting data structures and options.

     29:    Input Parameter:
     30: .  pc - the preconditioner context

     32:    Application Interface Routine: PCSetUp()

     34:    Note:
     35:    The interface routine PCSetUp() is not usually called directly by
     36:    the user, but instead is called by PCApply() if necessary.
     37: */
     38: static PetscErrorCode PCSetUp_SVD(PC pc)
     39: {
     40:   PC_SVD      *jac = (PC_SVD *)pc->data;
     41:   PetscScalar *a, *u, *v, *d, *work;
     42:   PetscBLASInt nb, lwork;
     43:   PetscInt     i, n;
     44:   PetscMPIInt  size;

     46:   PetscFunctionBegin;
     47:   PetscCall(MatDestroy(&jac->A));
     48:   PetscCallMPI(MPI_Comm_size(((PetscObject)pc->pmat)->comm, &size));
     49:   if (size > 1) {
     50:     Mat redmat;

     52:     PetscCall(MatCreateRedundantMatrix(pc->pmat, size, PETSC_COMM_SELF, MAT_INITIAL_MATRIX, &redmat));
     53:     PetscCall(MatConvert(redmat, MATSEQDENSE, MAT_INITIAL_MATRIX, &jac->A));
     54:     PetscCall(MatDestroy(&redmat));
     55:   } else {
     56:     PetscCall(MatConvert(pc->pmat, MATSEQDENSE, MAT_INITIAL_MATRIX, &jac->A));
     57:   }
     58:   if (!jac->diag) { /* assume square matrices */
     59:     PetscCall(MatCreateVecs(jac->A, &jac->diag, &jac->work));
     60:   }
     61:   if (!jac->U) {
     62:     PetscCall(MatDuplicate(jac->A, MAT_DO_NOT_COPY_VALUES, &jac->U));
     63:     PetscCall(MatDuplicate(jac->A, MAT_DO_NOT_COPY_VALUES, &jac->Vt));
     64:   }
     65:   PetscCall(MatGetSize(jac->A, &n, NULL));
     66:   if (!n) PetscFunctionReturn(PETSC_SUCCESS);

     68:   PetscCall(PetscBLASIntCast(n, &nb));
     69:   lwork = 5 * nb;
     70:   PetscCall(PetscMalloc1(lwork, &work));
     71:   PetscCall(MatDenseGetArray(jac->A, &a));
     72:   PetscCall(MatDenseGetArray(jac->U, &u));
     73:   PetscCall(MatDenseGetArray(jac->Vt, &v));
     74:   PetscCall(VecGetArray(jac->diag, &d));
     75: #if !defined(PETSC_USE_COMPLEX)
     76:   {
     77:     PetscBLASInt lierr;
     78:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     79:     PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &nb, &nb, a, &nb, d, u, &nb, v, &nb, work, &lwork, &lierr));
     80:     PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "gesvd() error %" PetscBLASInt_FMT, lierr);
     81:     PetscCall(PetscFPTrapPop());
     82:   }
     83: #else
     84:   {
     85:     PetscBLASInt lierr;
     86:     PetscReal   *rwork, *dd;
     87:     PetscCall(PetscMalloc1(5 * nb, &rwork));
     88:     PetscCall(PetscMalloc1(nb, &dd));
     89:     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     90:     PetscCallBLAS("LAPACKgesvd", LAPACKgesvd_("A", "A", &nb, &nb, a, &nb, dd, u, &nb, v, &nb, work, &lwork, rwork, &lierr));
     91:     PetscCheck(!lierr, PETSC_COMM_SELF, PETSC_ERR_LIB, "gesv() error %" PetscBLASInt_FMT, lierr);
     92:     PetscCall(PetscFree(rwork));
     93:     for (i = 0; i < n; i++) d[i] = dd[i];
     94:     PetscCall(PetscFree(dd));
     95:     PetscCall(PetscFPTrapPop());
     96:   }
     97: #endif
     98:   PetscCall(MatDenseRestoreArray(jac->A, &a));
     99:   PetscCall(MatDenseRestoreArray(jac->U, &u));
    100:   PetscCall(MatDenseRestoreArray(jac->Vt, &v));
    101:   for (i = n - 1; i >= 0; i--)
    102:     if (PetscRealPart(d[i]) > jac->zerosing) break;
    103:   jac->nzero = n - 1 - i;
    104:   if (jac->monitor) {
    105:     PetscCall(PetscViewerASCIIAddTab(jac->monitor, ((PetscObject)pc)->tablevel));
    106:     PetscCall(PetscViewerASCIIPrintf(jac->monitor, "    SVD: condition number %14.12e, %" PetscInt_FMT " of %" PetscInt_FMT " singular values are (nearly) zero\n", (double)PetscRealPart(d[0] / d[n - 1]), jac->nzero, n));
    107:     if (n < 10 || jac->monitorformat == PETSC_VIEWER_ALL) {
    108:       PetscCall(PetscViewerASCIIPrintf(jac->monitor, "    SVD: singular values:\n"));
    109:       for (i = 0; i < n; i++) {
    110:         if (i % 5 == 0) {
    111:           if (i != 0) PetscCall(PetscViewerASCIIPrintf(jac->monitor, "\n"));
    112:           PetscCall(PetscViewerASCIIPrintf(jac->monitor, "        "));
    113:         }
    114:         PetscCall(PetscViewerASCIIPrintf(jac->monitor, " %14.12e", (double)PetscRealPart(d[i])));
    115:       }
    116:       PetscCall(PetscViewerASCIIPrintf(jac->monitor, "\n"));
    117:     } else { /* print 5 smallest and 5 largest */
    118:       PetscCall(PetscViewerASCIIPrintf(jac->monitor, "    SVD: smallest singular values: %14.12e %14.12e %14.12e %14.12e %14.12e\n", (double)PetscRealPart(d[n - 1]), (double)PetscRealPart(d[n - 2]), (double)PetscRealPart(d[n - 3]), (double)PetscRealPart(d[n - 4]), (double)PetscRealPart(d[n - 5])));
    119:       PetscCall(PetscViewerASCIIPrintf(jac->monitor, "    SVD: largest singular values : %14.12e %14.12e %14.12e %14.12e %14.12e\n", (double)PetscRealPart(d[4]), (double)PetscRealPart(d[3]), (double)PetscRealPart(d[2]), (double)PetscRealPart(d[1]), (double)PetscRealPart(d[0])));
    120:     }
    121:     PetscCall(PetscViewerASCIISubtractTab(jac->monitor, ((PetscObject)pc)->tablevel));
    122:   }
    123:   PetscCall(PetscInfo(pc, "Largest and smallest singular values %14.12e %14.12e\n", (double)PetscRealPart(d[0]), (double)PetscRealPart(d[n - 1])));
    124:   for (i = 0; i < n - jac->nzero; i++) d[i] = 1.0 / d[i];
    125:   for (; i < n; i++) d[i] = 0.0;
    126:   if (jac->essrank > 0)
    127:     for (i = 0; i < n - jac->nzero - jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
    128:   PetscCall(PetscInfo(pc, "Number of zero or nearly singular values %" PetscInt_FMT "\n", jac->nzero));
    129:   PetscCall(VecRestoreArray(jac->diag, &d));
    130:   PetscCall(PetscFree(work));
    131:   PetscFunctionReturn(PETSC_SUCCESS);
    132: }

    134: static PetscErrorCode PCSVDGetVec(PC pc, PCSide side, AccessMode amode, Vec x, Vec *xred)
    135: {
    136:   PC_SVD     *jac = (PC_SVD *)pc->data;
    137:   PetscMPIInt size;

    139:   PetscFunctionBegin;
    140:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
    141:   *xred = NULL;
    142:   switch (side) {
    143:   case PC_LEFT:
    144:     if (size == 1) *xred = x;
    145:     else {
    146:       if (!jac->left2red) PetscCall(VecScatterCreateToAll(x, &jac->left2red, &jac->leftred));
    147:       if (amode & READ) {
    148:         PetscCall(VecScatterBegin(jac->left2red, x, jac->leftred, INSERT_VALUES, SCATTER_FORWARD));
    149:         PetscCall(VecScatterEnd(jac->left2red, x, jac->leftred, INSERT_VALUES, SCATTER_FORWARD));
    150:       }
    151:       *xred = jac->leftred;
    152:     }
    153:     break;
    154:   case PC_RIGHT:
    155:     if (size == 1) *xred = x;
    156:     else {
    157:       if (!jac->right2red) PetscCall(VecScatterCreateToAll(x, &jac->right2red, &jac->rightred));
    158:       if (amode & READ) {
    159:         PetscCall(VecScatterBegin(jac->right2red, x, jac->rightred, INSERT_VALUES, SCATTER_FORWARD));
    160:         PetscCall(VecScatterEnd(jac->right2red, x, jac->rightred, INSERT_VALUES, SCATTER_FORWARD));
    161:       }
    162:       *xred = jac->rightred;
    163:     }
    164:     break;
    165:   default:
    166:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Side must be LEFT or RIGHT");
    167:   }
    168:   PetscFunctionReturn(PETSC_SUCCESS);
    169: }

    171: static PetscErrorCode PCSVDRestoreVec(PC pc, PCSide side, AccessMode amode, Vec x, Vec *xred)
    172: {
    173:   PC_SVD     *jac = (PC_SVD *)pc->data;
    174:   PetscMPIInt size;

    176:   PetscFunctionBegin;
    177:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
    178:   switch (side) {
    179:   case PC_LEFT:
    180:     if (size != 1 && amode & WRITE) {
    181:       PetscCall(VecScatterBegin(jac->left2red, jac->leftred, x, INSERT_VALUES, SCATTER_REVERSE));
    182:       PetscCall(VecScatterEnd(jac->left2red, jac->leftred, x, INSERT_VALUES, SCATTER_REVERSE));
    183:     }
    184:     break;
    185:   case PC_RIGHT:
    186:     if (size != 1 && amode & WRITE) {
    187:       PetscCall(VecScatterBegin(jac->right2red, jac->rightred, x, INSERT_VALUES, SCATTER_REVERSE));
    188:       PetscCall(VecScatterEnd(jac->right2red, jac->rightred, x, INSERT_VALUES, SCATTER_REVERSE));
    189:     }
    190:     break;
    191:   default:
    192:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_PLIB, "Side must be LEFT or RIGHT");
    193:   }
    194:   *xred = NULL;
    195:   PetscFunctionReturn(PETSC_SUCCESS);
    196: }

    198: /*
    199:    PCApply_SVD - Applies the SVD preconditioner to a vector.

    201:    Input Parameters:
    202: .  pc - the preconditioner context
    203: .  x - input vector

    205:    Output Parameter:
    206: .  y - output vector

    208:    Application Interface Routine: PCApply()
    209:  */
    210: static PetscErrorCode PCApply_SVD(PC pc, Vec x, Vec y)
    211: {
    212:   PC_SVD *jac  = (PC_SVD *)pc->data;
    213:   Vec     work = jac->work, xred, yred;

    215:   PetscFunctionBegin;
    216:   PetscCall(PCSVDGetVec(pc, PC_RIGHT, READ, x, &xred));
    217:   PetscCall(PCSVDGetVec(pc, PC_LEFT, WRITE, y, &yred));
    218: #if !defined(PETSC_USE_COMPLEX)
    219:   PetscCall(MatMultTranspose(jac->U, xred, work));
    220: #else
    221:   PetscCall(MatMultHermitianTranspose(jac->U, xred, work));
    222: #endif
    223:   PetscCall(VecPointwiseMult(work, work, jac->diag));
    224: #if !defined(PETSC_USE_COMPLEX)
    225:   PetscCall(MatMultTranspose(jac->Vt, work, yred));
    226: #else
    227:   PetscCall(MatMultHermitianTranspose(jac->Vt, work, yred));
    228: #endif
    229:   PetscCall(PCSVDRestoreVec(pc, PC_RIGHT, READ, x, &xred));
    230:   PetscCall(PCSVDRestoreVec(pc, PC_LEFT, WRITE, y, &yred));
    231:   PetscFunctionReturn(PETSC_SUCCESS);
    232: }

    234: static PetscErrorCode PCMatApply_SVD(PC pc, Mat X, Mat Y)
    235: {
    236:   PC_SVD *jac = (PC_SVD *)pc->data;
    237:   Mat     W;

    239:   PetscFunctionBegin;
    240:   PetscCall(MatTransposeMatMult(jac->U, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &W));
    241:   PetscCall(MatDiagonalScale(W, jac->diag, NULL));
    242:   PetscCall(MatTransposeMatMult(jac->Vt, W, MAT_REUSE_MATRIX, PETSC_DETERMINE, &Y));
    243:   PetscCall(MatDestroy(&W));
    244:   PetscFunctionReturn(PETSC_SUCCESS);
    245: }

    247: static PetscErrorCode PCApplyTranspose_SVD(PC pc, Vec x, Vec y)
    248: {
    249:   PC_SVD *jac  = (PC_SVD *)pc->data;
    250:   Vec     work = jac->work, xred, yred;

    252:   PetscFunctionBegin;
    253:   PetscCall(PCSVDGetVec(pc, PC_LEFT, READ, x, &xred));
    254:   PetscCall(PCSVDGetVec(pc, PC_RIGHT, WRITE, y, &yred));
    255:   PetscCall(MatMult(jac->Vt, xred, work));
    256:   PetscCall(VecPointwiseMult(work, work, jac->diag));
    257:   PetscCall(MatMult(jac->U, work, yred));
    258:   PetscCall(PCSVDRestoreVec(pc, PC_LEFT, READ, x, &xred));
    259:   PetscCall(PCSVDRestoreVec(pc, PC_RIGHT, WRITE, y, &yred));
    260:   PetscFunctionReturn(PETSC_SUCCESS);
    261: }

    263: static PetscErrorCode PCReset_SVD(PC pc)
    264: {
    265:   PC_SVD *jac = (PC_SVD *)pc->data;

    267:   PetscFunctionBegin;
    268:   PetscCall(MatDestroy(&jac->A));
    269:   PetscCall(MatDestroy(&jac->U));
    270:   PetscCall(MatDestroy(&jac->Vt));
    271:   PetscCall(VecDestroy(&jac->diag));
    272:   PetscCall(VecDestroy(&jac->work));
    273:   PetscCall(VecScatterDestroy(&jac->right2red));
    274:   PetscCall(VecScatterDestroy(&jac->left2red));
    275:   PetscCall(VecDestroy(&jac->rightred));
    276:   PetscCall(VecDestroy(&jac->leftred));
    277:   PetscFunctionReturn(PETSC_SUCCESS);
    278: }

    280: /*
    281:    PCDestroy_SVD - Destroys the private context for the SVD preconditioner
    282:    that was created with PCCreate_SVD().

    284:    Input Parameter:
    285: .  pc - the preconditioner context

    287:    Application Interface Routine: PCDestroy()
    288: */
    289: static PetscErrorCode PCDestroy_SVD(PC pc)
    290: {
    291:   PC_SVD *jac = (PC_SVD *)pc->data;

    293:   PetscFunctionBegin;
    294:   PetscCall(PCReset_SVD(pc));
    295:   PetscCall(PetscViewerDestroy(&jac->monitor));
    296:   PetscCall(PetscFree(pc->data));
    297:   PetscFunctionReturn(PETSC_SUCCESS);
    298: }

    300: static PetscErrorCode PCSetFromOptions_SVD(PC pc, PetscOptionItems PetscOptionsObject)
    301: {
    302:   PC_SVD   *jac = (PC_SVD *)pc->data;
    303:   PetscBool flg;

    305:   PetscFunctionBegin;
    306:   PetscOptionsHeadBegin(PetscOptionsObject, "SVD options");
    307:   PetscCall(PetscOptionsReal("-pc_svd_zero_sing", "Singular values smaller than this treated as zero", "None", jac->zerosing, &jac->zerosing, NULL));
    308:   PetscCall(PetscOptionsInt("-pc_svd_ess_rank", "Essential rank of operator (0 to use entire operator)", "None", jac->essrank, &jac->essrank, NULL));
    309:   PetscCall(PetscOptionsViewer("-pc_svd_monitor", "Monitor the conditioning, and extremal singular values", "None", &jac->monitor, &jac->monitorformat, &flg));
    310:   PetscOptionsHeadEnd();
    311:   PetscFunctionReturn(PETSC_SUCCESS);
    312: }

    314: static PetscErrorCode PCView_SVD(PC pc, PetscViewer viewer)
    315: {
    316:   PC_SVD   *svd = (PC_SVD *)pc->data;
    317:   PetscBool isascii;

    319:   PetscFunctionBegin;
    320:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    321:   if (isascii) {
    322:     PetscCall(PetscViewerASCIIPrintf(viewer, "  All singular values smaller than %g treated as zero\n", (double)svd->zerosing));
    323:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Provided essential rank of the matrix %" PetscInt_FMT " (all other eigenvalues are zeroed)\n", svd->essrank));
    324:   }
    325:   PetscFunctionReturn(PETSC_SUCCESS);
    326: }

    328: /*
    329:    PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
    330:    and sets this as the private data within the generic preconditioning
    331:    context, PC, that was created within PCCreate().

    333:    Input Parameter:
    334: .  pc - the preconditioner context

    336:    Application Interface Routine: PCCreate()
    337: */

    339: /*MC
    340:    PCSVD - Use pseudo inverse defined by SVD of operator

    342:    Level: advanced

    344:   Options Database Keys:
    345: +  -pc_svd_zero_sing <rtol> - Singular values smaller than this are treated as zero
    346: -  -pc_svd_monitor          - Print information on the extreme singular values of the operator

    348:   Developer Note:
    349:   This implementation automatically creates a redundant copy of the
    350:   matrix on each process and uses a sequential SVD solve. Why does it do this instead
    351:   of using the composable `PCREDUNDANT` object?

    353: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCREDUNDANT`
    354: M*/

    356: PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
    357: {
    358:   PC_SVD     *jac;
    359:   PetscMPIInt size = 0;

    361:   PetscFunctionBegin;
    362:   /*
    363:      Creates the private data structure for this preconditioner and
    364:      attach it to the PC object.
    365:   */
    366:   PetscCall(PetscNew(&jac));
    367:   jac->zerosing = 1.e-12;
    368:   pc->data      = (void *)jac;

    370:   /*
    371:       Set the pointers for the functions that are provided above.
    372:       Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
    373:       are called, they will automatically call these functions.  Note we
    374:       choose not to provide a couple of these functions since they are
    375:       not needed.
    376:   */

    378: #if defined(PETSC_HAVE_COMPLEX)
    379:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
    380: #endif
    381:   if (size == 1) pc->ops->matapply = PCMatApply_SVD;
    382:   pc->ops->apply           = PCApply_SVD;
    383:   pc->ops->applytranspose  = PCApplyTranspose_SVD;
    384:   pc->ops->setup           = PCSetUp_SVD;
    385:   pc->ops->reset           = PCReset_SVD;
    386:   pc->ops->destroy         = PCDestroy_SVD;
    387:   pc->ops->setfromoptions  = PCSetFromOptions_SVD;
    388:   pc->ops->view            = PCView_SVD;
    389:   pc->ops->applyrichardson = NULL;
    390:   PetscFunctionReturn(PETSC_SUCCESS);
    391: }