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