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 implemented as an array of independent Vecs
12 : */
13 :
14 : #include <slepc/private/bvimpl.h>
15 :
16 : typedef struct {
17 : Vec *V;
18 : PetscInt vmip; /* Version of BVMultInPlace:
19 : 0: memory-efficient version, uses VecGetArray (default in CPU)
20 : 1: version that allocates (e-s) work vectors in every call (default in GPU) */
21 : } BV_VECS;
22 :
23 101 : static PetscErrorCode BVMult_Vecs(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
24 : {
25 101 : BV_VECS *y = (BV_VECS*)Y->data,*x = (BV_VECS*)X->data;
26 101 : PetscScalar *s=NULL;
27 101 : const PetscScalar *q;
28 101 : PetscInt i,j,ldq;
29 101 : PetscBool trivial=(alpha==1.0)?PETSC_TRUE:PETSC_FALSE;
30 :
31 101 : PetscFunctionBegin;
32 101 : if (Q) {
33 81 : PetscCall(MatDenseGetLDA(Q,&ldq));
34 81 : if (!trivial) {
35 81 : PetscCall(BVAllocateWork_Private(Y,X->k-X->l));
36 81 : s = Y->work;
37 : }
38 81 : PetscCall(MatDenseGetArrayRead(Q,&q));
39 546 : for (j=Y->l;j<Y->k;j++) {
40 465 : PetscCall(VecScale(y->V[Y->nc+j],beta));
41 465 : if (!trivial) {
42 2814 : for (i=X->l;i<X->k;i++) s[i-X->l] = alpha*q[i+j*ldq];
43 0 : } else s = (PetscScalar*)(q+j*ldq+X->l);
44 465 : PetscCall(VecMAXPY(y->V[Y->nc+j],X->k-X->l,s,x->V+X->nc+X->l));
45 : }
46 81 : PetscCall(MatDenseRestoreArrayRead(Q,&q));
47 : } else {
48 100 : for (j=0;j<Y->k-Y->l;j++) {
49 80 : PetscCall(VecScale(y->V[Y->nc+Y->l+j],beta));
50 80 : PetscCall(VecAXPY(y->V[Y->nc+Y->l+j],alpha,x->V[X->nc+X->l+j]));
51 : }
52 : }
53 101 : PetscFunctionReturn(PETSC_SUCCESS);
54 : }
55 :
56 979 : static PetscErrorCode BVMultVec_Vecs(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
57 : {
58 979 : BV_VECS *x = (BV_VECS*)X->data;
59 979 : PetscScalar *s=NULL,*qq=q;
60 979 : PetscInt i;
61 979 : PetscBool trivial=(alpha==1.0)?PETSC_TRUE:PETSC_FALSE;
62 :
63 979 : PetscFunctionBegin;
64 979 : if (!trivial) {
65 979 : PetscCall(BVAllocateWork_Private(X,X->k-X->l));
66 979 : s = X->work;
67 : }
68 979 : if (!q) PetscCall(VecGetArray(X->buffer,&qq));
69 979 : PetscCall(VecScale(y,beta));
70 979 : if (!trivial) {
71 7462 : for (i=0;i<X->k-X->l;i++) s[i] = alpha*qq[i];
72 0 : } else s = qq;
73 979 : PetscCall(VecMAXPY(y,X->k-X->l,s,x->V+X->nc+X->l));
74 979 : if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
75 979 : PetscFunctionReturn(PETSC_SUCCESS);
76 : }
77 :
78 : /*
79 : BVMultInPlace_Vecs_ME - V(:,s:e-1) = V*Q(:,s:e-1) for regular vectors.
80 :
81 : Memory-efficient version, uses VecGetArray (default in CPU)
82 :
83 : Writing V = [ V1 V2 V3 ] and Q(:,s:e-1) = [ Q1 Q2 Q3 ]', where V2
84 : corresponds to the columns s:e-1, the computation is done as
85 : V2 := V2*Q2 + V1*Q1 + V3*Q3
86 : */
87 117 : static PetscErrorCode BVMultInPlace_Vecs_ME(BV V,Mat Q,PetscInt s,PetscInt e)
88 : {
89 117 : BV_VECS *ctx = (BV_VECS*)V->data;
90 117 : const PetscScalar *q;
91 117 : PetscInt i,ldq;
92 :
93 117 : PetscFunctionBegin;
94 117 : PetscCall(MatDenseGetLDA(Q,&ldq));
95 117 : PetscCall(MatDenseGetArrayRead(Q,&q));
96 : /* V2 := V2*Q2 */
97 117 : PetscCall(BVMultInPlace_Vecs_Private(V,V->n,e-s,ldq,ctx->V+V->nc+s,q+s*ldq+s,PETSC_FALSE));
98 : /* V2 += V1*Q1 + V3*Q3 */
99 570 : for (i=s;i<e;i++) {
100 453 : if (PetscUnlikely(s>V->l)) PetscCall(VecMAXPY(ctx->V[V->nc+i],s-V->l,q+i*ldq+V->l,ctx->V+V->nc+V->l));
101 453 : if (V->k>e) PetscCall(VecMAXPY(ctx->V[V->nc+i],V->k-e,q+i*ldq+e,ctx->V+V->nc+e));
102 : }
103 117 : PetscCall(MatDenseRestoreArrayRead(Q,&q));
104 117 : PetscFunctionReturn(PETSC_SUCCESS);
105 : }
106 :
107 : /*
108 : BVMultInPlace_Vecs_Alloc - V(:,s:e-1) = V*Q(:,s:e-1) for regular vectors.
109 :
110 : Version that allocates (e-s) work vectors in every call (default in GPU)
111 : */
112 1 : static PetscErrorCode BVMultInPlace_Vecs_Alloc(BV V,Mat Q,PetscInt s,PetscInt e)
113 : {
114 1 : BV_VECS *ctx = (BV_VECS*)V->data;
115 1 : const PetscScalar *q;
116 1 : PetscInt i,ldq;
117 1 : Vec *W,z;
118 :
119 1 : PetscFunctionBegin;
120 1 : PetscCall(MatDenseGetLDA(Q,&ldq));
121 1 : PetscCall(MatDenseGetArrayRead(Q,&q));
122 1 : PetscCall(BVCreateVec(V,&z));
123 1 : PetscCall(VecDuplicateVecs(z,e-s,&W));
124 1 : PetscCall(VecDestroy(&z));
125 5 : for (i=s;i<e;i++) PetscCall(VecMAXPY(W[i-s],V->k-V->l,q+i*ldq+V->l,ctx->V+V->nc+V->l));
126 5 : for (i=s;i<e;i++) PetscCall(VecCopy(W[i-s],ctx->V[V->nc+i]));
127 1 : PetscCall(VecDestroyVecs(e-s,&W));
128 1 : PetscCall(MatDenseRestoreArrayRead(Q,&q));
129 1 : PetscFunctionReturn(PETSC_SUCCESS);
130 : }
131 :
132 : /*
133 : BVMultInPlaceHermitianTranspose_Vecs - V(:,s:e-1) = V*Q'(:,s:e-1) for regular vectors.
134 : */
135 1 : static PetscErrorCode BVMultInPlaceHermitianTranspose_Vecs(BV V,Mat Q,PetscInt s,PetscInt e)
136 : {
137 1 : BV_VECS *ctx = (BV_VECS*)V->data;
138 1 : const PetscScalar *q;
139 1 : PetscInt i,j,ldq,n;
140 :
141 1 : PetscFunctionBegin;
142 1 : PetscCall(MatGetSize(Q,NULL,&n));
143 1 : PetscCall(MatDenseGetLDA(Q,&ldq));
144 1 : PetscCall(MatDenseGetArrayRead(Q,&q));
145 : /* V2 := V2*Q2' */
146 1 : PetscCall(BVMultInPlace_Vecs_Private(V,V->n,e-s,ldq,ctx->V+V->nc+s,q+s*ldq+s,PETSC_TRUE));
147 : /* V2 += V1*Q1' + V3*Q3' */
148 5 : for (i=s;i<e;i++) {
149 8 : for (j=V->l;j<s;j++) PetscCall(VecAXPY(ctx->V[V->nc+i],q[i+j*ldq],ctx->V[V->nc+j]));
150 20 : for (j=e;j<n;j++) PetscCall(VecAXPY(ctx->V[V->nc+i],q[i+j*ldq],ctx->V[V->nc+j]));
151 : }
152 1 : PetscCall(MatDenseRestoreArrayRead(Q,&q));
153 1 : PetscFunctionReturn(PETSC_SUCCESS);
154 : }
155 :
156 156 : static PetscErrorCode BVDot_Vecs(BV X,BV Y,Mat M)
157 : {
158 156 : BV_VECS *x = (BV_VECS*)X->data,*y = (BV_VECS*)Y->data;
159 156 : PetscScalar *m;
160 156 : PetscInt j,ldm;
161 :
162 156 : PetscFunctionBegin;
163 156 : PetscCall(MatDenseGetLDA(M,&ldm));
164 156 : PetscCall(MatDenseGetArray(M,&m));
165 1043 : for (j=X->l;j<X->k;j++) PetscCall(VecMDot(x->V[X->nc+j],Y->k-Y->l,y->V+Y->nc+Y->l,m+j*ldm+Y->l));
166 156 : PetscCall(MatDenseRestoreArray(M,&m));
167 156 : PetscFunctionReturn(PETSC_SUCCESS);
168 : }
169 :
170 1371 : static PetscErrorCode BVDotVec_Vecs(BV X,Vec y,PetscScalar *q)
171 : {
172 1371 : BV_VECS *x = (BV_VECS*)X->data;
173 1371 : Vec z = y;
174 1371 : PetscScalar *qq=q;
175 :
176 1371 : PetscFunctionBegin;
177 1371 : if (PetscUnlikely(X->matrix)) {
178 57 : PetscCall(BV_IPMatMult(X,y));
179 57 : z = X->Bx;
180 : }
181 1371 : if (!q) PetscCall(VecGetArray(X->buffer,&qq));
182 1371 : PetscCall(VecMDot(z,X->k-X->l,x->V+X->nc+X->l,qq));
183 1371 : if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
184 1371 : PetscFunctionReturn(PETSC_SUCCESS);
185 : }
186 :
187 4 : static PetscErrorCode BVDotVec_Begin_Vecs(BV X,Vec y,PetscScalar *m)
188 : {
189 4 : BV_VECS *x = (BV_VECS*)X->data;
190 4 : Vec z = y;
191 :
192 4 : PetscFunctionBegin;
193 4 : if (PetscUnlikely(X->matrix)) {
194 0 : PetscCall(BV_IPMatMult(X,y));
195 0 : z = X->Bx;
196 : }
197 4 : PetscCall(VecMDotBegin(z,X->k-X->l,x->V+X->nc+X->l,m));
198 4 : PetscFunctionReturn(PETSC_SUCCESS);
199 : }
200 :
201 4 : static PetscErrorCode BVDotVec_End_Vecs(BV X,Vec y,PetscScalar *m)
202 : {
203 4 : BV_VECS *x = (BV_VECS*)X->data;
204 :
205 4 : PetscFunctionBegin;
206 4 : PetscCall(VecMDotEnd(y,X->k-X->l,x->V+X->nc+X->l,m));
207 4 : PetscFunctionReturn(PETSC_SUCCESS);
208 : }
209 :
210 642 : static PetscErrorCode BVScale_Vecs(BV bv,PetscInt j,PetscScalar alpha)
211 : {
212 642 : PetscInt i;
213 642 : BV_VECS *ctx = (BV_VECS*)bv->data;
214 :
215 642 : PetscFunctionBegin;
216 642 : if (PetscUnlikely(j<0)) {
217 42 : for (i=bv->l;i<bv->k;i++) PetscCall(VecScale(ctx->V[bv->nc+i],alpha));
218 637 : } else PetscCall(VecScale(ctx->V[bv->nc+j],alpha));
219 642 : PetscFunctionReturn(PETSC_SUCCESS);
220 : }
221 :
222 114 : static PetscErrorCode BVNorm_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
223 : {
224 114 : PetscInt i;
225 114 : PetscReal nrm;
226 114 : BV_VECS *ctx = (BV_VECS*)bv->data;
227 :
228 114 : PetscFunctionBegin;
229 114 : if (PetscUnlikely(j<0)) {
230 71 : PetscCheck(type==NORM_FROBENIUS,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
231 71 : *val = 0.0;
232 448 : for (i=bv->l;i<bv->k;i++) {
233 377 : PetscCall(VecNorm(ctx->V[bv->nc+i],NORM_2,&nrm));
234 377 : *val += nrm*nrm;
235 : }
236 71 : *val = PetscSqrtReal(*val);
237 43 : } else PetscCall(VecNorm(ctx->V[bv->nc+j],type,val));
238 114 : PetscFunctionReturn(PETSC_SUCCESS);
239 : }
240 :
241 2 : static PetscErrorCode BVNorm_Begin_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
242 : {
243 2 : BV_VECS *ctx = (BV_VECS*)bv->data;
244 :
245 2 : PetscFunctionBegin;
246 2 : PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
247 2 : PetscCall(VecNormBegin(ctx->V[bv->nc+j],type,val));
248 2 : PetscFunctionReturn(PETSC_SUCCESS);
249 : }
250 :
251 2 : static PetscErrorCode BVNorm_End_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
252 : {
253 2 : BV_VECS *ctx = (BV_VECS*)bv->data;
254 :
255 2 : PetscFunctionBegin;
256 2 : PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
257 2 : PetscCall(VecNormEnd(ctx->V[bv->nc+j],type,val));
258 2 : PetscFunctionReturn(PETSC_SUCCESS);
259 : }
260 :
261 6 : static PetscErrorCode BVNormalize_Vecs(BV bv,PetscScalar *eigi)
262 : {
263 6 : BV_VECS *ctx = (BV_VECS*)bv->data;
264 6 : PetscInt i;
265 :
266 6 : PetscFunctionBegin;
267 54 : for (i=bv->l;i<bv->k;i++) {
268 : #if !defined(PETSC_USE_COMPLEX)
269 48 : if (eigi && eigi[i] != 0.0) {
270 6 : PetscCall(VecNormalizeComplex(ctx->V[bv->nc+i],ctx->V[bv->nc+i+1],PETSC_TRUE,NULL));
271 6 : i++;
272 : } else
273 : #endif
274 : {
275 48 : PetscCall(VecNormalize(ctx->V[bv->nc+i],NULL));
276 : }
277 : }
278 6 : PetscFunctionReturn(PETSC_SUCCESS);
279 : }
280 :
281 98 : static PetscErrorCode BVMatMult_Vecs(BV V,Mat A,BV W)
282 : {
283 98 : BV_VECS *v = (BV_VECS*)V->data,*w = (BV_VECS*)W->data;
284 98 : PetscInt j;
285 98 : Mat Vmat,Wmat;
286 :
287 98 : PetscFunctionBegin;
288 98 : if (V->vmm) {
289 7 : PetscCall(BVGetMat(V,&Vmat));
290 7 : PetscCall(BVGetMat(W,&Wmat));
291 7 : PetscCall(MatProductCreateWithMat(A,Vmat,NULL,Wmat));
292 7 : PetscCall(MatProductSetType(Wmat,MATPRODUCT_AB));
293 7 : PetscCall(MatProductSetFromOptions(Wmat));
294 7 : PetscCall(MatProductSymbolic(Wmat));
295 7 : PetscCall(MatProductNumeric(Wmat));
296 7 : PetscCall(MatProductClear(Wmat));
297 7 : PetscCall(BVRestoreMat(V,&Vmat));
298 7 : PetscCall(BVRestoreMat(W,&Wmat));
299 : } else {
300 525 : for (j=0;j<V->k-V->l;j++) PetscCall(MatMult(A,v->V[V->nc+V->l+j],w->V[W->nc+W->l+j]));
301 : }
302 98 : PetscFunctionReturn(PETSC_SUCCESS);
303 : }
304 :
305 99 : static PetscErrorCode BVCopy_Vecs(BV V,BV W)
306 : {
307 99 : BV_VECS *v = (BV_VECS*)V->data,*w = (BV_VECS*)W->data;
308 99 : PetscInt j;
309 :
310 99 : PetscFunctionBegin;
311 742 : for (j=0;j<V->k-V->l;j++) PetscCall(VecCopy(v->V[V->nc+V->l+j],w->V[W->nc+W->l+j]));
312 99 : PetscFunctionReturn(PETSC_SUCCESS);
313 : }
314 :
315 72 : static PetscErrorCode BVCopyColumn_Vecs(BV V,PetscInt j,PetscInt i)
316 : {
317 72 : BV_VECS *v = (BV_VECS*)V->data;
318 :
319 72 : PetscFunctionBegin;
320 72 : PetscCall(VecCopy(v->V[V->nc+j],v->V[V->nc+i]));
321 72 : PetscFunctionReturn(PETSC_SUCCESS);
322 : }
323 :
324 9 : static PetscErrorCode BVResize_Vecs(BV bv,PetscInt m,PetscBool copy)
325 : {
326 9 : BV_VECS *ctx = (BV_VECS*)bv->data;
327 9 : Vec *newV,z;
328 9 : PetscInt j;
329 9 : char str[50];
330 :
331 9 : PetscFunctionBegin;
332 9 : PetscCall(BVCreateVec(bv,&z));
333 9 : PetscCall(VecDuplicateVecs(z,m,&newV));
334 9 : PetscCall(VecDestroy(&z));
335 9 : if (((PetscObject)bv)->name) {
336 111 : for (j=0;j<m;j++) {
337 102 : PetscCall(PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j));
338 102 : PetscCall(PetscObjectSetName((PetscObject)newV[j],str));
339 : }
340 : }
341 9 : if (copy) {
342 45 : for (j=0;j<PetscMin(m,bv->m);j++) PetscCall(VecCopy(ctx->V[j],newV[j]));
343 : }
344 9 : PetscCall(VecDestroyVecs(bv->m,&ctx->V));
345 9 : ctx->V = newV;
346 9 : PetscFunctionReturn(PETSC_SUCCESS);
347 : }
348 :
349 4770 : static PetscErrorCode BVGetColumn_Vecs(BV bv,PetscInt j,Vec *v)
350 : {
351 4770 : BV_VECS *ctx = (BV_VECS*)bv->data;
352 4770 : PetscInt l;
353 :
354 4770 : PetscFunctionBegin;
355 4770 : l = BVAvailableVec;
356 4770 : bv->cv[l] = ctx->V[bv->nc+j];
357 4770 : PetscFunctionReturn(PETSC_SUCCESS);
358 : }
359 :
360 4770 : static PetscErrorCode BVRestoreColumn_Vecs(BV bv,PetscInt j,Vec *v)
361 : {
362 4770 : PetscInt l;
363 :
364 4770 : PetscFunctionBegin;
365 4770 : l = (j==bv->ci[0])? 0: 1;
366 4770 : bv->cv[l] = NULL;
367 4770 : PetscFunctionReturn(PETSC_SUCCESS);
368 : }
369 :
370 67 : static PetscErrorCode BVGetArray_Vecs(BV bv,PetscScalar **a)
371 : {
372 67 : BV_VECS *ctx = (BV_VECS*)bv->data;
373 67 : PetscInt j;
374 67 : const PetscScalar *p;
375 :
376 67 : PetscFunctionBegin;
377 67 : PetscCall(PetscMalloc1((bv->nc+bv->m)*bv->n,a));
378 616 : for (j=0;j<bv->nc+bv->m;j++) {
379 549 : PetscCall(VecGetArrayRead(ctx->V[j],&p));
380 549 : PetscCall(PetscArraycpy(*a+j*bv->n,p,bv->n));
381 549 : PetscCall(VecRestoreArrayRead(ctx->V[j],&p));
382 : }
383 67 : PetscFunctionReturn(PETSC_SUCCESS);
384 : }
385 :
386 67 : static PetscErrorCode BVRestoreArray_Vecs(BV bv,PetscScalar **a)
387 : {
388 67 : BV_VECS *ctx = (BV_VECS*)bv->data;
389 67 : PetscInt j;
390 67 : PetscScalar *p;
391 :
392 67 : PetscFunctionBegin;
393 616 : for (j=0;j<bv->nc+bv->m;j++) {
394 549 : PetscCall(VecGetArray(ctx->V[j],&p));
395 549 : PetscCall(PetscArraycpy(p,*a+j*bv->n,bv->n));
396 549 : PetscCall(VecRestoreArray(ctx->V[j],&p));
397 : }
398 67 : PetscCall(PetscFree(*a));
399 67 : PetscFunctionReturn(PETSC_SUCCESS);
400 : }
401 :
402 8 : static PetscErrorCode BVGetArrayRead_Vecs(BV bv,const PetscScalar **a)
403 : {
404 8 : BV_VECS *ctx = (BV_VECS*)bv->data;
405 8 : PetscInt j;
406 8 : const PetscScalar *p;
407 :
408 8 : PetscFunctionBegin;
409 8 : PetscCall(PetscMalloc1((bv->nc+bv->m)*bv->n,(PetscScalar**)a));
410 54 : for (j=0;j<bv->nc+bv->m;j++) {
411 46 : PetscCall(VecGetArrayRead(ctx->V[j],&p));
412 46 : PetscCall(PetscArraycpy((PetscScalar*)*a+j*bv->n,p,bv->n));
413 46 : PetscCall(VecRestoreArrayRead(ctx->V[j],&p));
414 : }
415 8 : PetscFunctionReturn(PETSC_SUCCESS);
416 : }
417 :
418 8 : static PetscErrorCode BVRestoreArrayRead_Vecs(BV bv,const PetscScalar **a)
419 : {
420 8 : PetscFunctionBegin;
421 8 : PetscCall(PetscFree(*a));
422 8 : PetscFunctionReturn(PETSC_SUCCESS);
423 : }
424 :
425 : /*
426 : Sets the value of vmip flag and resets ops->multinplace accordingly
427 : */
428 662 : static inline PetscErrorCode BVVecsSetVmip(BV bv,PetscInt vmip)
429 : {
430 662 : typedef PetscErrorCode (*fmultinplace)(BV,Mat,PetscInt,PetscInt);
431 662 : fmultinplace multinplace[2] = {BVMultInPlace_Vecs_ME, BVMultInPlace_Vecs_Alloc};
432 662 : BV_VECS *ctx = (BV_VECS*)bv->data;
433 :
434 662 : PetscFunctionBegin;
435 662 : ctx->vmip = vmip;
436 662 : bv->ops->multinplace = multinplace[vmip];
437 662 : PetscFunctionReturn(PETSC_SUCCESS);
438 : }
439 :
440 120 : static PetscErrorCode BVSetFromOptions_Vecs(BV bv,PetscOptionItems *PetscOptionsObject)
441 : {
442 120 : BV_VECS *ctx = (BV_VECS*)bv->data;
443 :
444 120 : PetscFunctionBegin;
445 120 : PetscOptionsHeadBegin(PetscOptionsObject,"BV Vecs Options");
446 :
447 120 : PetscCall(PetscOptionsRangeInt("-bv_vecs_vmip","Version of BVMultInPlace operation","",ctx->vmip,&ctx->vmip,NULL,0,1));
448 120 : PetscCall(BVVecsSetVmip(bv,ctx->vmip));
449 :
450 120 : PetscOptionsHeadEnd();
451 120 : PetscFunctionReturn(PETSC_SUCCESS);
452 : }
453 :
454 12 : PetscErrorCode BVView_Vecs(BV bv,PetscViewer viewer)
455 : {
456 12 : BV_VECS *ctx = (BV_VECS*)bv->data;
457 12 : PetscInt j;
458 12 : PetscViewerFormat format;
459 12 : PetscBool isascii,ismatlab=PETSC_FALSE;
460 12 : const char *bvname,*name;
461 :
462 12 : PetscFunctionBegin;
463 12 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
464 12 : if (isascii) {
465 12 : PetscCall(PetscViewerGetFormat(viewer,&format));
466 12 : if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
467 10 : if (format == PETSC_VIEWER_ASCII_MATLAB) ismatlab = PETSC_TRUE;
468 : }
469 0 : if (ismatlab) {
470 0 : PetscCall(PetscObjectGetName((PetscObject)bv,&bvname));
471 0 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[];\n",bvname));
472 : }
473 48 : for (j=bv->nc;j<bv->nc+bv->m;j++) {
474 38 : PetscCall(VecView(ctx->V[j],viewer));
475 38 : if (ismatlab) {
476 0 : PetscCall(PetscObjectGetName((PetscObject)ctx->V[j],&name));
477 38 : PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[%s,%s];clear %s\n",bvname,bvname,name,name));
478 : }
479 : }
480 10 : PetscFunctionReturn(PETSC_SUCCESS);
481 : }
482 :
483 331 : static PetscErrorCode BVDestroy_Vecs(BV bv)
484 : {
485 331 : BV_VECS *ctx = (BV_VECS*)bv->data;
486 :
487 331 : PetscFunctionBegin;
488 331 : if (!bv->issplit) PetscCall(VecDestroyVecs(bv->nc+bv->m,&ctx->V));
489 331 : PetscCall(PetscFree(bv->data));
490 331 : PetscFunctionReturn(PETSC_SUCCESS);
491 : }
492 :
493 211 : static PetscErrorCode BVDuplicate_Vecs(BV V,BV W)
494 : {
495 211 : BV_VECS *ctx = (BV_VECS*)V->data;
496 :
497 211 : PetscFunctionBegin;
498 211 : PetscCall(BVVecsSetVmip(W,ctx->vmip));
499 211 : PetscFunctionReturn(PETSC_SUCCESS);
500 : }
501 :
502 331 : SLEPC_EXTERN PetscErrorCode BVCreate_Vecs(BV bv)
503 : {
504 331 : BV_VECS *ctx;
505 331 : PetscInt j,lsplit;
506 331 : PetscBool isgpu;
507 331 : char str[50];
508 331 : BV parent;
509 331 : Vec *Vpar,z;
510 :
511 331 : PetscFunctionBegin;
512 331 : PetscCall(PetscNew(&ctx));
513 331 : bv->data = (void*)ctx;
514 :
515 331 : if (PetscUnlikely(bv->issplit)) {
516 : /* split BV: share the Vecs of the parent BV */
517 80 : parent = bv->splitparent;
518 80 : lsplit = parent->lsplit;
519 80 : Vpar = ((BV_VECS*)parent->data)->V;
520 80 : ctx->V = (bv->issplit==1)? Vpar: Vpar+lsplit;
521 : } else {
522 : /* regular BV: create array of Vecs to store the BV columns */
523 251 : PetscCall(BVCreateVec(bv,&z));
524 251 : PetscCall(VecDuplicateVecs(z,bv->m,&ctx->V));
525 251 : PetscCall(VecDestroy(&z));
526 251 : if (((PetscObject)bv)->name) {
527 1032 : for (j=0;j<bv->m;j++) {
528 915 : PetscCall(PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j));
529 915 : PetscCall(PetscObjectSetName((PetscObject)ctx->V[j],str));
530 : }
531 : }
532 : }
533 331 : if (!bv->ld) bv->ld = bv->n;
534 331 : PetscCheck(bv->ld==bv->n,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"BVVECS does not support a user-defined leading dimension");
535 :
536 331 : if (PetscUnlikely(bv->Acreate)) {
537 18 : for (j=0;j<bv->m;j++) PetscCall(MatGetColumnVector(bv->Acreate,ctx->V[j],j));
538 2 : PetscCall(MatDestroy(&bv->Acreate));
539 : }
540 :
541 : /* Default version of BVMultInPlace */
542 331 : PetscCall(PetscStrcmpAny(bv->vtype,&isgpu,VECSEQCUDA,VECMPICUDA,""));
543 331 : ctx->vmip = isgpu? 1: 0;
544 :
545 : /* Default BVMatMult method */
546 331 : bv->vmm = BV_MATMULT_VECS;
547 :
548 : /* Deferred call to setfromoptions */
549 331 : if (bv->defersfo) {
550 3 : PetscObjectOptionsBegin((PetscObject)bv);
551 1 : PetscCall(BVSetFromOptions_Vecs(bv,PetscOptionsObject));
552 1 : PetscOptionsEnd();
553 : }
554 331 : PetscCall(BVVecsSetVmip(bv,ctx->vmip));
555 :
556 331 : bv->ops->mult = BVMult_Vecs;
557 331 : bv->ops->multvec = BVMultVec_Vecs;
558 331 : bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Vecs;
559 331 : bv->ops->dot = BVDot_Vecs;
560 331 : bv->ops->dotvec = BVDotVec_Vecs;
561 331 : bv->ops->dotvec_begin = BVDotVec_Begin_Vecs;
562 331 : bv->ops->dotvec_end = BVDotVec_End_Vecs;
563 331 : bv->ops->scale = BVScale_Vecs;
564 331 : bv->ops->norm = BVNorm_Vecs;
565 331 : bv->ops->norm_begin = BVNorm_Begin_Vecs;
566 331 : bv->ops->norm_end = BVNorm_End_Vecs;
567 331 : bv->ops->normalize = BVNormalize_Vecs;
568 331 : bv->ops->matmult = BVMatMult_Vecs;
569 331 : bv->ops->copy = BVCopy_Vecs;
570 331 : bv->ops->copycolumn = BVCopyColumn_Vecs;
571 331 : bv->ops->resize = BVResize_Vecs;
572 331 : bv->ops->getcolumn = BVGetColumn_Vecs;
573 331 : bv->ops->restorecolumn = BVRestoreColumn_Vecs;
574 331 : bv->ops->getarray = BVGetArray_Vecs;
575 331 : bv->ops->restorearray = BVRestoreArray_Vecs;
576 331 : bv->ops->getarrayread = BVGetArrayRead_Vecs;
577 331 : bv->ops->restorearrayread = BVRestoreArrayRead_Vecs;
578 331 : bv->ops->getmat = BVGetMat_Default;
579 331 : bv->ops->restoremat = BVRestoreMat_Default;
580 331 : bv->ops->destroy = BVDestroy_Vecs;
581 331 : bv->ops->duplicate = BVDuplicate_Vecs;
582 331 : bv->ops->setfromoptions = BVSetFromOptions_Vecs;
583 331 : bv->ops->view = BVView_Vecs;
584 331 : PetscFunctionReturn(PETSC_SUCCESS);
585 : }
|
