Actual source code: bvmat.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, 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: {
 24:   BV_MAT         *y = (BV_MAT*)Y->data,*x = (BV_MAT*)X->data;
 25:   PetscScalar    *px,*py,*q;
 26:   PetscInt       ldq;

 29:   MatDenseGetArray(x->A,&px);
 30:   MatDenseGetArray(y->A,&py);
 31:   if (Q) {
 32:     MatGetSize(Q,&ldq,NULL);
 33:     MatDenseGetArray(Q,&q);
 34:     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:     MatDenseRestoreArray(Q,&q);
 36:   } else {
 37:     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);
 38:   }
 39:   MatDenseRestoreArray(x->A,&px);
 40:   MatDenseRestoreArray(y->A,&py);
 41:   return(0);
 42: }

 44: PetscErrorCode BVMultVec_Mat(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
 45: {
 47:   BV_MAT         *x = (BV_MAT*)X->data;
 48:   PetscScalar    *px,*py,*qq=q;

 51:   MatDenseGetArray(x->A,&px);
 52:   VecGetArray(y,&py);
 53:   if (!q) { VecGetArray(X->buffer,&qq); }
 54:   BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,px+(X->nc+X->l)*X->n,qq,beta,py);
 55:   if (!q) { VecRestoreArray(X->buffer,&qq); }
 56:   MatDenseRestoreArray(x->A,&px);
 57:   VecRestoreArray(y,&py);
 58:   return(0);
 59: }

 61: PetscErrorCode BVMultInPlace_Mat(BV V,Mat Q,PetscInt s,PetscInt e)
 62: {
 64:   BV_MAT         *ctx = (BV_MAT*)V->data;
 65:   PetscScalar    *pv,*q;
 66:   PetscInt       ldq;

 69:   MatGetSize(Q,&ldq,NULL);
 70:   MatDenseGetArray(ctx->A,&pv);
 71:   MatDenseGetArray(Q,&q);
 72:   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);
 73:   MatDenseRestoreArray(Q,&q);
 74:   MatDenseRestoreArray(ctx->A,&pv);
 75:   return(0);
 76: }

 78: PetscErrorCode BVMultInPlaceTranspose_Mat(BV V,Mat Q,PetscInt s,PetscInt e)
 79: {
 81:   BV_MAT         *ctx = (BV_MAT*)V->data;
 82:   PetscScalar    *pv,*q;
 83:   PetscInt       ldq;

 86:   MatGetSize(Q,&ldq,NULL);
 87:   MatDenseGetArray(ctx->A,&pv);
 88:   MatDenseGetArray(Q,&q);
 89:   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);
 90:   MatDenseRestoreArray(Q,&q);
 91:   MatDenseRestoreArray(ctx->A,&pv);
 92:   return(0);
 93: }

 95: PetscErrorCode BVDot_Mat(BV X,BV Y,Mat M)
 96: {
 98:   BV_MAT         *x = (BV_MAT*)X->data,*y = (BV_MAT*)Y->data;
 99:   PetscScalar    *px,*py,*m;
100:   PetscInt       ldm;

103:   MatGetSize(M,&ldm,NULL);
104:   MatDenseGetArray(x->A,&px);
105:   MatDenseGetArray(y->A,&py);
106:   MatDenseGetArray(M,&m);
107:   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);
108:   MatDenseRestoreArray(M,&m);
109:   MatDenseRestoreArray(x->A,&px);
110:   MatDenseRestoreArray(y->A,&py);
111:   return(0);
112: }

114: PetscErrorCode BVDotVec_Mat(BV X,Vec y,PetscScalar *q)
115: {
116:   PetscErrorCode    ierr;
117:   BV_MAT            *x = (BV_MAT*)X->data;
118:   PetscScalar       *px,*qq=q;
119:   const PetscScalar *py;
120:   Vec               z = y;

123:   if (X->matrix) {
124:     BV_IPMatMult(X,y);
125:     z = X->Bx;
126:   }
127:   MatDenseGetArray(x->A,&px);
128:   VecGetArrayRead(z,&py);
129:   if (!q) { VecGetArray(X->buffer,&qq); }
130:   BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,qq,x->mpi);
131:   if (!q) { VecRestoreArray(X->buffer,&qq); }
132:   VecRestoreArrayRead(z,&py);
133:   MatDenseRestoreArray(x->A,&px);
134:   return(0);
135: }

137: PetscErrorCode BVDotVec_Local_Mat(BV X,Vec y,PetscScalar *m)
138: {
140:   BV_MAT         *x = (BV_MAT*)X->data;
141:   PetscScalar    *px,*py;
142:   Vec            z = y;

145:   if (X->matrix) {
146:     BV_IPMatMult(X,y);
147:     z = X->Bx;
148:   }
149:   MatDenseGetArray(x->A,&px);
150:   VecGetArray(z,&py);
151:   BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->n,py,m,PETSC_FALSE);
152:   VecRestoreArray(z,&py);
153:   MatDenseRestoreArray(x->A,&px);
154:   return(0);
155: }

