Actual source code: lrc.c

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

      4: typedef struct {
      5:   Mat         A;           /* sparse matrix */
      6:   Mat         U,V;         /* dense tall-skinny matrices */
      7:   Vec         c;           /* sequential vector containing the diagonal of C */
      8:   Vec         work1,work2; /* sequential vectors that hold partial products */
      9:   PetscMPIInt nwork;       /* length of work vectors */
     10:   Vec         xl,yl;       /* auxiliary sequential vectors for matmult operation */
     11: } Mat_LRC;

     14: PetscErrorCode MatMult_LRC(Mat N,Vec x,Vec y)
     15: {
     16:   Mat_LRC           *Na = (Mat_LRC*)N->data;
     17:   Mat               Uloc,Vloc;
     18:   PetscErrorCode    ierr;
     19:   PetscScalar       *w1,*w2;
     20:   const PetscScalar *a;

     23:   VecGetArrayRead(x,&a);
     24:   VecPlaceArray(Na->xl,a);
     25:   VecGetLocalVector(y,Na->yl);
     26:   MatDenseGetLocalMatrix(Na->U,&Uloc);
     27:   MatDenseGetLocalMatrix(Na->V,&Vloc);

     29:   /* multiply the local part of V with the local part of x */
     30: #if defined(PETSC_USE_COMPLEX)
     31:   MatMultHermitianTranspose(Vloc,Na->xl,Na->work1);
     32: #else
     33:   MatMultTranspose(Vloc,Na->xl,Na->work1);
     34: #endif

     36:   /* form the sum of all the local multiplies: this is work2 = V'*x =
     37:      sum_{all processors} work1 */
     38:   VecGetArray(Na->work1,&w1);
     39:   VecGetArray(Na->work2,&w2);
     40:   MPIU_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));
     41:   VecRestoreArray(Na->work1,&w1);
     42:   VecRestoreArray(Na->work2,&w2);

     44:   if (Na->c) {  /* work2 = C*work2 */
     45:     VecPointwiseMult(Na->work2,Na->c,Na->work2);
     46:   }

     48:   if (Na->A) {
     49:     /* form y = A*x */
     50:     MatMult(Na->A,x,y);
     51:     /* multiply-add y = y + U*work2 */
     52:     MatMultAdd(Uloc,Na->work2,Na->yl,Na->yl);
     53:   } else {
     54:     /* multiply y = U*work2 */
     55:     MatMult(Uloc,Na->work2,Na->yl);
     56:   }

     58:   VecRestoreArrayRead(x,&a);
     59:   VecResetArray(Na->xl);
     60:   VecRestoreLocalVector(y,Na->yl);
     61:   return(0);
     62: }

     64: PetscErrorCode MatDestroy_LRC(Mat N)
     65: {
     66:   Mat_LRC        *Na = (Mat_LRC*)N->data;

     70:   MatDestroy(&Na->A);
     71:   MatDestroy(&Na->U);
     72:   MatDestroy(&Na->V);
     73:   VecDestroy(&Na->c);
     74:   VecDestroy(&Na->work1);
     75:   VecDestroy(&Na->work2);
     76:   VecDestroy(&Na->xl);
     77:   VecDestroy(&Na->yl);
     78:   PetscFree(N->data);
     79:   PetscObjectComposeFunction((PetscObject)N,"MatLRCGetMats_C",0);
     80:   return(0);
     81: }

     83: PetscErrorCode MatLRCGetMats_LRC(Mat N,Mat *A,Mat *U,Vec *c,Mat *V)
     84: {
     85:   Mat_LRC *Na = (Mat_LRC*)N->data;

     88:   if (A) *A = Na->A;
     89:   if (U) *U = Na->U;
     90:   if (c) *c = Na->c;
     91:   if (V) *V = Na->V;
     92:   return(0);
     93: }

     95: /*@
     96:    MatLRCGetMats - Returns the constituents of an LRC matrix

     98:    Collective on Mat

    100:    Input Parameter:
    101: .  N    - matrix of type LRC

    103:    Output Parameters:
    104: +  A    - the (sparse) matrix
    105: .  U, V - two dense rectangular (tall and skinny) matrices
    106: -  c    - a sequential vector containing the diagonal of C

    108:    Note:
    109:    The returned matrices need not be destroyed by the caller.

    111:    Level: intermediate

    113: .seealso: MatCreateLRC()
    114: @*/
    115: PetscErrorCode  MatLRCGetMats(Mat N,Mat *A,Mat *U,Vec *c,Mat *V)
    116: {

    120:   PetscUseMethod(N,"MatLRCGetMats_C",(Mat,Mat*,Mat*,Vec*,Mat*),(N,A,U,c,V));
    121:   return(0);
    122: }

    124: /*@
    125:    MatCreateLRC - Creates a new matrix object that behaves like A + U*C*V'

    127:    Collective on Mat

    129:    Input Parameters:
    130: +  A    - the (sparse) matrix (can be NULL)
    131: .  U, V - two dense rectangular (tall and skinny) matrices
    132: -  c    - a sequential vector containing the diagonal of C (can be NULL)

    134:    Output Parameter:
    135: .  N    - the matrix that represents A + U*C*V'

    137:    Notes:
    138:    The matrix A + U*C*V' is not formed! Rather the new matrix
    139:    object performs the matrix-vector product by first multiplying by
    140:    A and then adding the other term.

    142:    C is a diagonal matrix (represented as a vector) of order k,
    143:    where k is the number of columns of both U and V.

    145:    If A is NULL then the new object behaves like a low-rank matrix U*C*V'.

    147:    Use V=U (or V=NULL) for a symmetric low-rank correction, A + U*C*U'.

    149:    If c is NULL then the low-rank correction is just U*V'.

    151:    Level: intermediate

    153: .seealso: MatLRCGetMats()
    154: @*/
    155: PetscErrorCode MatCreateLRC(Mat A,Mat U,Vec c,Mat V,Mat *N)
    156: {
    158:   PetscBool      match;
    159:   PetscInt       m,n,k,m1,n1,k1;
    160:   Mat_LRC        *Na;

    167:   else V=U;

    171:   PetscObjectTypeCompareAny((PetscObject)U,&match,MATSEQDENSE,MATMPIDENSE,"");
    172:   if (!match) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"Matrix U must be of type dense");
    173:   if (V) {
    174:     PetscObjectTypeCompareAny((PetscObject)V,&match,MATSEQDENSE,MATMPIDENSE,"");
    175:     if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Matrix V must be of type dense");
    176:   }

    178:   MatGetSize(U,NULL,&k);
    179:   MatGetSize(V,NULL,&k1);
    180:   if (k!=k1) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"U and V have different number of columns (%D vs %D)",k,k1);
    181:   MatGetLocalSize(U,&m,NULL);
    182:   MatGetLocalSize(V,&n,NULL);
    183:   if (A) {
    184:     MatGetLocalSize(A,&m1,&n1);
    185:     if (m!=m1) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_INCOMP,"Local dimensions of U %D and A %D do not match",m,m1);
    186:     if (n!=n1) SETERRQ2(PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Local dimensions of V %D and A %D do not match",n,n1);
    187:   }
    188:   if (c) {
    189:     VecGetSize(c,&k1);
    190:     if (k!=k1) SETERRQ2(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"The length of c %D does not match the number of columns of U and V (%D)",k1,k);
    191:     VecGetLocalSize(c,&k1);
    192:     if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"c must be a sequential vector");
    193:   }

    195:   MatCreate(PetscObjectComm((PetscObject)U),N);
    196:   MatSetSizes(*N,m,n,PETSC_DECIDE,PETSC_DECIDE);
    197:   PetscObjectChangeTypeName((PetscObject)*N,MATLRC);

    199:   PetscNewLog(*N,&Na);
    200:   (*N)->data = (void*)Na;
    201:   Na->A      = A;
    202:   Na->U      = U;
    203:   Na->c      = c;
    204:   Na->V      = V;

    206:   if (A) { PetscObjectReference((PetscObject)A); }
    207:   PetscObjectReference((PetscObject)Na->U);
    208:   PetscObjectReference((PetscObject)Na->V);
    209:   if (c) { PetscObjectReference((PetscObject)c); }

    211:   VecCreateSeq(PETSC_COMM_SELF,U->cmap->N,&Na->work1);
    212:   VecDuplicate(Na->work1,&Na->work2);
    213:   PetscMPIIntCast(U->cmap->N,&Na->nwork);

    215:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,V->rmap->n,NULL,&Na->xl);
    216:   VecCreateSeqWithArray(PETSC_COMM_SELF,1,U->rmap->n,NULL,&Na->yl);

    218:   (*N)->ops->destroy = MatDestroy_LRC;
    219:   (*N)->ops->mult    = MatMult_LRC;
    220:   (*N)->assembled    = PETSC_TRUE;
    221:   (*N)->preallocated = PETSC_TRUE;
    222:   (*N)->cmap->N      = V->rmap->N;
    223:   (*N)->rmap->N      = U->rmap->N;
    224:   (*N)->cmap->n      = V->rmap->n;
    225:   (*N)->rmap->n      = U->rmap->n;

    227:   PetscObjectComposeFunction((PetscObject)(*N),"MatLRCGetMats_C",MatLRCGetMats_LRC);
    228:   return(0);
    229: }