Actual source code: vecs.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: {
26: BV_VECS *y = (BV_VECS*)Y->data,*x = (BV_VECS*)X->data;
27: PetscScalar *q,*s=NULL;
28: PetscInt i,j,ldq;
31: if (Q) {
32: MatGetSize(Q,&ldq,NULL);
33: if (alpha!=1.0) {
34: BVAllocateWork_Private(Y,X->k-X->l);
35: s = Y->work;
36: }
37: MatDenseGetArray(Q,&q);
38: for (j=Y->l;j<Y->k;j++) {
39: VecScale(y->V[Y->nc+j],beta);
40: if (alpha!=1.0) {
41: for (i=X->l;i<X->k;i++) s[i-X->l] = alpha*q[i+j*ldq];
42: } else s = 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: MatDenseRestoreArray(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: {
58: BV_VECS *x = (BV_VECS*)X->data;
59: PetscScalar *s=NULL,*qq=q;
60: PetscInt i;
63: if (alpha!=1.0) {
64: BVAllocateWork_Private(X,X->k-X->l);
65: s = X->work;
66: }
67: if (!q) { VecGetArray(X->buffer,&qq); }
68: VecScale(y,beta);
69: if (alpha!=1.0) {
70: for (i=0;i<X->k-X->l;i++) s[i] = alpha*qq[i];
71: } else s = qq;
72: VecMAXPY(y,X->k-X->l,s,x->V+X->nc+X->l);
73: if (!q) { VecRestoreArray(X->buffer,&qq); }
74: return(0);
75: }
77: /*
78: BVMultInPlace_Vecs_ME - V(:,s:e-1) = V*Q(:,s:e-1) for regular vectors.
80: Memory-efficient version, uses VecGetArray (default in CPU)
82: Writing V = [ V1 V2 V3 ] and Q(:,s:e-1) = [ Q1 Q2 Q3 ]', where V2
83: corresponds to the columns s:e-1, the computation is done as
84: V2 := V2*Q2 + V1*Q1 + V3*Q3
85: */
86: PetscErrorCode BVMultInPlace_Vecs_ME(BV V,Mat Q,PetscInt s,PetscInt e)
87: {
89: BV_VECS *ctx = (BV_VECS*)V->data;
90: PetscScalar *q;
91: PetscInt i,ldq;
94: MatGetSize(Q,&ldq,NULL);
95: MatDenseGetArray(Q,&q);
96: /* V2 := V2*Q2 */
97: 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 (s>V->l) {
101: VecMAXPY(ctx->V[V->nc+i],s-V->l,q+i*ldq+V->l,ctx->V+V->nc+V->l);
102: }
103: if (V->k>e) {
104: VecMAXPY(ctx->V[V->nc+i],V->k-e,q+i*ldq+e,ctx->V+V->nc+e);
105: }
106: }
107: MatDenseRestoreArray(Q,&q);
108: return(0);
109: }
111: /*
112: BVMultInPlace_Vecs_Alloc - V(:,s:e-1) = V*Q(:,s:e-1) for regular vectors.
114: Version that allocates (e-s) work vectors in every call (default in GPU)
115: */
116: PetscErrorCode BVMultInPlace_Vecs_Alloc(BV V,Mat Q,PetscInt s,PetscInt e)
117: {
119: BV_VECS *ctx = (BV_VECS*)V->data;
120: PetscScalar *q;
121: PetscInt i,ldq;
122: Vec *W;
125: MatGetSize(Q,&ldq,NULL);
126: MatDenseGetArray(Q,&q);
127: VecDuplicateVecs(V->t,e-s,&W);
128: for (i=s;i<e;i++) {
129: VecMAXPY(W[i-s],V->k-V->l,q+i*ldq+V->l,ctx->V+V->nc+V->l);
130: }
131: for (i=s;i<e;i++) {
132: VecCopy(W[i-s],ctx->V[V->nc+i]);
133: }
134: VecDestroyVecs(e-s,&W);
135: MatDenseRestoreArray(Q,&q);
136: return(0);
137: }
139: /*
140: BVMultInPlaceTranspose_Vecs - V(:,s:e-1) = V*Q'(:,s:e-1) for regular vectors.
141: */
142: PetscErrorCode BVMultInPlaceTranspose_Vecs(BV V,Mat Q,PetscInt s,PetscInt e)
143: {
145: BV_VECS *ctx = (BV_VECS*)V->data;
146: PetscScalar *q;
147: PetscInt i,j,ldq,n;
150: MatGetSize(Q,&ldq,&n);
151: MatDenseGetArray(Q,&q);
152: /* V2 := V2*Q2' */
153: BVMultInPlace_Vecs_Private(V,V->n,e-s,ldq,ctx->V+V->nc+s,q+s*ldq+s,PETSC_TRUE);
154: /* V2 += V1*Q1' + V3*Q3' */
155: for (i=s;i<e;i++) {
156: for (j=V->l;j<s;j++) {
157: VecAXPY(ctx->V[V->nc+i],q[i+j*ldq],ctx->V[V->nc+j]);
158: }
159: for (j=e;j<n;j++) {
160: VecAXPY(ctx->V[V->nc+i],q[i+j*ldq],ctx->V[V->nc+j]);
161: }
162: }
163: MatDenseRestoreArray(Q,&q);
164: return(0);
165: }
167: PetscErrorCode BVDot_Vecs(BV X,BV Y,Mat M)
168: {
170: BV_VECS *x = (BV_VECS*)X->data,*y = (BV_VECS*)Y->data;
171: PetscScalar *m;
172: PetscInt j,ldm;
175: MatGetSize(M,&ldm,NULL);
176: MatDenseGetArray(M,&m);
177: for (j=X->l;j<X->k;j++) {
178: VecMDot(x->V[X->nc+j],Y->k-Y->l,y->V+Y->nc+Y->l,m+j*ldm+Y->l);
179: }
180: MatDenseRestoreArray(M,&m);
181: return(0);
182: }
184: PetscErrorCode BVDotVec_Vecs(BV X,Vec y,PetscScalar *q)
185: {
187: BV_VECS *x = (BV_VECS*)X->data;
188: Vec z = y;
189: PetscScalar *qq=q;
192: if (X->matrix) {
193: BV_IPMatMult(X,y);
194: z = X->Bx;
195: }
196: if (!q) { VecGetArray(X->buffer,&qq); }
197: VecMDot(z,X->k-X->l,x->V+X->nc+X->l,qq);
198: if (!q) { VecRestoreArray(X->buffer,&qq); }
199: return(0);
200: }
202: PetscErrorCode BVDotVec_Begin_Vecs(BV X,Vec y,PetscScalar *m)
203: {
205: BV_VECS *x = (BV_VECS*)X->data;
206: Vec z = y;
209: if (X->matrix) {
210: BV_IPMatMult(X,y);
211: z = X->Bx;
212: }
213: VecMDotBegin(z,X->k-X->l,x->V+X->nc+X->l,m);
214: return(0);
215: }
217: PetscErrorCode BVDotVec_End_Vecs(BV X,Vec y,PetscScalar *m)
218: {
220: BV_VECS *x = (BV_VECS*)X->data;
223: VecMDotEnd(y,X->k-X->l,x->V+X->nc+X->l,m);
224: return(0);
225: }
227: PetscErrorCode BVScale_Vecs(BV bv,PetscInt j,PetscScalar alpha)
228: {
230: PetscInt i;
231: BV_VECS *ctx = (BV_VECS*)bv->data;
234: if (j<0) {
235: for (i=bv->l;i<bv->k;i++) {
236: VecScale(ctx->V[bv->nc+i],alpha);
237: }
238: } else {
239: VecScale(ctx->V[bv->nc+j],alpha);
240: }
241: return(0);
242: }
244: PetscErrorCode BVNorm_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
245: {
247: PetscInt i;
248: PetscReal nrm;
249: BV_VECS *ctx = (BV_VECS*)bv->data;
252: if (j<0) {
253: if (type==NORM_FROBENIUS) {
254: *val = 0.0;
255: for (i=bv->l;i<bv->k;i++) {
256: VecNorm(ctx->V[bv->nc+i],NORM_2,&nrm);
257: *val += nrm*nrm;
258: }
259: *val = PetscSqrtReal(*val);
260: } else SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
261: } else {
262: VecNorm(ctx->V[bv->nc+j],type,val);
263: }
264: return(0);
265: }
267: PetscErrorCode BVNorm_Begin_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
268: {
270: BV_VECS *ctx = (BV_VECS*)bv->data;
273: if (j<0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
274: else {
275: VecNormBegin(ctx->V[bv->nc+j],type,val);
276: }
277: return(0);
278: }
280: PetscErrorCode BVNorm_End_Vecs(BV bv,PetscInt j,NormType type,PetscReal *val)
281: {
283: BV_VECS *ctx = (BV_VECS*)bv->data;
286: if (j<0) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Requested norm not implemented in BVVECS");
287: else {
288: VecNormEnd(ctx->V[bv->nc+j],type,val);
289: }
290: return(0);
291: }
293: PetscErrorCode BVMatMult_Vecs(BV V,Mat A,BV W)
294: {
296: BV_VECS *v = (BV_VECS*)V->data,*w = (BV_VECS*)W->data;
297: PetscInt j;
298: PetscBool flg;
299: Mat Vmat,Wmat;
302: MatHasOperation(A,MATOP_MAT_MULT,&flg);
303: if (V->vmm && flg) {
304: BVGetMat(V,&Vmat);
305: BVGetMat(W,&Wmat);
306: MatMatMult(A,Vmat,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Wmat);
307: BVRestoreMat(V,&Vmat);
308: BVRestoreMat(W,&Wmat);
309: } else {
310: for (j=0;j<V->k-V->l;j++) {
311: MatMult(A,v->V[V->nc+V->l+j],w->V[W->nc+W->l+j]);
312: }
313: }
314: return(0);
315: }
317: PetscErrorCode BVCopy_Vecs(BV V,BV W)
318: {
320: BV_VECS *v = (BV_VECS*)V->data,*w = (BV_VECS*)W->data;
321: PetscInt j;
324: for (j=0;j<V->k-V->l;j++) {
325: VecCopy(v->V[V->nc+V->l+j],w->V[W->nc+W->l+j]);
326: }
327: return(0);
328: }
330: PetscErrorCode BVCopyColumn_Vecs(BV V,PetscInt j,PetscInt i)
331: {
333: BV_VECS *v = (BV_VECS*)V->data;
336: VecCopy(v->V[V->nc+j],v->V[V->nc+i]);
337: return(0);
338: }
340: PetscErrorCode BVResize_Vecs(BV bv,PetscInt m,PetscBool copy)
341: {
343: BV_VECS *ctx = (BV_VECS*)bv->data;
344: Vec *newV;
345: PetscInt j,lsplit;
346: char str[50];
347: BV parent;
350: if (bv->issplit==2) {
351: parent = bv->splitparent;
352: lsplit = parent->lsplit;
353: ctx->V = ((BV_VECS*)parent->data)->V+lsplit;
354: } else if (!bv->issplit) {
355: VecDuplicateVecs(bv->t,m,&newV);
356: PetscLogObjectParents(bv,m,newV);
357: if (((PetscObject)bv)->name) {
358: for (j=0;j<m;j++) {
359: PetscSNPrintf(str,50,"%s_%D",((PetscObject)bv)->name,j);
360: PetscObjectSetName((PetscObject)newV[j],str);
361: }
362: }
363: if (copy) {
364: for (j=0;j<PetscMin(m,bv->m);j++) {
365: VecCopy(ctx->V[j],newV[j]);
366: }
367: }
368: VecDestroyVecs(bv->m,&ctx->V);
369: ctx->V = newV;
370: }
371: return(0);
372: }
374: PetscErrorCode BVGetColumn_Vecs(BV bv,PetscInt j,Vec *v)
375: {
376: BV_VECS *ctx = (BV_VECS*)bv->data;
377: PetscInt l;
380: l = BVAvailableVec;
381: bv->cv[l] = ctx->V[bv->nc+j];
382: return(0);
383: }
385: PetscErrorCode BVGetArray_Vecs(BV bv,PetscScalar **a)
386: {
387: PetscErrorCode ierr;
388: BV_VECS *ctx = (BV_VECS*)bv->data;
389: PetscInt j;
390: const PetscScalar *p;
393: PetscMalloc1((bv->nc+bv->m)*bv->n,a);
394: for (j=0;j<bv->nc+bv->m;j++) {
395: VecGetArrayRead(ctx->V[j],&p);
396: PetscMemcpy(*a+j*bv->n,p,bv->n*sizeof(PetscScalar));
397: VecRestoreArrayRead(ctx->V[j],&p);
398: }
399: return(0);
400: }
402: PetscErrorCode BVRestoreArray_Vecs(BV bv,PetscScalar **a)
403: {
405: BV_VECS *ctx = (BV_VECS*)bv->data;
406: PetscInt j;
407: PetscScalar *p;
410: for (j=0;j<bv->nc+bv->m;j++) {
411: VecGetArray(ctx->V[j],&p);
412: PetscMemcpy(p,*a+j*bv->n,bv->n*sizeof(PetscScalar));
413: VecRestoreArray(ctx->V[j],&p);
414: }
415: PetscFree(*a);
416: return(0);
417: }
419: PetscErrorCode BVGetArrayRead_Vecs(BV bv,const PetscScalar **a)
420: {
421: PetscErrorCode ierr;
422: BV_VECS *ctx = (BV_VECS*)bv->data;
423: PetscInt j;
424: const PetscScalar *p;
427: PetscMalloc1((bv->nc+bv->m)*bv->n,(PetscScalar**)a);
428: for (j=0;j<bv->nc+bv->m;j++) {
429: VecGetArrayRead(ctx->V[j],&p);
430: PetscMemcpy((PetscScalar*)*a+j*bv->n,p,bv->n*sizeof(PetscScalar));
431: VecRestoreArrayRead(ctx->V[j],&p);
432: }
433: return(0);
434: }
436: PetscErrorCode BVRestoreArrayRead_Vecs(BV bv,const PetscScalar **a)
437: {
441: PetscFree(*a);
442: return(0);
443: }
445: /*
446: Sets the value of vmip flag and resets ops->multinplace accordingly
447: */
448: PETSC_STATIC_INLINE PetscErrorCode BVVecsSetVmip(BV bv,PetscInt vmip)
449: {
450: typedef PetscErrorCode (*fmultinplace)(BV,Mat,PetscInt,PetscInt);
451: fmultinplace multinplace[2] = {BVMultInPlace_Vecs_ME, BVMultInPlace_Vecs_Alloc};
452: BV_VECS *ctx = (BV_VECS*)bv->data;
455: ctx->vmip = vmip;
456: bv->ops->multinplace = multinplace[vmip];
457: return(0);
458: }
460: PetscErrorCode BVSetFromOptions_Vecs(PetscOptionItems *PetscOptionsObject,BV bv)
461: {
463: BV_VECS *ctx = (BV_VECS*)bv->data;
466: PetscOptionsHead(PetscOptionsObject,"BV Vecs Options");
468: PetscOptionsInt("-bv_vecs_vmip","Version of BVMultInPlace operation","",ctx->vmip,&ctx->vmip,NULL);
469: if (ctx->vmip<0 || ctx->vmip>1) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_ARG_OUTOFRANGE,"Wrong version of BVMultInPlace");
470: BVVecsSetVmip(bv,ctx->vmip);
472: PetscOptionsTail();
473: return(0);
474: }
476: PetscErrorCode BVView_Vecs(BV bv,PetscViewer viewer)
477: {
478: PetscErrorCode ierr;
479: BV_VECS *ctx = (BV_VECS*)bv->data;
480: PetscInt j;
481: PetscViewerFormat format;
482: PetscBool isascii,ismatlab=PETSC_FALSE;
483: const char *bvname,*name;
486: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
487: if (isascii) {
488: PetscViewerGetFormat(viewer,&format);
489: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
490: if (format == PETSC_VIEWER_ASCII_MATLAB) ismatlab = PETSC_TRUE;
491: }
492: if (ismatlab) {
493: PetscObjectGetName((PetscObject)bv,&bvname);
494: PetscViewerASCIIPrintf(viewer,"%s=[];\n",bvname);
495: }
496: for (j=bv->nc;j<bv->nc+bv->m;j++) {
497: VecView(ctx->V[j],viewer);
498: if (ismatlab) {
499: PetscObjectGetName((PetscObject)ctx->V[j],&name);
500: PetscViewerASCIIPrintf(viewer,"%s=[%s,%s];clear %s\n",bvname,bvname,name,name);
501: }
502: }
503: return(0);
504: }
506: PetscErrorCode BVDestroy_Vecs(BV bv)
507: {
509: BV_VECS *ctx = (BV_VECS*)bv->data;
512: if (!bv->issplit) { VecDestroyVecs(bv->nc+bv->m,&ctx->V); }
513: PetscFree(bv->data);
514: return(0);
515: }
517: PetscErrorCode BVDuplicate_Vecs(BV V,BV W)
518: {
520: BV_VECS *ctx = (BV_VECS*)V->data;
523: BVVecsSetVmip(W,ctx->vmip);
524: return(0);
525: }
527: SLEPC_EXTERN PetscErrorCode BVCreate_Vecs(BV bv)
528: {
530: BV_VECS *ctx;
531: PetscInt j,lsplit;
532: PetscBool isgpu;
533: char str[50];
534: BV parent;
535: Vec *Vpar;
538: PetscNewLog(bv,&ctx);
539: bv->data = (void*)ctx;
541: if (bv->issplit) {
542: /* split BV: share the Vecs of the parent BV */
543: parent = bv->splitparent;
544: lsplit = parent->lsplit;
545: Vpar = ((BV_VECS*)parent->data)->V;
546: ctx->V = (bv->issplit==1)? Vpar: Vpar+lsplit;
547: } else {
548: /* regular BV: create array of Vecs to store the BV columns */
549: VecDuplicateVecs(bv->t,bv->m,&ctx->V);
550: PetscLogObjectParents(bv,bv->m,ctx->V);
551: if (((PetscObject)bv)->name) {
552: for (j=0;j<bv->m;j++) {
553: PetscSNPrintf(str,50,"%s_%D",((PetscObject)bv)->name,j);
554: PetscObjectSetName((PetscObject)ctx->V[j],str);
555: }
556: }
557: }
559: if (bv->Acreate) {
560: for (j=0;j<bv->m;j++) {
561: MatGetColumnVector(bv->Acreate,ctx->V[j],j);
562: }
563: MatDestroy(&bv->Acreate);
564: }
566: /* Default version of BVMultInPlace */
567: PetscObjectTypeCompareAny((PetscObject)bv->t,&isgpu,VECSEQCUDA,VECMPICUDA,"");
568: ctx->vmip = isgpu? 1: 0;
570: /* Default BVMatMult method */
571: bv->vmm = BV_MATMULT_VECS;
573: /* Deferred call to setfromoptions */
574: if (bv->defersfo) {
575: PetscObjectOptionsBegin((PetscObject)bv);
576: BVSetFromOptions_Vecs(PetscOptionsObject,bv);
577: PetscOptionsEnd();
578: }
579: BVVecsSetVmip(bv,ctx->vmip);
581: bv->ops->mult = BVMult_Vecs;
582: bv->ops->multvec = BVMultVec_Vecs;
583: bv->ops->multinplacetrans = BVMultInPlaceTranspose_Vecs;
584: bv->ops->dot = BVDot_Vecs;
585: bv->ops->dotvec = BVDotVec_Vecs;
586: bv->ops->dotvec_begin = BVDotVec_Begin_Vecs;
587: bv->ops->dotvec_end = BVDotVec_End_Vecs;
588: bv->ops->scale = BVScale_Vecs;
589: bv->ops->norm = BVNorm_Vecs;
590: bv->ops->norm_begin = BVNorm_Begin_Vecs;
591: bv->ops->norm_end = BVNorm_End_Vecs;
592: bv->ops->matmult = BVMatMult_Vecs;
593: bv->ops->copy = BVCopy_Vecs;
594: bv->ops->copycolumn = BVCopyColumn_Vecs;
595: bv->ops->resize = BVResize_Vecs;
596: bv->ops->getcolumn = BVGetColumn_Vecs;
597: bv->ops->getarray = BVGetArray_Vecs;
598: bv->ops->restorearray = BVRestoreArray_Vecs;
599: bv->ops->getarrayread = BVGetArrayRead_Vecs;
600: bv->ops->restorearrayread = BVRestoreArrayRead_Vecs;
601: bv->ops->destroy = BVDestroy_Vecs;
602: bv->ops->duplicate = BVDuplicate_Vecs;
603: bv->ops->setfromoptions = BVSetFromOptions_Vecs;
604: bv->ops->view = BVView_Vecs;
605: return(0);
606: }