Actual source code: contig.c

slepc-3.19.2 2023-09-05
Report Typos and Errors
  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:   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,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));
 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->n,beta,y->array+(Y->nc+Y->l)*Y->n));
 35:   PetscFunctionReturn(PETSC_SUCCESS);
 36: }

 38: 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->n,qq,beta,py));
 47:   if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
 48:   PetscCall(VecRestoreArray(y,&py));
 49:   PetscFunctionReturn(PETSC_SUCCESS);
 50: }

 52: 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:   PetscCall(MatDenseGetLDA(Q,&ldq));
 60:   PetscCall(MatDenseGetArrayRead(Q,&q));
 61:   PetscCall(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));
 62:   PetscCall(MatDenseRestoreArrayRead(Q,&q));
 63:   PetscFunctionReturn(PETSC_SUCCESS);
 64: }

 66: PetscErrorCode BVMultInPlaceHermitianTranspose_Contiguous(BV V,Mat Q,PetscInt s,PetscInt e)
 67: {
 68:   BV_CONTIGUOUS     *ctx = (BV_CONTIGUOUS*)V->data;
 69:   const PetscScalar *q;
 70:   PetscInt          ldq;

 72:   PetscFunctionBegin;
 73:   PetscCall(MatDenseGetLDA(Q,&ldq));
 74:   PetscCall(MatDenseGetArrayRead(Q,&q));
 75:   PetscCall(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));
 76:   PetscCall(MatDenseRestoreArrayRead(Q,&q));
 77:   PetscFunctionReturn(PETSC_SUCCESS);
 78: }

 80: PetscErrorCode BVDot_Contiguous(BV X,BV Y,Mat M)
 81: {
 82:   BV_CONTIGUOUS  *x = (BV_CONTIGUOUS*)X->data,*y = (BV_CONTIGUOUS*)Y->data;
 83:   PetscScalar    *m;
 84:   PetscInt       ldm;

 86:   PetscFunctionBegin;
 87:   PetscCall(MatDenseGetLDA(M,&ldm));
 88:   PetscCall(MatDenseGetArray(M,&m));
 89:   PetscCall(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));
 90:   PetscCall(MatDenseRestoreArray(M,&m));
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

 94: PetscErrorCode BVDotVec_Contiguous(BV X,Vec y,PetscScalar *q)
 95: {
 96:   BV_CONTIGUOUS     *x = (BV_CONTIGUOUS*)X->data;
 97:   const PetscScalar *py;
 98:   PetscScalar       *qq=q;
 99:   Vec               z = y;

101:   PetscFunctionBegin;
102:   if (PetscUnlikely(X->matrix)) {
103:     PetscCall(BV_IPMatMult(X,y));
104:     z = X->Bx;
105:   }
106:   PetscCall(VecGetArrayRead(z,&py));
107:   if (!q) PetscCall(VecGetArray(X->buffer,&qq));
108:   PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,x->array+(X->nc+X->l)*X->n,py,qq,x->mpi));
109:   if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
110:   PetscCall(VecRestoreArrayRead(z,&py));
111:   PetscFunctionReturn(PETSC_SUCCESS);
112: }

114: PetscErrorCode BVDotVec_Local_Contiguous(BV X,Vec y,PetscScalar *m)
115: {
116:   BV_CONTIGUOUS  *x = (BV_CONTIGUOUS*)X->data;
117:   PetscScalar    *py;
118:   Vec            z = y;

120:   PetscFunctionBegin;
121:   if (PetscUnlikely(X->matrix)) {
122:     PetscCall(BV_IPMatMult(X,y));
123:     z = X->Bx;
124:   }
125:   PetscCall(VecGetArray(z,&py));
126:   PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,x->array+(X->nc+X->l)*X->n,py,m,PETSC_FALSE));
127:   PetscCall(VecRestoreArray(z,&py));
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: PetscErrorCode BVScale_Contiguous(BV bv,PetscInt j,PetscScalar alpha)
132: {
133:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;

135:   PetscFunctionBegin;
136:   if (PetscUnlikely(j<0)) PetscCall(BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->n,ctx->array+(bv->nc+bv->l)*bv->n,alpha));
137:   else PetscCall(BVScale_BLAS_Private(bv,bv->n,ctx->array+(bv->nc+j)*bv->n,alpha));
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: PetscErrorCode BVNorm_Contiguous(BV bv,PetscInt j,NormType type,PetscReal *val)
142: {
143:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;

145:   PetscFunctionBegin;
146:   if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,type,val,ctx->mpi));
147:   else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,ctx->array+(bv->nc+j)*bv->n,type,val,ctx->mpi));
148:   PetscFunctionReturn(PETSC_SUCCESS);
149: }

