Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 : /*
11 : BV operations, except those involving global communication
12 : */
13 :
14 : #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
15 : #include <slepcds.h>
16 :
17 : /*@
18 : BVMult - Computes Y = beta*Y + alpha*X*Q.
19 :
20 : Logically Collective
21 :
22 : Input Parameters:
23 : + Y - first basis vectors context (modified on output)
24 : . alpha - first scalar
25 : . beta - second scalar
26 : . X - second basis vectors context
27 : - Q - (optional) sequential dense matrix
28 :
29 : Notes:
30 : X and Y must be different objects. The case X=Y can be addressed with
31 : BVMultInPlace().
32 :
33 : If matrix Q is NULL, then an AXPY operation Y = beta*Y + alpha*X is done
34 : (i.e. results as if Q = identity). If provided,
35 : the matrix Q must be a sequential dense Mat, with all entries equal on
36 : all processes (otherwise each process will compute a different update).
37 : The dimensions of Q must be at least m,n where m is the number of active
38 : columns of X and n is the number of active columns of Y.
39 :
40 : The leading columns of Y are not modified. Also, if X has leading
41 : columns specified, then these columns do not participate in the computation.
42 : Hence, only rows (resp. columns) of Q starting from lx (resp. ly) are used,
43 : where lx (resp. ly) is the number of leading columns of X (resp. Y).
44 :
45 : Level: intermediate
46 :
47 : .seealso: BVMultVec(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
48 : @*/
49 6863 : PetscErrorCode BVMult(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
50 : {
51 6863 : PetscInt m,n;
52 :
53 6863 : PetscFunctionBegin;
54 6863 : PetscValidHeaderSpecific(Y,BV_CLASSID,1);
55 27452 : PetscValidLogicalCollectiveScalar(Y,alpha,2);
56 27452 : PetscValidLogicalCollectiveScalar(Y,beta,3);
57 6863 : PetscValidHeaderSpecific(X,BV_CLASSID,4);
58 6863 : if (Q) PetscValidHeaderSpecific(Q,MAT_CLASSID,5);
59 6863 : PetscValidType(Y,1);
60 6863 : BVCheckSizes(Y,1);
61 6863 : BVCheckOp(Y,1,mult);
62 6863 : PetscValidType(X,4);
63 6863 : BVCheckSizes(X,4);
64 6863 : if (Q) PetscValidType(Q,5);
65 6863 : PetscCheckSameTypeAndComm(Y,1,X,4);
66 6863 : PetscCheck(X!=Y,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_WRONG,"X and Y arguments must be different");
67 6863 : if (Q) {
68 5385 : PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
69 5385 : PetscCall(MatGetSize(Q,&m,&n));
70 5385 : PetscCheck(m>=X->k,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,X->k);
71 5385 : PetscCheck(n>=Y->k,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,Y->k);
72 : }
73 6863 : PetscCheck(X->n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local dimension X %" PetscInt_FMT ", Y %" PetscInt_FMT,X->n,Y->n);
74 :
75 6863 : PetscCall(PetscLogEventBegin(BV_Mult,X,Y,0,0));
76 6863 : PetscUseTypeMethod(Y,mult,alpha,beta,X,Q);
77 6863 : PetscCall(PetscLogEventEnd(BV_Mult,X,Y,0,0));
78 6863 : PetscCall(PetscObjectStateIncrease((PetscObject)Y));
79 6863 : PetscFunctionReturn(PETSC_SUCCESS);
80 : }
81 :
82 : /*@
83 : BVMultVec - Computes y = beta*y + alpha*X*q.
84 :
85 : Logically Collective
86 :
87 : Input Parameters:
88 : + X - a basis vectors object
89 : . alpha - first scalar
90 : . beta - second scalar
91 : . y - a vector (modified on output)
92 : - q - an array of scalars
93 :
94 : Notes:
95 : This operation is the analogue of BVMult() but with a BV and a Vec,
96 : instead of two BV. Note that arguments are listed in different order
97 : with respect to BVMult().
98 :
99 : If X has leading columns specified, then these columns do not participate
100 : in the computation.
101 :
102 : The length of array q must be equal to the number of active columns of X
103 : minus the number of leading columns, i.e. the first entry of q multiplies
104 : the first non-leading column.
105 :
106 : Level: intermediate
107 :
108 : .seealso: BVMult(), BVMultColumn(), BVMultInPlace(), BVSetActiveColumns()
109 : @*/
110 469770 : PetscErrorCode BVMultVec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar q[])
111 : {
112 469770 : PetscInt n,N;
113 :
114 469770 : PetscFunctionBegin;
115 469770 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
116 1879080 : PetscValidLogicalCollectiveScalar(X,alpha,2);
117 1879080 : PetscValidLogicalCollectiveScalar(X,beta,3);
118 469770 : PetscValidHeaderSpecific(y,VEC_CLASSID,4);
119 469770 : PetscAssertPointer(q,5);
120 469770 : PetscValidType(X,1);
121 469770 : BVCheckSizes(X,1);
122 469770 : BVCheckOp(X,1,multvec);
123 469770 : PetscValidType(y,4);
124 469770 : PetscCheckSameComm(X,1,y,4);
125 :
126 469770 : PetscCall(VecGetSize(y,&N));
127 469770 : PetscCall(VecGetLocalSize(y,&n));
128 469770 : PetscCheck(N==X->N && n==X->n,PetscObjectComm((PetscObject)X),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,X->N,X->n);
129 :
130 469770 : PetscCall(PetscLogEventBegin(BV_MultVec,X,y,0,0));
131 469770 : PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
132 469770 : PetscCall(PetscLogEventEnd(BV_MultVec,X,y,0,0));
133 469770 : PetscFunctionReturn(PETSC_SUCCESS);
134 : }
135 :
136 : /*@
137 : BVMultColumn - Computes y = beta*y + alpha*X*q, where y is the j-th column
138 : of X.
139 :
140 : Logically Collective
141 :
142 : Input Parameters:
143 : + X - a basis vectors object
144 : . alpha - first scalar
145 : . beta - second scalar
146 : . j - the column index
147 : - q - an array of scalars
148 :
149 : Notes:
150 : This operation is equivalent to BVMultVec() but it uses column j of X
151 : rather than taking a Vec as an argument. The number of active columns of
152 : X is set to j before the computation, and restored afterwards.
153 : If X has leading columns specified, then these columns do not participate
154 : in the computation. Therefore, the length of array q must be equal to j
155 : minus the number of leading columns.
156 :
157 : Developer Notes:
158 : If q is NULL, then the coefficients are taken from position nc+l of the
159 : internal buffer vector, see BVGetBufferVec().
160 :
161 : Level: advanced
162 :
163 : .seealso: BVMult(), BVMultVec(), BVMultInPlace(), BVSetActiveColumns()
164 : @*/
165 254963 : PetscErrorCode BVMultColumn(BV X,PetscScalar alpha,PetscScalar beta,PetscInt j,PetscScalar *q)
166 : {
167 254963 : PetscInt ksave;
168 254963 : Vec y;
169 :
170 254963 : PetscFunctionBegin;
171 254963 : PetscValidHeaderSpecific(X,BV_CLASSID,1);
172 1019852 : PetscValidLogicalCollectiveScalar(X,alpha,2);
173 1019852 : PetscValidLogicalCollectiveScalar(X,beta,3);
174 1019852 : PetscValidLogicalCollectiveInt(X,j,4);
175 254963 : PetscValidType(X,1);
176 254963 : BVCheckSizes(X,1);
177 :
178 254963 : PetscCheck(j>=0,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
179 254963 : PetscCheck(j<X->m,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_OUTOFRANGE,"Index j=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j,X->m);
180 :
181 254963 : PetscCall(PetscLogEventBegin(BV_MultVec,X,0,0,0));
182 254963 : ksave = X->k;
183 254963 : X->k = j;
184 254963 : if (!q && !X->buffer) PetscCall(BVGetBufferVec(X,&X->buffer));
185 254963 : PetscCall(BVGetColumn(X,j,&y));
186 254963 : PetscUseTypeMethod(X,multvec,alpha,beta,y,q);
187 254963 : PetscCall(BVRestoreColumn(X,j,&y));
188 254963 : X->k = ksave;
189 254963 : PetscCall(PetscLogEventEnd(BV_MultVec,X,0,0,0));
190 254963 : PetscCall(PetscObjectStateIncrease((PetscObject)X));
191 254963 : PetscFunctionReturn(PETSC_SUCCESS);
192 : }
193 :
194 : /*@
195 : BVMultInPlace - Update a set of vectors as V(:,s:e-1) = V*Q(:,s:e-1).
196 :
197 : Logically Collective
198 :
199 : Input Parameters:
200 : + Q - a sequential dense matrix
201 : . s - first column of V to be overwritten
202 : - e - first column of V not to be overwritten
203 :
204 : Input/Output Parameter:
205 : . V - basis vectors
206 :
207 : Notes:
208 : The matrix Q must be a sequential dense Mat, with all entries equal on
209 : all processes (otherwise each process will compute a different update).
210 :
211 : This function computes V(:,s:e-1) = V*Q(:,s:e-1), that is, given a set of
212 : vectors V, columns from s to e-1 are overwritten with columns from s to
213 : e-1 of the matrix-matrix product V*Q. Only columns s to e-1 of Q are
214 : referenced.
215 :
216 : Level: intermediate
217 :
218 : .seealso: BVMult(), BVMultVec(), BVMultInPlaceHermitianTranspose(), BVSetActiveColumns()
219 : @*/
220 25272 : PetscErrorCode BVMultInPlace(BV V,Mat Q,PetscInt s,PetscInt e)
221 : {
222 25272 : PetscInt m,n;
223 :
224 25272 : PetscFunctionBegin;
225 25272 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
226 25272 : PetscValidHeaderSpecific(Q,MAT_CLASSID,2);
227 101088 : PetscValidLogicalCollectiveInt(V,s,3);
228 101088 : PetscValidLogicalCollectiveInt(V,e,4);
229 25272 : PetscValidType(V,1);
230 25272 : BVCheckSizes(V,1);
231 25272 : PetscValidType(Q,2);
232 25272 : PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
233 :
234 25272 : PetscCheck(s>=V->l && s<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument s has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,s,V->l,V->m);
235 25272 : PetscCheck(e>=V->l && e<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument e has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,e,V->l,V->m);
236 25272 : PetscCall(MatGetSize(Q,&m,&n));
237 25272 : PetscCheck(m>=V->k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " rows, should have at least %" PetscInt_FMT,m,V->k);
238 25272 : PetscCheck(e<=n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument only has %" PetscInt_FMT " columns, the requested value of e is larger: %" PetscInt_FMT,n,e);
239 :
240 25272 : PetscCall(PetscLogEventBegin(BV_MultInPlace,V,Q,0,0));
241 25272 : PetscUseTypeMethod(V,multinplace,Q,s,e);
242 25272 : PetscCall(PetscLogEventEnd(BV_MultInPlace,V,Q,0,0));
243 25272 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
244 25272 : PetscFunctionReturn(PETSC_SUCCESS);
245 : }
246 :
247 : /*@
248 : BVMultInPlaceHermitianTranspose - Update a set of vectors as V(:,s:e-1) = V*Q'(:,s:e-1).
249 :
250 : Logically Collective
251 :
252 : Input Parameters:
253 : + Q - a sequential dense matrix
254 : . s - first column of V to be overwritten
255 : - e - first column of V not to be overwritten
256 :
257 : Input/Output Parameter:
258 : . V - basis vectors
259 :
260 : Notes:
261 : This is a variant of BVMultInPlace() where the conjugate transpose
262 : of Q is used.
263 :
264 : Level: intermediate
265 :
266 : .seealso: BVMultInPlace()
267 : @*/
268 4 : PetscErrorCode BVMultInPlaceHermitianTranspose(BV V,Mat Q,PetscInt s,PetscInt e)
269 : {
270 4 : PetscInt m,n;
271 :
272 4 : PetscFunctionBegin;
273 4 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
274 4 : PetscValidHeaderSpecific(Q,MAT_CLASSID,2);
275 16 : PetscValidLogicalCollectiveInt(V,s,3);
276 16 : PetscValidLogicalCollectiveInt(V,e,4);
277 4 : PetscValidType(V,1);
278 4 : BVCheckSizes(V,1);
279 4 : PetscValidType(Q,2);
280 4 : PetscCheckTypeNames(Q,MATSEQDENSE,MATSEQDENSECUDA);
281 :
282 4 : PetscCheck(s>=V->l && s<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument s has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,s,V->l,V->m);
283 4 : PetscCheck(e>=V->l && e<=V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Argument e has wrong value %" PetscInt_FMT ", should be between %" PetscInt_FMT " and %" PetscInt_FMT,e,V->l,V->m);
284 4 : PetscCall(MatGetSize(Q,&m,&n));
285 4 : PetscCheck(n>=V->k,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument has %" PetscInt_FMT " columns, should have at least %" PetscInt_FMT,n,V->k);
286 4 : PetscCheck(e<=m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Mat argument only has %" PetscInt_FMT " rows, the requested value of e is larger: %" PetscInt_FMT,m,e);
287 :
288 4 : PetscCall(PetscLogEventBegin(BV_MultInPlace,V,Q,0,0));
289 4 : PetscUseTypeMethod(V,multinplacetrans,Q,s,e);
290 4 : PetscCall(PetscLogEventEnd(BV_MultInPlace,V,Q,0,0));
291 4 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
292 4 : PetscFunctionReturn(PETSC_SUCCESS);
293 : }
294 :
295 : /*@
296 : BVScale - Multiply the BV entries by a scalar value.
297 :
298 : Logically Collective
299 :
300 : Input Parameters:
301 : + bv - basis vectors
302 : - alpha - scaling factor
303 :
304 : Note:
305 : All active columns (except the leading ones) are scaled.
306 :
307 : Level: intermediate
308 :
309 : .seealso: BVScaleColumn(), BVSetActiveColumns()
310 : @*/
311 1632 : PetscErrorCode BVScale(BV bv,PetscScalar alpha)
312 : {
313 1632 : PetscFunctionBegin;
314 1632 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
315 6528 : PetscValidLogicalCollectiveScalar(bv,alpha,2);
316 1632 : PetscValidType(bv,1);
317 1632 : BVCheckSizes(bv,1);
318 1632 : if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
319 :
320 367 : PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
321 367 : PetscUseTypeMethod(bv,scale,-1,alpha);
322 367 : PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
323 367 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
324 367 : PetscFunctionReturn(PETSC_SUCCESS);
325 : }
326 :
327 : /*@
328 : BVScaleColumn - Scale one column of a BV.
329 :
330 : Logically Collective
331 :
332 : Input Parameters:
333 : + bv - basis vectors
334 : . j - column number to be scaled
335 : - alpha - scaling factor
336 :
337 : Level: intermediate
338 :
339 : .seealso: BVScale(), BVSetActiveColumns()
340 : @*/
341 97458 : PetscErrorCode BVScaleColumn(BV bv,PetscInt j,PetscScalar alpha)
342 : {
343 97458 : PetscFunctionBegin;
344 97458 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
345 389832 : PetscValidLogicalCollectiveInt(bv,j,2);
346 389832 : PetscValidLogicalCollectiveScalar(bv,alpha,3);
347 97458 : PetscValidType(bv,1);
348 97458 : BVCheckSizes(bv,1);
349 :
350 97458 : PetscCheck(j>=0 && j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %" PetscInt_FMT ", the number of columns is %" PetscInt_FMT,j,bv->m);
351 97458 : if (alpha == (PetscScalar)1.0) PetscFunctionReturn(PETSC_SUCCESS);
352 :
353 96164 : PetscCall(PetscLogEventBegin(BV_Scale,bv,0,0,0));
354 96164 : PetscUseTypeMethod(bv,scale,j,alpha);
355 96164 : PetscCall(PetscLogEventEnd(BV_Scale,bv,0,0,0));
356 96164 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
357 96164 : PetscFunctionReturn(PETSC_SUCCESS);
358 : }
359 :
360 2822 : static inline PetscErrorCode BVSetRandomColumn_Private(BV bv,PetscInt k)
361 : {
362 2822 : PetscInt i,low,high;
363 2822 : PetscScalar *px,t;
364 2822 : Vec x;
365 :
366 2822 : PetscFunctionBegin;
367 2822 : PetscCall(BVGetColumn(bv,k,&x));
368 2822 : if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
369 0 : PetscCall(VecGetOwnershipRange(x,&low,&high));
370 0 : PetscCall(VecGetArray(x,&px));
371 0 : for (i=0;i<bv->N;i++) {
372 0 : PetscCall(PetscRandomGetValue(bv->rand,&t));
373 0 : if (i>=low && i<high) px[i-low] = t;
374 : }
375 0 : PetscCall(VecRestoreArray(x,&px));
376 2822 : } else PetscCall(VecSetRandom(x,bv->rand));
377 2822 : PetscCall(BVRestoreColumn(bv,k,&x));
378 2822 : PetscFunctionReturn(PETSC_SUCCESS);
379 : }
380 :
381 106 : static inline PetscErrorCode BVSetRandomNormalColumn_Private(BV bv,PetscInt k,Vec w1,Vec w2)
382 : {
383 106 : PetscInt i,low,high;
384 106 : PetscScalar *px,s,t;
385 106 : Vec x;
386 :
387 106 : PetscFunctionBegin;
388 106 : PetscCall(BVGetColumn(bv,k,&x));
389 106 : if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
390 0 : PetscCall(VecGetOwnershipRange(x,&low,&high));
391 0 : PetscCall(VecGetArray(x,&px));
392 0 : for (i=0;i<bv->N;i++) {
393 0 : PetscCall(PetscRandomGetValue(bv->rand,&s));
394 0 : PetscCall(PetscRandomGetValue(bv->rand,&t));
395 0 : if (i>=low && i<high) {
396 : #if defined(PETSC_USE_COMPLEX)
397 0 : px[i-low] = PetscCMPLX(PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(s)))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(t)),PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(s)))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(t)));
398 : #else
399 : px[i-low] = PetscSqrtReal(-2.0*PetscLogReal(s))*PetscCosReal(2.0*PETSC_PI*t);
400 : #endif
401 : }
402 : }
403 0 : PetscCall(VecRestoreArray(x,&px));
404 106 : } else PetscCall(VecSetRandomNormal(x,bv->rand,w1,w2));
405 106 : PetscCall(BVRestoreColumn(bv,k,&x));
406 106 : PetscFunctionReturn(PETSC_SUCCESS);
407 : }
408 :
409 1094 : static inline PetscErrorCode BVSetRandomSignColumn_Private(BV bv,PetscInt k)
410 : {
411 1094 : PetscInt i,low,high;
412 1094 : PetscScalar *px,t;
413 1094 : Vec x;
414 :
415 1094 : PetscFunctionBegin;
416 1094 : PetscCall(BVGetColumn(bv,k,&x));
417 1094 : PetscCall(VecGetOwnershipRange(x,&low,&high));
418 1094 : if (bv->rrandom) { /* generate the same vector irrespective of number of processes */
419 0 : PetscCall(VecGetArray(x,&px));
420 0 : for (i=0;i<bv->N;i++) {
421 0 : PetscCall(PetscRandomGetValue(bv->rand,&t));
422 0 : if (i>=low && i<high) px[i-low] = (PetscRealPart(t)<0.5)? -1.0: 1.0;
423 : }
424 0 : PetscCall(VecRestoreArray(x,&px));
425 : } else {
426 1094 : PetscCall(VecSetRandom(x,bv->rand));
427 1094 : PetscCall(VecGetArray(x,&px));
428 231482 : for (i=low;i<high;i++) {
429 343402 : px[i-low] = (PetscRealPart(px[i-low])<0.5)? -1.0: 1.0;
430 : }
431 1094 : PetscCall(VecRestoreArray(x,&px));
432 : }
433 1094 : PetscCall(BVRestoreColumn(bv,k,&x));
434 1094 : PetscFunctionReturn(PETSC_SUCCESS);
435 : }
436 :
437 : /*@
438 : BVSetRandom - Set the columns of a BV to random numbers.
439 :
440 : Logically Collective
441 :
442 : Input Parameters:
443 : . bv - basis vectors
444 :
445 : Note:
446 : All active columns (except the leading ones) are modified.
447 :
448 : Level: advanced
449 :
450 : .seealso: BVSetRandomContext(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetRandomSign(), BVSetRandomCond(), BVSetActiveColumns()
451 : @*/
452 0 : PetscErrorCode BVSetRandom(BV bv)
453 : {
454 0 : PetscInt k;
455 :
456 0 : PetscFunctionBegin;
457 0 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
458 0 : PetscValidType(bv,1);
459 0 : BVCheckSizes(bv,1);
460 :
461 0 : PetscCall(BVGetRandomContext(bv,&bv->rand));
462 0 : PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
463 0 : for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomColumn_Private(bv,k));
464 0 : PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
465 0 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
466 0 : PetscFunctionReturn(PETSC_SUCCESS);
467 : }
468 :
469 : /*@
470 : BVSetRandomColumn - Set one column of a BV to random numbers.
471 :
472 : Logically Collective
473 :
474 : Input Parameters:
475 : + bv - basis vectors
476 : - j - column number to be set
477 :
478 : Level: advanced
479 :
480 : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomNormal(), BVSetRandomCond()
481 : @*/
482 2790 : PetscErrorCode BVSetRandomColumn(BV bv,PetscInt j)
483 : {
484 2790 : PetscFunctionBegin;
485 2790 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
486 11160 : PetscValidLogicalCollectiveInt(bv,j,2);
487 2790 : PetscValidType(bv,1);
488 2790 : BVCheckSizes(bv,1);
489 2790 : PetscCheck(j>=0 && j<bv->m,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Argument j has wrong value %" PetscInt_FMT ", the number of columns is %" PetscInt_FMT,j,bv->m);
490 :
491 2790 : PetscCall(BVGetRandomContext(bv,&bv->rand));
492 2790 : PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
493 2790 : PetscCall(BVSetRandomColumn_Private(bv,j));
494 2790 : PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
495 2790 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
496 2790 : PetscFunctionReturn(PETSC_SUCCESS);
497 : }
498 :
499 : /*@
500 : BVSetRandomNormal - Set the columns of a BV to random numbers with a normal
501 : distribution.
502 :
503 : Logically Collective
504 :
505 : Input Parameter:
506 : . bv - basis vectors
507 :
508 : Notes:
509 : All active columns (except the leading ones) are modified.
510 :
511 : Other functions such as BVSetRandom(), BVSetRandomColumn(), and BVSetRandomCond()
512 : produce random numbers with a uniform distribution. This function returns values
513 : that fit a normal distribution (Gaussian).
514 :
515 : Developer Notes:
516 : The current implementation obtains each of the columns by applying the Box-Muller
517 : transform on two random vectors with uniformly distributed entries.
518 :
519 : Level: advanced
520 :
521 : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomCond(), BVSetActiveColumns()
522 : @*/
523 10 : PetscErrorCode BVSetRandomNormal(BV bv)
524 : {
525 10 : PetscInt k;
526 10 : Vec w1=NULL,w2=NULL;
527 :
528 10 : PetscFunctionBegin;
529 10 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
530 10 : PetscValidType(bv,1);
531 10 : BVCheckSizes(bv,1);
532 :
533 10 : PetscCall(BVGetRandomContext(bv,&bv->rand));
534 10 : if (!bv->rrandom) {
535 10 : PetscCall(BVCreateVec(bv,&w1));
536 10 : PetscCall(BVCreateVec(bv,&w2));
537 : }
538 10 : PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
539 116 : for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomNormalColumn_Private(bv,k,w1,w2));
540 10 : PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
541 10 : if (!bv->rrandom) {
542 10 : PetscCall(VecDestroy(&w1));
543 10 : PetscCall(VecDestroy(&w2));
544 : }
545 10 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
546 10 : PetscFunctionReturn(PETSC_SUCCESS);
547 : }
548 :
549 : /*@
550 : BVSetRandomSign - Set the entries of a BV to values 1 or -1 with equal
551 : probability.
552 :
553 : Logically Collective
554 :
555 : Input Parameter:
556 : . bv - basis vectors
557 :
558 : Notes:
559 : All active columns (except the leading ones) are modified.
560 :
561 : This function is used, e.g., in contour integral methods when estimating
562 : the number of eigenvalues enclosed by the contour via an unbiased
563 : estimator of tr(f(A)) [Bai et al., JCAM 1996].
564 :
565 : Developer Notes:
566 : The current implementation obtains random numbers and then replaces them
567 : with -1 or 1 depending on the value being less than 0.5 or not.
568 :
569 : Level: advanced
570 :
571 : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetActiveColumns()
572 : @*/
573 93 : PetscErrorCode BVSetRandomSign(BV bv)
574 : {
575 93 : PetscScalar low,high;
576 93 : PetscInt k;
577 :
578 93 : PetscFunctionBegin;
579 93 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
580 93 : PetscValidType(bv,1);
581 93 : BVCheckSizes(bv,1);
582 :
583 93 : PetscCall(BVGetRandomContext(bv,&bv->rand));
584 93 : PetscCall(PetscRandomGetInterval(bv->rand,&low,&high));
585 93 : PetscCheck(PetscRealPart(low)==0.0 && PetscRealPart(high)==1.0,PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_WRONGSTATE,"The PetscRandom object in the BV must have interval [0,1]");
586 93 : PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
587 1187 : for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomSignColumn_Private(bv,k));
588 93 : PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
589 93 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
590 93 : PetscFunctionReturn(PETSC_SUCCESS);
591 : }
592 :
593 : /*@
594 : BVSetRandomCond - Set the columns of a BV to random numbers, in a way that
595 : the generated matrix has a given condition number.
596 :
597 : Logically Collective
598 :
599 : Input Parameters:
600 : + bv - basis vectors
601 : - condn - condition number
602 :
603 : Note:
604 : All active columns (except the leading ones) are modified.
605 :
606 : Level: advanced
607 :
608 : .seealso: BVSetRandomContext(), BVSetRandom(), BVSetRandomColumn(), BVSetRandomNormal(), BVSetActiveColumns()
609 : @*/
610 4 : PetscErrorCode BVSetRandomCond(BV bv,PetscReal condn)
611 : {
612 4 : PetscInt k,i;
613 4 : PetscScalar *eig,*d;
614 4 : DS ds;
615 4 : Mat A,X,Xt,M,G;
616 :
617 4 : PetscFunctionBegin;
618 4 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
619 4 : PetscValidType(bv,1);
620 4 : BVCheckSizes(bv,1);
621 :
622 4 : PetscCall(BVGetRandomContext(bv,&bv->rand));
623 4 : PetscCall(PetscLogEventBegin(BV_SetRandom,bv,0,0,0));
624 : /* B = rand(n,k) */
625 36 : for (k=bv->l;k<bv->k;k++) PetscCall(BVSetRandomColumn_Private(bv,k));
626 4 : PetscCall(DSCreate(PetscObjectComm((PetscObject)bv),&ds));
627 4 : PetscCall(DSSetType(ds,DSHEP));
628 4 : PetscCall(DSAllocate(ds,bv->m));
629 4 : PetscCall(DSSetDimensions(ds,bv->k,bv->l,bv->k));
630 : /* [V,S] = eig(B'*B) */
631 4 : PetscCall(DSGetMat(ds,DS_MAT_A,&A));
632 4 : PetscCall(BVDot(bv,bv,A));
633 4 : PetscCall(DSRestoreMat(ds,DS_MAT_A,&A));
634 4 : PetscCall(PetscMalloc1(bv->m,&eig));
635 4 : PetscCall(DSSolve(ds,eig,NULL));
636 4 : PetscCall(DSSynchronize(ds,eig,NULL));
637 4 : PetscCall(DSVectors(ds,DS_MAT_X,NULL,NULL));
638 : /* M = diag(linspace(1/condn,1,n)./sqrt(diag(S)))' */
639 4 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,bv->k,bv->k,NULL,&M));
640 4 : PetscCall(MatZeroEntries(M));
641 4 : PetscCall(MatDenseGetArray(M,&d));
642 36 : for (i=0;i<bv->k;i++) d[i+i*bv->m] = (1.0/condn+(1.0-1.0/condn)/(bv->k-1)*i)/PetscSqrtScalar(eig[i]);
643 4 : PetscCall(MatDenseRestoreArray(M,&d));
644 : /* G = X*M*X' */
645 4 : PetscCall(DSGetMat(ds,DS_MAT_X,&X));
646 4 : PetscCall(MatTranspose(X,MAT_INITIAL_MATRIX,&Xt));
647 4 : PetscCall(MatProductCreate(Xt,M,NULL,&G));
648 4 : PetscCall(MatProductSetType(G,MATPRODUCT_PtAP));
649 4 : PetscCall(MatProductSetFromOptions(G));
650 4 : PetscCall(MatProductSymbolic(G));
651 4 : PetscCall(MatProductNumeric(G));
652 4 : PetscCall(MatProductClear(G));
653 4 : PetscCall(DSRestoreMat(ds,DS_MAT_X,&X));
654 4 : PetscCall(MatDestroy(&Xt));
655 4 : PetscCall(MatDestroy(&M));
656 : /* B = B*G */
657 4 : PetscCall(BVMultInPlace(bv,G,bv->l,bv->k));
658 4 : PetscCall(MatDestroy(&G));
659 4 : PetscCall(PetscFree(eig));
660 4 : PetscCall(DSDestroy(&ds));
661 4 : PetscCall(PetscLogEventEnd(BV_SetRandom,bv,0,0,0));
662 4 : PetscCall(PetscObjectStateIncrease((PetscObject)bv));
663 4 : PetscFunctionReturn(PETSC_SUCCESS);
664 : }
665 :
666 : /*@
667 : BVMatMult - Computes the matrix-vector product for each column, Y=A*V.
668 :
669 : Neighbor-wise Collective
670 :
671 : Input Parameters:
672 : + V - basis vectors context
673 : - A - the matrix
674 :
675 : Output Parameter:
676 : . Y - the result
677 :
678 : Notes:
679 : Both V and Y must be distributed in the same manner. Only active columns
680 : (excluding the leading ones) are processed.
681 : In the result Y, columns are overwritten starting from the leading ones.
682 : The number of active columns in V and Y should match, although they need
683 : not be the same columns.
684 :
685 : It is possible to choose whether the computation is done column by column
686 : or as a Mat-Mat product, see BVSetMatMultMethod().
687 :
688 : Level: beginner
689 :
690 : .seealso: BVCopy(), BVSetActiveColumns(), BVMatMultColumn(), BVMatMultTranspose(), BVMatMultHermitianTranspose(), BVSetMatMultMethod()
691 : @*/
692 11316 : PetscErrorCode BVMatMult(BV V,Mat A,BV Y)
693 : {
694 11316 : PetscInt M,N,m,n;
695 :
696 11316 : PetscFunctionBegin;
697 11316 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
698 11316 : PetscValidType(V,1);
699 11316 : BVCheckSizes(V,1);
700 11316 : BVCheckOp(V,1,matmult);
701 11316 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
702 11316 : PetscValidType(A,2);
703 11316 : PetscValidHeaderSpecific(Y,BV_CLASSID,3);
704 11316 : PetscValidType(Y,3);
705 11316 : BVCheckSizes(Y,3);
706 11316 : PetscCheckSameComm(V,1,A,2);
707 11316 : PetscCheckSameTypeAndComm(V,1,Y,3);
708 :
709 11316 : PetscCall(MatGetSize(A,&M,&N));
710 11316 : PetscCall(MatGetLocalSize(A,&m,&n));
711 11316 : PetscCheck(M==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,M,Y->N);
712 11316 : PetscCheck(m==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,m,Y->n);
713 11316 : PetscCheck(N==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,N,V->N);
714 11316 : PetscCheck(n==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,n,V->n);
715 11316 : PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);
716 :
717 11316 : PetscCall(PetscLogEventBegin(BV_MatMult,V,A,Y,0));
718 11316 : PetscUseTypeMethod(V,matmult,A,Y);
719 11316 : PetscCall(PetscLogEventEnd(BV_MatMult,V,A,Y,0));
720 11316 : PetscCall(PetscObjectStateIncrease((PetscObject)Y));
721 11316 : PetscFunctionReturn(PETSC_SUCCESS);
722 : }
723 :
724 : /*@
725 : BVMatMultTranspose - Computes the matrix-vector product with the transpose
726 : of a matrix for each column, Y=A^T*V.
727 :
728 : Neighbor-wise Collective
729 :
730 : Input Parameters:
731 : + V - basis vectors context
732 : - A - the matrix
733 :
734 : Output Parameter:
735 : . Y - the result
736 :
737 : Notes:
738 : Both V and Y must be distributed in the same manner. Only active columns
739 : (excluding the leading ones) are processed.
740 : In the result Y, columns are overwritten starting from the leading ones.
741 : The number of active columns in V and Y should match, although they need
742 : not be the same columns.
743 :
744 : Currently implemented via MatCreateTranspose().
745 :
746 : Level: beginner
747 :
748 : .seealso: BVMatMult(), BVMatMultHermitianTranspose()
749 : @*/
750 0 : PetscErrorCode BVMatMultTranspose(BV V,Mat A,BV Y)
751 : {
752 0 : PetscInt M,N,m,n;
753 0 : Mat AT;
754 :
755 0 : PetscFunctionBegin;
756 0 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
757 0 : PetscValidType(V,1);
758 0 : BVCheckSizes(V,1);
759 0 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
760 0 : PetscValidType(A,2);
761 0 : PetscValidHeaderSpecific(Y,BV_CLASSID,3);
762 0 : PetscValidType(Y,3);
763 0 : BVCheckSizes(Y,3);
764 0 : PetscCheckSameComm(V,1,A,2);
765 0 : PetscCheckSameTypeAndComm(V,1,Y,3);
766 :
767 0 : PetscCall(MatGetSize(A,&M,&N));
768 0 : PetscCall(MatGetLocalSize(A,&m,&n));
769 0 : PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
770 0 : PetscCheck(m==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,m,V->n);
771 0 : PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
772 0 : PetscCheck(n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,n,Y->n);
773 0 : PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);
774 :
775 0 : PetscCall(MatCreateTranspose(A,&AT));
776 0 : PetscCall(BVMatMult(V,AT,Y));
777 0 : PetscCall(MatDestroy(&AT));
778 0 : PetscFunctionReturn(PETSC_SUCCESS);
779 : }
780 :
781 : /*@
782 : BVMatMultHermitianTranspose - Computes the matrix-vector product with the
783 : conjugate transpose of a matrix for each column, Y=A^H*V.
784 :
785 : Neighbor-wise Collective
786 :
787 : Input Parameters:
788 : + V - basis vectors context
789 : - A - the matrix
790 :
791 : Output Parameter:
792 : . Y - the result
793 :
794 : Note:
795 : Both V and Y must be distributed in the same manner. Only active columns
796 : (excluding the leading ones) are processed.
797 : In the result Y, columns are overwritten starting from the leading ones.
798 : The number of active columns in V and Y should match, although they need
799 : not be the same columns.
800 :
801 : Currently implemented via MatCreateHermitianTranspose().
802 :
803 : Level: beginner
804 :
805 : .seealso: BVMatMult(), BVMatMultTranspose()
806 : @*/
807 91 : PetscErrorCode BVMatMultHermitianTranspose(BV V,Mat A,BV Y)
808 : {
809 91 : PetscInt j,M,N,m,n;
810 91 : Vec v,y;
811 :
812 91 : PetscFunctionBegin;
813 91 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
814 91 : PetscValidType(V,1);
815 91 : BVCheckSizes(V,1);
816 91 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
817 91 : PetscValidType(A,2);
818 91 : PetscValidHeaderSpecific(Y,BV_CLASSID,3);
819 91 : PetscValidType(Y,3);
820 91 : BVCheckSizes(Y,3);
821 91 : PetscCheckSameComm(V,1,A,2);
822 91 : PetscCheckSameTypeAndComm(V,1,Y,3);
823 :
824 91 : PetscCall(MatGetSize(A,&M,&N));
825 91 : PetscCall(MatGetLocalSize(A,&m,&n));
826 91 : PetscCheck(M==V->N,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,M,V->N);
827 91 : PetscCheck(m==V->n,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_INCOMP,"Mismatching local row dimension A %" PetscInt_FMT ", V %" PetscInt_FMT,m,V->n);
828 91 : PetscCheck(N==Y->N,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,N,Y->N);
829 91 : PetscCheck(n==Y->n,PetscObjectComm((PetscObject)Y),PETSC_ERR_ARG_INCOMP,"Mismatching local column dimension A %" PetscInt_FMT ", Y %" PetscInt_FMT,n,Y->n);
830 91 : PetscCheck(V->k-V->l==Y->k-Y->l,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Y has %" PetscInt_FMT " active columns, should match %" PetscInt_FMT " active columns in V",Y->k-Y->l,V->k-V->l);
831 :
832 : /* TODO: recover this code if PETSc ever gets MATPRODUCT_AhB done
833 : PetscCall(MatCreateHermitianTranspose(A,&AH));
834 : PetscCall(BVMatMult(V,AH,Y));
835 : PetscCall(MatDestroy(&AH));
836 : */
837 402 : for (j=0;j<V->k-V->l;j++) {
838 311 : PetscCall(BVGetColumn(V,V->l+j,&v));
839 311 : PetscCall(BVGetColumn(Y,Y->l+j,&y));
840 311 : PetscCall(MatMultHermitianTranspose(A,v,y));
841 311 : PetscCall(BVRestoreColumn(V,V->l+j,&v));
842 311 : PetscCall(BVRestoreColumn(Y,Y->l+j,&y));
843 : }
844 91 : PetscFunctionReturn(PETSC_SUCCESS);
845 : }
846 :
847 : /*@
848 : BVMatMultColumn - Computes the matrix-vector product for a specified
849 : column, storing the result in the next column v_{j+1}=A*v_j.
850 :
851 : Neighbor-wise Collective
852 :
853 : Input Parameters:
854 : + V - basis vectors context
855 : . A - the matrix
856 : - j - the column
857 :
858 : Level: beginner
859 :
860 : .seealso: BVMatMult(), BVMatMultTransposeColumn(), BVMatMultHermitianTransposeColumn()
861 : @*/
862 61170 : PetscErrorCode BVMatMultColumn(BV V,Mat A,PetscInt j)
863 : {
864 61170 : Vec vj,vj1;
865 :
866 61170 : PetscFunctionBegin;
867 61170 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
868 61170 : PetscValidType(V,1);
869 61170 : BVCheckSizes(V,1);
870 61170 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
871 61170 : PetscCheckSameComm(V,1,A,2);
872 244680 : PetscValidLogicalCollectiveInt(V,j,3);
873 61170 : PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
874 61170 : PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);
875 :
876 61170 : PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
877 61170 : PetscCall(BVGetColumn(V,j,&vj));
878 61170 : PetscCall(BVGetColumn(V,j+1,&vj1));
879 61170 : PetscCall(MatMult(A,vj,vj1));
880 61170 : PetscCall(BVRestoreColumn(V,j,&vj));
881 61170 : PetscCall(BVRestoreColumn(V,j+1,&vj1));
882 61170 : PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
883 61170 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
884 61170 : PetscFunctionReturn(PETSC_SUCCESS);
885 : }
886 :
887 : /*@
888 : BVMatMultTransposeColumn - Computes the transpose matrix-vector product for a
889 : specified column, storing the result in the next column v_{j+1}=A^T*v_j.
890 :
891 : Neighbor-wise Collective
892 :
893 : Input Parameters:
894 : + V - basis vectors context
895 : . A - the matrix
896 : - j - the column
897 :
898 : Level: beginner
899 :
900 : .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultHermitianTransposeColumn()
901 : @*/
902 0 : PetscErrorCode BVMatMultTransposeColumn(BV V,Mat A,PetscInt j)
903 : {
904 0 : Vec vj,vj1;
905 :
906 0 : PetscFunctionBegin;
907 0 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
908 0 : PetscValidType(V,1);
909 0 : BVCheckSizes(V,1);
910 0 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
911 0 : PetscCheckSameComm(V,1,A,2);
912 0 : PetscValidLogicalCollectiveInt(V,j,3);
913 0 : PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
914 0 : PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);
915 :
916 0 : PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
917 0 : PetscCall(BVGetColumn(V,j,&vj));
918 0 : PetscCall(BVGetColumn(V,j+1,&vj1));
919 0 : PetscCall(MatMultTranspose(A,vj,vj1));
920 0 : PetscCall(BVRestoreColumn(V,j,&vj));
921 0 : PetscCall(BVRestoreColumn(V,j+1,&vj1));
922 0 : PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
923 0 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
924 0 : PetscFunctionReturn(PETSC_SUCCESS);
925 : }
926 :
927 : /*@
928 : BVMatMultHermitianTransposeColumn - Computes the conjugate-transpose matrix-vector
929 : product for a specified column, storing the result in the next column v_{j+1}=A^H*v_j.
930 :
931 : Neighbor-wise Collective
932 :
933 : Input Parameters:
934 : + V - basis vectors context
935 : . A - the matrix
936 : - j - the column
937 :
938 : Level: beginner
939 :
940 : .seealso: BVMatMult(), BVMatMultColumn(), BVMatMultTransposeColumn()
941 : @*/
942 0 : PetscErrorCode BVMatMultHermitianTransposeColumn(BV V,Mat A,PetscInt j)
943 : {
944 0 : Vec vj,vj1;
945 :
946 0 : PetscFunctionBegin;
947 0 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
948 0 : PetscValidType(V,1);
949 0 : BVCheckSizes(V,1);
950 0 : PetscValidHeaderSpecific(A,MAT_CLASSID,2);
951 0 : PetscCheckSameComm(V,1,A,2);
952 0 : PetscValidLogicalCollectiveInt(V,j,3);
953 0 : PetscCheck(j>=0,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Index j must be non-negative");
954 0 : PetscCheck(j+1<V->m,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_OUTOFRANGE,"Result should go in index j+1=%" PetscInt_FMT " but BV only has %" PetscInt_FMT " columns",j+1,V->m);
955 :
956 0 : PetscCall(PetscLogEventBegin(BV_MatMultVec,V,A,0,0));
957 0 : PetscCall(BVGetColumn(V,j,&vj));
958 0 : PetscCall(BVGetColumn(V,j+1,&vj1));
959 0 : PetscCall(MatMultHermitianTranspose(A,vj,vj1));
960 0 : PetscCall(BVRestoreColumn(V,j,&vj));
961 0 : PetscCall(BVRestoreColumn(V,j+1,&vj1));
962 0 : PetscCall(PetscLogEventEnd(BV_MatMultVec,V,A,0,0));
963 0 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
964 0 : PetscFunctionReturn(PETSC_SUCCESS);
965 : }
|
