Actual source code: lrc.c


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

  4: PETSC_EXTERN PetscErrorCode VecGetRootType_Private(Vec, VecType *);

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

 14: static PetscErrorCode MatMult_LRC_kernel(Mat N, Vec x, Vec y, PetscBool transpose)
 15: {
 16:   Mat_LRC    *Na = (Mat_LRC *)N->data;
 17:   PetscMPIInt size;
 18:   Mat         U, V;

 20:   U = transpose ? Na->V : Na->U;
 21:   V = transpose ? Na->U : Na->V;
 22:   MPI_Comm_size(PetscObjectComm((PetscObject)N), &size);
 23:   if (size == 1) {
 24:     MatMultHermitianTranspose(V, x, Na->work1);
 25:     if (Na->c) VecPointwiseMult(Na->work1, Na->c, Na->work1);
 26:     if (Na->A) {
 27:       if (transpose) {
 28:         MatMultTranspose(Na->A, x, y);
 29:       } else {
 30:         MatMult(Na->A, x, y);
 31:       }
 32:       MatMultAdd(U, Na->work1, y, y);
 33:     } else {
 34:       MatMult(U, Na->work1, y);
 35:     }
 36:   } else {
 37:     Mat                Uloc, Vloc;
 38:     Vec                yl, xl;
 39:     const PetscScalar *w1;
 40:     PetscScalar       *w2;
 41:     PetscInt           nwork;
 42:     PetscMPIInt        mpinwork;

 44:     xl = transpose ? Na->yl : Na->xl;
 45:     yl = transpose ? Na->xl : Na->yl;
 46:     VecGetLocalVector(y, yl);
 47:     MatDenseGetLocalMatrix(U, &Uloc);
 48:     MatDenseGetLocalMatrix(V, &Vloc);

 50:     /* multiply the local part of V with the local part of x */
 51:     VecGetLocalVectorRead(x, xl);
 52:     MatMultHermitianTranspose(Vloc, xl, Na->work1);
 53:     VecRestoreLocalVectorRead(x, xl);

 55:     /* form the sum of all the local multiplies: this is work2 = V'*x =
 56:        sum_{all processors} work1 */
 57:     VecGetArrayRead(Na->work1, &w1);
 58:     VecGetArrayWrite(Na->work2, &w2);
 59:     VecGetLocalSize(Na->work1, &nwork);
 60:     PetscMPIIntCast(nwork, &mpinwork);
 61:     MPIU_Allreduce(w1, w2, mpinwork, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N));
 62:     VecRestoreArrayRead(Na->work1, &w1);
 63:     VecRestoreArrayWrite(Na->work2, &w2);

 65:     if (Na->c) { /* work2 = C*work2 */
 66:       VecPointwiseMult(Na->work2, Na->c, Na->work2);
 67:     }

 69:     if (Na->A) {
 70:       /* form y = A*x or A^t*x */
 71:       if (transpose) {
 72:         MatMultTranspose(Na->A, x, y);
 73:       } else {
 74:         MatMult(Na->A, x, y);
 75:       }
 76:       /* multiply-add y = y + U*work2 */
 77:       MatMultAdd(Uloc, Na->work2, yl, yl);
 78:     } else {
 79:       /* multiply y = U*work2 */
 80:       MatMult(Uloc, Na->work2, yl);
 81:     }

 83:     VecRestoreLocalVector(y, yl);
 84:   }
 85:   return 0;
 86: }

 88: static PetscErrorCode MatMult_LRC(Mat N, Vec x, Vec y)
 89: {
 90:   MatMult_LRC_kernel(N, x, y, PETSC_FALSE);
 91:   return 0;
 92: }

 94: static PetscErrorCode MatMultTranspose_LRC(Mat N, Vec x, Vec y)
 95: {
 96:   MatMult_LRC_kernel(N, x, y, PETSC_TRUE);
 97:   return 0;
 98: }

