Actual source code: svec.c
slepc-3.17.2 2022-08-09
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 a single Vec
12: */
14: #include <slepc/private/bvimpl.h>
15: #include "svec.h"
17: PetscErrorCode BVMult_Svec(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
18: {
19: BV_SVEC *y = (BV_SVEC*)Y->data,*x = (BV_SVEC*)X->data;
20: const PetscScalar *px,*q;
21: PetscScalar *py;
22: PetscInt ldq;
24: VecGetArrayRead(x->v,&px);
25: VecGetArray(y->v,&py);
26: if (Q) {
27: MatDenseGetLDA(Q,&ldq);
28: MatDenseGetArrayRead(Q,&q);
29: 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);
30: MatDenseRestoreArrayRead(Q,&q);
31: } else 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);
32: VecRestoreArrayRead(x->v,&px);
33: VecRestoreArray(y->v,&py);
34: PetscFunctionReturn(0);
35: }
37: PetscErrorCode BVMultVec_Svec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
38: {
39: BV_SVEC *x = (BV_SVEC*)X->data;
40: PetscScalar *px,*py,*qq=q;
42: VecGetArray(x->v,&px);
43: VecGetArray(y,&py);
44: if (!q) VecGetArray(X->buffer,&qq);
45: BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,px+(X->nc+X->l)*X->n,qq,beta,py);
46: if (!q) VecRestoreArray(X->buffer,&qq);
47: VecRestoreArray(x->v,&px);
48: VecRestoreArray(y,&py);
49: PetscFunctionReturn(0);
50: }
52: PetscErrorCode BVMultInPlace_Svec(BV V,Mat Q,PetscInt s,PetscInt e)
53: {
54: BV_SVEC *ctx = (BV_SVEC*)V->data;
55: PetscScalar *pv;
56: const PetscScalar *q;
57: PetscInt ldq;
59: MatDenseGetLDA(Q,&ldq);
60: VecGetArray(ctx->v,&pv);
61: MatDenseGetArrayRead(Q,&q);
62: 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);
63: MatDenseRestoreArrayRead(Q,&q);
64: VecRestoreArray(ctx->v,&pv);
65: PetscFunctionReturn(0);
66: }
68: PetscErrorCode BVMultInPlaceHermitianTranspose_Svec(BV V,Mat Q,PetscInt s,PetscInt e)
69: {
70: BV_SVEC *ctx = (BV_SVEC*)V->data;
71: PetscScalar *pv;
72: const PetscScalar *q;
73: PetscInt ldq;
75: MatDenseGetLDA(Q,&ldq);
76: VecGetArray(ctx->v,&pv);
77: MatDenseGetArrayRead(Q,&q);
78: 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);
79: MatDenseRestoreArrayRead(Q,&q);
80: VecRestoreArray(ctx->v,&pv);
81: PetscFunctionReturn(0);
82: }
84: PetscErrorCode BVDot_Svec(BV X,BV Y,Mat M)
85: {
86: BV_SVEC *x = (BV_SVEC*)X->data,*y = (BV_SVEC*)Y->data;
87: const PetscScalar *px,*py;
88: PetscScalar *m;
89: PetscInt ldm;
91: MatDenseGetLDA(M,&ldm);
92: VecGetArrayRead(x->v,&px);
93: VecGetArrayRead(y->v,&py);
94: MatDenseGetArray(M,&m);
95: 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);
96: MatDenseRestoreArray(M,&m);
97: VecRestoreArrayRead(x->v,&px);
98: VecRestoreArrayRead(y->v,&py);
99: PetscFunctionReturn(0);
100: }
102: PetscErrorCode BVDotVec_Svec(BV X,Vec y,PetscScalar *q)
103: {
104: BV_SVEC *x = (BV_SVEC*)X->data;
105: const PetscScalar *px,*py;
106: PetscScalar *qq=q;
107: Vec z = y;
109: if (PetscUnlikely(X->matrix)) {
110: BV_IPMatMult(X,y);
111: z = X->Bx;
112: }
113: VecGetArrayRead(x->v,&px);
114: VecGetArrayRead(z,&py);
115: if (!q) VecGetArray(X->buffer,&qq);
116: BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,qq,x->mpi);
117: if (!q) VecRestoreArray(X->buffer,&qq);
118: VecRestoreArrayRead(z,&py);
119: VecRestoreArrayRead(x->v,&px);
120: PetscFunctionReturn(0);
121: }
123: PetscErrorCode BVDotVec_Local_Svec(BV X,Vec y,PetscScalar *m)
124: {
125: BV_SVEC *x = (BV_SVEC*)X->data;
126: PetscScalar *px,*py;
127: Vec z = y;
129: if (PetscUnlikely(X->matrix)) {
130: BV_IPMatMult(X,y);
131: z = X->Bx;
132: }
133: VecGetArray(x->v,&px);
134: VecGetArray(z,&py);
135: BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,m,PETSC_FALSE);
136: VecRestoreArray(z,&py);
137: VecRestoreArray(x->v,&px);
138: PetscFunctionReturn(0);
139: }
141: PetscErrorCode BVScale_Svec(BV bv,PetscInt j,PetscScalar alpha)
142: {
143: BV_SVEC *ctx = (BV_SVEC*)bv->data;
144: PetscScalar *array;
146: VecGetArray(ctx->v,&array);
147: if (PetscUnlikely(j<0)) BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->n,array+(bv->nc+bv->l)*bv->n,alpha);
148: else BVScale_BLAS_Private(bv,bv->n,array+(bv->nc+j)*bv->n,alpha);
149: VecRestoreArray(ctx->v,&array);
150: PetscFunctionReturn(0);
151: }
153: PetscErrorCode BVNorm_Svec(BV bv,PetscInt j,NormType type,PetscReal *val)
154: {
155: BV_SVEC *ctx = (BV_SVEC*)bv->data;
156: PetscScalar *array;
158: VecGetArray(ctx->v,&array);
159: if (PetscUnlikely(j<0)) BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,ctx->mpi);
160: else BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,ctx->mpi);
161: VecRestoreArray(ctx->v,&array);
162: PetscFunctionReturn(0);
163: }
165: PetscErrorCode BVNorm_Local_Svec(BV bv,PetscInt j,NormType type,PetscReal *val)
166: {
167: BV_SVEC *ctx = (BV_SVEC*)bv->data;
168: PetscScalar *array;
170: VecGetArray(ctx->v,&array);
171: if (PetscUnlikely(j<0)) BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,PETSC_FALSE);
172: else BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,PETSC_FALSE);
173: VecRestoreArray(ctx->v,&array);
174: PetscFunctionReturn(0);
175: }
177: PetscErrorCode BVNormalize_Svec(BV bv,PetscScalar *eigi)
178: {
179: BV_SVEC *ctx = (BV_SVEC*)bv->data;
180: PetscScalar *array,*wi=NULL;
182: VecGetArray(ctx->v,&array);
183: if (eigi) wi = eigi+bv->l;
184: BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,wi,ctx->mpi);
185: VecRestoreArray(ctx->v,&array);
186: PetscFunctionReturn(0);
187: }
189: PetscErrorCode BVMatMult_Svec(BV V,Mat A,BV W)
190: {
191: PetscInt j;
192: Mat Vmat,Wmat;
193: Vec vv,ww;
195: if (V->vmm) {
196: BVGetMat(V,&Vmat);
197: BVGetMat(W,&Wmat);
198: MatProductCreateWithMat(A,Vmat,NULL,Wmat);
199: MatProductSetType(Wmat,MATPRODUCT_AB);
200: MatProductSetFromOptions(Wmat);
201: MatProductSymbolic(Wmat);
202: MatProductNumeric(Wmat);
203: MatProductClear(Wmat);
204: BVRestoreMat(V,&Vmat);
205: BVRestoreMat(W,&Wmat);
206: } else {
207: for (j=0;j<V->k-V->l;j++) {
208: BVGetColumn(V,V->l+j,&vv);
209: BVGetColumn(W,W->l+j,&ww);
210: MatMult(A,vv,ww);
211: BVRestoreColumn(V,V->l+j,&vv);
212: BVRestoreColumn(W,W->l+j,&ww);
213: }
214: }
215: PetscFunctionReturn(0);
216: }
218: PetscErrorCode BVCopy_Svec(BV V,BV W)
219: {
220: BV_SVEC *v = (BV_SVEC*)V->data,*w = (BV_SVEC*)W->data;
221: PetscScalar *pv,*pw,*pvc,*pwc;
223: VecGetArray(v->v,&pv);
224: VecGetArray(w->v,&pw);
225: pvc = pv+(V->nc+V->l)*V->n;
226: pwc = pw+(W->nc+W->l)*W->n;
227: PetscArraycpy(pwc,pvc,(V->k-V->l)*V->n);
228: VecRestoreArray(v->v,&pv);
229: VecRestoreArray(w->v,&pw);
230: PetscFunctionReturn(0);
231: }
233: PetscErrorCode BVCopyColumn_Svec(BV V,PetscInt j,PetscInt i)
234: {
235: BV_SVEC *v = (BV_SVEC*)V->data;
236: PetscScalar *pv;
238: VecGetArray(v->v,&pv);
239: PetscArraycpy(pv+(V->nc+i)*V->n,pv+(V->nc+j)*V->n,V->n);
240: VecRestoreArray(v->v,&pv);
241: PetscFunctionReturn(0);
242: }
244: PetscErrorCode BVResize_Svec(BV bv,PetscInt m,PetscBool copy)
245: {
246: BV_SVEC *ctx = (BV_SVEC*)bv->data;
247: PetscScalar *pnew;
248: const PetscScalar *pv;
249: PetscInt bs;
250: Vec vnew;
251: char str[50];
253: VecGetBlockSize(bv->t,&bs);
254: VecCreate(PetscObjectComm((PetscObject)bv->t),&vnew);
255: VecSetType(vnew,((PetscObject)bv->t)->type_name);
256: VecSetSizes(vnew,m*bv->n,PETSC_DECIDE);
257: VecSetBlockSize(vnew,bs);
258: PetscLogObjectParent((PetscObject)bv,(PetscObject)vnew);
259: if (((PetscObject)bv)->name) {
260: PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);
261: PetscObjectSetName((PetscObject)vnew,str);
262: }
263: if (copy) {
264: VecGetArrayRead(ctx->v,&pv);
265: VecGetArray(vnew,&pnew);
266: PetscArraycpy(pnew,pv,PetscMin(m,bv->m)*bv->n);
267: VecRestoreArrayRead(ctx->v,&pv);
268: VecRestoreArray(vnew,&pnew);
269: }
270: VecDestroy(&ctx->v);
271: ctx->v = vnew;
272: PetscFunctionReturn(0);
273: }
275: PetscErrorCode BVGetColumn_Svec(BV bv,PetscInt j,Vec *v)
276: {
277: BV_SVEC *ctx = (BV_SVEC*)bv->data;
278: PetscScalar *pv;
279: PetscInt l;
281: l = BVAvailableVec;
282: VecGetArray(ctx->v,&pv);
283: VecPlaceArray(bv->cv[l],pv+(bv->nc+j)*bv->n);
284: PetscFunctionReturn(0);
285: }
287: PetscErrorCode BVRestoreColumn_Svec(BV bv,PetscInt j,Vec *v)
288: {
289: BV_SVEC *ctx = (BV_SVEC*)bv->data;
290: PetscInt l;
292: l = (j==bv->ci[0])? 0: 1;
293: VecResetArray(bv->cv[l]);
294: VecRestoreArray(ctx->v,NULL);
295: PetscFunctionReturn(0);
296: }
298: PetscErrorCode BVGetArray_Svec(BV bv,PetscScalar **a)
299: {
300: BV_SVEC *ctx = (BV_SVEC*)bv->data;
302: VecGetArray(ctx->v,a);
303: PetscFunctionReturn(0);
304: }
306: PetscErrorCode BVRestoreArray_Svec(BV bv,PetscScalar **a)
307: {
308: BV_SVEC *ctx = (BV_SVEC*)bv->data;
310: VecRestoreArray(ctx->v,a);
311: PetscFunctionReturn(0);
312: }
314: PetscErrorCode BVGetArrayRead_Svec(BV bv,const PetscScalar **a)
315: {
316: BV_SVEC *ctx = (BV_SVEC*)bv->data;
318: VecGetArrayRead(ctx->v,a);
319: PetscFunctionReturn(0);
320: }
322: PetscErrorCode BVRestoreArrayRead_Svec(BV bv,const PetscScalar **a)
323: {
324: BV_SVEC *ctx = (BV_SVEC*)bv->data;
326: VecRestoreArrayRead(ctx->v,a);
327: PetscFunctionReturn(0);
328: }
330: PetscErrorCode BVView_Svec(BV bv,PetscViewer viewer)
331: {
332: BV_SVEC *ctx = (BV_SVEC*)bv->data;
333: PetscViewerFormat format;
334: PetscBool isascii;
335: const char *bvname,*name;
337: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
338: if (isascii) {
339: PetscViewerGetFormat(viewer,&format);
340: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0);
341: VecView(ctx->v,viewer);
342: if (format == PETSC_VIEWER_ASCII_MATLAB) {
343: PetscObjectGetName((PetscObject)bv,&bvname);
344: PetscObjectGetName((PetscObject)ctx->v,&name);
345: PetscViewerASCIIPrintf(viewer,"%s=reshape(%s,%" PetscInt_FMT ",%" PetscInt_FMT ");clear %s\n",bvname,name,bv->N,bv->nc+bv->m,name);
346: if (bv->nc) PetscViewerASCIIPrintf(viewer,"%s=%s(:,%" PetscInt_FMT ":end);\n",bvname,bvname,bv->nc+1);
347: }
348: } else VecView(ctx->v,viewer);
349: PetscFunctionReturn(0);
350: }
352: PetscErrorCode BVDestroy_Svec(BV bv)
353: {
354: BV_SVEC *ctx = (BV_SVEC*)bv->data;
356: VecDestroy(&ctx->v);
357: VecDestroy(&bv->cv[0]);
358: VecDestroy(&bv->cv[1]);
359: PetscFree(bv->data);
360: bv->cuda = PETSC_FALSE;
361: PetscFunctionReturn(0);
362: }
364: SLEPC_EXTERN PetscErrorCode BVCreate_Svec(BV bv)
365: {
366: BV_SVEC *ctx;
367: PetscInt nloc,N,bs,tglobal=0,tlocal,lsplit;
368: PetscBool seq;
369: PetscScalar *vv;
370: const PetscScalar *aa,*array,*ptr;
371: char str[50];
372: BV parent;
373: Vec vpar;
374: #if defined(PETSC_HAVE_CUDA)
375: PetscScalar *gpuarray,*gptr;
376: #endif
378: PetscNewLog(bv,&ctx);
379: bv->data = (void*)ctx;
381: PetscObjectTypeCompareAny((PetscObject)bv->t,&bv->cuda,VECSEQCUDA,VECMPICUDA,"");
382: PetscObjectTypeCompareAny((PetscObject)bv->t,&ctx->mpi,VECMPI,VECMPICUDA,"");
384: PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq);
387: VecGetLocalSize(bv->t,&nloc);
388: VecGetSize(bv->t,&N);
389: VecGetBlockSize(bv->t,&bs);
390: tlocal = bv->m*nloc;
391: PetscIntMultError(bv->m,N,&tglobal);
393: if (PetscUnlikely(bv->issplit)) {
394: /* split BV: create Vec sharing the memory of the parent BV */
395: parent = bv->splitparent;
396: lsplit = parent->lsplit;
397: vpar = ((BV_SVEC*)parent->data)->v;
398: if (bv->cuda) {
399: #if defined(PETSC_HAVE_CUDA)
400: VecCUDAGetArray(vpar,&gpuarray);
401: gptr = (bv->issplit==1)? gpuarray: gpuarray+lsplit*nloc;
402: VecCUDARestoreArray(vpar,&gpuarray);
403: if (ctx->mpi) VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv->t),bs,tlocal,bv->m*N,NULL,&ctx->v);
404: else VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv->t),bs,tlocal,NULL,&ctx->v);
405: VecCUDAPlaceArray(ctx->v,gptr);
406: #endif
407: } else {
408: VecGetArrayRead(vpar,&array);
409: ptr = (bv->issplit==1)? array: array+lsplit*nloc;
410: VecRestoreArrayRead(vpar,&array);
411: if (ctx->mpi) VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,tlocal,bv->m*N,NULL,&ctx->v);
412: else VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,tlocal,NULL,&ctx->v);
413: VecPlaceArray(ctx->v,ptr);
414: }
415: } else {
416: /* regular BV: create Vec to store the BV entries */
417: VecCreate(PetscObjectComm((PetscObject)bv->t),&ctx->v);
418: VecSetType(ctx->v,((PetscObject)bv->t)->type_name);
419: VecSetSizes(ctx->v,tlocal,tglobal);
420: VecSetBlockSize(ctx->v,bs);
421: }
422: PetscLogObjectParent((PetscObject)bv,(PetscObject)ctx->v);
423: if (((PetscObject)bv)->name) {
424: PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);
425: PetscObjectSetName((PetscObject)ctx->v,str);
426: }
428: if (PetscUnlikely(bv->Acreate)) {
429: MatDenseGetArrayRead(bv->Acreate,&aa);
430: VecGetArray(ctx->v,&vv);
431: PetscArraycpy(vv,aa,tlocal);
432: VecRestoreArray(ctx->v,&vv);
433: MatDenseRestoreArrayRead(bv->Acreate,&aa);
434: MatDestroy(&bv->Acreate);
435: }
437: VecDuplicateEmpty(bv->t,&bv->cv[0]);
438: VecDuplicateEmpty(bv->t,&bv->cv[1]);
440: if (bv->cuda) {
441: #if defined(PETSC_HAVE_CUDA)
442: bv->ops->mult = BVMult_Svec_CUDA;
443: bv->ops->multvec = BVMultVec_Svec_CUDA;
444: bv->ops->multinplace = BVMultInPlace_Svec_CUDA;
445: bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Svec_CUDA;
446: bv->ops->dot = BVDot_Svec_CUDA;
447: bv->ops->dotvec = BVDotVec_Svec_CUDA;
448: bv->ops->dotvec_local = BVDotVec_Local_Svec_CUDA;
449: bv->ops->scale = BVScale_Svec_CUDA;
450: bv->ops->matmult = BVMatMult_Svec_CUDA;
451: bv->ops->copy = BVCopy_Svec_CUDA;
452: bv->ops->copycolumn = BVCopyColumn_Svec_CUDA;
453: bv->ops->resize = BVResize_Svec_CUDA;
454: bv->ops->getcolumn = BVGetColumn_Svec_CUDA;
455: bv->ops->restorecolumn = BVRestoreColumn_Svec_CUDA;
456: bv->ops->restoresplit = BVRestoreSplit_Svec_CUDA;
457: bv->ops->getmat = BVGetMat_Svec_CUDA;
458: bv->ops->restoremat = BVRestoreMat_Svec_CUDA;
459: #endif
460: } else {
461: bv->ops->mult = BVMult_Svec;
462: bv->ops->multvec = BVMultVec_Svec;
463: bv->ops->multinplace = BVMultInPlace_Svec;
464: bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Svec;
465: bv->ops->dot = BVDot_Svec;
466: bv->ops->dotvec = BVDotVec_Svec;
467: bv->ops->dotvec_local = BVDotVec_Local_Svec;
468: bv->ops->scale = BVScale_Svec;
469: bv->ops->matmult = BVMatMult_Svec;
470: bv->ops->copy = BVCopy_Svec;
471: bv->ops->copycolumn = BVCopyColumn_Svec;
472: bv->ops->resize = BVResize_Svec;
473: bv->ops->getcolumn = BVGetColumn_Svec;
474: bv->ops->restorecolumn = BVRestoreColumn_Svec;
475: }
476: bv->ops->norm = BVNorm_Svec;
477: bv->ops->norm_local = BVNorm_Local_Svec;
478: bv->ops->normalize = BVNormalize_Svec;
479: bv->ops->getarray = BVGetArray_Svec;
480: bv->ops->restorearray = BVRestoreArray_Svec;
481: bv->ops->getarrayread = BVGetArrayRead_Svec;
482: bv->ops->restorearrayread = BVRestoreArrayRead_Svec;
483: bv->ops->destroy = BVDestroy_Svec;
484: if (!ctx->mpi) bv->ops->view = BVView_Svec;
485: PetscFunctionReturn(0);
486: }