Actual source code: contig.c

slepc-3.14.1 2020-12-08
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2020, 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: {
 25:   BV_CONTIGUOUS  *y = (BV_CONTIGUOUS*)Y->data,*x = (BV_CONTIGUOUS*)X->data;
 26:   PetscScalar    *q;
 27:   PetscInt       ldq;

 30:   if (Q) {
 31:     MatGetSize(Q,&ldq,NULL);
 32:     MatDenseGetArray(Q,&q);
 33:     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);
 34:     MatDenseRestoreArray(Q,&q);
 35:   } else {
 36:     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);
 37:   }
 38:   return(0);
 39: }

 41: PetscErrorCode BVMultVec_Contiguous(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
 42: {
 44:   BV_CONTIGUOUS  *x = (BV_CONTIGUOUS*)X->data;
 45:   PetscScalar    *py,*qq=q;

 48:   VecGetArray(y,&py);
 49:   if (!q) { VecGetArray(X->buffer,&qq); }
 50:   BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,x->array+(X->nc+X->l)*X->n,qq,beta,py);
 51:   if (!q) { VecRestoreArray(X->buffer,&qq); }
 52:   VecRestoreArray(y,&py);
 53:   return(0);
 54: }

 56: PetscErrorCode BVMultInPlace_Contiguous(BV V,Mat Q,PetscInt s,PetscInt e)
 57: {
 59:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)V->data;
 60:   PetscScalar    *q;
 61:   PetscInt       ldq;

 64:   MatGetSize(Q,&ldq,NULL);
 65:   MatDenseGetArray(Q,&q);
 66:   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);
 67:   MatDenseRestoreArray(Q,&q);
 68:   return(0);
 69: }

 71: PetscErrorCode BVMultInPlaceTranspose_Contiguous(BV V,Mat Q,PetscInt s,PetscInt e)
 72: {
 74:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)V->data;
 75:   PetscScalar    *q;
 76:   PetscInt       ldq;

 79:   MatGetSize(Q,&ldq,NULL);
 80:   MatDenseGetArray(Q,&q);
 81:   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);
 82:   MatDenseRestoreArray(Q,&q);
 83:   return(0);
 84: }

 86: PetscErrorCode BVDot_Contiguous(BV X,BV Y,Mat M)
 87: {
 89:   BV_CONTIGUOUS  *x = (BV_CONTIGUOUS*)X->data,*y = (BV_CONTIGUOUS*)Y->data;
 90:   PetscScalar    *m;
 91:   PetscInt       ldm;

 94:   MatGetSize(M,&ldm,NULL);
 95:   MatDenseGetArray(M,&m);
 96:   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);
 97:   MatDenseRestoreArray(M,&m);
 98:   return(0);
 99: }

101: PetscErrorCode BVDotVec_Contiguous(BV X,Vec y,PetscScalar *q)
102: {
103:   PetscErrorCode    ierr;
104:   BV_CONTIGUOUS     *x = (BV_CONTIGUOUS*)X->data;
105:   const PetscScalar *py;
106:   PetscScalar       *qq=q;
107:   Vec               z = y;

110:   if (X->matrix) {
111:     BV_IPMatMult(X,y);
112:     z = X->Bx;
113:   }
114:   VecGetArrayRead(z,&py);
115:   if (!q) { VecGetArray(X->buffer,&qq); }
116:   BVDotVec_BLAS_Private(X,X->n,X->k-X->l,x->array+(X->nc+X->l)*X->n,py,qq,x->mpi);
117:   if (!q) { VecRestoreArray(X->buffer,&qq); }
118:   VecRestoreArrayRead(z,&py);
119:   return(0);
120: }

122: PetscErrorCode BVDotVec_Local_Contiguous(BV X,Vec y,PetscScalar *m)
123: {
125:   BV_CONTIGUOUS  *x = (BV_CONTIGUOUS*)X->data;
126:   PetscScalar    *py;
127:   Vec            z = y;

130:   if (X->matrix) {
131:     BV_IPMatMult(X,y);
132:     z = X->Bx;
133:   }
134:   VecGetArray(z,&py);
135:   BVDotVec_BLAS_Private(X,X->n,X->k-X->l,x->array+(X->nc+X->l)*X->n,py,m,PETSC_FALSE);
136:   VecRestoreArray(z,&py);
137:   return(0);
138: }

