Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 : /*
11 : BV operations involving global communication
12 : */
13 :
14 : #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
15 :
16 : /*
17 : BVDot for the particular case of non-standard inner product with
18 : matrix B, which is assumed to be symmetric (or complex Hermitian)
19 : */
20 68 : static inline PetscErrorCode BVDot_Private(BV X,BV Y,Mat M)
21 : {
22 68 : PetscObjectId idx,idy;
23 68 : PetscInt i,j,jend,m;
24 68 : PetscScalar *marray;
25 68 : PetscBool symm=PETSC_FALSE;
26 68 : Mat B;
27 68 : Vec z;
28 :
29 68 : PetscFunctionBegin;
30 68 : BVCheckOp(Y,1,dotvec);
31 68 : PetscCall(MatGetSize(M,&m,NULL));
32 68 : PetscCall(MatDenseGetArray(M,&marray));
33 68 : PetscCall(PetscObjectGetId((PetscObject)X,&idx));
34 68 : PetscCall(PetscObjectGetId((PetscObject)Y,&idy));
35 68 : B = Y->matrix;
36 68 : Y->matrix = NULL;
37 68 : if (idx==idy) symm=PETSC_TRUE; /* M=X'BX is symmetric */
38 68 : jend = X->k;
39 412 : for (j=X->l;j<jend;j++) {
40 344 : if (symm) Y->k = j+1;
41 344 : PetscCall(BVGetColumn(X->cached,j,&z));
42 344 : PetscUseTypeMethod(Y,dotvec,z,marray+j*m+Y->l);
43 344 : PetscCall(BVRestoreColumn(X->cached,j,&z));
44 344 : if (symm) {
45 1032 : for (i=X->l;i<j;i++)
46 736 : marray[j+i*m] = PetscConj(marray[i+j*m]);
47 : }
48 : }
49 68 : PetscCall(MatDenseRestoreArray(M,&marray));
50 68 : Y->matrix = B;
51 68 : PetscFunctionReturn(PETSC_SUCCESS);
52 : }
53 :
54 : /*@
55 : BVDot - Computes the 'block-dot' product of two basis vectors objects.
56 :
57 : Collective
58 :
59 : Input Parameters:
60 : + X - first basis vectors
61 : - Y - second basis vectors
62 :
63 : Output Parameter:
64 : . M - the resulting matrix
65 :
66 : Notes:
67 : This is the generalization of VecDot() for a collection of vectors, M = Y^H*X.
68 : The result is a matrix M whose entry m_ij is equal to y_i^H x_j (where y_i^H
69 : denotes the conjugate transpose of y_i).
70 :
71 : If a non-standard inner product has been specified with BVSetMatrix(),
72 : then the result is M = Y^H*B*X. In this case, both X and Y must have
73 : the same associated matrix.
74 :
75 : On entry, M must be a sequential dense Mat with dimensions m,n at least, where
76 : m is the number of active columns of Y and n is the number of active columns of X.
77 : Only rows (resp. columns) of M starting from ly (resp. lx) are computed,
78 : where ly (resp. lx) is the number of leading columns of Y (resp. X).
79 :
80 : X and Y need not be different objects.
81 :
82 : Level: intermediate
83 :
84 : .seealso: BVDotVec(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
85 : @*/
86 26943 : PetscErrorCode BVDot(BV X,BV Y,Mat M)
87 : {
88 26943 : PetscInt m,n;
89 :
90 26943 : PetscFunctionBegin;
91 26943 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
92 26943 : PetscValidHeaderSpecific(Y,BV_CLASSID,2);
93 26943 : PetscValidHeaderSpecific(M,MAT_CLASSID,3);
94 26943 : PetscValidType(X,1);
95 26943 : BVCheckSizes(X,1);
96 26943 : PetscValidType(Y,2);
97 26943 : BVCheckSizes(Y,2);
98 26943 : PetscValidType(M,3);
99 26943 : PetscCheckSameTypeAndComm(X,1,Y,2);
100 26943 : PetscCheckTypeNames(M,MATSEQDENSE,MATSEQDENSECUDA);
101 :
102 26943 : PetscCall(MatGetSize(M,&m,&n));
103 26943 : PetscCheck(m>=Y->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,Y->k);
104 26943 : PetscCheck(n>=X->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,X->k);
105 26943 : PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);
106 26943 : PetscCheck(X->matrix==Y->matrix,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"X and Y must have the same inner product matrix");
107 26943 : if (X->l==X->k || Y->l==Y->k) PetscFunctionReturn(PETSC_SUCCESS);
108 :
109 26845 : PetscCall(PetscLogEventBegin(BV_Dot,X,Y,0,0));
110 26845 : if (X->matrix) { /* non-standard inner product */
111 : /* compute BX first */
112 302 : PetscCall(BV_IPMatMultBV(X));
113 302 : if (X->vmm==BV_MATMULT_VECS) {
114 : /* perform computation column by column */
115 68 : PetscCall(BVDot_Private(X,Y,M));
116 234 : } else PetscUseTypeMethod(X->cached,dot,Y,M);
117 26543 : } else PetscUseTypeMethod(X,dot,Y,M);
118 26845 : PetscCall(PetscLogEventEnd(BV_Dot,X,Y,0,0));
119 26845 : PetscFunctionReturn(PETSC_SUCCESS);
120 : }
121 :
122 : /*@
123 : BVDotVec - Computes multiple dot products of a vector against all the
124 : column vectors of a BV.
125 :
126 : Collective
127 :
128 : Input Parameters:
129 : + X - basis vectors
130 : - y - a vector
131 :
132 : Output Parameter:
133 : . m - an array where the result must be placed
134 :
135 : Notes:
136 : This is analogue to VecMDot(), but using BV to represent a collection
137 : of vectors. The result is m = X^H*y, so m_i is equal to x_j^H y. Note
138 : that here X is transposed as opposed to BVDot().
139 :
140 : If a non-standard inner product has been specified with BVSetMatrix(),
141 : then the result is m = X^H*B*y.
142 :
143 : The length of array m must be equal to the number of active columns of X
144 : minus the number of leading columns, i.e. the first entry of m is the
145 : product of the first non-leading column with y.
146 :
147 : Level: intermediate
148 :
149 : .seealso: BVDot(), BVDotColumn(), BVSetActiveColumns(), BVSetMatrix()
150 : @*/
151 117180 : PetscErrorCode BVDotVec(BV X,Vec y,PetscScalar m[])
152 : {
153 117180 : PetscInt n;
154 :
155 117180 : PetscFunctionBegin;
156 117180 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
157 117180 : PetscValidHeaderSpecific(y,VEC_CLASSID,2);
158 117180 : PetscValidType(X,1);
159 117180 : BVCheckSizes(X,1);
160 117180 : BVCheckOp(X,1,dotvec);
161 117180 : PetscValidType(y,2);
162 117180 : PetscCheckSameTypeAndComm(X,1,y,2);
163 :
164 117180 : PetscCall(VecGetLocalSize(y,&n));
165 117180 : PetscCheck(X->n==n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", y %" PetscInt_FMT,X->n,n);
166 :
167 117180 : PetscCall(PetscLogEventBegin(BV_DotVec,X,y,0,0));
168 117180 : PetscUseTypeMethod(X,dotvec,y,m);
169 117180 : PetscCall(PetscLogEventEnd(BV_DotVec,X,y,0,0));
170 117180 : PetscFunctionReturn(PETSC_SUCCESS);
171 : }
172 :
173 : /*@
174 : BVDotVecBegin - Starts a split phase dot product computation.
175 :
176 : Input Parameters:
177 : + X - basis vectors
178 : . y - a vector
179 : - m - an array where the result will go (can be NULL)
180 :
181 : Note:
182 : Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().
183 :
184 : Level: advanced
185 :
186 : .seealso: BVDotVecEnd(), BVDotVec()
187 : @*/
188 112 : PetscErrorCode BVDotVecBegin(BV X,Vec y,PetscScalar *m)
189 : {
190 112 : PetscInt i,n,nv;
191 112 : PetscSplitReduction *sr;
192 112 : MPI_Comm comm;
193 :
194 112 : PetscFunctionBegin;
195 112 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
196 112 : PetscValidHeaderSpecific(y,VEC_CLASSID,2);
197 112 : PetscValidType(X,1);
198 112 : BVCheckSizes(X,1);
199 112 : PetscValidType(y,2);
200 112 : PetscCheckSameTypeAndComm(X,1,y,2);
201 :
202 112 : PetscCall(VecGetLocalSize(y,&n));
203 112 : PetscCheck(X->n==n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", y %" PetscInt_FMT,X->n,n);
204 :
205 112 : if (X->ops->dotvec_begin) PetscUseTypeMethod(X,dotvec_begin,y,m);
206 : else {
207 110 : BVCheckOp(X,1,dotvec_local);
208 110 : nv = X->k-X->l;
209 110 : PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
210 110 : PetscCall(PetscSplitReductionGet(comm,&sr));
211 110 : PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
212 1038 : for (i=0;i<nv;i++) {
213 928 : if (sr->numopsbegin+i >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
214 928 : sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
215 928 : sr->invecs[sr->numopsbegin+i] = (void*)X;
216 : }
217 110 : PetscCall(PetscLogEventBegin(BV_DotVec,X,y,0,0));
218 110 : PetscUseTypeMethod(X,dotvec_local,y,sr->lvalues+sr->numopsbegin);
219 110 : sr->numopsbegin += nv;
220 110 : PetscCall(PetscLogEventEnd(BV_DotVec,X,y,0,0));
221 : }
222 112 : PetscFunctionReturn(PETSC_SUCCESS);
223 : }
224 :
225 : /*@
226 : BVDotVecEnd - Ends a split phase dot product computation.
227 :
228 : Input Parameters:
229 : + X - basis vectors
230 : . y - a vector
231 : - m - an array where the result will go
232 :
233 : Note:
234 : Each call to BVDotVecBegin() should be paired with a call to BVDotVecEnd().
235 :
236 : Level: advanced
237 :
238 : .seealso: BVDotVecBegin(), BVDotVec()
239 : @*/
240 112 : PetscErrorCode BVDotVecEnd(BV X,Vec y,PetscScalar *m)
241 : {
242 112 : PetscInt i,nv;
243 112 : PetscSplitReduction *sr;
244 112 : MPI_Comm comm;
245 :
246 112 : PetscFunctionBegin;
247 112 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
248 112 : PetscValidType(X,1);
249 112 : BVCheckSizes(X,1);
250 :
251 112 : if (X->ops->dotvec_end) PetscUseTypeMethod(X,dotvec_end,y,m);
252 : else {
253 110 : nv = X->k-X->l;
254 110 : PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
255 110 : PetscCall(PetscSplitReductionGet(comm,&sr));
256 110 : PetscCall(PetscSplitReductionEnd(sr));
257 :
258 110 : PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times than VecxxxBegin()");
259 110 : PetscCheck((void*)X==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called BVxxxEnd() in a different order or with a different BV than BVxxxBegin()");
260 110 : PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_SUM,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Wrong type of reduction");
261 1038 : for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];
262 :
263 : /* Finished getting all the results so reset to no outstanding requests */
264 110 : if (sr->numopsend == sr->numopsbegin) {
265 104 : sr->state = STATE_BEGIN;
266 104 : sr->numopsend = 0;
267 104 : sr->numopsbegin = 0;
268 : }
269 : }
270 112 : PetscFunctionReturn(PETSC_SUCCESS);
271 : }
272 :
273 : /*@
274 : BVDotColumn - Computes multiple dot products of a column against all the
275 : previous columns of a BV.
276 :
277 : Collective
278 :
279 : Input Parameters:
280 : + X - basis vectors
281 : - j - the column index
282 :
283 : Output Parameter:
284 : . q - an array where the result must be placed
285 :
286 : Notes:
287 : This operation is equivalent to BVDotVec() but it uses column j of X
288 : rather than taking a Vec as an argument. The number of active columns of
289 : X is set to j before the computation, and restored afterwards.
290 : If X has leading columns specified, then these columns do not participate
291 : in the computation. Therefore, the length of array q must be equal to j
292 : minus the number of leading columns.
293 :
294 : Developer Notes:
295 : If q is NULL, then the result is written in position nc+l of the internal
296 : buffer vector, see BVGetBufferVec().
297 :
298 : Level: advanced
299 :
300 : .seealso: BVDot(), BVDotVec(), BVSetActiveColumns(), BVSetMatrix()
301 : @*/
302 44342 : PetscErrorCode BVDotColumn(BV X,PetscInt j,PetscScalar *q)
303 : {
304 44342 : PetscInt ksave;
305 44342 : Vec y;
306 :
307 44342 : PetscFunctionBegin;
308 44342 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
309 177368 : PetscValidLogicalCollectiveInt(X,j,2);
310 44342 : PetscValidType(X,1);
311 44342 : BVCheckSizes(X,1);
312 44342 : BVCheckOp(X,1,dotvec);
313 :
314 44342 : PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
315 44342 : PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
316 :
317 44342 : PetscCall(PetscLogEventBegin(BV_DotVec,X,0,0,0));
318 44342 : ksave = X->k;
319 44342 : X->k = j;
320 44342 : if (!q && !X->buffer) PetscCall(BVGetBufferVec(X,&X->buffer));
321 44342 : PetscCall(BVGetColumn(X,j,&y));
322 44342 : PetscUseTypeMethod(X,dotvec,y,q);
323 44342 : PetscCall(BVRestoreColumn(X,j,&y));
324 44342 : X->k = ksave;
325 44342 : PetscCall(PetscLogEventEnd(BV_DotVec,X,0,0,0));
326 44342 : PetscFunctionReturn(PETSC_SUCCESS);
327 : }
328 :
329 : /*@
330 : BVDotColumnBegin - Starts a split phase dot product computation.
331 :
332 : Input Parameters:
333 : + X - basis vectors
334 : . j - the column index
335 : - m - an array where the result will go (can be NULL)
336 :
337 : Note:
338 : Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().
339 :
340 : Level: advanced
341 :
342 : .seealso: BVDotColumnEnd(), BVDotColumn()
343 : @*/
344 2385 : PetscErrorCode BVDotColumnBegin(BV X,PetscInt j,PetscScalar *m)
345 : {
346 2385 : PetscInt i,nv,ksave;
347 2385 : PetscSplitReduction *sr;
348 2385 : MPI_Comm comm;
349 2385 : Vec y;
350 :
351 2385 : PetscFunctionBegin;
352 2385 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
353 9540 : PetscValidLogicalCollectiveInt(X,j,2);
354 2385 : PetscValidType(X,1);
355 2385 : BVCheckSizes(X,1);
356 :
357 2385 : PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
358 2385 : PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
359 2385 : ksave = X->k;
360 2385 : X->k = j;
361 2385 : PetscCall(BVGetColumn(X,j,&y));
362 :
363 2385 : if (X->ops->dotvec_begin) PetscUseTypeMethod(X,dotvec_begin,y,m);
364 : else {
365 2383 : BVCheckOp(X,1,dotvec_local);
366 2383 : nv = X->k-X->l;
367 2383 : PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
368 2383 : PetscCall(PetscSplitReductionGet(comm,&sr));
369 2383 : PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)X),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
370 24458 : for (i=0;i<nv;i++) {
371 22075 : if (sr->numopsbegin+i >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
372 22075 : sr->reducetype[sr->numopsbegin+i] = PETSC_SR_REDUCE_SUM;
373 22075 : sr->invecs[sr->numopsbegin+i] = (void*)X;
374 : }
375 2383 : PetscCall(PetscLogEventBegin(BV_DotVec,X,0,0,0));
376 2383 : PetscUseTypeMethod(X,dotvec_local,y,sr->lvalues+sr->numopsbegin);
377 2383 : sr->numopsbegin += nv;
378 2383 : PetscCall(PetscLogEventEnd(BV_DotVec,X,0,0,0));
379 : }
380 2385 : PetscCall(BVRestoreColumn(X,j,&y));
381 2385 : X->k = ksave;
382 2385 : PetscFunctionReturn(PETSC_SUCCESS);
383 : }
384 :
385 : /*@
386 : BVDotColumnEnd - Ends a split phase dot product computation.
387 :
388 : Input Parameters:
389 : + X - basis vectors
390 : . j - the column index
391 : - m - an array where the result will go
392 :
393 : Notes:
394 : Each call to BVDotColumnBegin() should be paired with a call to BVDotColumnEnd().
395 :
396 : Level: advanced
397 :
398 : .seealso: BVDotColumnBegin(), BVDotColumn()
399 : @*/
400 2385 : PetscErrorCode BVDotColumnEnd(BV X,PetscInt j,PetscScalar *m)
401 : {
402 2385 : PetscInt i,nv,ksave;
403 2385 : PetscSplitReduction *sr;
404 2385 : MPI_Comm comm;
405 2385 : Vec y;
406 :
407 2385 : PetscFunctionBegin;
408 2385 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
409 9540 : PetscValidLogicalCollectiveInt(X,j,2);
410 2385 : PetscValidType(X,1);
411 2385 : BVCheckSizes(X,1);
412 :
413 2385 : PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
414 2385 : PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
415 2385 : ksave = X->k;
416 2385 : X->k = j;
417 :
418 2385 : if (X->ops->dotvec_end) {
419 2 : PetscCall(BVGetColumn(X,j,&y));
420 2 : PetscUseTypeMethod(X,dotvec_end,y,m);
421 2 : PetscCall(BVRestoreColumn(X,j,&y));
422 : } else {
423 2383 : nv = X->k-X->l;
424 2383 : PetscCall(PetscObjectGetComm((PetscObject)X,&comm));
425 2383 : PetscCall(PetscSplitReductionGet(comm,&sr));
426 2383 : PetscCall(PetscSplitReductionEnd(sr));
427 :
428 2383 : PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times than VecxxxBegin()");
429 2383 : PetscCheck((void*)X==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Called BVxxxEnd() in a different order or with a different BV than BVxxxBegin()");
430 2383 : PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_SUM,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_WRONGSTATE,"Wrong type of reduction");
431 24458 : for (i=0;i<nv;i++) m[i] = sr->gvalues[sr->numopsend++];
432 :
433 : /* Finished getting all the results so reset to no outstanding requests */
434 2383 : if (sr->numopsend == sr->numopsbegin) {
435 623 : sr->state = STATE_BEGIN;
436 623 : sr->numopsend = 0;
437 623 : sr->numopsbegin = 0;
438 : }
439 : }
440 2385 : X->k = ksave;
441 2385 : PetscFunctionReturn(PETSC_SUCCESS);
442 : }
443 :
444 59139 : static inline PetscErrorCode BVNorm_Private(BV bv,Vec z,NormType type,PetscReal *val)
445 : {
446 59139 : PetscScalar p;
447 :
448 59139 : PetscFunctionBegin;
449 59139 : PetscCall(BV_IPMatMult(bv,z));
450 59139 : PetscCall(VecDot(bv->Bx,z,&p));
451 59139 : PetscCall(BV_SafeSqrt(bv,p,val));
452 59139 : PetscFunctionReturn(PETSC_SUCCESS);
453 : }
454 :
455 415 : static inline PetscErrorCode BVNorm_Begin_Private(BV bv,Vec z,NormType type,PetscReal *val)
456 : {
457 415 : PetscScalar p;
458 :
459 415 : PetscFunctionBegin;
460 415 : PetscCall(BV_IPMatMult(bv,z));
461 415 : PetscCall(VecDotBegin(bv->Bx,z,&p));
462 415 : PetscFunctionReturn(PETSC_SUCCESS);
463 : }
464 :
465 415 : static inline PetscErrorCode BVNorm_End_Private(BV bv,Vec z,NormType type,PetscReal *val)
466 : {
467 415 : PetscScalar p;
468 :
469 415 : PetscFunctionBegin;
470 415 : PetscCall(VecDotEnd(bv->Bx,z,&p));
471 415 : PetscCall(BV_SafeSqrt(bv,p,val));
472 415 : PetscFunctionReturn(PETSC_SUCCESS);
473 : }
474 :
475 : /*@
476 : BVNorm - Computes the matrix norm of the BV.
477 :
478 : Collective
479 :
480 : Input Parameters:
481 : + bv - basis vectors
482 : - type - the norm type
483 :
484 : Output Parameter:
485 : . val - the norm
486 :
487 : Notes:
488 : All active columns (except the leading ones) are considered as a matrix.
489 : The allowed norms are NORM_1, NORM_FROBENIUS, and NORM_INFINITY.
490 :
491 : This operation fails if a non-standard inner product has been
492 : specified with BVSetMatrix().
493 :
494 : Level: intermediate
495 :
496 : .seealso: BVNormVec(), BVNormColumn(), BVNormalize(), BVSetActiveColumns(), BVSetMatrix()
497 : @*/
498 291 : PetscErrorCode BVNorm(BV bv,NormType type,PetscReal *val)
499 : {
500 291 : PetscFunctionBegin;
501 291 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
502 1164 : PetscValidLogicalCollectiveEnum(bv,type,2);
503 291 : PetscAssertPointer(val,3);
504 291 : PetscValidType(bv,1);
505 291 : BVCheckSizes(bv,1);
506 :
507 291 : PetscCheck(type!=NORM_2 && type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
508 291 : PetscCheck(!bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Matrix norm not available for non-standard inner product");
509 :
510 291 : PetscCall(PetscLogEventBegin(BV_Norm,bv,0,0,0));
511 291 : PetscUseTypeMethod(bv,norm,-1,type,val);
512 291 : PetscCall(PetscLogEventEnd(BV_Norm,bv,0,0,0));
513 291 : PetscFunctionReturn(PETSC_SUCCESS);
514 : }
515 :
516 : /*@
517 : BVNormVec - Computes the norm of a given vector.
518 :
519 : Collective
520 :
521 : Input Parameters:
522 : + bv - basis vectors
523 : . v - the vector
524 : - type - the norm type
525 :
526 : Output Parameter:
527 : . val - the norm
528 :
529 : Notes:
530 : This is the analogue of BVNormColumn() but for a vector that is not in the BV.
531 : If a non-standard inner product has been specified with BVSetMatrix(),
532 : then the returned value is sqrt(v'*B*v), where B is the inner product
533 : matrix (argument 'type' is ignored). Otherwise, VecNorm() is called.
534 :
535 : Level: developer
536 :
537 : .seealso: BVNorm(), BVNormColumn(), BVSetMatrix()
538 : @*/
539 79697 : PetscErrorCode BVNormVec(BV bv,Vec v,NormType type,PetscReal *val)
540 : {
541 79697 : PetscInt n;
542 :
543 79697 : PetscFunctionBegin;
544 79697 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
545 79697 : PetscValidHeaderSpecific(v,VEC_CLASSID,2);
546 318788 : PetscValidLogicalCollectiveEnum(bv,type,3);
547 79697 : PetscAssertPointer(val,4);
548 79697 : PetscValidType(bv,1);
549 79697 : BVCheckSizes(bv,1);
550 79697 : PetscValidType(v,2);
551 79697 : PetscCheckSameComm(bv,1,v,2);
552 :
553 79697 : PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
554 :
555 79697 : PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
556 79697 : if (bv->matrix) { /* non-standard inner product */
557 6301 : PetscCall(VecGetLocalSize(v,&n));
558 6301 : PetscCheck(bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %" PetscInt_FMT ", v %" PetscInt_FMT,bv->n,n);
559 6301 : PetscCall(BVNorm_Private(bv,v,type,val));
560 73396 : } else PetscCall(VecNorm(v,type,val));
561 79697 : PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
562 79697 : PetscFunctionReturn(PETSC_SUCCESS);
563 : }
564 :
565 : /*@
566 : BVNormVecBegin - Starts a split phase norm computation.
567 :
568 : Input Parameters:
569 : + bv - basis vectors
570 : . v - the vector
571 : . type - the norm type
572 : - val - the norm
573 :
574 : Note:
575 : Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().
576 :
577 : Level: advanced
578 :
579 : .seealso: BVNormVecEnd(), BVNormVec()
580 : @*/
581 719 : PetscErrorCode BVNormVecBegin(BV bv,Vec v,NormType type,PetscReal *val)
582 : {
583 719 : PetscInt n;
584 :
585 719 : PetscFunctionBegin;
586 719 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
587 719 : PetscValidHeaderSpecific(v,VEC_CLASSID,2);
588 2876 : PetscValidLogicalCollectiveEnum(bv,type,3);
589 719 : PetscAssertPointer(val,4);
590 719 : PetscValidType(bv,1);
591 719 : BVCheckSizes(bv,1);
592 719 : PetscValidType(v,2);
593 719 : PetscCheckSameTypeAndComm(bv,1,v,2);
594 :
595 719 : PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
596 :
597 719 : PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
598 719 : if (bv->matrix) { /* non-standard inner product */
599 21 : PetscCall(VecGetLocalSize(v,&n));
600 21 : PetscCheck(bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension bv %" PetscInt_FMT ", v %" PetscInt_FMT,bv->n,n);
601 21 : PetscCall(BVNorm_Begin_Private(bv,v,type,val));
602 698 : } else PetscCall(VecNormBegin(v,type,val));
603 719 : PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
604 719 : PetscFunctionReturn(PETSC_SUCCESS);
605 : }
606 :
607 : /*@
608 : BVNormVecEnd - Ends a split phase norm computation.
609 :
610 : Input Parameters:
611 : + bv - basis vectors
612 : . v - the vector
613 : . type - the norm type
614 : - val - the norm
615 :
616 : Note:
617 : Each call to BVNormVecBegin() should be paired with a call to BVNormVecEnd().
618 :
619 : Level: advanced
620 :
621 : .seealso: BVNormVecBegin(), BVNormVec()
622 : @*/
623 719 : PetscErrorCode BVNormVecEnd(BV bv,Vec v,NormType type,PetscReal *val)
624 : {
625 719 : PetscFunctionBegin;
626 719 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
627 2876 : PetscValidLogicalCollectiveEnum(bv,type,3);
628 719 : PetscAssertPointer(val,4);
629 719 : PetscValidType(bv,1);
630 719 : BVCheckSizes(bv,1);
631 :
632 719 : PetscCheck(type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
633 :
634 719 : if (bv->matrix) PetscCall(BVNorm_End_Private(bv,v,type,val)); /* non-standard inner product */
635 698 : else PetscCall(VecNormEnd(v,type,val));
636 719 : PetscFunctionReturn(PETSC_SUCCESS);
637 : }
638 :
639 : /*@
640 : BVNormColumn - Computes the vector norm of a selected column.
641 :
642 : Collective
643 :
644 : Input Parameters:
645 : + bv - basis vectors
646 : . j - column number to be used
647 : - type - the norm type
648 :
649 : Output Parameter:
650 : . val - the norm
651 :
652 : Notes:
653 : The norm of V[j] is computed (NORM_1, NORM_2, or NORM_INFINITY).
654 : If a non-standard inner product has been specified with BVSetMatrix(),
655 : then the returned value is sqrt(V[j]'*B*V[j]),
656 : where B is the inner product matrix (argument 'type' is ignored).
657 :
658 : Level: intermediate
659 :
660 : .seealso: BVNorm(), BVNormVec(), BVNormalize(), BVSetActiveColumns(), BVSetMatrix()
661 : @*/
662 58080 : PetscErrorCode BVNormColumn(BV bv,PetscInt j,NormType type,PetscReal *val)
663 : {
664 58080 : Vec z;
665 :
666 58080 : PetscFunctionBegin;
667 58080 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
668 232320 : PetscValidLogicalCollectiveInt(bv,j,2);
669 232320 : PetscValidLogicalCollectiveEnum(bv,type,3);
670 58080 : PetscAssertPointer(val,4);
671 58080 : PetscValidType(bv,1);
672 58080 : BVCheckSizes(bv,1);
673 :
674 58080 : PetscCheck(j>=0 && j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %" PetscInt_FMT ", the number of columns is %" PetscInt_FMT,j,bv->m);
675 58080 : PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
676 :
677 58080 : PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
678 58080 : if (bv->matrix) { /* non-standard inner product */
679 52838 : PetscCall(BVGetColumn(bv,j,&z));
680 52838 : PetscCall(BVNorm_Private(bv,z,type,val));
681 52838 : PetscCall(BVRestoreColumn(bv,j,&z));
682 5242 : } else PetscUseTypeMethod(bv,norm,j,type,val);
683 58080 : PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
684 58080 : PetscFunctionReturn(PETSC_SUCCESS);
685 : }
686 :
687 : /*@
688 : BVNormColumnBegin - Starts a split phase norm computation.
689 :
690 : Input Parameters:
691 : + bv - basis vectors
692 : . j - column number to be used
693 : . type - the norm type
694 : - val - the norm
695 :
696 : Note:
697 : Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().
698 :
699 : Level: advanced
700 :
701 : .seealso: BVNormColumnEnd(), BVNormColumn()
702 : @*/
703 2708 : PetscErrorCode BVNormColumnBegin(BV bv,PetscInt j,NormType type,PetscReal *val)
704 : {
705 2708 : PetscSplitReduction *sr;
706 2708 : PetscReal lresult;
707 2708 : MPI_Comm comm;
708 2708 : Vec z;
709 :
710 2708 : PetscFunctionBegin;
711 2708 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
712 10832 : PetscValidLogicalCollectiveInt(bv,j,2);
713 10832 : PetscValidLogicalCollectiveEnum(bv,type,3);
714 2708 : PetscAssertPointer(val,4);
715 2708 : PetscValidType(bv,1);
716 2708 : BVCheckSizes(bv,1);
717 :
718 2708 : PetscCheck(j>=0 && j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %" PetscInt_FMT ", the number of columns is %" PetscInt_FMT,j,bv->m);
719 2708 : PetscCheck(type!=NORM_1_AND_2 || bv->matrix,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
720 :
721 2708 : PetscCall(PetscLogEventBegin(BV_NormVec,bv,0,0,0));
722 2708 : PetscCall(BVGetColumn(bv,j,&z));
723 2708 : if (bv->matrix) PetscCall(BVNorm_Begin_Private(bv,z,type,val)); /* non-standard inner product */
724 2314 : else if (bv->ops->norm_begin) PetscUseTypeMethod(bv,norm_begin,j,type,val);
725 : else {
726 2312 : BVCheckOp(bv,1,norm_local);
727 2312 : PetscCall(PetscObjectGetComm((PetscObject)z,&comm));
728 2312 : PetscCall(PetscSplitReductionGet(comm,&sr));
729 2312 : PetscCheck(sr->state==STATE_BEGIN,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Called before all VecxxxEnd() called");
730 2312 : if (sr->numopsbegin >= sr->maxops) PetscCall(PetscSplitReductionExtend(sr));
731 2312 : sr->invecs[sr->numopsbegin] = (void*)bv;
732 2312 : PetscUseTypeMethod(bv,norm_local,j,type,&lresult);
733 2312 : if (type == NORM_2) lresult = lresult*lresult;
734 2312 : if (type == NORM_MAX) sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_MAX;
735 2312 : else sr->reducetype[sr->numopsbegin] = PETSC_SR_REDUCE_SUM;
736 2312 : sr->lvalues[sr->numopsbegin++] = lresult;
737 : }
738 2708 : PetscCall(BVRestoreColumn(bv,j,&z));
739 2708 : PetscCall(PetscLogEventEnd(BV_NormVec,bv,0,0,0));
740 2708 : PetscFunctionReturn(PETSC_SUCCESS);
741 : }
742 :
743 : /*@
744 : BVNormColumnEnd - Ends a split phase norm computation.
745 :
746 : Input Parameters:
747 : + bv - basis vectors
748 : . j - column number to be used
749 : . type - the norm type
750 : - val - the norm
751 :
752 : Note:
753 : Each call to BVNormColumnBegin() should be paired with a call to BVNormColumnEnd().
754 :
755 : Level: advanced
756 :
757 : .seealso: BVNormColumnBegin(), BVNormColumn()
758 : @*/
759 2708 : PetscErrorCode BVNormColumnEnd(BV bv,PetscInt j,NormType type,PetscReal *val)
760 : {
761 2708 : PetscSplitReduction *sr;
762 2708 : MPI_Comm comm;
763 2708 : Vec z;
764 :
765 2708 : PetscFunctionBegin;
766 2708 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
767 10832 : PetscValidLogicalCollectiveInt(bv,j,2);
768 10832 : PetscValidLogicalCollectiveEnum(bv,type,3);
769 2708 : PetscAssertPointer(val,4);
770 2708 : PetscValidType(bv,1);
771 2708 : BVCheckSizes(bv,1);
772 :
773 2708 : PetscCheck(type!=NORM_1_AND_2,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not available");
774 :
775 2708 : PetscCall(BVGetColumn(bv,j,&z));
776 2708 : if (bv->matrix) PetscCall(BVNorm_End_Private(bv,z,type,val)); /* non-standard inner product */
777 2314 : else if (bv->ops->norm_end) PetscUseTypeMethod(bv,norm_end,j,type,val);
778 : else {
779 2312 : PetscCall(PetscObjectGetComm((PetscObject)z,&comm));
780 2312 : PetscCall(PetscSplitReductionGet(comm,&sr));
781 2312 : PetscCall(PetscSplitReductionEnd(sr));
782 :
783 2312 : PetscCheck(sr->numopsend<sr->numopsbegin,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() more times then VecxxxBegin()");
784 2312 : PetscCheck((void*)bv==sr->invecs[sr->numopsend],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called VecxxxEnd() in a different order or with a different vector than VecxxxBegin()");
785 2312 : PetscCheck(sr->reducetype[sr->numopsend]==PETSC_SR_REDUCE_MAX || type!=NORM_MAX,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Called BVNormEnd(,NORM_MAX,) on a reduction started with VecDotBegin() or NORM_1 or NORM_2");
786 2312 : *val = PetscRealPart(sr->gvalues[sr->numopsend++]);
787 2312 : if (type == NORM_2) *val = PetscSqrtReal(*val);
788 2312 : if (sr->numopsend == sr->numopsbegin) {
789 131 : sr->state = STATE_BEGIN;
790 131 : sr->numopsend = 0;
791 131 : sr->numopsbegin = 0;
792 : }
793 : }
794 2708 : PetscCall(BVRestoreColumn(bv,j,&z));
795 2708 : PetscFunctionReturn(PETSC_SUCCESS);
796 : }
797 :
798 : /*@
799 : BVNormalize - Normalize all columns (starting from the leading ones).
800 :
801 : Collective
802 :
803 : Input Parameters:
804 : + bv - basis vectors
805 : - eigi - (optional) imaginary parts of eigenvalues
806 :
807 : Notes:
808 : On output, all columns will have unit norm. The normalization is done with
809 : respect to the 2-norm, or to the B-norm if a non-standard inner product has
810 : been specified with BVSetMatrix(), see BVNormColumn().
811 :
812 : If the optional argument eigi is passed (taken into account only in real
813 : scalars) it is interpreted as the imaginary parts of the eigenvalues and
814 : the BV is supposed to contain the corresponding eigenvectors. Suppose the
815 : first three values are eigi = { 0, alpha, -alpha }, then the first column
816 : is normalized as usual, but the second and third ones are normalized assuming
817 : that they contain the real and imaginary parts of a complex conjugate pair of
818 : eigenvectors.
819 :
820 : If eigi is passed, the inner-product matrix is ignored.
821 :
822 : If there are leading columns, they are not modified (are assumed to be already
823 : normalized).
824 :
825 : Level: intermediate
826 :
827 : .seealso: BVNormColumn()
828 : @*/
829 3588 : PetscErrorCode BVNormalize(BV bv,PetscScalar *eigi)
830 : {
831 3588 : PetscReal norm;
832 3588 : PetscInt i;
833 :
834 3588 : PetscFunctionBegin;
835 3588 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
836 3588 : PetscValidType(bv,1);
837 3588 : BVCheckSizes(bv,1);
838 :
839 3588 : PetscCall(PetscLogEventBegin(BV_Normalize,bv,0,0,0));
840 3588 : if (bv->matrix && !eigi) {
841 1338 : for (i=bv->l;i<bv->k;i++) {
842 1252 : PetscCall(BVNormColumn(bv,i,NORM_2,&norm));
843 1252 : PetscCall(BVScaleColumn(bv,i,1.0/norm));
844 : }
845 3502 : } else PetscTryTypeMethod(bv,normalize,eigi);
846 3588 : PetscCall(PetscLogEventEnd(BV_Normalize,bv,0,0,0));
847 3588 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
848 3588 : PetscFunctionReturn(PETSC_SUCCESS);
849 : }
850 :
851 : /*
852 : Compute Y^H*A*X: right part column by column (with MatMult) and bottom
853 : part row by row (with MatMultHermitianTranspose); result placed in marray[*,ldm]
854 : */
855 18 : static inline PetscErrorCode BVMatProject_Vec(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)
856 : {
857 18 : PetscInt i,j,lx,ly,kx,ky,ulim;
858 18 : Vec z,f;
859 :
860 18 : PetscFunctionBegin;
861 18 : lx = X->l; kx = X->k;
862 18 : ly = Y->l; ky = Y->k;
863 18 : PetscCall(BVCreateVec(X,&f));
864 18 : BVCheckOp(Y,3,dotvec);
865 72 : for (j=lx;j<kx;j++) {
866 54 : PetscCall(BVGetColumn(X,j,&z));
867 54 : PetscCall(MatMult(A,z,f));
868 54 : PetscCall(BVRestoreColumn(X,j,&z));
869 54 : ulim = PetscMin(ly+(j-lx)+1,ky);
870 54 : Y->l = 0; Y->k = ulim;
871 54 : PetscUseTypeMethod(Y,dotvec,f,marray+j*ldm);
872 54 : if (symm) {
873 90 : for (i=0;i<j;i++) marray[j+i*ldm] = PetscConj(marray[i+j*ldm]);
874 : }
875 : }
876 18 : if (!symm) {
877 12 : BVCheckOp(X,1,dotvec);
878 12 : PetscCall(BV_AllocateCoeffs(Y));
879 48 : for (j=ly;j<ky;j++) {
880 36 : PetscCall(BVGetColumn(Y,j,&z));
881 36 : PetscCall(MatMultHermitianTranspose(A,z,f));
882 36 : PetscCall(BVRestoreColumn(Y,j,&z));
883 36 : ulim = PetscMin(lx+(j-ly),kx);
884 36 : X->l = 0; X->k = ulim;
885 36 : PetscUseTypeMethod(X,dotvec,f,Y->h);
886 180 : for (i=0;i<ulim;i++) marray[j+i*ldm] = PetscConj(Y->h[i]);
887 : }
888 : }
889 18 : PetscCall(VecDestroy(&f));
890 18 : X->l = lx; X->k = kx;
891 18 : Y->l = ly; Y->k = ky;
892 18 : PetscFunctionReturn(PETSC_SUCCESS);
893 : }
894 :
895 : /*
896 : Compute Y^H*A*X= [ -- | Y0'*W1 ]
897 : [ Y1'*W0 | Y1'*W1 ]
898 : Allocates auxiliary BV to store the result of A*X, then one BVDot
899 : call for top-right part and another one for bottom part;
900 : result placed in marray[*,ldm]
901 : */
902 2277 : static inline PetscErrorCode BVMatProject_MatMult(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm)
903 : {
904 2277 : PetscInt j,lx,ly,kx,ky;
905 2277 : const PetscScalar *harray;
906 2277 : Mat H;
907 2277 : BV W;
908 :
909 2277 : PetscFunctionBegin;
910 2277 : lx = X->l; kx = X->k;
911 2277 : ly = Y->l; ky = Y->k;
912 2277 : PetscCall(BVDuplicate(X,&W));
913 2277 : X->l = 0; X->k = kx;
914 2277 : W->l = 0; W->k = kx;
915 2277 : PetscCall(BVMatMult(X,A,W));
916 :
917 : /* top-right part, Y0'*AX1 */
918 2277 : if (ly>0 && lx<kx) {
919 0 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H));
920 0 : W->l = lx; W->k = kx;
921 0 : Y->l = 0; Y->k = ly;
922 0 : PetscCall(BVDot(W,Y,H));
923 0 : PetscCall(MatDenseGetArrayRead(H,&harray));
924 0 : for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+j*ly,ly));
925 0 : PetscCall(MatDenseRestoreArrayRead(H,&harray));
926 0 : PetscCall(MatDestroy(&H));
927 : }
928 :
929 : /* bottom part, Y1'*AX */
930 2277 : if (kx>0 && ly<ky) {
931 2277 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H));
932 2277 : W->l = 0; W->k = kx;
933 2277 : Y->l = ly; Y->k = ky;
934 2277 : PetscCall(BVDot(W,Y,H));
935 2277 : PetscCall(MatDenseGetArrayRead(H,&harray));
936 23112 : for (j=0;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm+ly,harray+j*ky+ly,ky-ly));
937 2277 : PetscCall(MatDenseRestoreArrayRead(H,&harray));
938 2277 : PetscCall(MatDestroy(&H));
939 : }
940 2277 : PetscCall(BVDestroy(&W));
941 2277 : X->l = lx; X->k = kx;
942 2277 : Y->l = ly; Y->k = ky;
943 2277 : PetscFunctionReturn(PETSC_SUCCESS);
944 : }
945 :
946 : /*
947 : Compute Y^H*A*X= [ -- | Y0'*W1 ]
948 : [ Y1'*W0 | Y1'*W1 ]
949 : First stage: allocate auxiliary BV to store A*X1, one BVDot for right part;
950 : Second stage: resize BV to accommodate A'*Y1, then call BVDot for transpose of
951 : bottom-left part; result placed in marray[*,ldm]
952 : */
953 586 : static inline PetscErrorCode BVMatProject_MatMult_2(BV X,Mat A,BV Y,PetscScalar *marray,PetscInt ldm,PetscBool symm)
954 : {
955 586 : PetscInt i,j,lx,ly,kx,ky;
956 586 : const PetscScalar *harray;
957 586 : Mat H;
958 586 : BV W;
959 :
960 586 : PetscFunctionBegin;
961 586 : lx = X->l; kx = X->k;
962 586 : ly = Y->l; ky = Y->k;
963 :
964 : /* right part, Y'*AX1 */
965 586 : PetscCall(BVDuplicateResize(X,kx-lx,&W));
966 586 : if (ky>0 && lx<kx) {
967 586 : PetscCall(BVMatMult(X,A,W));
968 586 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx-lx,NULL,&H));
969 586 : Y->l = 0; Y->k = ky;
970 586 : PetscCall(BVDot(W,Y,H));
971 586 : PetscCall(MatDenseGetArrayRead(H,&harray));
972 1248 : for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+(j-lx)*ky,ky));
973 586 : PetscCall(MatDenseRestoreArrayRead(H,&harray));
974 586 : PetscCall(MatDestroy(&H));
975 : }
976 :
977 : /* bottom-left part, Y1'*AX0 */
978 586 : if (lx>0 && ly<ky) {
979 511 : if (symm) {
980 : /* do not compute, just copy symmetric elements */
981 906 : for (i=ly;i<ky;i++) {
982 2469 : for (j=0;j<lx;j++) marray[i+j*ldm] = PetscConj(marray[j+i*ldm]);
983 : }
984 : } else {
985 67 : PetscCall(BVResize(W,ky-ly,PETSC_FALSE));
986 67 : Y->l = ly; Y->k = ky;
987 67 : PetscCall(BVMatMultHermitianTranspose(Y,A,W));
988 67 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,lx,ky-ly,NULL,&H));
989 67 : X->l = 0; X->k = lx;
990 67 : PetscCall(BVDot(W,X,H));
991 67 : PetscCall(MatDenseGetArrayRead(H,&harray));
992 174 : for (i=0;i<ky-ly;i++) {
993 428 : for (j=0;j<lx;j++) {
994 321 : marray[i+j*ldm+ly] = PetscConj(harray[j+i*(ky-ly)]);
995 : }
996 : }
997 67 : PetscCall(MatDenseRestoreArrayRead(H,&harray));
998 67 : PetscCall(MatDestroy(&H));
999 : }
1000 : }
1001 586 : PetscCall(BVDestroy(&W));
1002 586 : X->l = lx; X->k = kx;
1003 586 : Y->l = ly; Y->k = ky;
1004 586 : PetscFunctionReturn(PETSC_SUCCESS);
1005 : }
1006 :
1007 : /*
1008 : Compute Y^H*X = [ -- | Y0'*X1 ] (X contains A*X):
1009 : [ Y1'*X0 | Y1'*X1 ]
1010 : one BVDot call for top-right part and another one for bottom part;
1011 : result placed in marray[*,ldm]
1012 : */
1013 14325 : static inline PetscErrorCode BVMatProject_Dot(BV X,BV Y,PetscScalar *marray,PetscInt ldm)
1014 : {
1015 14325 : PetscInt j,lx,ly,kx,ky;
1016 14325 : const PetscScalar *harray;
1017 14325 : Mat H;
1018 :
1019 14325 : PetscFunctionBegin;
1020 14325 : lx = X->l; kx = X->k;
1021 14325 : ly = Y->l; ky = Y->k;
1022 :
1023 : /* top-right part, Y0'*X1 */
1024 14325 : if (ly>0 && lx<kx) {
1025 8097 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ly,kx,NULL,&H));
1026 8097 : X->l = lx; X->k = kx;
1027 8097 : Y->l = 0; Y->k = ly;
1028 8097 : PetscCall(BVDot(X,Y,H));
1029 8097 : PetscCall(MatDenseGetArrayRead(H,&harray));
1030 19871 : for (j=lx;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm,harray+j*ly,ly));
1031 8097 : PetscCall(MatDenseRestoreArrayRead(H,&harray));
1032 8097 : PetscCall(MatDestroy(&H));
1033 : }
1034 :
1035 : /* bottom part, Y1'*X */
1036 14325 : if (kx>0 && ly<ky) {
1037 14325 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,ky,kx,NULL,&H));
1038 14325 : X->l = 0; X->k = kx;
1039 14325 : Y->l = ly; Y->k = ky;
1040 14325 : PetscCall(BVDot(X,Y,H));
1041 14325 : PetscCall(MatDenseGetArrayRead(H,&harray));
1042 136744 : for (j=0;j<kx;j++) PetscCall(PetscArraycpy(marray+j*ldm+ly,harray+j*ky+ly,ky-ly));
1043 14325 : PetscCall(MatDenseRestoreArrayRead(H,&harray));
1044 14325 : PetscCall(MatDestroy(&H));
1045 : }
1046 14325 : X->l = lx; X->k = kx;
1047 14325 : Y->l = ly; Y->k = ky;
1048 14325 : PetscFunctionReturn(PETSC_SUCCESS);
1049 : }
1050 :
1051 : /*@
1052 : BVMatProject - Computes the projection of a matrix onto a subspace.
1053 :
1054 : Collective
1055 :
1056 : Input Parameters:
1057 : + X - basis vectors
1058 : . A - (optional) matrix to be projected
1059 : - Y - left basis vectors, can be equal to X
1060 :
1061 : Output Parameter:
1062 : . M - the resulting matrix
1063 :
1064 : Notes:
1065 : If A=NULL, then it is assumed that X already contains A*X.
1066 :
1067 : This operation is similar to BVDot(), with important differences.
1068 : The goal is to compute the matrix resulting from the orthogonal projection
1069 : of A onto the subspace spanned by the columns of X, M = X^H*A*X, or the
1070 : oblique projection onto X along Y, M = Y^H*A*X.
1071 :
1072 : A difference with respect to BVDot() is that the standard inner product
1073 : is always used, regardless of a non-standard inner product being specified
1074 : with BVSetMatrix().
1075 :
1076 : On entry, M must be a sequential dense Mat with dimensions ky,kx at least,
1077 : where ky (resp. kx) is the number of active columns of Y (resp. X).
1078 : Another difference with respect to BVDot() is that all entries of M are
1079 : computed except the leading ly,lx part, where ly (resp. lx) is the
1080 : number of leading columns of Y (resp. X). Hence, the leading columns of
1081 : X and Y participate in the computation, as opposed to BVDot().
1082 : The leading part of M is assumed to be already available from previous
1083 : computations.
1084 :
1085 : In the orthogonal projection case, Y=X, some computation can be saved if
1086 : A is real symmetric (or complex Hermitian). In order to exploit this
1087 : property, the symmetry flag of A must be set with MatSetOption().
1088 :
1089 : Level: intermediate
1090 :
1091 : .seealso: BVDot(), BVSetActiveColumns(), BVSetMatrix()
1092 : @*/
1093 17206 : PetscErrorCode BVMatProject(BV X,Mat A,BV Y,Mat M)
1094 : {
1095 17206 : PetscBool set,flg,symm=PETSC_FALSE;
1096 17206 : PetscInt m,n,ldm;
1097 17206 : PetscScalar *marray;
1098 17206 : Mat Xmatrix,Ymatrix;
1099 17206 : PetscObjectId idx,idy;
1100 :
1101 17206 : PetscFunctionBegin;
1102 17206 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
1103 17206 : if (A) PetscValidHeaderSpecific(A,MAT_CLASSID,2);
1104 17206 : PetscValidHeaderSpecific(Y,BV_CLASSID,3);
1105 17206 : PetscValidHeaderSpecific(M,MAT_CLASSID,4);
1106 17206 : PetscValidType(X,1);
1107 17206 : BVCheckSizes(X,1);
1108 17206 : if (A) {
1109 2881 : PetscValidType(A,2);
1110 2881 : PetscCheckSameComm(X,1,A,2);
1111 : }
1112 17206 : PetscValidType(Y,3);
1113 17206 : BVCheckSizes(Y,3);
1114 17206 : PetscValidType(M,4);
1115 17206 : PetscCheckSameTypeAndComm(X,1,Y,3);
1116 17206 : PetscCheckTypeNames(M,MATSEQDENSE,MATSEQDENSECUDA);
1117 :
1118 17206 : PetscCall(MatGetSize(M,&m,&n));
1119 17206 : PetscCall(MatDenseGetLDA(M,&ldm));
1120 17206 : PetscCheck(m>=Y->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,Y->k);
1121 17206 : PetscCheck(n>=X->k,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Matrix M has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,X->k);
1122 17206 : PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);
1123 :
1124 17206 : PetscCall(PetscLogEventBegin(BV_MatProject,X,A,Y,0));
1125 : /* temporarily set standard inner product */
1126 17206 : Xmatrix = X->matrix;
1127 17206 : Ymatrix = Y->matrix;
1128 17206 : X->matrix = Y->matrix = NULL;
1129 :
1130 17206 : PetscCall(PetscObjectGetId((PetscObject)X,&idx));
1131 17206 : PetscCall(PetscObjectGetId((PetscObject)Y,&idy));
1132 17206 : if (A && idx==idy) { /* check symmetry of M=X'AX */
1133 2866 : PetscCall(MatIsHermitianKnown(A,&set,&flg));
1134 2866 : symm = set? flg: PETSC_FALSE;
1135 : }
1136 :
1137 17206 : PetscCall(MatDenseGetArray(M,&marray));
1138 :
1139 17206 : if (A) {
1140 2881 : if (X->vmm==BV_MATMULT_VECS) {
1141 : /* perform computation column by column */
1142 18 : PetscCall(BVMatProject_Vec(X,A,Y,marray,ldm,symm));
1143 : } else {
1144 : /* use BVMatMult, then BVDot */
1145 2863 : PetscCall(MatHasOperation(A,MATOP_MULT_TRANSPOSE,&flg));
1146 2863 : if (symm || (flg && X->l>=X->k/2 && Y->l>=Y->k/2)) PetscCall(BVMatProject_MatMult_2(X,A,Y,marray,ldm,symm));
1147 2277 : else PetscCall(BVMatProject_MatMult(X,A,Y,marray,ldm));
1148 : }
1149 : } else {
1150 : /* use BVDot on subblocks */
1151 14325 : PetscCall(BVMatProject_Dot(X,Y,marray,ldm));
1152 : }
1153 :
1154 17206 : PetscCall(MatDenseRestoreArray(M,&marray));
1155 17206 : PetscCall(PetscLogEventEnd(BV_MatProject,X,A,Y,0));
1156 : /* restore non-standard inner product */
1157 17206 : X->matrix = Xmatrix;
1158 17206 : Y->matrix = Ymatrix;
1159 17206 : PetscFunctionReturn(PETSC_SUCCESS);
1160 : }
|
