Actual source code: lrc.c
petsc-3.9.4 2018-09-11
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: }