151: PetscErrorCode BVNorm_Local_Contiguous(BV bv,PetscInt j,NormType type,PetscReal *val)
152: {
153:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;

155:   PetscFunctionBegin;
156:   if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,type,val,PETSC_FALSE));
157:   else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,ctx->array+(bv->nc+j)*bv->n,type,val,PETSC_FALSE));
158:   PetscFunctionReturn(PETSC_SUCCESS);
159: }

161: PetscErrorCode BVNormalize_Contiguous(BV bv,PetscScalar *eigi)
162: {
163:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;
164:   PetscScalar    *wi=NULL;

166:   PetscFunctionBegin;
167:   if (eigi) wi = eigi+bv->l;
168:   PetscCall(BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,wi,ctx->mpi));
169:   PetscFunctionReturn(PETSC_SUCCESS);
170: }

172: PetscErrorCode BVMatMult_Contiguous(BV V,Mat A,BV W)
173: {
174:   BV_CONTIGUOUS  *v = (BV_CONTIGUOUS*)V->data,*w = (BV_CONTIGUOUS*)W->data;
175:   PetscInt       j;
176:   Mat            Vmat,Wmat;

178:   PetscFunctionBegin;
179:   if (V->vmm) {
180:     PetscCall(BVGetMat(V,&Vmat));
181:     PetscCall(BVGetMat(W,&Wmat));
182:     PetscCall(MatProductCreateWithMat(A,Vmat,NULL,Wmat));
183:     PetscCall(MatProductSetType(Wmat,MATPRODUCT_AB));
184:     PetscCall(MatProductSetFromOptions(Wmat));
185:     PetscCall(MatProductSymbolic(Wmat));
186:     PetscCall(MatProductNumeric(Wmat));
187:     PetscCall(MatProductClear(Wmat));
188:     PetscCall(BVRestoreMat(V,&Vmat));
189:     PetscCall(BVRestoreMat(W,&Wmat));
190:   } else {
191:     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]));
192:   }
193:   PetscFunctionReturn(PETSC_SUCCESS);
194: }