100: static PetscErrorCode MatDestroy_LRC(Mat N)
101: {
102:   Mat_LRC *Na = (Mat_LRC *)N->data;

104:   MatDestroy(&Na->A);
105:   MatDestroy(&Na->U);
106:   MatDestroy(&Na->V);
107:   VecDestroy(&Na->c);
108:   VecDestroy(&Na->work1);
109:   VecDestroy(&Na->work2);
110:   VecDestroy(&Na->xl);
111:   VecDestroy(&Na->yl);
112:   PetscFree(N->data);
113:   PetscObjectComposeFunction((PetscObject)N, "MatLRCGetMats_C", NULL);
114:   return 0;
115: }

117: static PetscErrorCode MatLRCGetMats_LRC(Mat N, Mat *A, Mat *U, Vec *c, Mat *V)
118: {
119:   Mat_LRC *Na = (Mat_LRC *)N->data;

121:   if (A) *A = Na->A;
122:   if (U) *U = Na->U;
123:   if (c) *c = Na->c;
124:   if (V) *V = Na->V;
125:   return 0;
126: }

128: /*@
129:    MatLRCGetMats - Returns the constituents of an LRC matrix

131:    Collective

133:    Input Parameter:
134: .  N - matrix of type `MATLRC`

136:    Output Parameters:
137: +  A - the (sparse) matrix
138: .  U - first dense rectangular (tall and skinny) matrix
139: .  c - a sequential vector containing the diagonal of C
140: -  V - second dense rectangular (tall and skinny) matrix

142:    Note:
143:    The returned matrices need not be destroyed by the caller.

145:    Level: intermediate

147: .seealso: `MATLRC`, `MatCreateLRC()`
148: @*/
149: PetscErrorCode MatLRCGetMats(Mat N, Mat *A, Mat *U, Vec *c, Mat *V)
150: {
151:   PetscUseMethod(N, "MatLRCGetMats_C", (Mat, Mat *, Mat *, Vec *, Mat *), (N, A, U, c, V));
152:   return 0;
153: }

155: /*MC
156:   MATLRC -  "lrc" - a matrix object that behaves like A + U*C*V'

158:   Note:
159:    The matrix A + U*C*V' is not formed! Rather the matrix  object performs the matrix-vector product `MatMult()`, by first multiplying by
160:    A and then adding the other term.

162:   Level: advanced

164: .seealso: `MatCreateLRC`
165: M*/