157: PetscErrorCode BVScale_Mat(BV bv,PetscInt j,PetscScalar alpha)
158: {
160:   BV_MAT         *ctx = (BV_MAT*)bv->data;
161:   PetscScalar    *array;

164:   MatDenseGetArray(ctx->A,&array);
165:   if (j<0) {
166:     BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->n,array+(bv->nc+bv->l)*bv->n,alpha);
167:   } else {
168:     BVScale_BLAS_Private(bv,bv->n,array+(bv->nc+j)*bv->n,alpha);
169:   }
170:   MatDenseRestoreArray(ctx->A,&array);
171:   return(0);
172: }

174: PetscErrorCode BVNorm_Mat(BV bv,PetscInt j,NormType type,PetscReal *val)
175: {
177:   BV_MAT         *ctx = (BV_MAT*)bv->data;
178:   PetscScalar    *array;

181:   MatDenseGetArray(ctx->A,&array);
182:   if (j<0) {
183:     BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,ctx->mpi);
184:   } else {
185:     BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,ctx->mpi);
186:   }
187:   MatDenseRestoreArray(ctx->A,&array);
188:   return(0);
189: }

191: PetscErrorCode BVNorm_Local_Mat(BV bv,PetscInt j,NormType type,PetscReal *val)
192: {
194:   BV_MAT         *ctx = (BV_MAT*)bv->data;
195:   PetscScalar    *array;

198:   MatDenseGetArray(ctx->A,&array);
199:   if (j<0) {
200:     BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->n,type,val,PETSC_FALSE);
201:   } else {
202:     BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->n,type,val,PETSC_FALSE);
203:   }
204:   MatDenseRestoreArray(ctx->A,&array);
205:   return(0);
206: }

208: PetscErrorCode BVMatMult_Mat(BV V,Mat A,BV W)
209: {
211:   BV_MAT         *v = (BV_MAT*)V->data,*w = (BV_MAT*)W->data;
212:   PetscScalar    *pv,*pw;
213:   PetscInt       j;
214:   PetscBool      flg;
215:   Mat            Vmat,Wmat;

218:   MatHasOperation(A,MATOP_MAT_MULT,&flg);
219:   if (V->vmm && flg) {
220:     BVGetMat(V,&Vmat);
221:     BVGetMat(W,&Wmat);
222:     MatMatMult(A,Vmat,MAT_REUSE_MATRIX,PETSC_DEFAULT,&Wmat);
223:     BVRestoreMat(V,&Vmat);
224:     BVRestoreMat(W,&Wmat);
225:   } else {
226:     MatDenseGetArray(v->A,&pv);
227:     MatDenseGetArray(w->A,&pw);
228:     for (j=0;j<V->k-V->l;j++) {
229:       VecPlaceArray(V->cv[1],pv+(V->nc+V->l+j)*V->n);
230:       VecPlaceArray(W->cv[1],pw+(W->nc+W->l+j)*W->n);
231:       MatMult(A,V->cv[1],W->cv[1]);
232:       VecResetArray(V->cv[1]);
233:       VecResetArray(W->cv[1]);
234:     }
235:     MatDenseRestoreArray(v->A,&pv);
236:     MatDenseRestoreArray(w->A,&pw);
237:   }
238:   return(0);
239: }

241: PetscErrorCode BVCopy_Mat(BV V,BV W)
242: {
244:   BV_MAT         *v = (BV_MAT*)V->data,*w = (BV_MAT*)W->data;
245:   PetscScalar    *pv,*pw,*pvc,*pwc;

248:   MatDenseGetArray(v->A,&pv);
249:   MatDenseGetArray(w->A,&pw);
250:   pvc = pv+(V->nc+V->l)*V->n;
251:   pwc = pw+(W->nc+W->l)*W->n;
252:   PetscMemcpy(pwc,pvc,(V->k-V->l)*V->n*sizeof(PetscScalar));
253:   MatDenseRestoreArray(v->A,&pv);
254:   MatDenseRestoreArray(w->A,&pw);
255:   return(0);
256: }

258: PetscErrorCode BVCopyColumn_Mat(BV V,PetscInt j,PetscInt i)
259: {
261:   BV_MAT         *v = (BV_MAT*)V->data;
262:   PetscScalar    *pv;

265:   MatDenseGetArray(v->A,&pv);
266:   PetscMemcpy(pv+(V->nc+i)*V->n,pv+(V->nc+j)*V->n,V->n*sizeof(PetscScalar));
267:   MatDenseRestoreArray(v->A,&pv);
268:   return(0);
269: }