196: PetscErrorCode BVCopy_Contiguous(BV V,BV W)
197: {
198:   BV_CONTIGUOUS  *v = (BV_CONTIGUOUS*)V->data,*w = (BV_CONTIGUOUS*)W->data;
199:   PetscScalar    *pvc,*pwc;

201:   PetscFunctionBegin;
202:   pvc = v->array+(V->nc+V->l)*V->n;
203:   pwc = w->array+(W->nc+W->l)*W->n;
204:   PetscCall(PetscArraycpy(pwc,pvc,(V->k-V->l)*V->n));
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: PetscErrorCode BVCopyColumn_Contiguous(BV V,PetscInt j,PetscInt i)
209: {
210:   BV_CONTIGUOUS  *v = (BV_CONTIGUOUS*)V->data;

212:   PetscFunctionBegin;
213:   PetscCall(PetscArraycpy(v->array+(V->nc+i)*V->n,v->array+(V->nc+j)*V->n,V->n));
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: PetscErrorCode BVResize_Contiguous(BV bv,PetscInt m,PetscBool copy)
218: {
219:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;
220:   PetscInt       j,bs;
221:   PetscScalar    *newarray;
222:   Vec            *newV;
223:   char           str[50];

225:   PetscFunctionBegin;
226:   PetscCall(VecGetBlockSize(bv->t,&bs));
227:   PetscCall(PetscMalloc1(m*bv->n,&newarray));
228:   PetscCall(PetscArrayzero(newarray,m*bv->n));
229:   PetscCall(PetscMalloc1(m,&newV));
230:   for (j=0;j<m;j++) {
231:     if (ctx->mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,PETSC_DECIDE,newarray+j*bv->n,newV+j));
232:     else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,newarray+j*bv->n,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->n));
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: 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: 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: 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: 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: 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: 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;
335:   PetscBool      seq;
336:   PetscScalar    *aa;
337:   char           str[50];
338:   PetscScalar    *array;
339:   BV             parent;
340:   Vec            *Vpar;

342:   PetscFunctionBegin;
343:   PetscCall(PetscNew(&ctx));
344:   bv->data = (void*)ctx;

346:   PetscCall(PetscObjectTypeCompare((PetscObject)bv->t,VECMPI,&ctx->mpi));
347:   if (!ctx->mpi) {
348:     PetscCall(PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq));
349:     PetscCheck(seq,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot create a contiguous BV from a non-standard template vector");
350:   }

352:   PetscCall(VecGetLocalSize(bv->t,&nloc));
353:   PetscCall(VecGetBlockSize(bv->t,&bs));

355:   if (PetscUnlikely(bv->issplit)) {
356:     /* split BV: share memory and Vecs of the parent BV */
357:     parent = bv->splitparent;
358:     lsplit = parent->lsplit;
359:     Vpar   = ((BV_CONTIGUOUS*)parent->data)->V;
360:     ctx->V = (bv->issplit==1)? Vpar: Vpar+lsplit;
361:     array  = ((BV_CONTIGUOUS*)parent->data)->array;
362:     ctx->array = (bv->issplit==1)? array: array+lsplit*nloc;
363:   } else {
364:     /* regular BV: allocate memory and Vecs for the BV entries */
365:     PetscCall(PetscCalloc1(bv->m*nloc,&ctx->array));
366:     PetscCall(PetscMalloc1(bv->m,&ctx->V));
367:     for (j=0;j<bv->m;j++) {
368:       if (ctx->mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,nloc,PETSC_DECIDE,ctx->array+j*nloc,ctx->V+j));
369:       else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,nloc,ctx->array+j*nloc,ctx->V+j));
370:     }
371:   }
372:   if (((PetscObject)bv)->name) {
373:     for (j=0;j<bv->m;j++) {
374:       PetscCall(PetscSNPrintf(str,sizeof(str),"%s_%" PetscInt_FMT,((PetscObject)bv)->name,j));
375:       PetscCall(PetscObjectSetName((PetscObject)ctx->V[j],str));
376:     }
377:   }

379:   if (PetscUnlikely(bv->Acreate)) {
380:     PetscCall(MatDenseGetArray(bv->Acreate,&aa));
381:     PetscCall(PetscArraycpy(ctx->array,aa,bv->m*nloc));
382:     PetscCall(MatDenseRestoreArray(bv->Acreate,&aa));
383:     PetscCall(MatDestroy(&bv->Acreate));
384:   }

386:   bv->ops->mult             = BVMult_Contiguous;
387:   bv->ops->multvec          = BVMultVec_Contiguous;
388:   bv->ops->multinplace      = BVMultInPlace_Contiguous;
389:   bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Contiguous;
390:   bv->ops->dot              = BVDot_Contiguous;
391:   bv->ops->dotvec           = BVDotVec_Contiguous;
392:   bv->ops->dotvec_local     = BVDotVec_Local_Contiguous;
393:   bv->ops->scale            = BVScale_Contiguous;
394:   bv->ops->norm             = BVNorm_Contiguous;
395:   bv->ops->norm_local       = BVNorm_Local_Contiguous;
396:   bv->ops->normalize        = BVNormalize_Contiguous;
397:   bv->ops->matmult          = BVMatMult_Contiguous;
398:   bv->ops->copy             = BVCopy_Contiguous;
399:   bv->ops->copycolumn       = BVCopyColumn_Contiguous;
400:   bv->ops->resize           = BVResize_Contiguous;
401:   bv->ops->getcolumn        = BVGetColumn_Contiguous;
402:   bv->ops->restorecolumn    = BVRestoreColumn_Contiguous;
403:   bv->ops->getarray         = BVGetArray_Contiguous;
404:   bv->ops->getarrayread     = BVGetArrayRead_Contiguous;
405:   bv->ops->getmat           = BVGetMat_Default;
406:   bv->ops->restoremat       = BVRestoreMat_Default;
407:   bv->ops->destroy          = BVDestroy_Contiguous;
408:   bv->ops->view             = BVView_Contiguous;
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }