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