271: PetscErrorCode BVResize_Mat(BV bv,PetscInt m,PetscBool copy)
272: {
274:   BV_MAT         *ctx = (BV_MAT*)bv->data;
275:   PetscScalar    *pA,*pnew;
276:   PetscInt       lsplit;
277:   Mat            A;
278:   char           str[50];
279:   PetscScalar    *array;
280:   BV             parent;

283:   if (bv->issplit==2) {
284:     parent = bv->splitparent;
285:     lsplit = parent->lsplit;
286:     MatDenseGetArray(((BV_MAT*)parent->data)->A,&array);
287:     MatDensePlaceArray(ctx->A,array+lsplit*bv->n);
288:     MatDenseRestoreArray(((BV_MAT*)parent->data)->A,&array);
289:   } else if (!bv->issplit) {
290:     MatCreateDense(PetscObjectComm((PetscObject)bv->t),bv->n,PETSC_DECIDE,PETSC_DECIDE,m,NULL,&A);
291:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
292:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
293:     PetscLogObjectParent((PetscObject)bv,(PetscObject)A);
294:     if (((PetscObject)bv)->name) {
295:       PetscSNPrintf(str,50,"%s_0",((PetscObject)bv)->name);
296:       PetscObjectSetName((PetscObject)A,str);
297:     }
298:     if (copy) {
299:       MatDenseGetArray(ctx->A,&pA);
300:       MatDenseGetArray(A,&pnew);
301:       PetscMemcpy(pnew,pA,PetscMin(m,bv->m)*bv->n*sizeof(PetscScalar));
302:       MatDenseRestoreArray(ctx->A,&pA);
303:       MatDenseRestoreArray(A,&pnew);
304:     }
305:     MatDestroy(&ctx->A);
306:     ctx->A = A;
307:   }
308:   return(0);
309: }

311: PetscErrorCode BVGetColumn_Mat(BV bv,PetscInt j,Vec *v)
312: {
314:   BV_MAT         *ctx = (BV_MAT*)bv->data;
315:   PetscScalar    *pA;
316:   PetscInt       l;

319:   l = BVAvailableVec;
320:   MatDenseGetArray(ctx->A,&pA);
321:   VecPlaceArray(bv->cv[l],pA+(bv->nc+j)*bv->n);
322:   return(0);
323: }

325: PetscErrorCode BVRestoreColumn_Mat(BV bv,PetscInt j,Vec *v)
326: {
328:   BV_MAT         *ctx = (BV_MAT*)bv->data;
329:   PetscScalar    *pA;
330:   PetscInt       l;

333:   l = (j==bv->ci[0])? 0: 1;
334:   VecResetArray(bv->cv[l]);
335:   MatDenseRestoreArray(ctx->A,&pA);
336:   return(0);
337: }

339: PetscErrorCode BVGetArray_Mat(BV bv,PetscScalar **a)
340: {
342:   BV_MAT         *ctx = (BV_MAT*)bv->data;

345:   MatDenseGetArray(ctx->A,a);
346:   return(0);
347: }

349: PetscErrorCode BVRestoreArray_Mat(BV bv,PetscScalar **a)
350: {
352:   BV_MAT         *ctx = (BV_MAT*)bv->data;

355:   if (a) { MatDenseRestoreArray(ctx->A,a); }
356:   return(0);
357: }

359: PetscErrorCode BVGetArrayRead_Mat(BV bv,const PetscScalar **a)
360: {
362:   BV_MAT         *ctx = (BV_MAT*)bv->data;

365:   MatDenseGetArray(ctx->A,(PetscScalar**)a);
366:   return(0);
367: }

369: PetscErrorCode BVRestoreArrayRead_Mat(BV bv,const PetscScalar **a)
370: {
372:   BV_MAT         *ctx = (BV_MAT*)bv->data;

375:   if (a) { MatDenseRestoreArray(ctx->A,(PetscScalar**)a); }
376:   return(0);
377: }

379: PetscErrorCode BVView_Mat(BV bv,PetscViewer viewer)
380: {
381:   PetscErrorCode    ierr;
382:   BV_MAT            *ctx = (BV_MAT*)bv->data;
383:   PetscViewerFormat format;
384:   PetscBool         isascii;
385:   const char        *bvname,*name;

388:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
389:   if (isascii) {
390:     PetscViewerGetFormat(viewer,&format);
391:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
392:     MatView(ctx->A,viewer);
393:     if (format == PETSC_VIEWER_ASCII_MATLAB) {
394:       PetscObjectGetName((PetscObject)bv,&bvname);
395:       PetscObjectGetName((PetscObject)ctx->A,&name);
396:       PetscViewerASCIIPrintf(viewer,"%s=%s;clear %s\n",bvname,name,name);
397:       if (bv->nc) {
398:         PetscViewerASCIIPrintf(viewer,"%s=%s(:,%D:end);\n",bvname,bvname,bv->nc+1);
399:       }
400:     }
401:   } else {
402:     MatView(ctx->A,viewer);
403:   }
404:   return(0);
405: }

