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