Actual source code: vecs.c
slepc-3.21.1 2024-04-26
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,z;
119: PetscFunctionBegin;
120: PetscCall(MatDenseGetLDA(Q,&ldq));
121: PetscCall(MatDenseGetArrayRead(Q,&q));
122: PetscCall(BVCreateVec(V,&z));
123: PetscCall(VecDuplicateVecs(z,e-s,&W));
124: PetscCall(VecDestroy(&z));
125: 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: for (i=s;i<e;i++) PetscCall(VecCopy(W[i-s],ctx->V[V->nc+i]));
127: PetscCall(VecDestroyVecs(e-s,&W));
128: PetscCall(MatDenseRestoreArrayRead(Q,&q));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /*
133: BVMultInPlaceHermitianTranspose_Vecs - V(:,s:e-1) = V*Q'(:,s:e-1) for regular vectors.
134: */
135: static PetscErrorCode BVMultInPlaceHermitianTranspose_Vecs(BV V,Mat Q,PetscInt s,PetscInt e)
136: {
137: BV_VECS *ctx = (BV_VECS*)V->data;
138: const PetscScalar *q;
139: PetscInt i,j,ldq,n;
141: PetscFunctionBegin;
142: PetscCall(MatGetSize(Q,NULL,&n));
143: PetscCall(MatDenseGetLDA(Q,&ldq));
144: PetscCall(MatDenseGetArrayRead(Q,&q));
145: /* V2 := V2*Q2' */
146: 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: for (i=s;i<e;i++) {
149: for (j=V->l;j<s;j++) PetscCall(VecAXPY(ctx->V[V->nc+i],q[i+j*ldq],ctx->V[V->nc+j]));
150: for (j=e;j<n;j++) PetscCall(VecAXPY(ctx->V[V->nc+i],q[i+j*ldq],ctx->V[V->nc+j]));
151: }
152: PetscCall(MatDenseRestoreArrayRead(Q,&q));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode BVDot_Vecs(BV X,BV Y,Mat M)
157: {
158: BV_VECS *x = (BV_VECS*)X->data,*y = (BV_VECS*)Y->data;
159: PetscScalar *m;
160: PetscInt j,ldm;
162: PetscFunctionBegin;
163: PetscCall(MatDenseGetLDA(M,&ldm));
164: PetscCall(MatDenseGetArray(M,&m));
165: 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: PetscCall(MatDenseRestoreArray(M,&m));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: static PetscErrorCode BVDotVec_Vecs(BV X,Vec y,PetscScalar *q)
171: {
172: BV_VECS *x = (BV_VECS*)X->data;
173: Vec z = y;
174: PetscScalar *qq=q;
176: PetscFunctionBegin;
177: if (PetscUnlikely(X->matrix)) {
178: PetscCall(BV_IPMatMult(X,y));
179: z = X->Bx;
180: }
181: if (!q) PetscCall(VecGetArray(X->buffer,&qq));
182: PetscCall(VecMDot(z,X->k-X->l,x->V+X->nc+X->l,qq));
183: if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode BVDotVec_Begin_Vecs(BV X,Vec y,PetscScalar *m)
188: {
189: BV_VECS *x = (BV_VECS*)X->data;
190: Vec z = y;
192: PetscFunctionBegin;
193: if (PetscUnlikely(X->matrix)) {
194: PetscCall(BV_IPMatMult(X,y));
195: z = X->Bx;
196: }
197: PetscCall(VecMDotBegin(z,X->k-X->l,x->V+X->nc+X->l,m));
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: static PetscErrorCode BVDotVec_End_Vecs(BV X,Vec y,PetscScalar *m)
202: {
203: BV_VECS *x = (BV_VECS*)X->data;
205: PetscFunctionBegin;
206: PetscCall(VecMDotEnd(y,X->k-X->l,x->V+X->nc+X->l,m));
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: static PetscErrorCode BVScale_Vecs(BV bv,PetscInt j,PetscScalar alpha)
211: {
212: PetscInt i;
213: BV_VECS *ctx = (BV_VECS*)bv->data;
215: PetscFunctionBegin;
216: if (PetscUnlikely(j<0)) {
217: for (i=bv->l;i<bv->k;i++) PetscCall(VecScale(ctx->V[bv->nc+i],alpha));
218: } else PetscCall(VecScale(ctx->V[bv->nc+j],alpha));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode BVNorm_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
223: {
224: PetscInt i;
225: PetscReal nrm;
226: BV_VECS *ctx = (BV_VECS*)bv->data;
228: PetscFunctionBegin;
229: if (PetscUnlikely(j<0)) {
230: PetscCheck(type==NORM_FROBENIUS,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
231: *val = 0.0;
232: for (i=bv->l;i<bv->k;i++) {
233: PetscCall(VecNorm(ctx->V[bv->nc+i],NORM_2,&nrm));
234: *val += nrm*nrm;
235: }
236: *val = PetscSqrtReal(*val);
237: } else PetscCall(VecNorm(ctx->V[bv->nc+j],type,val));
238: PetscFunctionReturn(PETSC_SUCCESS);
239: }
241: static PetscErrorCode BVNorm_Begin_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
242: {
243: BV_VECS *ctx = (BV_VECS*)bv->data;
245: PetscFunctionBegin;
246: PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
247: PetscCall(VecNormBegin(ctx->V[bv->nc+j],type,val));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: static PetscErrorCode BVNorm_End_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
252: {
253: BV_VECS *ctx = (BV_VECS*)bv->data;
255: PetscFunctionBegin;
256: PetscCheck(j>=0,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
257: PetscCall(VecNormEnd(ctx->V[bv->nc+j],type,val));
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: static PetscErrorCode BVNormalize_Vecs(BV bv,PetscScalar *eigi)
262: {
263: BV_VECS *ctx = (BV_VECS*)bv->data;
264: PetscInt i;
266: PetscFunctionBegin;
267: for (i=bv->l;i<bv->k;i++) {
268: #if !defined(PETSC_USE_COMPLEX)
269: if (eigi && eigi[i] != 0.0) {
270: PetscCall(VecNormalizeComplex(ctx->V[bv->nc+i],ctx->V[bv->nc+i+1],PETSC_TRUE,NULL));
271: i++;
272: } else
273: #endif
274: {
275: PetscCall(VecNormalize(ctx->V[bv->nc+i],NULL));
276: }
277: }
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: static PetscErrorCode BVMatMult_Vecs(BV V,Mat A,BV W)
282: {
283: BV_VECS *v = (BV_VECS*)V->data,*w = (BV_VECS*)W->data;
284: PetscInt j;
285: Mat Vmat,Wmat;
287: PetscFunctionBegin;
288: if (V->vmm) {
289: PetscCall(BVGetMat(V,&Vmat));
290: PetscCall(BVGetMat(W,&Wmat));
291: PetscCall(MatProductCreateWithMat(A,Vmat,NULL,Wmat));
292: PetscCall(MatProductSetType(Wmat,MATPRODUCT_AB));
293: PetscCall(MatProductSetFromOptions(Wmat));
294: PetscCall(MatProductSymbolic(Wmat));
295: PetscCall(MatProductNumeric(Wmat));
296: PetscCall(MatProductClear(Wmat));
297: PetscCall(BVRestoreMat(V,&Vmat));
298: PetscCall(BVRestoreMat(W,&Wmat));
299: } else {
300: 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: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode BVCopy_Vecs(BV V,BV W)
306: {
307: BV_VECS *v = (BV_VECS*)V->data,*w = (BV_VECS*)W->data;
308: PetscInt j;
310: PetscFunctionBegin;
311: 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: PetscFunctionReturn(PETSC_SUCCESS);
313: }
315: static PetscErrorCode BVCopyColumn_Vecs(BV V,PetscInt j,PetscInt i)
316: {
317: BV_VECS *v = (BV_VECS*)V->data;
319: PetscFunctionBegin;
320: PetscCall(VecCopy(v->V[V->nc+j],v->V[V->nc+i]));
321: PetscFunctionReturn(PETSC_SUCCESS);
322: }
324: static PetscErrorCode BVResize_Vecs(BV bv,PetscInt m,PetscBool copy)
325: {
326: BV_VECS *ctx = (BV_VECS*)bv->data;
327: Vec *newV,z;
328: PetscInt j;
329: char str[50];
331: PetscFunctionBegin;
332: PetscCall(BVCreateVec(bv,&z));
333: PetscCall(VecDuplicateVecs(z,m,&newV));
334: PetscCall(VecDestroy(&z));
335: if (((PetscObject)bv)->name) {
336: for (j=0;j<m;j++) {
337: PetscCall(PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j));
338: PetscCall(PetscObjectSetName((PetscObject)newV[j],str));
339: }
340: }
341: if (copy) {
342: for (j=0;j<PetscMin(m,bv->m);j++) PetscCall(VecCopy(ctx->V[j],newV[j]));
343: }
344: PetscCall(VecDestroyVecs(bv->m,&ctx->V));
345: ctx->V = newV;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: static PetscErrorCode BVGetColumn_Vecs(BV bv,PetscInt j,Vec *v)
350: {
351: BV_VECS *ctx = (BV_VECS*)bv->data;
352: PetscInt l;
354: PetscFunctionBegin;
355: l = BVAvailableVec;
356: bv->cv[l] = ctx->V[bv->nc+j];
357: PetscFunctionReturn(PETSC_SUCCESS);
358: }
360: static PetscErrorCode BVRestoreColumn_Vecs(BV bv,PetscInt j,Vec *v)
361: {
362: PetscInt l;
364: PetscFunctionBegin;
365: l = (j==bv->ci[0])? 0: 1;
366: bv->cv[l] = NULL;
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode BVGetArray_Vecs(BV bv,PetscScalar **a)
371: {
372: BV_VECS *ctx = (BV_VECS*)bv->data;
373: PetscInt j;
374: const PetscScalar *p;
376: PetscFunctionBegin;
377: PetscCall(PetscMalloc1((bv->nc+bv->m)*bv->n,a));
378: for (j=0;j<bv->nc+bv->m;j++) {
379: PetscCall(VecGetArrayRead(ctx->V[j],&p));
380: PetscCall(PetscArraycpy(*a+j*bv->n,p,bv->n));
381: PetscCall(VecRestoreArrayRead(ctx->V[j],&p));
382: }
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: static PetscErrorCode BVRestoreArray_Vecs(BV bv,PetscScalar **a)
387: {
388: BV_VECS *ctx = (BV_VECS*)bv->data;
389: PetscInt j;
390: PetscScalar *p;
392: PetscFunctionBegin;
393: for (j=0;j<bv->nc+bv->m;j++) {
394: PetscCall(VecGetArray(ctx->V[j],&p));
395: PetscCall(PetscArraycpy(p,*a+j*bv->n,bv->n));
396: PetscCall(VecRestoreArray(ctx->V[j],&p));
397: }
398: PetscCall(PetscFree(*a));
399: PetscFunctionReturn(PETSC_SUCCESS);
400: }
402: static PetscErrorCode BVGetArrayRead_Vecs(BV bv,const PetscScalar **a)
403: {
404: BV_VECS *ctx = (BV_VECS*)bv->data;
405: PetscInt j;
406: const PetscScalar *p;
408: PetscFunctionBegin;
409: PetscCall(PetscMalloc1((bv->nc+bv->m)*bv->n,(PetscScalar**)a));
410: for (j=0;j<bv->nc+bv->m;j++) {
411: PetscCall(VecGetArrayRead(ctx->V[j],&p));
412: PetscCall(PetscArraycpy((PetscScalar*)*a+j*bv->n,p,bv->n));
413: PetscCall(VecRestoreArrayRead(ctx->V[j],&p));
414: }
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: static PetscErrorCode BVRestoreArrayRead_Vecs(BV bv,const PetscScalar **a)
419: {
420: PetscFunctionBegin;
421: PetscCall(PetscFree(*a));
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: /*
426: Sets the value of vmip flag and resets ops->multinplace accordingly
427: */
428: static inline PetscErrorCode BVVecsSetVmip(BV bv,PetscInt vmip)
429: {
430: typedef PetscErrorCode (*fmultinplace)(BV,Mat,PetscInt,PetscInt);
431: fmultinplace multinplace[2] = {BVMultInPlace_Vecs_ME, BVMultInPlace_Vecs_Alloc};
432: BV_VECS *ctx = (BV_VECS*)bv->data;
434: PetscFunctionBegin;
435: ctx->vmip = vmip;
436: bv->ops->multinplace = multinplace[vmip];
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
440: static PetscErrorCode BVSetFromOptions_Vecs(BV bv,PetscOptionItems *PetscOptionsObject)
441: {
442: BV_VECS *ctx = (BV_VECS*)bv->data;
444: PetscFunctionBegin;
445: PetscOptionsHeadBegin(PetscOptionsObject,"BV Vecs Options");
447: PetscCall(PetscOptionsRangeInt("-bv_vecs_vmip","Version of BVMultInPlace operation","",ctx->vmip,&ctx->vmip,NULL,0,1));
448: PetscCall(BVVecsSetVmip(bv,ctx->vmip));
450: PetscOptionsHeadEnd();
451: PetscFunctionReturn(PETSC_SUCCESS);
452: }
454: PetscErrorCode BVView_Vecs(BV bv,PetscViewer viewer)
455: {
456: BV_VECS *ctx = (BV_VECS*)bv->data;
457: PetscInt j;
458: PetscViewerFormat format;
459: PetscBool isascii,ismatlab=PETSC_FALSE;
460: const char *bvname,*name;
462: PetscFunctionBegin;
463: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
464: if (isascii) {
465: PetscCall(PetscViewerGetFormat(viewer,&format));
466: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
467: if (format == PETSC_VIEWER_ASCII_MATLAB) ismatlab = PETSC_TRUE;
468: }
469: if (ismatlab) {
470: PetscCall(PetscObjectGetName((PetscObject)bv,&bvname));
471: PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[];\n",bvname));
472: }
473: for (j=bv->nc;j<bv->nc+bv->m;j++) {
474: PetscCall(VecView(ctx->V[j],viewer));
475: if (ismatlab) {
476: PetscCall(PetscObjectGetName((PetscObject)ctx->V[j],&name));
477: PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[%s,%s];clear %s\n",bvname,bvname,name,name));
478: }
479: }
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: static PetscErrorCode BVDestroy_Vecs(BV bv)
484: {
485: BV_VECS *ctx = (BV_VECS*)bv->data;
487: PetscFunctionBegin;
488: if (!bv->issplit) PetscCall(VecDestroyVecs(bv->nc+bv->m,&ctx->V));
489: PetscCall(PetscFree(bv->data));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: static PetscErrorCode BVDuplicate_Vecs(BV V,BV W)
494: {
495: BV_VECS *ctx = (BV_VECS*)V->data;
497: PetscFunctionBegin;
498: PetscCall(BVVecsSetVmip(W,ctx->vmip));
499: PetscFunctionReturn(PETSC_SUCCESS);
500: }
502: SLEPC_EXTERN PetscErrorCode BVCreate_Vecs(BV bv)
503: {
504: BV_VECS *ctx;
505: PetscInt j,lsplit;
506: PetscBool isgpu;
507: char str[50];
508: BV parent;
509: Vec *Vpar,z;
511: PetscFunctionBegin;
512: PetscCall(PetscNew(&ctx));
513: bv->data = (void*)ctx;
515: if (PetscUnlikely(bv->issplit)) {
516: /* split BV: share the Vecs of the parent BV */
517: parent = bv->splitparent;
518: lsplit = parent->lsplit;
519: Vpar = ((BV_VECS*)parent->data)->V;
520: ctx->V = (bv->issplit==1)? Vpar: Vpar+lsplit;
521: } else {
522: /* regular BV: create array of Vecs to store the BV columns */
523: PetscCall(BVCreateVec(bv,&z));
524: PetscCall(VecDuplicateVecs(z,bv->m,&ctx->V));
525: PetscCall(VecDestroy(&z));
526: if (((PetscObject)bv)->name) {
527: for (j=0;j<bv->m;j++) {
528: PetscCall(PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j));
529: PetscCall(PetscObjectSetName((PetscObject)ctx->V[j],str));
530: }
531: }
532: }
533: if (!bv->ld) bv->ld = bv->n;
534: PetscCheck(bv->ld==bv->n,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"BVVECS does not support a user-defined leading dimension");
536: if (PetscUnlikely(bv->Acreate)) {
537: for (j=0;j<bv->m;j++) PetscCall(MatGetColumnVector(bv->Acreate,ctx->V[j],j));
538: PetscCall(MatDestroy(&bv->Acreate));
539: }
541: /* Default version of BVMultInPlace */
542: PetscCall(PetscStrcmpAny(bv->vtype,&isgpu,VECSEQCUDA,VECMPICUDA,""));
543: ctx->vmip = isgpu? 1: 0;
545: /* Default BVMatMult method */
546: bv->vmm = BV_MATMULT_VECS;
548: /* Deferred call to setfromoptions */
549: if (bv->defersfo) {
550: PetscObjectOptionsBegin((PetscObject)bv);
551: PetscCall(BVSetFromOptions_Vecs(bv,PetscOptionsObject));
552: PetscOptionsEnd();
553: }
554: PetscCall(BVVecsSetVmip(bv,ctx->vmip));
556: bv->ops->mult = BVMult_Vecs;
557: bv->ops->multvec = BVMultVec_Vecs;
558: bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Vecs;
559: bv->ops->dot = BVDot_Vecs;
560: bv->ops->dotvec = BVDotVec_Vecs;
561: bv->ops->dotvec_begin = BVDotVec_Begin_Vecs;
562: bv->ops->dotvec_end = BVDotVec_End_Vecs;
563: bv->ops->scale = BVScale_Vecs;
564: bv->ops->norm = BVNorm_Vecs;
565: bv->ops->norm_begin = BVNorm_Begin_Vecs;
566: bv->ops->norm_end = BVNorm_End_Vecs;
567: bv->ops->normalize = BVNormalize_Vecs;
568: bv->ops->matmult = BVMatMult_Vecs;
569: bv->ops->copy = BVCopy_Vecs;
570: bv->ops->copycolumn = BVCopyColumn_Vecs;
571: bv->ops->resize = BVResize_Vecs;
572: bv->ops->getcolumn = BVGetColumn_Vecs;
573: bv->ops->restorecolumn = BVRestoreColumn_Vecs;
574: bv->ops->getarray = BVGetArray_Vecs;
575: bv->ops->restorearray = BVRestoreArray_Vecs;
576: bv->ops->getarrayread = BVGetArrayRead_Vecs;
577: bv->ops->restorearrayread = BVRestoreArrayRead_Vecs;
578: bv->ops->getmat = BVGetMat_Default;
579: bv->ops->restoremat = BVRestoreMat_Default;
580: bv->ops->destroy = BVDestroy_Vecs;
581: bv->ops->duplicate = BVDuplicate_Vecs;
582: bv->ops->setfromoptions = BVSetFromOptions_Vecs;
583: bv->ops->view = BVView_Vecs;
584: PetscFunctionReturn(PETSC_SUCCESS);
585: }