140: PetscErrorCode BVScale_Contiguous(BV bv,PetscInt j,PetscScalar alpha)
141: {
143:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;

146:   if (j<0) {
147:     BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->n,ctx->array+(bv->nc+bv->l)*bv->n,alpha);
148:   } else {
149:     BVScale_BLAS_Private(bv,bv->n,ctx->array+(bv->nc+j)*bv->n,alpha);
150:   }
151:   return(0);
152: }

154: PetscErrorCode BVNorm_Contiguous(BV bv,PetscInt j,NormType type,PetscReal *val)
155: {
157:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;

160:   if (j<0) {
161:     BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,type,val,ctx->mpi);
162:   } else {
163:     BVNorm_LAPACK_Private(bv,bv->n,1,ctx->array+(bv->nc+j)*bv->n,type,val,ctx->mpi);
164:   }
165:   return(0);
166: }

168: PetscErrorCode BVNorm_Local_Contiguous(BV bv,PetscInt j,NormType type,PetscReal *val)
169: {
171:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;

174:   if (j<0) {
175:     BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,type,val,PETSC_FALSE);
176:   } else {
177:     BVNorm_LAPACK_Private(bv,bv->n,1,ctx->array+(bv->nc+j)*bv->n,type,val,PETSC_FALSE);
178:   }
179:   return(0);
180: }

182: PetscErrorCode BVNormalize_Contiguous(BV bv,PetscScalar *eigi)
183: {
185:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;
186:   PetscScalar    *wi=NULL;

189:   if (eigi) wi = eigi+bv->l;
190:   BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,ctx->array+(bv->nc+bv->l)*bv->n,wi,ctx->mpi);
191:   return(0);
192: }

194: PetscErrorCode BVMatMult_Contiguous(BV V,Mat A,BV W)
195: {
197:   BV_CONTIGUOUS  *v = (BV_CONTIGUOUS*)V->data,*w = (BV_CONTIGUOUS*)W->data;
198:   PetscInt       j;
199:   Mat            Vmat,Wmat;

202:   if (V->vmm) {
203:     BVGetMat(V,&Vmat);
204:     BVGetMat(W,&Wmat);
205:     MatProductCreateWithMat(A,Vmat,NULL,Wmat);
206:     MatProductSetType(Wmat,MATPRODUCT_AB);
207:     MatProductSetFromOptions(Wmat);
208:     MatProductSymbolic(Wmat);
209:     MatProductNumeric(Wmat);
210:     MatProductClear(Wmat);
211:     BVRestoreMat(V,&Vmat);
212:     BVRestoreMat(W,&Wmat);
213:   } else {
214:     for (j=0;j<V->k-V->l;j++) {
215:       MatMult(A,v->V[V->nc+V->l+j],w->V[W->nc+W->l+j]);
216:     }
217:   }
218:   return(0);
219: }

221: PetscErrorCode BVCopy_Contiguous(BV V,BV W)
222: {
224:   BV_CONTIGUOUS  *v = (BV_CONTIGUOUS*)V->data,*w = (BV_CONTIGUOUS*)W->data;
225:   PetscScalar    *pvc,*pwc;

228:   pvc = v->array+(V->nc+V->l)*V->n;
229:   pwc = w->array+(W->nc+W->l)*W->n;
230:   PetscArraycpy(pwc,pvc,(V->k-V->l)*V->n);
231:   return(0);
232: }

234: PetscErrorCode BVCopyColumn_Contiguous(BV V,PetscInt j,PetscInt i)
235: {
237:   BV_CONTIGUOUS  *v = (BV_CONTIGUOUS*)V->data;

240:   PetscArraycpy(v->array+(V->nc+i)*V->n,v->array+(V->nc+j)*V->n,V->n);
241:   return(0);
242: }

