1
2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
3
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;
12
13
MatMult_LRC(Mat N,Vec x,Vec y)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;
21
22 PetscFunctionBegin;
23 ierr = VecGetArrayRead(x,&a);CHKERRQ(ierr);
24 ierr = VecPlaceArray(Na->xl,a);CHKERRQ(ierr);
25 ierr = VecGetLocalVector(y,Na->yl);CHKERRQ(ierr);
26 ierr = MatDenseGetLocalMatrix(Na->U,&Uloc);CHKERRQ(ierr);
27 ierr = MatDenseGetLocalMatrix(Na->V,&Vloc);CHKERRQ(ierr);
28
29 /* multiply the local part of V with the local part of x */
30 #if defined(PETSC_USE_COMPLEX)
31 ierr = MatMultHermitianTranspose(Vloc,Na->xl,Na->work1);CHKERRQ(ierr);
32 #else
33 ierr = MatMultTranspose(Vloc,Na->xl,Na->work1);CHKERRQ(ierr);
34 #endif
35
36 /* form the sum of all the local multiplies: this is work2 = V'*x =
37 sum_{all processors} work1 */
38 ierr = VecGetArray(Na->work1,&w1);CHKERRQ(ierr);
39 ierr = VecGetArray(Na->work2,&w2);CHKERRQ(ierr);
40 ierr = MPIU_Allreduce(w1,w2,Na->nwork,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)N));CHKERRQ(ierr);
41 ierr = VecRestoreArray(Na->work1,&w1);CHKERRQ(ierr);
42 ierr = VecRestoreArray(Na->work2,&w2);CHKERRQ(ierr);
43
44 if (Na->c) { /* work2 = C*work2 */
45 ierr = VecPointwiseMult(Na->work2,Na->c,Na->work2);CHKERRQ(ierr);
46 }
47
48 if (Na->A) {
49 /* form y = A*x */
50 ierr = MatMult(Na->A,x,y);CHKERRQ(ierr);
51 /* multiply-add y = y + U*work2 */
52 ierr = MatMultAdd(Uloc,Na->work2,Na->yl,Na->yl);CHKERRQ(ierr);
53 } else {
54 /* multiply y = U*work2 */
55 ierr = MatMult(Uloc,Na->work2,Na->yl);CHKERRQ(ierr);
56 }
57
58 ierr = VecRestoreArrayRead(x,&a);CHKERRQ(ierr);
59 ierr = VecResetArray(Na->xl);CHKERRQ(ierr);
60 ierr = VecRestoreLocalVector(y,Na->yl);CHKERRQ(ierr);
61 PetscFunctionReturn(0);
62 }
63
MatDestroy_LRC(Mat N)64 PetscErrorCode MatDestroy_LRC(Mat N)
65 {
66 Mat_LRC *Na = (Mat_LRC*)N->data;
67 PetscErrorCode ierr;
68
69 PetscFunctionBegin;
70 ierr = MatDestroy(&Na->A);CHKERRQ(ierr);
71 ierr = MatDestroy(&Na->U);CHKERRQ(ierr);
72 ierr = MatDestroy(&Na->V);CHKERRQ(ierr);
73 ierr = VecDestroy(&Na->c);CHKERRQ(ierr);
74 ierr = VecDestroy(&Na->work1);CHKERRQ(ierr);
75 ierr = VecDestroy(&Na->work2);CHKERRQ(ierr);
76 ierr = VecDestroy(&Na->xl);CHKERRQ(ierr);
77 ierr = VecDestroy(&Na->yl);CHKERRQ(ierr);
78 ierr = PetscFree(N->data);CHKERRQ(ierr);
79 ierr = PetscObjectComposeFunction((PetscObject)N,"MatLRCGetMats_C",0);CHKERRQ(ierr);
80 PetscFunctionReturn(0);
81 }
82
MatLRCGetMats_LRC(Mat N,Mat * A,Mat * U,Vec * c,Mat * V)83 PetscErrorCode MatLRCGetMats_LRC(Mat N,Mat *A,Mat *U,Vec *c,Mat *V)
84 {
85 Mat_LRC *Na = (Mat_LRC*)N->data;
86
87 PetscFunctionBegin;
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 PetscFunctionReturn(0);
93 }
94
95 /*@
96 MatLRCGetMats - Returns the constituents of an LRC matrix
97
98 Collective on Mat
99
100 Input Parameter:
101 . N - matrix of type LRC
102
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
107
108 Note:
109 The returned matrices need not be destroyed by the caller.
110
111 Level: intermediate
112
113 .seealso: MatCreateLRC()
114 @*/
MatLRCGetMats(Mat N,Mat * A,Mat * U,Vec * c,Mat * V)115 PetscErrorCode MatLRCGetMats(Mat N,Mat *A,Mat *U,Vec *c,Mat *V)
116 {
117 PetscErrorCode ierr;
118
119 PetscFunctionBegin;
120 ierr = PetscUseMethod(N,"MatLRCGetMats_C",(Mat,Mat*,Mat*,Vec*,Mat*),(N,A,U,c,V));CHKERRQ(ierr);
121 PetscFunctionReturn(0);
122 }
123
124 /*@
125 MatCreateLRC - Creates a new matrix object that behaves like A + U*C*V'
126
127 Collective on Mat
128
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)
133
134 Output Parameter:
135 . N - the matrix that represents A + U*C*V'
136
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.
141
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.
144
145 If A is NULL then the new object behaves like a low-rank matrix U*C*V'.
146
147 Use V=U (or V=NULL) for a symmetric low-rank correction, A + U*C*U'.
148
149 If c is NULL then the low-rank correction is just U*V'.
150
151 Level: intermediate
152
153 .seealso: MatLRCGetMats()
154 @*/
MatCreateLRC(Mat A,Mat U,Vec c,Mat V,Mat * N)155 PetscErrorCode MatCreateLRC(Mat A,Mat U,Vec c,Mat V,Mat *N)
156 {
157 PetscErrorCode ierr;
158 PetscBool match;
159 PetscInt m,n,k,m1,n1,k1;
160 Mat_LRC *Na;
161
162 PetscFunctionBegin;
163 if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,1);
164 PetscValidHeaderSpecific(U,MAT_CLASSID,2);
165 if (c) PetscValidHeaderSpecific(c,VEC_CLASSID,3);
166 if (V) PetscValidHeaderSpecific(V,MAT_CLASSID,4);
167 else V=U;
168 if (A) PetscCheckSameComm(A,1,U,2);
169 PetscCheckSameComm(U,2,V,4);
170
171 ierr = PetscObjectTypeCompareAny((PetscObject)U,&match,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr);
172 if (!match) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"Matrix U must be of type dense");
173 if (V) {
174 ierr = PetscObjectTypeCompareAny((PetscObject)V,&match,MATSEQDENSE,MATMPIDENSE,"");CHKERRQ(ierr);
175 if (!match) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Matrix V must be of type dense");
176 }
177
178 ierr = MatGetSize(U,NULL,&k);CHKERRQ(ierr);
179 ierr = MatGetSize(V,NULL,&k1);CHKERRQ(ierr);
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 ierr = MatGetLocalSize(U,&m,NULL);CHKERRQ(ierr);
182 ierr = MatGetLocalSize(V,&n,NULL);CHKERRQ(ierr);
183 if (A) {
184 ierr = MatGetLocalSize(A,&m1,&n1);CHKERRQ(ierr);
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 ierr = VecGetSize(c,&k1);CHKERRQ(ierr);
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 ierr = VecGetLocalSize(c,&k1);CHKERRQ(ierr);
192 if (k!=k1) SETERRQ(PetscObjectComm((PetscObject)c),PETSC_ERR_ARG_INCOMP,"c must be a sequential vector");
193 }
194
195 ierr = MatCreate(PetscObjectComm((PetscObject)U),N);CHKERRQ(ierr);
196 ierr = MatSetSizes(*N,m,n,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr);
197 ierr = PetscObjectChangeTypeName((PetscObject)*N,MATLRC);CHKERRQ(ierr);
198
199 ierr = PetscNewLog(*N,&Na);CHKERRQ(ierr);
200 (*N)->data = (void*)Na;
201 Na->A = A;
202 Na->U = U;
203 Na->c = c;
204 Na->V = V;
205
206 if (A) { ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); }
207 ierr = PetscObjectReference((PetscObject)Na->U);CHKERRQ(ierr);
208 ierr = PetscObjectReference((PetscObject)Na->V);CHKERRQ(ierr);
209 if (c) { ierr = PetscObjectReference((PetscObject)c);CHKERRQ(ierr); }
210
211 ierr = VecCreateSeq(PETSC_COMM_SELF,U->cmap->N,&Na->work1);CHKERRQ(ierr);
212 ierr = VecDuplicate(Na->work1,&Na->work2);CHKERRQ(ierr);
213 ierr = PetscMPIIntCast(U->cmap->N,&Na->nwork);CHKERRQ(ierr);
214
215 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,V->rmap->n,NULL,&Na->xl);CHKERRQ(ierr);
216 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,U->rmap->n,NULL,&Na->yl);CHKERRQ(ierr);
217
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;
226
227 ierr = PetscObjectComposeFunction((PetscObject)(*N),"MatLRCGetMats_C",MatLRCGetMats_LRC);CHKERRQ(ierr);
228 PetscFunctionReturn(0);
229 }
230