Actual source code: contig.c
slepc-3.18.1 2022-11-02
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 an array of Vecs sharing a contiguous array for elements
12: */
14: #include <slepc/private/bvimpl.h>
16: typedef struct {
17: Vec *V;
18: PetscScalar *array;
19: PetscBool mpi;
20: } BV_CONTIGUOUS;
22: PetscErrorCode BVMult_Contiguous(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
23: {
24: BV_CONTIGUOUS *y = (BV_CONTIGUOUS*)Y->data,*x = (BV_CONTIGUOUS*)X->data;
25: const PetscScalar *q;
26: PetscInt ldq;
28: if (Q) {
29: MatDenseGetLDA(Q,&ldq);
30: MatDenseGetArrayRead(Q,&q);
31: BVMult_BLAS_Private(Y,Y->n,Y->k-Y->l,X->k-X->l,ldq,alpha,x->array+(X->nc+X->l)*X->n,q+Y->l*ldq+X->l,beta,y->array+(Y->nc+Y->l)*Y->n);
32: MatDenseRestoreArrayRead(Q,&q);
33: } else BVAXPY_BLAS_Private(Y,Y->n,Y->k-Y->l,alpha,x->array+(X->nc+X->l)*X->n,beta,y->array+(Y->nc+Y->l)*Y->n);
34: return 0;
35: }
37: PetscErrorCode BVMultVec_Contiguous(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
38: {
39: BV_CONTIGUOUS *x = (BV_CONTIGUOUS*)X->data;
40: PetscScalar *py,*qq=q;
42: VecGetArray(y,&py);
43: if (!q) VecGetArray(X->buffer,&qq);
44: BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,x->array+(X->nc+X->l)*X->n,qq,beta,py);
45: if (!q) VecRestoreArray(X->buffer,&qq);
46: VecRestoreArray(y,&py);
47: return 0;
48: }
50: PetscErrorCode BVMultInPlace_Contiguous(BV V,Mat Q,PetscInt s,PetscInt e)
51: {
52: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)V->data;
53: const PetscScalar *q;
54: PetscInt ldq;
56: MatDenseGetLDA(Q,&ldq);
57: MatDenseGetArrayRead(Q,&q);
58: BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,ldq,s-V->l,e-V->l,ctx->array+(V->nc+V->l)*V->n,q+V->l*ldq+V->l,PETSC_FALSE);
59: MatDenseRestoreArrayRead(Q,&q);
60: return 0;
61: }
63: PetscErrorCode BVMultInPlaceHermitianTranspose_Contiguous(BV V,Mat Q,PetscInt s,PetscInt e)
64: {
65: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)V->data;
66: const PetscScalar *q;
67: PetscInt ldq;
69: MatDenseGetLDA(Q,&ldq);
70: MatDenseGetArrayRead(Q,&q);
71: BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,ldq,s-V->l,e-V->l,ctx->array+(V->nc+V->l)*V->n,q+V->l*ldq+V->l,PETSC_TRUE);
72: MatDenseRestoreArrayRead(Q,&q);
73: return 0;
74: }
76: PetscErrorCode BVDot_Contiguous(BV X,BV Y,Mat M)
77: {
78: BV_CONTIGUOUS *x = (BV_CONTIGUOUS*)X->data,*y = (BV_CONTIGUOUS*)Y->data;
79: PetscScalar *m;
80: PetscInt ldm;
82: MatDenseGetLDA(M,&ldm);
83: MatDenseGetArray(M,&m);
84: BVDot_BLAS_Private(X,Y->k-Y->l,X->k-X->l,X->n,ldm,y->array+(Y->nc+Y->l)*Y->n,x->array+(X->nc+X->l)*X->n,m+X->l*ldm+Y->l,x->mpi);
85: MatDenseRestoreArray(M,&m);
86: return 0;
87: }
89: PetscErrorCode BVDotVec_Contiguous(BV X,Vec y,PetscScalar *q)
90: {
91: BV_CONTIGUOUS *x = (BV_CONTIGUOUS*)X->data;
92: const PetscScalar *py;
93: PetscScalar *qq=q;
94: Vec z = y;
96: if (PetscUnlikely(X->matrix)) {
97: BV_IPMatMult(X,y);
98: z = X->Bx;
99: }
100: VecGetArrayRead(z,&py);
101: if (!q) VecGetArray(X->buffer,&qq);
102: BVDotVec_BLAS_Private(X,X->n,X->k-X->l,x->array+(X->nc+X->l)*X->n,py,qq,x->mpi);
103: if (!q) VecRestoreArray(X->buffer,&qq);
104: VecRestoreArrayRead(z,&py);
105: return 0;
106: }
108: PetscErrorCode BVDotVec_Local_Contiguous(BV X,Vec y,PetscScalar *m)
109: {
110: BV_CONTIGUOUS *x = (BV_CONTIGUOUS*)X->data;
111: PetscScalar *py;
112: Vec z = y;
114: if (PetscUnlikely(X->matrix)) {
115: BV_IPMatMult(X,y);
116: z = X->Bx;
117: }
118: VecGetArray(z,&py);
119: BVDotVec_BLAS_Private(X,X->n,X->k-X->l,x->array+(X->nc+X->l)*X->n,py,m,PETSC_FALSE);
120: VecRestoreArray(z,&py);
121: return 0;
122: }
124: PetscErrorCode BVScale_Contiguous(BV bv,PetscInt j,PetscScalar alpha)
125: {
126: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
128: if (PetscUnlikely(j<0)) BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->n,ctx->array+(bv->nc+bv->l)*bv->n,alpha);
129: else BVScale_BLAS_Private(bv,bv->n,ctx->array+(bv->nc+j)*bv->n,alpha);
130: return 0;
131: }
133: PetscErrorCode BVNorm_Contiguous(BV bv,PetscInt j,NormType type,PetscReal *val)
134: {
135: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
137: if (PetscUnlikely(j<0)) BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,type,val,ctx->mpi);
138: else BVNorm_LAPACK_Private(bv,bv->n,1,ctx->array+(bv->nc+j)*bv->n,type,val,ctx->mpi);
139: return 0;
140: }
142: PetscErrorCode BVNorm_Local_Contiguous(BV bv,PetscInt j,NormType type,PetscReal *val)
143: {
144: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
146: if (PetscUnlikely(j<0)) BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,type,val,PETSC_FALSE);
147: else BVNorm_LAPACK_Private(bv,bv->n,1,ctx->array+(bv->nc+j)*bv->n,type,val,PETSC_FALSE);
148: return 0;
149: }
151: PetscErrorCode BVNormalize_Contiguous(BV bv,PetscScalar *eigi)
152: {
153: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
154: PetscScalar *wi=NULL;
156: if (eigi) wi = eigi+bv->l;
157: BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,wi,ctx->mpi);
158: return 0;
159: }
161: PetscErrorCode BVMatMult_Contiguous(BV V,Mat A,BV W)
162: {
163: BV_CONTIGUOUS *v = (BV_CONTIGUOUS*)V->data,*w = (BV_CONTIGUOUS*)W->data;
164: PetscInt j;
165: Mat Vmat,Wmat;
167: if (V->vmm) {
168: BVGetMat(V,&Vmat);
169: BVGetMat(W,&Wmat);
170: MatProductCreateWithMat(A,Vmat,NULL,Wmat);
171: MatProductSetType(Wmat,MATPRODUCT_AB);
172: MatProductSetFromOptions(Wmat);
173: MatProductSymbolic(Wmat);
174: MatProductNumeric(Wmat);
175: MatProductClear(Wmat);
176: BVRestoreMat(V,&Vmat);
177: BVRestoreMat(W,&Wmat);
178: } else {
179: for (j=0;j<V->k-V->l;j++) MatMult(A,v->V[V->nc+V->l+j],w->V[W->nc+W->l+j]);
180: }
181: return 0;
182: }
184: PetscErrorCode BVCopy_Contiguous(BV V,BV W)
185: {
186: BV_CONTIGUOUS *v = (BV_CONTIGUOUS*)V->data,*w = (BV_CONTIGUOUS*)W->data;
187: PetscScalar *pvc,*pwc;
189: pvc = v->array+(V->nc+V->l)*V->n;
190: pwc = w->array+(W->nc+W->l)*W->n;
191: PetscArraycpy(pwc,pvc,(V->k-V->l)*V->n);
192: return 0;
193: }
195: PetscErrorCode BVCopyColumn_Contiguous(BV V,PetscInt j,PetscInt i)
196: {
197: BV_CONTIGUOUS *v = (BV_CONTIGUOUS*)V->data;
199: PetscArraycpy(v->array+(V->nc+i)*V->n,v->array+(V->nc+j)*V->n,V->n);
200: return 0;
201: }
203: PetscErrorCode BVResize_Contiguous(BV bv,PetscInt m,PetscBool copy)
204: {
205: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
206: PetscInt j,bs;
207: PetscScalar *newarray;
208: Vec *newV;
209: char str[50];
211: VecGetBlockSize(bv->t,&bs);
212: PetscMalloc1(m*bv->n,&newarray);
213: PetscArrayzero(newarray,m*bv->n);
214: PetscMalloc1(m,&newV);
215: for (j=0;j<m;j++) {
216: if (ctx->mpi) VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,PETSC_DECIDE,newarray+j*bv->n,newV+j);
217: else VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,newarray+j*bv->n,newV+j);
218: }
219: if (((PetscObject)bv)->name) {
220: for (j=0;j<m;j++) {
221: PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j);
222: PetscObjectSetName((PetscObject)newV[j],str);
223: }
224: }
225: if (copy) PetscArraycpy(newarray,ctx->array,PetscMin(m,bv->m)*bv->n);
226: VecDestroyVecs(bv->m,&ctx->V);
227: ctx->V = newV;
228: PetscFree(ctx->array);
229: ctx->array = newarray;
230: return 0;
231: }
233: PetscErrorCode BVGetColumn_Contiguous(BV bv,PetscInt j,Vec *v)
234: {
235: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
236: PetscInt l;
238: l = BVAvailableVec;
239: bv->cv[l] = ctx->V[bv->nc+j];
240: return 0;
241: }
243: PetscErrorCode BVRestoreColumn_Contiguous(BV bv,PetscInt j,Vec *v)
244: {
245: PetscInt l;
247: l = (j==bv->ci[0])? 0: 1;
248: bv->cv[l] = NULL;
249: return 0;
250: }
252: PetscErrorCode BVGetArray_Contiguous(BV bv,PetscScalar **a)
253: {
254: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
256: *a = ctx->array;
257: return 0;
258: }
260: PetscErrorCode BVGetArrayRead_Contiguous(BV bv,const PetscScalar **a)
261: {
262: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
264: *a = ctx->array;
265: return 0;
266: }
268: PetscErrorCode BVDestroy_Contiguous(BV bv)
269: {
270: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
272: if (!bv->issplit) {
273: VecDestroyVecs(bv->nc+bv->m,&ctx->V);
274: PetscFree(ctx->array);
275: }
276: PetscFree(bv->data);
277: return 0;
278: }
280: PetscErrorCode BVView_Contiguous(BV bv,PetscViewer viewer)
281: {
282: PetscInt j;
283: Vec v;
284: PetscViewerFormat format;
285: PetscBool isascii,ismatlab=PETSC_FALSE;
286: const char *bvname,*name;
288: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
289: if (isascii) {
290: PetscViewerGetFormat(viewer,&format);
291: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return 0;
292: if (format == PETSC_VIEWER_ASCII_MATLAB) ismatlab = PETSC_TRUE;
293: }
294: if (ismatlab) {
295: PetscObjectGetName((PetscObject)bv,&bvname);
296: PetscViewerASCIIPrintf(viewer,"%s=[];\n",bvname);
297: }
298: for (j=0;j<bv->m;j++) {
299: BVGetColumn(bv,j,&v);
300: VecView(v,viewer);
301: if (ismatlab) {
302: PetscObjectGetName((PetscObject)v,&name);
303: PetscViewerASCIIPrintf(viewer,"%s=[%s,%s];clear %s\n",bvname,bvname,name,name);
304: }
305: BVRestoreColumn(bv,j,&v);
306: }
307: return 0;
308: }
310: SLEPC_EXTERN PetscErrorCode BVCreate_Contiguous(BV bv)
311: {
312: BV_CONTIGUOUS *ctx;
313: PetscInt j,nloc,bs,lsplit;
314: PetscBool seq;
315: PetscScalar *aa;
316: char str[50];
317: PetscScalar *array;
318: BV parent;
319: Vec *Vpar;
321: PetscNew(&ctx);
322: bv->data = (void*)ctx;
324: PetscObjectTypeCompare((PetscObject)bv->t,VECMPI,&ctx->mpi);
325: if (!ctx->mpi) {
326: PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq);
328: }
330: VecGetLocalSize(bv->t,&nloc);
331: VecGetBlockSize(bv->t,&bs);
333: if (PetscUnlikely(bv->issplit)) {
334: /* split BV: share memory and Vecs of the parent BV */
335: parent = bv->splitparent;
336: lsplit = parent->lsplit;
337: Vpar = ((BV_CONTIGUOUS*)parent->data)->V;
338: ctx->V = (bv->issplit==1)? Vpar: Vpar+lsplit;
339: array = ((BV_CONTIGUOUS*)parent->data)->array;
340: ctx->array = (bv->issplit==1)? array: array+lsplit*nloc;
341: } else {
342: /* regular BV: allocate memory and Vecs for the BV entries */
343: PetscCalloc1(bv->m*nloc,&ctx->array);
344: PetscMalloc1(bv->m,&ctx->V);
345: for (j=0;j<bv->m;j++) {
346: if (ctx->mpi) VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,nloc,PETSC_DECIDE,ctx->array+j*nloc,ctx->V+j);
347: else VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,nloc,ctx->array+j*nloc,ctx->V+j);
348: }
349: }
350: if (((PetscObject)bv)->name) {
351: for (j=0;j<bv->m;j++) {
352: PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j);
353: PetscObjectSetName((PetscObject)ctx->V[j],str);
354: }
355: }
357: if (PetscUnlikely(bv->Acreate)) {
358: MatDenseGetArray(bv->Acreate,&aa);
359: PetscArraycpy(ctx->array,aa,bv->m*nloc);
360: MatDenseRestoreArray(bv->Acreate,&aa);
361: MatDestroy(&bv->Acreate);
362: }
364: bv->ops->mult = BVMult_Contiguous;
365: bv->ops->multvec = BVMultVec_Contiguous;
366: bv->ops->multinplace = BVMultInPlace_Contiguous;
367: bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Contiguous;
368: bv->ops->dot = BVDot_Contiguous;
369: bv->ops->dotvec = BVDotVec_Contiguous;
370: bv->ops->dotvec_local = BVDotVec_Local_Contiguous;
371: bv->ops->scale = BVScale_Contiguous;
372: bv->ops->norm = BVNorm_Contiguous;
373: bv->ops->norm_local = BVNorm_Local_Contiguous;
374: bv->ops->normalize = BVNormalize_Contiguous;
375: bv->ops->matmult = BVMatMult_Contiguous;
376: bv->ops->copy = BVCopy_Contiguous;
377: bv->ops->copycolumn = BVCopyColumn_Contiguous;
378: bv->ops->resize = BVResize_Contiguous;
379: bv->ops->getcolumn = BVGetColumn_Contiguous;
380: bv->ops->restorecolumn = BVRestoreColumn_Contiguous;
381: bv->ops->getarray = BVGetArray_Contiguous;
382: bv->ops->getarrayread = BVGetArrayRead_Contiguous;
383: bv->ops->getmat = BVGetMat_Default;
384: bv->ops->restoremat = BVRestoreMat_Default;
385: bv->ops->destroy = BVDestroy_Contiguous;
386: bv->ops->view = BVView_Contiguous;
387: return 0;
388: }