244: PetscErrorCode BVResize_Contiguous(BV bv,PetscInt m,PetscBool copy)
245: {
247:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;
248:   PetscInt       j,bs;
249:   PetscScalar    *newarray;
250:   Vec            *newV;
251:   char           str[50];

254:   VecGetBlockSize(bv->t,&bs);
255:   PetscMalloc1(m*bv->n,&newarray);
256:   PetscArrayzero(newarray,m*bv->n);
257:   PetscMalloc1(m,&newV);
258:   for (j=0;j<m;j++) {
259:     if (ctx->mpi) {
260:       VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,PETSC_DECIDE,newarray+j*bv->n,newV+j);
261:     } else {
262:       VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv->t),bs,bv->n,newarray+j*bv->n,newV+j);
263:     }
264:   }
265:   PetscLogObjectParents(bv,m,newV);
266:   if (((PetscObject)bv)->name) {
267:     for (j=0;j<m;j++) {
268:       PetscSNPrintf(str,sizeof(str),"%s_%d",((PetscObject)bv)->name,(int)j);
269:       PetscObjectSetName((PetscObject)newV[j],str);
270:     }
271:   }
272:   if (copy) {
273:     PetscArraycpy(newarray,ctx->array,PetscMin(m,bv->m)*bv->n);
274:   }
275:   VecDestroyVecs(bv->m,&ctx->V);
276:   ctx->V = newV;
277:   PetscFree(ctx->array);
278:   ctx->array = newarray;
279:   return(0);
280: }

282: PetscErrorCode BVGetColumn_Contiguous(BV bv,PetscInt j,Vec *v)
283: {
284:   BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;
285:   PetscInt      l;

288:   l = BVAvailableVec;
289:   bv->cv[l] = ctx->V[bv->nc+j];
290:   return(0);
291: }

293: PetscErrorCode BVGetArray_Contiguous(BV bv,PetscScalar **a)
294: {
295:   BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;

298:   *a = ctx->array;
299:   return(0);
300: }

302: PetscErrorCode BVGetArrayRead_Contiguous(BV bv,const PetscScalar **a)
303: {
304:   BV_CONTIGUOUS *ctx = (BV_CONTIGUOUS*)bv->data;

307:   *a = ctx->array;
308:   return(0);
309: }

311: PetscErrorCode BVDestroy_Contiguous(BV bv)
312: {
314:   BV_CONTIGUOUS  *ctx = (BV_CONTIGUOUS*)bv->data;

317:   if (!bv->issplit) {
318:     VecDestroyVecs(bv->nc+bv->m,&ctx->V);
319:     PetscFree(ctx->array);
320:   }
321:   PetscFree(bv->data);
322:   return(0);
323: }

325: SLEPC_EXTERN PetscErrorCode BVCreate_Contiguous(BV bv)
326: {
328:   BV_CONTIGUOUS  *ctx;
329:   PetscInt       j,nloc,bs,lsplit;
330:   PetscBool      seq;
331:   PetscScalar    *aa;
332:   char           str[50];
333:   PetscScalar    *array;
334:   BV             parent;
335:   Vec            *Vpar;

338:   PetscNewLog(bv,&ctx);
339:   bv->data = (void*)ctx;

341:   PetscObjectTypeCompare((PetscObject)bv->t,VECMPI,&ctx->mpi);
342:   if (!ctx->mpi) {
343:     PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq);
344:     if (!seq) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot create a contiguous BV from a non-standard template vector");
345:   }

347:   VecGetLocalSize(bv->t,&nloc);
348:   VecGetBlockSize(bv->t,&bs);

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

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

385:   bv->ops->mult             = BVMult_Contiguous;
386:   bv->ops->multvec          = BVMultVec_Contiguous;
387:   bv->ops->multinplace      = BVMultInPlace_Contiguous;
388:   bv->ops->multinplacetrans = BVMultInPlaceTranspose_Contiguous;
389:   bv->ops->dot              = BVDot_Contiguous;
390:   bv->ops->dotvec           = BVDotVec_Contiguous;
391:   bv->ops->dotvec_local     = BVDotVec_Local_Contiguous;
392:   bv->ops->scale            = BVScale_Contiguous;
393:   bv->ops->norm             = BVNorm_Contiguous;
394:   bv->ops->norm_local       = BVNorm_Local_Contiguous;
395:   bv->ops->normalize        = BVNormalize_Contiguous;
396:   bv->ops->matmult          = BVMatMult_Contiguous;
397:   bv->ops->copy             = BVCopy_Contiguous;
398:   bv->ops->copycolumn       = BVCopyColumn_Contiguous;
399:   bv->ops->resize           = BVResize_Contiguous;
400:   bv->ops->getcolumn        = BVGetColumn_Contiguous;
401:   bv->ops->getarray         = BVGetArray_Contiguous;
402:   bv->ops->getarrayread     = BVGetArrayRead_Contiguous;
403:   bv->ops->destroy          = BVDestroy_Contiguous;
404:   return(0);
405: }