Actual source code: bvmat.c
slepc-3.14.2 2021-02-01
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2020, 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 with a dense Mat
12: */
14: #include <slepc/private/bvimpl.h>
16: typedef struct {
17: Mat A;
18: PetscBool mpi;
19: } BV_MAT;
21: PetscErrorCode BVMult_Mat(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
22: {
24: BV_MAT *y = (BV_MAT*)Y->data,*x = (BV_MAT*)X->data;
25: PetscScalar *px,*py,*q;
26: PetscInt ldq;
29: MatDenseGetArray(x->A,&px);
30: MatDenseGetArray(y->A,&py);
31: if (Q) {
32: MatGetSize(Q,&ldq,NULL);
33: MatDenseGetArray(Q,&q);
34: 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);
35: MatDenseRestoreArray(Q,&q);
36: } else {
37: 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);
38: }
39: MatDenseRestoreArray(x->A,&px);
40: MatDenseRestoreArray(y->A,&py);
41: return(0);
42: }
44: PetscErrorCode BVMultVec_Mat(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
45: {
47: BV_MAT *x = (BV_MAT*)X->data;
48: PetscScalar *px,*py,*qq=q;
51: MatDenseGetArray(x->A,&px);
52: VecGetArray(y,&py);
53: if (!q) { VecGetArray(X->buffer,&qq); }
54: BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,px+(X->nc+X->l)*X->n,qq,beta,py);
55: if (!q) { VecRestoreArray(X->buffer,&qq); }
56: MatDenseRestoreArray(x->A,&px);
57: VecRestoreArray(y,&py);
58: return(0);
59: }
61: PetscErrorCode BVMultInPlace_Mat(BV V,Mat Q,PetscInt s,PetscInt e)
62: {
64: BV_MAT *ctx = (BV_MAT*)V->data;
65: PetscScalar *pv,*q;
66: PetscInt ldq;
69: MatGetSize(Q,&ldq,NULL);
70: MatDenseGetArray(ctx->A,&pv);
71: MatDenseGetArray(Q,&q);
72: 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);
73: MatDenseRestoreArray(Q,&q);
74: MatDenseRestoreArray(ctx->A,&pv);
75: return(0);
76: }
78: PetscErrorCode BVMultInPlaceTranspose_Mat(BV V,Mat Q,PetscInt s,PetscInt e)
79: {
81: BV_MAT *ctx = (BV_MAT*)V->data;
82: PetscScalar *pv,*q;
83: PetscInt ldq;
86: MatGetSize(Q,&ldq,NULL);
87: MatDenseGetArray(ctx->A,&pv);
88: MatDenseGetArray(Q,&q);
89: 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);
90: MatDenseRestoreArray(Q,&q);
91: MatDenseRestoreArray(ctx->A,&pv);
92: return(0);
93: }
95: PetscErrorCode BVDot_Mat(BV X,BV Y,Mat M)
96: {
98: BV_MAT *x = (BV_MAT*)X->data,*y = (BV_MAT*)Y->data;
99: PetscScalar *px,*py,*m;
100: PetscInt ldm;
103: MatGetSize(M,&ldm,NULL);
104: MatDenseGetArray(x->A,&px);
105: MatDenseGetArray(y->A,&py);
106: MatDenseGetArray(M,&m);
107: 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);
108: MatDenseRestoreArray(M,&m);
109: MatDenseRestoreArray(x->A,&px);
110: MatDenseRestoreArray(y->A,&py);
111: return(0);
112: }
114: PetscErrorCode BVDotVec_Mat(BV X,Vec y,PetscScalar *q)
115: {
116: PetscErrorCode ierr;
117: BV_MAT *x = (BV_MAT*)X->data;
118: PetscScalar *px,*qq=q;
119: const PetscScalar *py;
120: Vec z = y;
123: if (X->matrix) {
124: BV_IPMatMult(X,y);
125: z = X->Bx;
126: }
127: MatDenseGetArray(x->A,&px);
128: VecGetArrayRead(z,&py);
129: if (!q) { VecGetArray(X->buffer,&qq); }
130: BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,qq,x->mpi);
131: if (!q) { VecRestoreArray(X->buffer,&qq); }
132: VecRestoreArrayRead(z,&py);
133: MatDenseRestoreArray(x->A,&px);
134: return(0);
135: }
137: PetscErrorCode BVDotVec_Local_Mat(BV X,Vec y,PetscScalar *m)
138: {
140: BV_MAT *x = (BV_MAT*)X->data;
141: PetscScalar *px,*py;
142: Vec z = y;
145: if (X->matrix) {
146: BV_IPMatMult(X,y);
147: z = X->Bx;
148: }
149: MatDenseGetArray(x->A,&px);
150: VecGetArray(z,&py);
151: BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,m,PETSC_FALSE);
152: VecRestoreArray(z,&py);
153: MatDenseRestoreArray(x->A,&px);
154: return(0);
155: }
157: PetscErrorCode BVScale_Mat(BV bv,PetscInt j,PetscScalar alpha)
158: {
160: BV_MAT *ctx = (BV_MAT*)bv->data;
161: PetscScalar *array;
164: MatDenseGetArray(ctx->A,&array);
165: if (j<0) {
166: BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->n,array+(bv->nc+bv->l)*bv->n,alpha);
167: } else {
168: BVScale_BLAS_Private(bv,bv->n,array+(bv->nc+j)*bv->n,alpha);
169: }
170: MatDenseRestoreArray(ctx->A,&array);
171: return(0);
172: }
174: PetscErrorCode BVNorm_Mat(BV bv,PetscInt j,NormType type,PetscReal *val)
175: {
177: BV_MAT *ctx = (BV_MAT*)bv->data;
178: PetscScalar *array;
181: MatDenseGetArray(ctx->A,&array);
182: if (j<0) {
183: BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,ctx->mpi);
184: } else {
185: BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,ctx->mpi);
186: }
187: MatDenseRestoreArray(ctx->A,&array);
188: return(0);
189: }
191: PetscErrorCode BVNorm_Local_Mat(BV bv,PetscInt j,NormType type,PetscReal *val)
192: {
194: BV_MAT *ctx = (BV_MAT*)bv->data;
195: PetscScalar *array;
198: MatDenseGetArray(ctx->A,&array);
199: if (j<0) {
200: BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,PETSC_FALSE);
201: } else {
202: BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,PETSC_FALSE);
203: }
204: MatDenseRestoreArray(ctx->A,&array);
205: return(0);
206: }
208: PetscErrorCode BVNormalize_Mat(BV bv,PetscScalar *eigi)
209: {
211: BV_MAT *ctx = (BV_MAT*)bv->data;
212: PetscScalar *array,*wi=NULL;
215: MatDenseGetArray(ctx->A,&array);
216: if (eigi) wi = eigi+bv->l;
217: BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,wi,ctx->mpi);
218: MatDenseRestoreArray(ctx->A,&array);
219: return(0);
220: }
222: PetscErrorCode BVMatMult_Mat(BV V,Mat A,BV W)
223: {
225: PetscInt j;
226: Mat Vmat,Wmat;
227: Vec vv,ww;
230: if (V->vmm) {
231: BVGetMat(V,&Vmat);
232: BVGetMat(W,&Wmat);
233: MatProductCreateWithMat(A,Vmat,NULL,Wmat);
234: MatProductSetType(Wmat,MATPRODUCT_AB);
235: MatProductSetFromOptions(Wmat);
236: MatProductSymbolic(Wmat);
237: MatProductNumeric(Wmat);
238: MatProductClear(Wmat);
239: BVRestoreMat(V,&Vmat);
240: BVRestoreMat(W,&Wmat);
241: } else {
242: for (j=0;j<V->k-V->l;j++) {
243: BVGetColumn(V,V->l+j,&vv);
244: BVGetColumn(W,W->l+j,&ww);
245: MatMult(A,vv,ww);
246: BVRestoreColumn(V,V->l+j,&vv);
247: BVRestoreColumn(W,W->l+j,&ww);
248: }
249: }
250: return(0);
251: }
253: PetscErrorCode BVCopy_Mat(BV V,BV W)
254: {
256: BV_MAT *v = (BV_MAT*)V->data,*w = (BV_MAT*)W->data;
257: PetscScalar *pv,*pw,*pvc,*pwc;
260: MatDenseGetArray(v->A,&pv);
261: MatDenseGetArray(w->A,&pw);
262: pvc = pv+(V->nc+V->l)*V->n;
263: pwc = pw+(W->nc+W->l)*W->n;
264: PetscArraycpy(pwc,pvc,(V->k-V->l)*V->n);
265: MatDenseRestoreArray(v->A,&pv);
266: MatDenseRestoreArray(w->A,&pw);
267: return(0);
268: }
270: PetscErrorCode BVCopyColumn_Mat(BV V,PetscInt j,PetscInt i)
271: {
273: BV_MAT *v = (BV_MAT*)V->data;
274: PetscScalar *pv;
277: MatDenseGetArray(v->A,&pv);
278: PetscArraycpy(pv+(V->nc+i)*V->n,pv+(V->nc+j)*V->n,V->n);
279: MatDenseRestoreArray(v->A,&pv);
280: return(0);
281: }
283: PetscErrorCode BVResize_Mat(BV bv,PetscInt m,PetscBool copy)
284: {
286: BV_MAT *ctx = (BV_MAT*)bv->data;
287: PetscScalar *pA,*pnew;
288: Mat A;
289: char str[50];
292: MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,PETSC_DECIDE,m,NULL,&A);
293: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
294: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
295: PetscLogObjectParent((PetscObject)bv,(PetscObject)A);
296: if (((PetscObject)bv)->name) {
297: PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);
298: PetscObjectSetName((PetscObject)A,str);
299: }
300: if (copy) {
301: MatDenseGetArray(ctx->A,&pA);
302: MatDenseGetArray(A,&pnew);
303: PetscArraycpy(pnew,pA,PetscMin(m,bv->m)*bv->n);
304: MatDenseRestoreArray(ctx->A,&pA);
305: MatDenseRestoreArray(A,&pnew);
306: }
307: MatDestroy(&ctx->A);
308: ctx->A = A;
309: return(0);
310: }
312: PetscErrorCode BVGetColumn_Mat(BV bv,PetscInt j,Vec *v)
313: {
315: BV_MAT *ctx = (BV_MAT*)bv->data;
316: PetscScalar *pA;
317: PetscInt l;
320: l = BVAvailableVec;
321: MatDenseGetArray(ctx->A,&pA);
322: VecPlaceArray(bv->cv[l],pA+(bv->nc+j)*bv->n);
323: return(0);
324: }
326: PetscErrorCode BVRestoreColumn_Mat(BV bv,PetscInt j,Vec *v)
327: {
329: BV_MAT *ctx = (BV_MAT*)bv->data;
330: PetscScalar *pA;
331: PetscInt l;
334: l = (j==bv->ci[0])? 0: 1;
335: VecResetArray(bv->cv[l]);
336: MatDenseRestoreArray(ctx->A,&pA);
337: return(0);
338: }
340: PetscErrorCode BVGetArray_Mat(BV bv,PetscScalar **a)
341: {
343: BV_MAT *ctx = (BV_MAT*)bv->data;
346: MatDenseGetArray(ctx->A,a);
347: return(0);
348: }
350: PetscErrorCode BVRestoreArray_Mat(BV bv,PetscScalar **a)
351: {
353: BV_MAT *ctx = (BV_MAT*)bv->data;
356: if (a) { MatDenseRestoreArray(ctx->A,a); }
357: return(0);
358: }
360: PetscErrorCode BVGetArrayRead_Mat(BV bv,const PetscScalar **a)
361: {
363: BV_MAT *ctx = (BV_MAT*)bv->data;
366: MatDenseGetArray(ctx->A,(PetscScalar**)a);
367: return(0);
368: }
370: PetscErrorCode BVRestoreArrayRead_Mat(BV bv,const PetscScalar **a)
371: {
373: BV_MAT *ctx = (BV_MAT*)bv->data;
376: if (a) { MatDenseRestoreArray(ctx->A,(PetscScalar**)a); }
377: return(0);
378: }
380: PetscErrorCode BVView_Mat(BV bv,PetscViewer viewer)
381: {
382: PetscErrorCode ierr;
383: BV_MAT *ctx = (BV_MAT*)bv->data;
384: PetscViewerFormat format;
385: PetscBool isascii;
386: const char *bvname,*name;
389: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
390: if (isascii) {
391: PetscViewerGetFormat(viewer,&format);
392: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
393: MatView(ctx->A,viewer);
394: if (format == PETSC_VIEWER_ASCII_MATLAB) {
395: PetscObjectGetName((PetscObject)bv,&bvname);
396: PetscObjectGetName((PetscObject)ctx->A,&name);
397: PetscViewerASCIIPrintf(viewer,"%s=%s;clear %s\n",bvname,name,name);
398: if (bv->nc) {
399: PetscViewerASCIIPrintf(viewer,"%s=%s(:,%D:end);\n",bvname,bvname,bv->nc+1);
400: }
401: }
402: } else {
403: MatView(ctx->A,viewer);
404: }
405: return(0);
406: }
408: PetscErrorCode BVDestroy_Mat(BV bv)
409: {
411: BV_MAT *ctx = (BV_MAT*)bv->data;
414: MatDestroy(&ctx->A);
415: VecDestroy(&bv->cv[0]);
416: VecDestroy(&bv->cv[1]);
417: PetscFree(bv->data);
418: return(0);
419: }
421: SLEPC_EXTERN PetscErrorCode BVCreate_Mat(BV bv)
422: {
424: BV_MAT *ctx;
425: PetscInt nloc,bs,lsplit;
426: PetscBool seq;
427: char str[50];
428: PetscScalar *array,*ptr;
429: BV parent;
432: PetscNewLog(bv,&ctx);
433: bv->data = (void*)ctx;
435: PetscObjectTypeCompare((PetscObject)bv->t,VECMPI,&ctx->mpi);
436: if (!ctx->mpi) {
437: PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq);
438: if (!seq) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot create a BVMAT from a non-standard template vector");
439: }
441: VecGetLocalSize(bv->t,&nloc);
442: VecGetBlockSize(bv->t,&bs);
444: if (bv->issplit) {
445: /* split BV: share the memory of the parent BV */
446: parent = bv->splitparent;
447: lsplit = parent->lsplit;
448: MatDenseGetArray(((BV_MAT*)parent->data)->A,&array);
449: ptr = (bv->issplit==1)? array: array+lsplit*nloc;
450: MatDenseRestoreArray(((BV_MAT*)parent->data)->A,&array);
451: } else {
452: /* regular BV: allocate memory for the BV entries */
453: ptr = NULL;
454: }
455: MatCreateDense(PetscObjectComm((PetscObject)bv->t),nloc,PETSC_DECIDE,PETSC_DECIDE,bv->m,ptr,&ctx->A);
456: MatAssemblyBegin(ctx->A,MAT_FINAL_ASSEMBLY);
457: MatAssemblyEnd(ctx->A,MAT_FINAL_ASSEMBLY);
458: PetscLogObjectParent((PetscObject)bv,(PetscObject)ctx->A);
459: if (((PetscObject)bv)->name) {
460: PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name);
461: PetscObjectSetName((PetscObject)ctx->A,str);
462: }
464: if (bv->Acreate) {
465: MatCopy(bv->Acreate,ctx->A,SAME_NONZERO_PATTERN);
466: MatDestroy(&bv->Acreate);
467: }
469: VecDuplicateEmpty(bv->t,&bv->cv[0]);
470: VecDuplicateEmpty(bv->t,&bv->cv[1]);
472: bv->ops->mult = BVMult_Mat;
473: bv->ops->multvec = BVMultVec_Mat;
474: bv->ops->multinplace = BVMultInPlace_Mat;
475: bv->ops->multinplacetrans = BVMultInPlaceTranspose_Mat;
476: bv->ops->dot = BVDot_Mat;
477: bv->ops->dotvec = BVDotVec_Mat;
478: bv->ops->dotvec_local = BVDotVec_Local_Mat;
479: bv->ops->scale = BVScale_Mat;
480: bv->ops->norm = BVNorm_Mat;
481: bv->ops->norm_local = BVNorm_Local_Mat;
482: bv->ops->normalize = BVNormalize_Mat;
483: bv->ops->matmult = BVMatMult_Mat;
484: bv->ops->copy = BVCopy_Mat;
485: bv->ops->copycolumn = BVCopyColumn_Mat;
486: bv->ops->resize = BVResize_Mat;
487: bv->ops->getcolumn = BVGetColumn_Mat;
488: bv->ops->restorecolumn = BVRestoreColumn_Mat;
489: bv->ops->getarray = BVGetArray_Mat;
490: bv->ops->restorearray = BVRestoreArray_Mat;
491: bv->ops->getarrayread = BVGetArrayRead_Mat;
492: bv->ops->restorearrayread = BVRestoreArrayRead_Mat;
493: bv->ops->destroy = BVDestroy_Mat;
494: if (!ctx->mpi) bv->ops->view = BVView_Mat;
495: return(0);
496: }