Actual source code: bvmat.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 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 BVMatMult_Mat(BV V,Mat A,BV W)
209: {
211: BV_MAT *v = (BV_MAT*)V->data,*w = (BV_MAT*)W->data;
212: PetscScalar *pv,*pw;
213: PetscInt j;
214: PetscBool flg;
215: Mat Vmat,Wmat;
218: MatHasOperation(A,MATOP_MAT_MULT,&flg);
219: if (V->vmm && flg) {
220: BVGetMat(V,&Vmat);
221: BVGetMat(W,&Wmat);
222: MatMatMult(A,Vmat,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Wmat);
223: BVRestoreMat(V,&Vmat);
224: BVRestoreMat(W,&Wmat);
225: } else {
226: MatDenseGetArray(v->A,&pv);
227: MatDenseGetArray(w->A,&pw);
228: for (j=0;j<V->k-V->l;j++) {
229: VecPlaceArray(V->cv[1],pv+(V->nc+V->l+j)*V->n);
230: VecPlaceArray(W->cv[1],pw+(W->nc+W->l+j)*W->n);
231: MatMult(A,V->cv[1],W->cv[1]);
232: VecResetArray(V->cv[1]);
233: VecResetArray(W->cv[1]);
234: }
235: MatDenseRestoreArray(v->A,&pv);
236: MatDenseRestoreArray(w->A,&pw);
237: }
238: return(0);
239: }
241: PetscErrorCode BVCopy_Mat(BV V,BV W)
242: {
244: BV_MAT *v = (BV_MAT*)V->data,*w = (BV_MAT*)W->data;
245: PetscScalar *pv,*pw,*pvc,*pwc;
248: MatDenseGetArray(v->A,&pv);
249: MatDenseGetArray(w->A,&pw);
250: pvc = pv+(V->nc+V->l)*V->n;
251: pwc = pw+(W->nc+W->l)*W->n;
252: PetscMemcpy(pwc,pvc,(V->k-V->l)*V->n*sizeof(PetscScalar));
253: MatDenseRestoreArray(v->A,&pv);
254: MatDenseRestoreArray(w->A,&pw);
255: return(0);
256: }
258: PetscErrorCode BVCopyColumn_Mat(BV V,PetscInt j,PetscInt i)
259: {
261: BV_MAT *v = (BV_MAT*)V->data;
262: PetscScalar *pv;
265: MatDenseGetArray(v->A,&pv);
266: PetscMemcpy(pv+(V->nc+i)*V->n,pv+(V->nc+j)*V->n,V->n*sizeof(PetscScalar));
267: MatDenseRestoreArray(v->A,&pv);
268: return(0);
269: }
271: PetscErrorCode BVResize_Mat(BV bv,PetscInt m,PetscBool copy)
272: {
274: BV_MAT *ctx = (BV_MAT*)bv->data;
275: PetscScalar *pA,*pnew;
276: PetscInt lsplit;
277: Mat A;
278: char str[50];
279: PetscScalar *array;
280: BV parent;
283: if (bv->issplit==2) {
284: parent = bv->splitparent;
285: lsplit = parent->lsplit;
286: MatDenseGetArray(((BV_MAT*)parent->data)->A,&array);
287: MatDensePlaceArray(ctx->A,array+lsplit*bv->n);
288: MatDenseRestoreArray(((BV_MAT*)parent->data)->A,&array);
289: } else if (!bv->issplit) {
290: MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,PETSC_DECIDE,m,NULL,&A);
291: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
292: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
293: PetscLogObjectParent((PetscObject)bv,(PetscObject)A);
294: if (((PetscObject)bv)->name) {
295: PetscSNPrintf(str,50,"%s_0",((PetscObject)bv)->name);
296: PetscObjectSetName((PetscObject)A,str);
297: }
298: if (copy) {
299: MatDenseGetArray(ctx->A,&pA);
300: MatDenseGetArray(A,&pnew);
301: PetscMemcpy(pnew,pA,PetscMin(m,bv->m)*bv->n*sizeof(PetscScalar));
302: MatDenseRestoreArray(ctx->A,&pA);
303: MatDenseRestoreArray(A,&pnew);
304: }
305: MatDestroy(&ctx->A);
306: ctx->A = A;
307: }
308: return(0);
309: }
311: PetscErrorCode BVGetColumn_Mat(BV bv,PetscInt j,Vec *v)
312: {
314: BV_MAT *ctx = (BV_MAT*)bv->data;
315: PetscScalar *pA;
316: PetscInt l;
319: l = BVAvailableVec;
320: MatDenseGetArray(ctx->A,&pA);
321: VecPlaceArray(bv->cv[l],pA+(bv->nc+j)*bv->n);
322: return(0);
323: }
325: PetscErrorCode BVRestoreColumn_Mat(BV bv,PetscInt j,Vec *v)
326: {
328: BV_MAT *ctx = (BV_MAT*)bv->data;
329: PetscScalar *pA;
330: PetscInt l;
333: l = (j==bv->ci[0])? 0: 1;
334: VecResetArray(bv->cv[l]);
335: MatDenseRestoreArray(ctx->A,&pA);
336: return(0);
337: }
339: PetscErrorCode BVGetArray_Mat(BV bv,PetscScalar **a)
340: {
342: BV_MAT *ctx = (BV_MAT*)bv->data;
345: MatDenseGetArray(ctx->A,a);
346: return(0);
347: }
349: PetscErrorCode BVRestoreArray_Mat(BV bv,PetscScalar **a)
350: {
352: BV_MAT *ctx = (BV_MAT*)bv->data;
355: if (a) { MatDenseRestoreArray(ctx->A,a); }
356: return(0);
357: }
359: PetscErrorCode BVGetArrayRead_Mat(BV bv,const PetscScalar **a)
360: {
362: BV_MAT *ctx = (BV_MAT*)bv->data;
365: MatDenseGetArray(ctx->A,(PetscScalar**)a);
366: return(0);
367: }
369: PetscErrorCode BVRestoreArrayRead_Mat(BV bv,const PetscScalar **a)
370: {
372: BV_MAT *ctx = (BV_MAT*)bv->data;
375: if (a) { MatDenseRestoreArray(ctx->A,(PetscScalar**)a); }
376: return(0);
377: }
379: PetscErrorCode BVView_Mat(BV bv,PetscViewer viewer)
380: {
381: PetscErrorCode ierr;
382: BV_MAT *ctx = (BV_MAT*)bv->data;
383: PetscViewerFormat format;
384: PetscBool isascii;
385: const char *bvname,*name;
388: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
389: if (isascii) {
390: PetscViewerGetFormat(viewer,&format);
391: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
392: MatView(ctx->A,viewer);
393: if (format == PETSC_VIEWER_ASCII_MATLAB) {
394: PetscObjectGetName((PetscObject)bv,&bvname);
395: PetscObjectGetName((PetscObject)ctx->A,&name);
396: PetscViewerASCIIPrintf(viewer,"%s=%s;clear %s\n",bvname,name,name);
397: if (bv->nc) {
398: PetscViewerASCIIPrintf(viewer,"%s=%s(:,%D:end);\n",bvname,bvname,bv->nc+1);
399: }
400: }
401: } else {
402: MatView(ctx->A,viewer);
403: }
404: return(0);
405: }
407: PetscErrorCode BVDestroy_Mat(BV bv)
408: {
410: BV_MAT *ctx = (BV_MAT*)bv->data;
413: MatDestroy(&ctx->A);
414: VecDestroy(&bv->cv[0]);
415: VecDestroy(&bv->cv[1]);
416: PetscFree(bv->data);
417: return(0);
418: }
420: SLEPC_EXTERN PetscErrorCode BVCreate_Mat(BV bv)
421: {
423: BV_MAT *ctx;
424: PetscInt nloc,bs,lsplit;
425: PetscBool seq;
426: char str[50];
427: PetscScalar *array,*ptr;
428: BV parent;
431: PetscNewLog(bv,&ctx);
432: bv->data = (void*)ctx;
434: PetscObjectTypeCompare((PetscObject)bv->t,VECMPI,&ctx->mpi);
435: if (!ctx->mpi) {
436: PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq);
437: if (!seq) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot create a BVMAT from a non-standard template vector");
438: }
440: VecGetLocalSize(bv->t,&nloc);
441: VecGetBlockSize(bv->t,&bs);
443: if (bv->issplit) {
444: /* split BV: share the memory of the parent BV */
445: parent = bv->splitparent;
446: lsplit = parent->lsplit;
447: MatDenseGetArray(((BV_MAT*)parent->data)->A,&array);
448: ptr = (bv->issplit==1)? array: array+lsplit*nloc;
449: MatDenseRestoreArray(((BV_MAT*)parent->data)->A,&array);
450: } else {
451: /* regular BV: allocate memory for the BV entries */
452: ptr = NULL;
453: }
454: MatCreateDense(PetscObjectComm((PetscObject)bv->t),nloc,PETSC_DECIDE,PETSC_DECIDE,bv->m,ptr,&ctx->A);
455: MatAssemblyBegin(ctx->A,MAT_FINAL_ASSEMBLY);
456: MatAssemblyEnd(ctx->A,MAT_FINAL_ASSEMBLY);
457: PetscLogObjectParent((PetscObject)bv,(PetscObject)ctx->A);
458: if (((PetscObject)bv)->name) {
459: PetscSNPrintf(str,50,"%s_0",((PetscObject)bv)->name);
460: PetscObjectSetName((PetscObject)ctx->A,str);
461: }
463: if (bv->Acreate) {
464: MatCopy(bv->Acreate,ctx->A,SAME_NONZERO_PATTERN);
465: MatDestroy(&bv->Acreate);
466: }
468: VecDuplicateEmpty(bv->t,&bv->cv[0]);
469: VecDuplicateEmpty(bv->t,&bv->cv[1]);
471: bv->ops->mult = BVMult_Mat;
472: bv->ops->multvec = BVMultVec_Mat;
473: bv->ops->multinplace = BVMultInPlace_Mat;
474: bv->ops->multinplacetrans = BVMultInPlaceTranspose_Mat;
475: bv->ops->dot = BVDot_Mat;
476: bv->ops->dotvec = BVDotVec_Mat;
477: bv->ops->dotvec_local = BVDotVec_Local_Mat;
478: bv->ops->scale = BVScale_Mat;
479: bv->ops->norm = BVNorm_Mat;
480: bv->ops->norm_local = BVNorm_Local_Mat;
481: bv->ops->matmult = BVMatMult_Mat;
482: bv->ops->copy = BVCopy_Mat;
483: bv->ops->copycolumn = BVCopyColumn_Mat;
484: bv->ops->resize = BVResize_Mat;
485: bv->ops->getcolumn = BVGetColumn_Mat;
486: bv->ops->restorecolumn = BVRestoreColumn_Mat;
487: bv->ops->getarray = BVGetArray_Mat;
488: bv->ops->restorearray = BVRestoreArray_Mat;
489: bv->ops->getarrayread = BVGetArrayRead_Mat;
490: bv->ops->restorearrayread = BVRestoreArrayRead_Mat;
491: bv->ops->destroy = BVDestroy_Mat;
492: if (!ctx->mpi) bv->ops->view = BVView_Mat;
493: return(0);
494: }