Actual source code: vecs.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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: }