Actual source code: svec.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 a single Vec
12: */
14: #include <slepc/private/bvimpl.h>
15: #include "./svecimpl.h"
17: PetscErrorCode BVMult_Svec(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
18: {
19: PetscErrorCode ierr;
20: BV_SVEC *y = (BV_SVEC*)Y->data,*x = (BV_SVEC*)X->data;
21: const PetscScalar *px;
22: PetscScalar *py,*q;
23: PetscInt ldq;
26: VecGetArrayRead(x->v,&px);
27: VecGetArray(y->v,&py);
28: if (Q) {
29: MatGetSize(Q,&ldq,NULL);
30: MatDenseGetArray(Q,&q);
31: BVMult_BLAS_Private(Y,Y->n,Y->k-Y->l,X->k-X->l,ldq,alpha,px+(X->nc+X->l)*X->n,q+Y->l*ldq+X->l,beta,py+(Y->nc+Y->l)*Y->n);
32: MatDenseRestoreArray(Q,&q);
33: } else {
34: BVAXPY_BLAS_Private(Y,Y->n,Y->k-Y->l,alpha,px+(X->nc+X->l)*X->n,beta,py+(Y->nc+Y->l)*Y->n);
35: }
36: VecRestoreArrayRead(x->v,&px);
37: VecRestoreArray(y->v,&py);
38: return(0);
39: }
41: PetscErrorCode BVMultVec_Svec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
42: {
44: BV_SVEC *x = (BV_SVEC*)X->data;
45: PetscScalar *px,*py,*qq=q;
48: VecGetArray(x->v,&px);
49: VecGetArray(y,&py);
50: if (!q) { VecGetArray(X->buffer,&qq); }
51: BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,px+(X->nc+X->l)*X->n,qq,beta,py);
52: if (!q) { VecRestoreArray(X->buffer,&qq); }
53: VecRestoreArray(x->v,&px);
54: VecRestoreArray(y,&py);
55: return(0);
56: }
58: PetscErrorCode BVMultInPlace_Svec(BV V,Mat Q,PetscInt s,PetscInt e)
59: {
61: BV_SVEC *ctx = (BV_SVEC*)V->data;
62: PetscScalar *pv,*q;
63: PetscInt ldq;
66: MatGetSize(Q,&ldq,NULL);
67: VecGetArray(ctx->v,&pv);
68: MatDenseGetArray(Q,&q);
69: BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,ldq,s-V->l,e-V->l,pv+(V->nc+V->l)*V->n,q+V->l*ldq+V->l,PETSC_FALSE);
70: MatDenseRestoreArray(Q,&q);
71: VecRestoreArray(ctx->v,&pv);
72: return(0);
73: }
75: PetscErrorCode BVMultInPlaceTranspose_Svec(BV V,Mat Q,PetscInt s,PetscInt e)
76: {
78: BV_SVEC *ctx = (BV_SVEC*)V->data;
79: PetscScalar *pv,*q;
80: PetscInt ldq;
83: MatGetSize(Q,&ldq,NULL);
84: VecGetArray(ctx->v,&pv);
85: MatDenseGetArray(Q,&q);
86: BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,ldq,s-V->l,e-V->l,pv+(V->nc+V->l)*V->n,q+V->l*ldq+V->l,PETSC_TRUE);
87: MatDenseRestoreArray(Q,&q);
88: VecRestoreArray(ctx->v,&pv);
89: return(0);
90: }
92: PetscErrorCode BVDot_Svec(BV X,BV Y,Mat M)
93: {
94: PetscErrorCode ierr;
95: BV_SVEC *x = (BV_SVEC*)X->data,*y = (BV_SVEC*)Y->data;
96: const PetscScalar *px,*py;
97: PetscScalar *m;
98: PetscInt ldm;
101: MatGetSize(M,&ldm,NULL);
102: VecGetArrayRead(x->v,&px);
103: VecGetArrayRead(y->v,&py);
104: MatDenseGetArray(M,&m);
105: BVDot_BLAS_Private(X,Y->k-Y->l,X->k-X->l,X->n,ldm,py+(Y->nc+Y->l)*Y->n,px+(X->nc+X->l)*X->n,m+X->l*ldm+Y->l,x->mpi);
106: MatDenseRestoreArray(M,&m);
107: VecRestoreArrayRead(x->v,&px);
108: VecRestoreArrayRead(y->v,&py);
109: return(0);
110: }
112: PetscErrorCode BVDotVec_Svec(BV X,Vec y,PetscScalar *q)
113: {
114: PetscErrorCode ierr;
115: BV_SVEC *x = (BV_SVEC*)X->data;
116: const PetscScalar *px,*py;
117: PetscScalar *qq=q;
118: Vec z = y;
121: if (X->matrix) {
122: BV_IPMatMult(X,y);
123: z = X->Bx;
124: }
125: VecGetArrayRead(x->v,&px);
126: VecGetArrayRead(z,&py);
127: if (!q) { VecGetArray(X->buffer,&qq); }
128: BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,qq,x->mpi);
129: if (!q) { VecRestoreArray(X->buffer,&qq); }
130: VecRestoreArrayRead(z,&py);
131: VecRestoreArrayRead(x->v,&px);
132: return(0);
133: }
135: PetscErrorCode BVDotVec_Local_Svec(BV X,Vec y,PetscScalar *m)
136: {
138: BV_SVEC *x = (BV_SVEC*)X->data;
139: PetscScalar *px,*py;
140: Vec z = y;
143: if (X->matrix) {
144: BV_IPMatMult(X,y);
145: z = X->Bx;
146: }
147: VecGetArray(x->v,&px);
148: VecGetArray(z,&py);
149: BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,m,PETSC_FALSE);
150: VecRestoreArray(z,&py);
151: VecRestoreArray(x->v,&px);
152: return(0);
153: }
155: PetscErrorCode BVScale_Svec(BV bv,PetscInt j,PetscScalar alpha)
156: {
158: BV_SVEC *ctx = (BV_SVEC*)bv->data;
159: PetscScalar *array;
162: VecGetArray(ctx->v,&array);
163: if (j<0) {
164: BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->n,array+(bv->nc+bv->l)*bv->n,alpha);
165: } else {
166: BVScale_BLAS_Private(bv,bv->n,array+(bv->nc+j)*bv->n,alpha);
167: }
168: VecRestoreArray(ctx->v,&array);
169: return(0);
170: }
172: PetscErrorCode BVNorm_Svec(BV bv,PetscInt j,NormType type,PetscReal *val)
173: {
175: BV_SVEC *ctx = (BV_SVEC*)bv->data;
176: PetscScalar *array;
179: VecGetArray(ctx->v,&array);
180: if (j<0) {
181: BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,ctx->mpi);
182: } else {
183: BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,ctx->mpi);
184: }
185: VecRestoreArray(ctx->v,&array);
186: return(0);
187: }
189: PetscErrorCode BVNorm_Local_Svec(BV bv,PetscInt j,NormType type,PetscReal *val)
190: {
192: BV_SVEC *ctx = (BV_SVEC*)bv->data;
193: PetscScalar *array;
196: VecGetArray(ctx->v,&array);
197: if (j<0) {
198: BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,PETSC_FALSE);
199: } else {
200: BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,PETSC_FALSE);
201: }
202: VecRestoreArray(ctx->v,&array);
203: return(0);
204: }
206: PetscErrorCode BVMatMult_Svec(BV V,Mat A,BV W)
207: {
208: PetscErrorCode ierr;
209: BV_SVEC *v = (BV_SVEC*)V->data,*w = (BV_SVEC*)W->data;
210: PetscScalar *pv,*pw;
211: PetscInt j;
212: PetscBool flg;
213: Mat Vmat,Wmat,aux;
214: PetscObjectState Astate;
217: MatHasOperation(A,MATOP_MAT_MULT,&flg);
218: if (V->vmm && flg) {
219: BVGetMat(V,&Vmat);
220: BVGetMat(W,&Wmat);
221: PetscObjectStateGet((PetscObject)A,&Astate);
222: if (V->Amult==A && V->Amultstate==Astate) {
223: MatMatMult(A,Vmat,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Wmat);
224: } else {
225: MatMatMult(A,Vmat,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&aux);
226: MatCopy(aux,Wmat,SAME_NONZERO_PATTERN);
227: MatDestroy(&aux);
228: V->Amult = A;
229: V->Amultstate = Astate;
230: }
231: BVRestoreMat(V,&Vmat);
232: BVRestoreMat(W,&Wmat);
233: } else {
234: VecGetArray(v->v,&pv);
235: VecGetArray(w->v,&pw);
236: for (j=0;j<V->k-V->l;j++) {
237: VecPlaceArray(V->cv[1],pv+(V->nc+V->l+j)*V->n);
238: VecPlaceArray(W->cv[1],pw+(W->nc+W->l+j)*W->n);
239: MatMult(A,V->cv[1],W->cv[1]);
240: VecResetArray(V->cv[1]);
241: VecResetArray(W->cv[1]);
242: }
243: VecRestoreArray(v->v,&pv);
244: VecRestoreArray(w->v,&pw);
245: }
246: return(0);
247: }
249: PetscErrorCode BVCopy_Svec(BV V,BV W)
250: {
252: BV_SVEC *v = (BV_SVEC*)V->data,*w = (BV_SVEC*)W->data;
253: PetscScalar *pv,*pw,*pvc,*pwc;
256: VecGetArray(v->v,&pv);
257: VecGetArray(w->v,&pw);
258: pvc = pv+(V->nc+V->l)*V->n;
259: pwc = pw+(W->nc+W->l)*W->n;
260: PetscMemcpy(pwc,pvc,(V->k-V->l)*V->n*sizeof(PetscScalar));
261: VecRestoreArray(v->v,&pv);
262: VecRestoreArray(w->v,&pw);
263: return(0);
264: }
266: PetscErrorCode BVCopyColumn_Svec(BV V,PetscInt j,PetscInt i)
267: {
269: BV_SVEC *v = (BV_SVEC*)V->data;
270: PetscScalar *pv;
273: VecGetArray(v->v,&pv);
274: PetscMemcpy(pv+(V->nc+i)*V->n,pv+(V->nc+j)*V->n,V->n*sizeof(PetscScalar));
275: VecRestoreArray(v->v,&pv);
276: return(0);
277: }
279: PetscErrorCode BVResize_Svec(BV bv,PetscInt m,PetscBool copy)
280: {
281: PetscErrorCode ierr;
282: BV_SVEC *ctx = (BV_SVEC*)bv->data;
283: PetscScalar *pnew;
284: const PetscScalar *pv;
285: PetscInt bs,lsplit;
286: Vec vnew,vpar;
287: char str[50];
288: BV parent;
291: if (bv->issplit==2) {
292: parent = bv->splitparent;
293: lsplit = parent->lsplit;
294: vpar = ((BV_SVEC*)parent->data)->v;
295: VecGetArrayRead(vpar,&pv);
296: VecResetArray(ctx->v);
297: VecPlaceArray(ctx->v,pv+lsplit*bv->n);
298: VecRestoreArrayRead(vpar,&pv);
299: } else if (!bv->issplit) {
300: VecGetBlockSize(bv->t,&bs);
301: VecCreate(PetscObjectComm((PetscObject)bv->t),&vnew);
302: VecSetType(vnew,((PetscObject)bv->t)->type_name);
303: VecSetSizes(vnew,m*bv->n,PETSC_DECIDE);
304: VecSetBlockSize(vnew,bs);
305: PetscLogObjectParent((PetscObject)bv,(PetscObject)vnew);
306: if (((PetscObject)bv)->name) {
307: PetscSNPrintf(str,50,"%s_0",((PetscObject)bv)->name);
308: PetscObjectSetName((PetscObject)vnew,str);
309: }
310: if (copy) {
311: VecGetArrayRead(ctx->v,&pv);
312: VecGetArray(vnew,&pnew);
313: PetscMemcpy(pnew,pv,PetscMin(m,bv->m)*bv->n*sizeof(PetscScalar));
314: VecRestoreArrayRead(ctx->v,&pv);
315: VecRestoreArray(vnew,&pnew);
316: }
317: VecDestroy(&ctx->v);
318: ctx->v = vnew;
319: }
320: return(0);
321: }
323: PetscErrorCode BVGetColumn_Svec(BV bv,PetscInt j,Vec *v)
324: {
326: BV_SVEC *ctx = (BV_SVEC*)bv->data;
327: PetscScalar *pv;
328: PetscInt l;
331: l = BVAvailableVec;
332: VecGetArray(ctx->v,&pv);
333: VecPlaceArray(bv->cv[l],pv+(bv->nc+j)*bv->n);
334: return(0);
335: }
337: PetscErrorCode BVRestoreColumn_Svec(BV bv,PetscInt j,Vec *v)
338: {
340: BV_SVEC *ctx = (BV_SVEC*)bv->data;
341: PetscInt l;
344: l = (j==bv->ci[0])? 0: 1;
345: VecResetArray(bv->cv[l]);
346: VecRestoreArray(ctx->v,NULL);
347: return(0);
348: }
350: PetscErrorCode BVGetArray_Svec(BV bv,PetscScalar **a)
351: {
353: BV_SVEC *ctx = (BV_SVEC*)bv->data;
356: VecGetArray(ctx->v,a);
357: return(0);
358: }
360: PetscErrorCode BVRestoreArray_Svec(BV bv,PetscScalar **a)
361: {
363: BV_SVEC *ctx = (BV_SVEC*)bv->data;
366: VecRestoreArray(ctx->v,a);
367: return(0);
368: }
370: PetscErrorCode BVGetArrayRead_Svec(BV bv,const PetscScalar **a)
371: {
373: BV_SVEC *ctx = (BV_SVEC*)bv->data;
376: VecGetArrayRead(ctx->v,a);
377: return(0);
378: }
380: PetscErrorCode BVRestoreArrayRead_Svec(BV bv,const PetscScalar **a)
381: {
383: BV_SVEC *ctx = (BV_SVEC*)bv->data;
386: VecRestoreArrayRead(ctx->v,a);
387: return(0);
388: }
390: PetscErrorCode BVView_Svec(BV bv,PetscViewer viewer)
391: {
392: PetscErrorCode ierr;
393: BV_SVEC *ctx = (BV_SVEC*)bv->data;
394: PetscViewerFormat format;
395: PetscBool isascii;
396: const char *bvname,*name;
399: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
400: if (isascii) {
401: PetscViewerGetFormat(viewer,&format);
402: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
403: VecView(ctx->v,viewer);
404: if (format == PETSC_VIEWER_ASCII_MATLAB) {
405: PetscObjectGetName((PetscObject)bv,&bvname);
406: PetscObjectGetName((PetscObject)ctx->v,&name);
407: PetscViewerASCIIPrintf(viewer,"%s=reshape(%s,%D,%D);clear %s\n",bvname,name,bv->N,bv->nc+bv->m,name);
408: if (bv->nc) {
409: PetscViewerASCIIPrintf(viewer,"%s=%s(:,%D:end);\n",bvname,bvname,bv->nc+1);
410: }
411: }
412: } else {
413: VecView(ctx->v,viewer);
414: }
415: return(0);
416: }
418: PetscErrorCode BVDestroy_Svec(BV bv)
419: {
421: BV_SVEC *ctx = (BV_SVEC*)bv->data;
424: VecDestroy(&ctx->v);
425: VecDestroy(&bv->cv[0]);
426: VecDestroy(&bv->cv[1]);
427: PetscFree(bv->data);
428: bv->cuda = PETSC_FALSE;
429: return(0);
430: }
432: SLEPC_EXTERN PetscErrorCode BVCreate_Svec(BV bv)
433: {
434: PetscErrorCode ierr;
435: BV_SVEC *ctx;
436: PetscInt nloc,N,bs,tglobal,tlocal,lsplit;
437: PetscBool seq;
438: PetscScalar *aa,*vv;
439: const PetscScalar *array,*ptr;
440: char str[50];
441: BV parent;
442: Vec vpar;
443: #if defined(PETSC_HAVE_CUDA)
444: PetscScalar *gpuarray,*gptr;
445: #endif
448: PetscNewLog(bv,&ctx);
449: bv->data = (void*)ctx;
451: PetscObjectTypeCompareAny((PetscObject)bv->t,&bv->cuda,VECSEQCUDA,VECMPICUDA,"");
452: PetscObjectTypeCompareAny((PetscObject)bv->t,&ctx->mpi,VECMPI,VECMPICUDA,"");
454: PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq);
455: if (!seq && !ctx->mpi && !bv->cuda) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"BVSVEC does not support the type of the provided template vector");
457: VecGetLocalSize(bv->t,&nloc);
458: VecGetSize(bv->t,&N);
459: VecGetBlockSize(bv->t,&bs);
460: tlocal = bv->m*nloc;
461: tglobal = bv->m*N;
462: if (tglobal<0) SETERRQ2(PetscObjectComm((PetscObject)bv),1,"The product %D times %D overflows the size of PetscInt; consider reducing the number of columns, or use BVVECS instead",bv->m,N);
464: if (bv->issplit) {
465: /* split BV: create Vec sharing the memory of the parent BV */
466: parent = bv->splitparent;
467: lsplit = parent->lsplit;
468: vpar = ((BV_SVEC*)parent->data)->v;
469: if (bv->cuda) {
470: #if defined(PETSC_HAVE_CUDA)
471: VecCUDAGetArrayReadWrite(vpar,&gpuarray);
472: gptr = (bv->issplit==1)? gpuarray: gpuarray+lsplit*nloc;
473: VecCUDARestoreArrayReadWrite(vpar,&gpuarray);
474: if (ctx->mpi) {
475: VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv->t),bs,tlocal,bv->m*N,NULL,&ctx->v);
476: } else {
477: VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv->t),bs,tlocal,NULL,&ctx->v);
478: }
479: VecCUDAPlaceArray(ctx->v,gptr);
480: #endif
481: } else {
482: VecGetArrayRead(vpar,&array);
483: ptr = (bv->issplit==1)? array: array+lsplit*nloc;
484: VecRestoreArrayRead(vpar,&array);
485: if (ctx->mpi) {
486: VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,tlocal,bv->m*N,NULL,&ctx->v);
487: } else {
488: VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,tlocal,NULL,&ctx->v);
489: }
490: VecPlaceArray(ctx->v,ptr);
491: }
492: } else {
493: /* regular BV: create Vec to store the BV entries */
494: VecCreate(PetscObjectComm((PetscObject)bv->t),&ctx->v);
495: VecSetType(ctx->v,((PetscObject)bv->t)->type_name);
496: VecSetSizes(ctx->v,tlocal,tglobal);
497: VecSetBlockSize(ctx->v,bs);
498: }
499: PetscLogObjectParent((PetscObject)bv,(PetscObject)ctx->v);
500: if (((PetscObject)bv)->name) {
501: PetscSNPrintf(str,50,"%s_0",((PetscObject)bv)->name);
502: PetscObjectSetName((PetscObject)ctx->v,str);
503: }
505: if (bv->Acreate) {
506: MatDenseGetArray(bv->Acreate,&aa);
507: VecGetArray(ctx->v,&vv);
508: PetscMemcpy(vv,aa,tlocal*sizeof(PetscScalar));
509: VecRestoreArray(ctx->v,&vv);
510: MatDenseRestoreArray(bv->Acreate,&aa);
511: MatDestroy(&bv->Acreate);
512: }
514: VecDuplicateEmpty(bv->t,&bv->cv[0]);
515: VecDuplicateEmpty(bv->t,&bv->cv[1]);
517: if (bv->cuda) {
518: #if defined(PETSC_HAVE_CUDA)
519: bv->ops->mult = BVMult_Svec_CUDA;
520: bv->ops->multvec = BVMultVec_Svec_CUDA;
521: bv->ops->multinplace = BVMultInPlace_Svec_CUDA;
522: bv->ops->multinplacetrans = BVMultInPlaceTranspose_Svec_CUDA;
523: bv->ops->dot = BVDot_Svec_CUDA;
524: bv->ops->dotvec = BVDotVec_Svec_CUDA;
525: bv->ops->dotvec_local = BVDotVec_Local_Svec_CUDA;
526: bv->ops->scale = BVScale_Svec_CUDA;
527: bv->ops->matmult = BVMatMult_Svec_CUDA;
528: bv->ops->copy = BVCopy_Svec_CUDA;
529: bv->ops->copycolumn = BVCopyColumn_Svec_CUDA;
530: bv->ops->resize = BVResize_Svec_CUDA;
531: bv->ops->getcolumn = BVGetColumn_Svec_CUDA;
532: bv->ops->restorecolumn = BVRestoreColumn_Svec_CUDA;
533: bv->ops->restoresplit = BVRestoreSplit_Svec_CUDA;
534: #endif
535: } else {
536: bv->ops->mult = BVMult_Svec;
537: bv->ops->multvec = BVMultVec_Svec;
538: bv->ops->multinplace = BVMultInPlace_Svec;
539: bv->ops->multinplacetrans = BVMultInPlaceTranspose_Svec;
540: bv->ops->dot = BVDot_Svec;
541: bv->ops->dotvec = BVDotVec_Svec;
542: bv->ops->dotvec_local = BVDotVec_Local_Svec;
543: bv->ops->scale = BVScale_Svec;
544: bv->ops->matmult = BVMatMult_Svec;
545: bv->ops->copy = BVCopy_Svec;
546: bv->ops->copycolumn = BVCopyColumn_Svec;
547: bv->ops->resize = BVResize_Svec;
548: bv->ops->getcolumn = BVGetColumn_Svec;
549: bv->ops->restorecolumn = BVRestoreColumn_Svec;
550: }
551: bv->ops->norm = BVNorm_Svec;
552: bv->ops->norm_local = BVNorm_Local_Svec;
553: bv->ops->getarray = BVGetArray_Svec;
554: bv->ops->restorearray = BVRestoreArray_Svec;
555: bv->ops->getarrayread = BVGetArrayRead_Svec;
556: bv->ops->restorearrayread = BVRestoreArrayRead_Svec;
557: bv->ops->destroy = BVDestroy_Svec;
558: if (!ctx->mpi) bv->ops->view = BVView_Svec;
559: return(0);
560: }