Actual source code: contig.c
slepc-3.21.1 2024-04-26
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: static 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: PetscFunctionBegin;
29: if (Q) {
30: PetscCall(MatDenseGetLDA(Q,&ldq));
31: PetscCall(MatDenseGetArrayRead(Q,&q));
32: PetscCall(BVMult_BLAS_Private(Y,Y->n,Y->k-Y->l,X->k-X->l,alpha,x->array+(X->nc+X->l)*X->ld,X->ld,q+Y->l*ldq+X->l,ldq,beta,y->array+(Y->nc+Y->l)*Y->ld,Y->ld));
33: PetscCall(MatDenseRestoreArrayRead(Q,&q));
34: } else PetscCall(BVAXPY_BLAS_Private(Y,Y->n,Y->k-Y->l,alpha,x->array+(X->nc+X->l)*X->ld,X->ld,beta,y->array+(Y->nc+Y->l)*Y->ld,Y->ld));
35: PetscFunctionReturn(PETSC_SUCCESS);
36: }
38: static PetscErrorCode BVMultVec_Contiguous(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
39: {
40: BV_CONTIGUOUS *x = (BV_CONTIGUOUS*)X->data;
41: PetscScalar *py,*qq=q;
43: PetscFunctionBegin;
44: PetscCall(VecGetArray(y,&py));
45: if (!q) PetscCall(VecGetArray(X->buffer,&qq));
46: PetscCall(BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,x->array+(X->nc+X->l)*X->ld,X->ld,qq,beta,py));
47: if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
48: PetscCall(VecRestoreArray(y,&py));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: static PetscErrorCode BVMultInPlace_Contiguous(BV V,Mat Q,PetscInt s,PetscInt e)
53: {
54: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)V->data;
55: const PetscScalar *q;
56: PetscInt ldq;
58: PetscFunctionBegin;
59: if (s>=e || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
60: PetscCall(MatDenseGetLDA(Q,&ldq));
61: PetscCall(MatDenseGetArrayRead(Q,&q));
62: PetscCall(BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,s-V->l,e-V->l,ctx->array+(V->nc+V->l)*V->ld,V->ld,q+V->l*ldq+V->l,ldq,PETSC_FALSE));
63: PetscCall(MatDenseRestoreArrayRead(Q,&q));
64: PetscFunctionReturn(PETSC_SUCCESS);
65: }
67: static PetscErrorCode BVMultInPlaceHermitianTranspose_Contiguous(BV V,Mat Q,PetscInt s,PetscInt e)
68: {
69: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)V->data;
70: const PetscScalar *q;
71: PetscInt ldq;
73: PetscFunctionBegin;
74: if (s>=e || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
75: PetscCall(MatDenseGetLDA(Q,&ldq));
76: PetscCall(MatDenseGetArrayRead(Q,&q));
77: PetscCall(BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,s-V->l,e-V->l,ctx->array+(V->nc+V->l)*V->ld,V->ld,q+V->l*ldq+V->l,ldq,PETSC_TRUE));
78: PetscCall(MatDenseRestoreArrayRead(Q,&q));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: static PetscErrorCode BVDot_Contiguous(BV X,BV Y,Mat M)
83: {
84: BV_CONTIGUOUS *x = (BV_CONTIGUOUS*)X->data,*y = (BV_CONTIGUOUS*)Y->data;
85: PetscScalar *m;
86: PetscInt ldm;
88: PetscFunctionBegin;
89: PetscCall(MatDenseGetLDA(M,&ldm));
90: PetscCall(MatDenseGetArray(M,&m));
91: PetscCall(BVDot_BLAS_Private(X,Y->k-Y->l,X->k-X->l,X->n,y->array+(Y->nc+Y->l)*Y->ld,Y->ld,x->array+(X->nc+X->l)*X->ld,X->ld,m+X->l*ldm+Y->l,ldm,x->mpi));
92: PetscCall(MatDenseRestoreArray(M,&m));
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode BVDotVec_Contiguous(BV X,Vec y,PetscScalar *q)
97: {
98: BV_CONTIGUOUS *x = (BV_CONTIGUOUS*)X->data;
99: const PetscScalar *py;
100: PetscScalar *qq=q;
101: Vec z = y;
103: PetscFunctionBegin;
104: if (PetscUnlikely(X->matrix)) {
105: PetscCall(BV_IPMatMult(X,y));
106: z = X->Bx;
107: }
108: PetscCall(VecGetArrayRead(z,&py));
109: if (!q) PetscCall(VecGetArray(X->buffer,&qq));
110: PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,x->array+(X->nc+X->l)*X->ld,X->ld,py,qq,x->mpi));
111: if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
112: PetscCall(VecRestoreArrayRead(z,&py));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode BVDotVec_Local_Contiguous(BV X,Vec y,PetscScalar *m)
117: {
118: BV_CONTIGUOUS *x = (BV_CONTIGUOUS*)X->data;
119: PetscScalar *py;
120: Vec z = y;
122: PetscFunctionBegin;
123: if (PetscUnlikely(X->matrix)) {
124: PetscCall(BV_IPMatMult(X,y));
125: z = X->Bx;
126: }
127: PetscCall(VecGetArray(z,&py));
128: PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,x->array+(X->nc+X->l)*X->ld,X->ld,py,m,PETSC_FALSE));
129: PetscCall(VecRestoreArray(z,&py));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: static PetscErrorCode BVScale_Contiguous(BV bv,PetscInt j,PetscScalar alpha)
134: {
135: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
137: PetscFunctionBegin;
138: if (!bv->n) PetscFunctionReturn(PETSC_SUCCESS);
139: if (PetscUnlikely(j<0)) PetscCall(BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->ld,ctx->array+(bv->nc+bv->l)*bv->ld,alpha));
140: else PetscCall(BVScale_BLAS_Private(bv,bv->n,ctx->array+(bv->nc+j)*bv->ld,alpha));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode BVNorm_Contiguous(BV bv,PetscInt j,NormType type,PetscReal *val)
145: {
146: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
148: PetscFunctionBegin;
149: if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->ld,bv->ld,type,val,ctx->mpi));
150: else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,ctx->array+(bv->nc+j)*bv->ld,bv->ld,type,val,ctx->mpi));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode BVNorm_Local_Contiguous(BV bv,PetscInt j,NormType type,PetscReal *val)
155: {
156: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
158: PetscFunctionBegin;
159: if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->ld,bv->ld,type,val,PETSC_FALSE));
160: else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,ctx->array+(bv->nc+j)*bv->ld,bv->ld,type,val,PETSC_FALSE));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: static PetscErrorCode BVNormalize_Contiguous(BV bv,PetscScalar *eigi)
165: {
166: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
167: PetscScalar *wi=NULL;
169: PetscFunctionBegin;
170: if (eigi) wi = eigi+bv->l;
171: PetscCall(BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->ld,bv->ld,wi,ctx->mpi));
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: static PetscErrorCode BVMatMult_Contiguous(BV V,Mat A,BV W)
176: {
177: BV_CONTIGUOUS *v = (BV_CONTIGUOUS*)V->data,*w = (BV_CONTIGUOUS*)W->data;
178: PetscInt j;
179: Mat Vmat,Wmat;
181: PetscFunctionBegin;
182: if (V->vmm) {
183: PetscCall(BVGetMat(V,&Vmat));
184: PetscCall(BVGetMat(W,&Wmat));
185: PetscCall(MatProductCreateWithMat(A,Vmat,NULL,Wmat));
186: PetscCall(MatProductSetType(Wmat,MATPRODUCT_AB));
187: PetscCall(MatProductSetFromOptions(Wmat));
188: PetscCall(MatProductSymbolic(Wmat));
189: PetscCall(MatProductNumeric(Wmat));
190: PetscCall(MatProductClear(Wmat));
191: PetscCall(BVRestoreMat(V,&Vmat));
192: PetscCall(BVRestoreMat(W,&Wmat));
193: } else {
194: for (j=0;j<V->k-V->l;j++) PetscCall(MatMult(A,v->V[V->nc+V->l+j],w->V[W->nc+W->l+j]));
195: }
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: static PetscErrorCode BVCopy_Contiguous(BV V,BV W)
200: {
201: BV_CONTIGUOUS *v = (BV_CONTIGUOUS*)V->data,*w = (BV_CONTIGUOUS*)W->data;
202: PetscInt j;
204: PetscFunctionBegin;
205: for (j=0;j<V->k-V->l;j++) PetscCall(PetscArraycpy(w->array+(W->nc+W->l+j)*W->ld,v->array+(V->nc+V->l+j)*V->ld,V->n));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode BVCopyColumn_Contiguous(BV V,PetscInt j,PetscInt i)
210: {
211: BV_CONTIGUOUS *v = (BV_CONTIGUOUS*)V->data;
213: PetscFunctionBegin;
214: PetscCall(PetscArraycpy(v->array+(V->nc+i)*V->ld,v->array+(V->nc+j)*V->ld,V->n));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: static PetscErrorCode BVResize_Contiguous(BV bv,PetscInt m,PetscBool copy)
219: {
220: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
221: PetscInt j,bs;
222: PetscScalar *newarray;
223: Vec *newV;
224: char str[50];
226: PetscFunctionBegin;
227: PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
228: PetscCall(PetscCalloc1(m*bv->ld,&newarray));
229: PetscCall(PetscMalloc1(m,&newV));
230: for (j=0;j<m;j++) {
231: if (ctx->mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,bv->n,PETSC_DECIDE,newarray+j*bv->ld,newV+j));
232: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,bv->n,newarray+j*bv->ld,newV+j));
233: }
234: if (((PetscObject)bv)->name) {
235: for (j=0;j<m;j++) {
236: PetscCall(PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j));
237: PetscCall(PetscObjectSetName((PetscObject)newV[j],str));
238: }
239: }
240: if (copy) PetscCall(PetscArraycpy(newarray,ctx->array,PetscMin(m,bv->m)*bv->ld));
241: PetscCall(VecDestroyVecs(bv->m,&ctx->V));
242: ctx->V = newV;
243: PetscCall(PetscFree(ctx->array));
244: ctx->array = newarray;
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode BVGetColumn_Contiguous(BV bv,PetscInt j,Vec *v)
249: {
250: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
251: PetscInt l;
253: PetscFunctionBegin;
254: l = BVAvailableVec;
255: bv->cv[l] = ctx->V[bv->nc+j];
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: static PetscErrorCode BVRestoreColumn_Contiguous(BV bv,PetscInt j,Vec *v)
260: {
261: PetscInt l;
263: PetscFunctionBegin;
264: l = (j==bv->ci[0])? 0: 1;
265: bv->cv[l] = NULL;
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: static PetscErrorCode BVGetArray_Contiguous(BV bv,PetscScalar **a)
270: {
271: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
273: PetscFunctionBegin;
274: *a = ctx->array;
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
278: static PetscErrorCode BVGetArrayRead_Contiguous(BV bv,const PetscScalar **a)
279: {
280: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
282: PetscFunctionBegin;
283: *a = ctx->array;
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode BVDestroy_Contiguous(BV bv)
288: {
289: BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
291: PetscFunctionBegin;
292: if (!bv->issplit) {
293: PetscCall(VecDestroyVecs(bv->nc+bv->m,&ctx->V));
294: PetscCall(PetscFree(ctx->array));
295: }
296: PetscCall(PetscFree(bv->data));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode BVView_Contiguous(BV bv,PetscViewer viewer)
301: {
302: PetscInt j;
303: Vec v;
304: PetscViewerFormat format;
305: PetscBool isascii,ismatlab=PETSC_FALSE;
306: const char *bvname,*name;
308: PetscFunctionBegin;
309: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
310: if (isascii) {
311: PetscCall(PetscViewerGetFormat(viewer,&format));
312: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
313: if (format == PETSC_VIEWER_ASCII_MATLAB) ismatlab = PETSC_TRUE;
314: }
315: if (ismatlab) {
316: PetscCall(PetscObjectGetName((PetscObject)bv,&bvname));
317: PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[];\n",bvname));
318: }
319: for (j=0;j<bv->m;j++) {
320: PetscCall(BVGetColumn(bv,j,&v));
321: PetscCall(VecView(v,viewer));
322: if (ismatlab) {
323: PetscCall(PetscObjectGetName((PetscObject)v,&name));
324: PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[%s,%s];clear %s\n",bvname,bvname,name,name));
325: }
326: PetscCall(BVRestoreColumn(bv,j,&v));
327: }
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: SLEPC_EXTERN PetscErrorCode BVCreate_Contiguous(BV bv)
332: {
333: BV_CONTIGUOUS *ctx;
334: PetscInt j,nloc,bs,lsplit,lda;
335: PetscBool seq,isdense;
336: PetscScalar *aa;
337: char str[50];
338: PetscScalar *array;
339: BV parent;
340: Vec *Vpar;
341: MatType mtype;
343: PetscFunctionBegin;
344: PetscCall(PetscNew(&ctx));
345: bv->data = (void*)ctx;
347: PetscCall(PetscStrcmp(bv->vtype,VECMPI,&ctx->mpi));
348: if (!ctx->mpi) {
349: PetscCall(PetscStrcmp(bv->vtype,VECSEQ,&seq));
350: PetscCheck(seq,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot create a contiguous BV from a non-standard vector type: %s",bv->vtype);
351: }
353: PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
354: PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
355: PetscCall(BV_SetDefaultLD(bv,nloc));
357: if (PetscUnlikely(bv->issplit)) {
358: /* split BV: share memory and Vecs of the parent BV */
359: parent = bv->splitparent;
360: lsplit = parent->lsplit;
361: Vpar = ((BV_CONTIGUOUS*)parent->data)->V;
362: ctx->V = (bv->issplit==1)? Vpar: Vpar+lsplit;
363: array = ((BV_CONTIGUOUS*)parent->data)->array;
364: ctx->array = (bv->issplit==1)? array: array+lsplit*bv->ld;
365: } else {
366: /* regular BV: allocate memory and Vecs for the BV entries */
367: PetscCall(PetscCalloc1(bv->m*bv->ld,&ctx->array));
368: PetscCall(PetscMalloc1(bv->m,&ctx->V));
369: for (j=0;j<bv->m;j++) {
370: if (ctx->mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,PETSC_DECIDE,ctx->array+j*bv->ld,ctx->V+j));
371: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,nloc,ctx->array+j*bv->ld,ctx->V+j));
372: }
373: }
374: if (((PetscObject)bv)->name) {
375: for (j=0;j<bv->m;j++) {
376: PetscCall(PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j));
377: PetscCall(PetscObjectSetName((PetscObject)ctx->V[j],str));
378: }
379: }
381: if (PetscUnlikely(bv->Acreate)) {
382: PetscCall(MatGetType(bv->Acreate,&mtype));
383: PetscCall(PetscStrcmpAny(mtype,&isdense,MATSEQDENSE,MATMPIDENSE,""));
384: PetscCheck(isdense,PetscObjectComm((PetscObject)bv->Acreate),PETSC_ERR_SUP,"BVCONTIGUOUS requires a dense matrix in BVCreateFromMat()");
385: PetscCall(MatDenseGetArray(bv->Acreate,&aa));
386: PetscCall(MatDenseGetLDA(bv->Acreate,&lda));
387: for (j=0;j<bv->m;j++) PetscCall(PetscArraycpy(ctx->array+j*bv->ld,aa+j*lda,bv->n));
388: PetscCall(MatDenseRestoreArray(bv->Acreate,&aa));
389: PetscCall(MatDestroy(&bv->Acreate));
390: }
392: bv->ops->mult = BVMult_Contiguous;
393: bv->ops->multvec = BVMultVec_Contiguous;
394: bv->ops->multinplace = BVMultInPlace_Contiguous;
395: bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Contiguous;
396: bv->ops->dot = BVDot_Contiguous;
397: bv->ops->dotvec = BVDotVec_Contiguous;
398: bv->ops->dotvec_local = BVDotVec_Local_Contiguous;
399: bv->ops->scale = BVScale_Contiguous;
400: bv->ops->norm = BVNorm_Contiguous;
401: bv->ops->norm_local = BVNorm_Local_Contiguous;
402: bv->ops->normalize = BVNormalize_Contiguous;
403: bv->ops->matmult = BVMatMult_Contiguous;
404: bv->ops->copy = BVCopy_Contiguous;
405: bv->ops->copycolumn = BVCopyColumn_Contiguous;
406: bv->ops->resize = BVResize_Contiguous;
407: bv->ops->getcolumn = BVGetColumn_Contiguous;
408: bv->ops->restorecolumn = BVRestoreColumn_Contiguous;
409: bv->ops->getarray = BVGetArray_Contiguous;
410: bv->ops->getarrayread = BVGetArrayRead_Contiguous;
411: bv->ops->getmat = BVGetMat_Default;
412: bv->ops->restoremat = BVRestoreMat_Default;
413: bv->ops->destroy = BVDestroy_Contiguous;
414: bv->ops->view = BVView_Contiguous;
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }