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 : Basic BV routines
12 : */
13 :
14 : #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
15 :
16 : PetscBool BVRegisterAllCalled = PETSC_FALSE;
17 : PetscFunctionList BVList = NULL;
18 :
19 : /*@
20 : BVSetType - Selects the type for the BV object.
21 :
22 : Logically Collective
23 :
24 : Input Parameters:
25 : + bv - the basis vectors context
26 : - type - a known type
27 :
28 : Options Database Key:
29 : . -bv_type <type> - Sets BV type
30 :
31 : Level: intermediate
32 :
33 : .seealso: BVGetType()
34 : @*/
35 14220 : PetscErrorCode BVSetType(BV bv,BVType type)
36 : {
37 14220 : PetscErrorCode (*r)(BV);
38 14220 : PetscBool match;
39 :
40 14220 : PetscFunctionBegin;
41 14220 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
42 14220 : PetscAssertPointer(type,2);
43 :
44 14220 : PetscCall(PetscObjectTypeCompare((PetscObject)bv,type,&match));
45 14220 : if (match) PetscFunctionReturn(PETSC_SUCCESS);
46 14217 : PetscCall(PetscStrcmp(type,BVTENSOR,&match));
47 14217 : PetscCheck(!match,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Use BVCreateTensor() to create a BV of type tensor");
48 :
49 14217 : PetscCall(PetscFunctionListFind(BVList,type,&r));
50 14217 : PetscCheck(r,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unable to find requested BV type %s",type);
51 :
52 14217 : PetscTryTypeMethod(bv,destroy);
53 14217 : PetscCall(PetscMemzero(bv->ops,sizeof(struct _BVOps)));
54 :
55 14217 : PetscCall(PetscObjectChangeTypeName((PetscObject)bv,type));
56 14217 : if (bv->n < 0 && bv->N < 0) {
57 1639 : bv->ops->create = r;
58 : } else {
59 12578 : PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
60 12578 : PetscCall((*r)(bv));
61 12578 : PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
62 : }
63 14217 : PetscFunctionReturn(PETSC_SUCCESS);
64 : }
65 :
66 : /*@
67 : BVGetType - Gets the BV type name (as a string) from the BV context.
68 :
69 : Not Collective
70 :
71 : Input Parameter:
72 : . bv - the basis vectors context
73 :
74 : Output Parameter:
75 : . type - name of the type of basis vectors
76 :
77 : Level: intermediate
78 :
79 : .seealso: BVSetType()
80 : @*/
81 160 : PetscErrorCode BVGetType(BV bv,BVType *type)
82 : {
83 160 : PetscFunctionBegin;
84 160 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
85 160 : PetscAssertPointer(type,2);
86 160 : *type = ((PetscObject)bv)->type_name;
87 160 : PetscFunctionReturn(PETSC_SUCCESS);
88 : }
89 :
90 : /*@
91 : BVSetSizes - Sets the local and global sizes, and the number of columns.
92 :
93 : Collective
94 :
95 : Input Parameters:
96 : + bv - the basis vectors
97 : . n - the local size (or PETSC_DECIDE to have it set)
98 : . N - the global size (or PETSC_DECIDE)
99 : - m - the number of columns
100 :
101 : Notes:
102 : n and N cannot be both PETSC_DECIDE.
103 : If one processor calls this with N of PETSC_DECIDE then all processors must,
104 : otherwise the program will hang.
105 :
106 : Level: beginner
107 :
108 : .seealso: BVSetSizesFromVec(), BVGetSizes(), BVResize()
109 : @*/
110 451 : PetscErrorCode BVSetSizes(BV bv,PetscInt n,PetscInt N,PetscInt m)
111 : {
112 451 : PetscInt ma;
113 451 : PetscMPIInt size;
114 :
115 451 : PetscFunctionBegin;
116 451 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
117 969 : if (N >= 0) PetscValidLogicalCollectiveInt(bv,N,3);
118 1353 : PetscValidLogicalCollectiveInt(bv,m,4);
119 451 : PetscCheck(N<0 || n<=N,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local size %" PetscInt_FMT " cannot be larger than global size %" PetscInt_FMT,n,N);
120 451 : PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
121 451 : PetscCheck((bv->n<0 && bv->N<0) || (bv->n==n && bv->N==N),PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change/reset vector sizes to %" PetscInt_FMT " local %" PetscInt_FMT " global after previously setting them to %" PetscInt_FMT " local %" PetscInt_FMT " global",n,N,bv->n,bv->N);
122 451 : PetscCheck(bv->m<=0 || bv->m==m,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot change the number of columns to %" PetscInt_FMT " after previously setting it to %" PetscInt_FMT "; use BVResize()",m,bv->m);
123 451 : PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
124 451 : bv->n = n;
125 451 : bv->N = N;
126 451 : bv->m = m;
127 451 : bv->k = m;
128 : /* create layout and get actual dimensions */
129 451 : PetscCall(PetscLayoutCreate(PetscObjectComm((PetscObject)bv),&bv->map));
130 451 : PetscCall(PetscLayoutSetSize(bv->map,bv->N));
131 451 : PetscCall(PetscLayoutSetLocalSize(bv->map,bv->n));
132 451 : PetscCall(PetscLayoutSetUp(bv->map));
133 451 : PetscCall(PetscLayoutGetSize(bv->map,&bv->N));
134 451 : PetscCall(PetscLayoutGetLocalSize(bv->map,&bv->n));
135 451 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
136 541 : PetscCall(BVSetVecType(bv,(size==1)?VECSEQ:VECMPI));
137 451 : if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
138 0 : PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
139 0 : PetscCheck(bv->n==ma,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %" PetscInt_FMT " does not match that of matrix given at BVSetMatrix %" PetscInt_FMT,bv->n,ma);
140 : }
141 451 : if (bv->ops->create) {
142 134 : PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
143 134 : PetscUseTypeMethod(bv,create);
144 134 : PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
145 134 : bv->ops->create = NULL;
146 134 : bv->defersfo = PETSC_FALSE;
147 : }
148 451 : PetscFunctionReturn(PETSC_SUCCESS);
149 : }
150 :
151 : /*@
152 : BVSetSizesFromVec - Sets the local and global sizes, and the number of columns.
153 : Local and global sizes are specified indirectly by passing a template vector.
154 :
155 : Collective
156 :
157 : Input Parameters:
158 : + bv - the basis vectors
159 : . t - the template vector
160 : - m - the number of columns
161 :
162 : Level: beginner
163 :
164 : .seealso: BVSetSizes(), BVGetSizes(), BVResize()
165 : @*/
166 1978 : PetscErrorCode BVSetSizesFromVec(BV bv,Vec t,PetscInt m)
167 : {
168 1978 : PetscInt ma;
169 1978 : PetscLayout map;
170 1978 : VecType vtype;
171 :
172 1978 : PetscFunctionBegin;
173 1978 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
174 1978 : PetscValidHeaderSpecific(t,VEC_CLASSID,2);
175 1978 : PetscCheckSameComm(bv,1,t,2);
176 5934 : PetscValidLogicalCollectiveInt(bv,m,3);
177 1978 : PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
178 1978 : PetscCheck(!bv->map,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Vector layout was already defined by a previous call to BVSetSizes/FromVec");
179 1978 : PetscCall(VecGetType(t,&vtype));
180 1978 : PetscCall(BVSetVecType(bv,vtype));
181 1978 : PetscCall(VecGetLayout(t,&map));
182 1978 : PetscCall(PetscLayoutReference(map,&bv->map));
183 1978 : PetscCall(VecGetSize(t,&bv->N));
184 1978 : PetscCall(VecGetLocalSize(t,&bv->n));
185 1978 : if (bv->matrix) { /* check compatible dimensions of user-provided matrix */
186 5 : PetscCall(MatGetLocalSize(bv->matrix,&ma,NULL));
187 5 : PetscCheck(bv->n==ma,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Local dimension %" PetscInt_FMT " does not match that of matrix given at BVSetMatrix %" PetscInt_FMT,bv->n,ma);
188 : }
189 1978 : bv->m = m;
190 1978 : bv->k = m;
191 1978 : if (bv->ops->create) {
192 1498 : PetscUseTypeMethod(bv,create);
193 1498 : bv->ops->create = NULL;
194 1498 : bv->defersfo = PETSC_FALSE;
195 : }
196 1978 : PetscFunctionReturn(PETSC_SUCCESS);
197 : }
198 :
199 : /*@
200 : BVGetSizes - Returns the local and global sizes, and the number of columns.
201 :
202 : Not Collective
203 :
204 : Input Parameter:
205 : . bv - the basis vectors
206 :
207 : Output Parameters:
208 : + n - the local size
209 : . N - the global size
210 : - m - the number of columns
211 :
212 : Note:
213 : Normal usage requires that bv has already been given its sizes, otherwise
214 : the call fails. However, this function can also be used to determine if
215 : a BV object has been initialized completely (sizes and type). For this,
216 : call with n=NULL and N=NULL, then a return value of m=0 indicates that
217 : the BV object is not ready for use yet.
218 :
219 : Level: beginner
220 :
221 : .seealso: BVSetSizes(), BVSetSizesFromVec()
222 : @*/
223 7373 : PetscErrorCode BVGetSizes(BV bv,PetscInt *n,PetscInt *N,PetscInt *m)
224 : {
225 7373 : PetscFunctionBegin;
226 7373 : if (!bv) {
227 93 : if (m && !n && !N) *m = 0;
228 93 : PetscFunctionReturn(PETSC_SUCCESS);
229 : }
230 7280 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
231 7280 : if (n || N) BVCheckSizes(bv,1);
232 7280 : if (n) *n = bv->n;
233 7280 : if (N) *N = bv->N;
234 7280 : if (m) *m = bv->m;
235 7280 : if (m && !n && !N && !((PetscObject)bv)->type_name) *m = 0;
236 7280 : PetscFunctionReturn(PETSC_SUCCESS);
237 : }
238 :
239 : /*@
240 : BVSetNumConstraints - Set the number of constraints.
241 :
242 : Logically Collective
243 :
244 : Input Parameters:
245 : + V - basis vectors
246 : - nc - number of constraints
247 :
248 : Notes:
249 : This function sets the number of constraints to nc and marks all remaining
250 : columns as regular. Normal user would call BVInsertConstraints() instead.
251 :
252 : If nc is smaller than the previously set value, then some of the constraints
253 : are discarded. In particular, using nc=0 removes all constraints preserving
254 : the content of regular columns.
255 :
256 : Level: developer
257 :
258 : .seealso: BVInsertConstraints()
259 : @*/
260 62 : PetscErrorCode BVSetNumConstraints(BV V,PetscInt nc)
261 : {
262 62 : PetscInt total,diff,i;
263 62 : Vec x,y;
264 :
265 62 : PetscFunctionBegin;
266 62 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
267 186 : PetscValidLogicalCollectiveInt(V,nc,2);
268 62 : PetscCheck(nc>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Number of constraints (given %" PetscInt_FMT ") cannot be negative",nc);
269 62 : PetscValidType(V,1);
270 62 : BVCheckSizes(V,1);
271 62 : PetscCheck(V->ci[0]==-V->nc-1 && V->ci[1]==-V->nc-1,PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Cannot call BVSetNumConstraints after BVGetColumn");
272 :
273 62 : diff = nc-V->nc;
274 62 : if (!diff) PetscFunctionReturn(PETSC_SUCCESS);
275 26 : total = V->nc+V->m;
276 26 : PetscCheck(total-nc>0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Not enough columns for the given nc value");
277 26 : if (diff<0) { /* lessen constraints, shift contents of BV */
278 269 : for (i=0;i<V->m;i++) {
279 255 : PetscCall(BVGetColumn(V,i,&x));
280 255 : PetscCall(BVGetColumn(V,i+diff,&y));
281 255 : PetscCall(VecCopy(x,y));
282 255 : PetscCall(BVRestoreColumn(V,i,&x));
283 255 : PetscCall(BVRestoreColumn(V,i+diff,&y));
284 : }
285 : }
286 26 : V->nc = nc;
287 26 : V->ci[0] = -V->nc-1;
288 26 : V->ci[1] = -V->nc-1;
289 26 : V->m = total-nc;
290 26 : V->l = PetscMin(V->l,V->m);
291 26 : V->k = PetscMin(V->k,V->m);
292 26 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
293 26 : PetscFunctionReturn(PETSC_SUCCESS);
294 : }
295 :
296 : /*@
297 : BVGetNumConstraints - Returns the number of constraints.
298 :
299 : Not Collective
300 :
301 : Input Parameter:
302 : . bv - the basis vectors
303 :
304 : Output Parameters:
305 : . nc - the number of constraints
306 :
307 : Level: advanced
308 :
309 : .seealso: BVGetSizes(), BVInsertConstraints()
310 : @*/
311 48 : PetscErrorCode BVGetNumConstraints(BV bv,PetscInt *nc)
312 : {
313 48 : PetscFunctionBegin;
314 48 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
315 48 : PetscAssertPointer(nc,2);
316 48 : *nc = bv->nc;
317 48 : PetscFunctionReturn(PETSC_SUCCESS);
318 : }
319 :
320 : /*@
321 : BVResize - Change the number of columns.
322 :
323 : Collective
324 :
325 : Input Parameters:
326 : + bv - the basis vectors
327 : . m - the new number of columns
328 : - copy - a flag indicating whether current values should be kept
329 :
330 : Note:
331 : Internal storage is reallocated. If the copy flag is set to true, then
332 : the contents are copied to the leading part of the new space.
333 :
334 : Level: advanced
335 :
336 : .seealso: BVSetSizes(), BVSetSizesFromVec()
337 : @*/
338 497 : PetscErrorCode BVResize(BV bv,PetscInt m,PetscBool copy)
339 : {
340 497 : PetscScalar *array;
341 497 : const PetscScalar *omega;
342 497 : Vec v;
343 :
344 497 : PetscFunctionBegin;
345 497 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
346 1491 : PetscValidLogicalCollectiveInt(bv,m,2);
347 1491 : PetscValidLogicalCollectiveBool(bv,copy,3);
348 497 : PetscValidType(bv,1);
349 497 : PetscCheck(m>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Number of columns %" PetscInt_FMT " must be positive",m);
350 497 : PetscCheck(!bv->nc || bv->issplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot resize a BV with constraints");
351 497 : if (bv->m == m) PetscFunctionReturn(PETSC_SUCCESS);
352 153 : BVCheckOp(bv,1,resize);
353 :
354 153 : PetscCall(PetscLogEventBegin(BV_Create,bv,0,0,0));
355 153 : PetscUseTypeMethod(bv,resize,m,copy);
356 153 : PetscCall(VecDestroy(&bv->buffer));
357 153 : PetscCall(BVDestroy(&bv->cached));
358 153 : PetscCall(PetscFree2(bv->h,bv->c));
359 153 : if (bv->omega) {
360 2 : if (bv->cuda) {
361 : #if defined(PETSC_HAVE_CUDA)
362 : PetscCall(VecCreateSeqCUDA(PETSC_COMM_SELF,m,&v));
363 : #else
364 0 : SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
365 : #endif
366 2 : } else if (bv->hip) {
367 : #if defined(PETSC_HAVE_HIP)
368 : PetscCall(VecCreateSeqHIP(PETSC_COMM_SELF,m,&v));
369 : #else
370 0 : SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_PLIB,"Something wrong happened");
371 : #endif
372 2 : } else PetscCall(VecCreateSeq(PETSC_COMM_SELF,m,&v));
373 2 : if (copy) {
374 2 : PetscCall(VecGetArray(v,&array));
375 2 : PetscCall(VecGetArrayRead(bv->omega,&omega));
376 2 : PetscCall(PetscArraycpy(array,omega,PetscMin(m,bv->m)));
377 2 : PetscCall(VecRestoreArrayRead(bv->omega,&omega));
378 2 : PetscCall(VecRestoreArray(v,&array));
379 0 : } else PetscCall(VecSet(v,1.0));
380 2 : PetscCall(VecDestroy(&bv->omega));
381 2 : bv->omega = v;
382 : }
383 153 : bv->m = m;
384 153 : bv->k = PetscMin(bv->k,m);
385 153 : bv->l = PetscMin(bv->l,m);
386 153 : PetscCall(PetscLogEventEnd(BV_Create,bv,0,0,0));
387 153 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
388 153 : PetscFunctionReturn(PETSC_SUCCESS);
389 : }
390 :
391 : /*@
392 : BVSetActiveColumns - Specify the columns that will be involved in operations.
393 :
394 : Logically Collective
395 :
396 : Input Parameters:
397 : + bv - the basis vectors context
398 : . l - number of leading columns
399 : - k - number of active columns
400 :
401 : Notes:
402 : In operations such as BVMult() or BVDot(), only the first k columns are
403 : considered. This is useful when the BV is filled from left to right, so
404 : the last m-k columns do not have relevant information.
405 :
406 : Also in operations such as BVMult() or BVDot(), the first l columns are
407 : normally not included in the computation. See the manpage of each
408 : operation.
409 :
410 : In orthogonalization operations, the first l columns are treated
411 : differently, they participate in the orthogonalization but the computed
412 : coefficients are not stored.
413 :
414 : Use PETSC_CURRENT to leave any of the values unchanged. Use PETSC_DETERMINE
415 : to set l to the minimum value (0) and k to the maximum (m).
416 :
417 : Level: intermediate
418 :
419 : .seealso: BVGetActiveColumns(), BVSetSizes()
420 : @*/
421 384879 : PetscErrorCode BVSetActiveColumns(BV bv,PetscInt l,PetscInt k)
422 : {
423 384879 : PetscFunctionBegin;
424 384879 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
425 1154637 : PetscValidLogicalCollectiveInt(bv,l,2);
426 1154637 : PetscValidLogicalCollectiveInt(bv,k,3);
427 384879 : BVCheckSizes(bv,1);
428 384879 : if (PetscUnlikely(k == PETSC_DETERMINE)) {
429 0 : bv->k = bv->m;
430 384879 : } else if (k != PETSC_CURRENT) {
431 384879 : PetscCheck(k>=0 && k<=bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of k (%" PetscInt_FMT "). Must be between 0 and m (%" PetscInt_FMT ")",k,bv->m);
432 384879 : bv->k = k;
433 : }
434 384879 : if (PetscUnlikely(l == PETSC_DETERMINE)) {
435 0 : bv->l = 0;
436 384879 : } else if (l != PETSC_CURRENT) {
437 384879 : PetscCheck(l>=0 && l<=bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of l (%" PetscInt_FMT "). Must be between 0 and k (%" PetscInt_FMT ")",l,bv->k);
438 384879 : bv->l = l;
439 : }
440 384879 : PetscFunctionReturn(PETSC_SUCCESS);
441 : }
442 :
443 : /*@
444 : BVGetActiveColumns - Returns the current active dimensions.
445 :
446 : Not Collective
447 :
448 : Input Parameter:
449 : . bv - the basis vectors context
450 :
451 : Output Parameters:
452 : + l - number of leading columns
453 : - k - number of active columns
454 :
455 : Level: intermediate
456 :
457 : .seealso: BVSetActiveColumns()
458 : @*/
459 105330 : PetscErrorCode BVGetActiveColumns(BV bv,PetscInt *l,PetscInt *k)
460 : {
461 105330 : PetscFunctionBegin;
462 105330 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
463 105330 : if (l) *l = bv->l;
464 105330 : if (k) *k = bv->k;
465 105330 : PetscFunctionReturn(PETSC_SUCCESS);
466 : }
467 :
468 : /*@
469 : BVSetMatrix - Specifies the inner product to be used in orthogonalization.
470 :
471 : Collective
472 :
473 : Input Parameters:
474 : + bv - the basis vectors context
475 : . B - a symmetric matrix (may be NULL)
476 : - indef - a flag indicating if the matrix is indefinite
477 :
478 : Notes:
479 : This is used to specify a non-standard inner product, whose matrix
480 : representation is given by B. Then, all inner products required during
481 : orthogonalization are computed as (x,y)_B=y^H*B*x rather than the
482 : standard form (x,y)=y^H*x.
483 :
484 : Matrix B must be real symmetric (or complex Hermitian). A genuine inner
485 : product requires that B is also positive (semi-)definite. However, we
486 : also allow for an indefinite B (setting indef=PETSC_TRUE), in which
487 : case the orthogonalization uses an indefinite inner product.
488 :
489 : This affects operations BVDot(), BVNorm(), BVOrthogonalize(), and variants.
490 :
491 : Setting B=NULL has the same effect as if the identity matrix was passed.
492 :
493 : Level: advanced
494 :
495 : .seealso: BVGetMatrix(), BVDot(), BVNorm(), BVOrthogonalize(), BVSetDefiniteTolerance()
496 : @*/
497 5081 : PetscErrorCode BVSetMatrix(BV bv,Mat B,PetscBool indef)
498 : {
499 5081 : PetscInt m,n;
500 :
501 5081 : PetscFunctionBegin;
502 5081 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
503 15243 : PetscValidLogicalCollectiveBool(bv,indef,3);
504 5081 : if (B!=bv->matrix || (B && ((PetscObject)B)->id!=((PetscObject)bv->matrix)->id) || indef!=bv->indef) {
505 2408 : if (B) {
506 1286 : PetscValidHeaderSpecific(B,MAT_CLASSID,2);
507 1286 : PetscCall(MatGetLocalSize(B,&m,&n));
508 1286 : PetscCheck(m==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Matrix must be square");
509 1286 : PetscCheck(!bv->m || bv->n==n,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension BV %" PetscInt_FMT ", Mat %" PetscInt_FMT,bv->n,n);
510 : }
511 2408 : if (B) PetscCall(PetscObjectReference((PetscObject)B));
512 2408 : PetscCall(MatDestroy(&bv->matrix));
513 2408 : bv->matrix = B;
514 2408 : bv->indef = indef;
515 2408 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
516 2408 : if (bv->Bx) PetscCall(PetscObjectStateIncrease((PetscObject)bv->Bx));
517 2408 : if (bv->cached) PetscCall(PetscObjectStateIncrease((PetscObject)bv->cached));
518 : }
519 5081 : PetscFunctionReturn(PETSC_SUCCESS);
520 : }
521 :
522 : /*@
523 : BVGetMatrix - Retrieves the matrix representation of the inner product.
524 :
525 : Not Collective
526 :
527 : Input Parameter:
528 : . bv - the basis vectors context
529 :
530 : Output Parameters:
531 : + B - the matrix of the inner product (may be NULL)
532 : - indef - the flag indicating if the matrix is indefinite
533 :
534 : Level: advanced
535 :
536 : .seealso: BVSetMatrix()
537 : @*/
538 1750 : PetscErrorCode BVGetMatrix(BV bv,Mat *B,PetscBool *indef)
539 : {
540 1750 : PetscFunctionBegin;
541 1750 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
542 1750 : if (B) *B = bv->matrix;
543 1750 : if (indef) *indef = bv->indef;
544 1750 : PetscFunctionReturn(PETSC_SUCCESS);
545 : }
546 :
547 : /*@
548 : BVApplyMatrix - Multiplies a vector by the matrix representation of the
549 : inner product.
550 :
551 : Neighbor-wise Collective
552 :
553 : Input Parameters:
554 : + bv - the basis vectors context
555 : - x - the vector
556 :
557 : Output Parameter:
558 : . y - the result
559 :
560 : Note:
561 : If no matrix was specified this function copies the vector.
562 :
563 : Level: advanced
564 :
565 : .seealso: BVSetMatrix(), BVApplyMatrixBV()
566 : @*/
567 595 : PetscErrorCode BVApplyMatrix(BV bv,Vec x,Vec y)
568 : {
569 595 : PetscFunctionBegin;
570 595 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
571 595 : PetscValidHeaderSpecific(x,VEC_CLASSID,2);
572 595 : PetscValidHeaderSpecific(y,VEC_CLASSID,3);
573 595 : if (bv->matrix) {
574 595 : PetscCall(BV_IPMatMult(bv,x));
575 595 : PetscCall(VecCopy(bv->Bx,y));
576 0 : } else PetscCall(VecCopy(x,y));
577 595 : PetscFunctionReturn(PETSC_SUCCESS);
578 : }
579 :
580 : /*@
581 : BVApplyMatrixBV - Multiplies the BV vectors by the matrix representation
582 : of the inner product.
583 :
584 : Neighbor-wise Collective
585 :
586 : Input Parameter:
587 : . X - the basis vectors context
588 :
589 : Output Parameter:
590 : . Y - the basis vectors to store the result (optional)
591 :
592 : Note:
593 : This function computes Y = B*X, where B is the matrix given with
594 : BVSetMatrix(). This operation is computed as in BVMatMult().
595 : If no matrix was specified, then it just copies Y = X.
596 :
597 : If no Y is given, the result is stored internally in the cached BV.
598 :
599 : Level: developer
600 :
601 : .seealso: BVSetMatrix(), BVApplyMatrix(), BVMatMult(), BVGetCachedBV()
602 : @*/
603 8 : PetscErrorCode BVApplyMatrixBV(BV X,BV Y)
604 : {
605 8 : PetscFunctionBegin;
606 8 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
607 8 : if (Y) {
608 8 : PetscValidHeaderSpecific(Y,BV_CLASSID,2);
609 8 : if (X->matrix) PetscCall(BVMatMult(X,X->matrix,Y));
610 0 : else PetscCall(BVCopy(X,Y));
611 0 : } else PetscCall(BV_IPMatMultBV(X));
612 8 : PetscFunctionReturn(PETSC_SUCCESS);
613 : }
614 :
615 : /*@
616 : BVSetSignature - Sets the signature matrix to be used in orthogonalization.
617 :
618 : Logically Collective
619 :
620 : Input Parameters:
621 : + bv - the basis vectors context
622 : - omega - a vector representing the diagonal of the signature matrix
623 :
624 : Note:
625 : The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
626 :
627 : Level: developer
628 :
629 : .seealso: BVSetMatrix(), BVGetSignature()
630 : @*/
631 3222 : PetscErrorCode BVSetSignature(BV bv,Vec omega)
632 : {
633 3222 : PetscInt i,n;
634 3222 : const PetscScalar *pomega;
635 3222 : PetscScalar *intern;
636 :
637 3222 : PetscFunctionBegin;
638 3222 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
639 3222 : BVCheckSizes(bv,1);
640 3222 : PetscValidHeaderSpecific(omega,VEC_CLASSID,2);
641 3222 : PetscValidType(omega,2);
642 :
643 3222 : PetscCall(VecGetSize(omega,&n));
644 3222 : PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
645 3222 : PetscCall(BV_AllocateSignature(bv));
646 3222 : if (bv->indef) {
647 3218 : PetscCall(VecGetArrayRead(omega,&pomega));
648 3218 : PetscCall(VecGetArray(bv->omega,&intern));
649 53316 : for (i=0;i<n;i++) intern[bv->nc+i] = pomega[i];
650 3218 : PetscCall(VecRestoreArray(bv->omega,&intern));
651 3218 : PetscCall(VecRestoreArrayRead(omega,&pomega));
652 4 : } else PetscCall(PetscInfo(bv,"Ignoring signature because BV is not indefinite\n"));
653 3222 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
654 3222 : PetscFunctionReturn(PETSC_SUCCESS);
655 : }
656 :
657 : /*@
658 : BVGetSignature - Retrieves the signature matrix from last orthogonalization.
659 :
660 : Not Collective
661 :
662 : Input Parameter:
663 : . bv - the basis vectors context
664 :
665 : Output Parameter:
666 : . omega - a vector representing the diagonal of the signature matrix
667 :
668 : Note:
669 : The signature matrix Omega = V'*B*V is relevant only for an indefinite B.
670 :
671 : Level: developer
672 :
673 : .seealso: BVSetMatrix(), BVSetSignature()
674 : @*/
675 87 : PetscErrorCode BVGetSignature(BV bv,Vec omega)
676 : {
677 87 : PetscInt i,n;
678 87 : PetscScalar *pomega;
679 87 : const PetscScalar *intern;
680 :
681 87 : PetscFunctionBegin;
682 87 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
683 87 : BVCheckSizes(bv,1);
684 87 : PetscValidHeaderSpecific(omega,VEC_CLASSID,2);
685 87 : PetscValidType(omega,2);
686 :
687 87 : PetscCall(VecGetSize(omega,&n));
688 87 : PetscCheck(n==bv->k,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Vec argument has %" PetscInt_FMT " elements, should be %" PetscInt_FMT,n,bv->k);
689 87 : if (bv->indef && bv->omega) {
690 87 : PetscCall(VecGetArray(omega,&pomega));
691 87 : PetscCall(VecGetArrayRead(bv->omega,&intern));
692 1063 : for (i=0;i<n;i++) pomega[i] = intern[bv->nc+i];
693 87 : PetscCall(VecRestoreArrayRead(bv->omega,&intern));
694 87 : PetscCall(VecRestoreArray(omega,&pomega));
695 0 : } else PetscCall(VecSet(omega,1.0));
696 87 : PetscFunctionReturn(PETSC_SUCCESS);
697 : }
698 :
699 : /*@
700 : BVSetBufferVec - Attach a vector object to be used as buffer space for
701 : several operations.
702 :
703 : Collective
704 :
705 : Input Parameters:
706 : + bv - the basis vectors context)
707 : - buffer - the vector
708 :
709 : Notes:
710 : Use BVGetBufferVec() to retrieve the vector (for example, to free it
711 : at the end of the computations).
712 :
713 : The vector must be sequential of length (nc+m)*m, where m is the number
714 : of columns of bv and nc is the number of constraints.
715 :
716 : Level: developer
717 :
718 : .seealso: BVGetBufferVec(), BVSetSizes(), BVGetNumConstraints()
719 : @*/
720 0 : PetscErrorCode BVSetBufferVec(BV bv,Vec buffer)
721 : {
722 0 : PetscInt ld,n;
723 0 : PetscMPIInt size;
724 :
725 0 : PetscFunctionBegin;
726 0 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
727 0 : PetscValidHeaderSpecific(buffer,VEC_CLASSID,2);
728 0 : BVCheckSizes(bv,1);
729 0 : PetscCall(VecGetSize(buffer,&n));
730 0 : ld = bv->m+bv->nc;
731 0 : PetscCheck(n==ld*bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_SIZ,"Buffer size must be %" PetscInt_FMT,ld*bv->m);
732 0 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)buffer),&size));
733 0 : PetscCheck(size==1,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Buffer must be a sequential vector");
734 :
735 0 : PetscCall(PetscObjectReference((PetscObject)buffer));
736 0 : PetscCall(VecDestroy(&bv->buffer));
737 0 : bv->buffer = buffer;
738 0 : PetscFunctionReturn(PETSC_SUCCESS);
739 : }
740 :
741 : /*@
742 : BVGetBufferVec - Obtain the buffer vector associated with the BV object.
743 :
744 : Collective
745 :
746 : Input Parameters:
747 : . bv - the basis vectors context
748 :
749 : Output Parameter:
750 : . buffer - vector
751 :
752 : Notes:
753 : The vector is created if not available previously. It is a sequential vector
754 : of length (nc+m)*m, where m is the number of columns of bv and nc is the number
755 : of constraints.
756 :
757 : Developer Notes:
758 : The buffer vector is viewed as a column-major matrix with leading dimension
759 : ld=nc+m, and m columns at most. In the most common usage, it has the structure
760 : .vb
761 : | | C |
762 : |s|---|
763 : | | H |
764 : .ve
765 : where H is an upper Hessenberg matrix of order m x (m-1), C contains coefficients
766 : related to orthogonalization against constraints (first nc rows), and s is the
767 : first column that contains scratch values computed during Gram-Schmidt
768 : orthogonalization. In particular, BVDotColumn() and BVMultColumn() use s to
769 : store the coefficients.
770 :
771 : Level: developer
772 :
773 : .seealso: BVSetBufferVec(), BVSetSizes(), BVGetNumConstraints(), BVDotColumn(), BVMultColumn()
774 : @*/
775 6800 : PetscErrorCode BVGetBufferVec(BV bv,Vec *buffer)
776 : {
777 6800 : PetscInt ld;
778 :
779 6800 : PetscFunctionBegin;
780 6800 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
781 6800 : PetscAssertPointer(buffer,2);
782 6800 : BVCheckSizes(bv,1);
783 6800 : if (!bv->buffer) {
784 1820 : ld = bv->m+bv->nc;
785 1820 : PetscCall(VecCreate(PETSC_COMM_SELF,&bv->buffer));
786 1820 : PetscCall(VecSetSizes(bv->buffer,PETSC_DECIDE,ld*bv->m));
787 1820 : PetscCall(VecSetType(bv->buffer,bv->vtype));
788 : }
789 6800 : *buffer = bv->buffer;
790 6800 : PetscFunctionReturn(PETSC_SUCCESS);
791 : }
792 :
793 : /*@
794 : BVSetRandomContext - Sets the PetscRandom object associated with the BV,
795 : to be used in operations that need random numbers.
796 :
797 : Collective
798 :
799 : Input Parameters:
800 : + bv - the basis vectors context
801 : - rand - the random number generator context
802 :
803 : Level: advanced
804 :
805 : .seealso: BVGetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
806 : @*/
807 0 : PetscErrorCode BVSetRandomContext(BV bv,PetscRandom rand)
808 : {
809 0 : PetscFunctionBegin;
810 0 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
811 0 : PetscValidHeaderSpecific(rand,PETSC_RANDOM_CLASSID,2);
812 0 : PetscCheckSameComm(bv,1,rand,2);
813 0 : PetscCall(PetscObjectReference((PetscObject)rand));
814 0 : PetscCall(PetscRandomDestroy(&bv->rand));
815 0 : bv->rand = rand;
816 0 : PetscFunctionReturn(PETSC_SUCCESS);
817 : }
818 :
819 : /*@
820 : BVGetRandomContext - Gets the PetscRandom object associated with the BV.
821 :
822 : Collective
823 :
824 : Input Parameter:
825 : . bv - the basis vectors context
826 :
827 : Output Parameter:
828 : . rand - the random number generator context
829 :
830 : Level: advanced
831 :
832 : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomColumn(), BVSetRandomCond()
833 : @*/
834 3079 : PetscErrorCode BVGetRandomContext(BV bv,PetscRandom* rand)
835 : {
836 3079 : PetscFunctionBegin;
837 3079 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
838 3079 : PetscAssertPointer(rand,2);
839 3079 : if (!bv->rand) {
840 934 : PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)bv),&bv->rand));
841 934 : if (bv->cuda) PetscCall(PetscRandomSetType(bv->rand,PETSCCURAND));
842 934 : if (bv->sfocalled) PetscCall(PetscRandomSetFromOptions(bv->rand));
843 934 : if (bv->rrandom) {
844 16 : PetscCall(PetscRandomSetSeed(bv->rand,0x12345678));
845 16 : PetscCall(PetscRandomSeed(bv->rand));
846 : }
847 : }
848 3079 : *rand = bv->rand;
849 3079 : PetscFunctionReturn(PETSC_SUCCESS);
850 : }
851 :
852 : /*@
853 : BVSetFromOptions - Sets BV options from the options database.
854 :
855 : Collective
856 :
857 : Input Parameter:
858 : . bv - the basis vectors context
859 :
860 : Level: beginner
861 :
862 : .seealso: BVSetOptionsPrefix()
863 : @*/
864 1936 : PetscErrorCode BVSetFromOptions(BV bv)
865 : {
866 1936 : char type[256];
867 1936 : PetscBool flg1,flg2,flg3,flg4;
868 1936 : PetscReal r;
869 1936 : BVOrthogType otype;
870 1936 : BVOrthogRefineType orefine;
871 1936 : BVOrthogBlockType oblock;
872 :
873 1936 : PetscFunctionBegin;
874 1936 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
875 1936 : PetscCall(BVRegisterAll());
876 5808 : PetscObjectOptionsBegin((PetscObject)bv);
877 1961 : PetscCall(PetscOptionsFList("-bv_type","Basis Vectors type","BVSetType",BVList,(char*)(((PetscObject)bv)->type_name?((PetscObject)bv)->type_name:BVMAT),type,sizeof(type),&flg1));
878 1936 : if (flg1) PetscCall(BVSetType(bv,type));
879 1454 : else if (!((PetscObject)bv)->type_name) PetscCall(BVSetType(bv,BVMAT));
880 :
881 1936 : otype = bv->orthog_type;
882 1936 : PetscCall(PetscOptionsEnum("-bv_orthog_type","Orthogonalization method","BVSetOrthogonalization",BVOrthogTypes,(PetscEnum)otype,(PetscEnum*)&otype,&flg1));
883 1936 : orefine = bv->orthog_ref;
884 1936 : PetscCall(PetscOptionsEnum("-bv_orthog_refine","Iterative refinement mode during orthogonalization","BVSetOrthogonalization",BVOrthogRefineTypes,(PetscEnum)orefine,(PetscEnum*)&orefine,&flg2));
885 1936 : r = bv->orthog_eta;
886 1936 : PetscCall(PetscOptionsReal("-bv_orthog_eta","Parameter of iterative refinement during orthogonalization","BVSetOrthogonalization",r,&r,&flg3));
887 1936 : oblock = bv->orthog_block;
888 1936 : PetscCall(PetscOptionsEnum("-bv_orthog_block","Block orthogonalization method","BVSetOrthogonalization",BVOrthogBlockTypes,(PetscEnum)oblock,(PetscEnum*)&oblock,&flg4));
889 1936 : if (flg1 || flg2 || flg3 || flg4) PetscCall(BVSetOrthogonalization(bv,otype,orefine,r,oblock));
890 :
891 1936 : PetscCall(PetscOptionsEnum("-bv_matmult","Method for BVMatMult","BVSetMatMultMethod",BVMatMultTypes,(PetscEnum)bv->vmm,(PetscEnum*)&bv->vmm,NULL));
892 :
893 1936 : PetscCall(PetscOptionsReal("-bv_definite_tol","Tolerance for checking a definite inner product","BVSetDefiniteTolerance",r,&r,&flg1));
894 1936 : if (flg1) PetscCall(BVSetDefiniteTolerance(bv,r));
895 :
896 : /* undocumented option to generate random vectors that are independent of the number of processes */
897 1936 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-bv_reproducible_random",&bv->rrandom,NULL));
898 :
899 1936 : if (bv->ops->create) bv->defersfo = PETSC_TRUE; /* defer call to setfromoptions */
900 592 : else PetscTryTypeMethod(bv,setfromoptions,PetscOptionsObject);
901 1936 : PetscCall(PetscObjectProcessOptionsHandlers((PetscObject)bv,PetscOptionsObject));
902 1936 : PetscOptionsEnd();
903 1936 : bv->sfocalled = PETSC_TRUE;
904 1936 : PetscFunctionReturn(PETSC_SUCCESS);
905 : }
906 :
907 : /*@
908 : BVSetOrthogonalization - Specifies the method used for the orthogonalization of
909 : vectors (classical or modified Gram-Schmidt with or without refinement), and
910 : for the block-orthogonalization (simultaneous orthogonalization of a set of
911 : vectors).
912 :
913 : Logically Collective
914 :
915 : Input Parameters:
916 : + bv - the basis vectors context
917 : . type - the method of vector orthogonalization
918 : . refine - type of refinement
919 : . eta - parameter for selective refinement
920 : - block - the method of block orthogonalization
921 :
922 : Options Database Keys:
923 : + -bv_orthog_type <type> - Where <type> is cgs for Classical Gram-Schmidt orthogonalization
924 : (default) or mgs for Modified Gram-Schmidt orthogonalization
925 : . -bv_orthog_refine <ref> - Where <ref> is one of never, ifneeded (default) or always
926 : . -bv_orthog_eta <eta> - For setting the value of eta
927 : - -bv_orthog_block <block> - Where <block> is the block-orthogonalization method
928 :
929 : Notes:
930 : The default settings work well for most problems.
931 :
932 : The parameter eta should be a real value between 0 and 1, that is used only when
933 : the refinement type is "ifneeded". Use PETSC_DETERMINE to set a reasonable
934 : default value. Use PETSC_CURRENT to leave the current value unchanged.
935 :
936 : When using several processors, MGS is likely to result in bad scalability.
937 :
938 : If the method set for block orthogonalization is GS, then the computation
939 : is done column by column with the vector orthogonalization.
940 :
941 : Level: advanced
942 :
943 : .seealso: BVOrthogonalizeColumn(), BVGetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
944 : @*/
945 448 : PetscErrorCode BVSetOrthogonalization(BV bv,BVOrthogType type,BVOrthogRefineType refine,PetscReal eta,BVOrthogBlockType block)
946 : {
947 448 : PetscFunctionBegin;
948 448 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
949 1344 : PetscValidLogicalCollectiveEnum(bv,type,2);
950 1344 : PetscValidLogicalCollectiveEnum(bv,refine,3);
951 1344 : PetscValidLogicalCollectiveReal(bv,eta,4);
952 1344 : PetscValidLogicalCollectiveEnum(bv,block,5);
953 448 : switch (type) {
954 448 : case BV_ORTHOG_CGS:
955 : case BV_ORTHOG_MGS:
956 448 : bv->orthog_type = type;
957 448 : break;
958 0 : default:
959 0 : SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown orthogonalization type");
960 : }
961 448 : switch (refine) {
962 448 : case BV_ORTHOG_REFINE_NEVER:
963 : case BV_ORTHOG_REFINE_IFNEEDED:
964 : case BV_ORTHOG_REFINE_ALWAYS:
965 448 : bv->orthog_ref = refine;
966 448 : break;
967 0 : default:
968 0 : SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown refinement type");
969 : }
970 448 : if (eta == (PetscReal)PETSC_DETERMINE) {
971 0 : bv->orthog_eta = 0.7071;
972 448 : } else if (eta != (PetscReal)PETSC_CURRENT) {
973 432 : PetscCheck(eta>0.0 && eta<=1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Invalid eta value");
974 432 : bv->orthog_eta = eta;
975 : }
976 448 : switch (block) {
977 448 : case BV_ORTHOG_BLOCK_GS:
978 : case BV_ORTHOG_BLOCK_CHOL:
979 : case BV_ORTHOG_BLOCK_TSQR:
980 : case BV_ORTHOG_BLOCK_TSQRCHOL:
981 : case BV_ORTHOG_BLOCK_SVQB:
982 448 : bv->orthog_block = block;
983 448 : break;
984 0 : default:
985 0 : SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown block orthogonalization type");
986 : }
987 448 : PetscFunctionReturn(PETSC_SUCCESS);
988 : }
989 :
990 : /*@
991 : BVGetOrthogonalization - Gets the orthogonalization settings from the BV object.
992 :
993 : Not Collective
994 :
995 : Input Parameter:
996 : . bv - basis vectors context
997 :
998 : Output Parameters:
999 : + type - the method of vector orthogonalization
1000 : . refine - type of refinement
1001 : . eta - parameter for selective refinement
1002 : - block - the method of block orthogonalization
1003 :
1004 : Level: advanced
1005 :
1006 : .seealso: BVOrthogonalizeColumn(), BVSetOrthogonalization(), BVOrthogType, BVOrthogRefineType, BVOrthogBlockType
1007 : @*/
1008 433 : PetscErrorCode BVGetOrthogonalization(BV bv,BVOrthogType *type,BVOrthogRefineType *refine,PetscReal *eta,BVOrthogBlockType *block)
1009 : {
1010 433 : PetscFunctionBegin;
1011 433 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1012 433 : if (type) *type = bv->orthog_type;
1013 433 : if (refine) *refine = bv->orthog_ref;
1014 433 : if (eta) *eta = bv->orthog_eta;
1015 433 : if (block) *block = bv->orthog_block;
1016 433 : PetscFunctionReturn(PETSC_SUCCESS);
1017 : }
1018 :
1019 : /*@
1020 : BVSetMatMultMethod - Specifies the method used for the BVMatMult() operation.
1021 :
1022 : Logically Collective
1023 :
1024 : Input Parameters:
1025 : + bv - the basis vectors context
1026 : - method - the method for the BVMatMult() operation
1027 :
1028 : Options Database Keys:
1029 : . -bv_matmult <meth> - choose one of the methods: vecs, mat
1030 :
1031 : Notes:
1032 : Allowed values are
1033 : + BV_MATMULT_VECS - perform a matrix-vector multiply per each column
1034 : . BV_MATMULT_MAT - carry out a Mat-Mat product with a dense matrix
1035 : - BV_MATMULT_MAT_SAVE - this case is deprecated
1036 :
1037 : The default is BV_MATMULT_MAT except in the case of BVVECS.
1038 :
1039 : Level: advanced
1040 :
1041 : .seealso: BVMatMult(), BVGetMatMultMethod(), BVMatMultType
1042 : @*/
1043 80 : PetscErrorCode BVSetMatMultMethod(BV bv,BVMatMultType method)
1044 : {
1045 80 : PetscFunctionBegin;
1046 80 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1047 240 : PetscValidLogicalCollectiveEnum(bv,method,2);
1048 80 : switch (method) {
1049 80 : case BV_MATMULT_VECS:
1050 : case BV_MATMULT_MAT:
1051 80 : bv->vmm = method;
1052 80 : break;
1053 0 : case BV_MATMULT_MAT_SAVE:
1054 0 : PetscCall(PetscInfo(bv,"BV_MATMULT_MAT_SAVE is deprecated, using BV_MATMULT_MAT\n"));
1055 0 : bv->vmm = BV_MATMULT_MAT;
1056 0 : break;
1057 0 : default:
1058 0 : SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Unknown matmult method");
1059 : }
1060 80 : PetscFunctionReturn(PETSC_SUCCESS);
1061 : }
1062 :
1063 : /*@
1064 : BVGetMatMultMethod - Gets the method used for the BVMatMult() operation.
1065 :
1066 : Not Collective
1067 :
1068 : Input Parameter:
1069 : . bv - basis vectors context
1070 :
1071 : Output Parameter:
1072 : . method - the method for the BVMatMult() operation
1073 :
1074 : Level: advanced
1075 :
1076 : .seealso: BVMatMult(), BVSetMatMultMethod(), BVMatMultType
1077 : @*/
1078 40 : PetscErrorCode BVGetMatMultMethod(BV bv,BVMatMultType *method)
1079 : {
1080 40 : PetscFunctionBegin;
1081 40 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1082 40 : PetscAssertPointer(method,2);
1083 40 : *method = bv->vmm;
1084 40 : PetscFunctionReturn(PETSC_SUCCESS);
1085 : }
1086 :
1087 : /*@
1088 : BVGetColumn - Returns a Vec object that contains the entries of the
1089 : requested column of the basis vectors object.
1090 :
1091 : Logically Collective
1092 :
1093 : Input Parameters:
1094 : + bv - the basis vectors context
1095 : - j - the index of the requested column
1096 :
1097 : Output Parameter:
1098 : . v - vector containing the jth column
1099 :
1100 : Notes:
1101 : The returned Vec must be seen as a reference (not a copy) of the BV
1102 : column, that is, modifying the Vec will change the BV entries as well.
1103 :
1104 : The returned Vec must not be destroyed. BVRestoreColumn() must be
1105 : called when it is no longer needed. At most, two columns can be fetched,
1106 : that is, this function can only be called twice before the corresponding
1107 : BVRestoreColumn() is invoked.
1108 :
1109 : A negative index j selects the i-th constraint, where i=-j. Constraints
1110 : should not be modified.
1111 :
1112 : Level: beginner
1113 :
1114 : .seealso: BVRestoreColumn(), BVInsertConstraints()
1115 : @*/
1116 1579480 : PetscErrorCode BVGetColumn(BV bv,PetscInt j,Vec *v)
1117 : {
1118 1579480 : PetscInt l;
1119 :
1120 1579480 : PetscFunctionBegin;
1121 1579480 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1122 1579480 : PetscValidType(bv,1);
1123 1579480 : BVCheckSizes(bv,1);
1124 1579480 : BVCheckOp(bv,1,getcolumn);
1125 4738440 : PetscValidLogicalCollectiveInt(bv,j,2);
1126 1579480 : PetscCheck(j>=0 || -j<=bv->nc,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %" PetscInt_FMT " but only %" PetscInt_FMT " are available",-j,bv->nc);
1127 1579480 : PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %" PetscInt_FMT " but only %" PetscInt_FMT " are available",j,bv->m);
1128 1579480 : PetscCheck(j!=bv->ci[0] && j!=bv->ci[1],PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Column %" PetscInt_FMT " already fetched in a previous call to BVGetColumn",j);
1129 1579480 : l = BVAvailableVec;
1130 1579480 : PetscCheck(l!=-1,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Too many requested columns; you must call BVRestoreColumn for one of the previously fetched columns");
1131 1579480 : PetscUseTypeMethod(bv,getcolumn,j,v);
1132 1579480 : bv->ci[l] = j;
1133 1579480 : PetscCall(VecGetState(bv->cv[l],&bv->st[l]));
1134 1579480 : PetscCall(PetscObjectGetId((PetscObject)bv->cv[l],&bv->id[l]));
1135 1579480 : *v = bv->cv[l];
1136 1579480 : PetscFunctionReturn(PETSC_SUCCESS);
1137 : }
1138 :
1139 : /*@
1140 : BVRestoreColumn - Restore a column obtained with BVGetColumn().
1141 :
1142 : Logically Collective
1143 :
1144 : Input Parameters:
1145 : + bv - the basis vectors context
1146 : . j - the index of the column
1147 : - v - vector obtained with BVGetColumn()
1148 :
1149 : Note:
1150 : The arguments must match the corresponding call to BVGetColumn().
1151 :
1152 : Level: beginner
1153 :
1154 : .seealso: BVGetColumn()
1155 : @*/
1156 1579480 : PetscErrorCode BVRestoreColumn(BV bv,PetscInt j,Vec *v)
1157 : {
1158 1579480 : PetscObjectId id;
1159 1579480 : PetscObjectState st;
1160 1579480 : PetscInt l;
1161 :
1162 1579480 : PetscFunctionBegin;
1163 1579480 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1164 1579480 : PetscValidType(bv,1);
1165 1579480 : BVCheckSizes(bv,1);
1166 4738440 : PetscValidLogicalCollectiveInt(bv,j,2);
1167 1579480 : PetscAssertPointer(v,3);
1168 1579480 : PetscValidHeaderSpecific(*v,VEC_CLASSID,3);
1169 1579480 : PetscCheck(j>=0 || -j<=bv->nc,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested constraint %" PetscInt_FMT " but only %" PetscInt_FMT " are available",-j,bv->nc);
1170 1579480 : PetscCheck(j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"You requested column %" PetscInt_FMT " but only %" PetscInt_FMT " are available",j,bv->m);
1171 1579480 : PetscCheck(j==bv->ci[0] || j==bv->ci[1],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Column %" PetscInt_FMT " has not been fetched with a call to BVGetColumn",j);
1172 1579480 : l = (j==bv->ci[0])? 0: 1;
1173 1579480 : PetscCall(PetscObjectGetId((PetscObject)*v,&id));
1174 1579480 : PetscCheck(id==bv->id[l],PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same Vec that was obtained with BVGetColumn");
1175 1579480 : PetscCall(VecGetState(*v,&st));
1176 1579480 : if (st!=bv->st[l]) PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1177 1579480 : PetscUseTypeMethod(bv,restorecolumn,j,v);
1178 1579480 : bv->ci[l] = -bv->nc-1;
1179 1579480 : bv->st[l] = -1;
1180 1579480 : bv->id[l] = 0;
1181 1579480 : *v = NULL;
1182 1579480 : PetscFunctionReturn(PETSC_SUCCESS);
1183 : }
1184 :
1185 : /*@C
1186 : BVGetArray - Returns a pointer to a contiguous array that contains this
1187 : processor's portion of the BV data.
1188 :
1189 : Logically Collective
1190 :
1191 : Input Parameters:
1192 : . bv - the basis vectors context
1193 :
1194 : Output Parameter:
1195 : . a - location to put pointer to the array
1196 :
1197 : Notes:
1198 : BVRestoreArray() must be called when access to the array is no longer needed.
1199 : This operation may imply a data copy, for BV types that do not store
1200 : data contiguously in memory.
1201 :
1202 : The pointer will normally point to the first entry of the first column,
1203 : but if the BV has constraints then these go before the regular columns.
1204 :
1205 : Note that for manipulating the pointer to the BV array, one must take into
1206 : account the leading dimension, which might be different from the local
1207 : number of rows, see BVGetLeadingDimension().
1208 :
1209 : Use BVGetArrayRead() for read-only access.
1210 :
1211 : Level: advanced
1212 :
1213 : .seealso: BVRestoreArray(), BVInsertConstraints(), BVGetLeadingDimension(), BVGetArrayRead()
1214 : @*/
1215 27738 : PetscErrorCode BVGetArray(BV bv,PetscScalar **a)
1216 : {
1217 27738 : PetscFunctionBegin;
1218 27738 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1219 27738 : PetscValidType(bv,1);
1220 27738 : BVCheckSizes(bv,1);
1221 27738 : BVCheckOp(bv,1,getarray);
1222 27738 : PetscUseTypeMethod(bv,getarray,a);
1223 27738 : PetscFunctionReturn(PETSC_SUCCESS);
1224 : }
1225 :
1226 : /*@C
1227 : BVRestoreArray - Restore the BV object after BVGetArray() has been called.
1228 :
1229 : Logically Collective
1230 :
1231 : Input Parameters:
1232 : + bv - the basis vectors context
1233 : - a - location of pointer to array obtained from BVGetArray()
1234 :
1235 : Note:
1236 : This operation may imply a data copy, for BV types that do not store
1237 : data contiguously in memory.
1238 :
1239 : Level: advanced
1240 :
1241 : .seealso: BVGetColumn()
1242 : @*/
1243 27738 : PetscErrorCode BVRestoreArray(BV bv,PetscScalar **a)
1244 : {
1245 27738 : PetscFunctionBegin;
1246 27738 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1247 27738 : PetscValidType(bv,1);
1248 27738 : BVCheckSizes(bv,1);
1249 27738 : PetscTryTypeMethod(bv,restorearray,a);
1250 27738 : if (a) *a = NULL;
1251 27738 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
1252 27738 : PetscFunctionReturn(PETSC_SUCCESS);
1253 : }
1254 :
1255 : /*@C
1256 : BVGetArrayRead - Returns a read-only pointer to a contiguous array that
1257 : contains this processor's portion of the BV data.
1258 :
1259 : Not Collective
1260 :
1261 : Input Parameters:
1262 : . bv - the basis vectors context
1263 :
1264 : Output Parameter:
1265 : . a - location to put pointer to the array
1266 :
1267 : Notes:
1268 : BVRestoreArrayRead() must be called when access to the array is no
1269 : longer needed. This operation may imply a data copy, for BV types that
1270 : do not store data contiguously in memory.
1271 :
1272 : The pointer will normally point to the first entry of the first column,
1273 : but if the BV has constraints then these go before the regular columns.
1274 :
1275 : Level: advanced
1276 :
1277 : .seealso: BVRestoreArray(), BVInsertConstraints(), BVGetLeadingDimension(), BVGetArray()
1278 : @*/
1279 190 : PetscErrorCode BVGetArrayRead(BV bv,const PetscScalar **a)
1280 : {
1281 190 : PetscFunctionBegin;
1282 190 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1283 190 : PetscValidType(bv,1);
1284 190 : BVCheckSizes(bv,1);
1285 190 : BVCheckOp(bv,1,getarrayread);
1286 190 : PetscUseTypeMethod(bv,getarrayread,a);
1287 190 : PetscFunctionReturn(PETSC_SUCCESS);
1288 : }
1289 :
1290 : /*@C
1291 : BVRestoreArrayRead - Restore the BV object after BVGetArrayRead() has
1292 : been called.
1293 :
1294 : Not Collective
1295 :
1296 : Input Parameters:
1297 : + bv - the basis vectors context
1298 : - a - location of pointer to array obtained from BVGetArrayRead()
1299 :
1300 : Level: advanced
1301 :
1302 : .seealso: BVGetColumn()
1303 : @*/
1304 190 : PetscErrorCode BVRestoreArrayRead(BV bv,const PetscScalar **a)
1305 : {
1306 190 : PetscFunctionBegin;
1307 190 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1308 190 : PetscValidType(bv,1);
1309 190 : BVCheckSizes(bv,1);
1310 190 : PetscTryTypeMethod(bv,restorearrayread,a);
1311 190 : if (a) *a = NULL;
1312 190 : PetscFunctionReturn(PETSC_SUCCESS);
1313 : }
1314 :
1315 : /*@
1316 : BVCreateVec - Creates a new Vec object with the same type and dimensions
1317 : as the columns of the basis vectors object.
1318 :
1319 : Collective
1320 :
1321 : Input Parameter:
1322 : . bv - the basis vectors context
1323 :
1324 : Output Parameter:
1325 : . v - the new vector
1326 :
1327 : Note:
1328 : The user is responsible of destroying the returned vector.
1329 :
1330 : Level: beginner
1331 :
1332 : .seealso: BVCreateMat(), BVCreateVecEmpty()
1333 : @*/
1334 1283 : PetscErrorCode BVCreateVec(BV bv,Vec *v)
1335 : {
1336 1283 : PetscFunctionBegin;
1337 1283 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1338 1283 : BVCheckSizes(bv,1);
1339 1283 : PetscAssertPointer(v,2);
1340 1283 : PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),v));
1341 1283 : PetscCall(VecSetLayout(*v,bv->map));
1342 1283 : PetscCall(VecSetType(*v,bv->vtype));
1343 1283 : PetscCall(VecSetUp(*v));
1344 1283 : PetscFunctionReturn(PETSC_SUCCESS);
1345 : }
1346 :
1347 : /*@
1348 : BVCreateVecEmpty - Creates a new Vec object with the same type and dimensions
1349 : as the columns of the basis vectors object, but without internal array.
1350 :
1351 : Collective
1352 :
1353 : Input Parameter:
1354 : . bv - the basis vectors context
1355 :
1356 : Output Parameter:
1357 : . v - the new vector
1358 :
1359 : Note:
1360 : This works as BVCreateVec(), but the new vector does not have the array allocated,
1361 : so the intended usage is with VecPlaceArray().
1362 :
1363 : Level: developer
1364 :
1365 : .seealso: BVCreateVec()
1366 : @*/
1367 27088 : PetscErrorCode BVCreateVecEmpty(BV bv,Vec *v)
1368 : {
1369 27088 : PetscBool standard,cuda,hip,mpi;
1370 27088 : PetscInt N,nloc,bs;
1371 :
1372 27088 : PetscFunctionBegin;
1373 27088 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1374 27088 : BVCheckSizes(bv,1);
1375 27088 : PetscAssertPointer(v,2);
1376 :
1377 27088 : PetscCall(PetscStrcmpAny(bv->vtype,&standard,VECSEQ,VECMPI,""));
1378 27088 : PetscCall(PetscStrcmpAny(bv->vtype,&cuda,VECSEQCUDA,VECMPICUDA,""));
1379 27088 : PetscCall(PetscStrcmpAny(bv->vtype,&hip,VECSEQHIP,VECMPIHIP,""));
1380 27088 : if (standard || cuda || hip) {
1381 27088 : PetscCall(PetscStrcmpAny(bv->vtype,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
1382 27088 : PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
1383 27088 : PetscCall(PetscLayoutGetSize(bv->map,&N));
1384 27088 : PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
1385 27088 : if (cuda) {
1386 : #if defined(PETSC_HAVE_CUDA)
1387 : if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1388 : else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1389 : #endif
1390 27088 : } else if (hip) {
1391 : #if defined(PETSC_HAVE_HIP)
1392 : if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1393 : else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1394 : #endif
1395 : } else {
1396 27088 : if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,N,NULL,v));
1397 22730 : else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,N,NULL,v));
1398 : }
1399 0 : } else PetscCall(BVCreateVec(bv,v)); /* standard duplicate, with internal array */
1400 27088 : PetscFunctionReturn(PETSC_SUCCESS);
1401 : }
1402 :
1403 : /*@
1404 : BVSetVecType - Set the vector type to be used when creating vectors via BVCreateVec().
1405 :
1406 : Collective
1407 :
1408 : Input Parameters:
1409 : + bv - the basis vectors context
1410 : - vtype - the vector type
1411 :
1412 : Level: advanced
1413 :
1414 : Note:
1415 : This is not needed if the BV object is set up with BVSetSizesFromVec(), but may be
1416 : required in the case of BVSetSizes() if one wants to work with non-standard vectors.
1417 :
1418 : .seealso: BVGetVecType(), BVSetSizesFromVec(), BVSetSizes()
1419 : @*/
1420 14806 : PetscErrorCode BVSetVecType(BV bv,VecType vtype)
1421 : {
1422 14806 : PetscBool std;
1423 14806 : PetscMPIInt size;
1424 :
1425 14806 : PetscFunctionBegin;
1426 14806 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1427 14806 : PetscCall(PetscFree(bv->vtype));
1428 14806 : PetscCall(PetscStrcmp(vtype,VECSTANDARD,&std));
1429 14806 : if (std) {
1430 74 : PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)bv),&size));
1431 82 : PetscCall(PetscStrallocpy((size==1)?VECSEQ:VECMPI,(char**)&bv->vtype));
1432 14732 : } else PetscCall(PetscStrallocpy(vtype,(char**)&bv->vtype));
1433 14806 : PetscFunctionReturn(PETSC_SUCCESS);
1434 : }
1435 :
1436 : /*@
1437 : BVGetVecType - Get the vector type to be used when creating vectors via BVCreateVec().
1438 :
1439 : Not Collective
1440 :
1441 : Input Parameter:
1442 : . bv - the basis vectors context
1443 :
1444 : Output Parameter:
1445 : . vtype - the vector type
1446 :
1447 : Level: advanced
1448 :
1449 : .seealso: BVSetVecType()
1450 : @*/
1451 290 : PetscErrorCode BVGetVecType(BV bv,VecType *vtype)
1452 : {
1453 290 : PetscFunctionBegin;
1454 290 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1455 290 : PetscAssertPointer(vtype,2);
1456 290 : *vtype = bv->vtype;
1457 290 : PetscFunctionReturn(PETSC_SUCCESS);
1458 : }
1459 :
1460 : /*@
1461 : BVCreateMat - Creates a new Mat object of dense type and copies the contents
1462 : of the BV object.
1463 :
1464 : Collective
1465 :
1466 : Input Parameter:
1467 : . bv - the basis vectors context
1468 :
1469 : Output Parameter:
1470 : . A - the new matrix
1471 :
1472 : Notes:
1473 : The user is responsible of destroying the returned matrix.
1474 :
1475 : The matrix contains all columns of the BV, not just the active columns.
1476 :
1477 : Level: intermediate
1478 :
1479 : .seealso: BVCreateFromMat(), BVCreateVec(), BVGetMat()
1480 : @*/
1481 39 : PetscErrorCode BVCreateMat(BV bv,Mat *A)
1482 : {
1483 39 : PetscInt ksave,lsave;
1484 39 : Mat B;
1485 :
1486 39 : PetscFunctionBegin;
1487 39 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1488 39 : BVCheckSizes(bv,1);
1489 39 : PetscAssertPointer(A,2);
1490 :
1491 39 : PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,bv->m,bv->ld,NULL,A));
1492 39 : lsave = bv->l;
1493 39 : ksave = bv->k;
1494 39 : bv->l = 0;
1495 39 : bv->k = bv->m;
1496 39 : PetscCall(BVGetMat(bv,&B));
1497 39 : PetscCall(MatCopy(B,*A,SAME_NONZERO_PATTERN));
1498 39 : PetscCall(BVRestoreMat(bv,&B));
1499 39 : bv->l = lsave;
1500 39 : bv->k = ksave;
1501 39 : PetscFunctionReturn(PETSC_SUCCESS);
1502 : }
1503 :
1504 27393 : PetscErrorCode BVGetMat_Default(BV bv,Mat *A)
1505 : {
1506 27393 : PetscScalar *vv,*aa;
1507 27393 : PetscBool create=PETSC_FALSE;
1508 27393 : PetscInt m,cols;
1509 :
1510 27393 : PetscFunctionBegin;
1511 27393 : m = bv->k-bv->l;
1512 27393 : if (!bv->Aget) create=PETSC_TRUE;
1513 : else {
1514 23588 : PetscCall(MatDenseGetArray(bv->Aget,&aa));
1515 23588 : PetscCheck(!aa,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVGetMat already called on this BV");
1516 23588 : PetscCall(MatGetSize(bv->Aget,NULL,&cols));
1517 23588 : if (cols!=m) {
1518 2498 : PetscCall(MatDestroy(&bv->Aget));
1519 : create=PETSC_TRUE;
1520 : }
1521 : }
1522 27393 : PetscCall(BVGetArray(bv,&vv));
1523 27393 : if (create) {
1524 6303 : PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)bv),bv->vtype,bv->n,PETSC_DECIDE,bv->N,m,bv->ld,vv,&bv->Aget)); /* pass a pointer to avoid allocation of storage */
1525 6303 : PetscCall(MatDenseReplaceArray(bv->Aget,NULL)); /* replace with a null pointer, the value after BVRestoreMat */
1526 : }
1527 27393 : PetscCall(MatDensePlaceArray(bv->Aget,vv+(bv->nc+bv->l)*bv->ld)); /* set the actual pointer */
1528 27393 : *A = bv->Aget;
1529 27393 : PetscFunctionReturn(PETSC_SUCCESS);
1530 : }
1531 :
1532 : /*@
1533 : BVGetMat - Returns a Mat object of dense type that shares the memory of
1534 : the BV object.
1535 :
1536 : Collective
1537 :
1538 : Input Parameter:
1539 : . bv - the basis vectors context
1540 :
1541 : Output Parameter:
1542 : . A - the matrix
1543 :
1544 : Notes:
1545 : The returned matrix contains only the active columns. If the content of
1546 : the Mat is modified, these changes are also done in the BV object. The
1547 : user must call BVRestoreMat() when no longer needed.
1548 :
1549 : This operation implies a call to BVGetArray(), which may result in data
1550 : copies.
1551 :
1552 : Level: advanced
1553 :
1554 : .seealso: BVRestoreMat(), BVCreateMat(), BVGetArray()
1555 : @*/
1556 27393 : PetscErrorCode BVGetMat(BV bv,Mat *A)
1557 : {
1558 27393 : char name[64];
1559 :
1560 27393 : PetscFunctionBegin;
1561 27393 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1562 27393 : BVCheckSizes(bv,1);
1563 27393 : PetscAssertPointer(A,2);
1564 27393 : PetscUseTypeMethod(bv,getmat,A);
1565 27393 : if (((PetscObject)bv)->name) { /* set A's name based on BV name */
1566 401 : PetscCall(PetscStrncpy(name,"Mat_",sizeof(name)));
1567 401 : PetscCall(PetscStrlcat(name,((PetscObject)bv)->name,sizeof(name)));
1568 401 : PetscCall(PetscObjectSetName((PetscObject)*A,name));
1569 : }
1570 27393 : PetscFunctionReturn(PETSC_SUCCESS);
1571 : }
1572 :
1573 27393 : PetscErrorCode BVRestoreMat_Default(BV bv,Mat *A)
1574 : {
1575 27393 : PetscScalar *vv,*aa;
1576 :
1577 27393 : PetscFunctionBegin;
1578 27393 : PetscCall(MatDenseGetArray(bv->Aget,&aa));
1579 27393 : vv = aa-(bv->nc+bv->l)*bv->ld;
1580 27393 : PetscCall(MatDenseResetArray(bv->Aget));
1581 27393 : PetscCall(BVRestoreArray(bv,&vv));
1582 27393 : *A = NULL;
1583 27393 : PetscFunctionReturn(PETSC_SUCCESS);
1584 : }
1585 :
1586 : /*@
1587 : BVRestoreMat - Restores the Mat obtained with BVGetMat().
1588 :
1589 : Logically Collective
1590 :
1591 : Input Parameters:
1592 : + bv - the basis vectors context
1593 : - A - the fetched matrix
1594 :
1595 : Note:
1596 : A call to this function must match a previous call of BVGetMat().
1597 : The effect is that the contents of the Mat are copied back to the
1598 : BV internal data structures.
1599 :
1600 : Level: advanced
1601 :
1602 : .seealso: BVGetMat(), BVRestoreArray()
1603 : @*/
1604 27393 : PetscErrorCode BVRestoreMat(BV bv,Mat *A)
1605 : {
1606 27393 : PetscFunctionBegin;
1607 27393 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1608 27393 : BVCheckSizes(bv,1);
1609 27393 : PetscAssertPointer(A,2);
1610 27393 : PetscCheck(bv->Aget,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"BVRestoreMat must match a previous call to BVGetMat");
1611 27393 : PetscCheck(bv->Aget==*A,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Mat argument is not the same as the one obtained with BVGetMat");
1612 27393 : PetscUseTypeMethod(bv,restoremat,A);
1613 27393 : PetscFunctionReturn(PETSC_SUCCESS);
1614 : }
1615 :
1616 : /*
1617 : Copy all user-provided attributes of V to another BV object W
1618 : */
1619 11876 : static inline PetscErrorCode BVDuplicate_Private(BV V,BV W)
1620 : {
1621 11876 : PetscFunctionBegin;
1622 11876 : PetscCall(PetscLayoutReference(V->map,&W->map));
1623 11876 : PetscCall(BVSetVecType(W,V->vtype));
1624 11876 : W->ld = V->ld;
1625 11876 : PetscCall(BVSetType(W,((PetscObject)V)->type_name));
1626 11876 : W->orthog_type = V->orthog_type;
1627 11876 : W->orthog_ref = V->orthog_ref;
1628 11876 : W->orthog_eta = V->orthog_eta;
1629 11876 : W->orthog_block = V->orthog_block;
1630 11876 : if (V->matrix) PetscCall(PetscObjectReference((PetscObject)V->matrix));
1631 11876 : W->matrix = V->matrix;
1632 11876 : W->indef = V->indef;
1633 11876 : W->vmm = V->vmm;
1634 11876 : W->rrandom = V->rrandom;
1635 11876 : W->deftol = V->deftol;
1636 11876 : if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
1637 11876 : W->rand = V->rand;
1638 11876 : W->sfocalled = V->sfocalled;
1639 11876 : PetscTryTypeMethod(V,duplicate,W);
1640 11876 : PetscCall(PetscObjectStateIncrease((PetscObject)W));
1641 11876 : PetscFunctionReturn(PETSC_SUCCESS);
1642 : }
1643 :
1644 : /*@
1645 : BVDuplicate - Creates a new basis vector object of the same type and
1646 : dimensions as an existing one.
1647 :
1648 : Collective
1649 :
1650 : Input Parameter:
1651 : . V - basis vectors context
1652 :
1653 : Output Parameter:
1654 : . W - location to put the new BV
1655 :
1656 : Notes:
1657 : The new BV has the same type and dimensions as V. Also, the inner
1658 : product matrix and orthogonalization options are copied.
1659 :
1660 : BVDuplicate() DOES NOT COPY the entries, but rather allocates storage
1661 : for the new basis vectors. Use BVCopy() to copy the contents.
1662 :
1663 : Level: intermediate
1664 :
1665 : .seealso: BVDuplicateResize(), BVCreate(), BVSetSizesFromVec(), BVCopy()
1666 : @*/
1667 2905 : PetscErrorCode BVDuplicate(BV V,BV *W)
1668 : {
1669 2905 : PetscFunctionBegin;
1670 2905 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
1671 2905 : PetscValidType(V,1);
1672 2905 : BVCheckSizes(V,1);
1673 2905 : PetscAssertPointer(W,2);
1674 2905 : PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1675 2905 : (*W)->N = V->N;
1676 2905 : (*W)->n = V->n;
1677 2905 : (*W)->m = V->m;
1678 2905 : (*W)->k = V->m;
1679 2905 : PetscCall(BVDuplicate_Private(V,*W));
1680 2905 : PetscFunctionReturn(PETSC_SUCCESS);
1681 : }
1682 :
1683 : /*@
1684 : BVDuplicateResize - Creates a new basis vector object of the same type and
1685 : dimensions as an existing one, but with possibly different number of columns.
1686 :
1687 : Collective
1688 :
1689 : Input Parameters:
1690 : + V - basis vectors context
1691 : - m - the new number of columns
1692 :
1693 : Output Parameter:
1694 : . W - location to put the new BV
1695 :
1696 : Note:
1697 : This is equivalent of a call to BVDuplicate() followed by BVResize(). The
1698 : contents of V are not copied to W.
1699 :
1700 : Level: intermediate
1701 :
1702 : .seealso: BVDuplicate(), BVResize()
1703 : @*/
1704 8541 : PetscErrorCode BVDuplicateResize(BV V,PetscInt m,BV *W)
1705 : {
1706 8541 : PetscFunctionBegin;
1707 8541 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
1708 8541 : PetscValidType(V,1);
1709 8541 : BVCheckSizes(V,1);
1710 25623 : PetscValidLogicalCollectiveInt(V,m,2);
1711 8541 : PetscAssertPointer(W,3);
1712 8541 : PetscCall(BVCreate(PetscObjectComm((PetscObject)V),W));
1713 8541 : (*W)->N = V->N;
1714 8541 : (*W)->n = V->n;
1715 8541 : (*W)->m = m;
1716 8541 : (*W)->k = m;
1717 8541 : PetscCall(BVDuplicate_Private(V,*W));
1718 8541 : PetscFunctionReturn(PETSC_SUCCESS);
1719 : }
1720 :
1721 : /*@
1722 : BVGetCachedBV - Returns a BV object stored internally that holds the
1723 : result of B*X after a call to BVApplyMatrixBV().
1724 :
1725 : Collective
1726 :
1727 : Input Parameter:
1728 : . bv - the basis vectors context
1729 :
1730 : Output Parameter:
1731 : . cached - the cached BV
1732 :
1733 : Note:
1734 : The cached BV is created if not available previously.
1735 :
1736 : Level: developer
1737 :
1738 : .seealso: BVApplyMatrixBV()
1739 : @*/
1740 2466 : PetscErrorCode BVGetCachedBV(BV bv,BV *cached)
1741 : {
1742 2466 : PetscFunctionBegin;
1743 2466 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1744 2466 : PetscAssertPointer(cached,2);
1745 2466 : BVCheckSizes(bv,1);
1746 2466 : if (!bv->cached) {
1747 110 : PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),&bv->cached));
1748 110 : bv->cached->N = bv->N;
1749 110 : bv->cached->n = bv->n;
1750 110 : bv->cached->m = bv->m;
1751 110 : bv->cached->k = bv->m;
1752 110 : PetscCall(BVDuplicate_Private(bv,bv->cached));
1753 : }
1754 2466 : *cached = bv->cached;
1755 2466 : PetscFunctionReturn(PETSC_SUCCESS);
1756 : }
1757 :
1758 : /*@
1759 : BVCopy - Copies a basis vector object into another one, W <- V.
1760 :
1761 : Logically Collective
1762 :
1763 : Input Parameter:
1764 : . V - basis vectors context
1765 :
1766 : Output Parameter:
1767 : . W - the copy
1768 :
1769 : Note:
1770 : Both V and W must be distributed in the same manner; local copies are
1771 : done. Only active columns (excluding the leading ones) are copied.
1772 : In the destination W, columns are overwritten starting from the leading ones.
1773 : Constraints are not copied.
1774 :
1775 : Level: beginner
1776 :
1777 : .seealso: BVCopyVec(), BVCopyColumn(), BVDuplicate(), BVSetActiveColumns()
1778 : @*/
1779 13773 : PetscErrorCode BVCopy(BV V,BV W)
1780 : {
1781 13773 : PetscScalar *womega;
1782 13773 : const PetscScalar *vomega;
1783 :
1784 13773 : PetscFunctionBegin;
1785 13773 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
1786 13773 : PetscValidType(V,1);
1787 13773 : BVCheckSizes(V,1);
1788 13773 : BVCheckOp(V,1,copy);
1789 13773 : PetscValidHeaderSpecific(W,BV_CLASSID,2);
1790 13773 : PetscValidType(W,2);
1791 13773 : BVCheckSizes(W,2);
1792 13773 : PetscCheckSameTypeAndComm(V,1,W,2);
1793 13773 : PetscCheck(V->n==W->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension V %" PetscInt_FMT ", W %" PetscInt_FMT,V->n,W->n);
1794 13773 : PetscCheck(V->k-V->l<=W->m-W->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"W has %" PetscInt_FMT " non-leading columns, not enough to store %" PetscInt_FMT " columns",W->m-W->l,V->k-V->l);
1795 13773 : if (V==W || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
1796 :
1797 13773 : PetscCall(PetscLogEventBegin(BV_Copy,V,W,0,0));
1798 13773 : if (V->indef && V->matrix && V->indef==W->indef && V->matrix==W->matrix) {
1799 : /* copy signature */
1800 8 : PetscCall(BV_AllocateSignature(W));
1801 8 : PetscCall(VecGetArrayRead(V->omega,&vomega));
1802 8 : PetscCall(VecGetArray(W->omega,&womega));
1803 8 : PetscCall(PetscArraycpy(womega+W->nc+W->l,vomega+V->nc+V->l,V->k-V->l));
1804 8 : PetscCall(VecRestoreArray(W->omega,&womega));
1805 8 : PetscCall(VecRestoreArrayRead(V->omega,&vomega));
1806 : }
1807 13773 : PetscUseTypeMethod(V,copy,W);
1808 13773 : PetscCall(PetscLogEventEnd(BV_Copy,V,W,0,0));
1809 13773 : PetscCall(PetscObjectStateIncrease((PetscObject)W));
1810 13773 : PetscFunctionReturn(PETSC_SUCCESS);
1811 : }
1812 :
1813 : /*@
1814 : BVCopyVec - Copies one of the columns of a basis vectors object into a Vec.
1815 :
1816 : Logically Collective
1817 :
1818 : Input Parameters:
1819 : + V - basis vectors context
1820 : - j - the column number to be copied
1821 :
1822 : Output Parameter:
1823 : . w - the copied column
1824 :
1825 : Note:
1826 : Both V and w must be distributed in the same manner; local copies are done.
1827 :
1828 : Level: beginner
1829 :
1830 : .seealso: BVCopy(), BVCopyColumn(), BVInsertVec()
1831 : @*/
1832 34706 : PetscErrorCode BVCopyVec(BV V,PetscInt j,Vec w)
1833 : {
1834 34706 : PetscInt n,N;
1835 34706 : Vec z;
1836 :
1837 34706 : PetscFunctionBegin;
1838 34706 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
1839 34706 : PetscValidType(V,1);
1840 34706 : BVCheckSizes(V,1);
1841 104118 : PetscValidLogicalCollectiveInt(V,j,2);
1842 34706 : PetscValidHeaderSpecific(w,VEC_CLASSID,3);
1843 34706 : PetscCheckSameComm(V,1,w,3);
1844 :
1845 34706 : PetscCall(VecGetSize(w,&N));
1846 34706 : PetscCall(VecGetLocalSize(w,&n));
1847 34706 : PetscCheck(N==V->N && n==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Vec sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ") do not match BV sizes (global %" PetscInt_FMT ", local %" PetscInt_FMT ")",N,n,V->N,V->n);
1848 :
1849 34706 : PetscCall(PetscLogEventBegin(BV_Copy,V,w,0,0));
1850 34706 : PetscCall(BVGetColumn(V,j,&z));
1851 34706 : PetscCall(VecCopy(z,w));
1852 34706 : PetscCall(BVRestoreColumn(V,j,&z));
1853 34706 : PetscCall(PetscLogEventEnd(BV_Copy,V,w,0,0));
1854 34706 : PetscFunctionReturn(PETSC_SUCCESS);
1855 : }
1856 :
1857 : /*@
1858 : BVCopyColumn - Copies the values from one of the columns to another one.
1859 :
1860 : Logically Collective
1861 :
1862 : Input Parameters:
1863 : + V - basis vectors context
1864 : . j - the number of the source column
1865 : - i - the number of the destination column
1866 :
1867 : Level: beginner
1868 :
1869 : .seealso: BVCopy(), BVCopyVec()
1870 : @*/
1871 10906 : PetscErrorCode BVCopyColumn(BV V,PetscInt j,PetscInt i)
1872 : {
1873 10906 : PetscScalar *omega;
1874 :
1875 10906 : PetscFunctionBegin;
1876 10906 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
1877 10906 : PetscValidType(V,1);
1878 10906 : BVCheckSizes(V,1);
1879 32718 : PetscValidLogicalCollectiveInt(V,j,2);
1880 32718 : PetscValidLogicalCollectiveInt(V,i,3);
1881 10906 : if (j==i) PetscFunctionReturn(PETSC_SUCCESS);
1882 :
1883 10875 : PetscCall(PetscLogEventBegin(BV_Copy,V,0,0,0));
1884 10875 : if (V->omega) {
1885 2731 : PetscCall(VecGetArray(V->omega,&omega));
1886 2731 : omega[i] = omega[j];
1887 2731 : PetscCall(VecRestoreArray(V->omega,&omega));
1888 : }
1889 10875 : PetscUseTypeMethod(V,copycolumn,j,i);
1890 10875 : PetscCall(PetscLogEventEnd(BV_Copy,V,0,0,0));
1891 10875 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
1892 10875 : PetscFunctionReturn(PETSC_SUCCESS);
1893 : }
1894 :
1895 344 : static PetscErrorCode BVGetSplit_Private(BV bv,PetscBool left,BV *split)
1896 : {
1897 344 : PetscInt ncols;
1898 :
1899 344 : PetscFunctionBegin;
1900 344 : ncols = left? bv->nc+bv->l: bv->m-bv->l;
1901 344 : if (*split && (ncols!=(*split)->m || bv->N!=(*split)->N)) PetscCall(BVDestroy(split));
1902 344 : if (!*split) {
1903 320 : PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
1904 320 : (*split)->issplit = left? 1: 2;
1905 320 : (*split)->splitparent = bv;
1906 320 : (*split)->N = bv->N;
1907 320 : (*split)->n = bv->n;
1908 320 : (*split)->m = bv->m;
1909 320 : (*split)->k = bv->m;
1910 320 : PetscCall(BVDuplicate_Private(bv,*split));
1911 : }
1912 344 : (*split)->l = 0;
1913 344 : (*split)->k = left? bv->l: bv->k-bv->l;
1914 344 : (*split)->nc = left? bv->nc: 0;
1915 344 : (*split)->m = ncols-(*split)->nc;
1916 344 : if ((*split)->nc) {
1917 24 : (*split)->ci[0] = -(*split)->nc-1;
1918 24 : (*split)->ci[1] = -(*split)->nc-1;
1919 : }
1920 344 : if (left) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
1921 172 : else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
1922 344 : PetscFunctionReturn(PETSC_SUCCESS);
1923 : }
1924 :
1925 : /*@
1926 : BVGetSplit - Splits the BV object into two BV objects that share the
1927 : internal data, one of them containing the leading columns and the other
1928 : one containing the remaining columns.
1929 :
1930 : Collective
1931 :
1932 : Input Parameter:
1933 : . bv - the basis vectors context
1934 :
1935 : Output Parameters:
1936 : + L - left BV containing leading columns (can be NULL)
1937 : - R - right BV containing remaining columns (can be NULL)
1938 :
1939 : Notes:
1940 : The columns are split in two sets. The leading columns (including the
1941 : constraints) are assigned to the left BV and the remaining columns
1942 : are assigned to the right BV. The number of leading columns, as
1943 : specified with BVSetActiveColumns(), must be between 1 and m-1 (to
1944 : guarantee that both L and R have at least one column).
1945 :
1946 : The returned BV's must be seen as references (not copies) of the input
1947 : BV, that is, modifying them will change the entries of bv as well.
1948 : The returned BV's must not be destroyed. BVRestoreSplit() must be called
1949 : when they are no longer needed.
1950 :
1951 : Pass NULL for any of the output BV's that is not needed.
1952 :
1953 : Level: advanced
1954 :
1955 : .seealso: BVRestoreSplit(), BVSetActiveColumns(), BVSetNumConstraints(), BVGetSplitRows()
1956 : @*/
1957 172 : PetscErrorCode BVGetSplit(BV bv,BV *L,BV *R)
1958 : {
1959 172 : PetscFunctionBegin;
1960 172 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1961 172 : PetscValidType(bv,1);
1962 172 : BVCheckSizes(bv,1);
1963 172 : PetscCheck(bv->l,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must indicate the number of leading columns with BVSetActiveColumns()");
1964 172 : PetscCheck(bv->lsplit>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplit() after BVGetSplitRows()");
1965 172 : PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplit()");
1966 172 : bv->lsplit = bv->nc+bv->l;
1967 172 : PetscCall(BVGetSplit_Private(bv,PETSC_TRUE,&bv->L));
1968 172 : PetscCall(BVGetSplit_Private(bv,PETSC_FALSE,&bv->R));
1969 172 : if (L) *L = bv->L;
1970 172 : if (R) *R = bv->R;
1971 172 : PetscFunctionReturn(PETSC_SUCCESS);
1972 : }
1973 :
1974 : /*@
1975 : BVRestoreSplit - Restore the BV objects obtained with BVGetSplit().
1976 :
1977 : Logically Collective
1978 :
1979 : Input Parameters:
1980 : + bv - the basis vectors context
1981 : . L - left BV obtained with BVGetSplit()
1982 : - R - right BV obtained with BVGetSplit()
1983 :
1984 : Note:
1985 : The arguments must match the corresponding call to BVGetSplit().
1986 :
1987 : Level: advanced
1988 :
1989 : .seealso: BVGetSplit()
1990 : @*/
1991 172 : PetscErrorCode BVRestoreSplit(BV bv,BV *L,BV *R)
1992 : {
1993 172 : PetscFunctionBegin;
1994 172 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
1995 172 : PetscValidType(bv,1);
1996 172 : BVCheckSizes(bv,1);
1997 172 : PetscCheck(bv->lsplit>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplit first");
1998 172 : PetscCheck(!L || *L==bv->L,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 2 is not the same BV that was obtained with BVGetSplit");
1999 172 : PetscCheck(!R || *R==bv->R,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Argument 3 is not the same BV that was obtained with BVGetSplit");
2000 172 : PetscCheck(!L || ((*L)->ci[0]<=(*L)->nc-1 && (*L)->ci[1]<=(*L)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 2 has unrestored columns, use BVRestoreColumn()");
2001 172 : PetscCheck(!R || ((*R)->ci[0]<=(*R)->nc-1 && (*R)->ci[1]<=(*R)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Argument 3 has unrestored columns, use BVRestoreColumn()");
2002 :
2003 172 : PetscTryTypeMethod(bv,restoresplit,L,R);
2004 172 : bv->lsplit = 0;
2005 172 : if (L) *L = NULL;
2006 172 : if (R) *R = NULL;
2007 172 : PetscFunctionReturn(PETSC_SUCCESS);
2008 : }
2009 :
2010 : /*
2011 : Copy all user-provided attributes of V to another BV object W with different layout
2012 : */
2013 70 : static inline PetscErrorCode BVDuplicateNewLayout_Private(BV V,BV W)
2014 : {
2015 70 : PetscFunctionBegin;
2016 70 : PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)V),W->n,W->N,1,&W->map));
2017 70 : PetscCall(BVSetVecType(W,V->vtype));
2018 70 : W->ld = V->ld;
2019 70 : PetscCall(BVSetType(W,((PetscObject)V)->type_name));
2020 70 : W->orthog_type = V->orthog_type;
2021 70 : W->orthog_ref = V->orthog_ref;
2022 70 : W->orthog_eta = V->orthog_eta;
2023 70 : W->orthog_block = V->orthog_block;
2024 70 : W->vmm = V->vmm;
2025 70 : W->rrandom = V->rrandom;
2026 70 : W->deftol = V->deftol;
2027 70 : if (V->rand) PetscCall(PetscObjectReference((PetscObject)V->rand));
2028 70 : W->rand = V->rand;
2029 70 : W->sfocalled = V->sfocalled;
2030 70 : PetscTryTypeMethod(V,duplicate,W);
2031 70 : PetscCall(PetscObjectStateIncrease((PetscObject)W));
2032 70 : PetscFunctionReturn(PETSC_SUCCESS);
2033 : }
2034 :
2035 110 : static PetscErrorCode BVGetSplitRows_Private(BV bv,PetscBool top,IS is,BV *split)
2036 : {
2037 110 : PetscInt rstart,rend,lstart;
2038 110 : PetscInt N,n;
2039 110 : PetscBool contig;
2040 :
2041 110 : PetscFunctionBegin;
2042 110 : PetscCall(PetscLayoutGetRange(bv->map,&rstart,&rend));
2043 110 : PetscCall(ISContiguousLocal(is,rstart,rend,&lstart,&contig));
2044 110 : PetscCheck(contig,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"%s index set is not contiguous",(top==PETSC_TRUE)?"Upper":"Lower");
2045 110 : if (top) PetscCheck(lstart==0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONG,"Upper index set should start at first local row");
2046 55 : else bv->lsplit = -lstart;
2047 110 : PetscCall(ISGetSize(is,&N));
2048 110 : PetscCall(ISGetLocalSize(is,&n));
2049 110 : if (*split && (bv->m!=(*split)->m || N!=(*split)->N)) PetscCall(BVDestroy(split));
2050 110 : if (!*split) {
2051 70 : PetscCall(BVCreate(PetscObjectComm((PetscObject)bv),split));
2052 70 : (*split)->issplit = top? -1: -2;
2053 70 : (*split)->splitparent = bv;
2054 70 : (*split)->N = N;
2055 70 : (*split)->n = n;
2056 70 : (*split)->m = bv->m;
2057 70 : PetscCall(BVDuplicateNewLayout_Private(bv,*split));
2058 : }
2059 110 : (*split)->k = bv->k;
2060 110 : (*split)->l = bv->l;
2061 110 : (*split)->nc = bv->nc;
2062 110 : if ((*split)->nc) {
2063 12 : (*split)->ci[0] = -(*split)->nc-1;
2064 12 : (*split)->ci[1] = -(*split)->nc-1;
2065 : }
2066 110 : if (top) PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->rstate));
2067 55 : else PetscCall(PetscObjectStateGet((PetscObject)*split,&bv->lstate));
2068 110 : PetscFunctionReturn(PETSC_SUCCESS);
2069 : }
2070 :
2071 : /*@
2072 : BVGetSplitRows - Splits the BV object into two BV objects that share the
2073 : internal data, using a disjoint horizontal splitting.
2074 :
2075 : Collective
2076 :
2077 : Input Parameters:
2078 : + bv - the basis vectors context
2079 : . isup - the index set that defines the upper part of the horizontal splitting
2080 : - islo - the index set that defines the lower part of the horizontal splitting
2081 :
2082 : Output Parameters:
2083 : + U - the resulting BV containing the upper rows
2084 : - L - the resulting BV containing the lower rows
2085 :
2086 : Notes:
2087 : The index sets must be such that every MPI process can extract the selected
2088 : rows from its local part of the input BV, and this part must be contiguous.
2089 : With one process, isup will list contiguous indices starting from 0, and islo
2090 : will contain the remaining indices, hence we refer to upper and lower part.
2091 : However, with several processes the indices will be interleaved because
2092 : isup will refer to the upper part of the local array.
2093 :
2094 : The intended use of this function is with matrices of MATNEST type, where
2095 : MatNestGetISs() will return the appropriate index sets.
2096 :
2097 : The returned BV's must be seen as references (not copies) of the input
2098 : BV, that is, modifying them will change the entries of bv as well.
2099 : The returned BV's must not be destroyed. BVRestoreSplitRows() must be called
2100 : when they are no longer needed.
2101 :
2102 : Pass NULL for any of the output BV's that is not needed.
2103 :
2104 : Level: advanced
2105 :
2106 : .seealso: BVRestoreSplitRows(), BVGetSplit()
2107 : @*/
2108 55 : PetscErrorCode BVGetSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2109 : {
2110 55 : PetscFunctionBegin;
2111 55 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
2112 55 : PetscValidType(bv,1);
2113 55 : BVCheckSizes(bv,1);
2114 55 : if (U) PetscValidHeaderSpecific(isup,IS_CLASSID,2);
2115 55 : if (L) PetscValidHeaderSpecific(islo,IS_CLASSID,3);
2116 55 : PetscCheck(bv->lsplit<=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot call BVGetSplitRows() after BVGetSplit()");
2117 55 : PetscCheck(!bv->lsplit,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Cannot get the split BV's twice before restoring them with BVRestoreSplitRows()");
2118 55 : if (U) {
2119 55 : PetscCall(BVGetSplitRows_Private(bv,PETSC_TRUE,isup,&bv->R));
2120 55 : *U = bv->R;
2121 : }
2122 55 : if (L) {
2123 55 : PetscCall(BVGetSplitRows_Private(bv,PETSC_FALSE,islo,&bv->L));
2124 55 : *L = bv->L;
2125 : }
2126 55 : PetscFunctionReturn(PETSC_SUCCESS);
2127 : }
2128 :
2129 : /*@
2130 : BVRestoreSplitRows - Restore the BV objects obtained with BVGetSplitRows().
2131 :
2132 : Logically Collective
2133 :
2134 : Input Parameters:
2135 : + bv - the basis vectors context
2136 : . isup - the index set that defines the upper part of the horizontal splitting
2137 : . islo - the index set that defines the lower part of the horizontal splitting
2138 : . U - the BV containing the upper rows
2139 : - L - the BV containing the lower rows
2140 :
2141 : Note:
2142 : The arguments must match the corresponding call to BVGetSplitRows().
2143 :
2144 : Level: advanced
2145 :
2146 : .seealso: BVGetSplitRows()
2147 : @*/
2148 55 : PetscErrorCode BVRestoreSplitRows(BV bv,IS isup,IS islo,BV *U,BV *L)
2149 : {
2150 55 : PetscFunctionBegin;
2151 55 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
2152 55 : PetscValidType(bv,1);
2153 55 : BVCheckSizes(bv,1);
2154 55 : PetscCheck(bv->lsplit<0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"Must call BVGetSplitRows first");
2155 55 : PetscCheck(!U || ((*U)->ci[0]<=(*U)->nc-1 && (*U)->ci[1]<=(*U)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The upper BV has unrestored columns, use BVRestoreColumn()");
2156 55 : PetscCheck(!L || ((*L)->ci[0]<=(*L)->nc-1 && (*L)->ci[1]<=(*L)->nc-1),PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The lower BV has unrestored columns, use BVRestoreColumn()");
2157 55 : PetscTryTypeMethod(bv,restoresplitrows,isup,islo,U,L);
2158 55 : bv->lsplit = 0;
2159 55 : PetscFunctionReturn(PETSC_SUCCESS);
2160 : }
2161 :
2162 : /*@
2163 : BVSetDefiniteTolerance - Set the tolerance to be used when checking a
2164 : definite inner product.
2165 :
2166 : Logically Collective
2167 :
2168 : Input Parameters:
2169 : + bv - basis vectors
2170 : - deftol - the tolerance
2171 :
2172 : Options Database Key:
2173 : . -bv_definite_tol <deftol> - the tolerance
2174 :
2175 : Notes:
2176 : When using a non-standard inner product, see BVSetMatrix(), the solver needs
2177 : to compute sqrt(z'*B*z) for various vectors z. If the inner product has not
2178 : been declared indefinite, the value z'*B*z must be positive, but due to
2179 : rounding error a tiny value may become negative. A tolerance is used to
2180 : detect this situation. Likewise, in complex arithmetic z'*B*z should be
2181 : real, and we use the same tolerance to check whether a nonzero imaginary part
2182 : can be considered negligible.
2183 :
2184 : This function sets this tolerance, which defaults to 10*eps, where eps is
2185 : the machine epsilon. The default value should be good for most applications.
2186 :
2187 : Level: advanced
2188 :
2189 : .seealso: BVSetMatrix()
2190 : @*/
2191 5 : PetscErrorCode BVSetDefiniteTolerance(BV bv,PetscReal deftol)
2192 : {
2193 5 : PetscFunctionBegin;
2194 5 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
2195 15 : PetscValidLogicalCollectiveReal(bv,deftol,2);
2196 5 : if (deftol == (PetscReal)PETSC_DEFAULT || deftol == (PetscReal)PETSC_DECIDE) bv->deftol = 10*PETSC_MACHINE_EPSILON;
2197 : else {
2198 5 : PetscCheck(deftol>0.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of deftol. Must be > 0");
2199 5 : bv->deftol = deftol;
2200 : }
2201 5 : PetscFunctionReturn(PETSC_SUCCESS);
2202 : }
2203 :
2204 : /*@
2205 : BVGetDefiniteTolerance - Returns the tolerance for checking a definite
2206 : inner product.
2207 :
2208 : Not Collective
2209 :
2210 : Input Parameter:
2211 : . bv - the basis vectors
2212 :
2213 : Output Parameter:
2214 : . deftol - the tolerance
2215 :
2216 : Level: advanced
2217 :
2218 : .seealso: BVSetDefiniteTolerance()
2219 : @*/
2220 0 : PetscErrorCode BVGetDefiniteTolerance(BV bv,PetscReal *deftol)
2221 : {
2222 0 : PetscFunctionBegin;
2223 0 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
2224 0 : PetscAssertPointer(deftol,2);
2225 0 : *deftol = bv->deftol;
2226 0 : PetscFunctionReturn(PETSC_SUCCESS);
2227 : }
2228 :
2229 : /*@
2230 : BVSetLeadingDimension - Set the leading dimension to be used for storing the BV data.
2231 :
2232 : Not Collective
2233 :
2234 : Input Parameters:
2235 : + bv - basis vectors
2236 : - ld - the leading dimension
2237 :
2238 : Notes:
2239 : This parameter is relevant for BVMAT, though it might be employed in other types
2240 : as well.
2241 :
2242 : When the internal data of the BV is stored as a dense matrix, the leading dimension
2243 : has the same meaning as in MatDenseSetLDA(), i.e., the distance in number of
2244 : elements from one entry of the matrix to the one in the next column at the same
2245 : row. The leading dimension refers to the local array, and hence can be different
2246 : in different processes.
2247 :
2248 : The user does not need to change this parameter. The default value is equal to the
2249 : number of local rows, but this value may be increased a little to guarantee alignment
2250 : (especially in the case of GPU storage).
2251 :
2252 : Level: advanced
2253 :
2254 : .seealso: BVGetLeadingDimension()
2255 : @*/
2256 0 : PetscErrorCode BVSetLeadingDimension(BV bv,PetscInt ld)
2257 : {
2258 0 : PetscFunctionBegin;
2259 0 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
2260 0 : PetscValidLogicalCollectiveInt(bv,ld,2);
2261 0 : PetscCheck((bv->n<0 && bv->N<0) || !bv->ops->create,PetscObjectComm((PetscObject)bv),PETSC_ERR_ORDER,"Must call BVSetLeadingDimension() before setting the BV type and sizes");
2262 0 : if (ld == PETSC_DEFAULT || ld == PETSC_DECIDE) bv->ld = 0;
2263 : else {
2264 0 : PetscCheck(ld>0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Illegal value of ld. Must be > 0");
2265 0 : bv->ld = ld;
2266 : }
2267 0 : PetscFunctionReturn(PETSC_SUCCESS);
2268 : }
2269 :
2270 : /*@
2271 : BVGetLeadingDimension - Returns the leading dimension of the BV.
2272 :
2273 : Not Collective
2274 :
2275 : Input Parameter:
2276 : . bv - the basis vectors
2277 :
2278 : Output Parameter:
2279 : . ld - the leading dimension
2280 :
2281 : Level: advanced
2282 :
2283 : Notes:
2284 : The returned value may be different in different processes.
2285 :
2286 : The leading dimension must be used when accessing the internal array via
2287 : BVGetArray() or BVGetArrayRead().
2288 :
2289 : .seealso: BVSetLeadingDimension(), BVGetArray(), BVGetArrayRead()
2290 : @*/
2291 222 : PetscErrorCode BVGetLeadingDimension(BV bv,PetscInt *ld)
2292 : {
2293 222 : PetscFunctionBegin;
2294 222 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
2295 222 : PetscAssertPointer(ld,2);
2296 222 : *ld = bv->ld;
2297 222 : PetscFunctionReturn(PETSC_SUCCESS);
2298 : }
|
