Actual source code: bvmat.c
slepc-3.19.2 2023-09-05
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 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: {
23: BV_MAT *y = (BV_MAT*)Y->data,*x = (BV_MAT*)X->data;
24: PetscScalar *py;
25: const PetscScalar *px,*q;
26: PetscInt ldq;
28: PetscFunctionBegin;
29: PetscCall(MatDenseGetArrayRead(x->A,&px));
30: PetscCall(MatDenseGetArray(y->A,&py));
31: if (Q) {
32: PetscCall(MatDenseGetLDA(Q,&ldq));
33: PetscCall(MatDenseGetArrayRead(Q,&q));
34: PetscCall(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: PetscCall(MatDenseRestoreArrayRead(Q,&q));
36: } else PetscCall(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));
37: PetscCall(MatDenseRestoreArrayRead(x->A,&px));
38: PetscCall(MatDenseRestoreArray(y->A,&py));
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: PetscErrorCode BVMultVec_Mat(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
43: {
44: BV_MAT *x = (BV_MAT*)X->data;
45: PetscScalar *py,*qq=q;
46: const PetscScalar *px;
48: PetscFunctionBegin;
49: PetscCall(MatDenseGetArrayRead(x->A,&px));
50: PetscCall(VecGetArray(y,&py));
51: if (!q) PetscCall(VecGetArray(X->buffer,&qq));
52: PetscCall(BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,px+(X->nc+X->l)*X->n,qq,beta,py));
53: if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
54: PetscCall(MatDenseRestoreArrayRead(x->A,&px));
55: PetscCall(VecRestoreArray(y,&py));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: PetscErrorCode BVMultInPlace_Mat(BV V,Mat Q,PetscInt s,PetscInt e)
60: {
61: BV_MAT *ctx = (BV_MAT*)V->data;
62: PetscScalar *pv;
63: const PetscScalar *q;
64: PetscInt ldq;
66: PetscFunctionBegin;
67: PetscCall(MatDenseGetLDA(Q,&ldq));
68: PetscCall(MatDenseGetArray(ctx->A,&pv));
69: PetscCall(MatDenseGetArrayRead(Q,&q));
70: PetscCall(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));
71: PetscCall(MatDenseRestoreArrayRead(Q,&q));
72: PetscCall(MatDenseRestoreArray(ctx->A,&pv));
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: PetscErrorCode BVMultInPlaceHermitianTranspose_Mat(BV V,Mat Q,PetscInt s,PetscInt e)
77: {
78: BV_MAT *ctx = (BV_MAT*)V->data;
79: PetscScalar *pv;
80: const PetscScalar *q;
81: PetscInt ldq;
83: PetscFunctionBegin;
84: PetscCall(MatDenseGetLDA(Q,&ldq));
85: PetscCall(MatDenseGetArray(ctx->A,&pv));
86: PetscCall(MatDenseGetArrayRead(Q,&q));
87: PetscCall(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));
88: PetscCall(MatDenseRestoreArrayRead(Q,&q));
89: PetscCall(MatDenseRestoreArray(ctx->A,&pv));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: PetscErrorCode BVDot_Mat(BV X,BV Y,Mat M)
94: {
95: BV_MAT *x = (BV_MAT*)X->data,*y = (BV_MAT*)Y->data;
96: PetscScalar *m;
97: const PetscScalar *px,*py;
98: PetscInt ldm;
100: PetscFunctionBegin;
101: PetscCall(MatDenseGetLDA(M,&ldm));
102: PetscCall(MatDenseGetArrayRead(x->A,&px));
103: PetscCall(MatDenseGetArrayRead(y->A,&py));
104: PetscCall(MatDenseGetArray(M,&m));
105: PetscCall(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: PetscCall(MatDenseRestoreArray(M,&m));
107: PetscCall(MatDenseRestoreArrayRead(x->A,&px));
108: PetscCall(MatDenseRestoreArrayRead(y->A,&py));
109: PetscFunctionReturn(PETSC_SUCCESS);
110: }
112: PetscErrorCode BVDotVec_Mat(BV X,Vec y,PetscScalar *q)
113: {
114: BV_MAT *x = (BV_MAT*)X->data;
115: PetscScalar *qq=q;
116: const PetscScalar *px,*py;
117: Vec z = y;
119: PetscFunctionBegin;
120: if (PetscUnlikely(X->matrix)) {
121: PetscCall(BV_IPMatMult(X,y));
122: z = X->Bx;
123: }
124: PetscCall(MatDenseGetArrayRead(x->A,&px));
125: PetscCall(VecGetArrayRead(z,&py));
126: if (!q) PetscCall(VecGetArray(X->buffer,&qq));
127: PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,qq,x->mpi));
128: if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
129: PetscCall(VecRestoreArrayRead(z,&py));
130: PetscCall(MatDenseRestoreArrayRead(x->A,&px));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: PetscErrorCode BVDotVec_Local_Mat(BV X,Vec y,PetscScalar *m)
135: {
136: BV_MAT *x = (BV_MAT*)X->data;
137: const PetscScalar *px,*py;
138: Vec z = y;
140: PetscFunctionBegin;
141: if (PetscUnlikely(X->matrix)) {
142: PetscCall(BV_IPMatMult(X,y));
143: z = X->Bx;
144: }
145: PetscCall(MatDenseGetArrayRead(x->A,&px));
146: PetscCall(VecGetArrayRead(z,&py));
147: PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,m,PETSC_FALSE));
148: PetscCall(VecRestoreArrayRead(z,&py));
149: PetscCall(MatDenseRestoreArrayRead(x->A,&px));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: PetscErrorCode BVScale_Mat(BV bv,PetscInt j,PetscScalar alpha)
154: {
155: BV_MAT *ctx = (BV_MAT*)bv->data;
156: PetscScalar *array;
158: PetscFunctionBegin;
159: PetscCall(MatDenseGetArray(ctx->A,&array));
160: if (PetscUnlikely(j<0)) PetscCall(BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->n,array+(bv->nc+bv->l)*bv->n,alpha));
161: else PetscCall(BVScale_BLAS_Private(bv,bv->n,array+(bv->nc+j)*bv->n,alpha));
162: PetscCall(MatDenseRestoreArray(ctx->A,&array));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: PetscErrorCode BVNorm_Mat(BV bv,PetscInt j,NormType type,PetscReal *val)
167: {
168: BV_MAT *ctx = (BV_MAT*)bv->data;
169: const PetscScalar *array;
171: PetscFunctionBegin;
172: PetscCall(MatDenseGetArrayRead(ctx->A,&array));
173: if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,ctx->mpi));
174: else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,ctx->mpi));
175: PetscCall(MatDenseRestoreArrayRead(ctx->A,&array));
176: PetscFunctionReturn(PETSC_SUCCESS);
177: }
179: PetscErrorCode BVNorm_Local_Mat(BV bv,PetscInt j,NormType type,PetscReal *val)
180: {
181: BV_MAT *ctx = (BV_MAT*)bv->data;
182: const PetscScalar *array;
184: PetscFunctionBegin;
185: PetscCall(MatDenseGetArrayRead(ctx->A,&array));
186: if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,PETSC_FALSE));
187: else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,PETSC_FALSE));
188: PetscCall(MatDenseRestoreArrayRead(ctx->A,&array));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: PetscErrorCode BVNormalize_Mat(BV bv,PetscScalar *eigi)
193: {
194: BV_MAT *ctx = (BV_MAT*)bv->data;
195: PetscScalar *array,*wi=NULL;
197: PetscFunctionBegin;
198: PetscCall(MatDenseGetArray(ctx->A,&array));
199: if (eigi) wi = eigi+bv->l;
200: PetscCall(BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,wi,ctx->mpi));
201: PetscCall(MatDenseRestoreArray(ctx->A,&array));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: PetscErrorCode BVMatMult_Mat(BV V,Mat A,BV W)
206: {
207: PetscInt j;
208: Mat Vmat,Wmat;
209: Vec vv,ww;
211: PetscFunctionBegin;
212: if (V->vmm) {
213: PetscCall(BVGetMat(V,&Vmat));
214: PetscCall(BVGetMat(W,&Wmat));
215: PetscCall(MatProductCreateWithMat(A,Vmat,NULL,Wmat));
216: PetscCall(MatProductSetType(Wmat,MATPRODUCT_AB));
217: PetscCall(MatProductSetFromOptions(Wmat));
218: PetscCall(MatProductSymbolic(Wmat));
219: PetscCall(MatProductNumeric(Wmat));
220: PetscCall(MatProductClear(Wmat));
221: PetscCall(BVRestoreMat(V,&Vmat));
222: PetscCall(BVRestoreMat(W,&Wmat));
223: } else {
224: for (j=0;j<V->k-V->l;j++) {
225: PetscCall(BVGetColumn(V,V->l+j,&vv));
226: PetscCall(BVGetColumn(W,W->l+j,&ww));
227: PetscCall(MatMult(A,vv,ww));
228: PetscCall(BVRestoreColumn(V,V->l+j,&vv));
229: PetscCall(BVRestoreColumn(W,W->l+j,&ww));
230: }
231: }
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: PetscErrorCode BVCopy_Mat(BV V,BV W)
236: {
237: BV_MAT *v = (BV_MAT*)V->data,*w = (BV_MAT*)W->data;
238: PetscScalar *pw,*pwc;
239: const PetscScalar *pv,*pvc;
241: PetscFunctionBegin;
242: PetscCall(MatDenseGetArrayRead(v->A,&pv));
243: PetscCall(MatDenseGetArray(w->A,&pw));
244: pvc = pv+(V->nc+V->l)*V->n;
245: pwc = pw+(W->nc+W->l)*W->n;
246: PetscCall(PetscArraycpy(pwc,pvc,(V->k-V->l)*V->n));
247: PetscCall(MatDenseRestoreArrayRead(v->A,&pv));
248: PetscCall(MatDenseRestoreArray(w->A,&pw));
249: PetscFunctionReturn(PETSC_SUCCESS);
250: }
252: PetscErrorCode BVCopyColumn_Mat(BV V,PetscInt j,PetscInt i)
253: {
254: BV_MAT *v = (BV_MAT*)V->data;
255: PetscScalar *pv;
257: PetscFunctionBegin;
258: PetscCall(MatDenseGetArray(v->A,&pv));
259: PetscCall(PetscArraycpy(pv+(V->nc+i)*V->n,pv+(V->nc+j)*V->n,V->n));
260: PetscCall(MatDenseRestoreArray(v->A,&pv));
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: PetscErrorCode BVResize_Mat(BV bv,PetscInt m,PetscBool copy)
265: {
266: BV_MAT *ctx = (BV_MAT*)bv->data;
267: PetscScalar *pnew;
268: const PetscScalar *pA;
269: Mat A;
270: char str[50];
272: PetscFunctionBegin;
273: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,PETSC_DECIDE,m,NULL,&A));
274: if (((PetscObject)bv)->name) {
275: PetscCall(PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name));
276: PetscCall(PetscObjectSetName((PetscObject)A,str));
277: }
278: if (copy) {
279: PetscCall(MatDenseGetArrayRead(ctx->A,&pA));
280: PetscCall(MatDenseGetArrayWrite(A,&pnew));
281: PetscCall(PetscArraycpy(pnew,pA,PetscMin(m,bv->m)*bv->n));
282: PetscCall(MatDenseRestoreArrayRead(ctx->A,&pA));
283: PetscCall(MatDenseRestoreArrayWrite(A,&pnew));
284: }
285: PetscCall(MatDestroy(&ctx->A));
286: ctx->A = A;
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: PetscErrorCode BVGetColumn_Mat(BV bv,PetscInt j,Vec *v)
291: {
292: BV_MAT *ctx = (BV_MAT*)bv->data;
293: PetscScalar *pA;
294: PetscInt l;
296: PetscFunctionBegin;
297: l = BVAvailableVec;
298: PetscCall(MatDenseGetArray(ctx->A,&pA));
299: PetscCall(VecPlaceArray(bv->cv[l],pA+(bv->nc+j)*bv->n));
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: PetscErrorCode BVRestoreColumn_Mat(BV bv,PetscInt j,Vec *v)
304: {
305: BV_MAT *ctx = (BV_MAT*)bv->data;
306: PetscScalar *pA;
307: PetscInt l;
309: PetscFunctionBegin;
310: l = (j==bv->ci[0])? 0: 1;
311: PetscCall(VecResetArray(bv->cv[l]));
312: PetscCall(MatDenseRestoreArray(ctx->A,&pA));
313: PetscFunctionReturn(PETSC_SUCCESS);
314: }
316: PetscErrorCode BVGetArray_Mat(BV bv,PetscScalar **a)
317: {
318: BV_MAT *ctx = (BV_MAT*)bv->data;
320: PetscFunctionBegin;
321: PetscCall(MatDenseGetArray(ctx->A,a));
322: PetscFunctionReturn(PETSC_SUCCESS);
323: }
325: PetscErrorCode BVRestoreArray_Mat(BV bv,PetscScalar **a)
326: {
327: BV_MAT *ctx = (BV_MAT*)bv->data;
329: PetscFunctionBegin;
330: if (a) PetscCall(MatDenseRestoreArray(ctx->A,a));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: PetscErrorCode BVGetArrayRead_Mat(BV bv,const PetscScalar **a)
335: {
336: BV_MAT *ctx = (BV_MAT*)bv->data;
338: PetscFunctionBegin;
339: PetscCall(MatDenseGetArrayRead(ctx->A,a));
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: PetscErrorCode BVRestoreArrayRead_Mat(BV bv,const PetscScalar **a)
344: {
345: BV_MAT *ctx = (BV_MAT*)bv->data;
347: PetscFunctionBegin;
348: if (a) PetscCall(MatDenseRestoreArrayRead(ctx->A,a));
349: PetscFunctionReturn(PETSC_SUCCESS);
350: }
352: PetscErrorCode BVView_Mat(BV bv,PetscViewer viewer)
353: {
354: BV_MAT *ctx = (BV_MAT*)bv->data;
355: PetscViewerFormat format;
356: PetscBool isascii;
357: const char *bvname,*name;
359: PetscFunctionBegin;
360: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
361: if (isascii) {
362: PetscCall(PetscViewerGetFormat(viewer,&format));
363: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
364: PetscCall(MatView(ctx->A,viewer));
365: if (format == PETSC_VIEWER_ASCII_MATLAB) {
366: PetscCall(PetscObjectGetName((PetscObject)bv,&bvname));
367: PetscCall(PetscObjectGetName((PetscObject)ctx->A,&name));
368: PetscCall(PetscViewerASCIIPrintf(viewer,"%s=%s;clear %s\n",bvname,name,name));
369: if (bv->nc) PetscCall(PetscViewerASCIIPrintf(viewer,"%s=%s(:,%" PetscInt_FMT ":end);\n",bvname,bvname,bv->nc+1));
370: }
371: } else PetscCall(MatView(ctx->A,viewer));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: PetscErrorCode BVDestroy_Mat(BV bv)
376: {
377: BV_MAT *ctx = (BV_MAT*)bv->data;
379: PetscFunctionBegin;
380: PetscCall(MatDestroy(&ctx->A));
381: PetscCall(VecDestroy(&bv->cv[0]));
382: PetscCall(VecDestroy(&bv->cv[1]));
383: PetscCall(PetscFree(bv->data));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: SLEPC_EXTERN PetscErrorCode BVCreate_Mat(BV bv)
388: {
389: BV_MAT *ctx;
390: PetscInt nloc,bs,lsplit;
391: PetscBool seq;
392: char str[50];
393: PetscScalar *array,*ptr;
394: BV parent;
396: PetscFunctionBegin;
397: PetscCall(PetscNew(&ctx));
398: bv->data = (void*)ctx;
400: PetscCall(PetscObjectTypeCompare((PetscObject)bv->t,VECMPI,&ctx->mpi));
401: if (!ctx->mpi) {
402: PetscCall(PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq));
403: PetscCheck(seq,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot create a BVMAT from a non-standard template vector");
404: }
406: PetscCall(VecGetLocalSize(bv->t,&nloc));
407: PetscCall(VecGetBlockSize(bv->t,&bs));
409: if (PetscUnlikely(bv->issplit)) {
410: /* split BV: share the memory of the parent BV */
411: parent = bv->splitparent;
412: lsplit = parent->lsplit;
413: PetscCall(MatDenseGetArray(((BV_MAT*)parent->data)->A,&array));
414: ptr = (bv->issplit==1)? array: array+lsplit*nloc;
415: PetscCall(MatDenseRestoreArray(((BV_MAT*)parent->data)->A,&array));
416: } else {
417: /* regular BV: allocate memory for the BV entries */
418: ptr = NULL;
419: }
420: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)bv->t),nloc,PETSC_DECIDE,PETSC_DECIDE,bv->m,ptr,&ctx->A));
421: if (((PetscObject)bv)->name) {
422: PetscCall(PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name));
423: PetscCall(PetscObjectSetName((PetscObject)ctx->A,str));
424: }
426: if (PetscUnlikely(bv->Acreate)) {
427: PetscCall(MatCopy(bv->Acreate,ctx->A,SAME_NONZERO_PATTERN));
428: PetscCall(MatDestroy(&bv->Acreate));
429: }
431: PetscCall(VecDuplicateEmpty(bv->t,&bv->cv[0]));
432: PetscCall(VecDuplicateEmpty(bv->t,&bv->cv[1]));
434: bv->ops->mult = BVMult_Mat;
435: bv->ops->multvec = BVMultVec_Mat;
436: bv->ops->multinplace = BVMultInPlace_Mat;
437: bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Mat;
438: bv->ops->dot = BVDot_Mat;
439: bv->ops->dotvec = BVDotVec_Mat;
440: bv->ops->dotvec_local = BVDotVec_Local_Mat;
441: bv->ops->scale = BVScale_Mat;
442: bv->ops->norm = BVNorm_Mat;
443: bv->ops->norm_local = BVNorm_Local_Mat;
444: bv->ops->normalize = BVNormalize_Mat;
445: bv->ops->matmult = BVMatMult_Mat;
446: bv->ops->copy = BVCopy_Mat;
447: bv->ops->copycolumn = BVCopyColumn_Mat;
448: bv->ops->resize = BVResize_Mat;
449: bv->ops->getcolumn = BVGetColumn_Mat;
450: bv->ops->restorecolumn = BVRestoreColumn_Mat;
451: bv->ops->getarray = BVGetArray_Mat;
452: bv->ops->restorearray = BVRestoreArray_Mat;
453: bv->ops->getarrayread = BVGetArrayRead_Mat;
454: bv->ops->restorearrayread = BVRestoreArrayRead_Mat;
455: bv->ops->getmat = BVGetMat_Default;
456: bv->ops->restoremat = BVRestoreMat_Default;
457: bv->ops->destroy = BVDestroy_Mat;
458: if (!ctx->mpi) bv->ops->view = BVView_Mat;
459: PetscFunctionReturn(PETSC_SUCCESS);
460: }