167: /*@
168:    MatCreateLRC - Creates a new matrix object that behaves like A + U*C*V' of type `MATLRC`

170:    Collective

172:    Input Parameters:
173: +  A    - the (sparse) matrix (can be NULL)
174: .  U, V - two dense rectangular (tall and skinny) matrices
175: -  c    - a vector containing the diagonal of C (can be NULL)

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

180:    Notes:
181:    The matrix A + U*C*V' is not formed! Rather the new matrix
182:    object performs the matrix-vector product `MatMult()`, by first multiplying by
183:    A and then adding the other term.

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

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

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

192:    If c is NULL then the low-rank correction is just U*V'.
193:    If a sequential c vector is used for a parallel matrix,
194:    PETSc assumes that the values of the vector are consistently set across processors.

196:    Level: intermediate

198: .seealso: `MATLRC`, `MatLRCGetMats()`
199: @*/
200: PetscErrorCode MatCreateLRC(Mat A, Mat U, Vec c, Mat V, Mat *N)
201: {
202:   PetscBool   match;
203:   PetscInt    m, n, k, m1, n1, k1;
204:   Mat_LRC    *Na;
205:   Mat         Uloc;
206:   PetscMPIInt size, csize = 0;

211:   if (V) {
214:   }

217:   if (!V) V = U;
218:   PetscObjectBaseTypeCompareAny((PetscObject)U, &match, MATSEQDENSE, MATMPIDENSE, "");
220:   PetscObjectBaseTypeCompareAny((PetscObject)V, &match, MATSEQDENSE, MATMPIDENSE, "");
222:   PetscStrcmp(U->defaultvectype, V->defaultvectype, &match);
224:   if (A) {
225:     PetscStrcmp(A->defaultvectype, U->defaultvectype, &match);
227:   }

229:   MPI_Comm_size(PetscObjectComm((PetscObject)U), &size);
230:   MatGetSize(U, NULL, &k);
231:   MatGetSize(V, NULL, &k1);
233:   MatGetLocalSize(U, &m, NULL);
234:   MatGetLocalSize(V, &n, NULL);
235:   if (A) {
236:     MatGetLocalSize(A, &m1, &n1);
239:   }
240:   if (c) {
241:     MPI_Comm_size(PetscObjectComm((PetscObject)c), &csize);
242:     VecGetSize(c, &k1);
245:   }

247:   MatCreate(PetscObjectComm((PetscObject)U), N);
248:   MatSetSizes(*N, m, n, PETSC_DECIDE, PETSC_DECIDE);
249:   MatSetVecType(*N, U->defaultvectype);
250:   PetscObjectChangeTypeName((PetscObject)*N, MATLRC);
251:   /* Flag matrix as symmetric if A is symmetric and U == V */
252:   MatSetOption(*N, MAT_SYMMETRIC, (PetscBool)((A ? A->symmetric == PETSC_BOOL3_TRUE : PETSC_TRUE) && U == V));

254:   PetscNew(&Na);
255:   (*N)->data = (void *)Na;
256:   Na->A      = A;
257:   Na->U      = U;
258:   Na->c      = c;
259:   Na->V      = V;

261:   PetscObjectReference((PetscObject)A);
262:   PetscObjectReference((PetscObject)Na->U);
263:   PetscObjectReference((PetscObject)Na->V);
264:   PetscObjectReference((PetscObject)c);

266:   MatDenseGetLocalMatrix(Na->U, &Uloc);
267:   MatCreateVecs(Uloc, &Na->work1, NULL);
268:   if (size != 1) {
269:     Mat Vloc;

271:     if (Na->c && csize != 1) { /* scatter parallel vector to sequential */
272:       VecScatter sct;

274:       VecScatterCreateToAll(Na->c, &sct, &c);
275:       VecScatterBegin(sct, Na->c, c, INSERT_VALUES, SCATTER_FORWARD);
276:       VecScatterEnd(sct, Na->c, c, INSERT_VALUES, SCATTER_FORWARD);
277:       VecScatterDestroy(&sct);
278:       VecDestroy(&Na->c);
279:       Na->c = c;
280:     }
281:     MatDenseGetLocalMatrix(Na->V, &Vloc);
282:     VecDuplicate(Na->work1, &Na->work2);
283:     MatCreateVecs(Vloc, NULL, &Na->xl);
284:     MatCreateVecs(Uloc, NULL, &Na->yl);
285:   }

287:   /* Internally create a scaling vector if roottypes do not match */
288:   if (Na->c) {
289:     VecType rt1, rt2;

291:     VecGetRootType_Private(Na->work1, &rt1);
292:     VecGetRootType_Private(Na->c, &rt2);
293:     PetscStrcmp(rt1, rt2, &match);
294:     if (!match) {
295:       VecDuplicate(Na->c, &c);
296:       VecCopy(Na->c, c);
297:       VecDestroy(&Na->c);
298:       Na->c = c;
299:     }
300:   }

302:   (*N)->ops->destroy       = MatDestroy_LRC;
303:   (*N)->ops->mult          = MatMult_LRC;
304:   (*N)->ops->multtranspose = MatMultTranspose_LRC;

306:   (*N)->assembled    = PETSC_TRUE;
307:   (*N)->preallocated = PETSC_TRUE;

309:   PetscObjectComposeFunction((PetscObject)(*N), "MatLRCGetMats_C", MatLRCGetMats_LRC);
310:   MatSetUp(*N);
311:   return 0;
312: }