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: }