407: PetscErrorCode BVDestroy_Mat(BV bv)
408: {
410:   BV_MAT         *ctx = (BV_MAT*)bv->data;

413:   MatDestroy(&ctx->A);
414:   VecDestroy(&bv->cv[0]);
415:   VecDestroy(&bv->cv[1]);
416:   PetscFree(bv->data);
417:   return(0);
418: }

420: SLEPC_EXTERN PetscErrorCode BVCreate_Mat(BV bv)
421: {
423:   BV_MAT         *ctx;
424:   PetscInt       nloc,bs,lsplit;
425:   PetscBool      seq;
426:   char           str[50];
427:   PetscScalar    *array,*ptr;
428:   BV             parent;

431:   PetscNewLog(bv,&ctx);
432:   bv->data = (void*)ctx;

434:   PetscObjectTypeCompare((PetscObject)bv->t,VECMPI,&ctx->mpi);
435:   if (!ctx->mpi) {
436:     PetscObjectTypeCompare((PetscObject)bv->t,VECSEQ,&seq);
437:     if (!seq) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Cannot create a BVMAT from a non-standard template vector");
438:   }

440:   VecGetLocalSize(bv->t,&nloc);
441:   VecGetBlockSize(bv->t,&bs);

443:   if (bv->issplit) {
444:     /* split BV: share the memory of the parent BV */
445:     parent = bv->splitparent;
446:     lsplit = parent->lsplit;
447:     MatDenseGetArray(((BV_MAT*)parent->data)->A,&array);
448:     ptr = (bv->issplit==1)? array: array+lsplit*nloc;
449:     MatDenseRestoreArray(((BV_MAT*)parent->data)->A,&array);
450:   } else {
451:     /* regular BV: allocate memory for the BV entries */
452:     ptr = NULL;
453:   }
454:   MatCreateDense(PetscObjectComm((PetscObject)bv->t),nloc,PETSC_DECIDE,PETSC_DECIDE,bv->m,ptr,&ctx->A);
455:   MatAssemblyBegin(ctx->A,MAT_FINAL_ASSEMBLY);
456:   MatAssemblyEnd(ctx->A,MAT_FINAL_ASSEMBLY);
457:   PetscLogObjectParent((PetscObject)bv,(PetscObject)ctx->A);
458:   if (((PetscObject)bv)->name) {
459:     PetscSNPrintf(str,50,"%s_0",((PetscObject)bv)->name);
460:     PetscObjectSetName((PetscObject)ctx->A,str);
461:   }

463:   if (bv->Acreate) {
464:     MatCopy(bv->Acreate,ctx->A,SAME_NONZERO_PATTERN);
465:     MatDestroy(&bv->Acreate);
466:   }

468:   VecDuplicateEmpty(bv->t,&bv->cv[0]);
469:   VecDuplicateEmpty(bv->t,&bv->cv[1]);

471:   bv->ops->mult             = BVMult_Mat;
472:   bv->ops->multvec          = BVMultVec_Mat;
473:   bv->ops->multinplace      = BVMultInPlace_Mat;
474:   bv->ops->multinplacetrans = BVMultInPlaceTranspose_Mat;
475:   bv->ops->dot              = BVDot_Mat;
476:   bv->ops->dotvec           = BVDotVec_Mat;
477:   bv->ops->dotvec_local     = BVDotVec_Local_Mat;
478:   bv->ops->scale            = BVScale_Mat;
479:   bv->ops->norm             = BVNorm_Mat;
480:   bv->ops->norm_local       = BVNorm_Local_Mat;
481:   bv->ops->matmult          = BVMatMult_Mat;
482:   bv->ops->copy             = BVCopy_Mat;
483:   bv->ops->copycolumn       = BVCopyColumn_Mat;
484:   bv->ops->resize           = BVResize_Mat;
485:   bv->ops->getcolumn        = BVGetColumn_Mat;
486:   bv->ops->restorecolumn    = BVRestoreColumn_Mat;
487:   bv->ops->getarray         = BVGetArray_Mat;
488:   bv->ops->restorearray     = BVRestoreArray_Mat;
489:   bv->ops->getarrayread     = BVGetArrayRead_Mat;
490:   bv->ops->restorearrayread = BVRestoreArrayRead_Mat;
491:   bv->ops->destroy          = BVDestroy_Mat;
492:   if (!ctx->mpi) bv->ops->view = BVView_Mat;
493:   return(0);
494: }