Line data    Source code

       1             : /*
       2             :    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
       3             :    SLEPc - Scalable Library for Eigenvalue Problem Computations
       4             :    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
       5             : 
       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 a single Vec
      12             : */
      13             : 
      14             : #include <slepc/private/bvimpl.h>
      15             : #include "svec.h"
      16             : 
      17          98 : static PetscErrorCode BVMult_Svec(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
      18             : {
      19          98 :   BV_SVEC           *y = (BV_SVEC*)Y->data,*x = (BV_SVEC*)X->data;
      20          98 :   const PetscScalar *px,*q;
      21          98 :   PetscScalar       *py;
      22          98 :   PetscInt          ldq;
      23             : 
      24          98 :   PetscFunctionBegin;
      25          98 :   PetscCall(VecGetArrayRead(x->v,&px));
      26          98 :   PetscCall(VecGetArray(y->v,&py));
      27          98 :   if (Q) {
      28          78 :     PetscCall(MatDenseGetLDA(Q,&ldq));
      29          78 :     PetscCall(MatDenseGetArrayRead(Q,&q));
      30          78 :     PetscCall(BVMult_BLAS_Private(Y,Y->n,Y->k-Y->l,X->k-X->l,alpha,px+(X->nc+X->l)*X->ld,X->ld,q+Y->l*ldq+X->l,ldq,beta,py+(Y->nc+Y->l)*Y->ld,Y->ld));
      31          78 :     PetscCall(MatDenseRestoreArrayRead(Q,&q));
      32          20 :   } else PetscCall(BVAXPY_BLAS_Private(Y,Y->n,Y->k-Y->l,alpha,px+(X->nc+X->l)*X->ld,X->ld,beta,py+(Y->nc+Y->l)*Y->ld,Y->ld));
      33          98 :   PetscCall(VecRestoreArrayRead(x->v,&px));
      34          98 :   PetscCall(VecRestoreArray(y->v,&py));
      35          98 :   PetscFunctionReturn(PETSC_SUCCESS);
      36             : }
      37             : 
      38         357 : static PetscErrorCode BVMultVec_Svec(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
      39             : {
      40         357 :   BV_SVEC        *x = (BV_SVEC*)X->data;
      41         357 :   PetscScalar    *px,*py,*qq=q;
      42             : 
      43         357 :   PetscFunctionBegin;
      44         357 :   PetscCall(VecGetArray(x->v,&px));
      45         357 :   PetscCall(VecGetArray(y,&py));
      46         357 :   if (!q) PetscCall(VecGetArray(X->buffer,&qq));
      47         357 :   PetscCall(BVMultVec_BLAS_Private(X,X->n,X->k-X->l,alpha,px+(X->nc+X->l)*X->ld,X->ld,qq,beta,py));
      48         357 :   if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
      49         357 :   PetscCall(VecRestoreArray(x->v,&px));
      50         357 :   PetscCall(VecRestoreArray(y,&py));
      51         357 :   PetscFunctionReturn(PETSC_SUCCESS);
      52             : }
      53             : 
      54          46 : static PetscErrorCode BVMultInPlace_Svec(BV V,Mat Q,PetscInt s,PetscInt e)
      55             : {
      56          46 :   BV_SVEC           *ctx = (BV_SVEC*)V->data;
      57          46 :   PetscScalar       *pv;
      58          46 :   const PetscScalar *q;
      59          46 :   PetscInt          ldq;
      60             : 
      61          46 :   PetscFunctionBegin;
      62          46 :   if (s>=e || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
      63          46 :   PetscCall(MatDenseGetLDA(Q,&ldq));
      64          46 :   PetscCall(VecGetArray(ctx->v,&pv));
      65          46 :   PetscCall(MatDenseGetArrayRead(Q,&q));
      66          46 :   PetscCall(BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,s-V->l,e-V->l,pv+(V->nc+V->l)*V->ld,V->ld,q+V->l*ldq+V->l,ldq,PETSC_FALSE));
      67          46 :   PetscCall(MatDenseRestoreArrayRead(Q,&q));
      68          46 :   PetscCall(VecRestoreArray(ctx->v,&pv));
      69          46 :   PetscFunctionReturn(PETSC_SUCCESS);
      70             : }
      71             : 
      72           1 : static PetscErrorCode BVMultInPlaceHermitianTranspose_Svec(BV V,Mat Q,PetscInt s,PetscInt e)
      73             : {
      74           1 :   BV_SVEC           *ctx = (BV_SVEC*)V->data;
      75           1 :   PetscScalar       *pv;
      76           1 :   const PetscScalar *q;
      77           1 :   PetscInt          ldq;
      78             : 
      79           1 :   PetscFunctionBegin;
      80           1 :   if (s>=e || !V->n) PetscFunctionReturn(PETSC_SUCCESS);
      81           1 :   PetscCall(MatDenseGetLDA(Q,&ldq));
      82           1 :   PetscCall(VecGetArray(ctx->v,&pv));
      83           1 :   PetscCall(MatDenseGetArrayRead(Q,&q));
      84           1 :   PetscCall(BVMultInPlace_BLAS_Private(V,V->n,V->k-V->l,s-V->l,e-V->l,pv+(V->nc+V->l)*V->ld,V->ld,q+V->l*ldq+V->l,ldq,PETSC_TRUE));
      85           1 :   PetscCall(MatDenseRestoreArrayRead(Q,&q));
      86           1 :   PetscCall(VecRestoreArray(ctx->v,&pv));
      87           1 :   PetscFunctionReturn(PETSC_SUCCESS);
      88             : }
      89             : 
      90         241 : static PetscErrorCode BVDot_Svec(BV X,BV Y,Mat M)
      91             : {
      92         241 :   BV_SVEC           *x = (BV_SVEC*)X->data,*y = (BV_SVEC*)Y->data;
      93         241 :   const PetscScalar *px,*py;
      94         241 :   PetscScalar       *m;
      95         241 :   PetscInt          ldm;
      96             : 
      97         241 :   PetscFunctionBegin;
      98         241 :   PetscCall(MatDenseGetLDA(M,&ldm));
      99         241 :   PetscCall(VecGetArrayRead(x->v,&px));
     100         241 :   PetscCall(VecGetArrayRead(y->v,&py));
     101         241 :   PetscCall(MatDenseGetArray(M,&m));
     102         241 :   PetscCall(BVDot_BLAS_Private(X,Y->k-Y->l,X->k-X->l,X->n,py+(Y->nc+Y->l)*Y->ld,Y->ld,px+(X->nc+X->l)*X->ld,X->ld,m+X->l*ldm+Y->l,ldm,x->mpi));
     103         241 :   PetscCall(MatDenseRestoreArray(M,&m));
     104         241 :   PetscCall(VecRestoreArrayRead(x->v,&px));
     105         241 :   PetscCall(VecRestoreArrayRead(y->v,&py));
     106         241 :   PetscFunctionReturn(PETSC_SUCCESS);
     107             : }
     108             : 
     109         415 : static PetscErrorCode BVDotVec_Svec(BV X,Vec y,PetscScalar *q)
     110             : {
     111         415 :   BV_SVEC           *x = (BV_SVEC*)X->data;
     112         415 :   const PetscScalar *px,*py;
     113         415 :   PetscScalar       *qq=q;
     114         415 :   Vec               z = y;
     115             : 
     116         415 :   PetscFunctionBegin;
     117         415 :   if (PetscUnlikely(X->matrix)) {
     118          62 :     PetscCall(BV_IPMatMult(X,y));
     119          62 :     z = X->Bx;
     120             :   }
     121         415 :   PetscCall(VecGetArrayRead(x->v,&px));
     122         415 :   PetscCall(VecGetArrayRead(z,&py));
     123         415 :   if (!q) PetscCall(VecGetArray(X->buffer,&qq));
     124         415 :   PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->ld,X->ld,py,qq,x->mpi));
     125         415 :   if (!q) PetscCall(VecRestoreArray(X->buffer,&qq));
     126         415 :   PetscCall(VecRestoreArrayRead(z,&py));
     127         415 :   PetscCall(VecRestoreArrayRead(x->v,&px));
     128         415 :   PetscFunctionReturn(PETSC_SUCCESS);
     129             : }
     130             : 
     131           4 : static PetscErrorCode BVDotVec_Local_Svec(BV X,Vec y,PetscScalar *m)
     132             : {
     133           4 :   BV_SVEC        *x = (BV_SVEC*)X->data;
     134           4 :   PetscScalar    *px,*py;
     135           4 :   Vec            z = y;
     136             : 
     137           4 :   PetscFunctionBegin;
     138           4 :   if (PetscUnlikely(X->matrix)) {
     139           0 :     PetscCall(BV_IPMatMult(X,y));
     140           0 :     z = X->Bx;
     141             :   }
     142           4 :   PetscCall(VecGetArray(x->v,&px));
     143           4 :   PetscCall(VecGetArray(z,&py));
     144           4 :   PetscCall(BVDotVec_BLAS_Private(X,X->n,X->k-X->l,px+(X->nc+X->l)*X->ld,X->ld,py,m,PETSC_FALSE));
     145           4 :   PetscCall(VecRestoreArray(z,&py));
     146           4 :   PetscCall(VecRestoreArray(x->v,&px));
     147           4 :   PetscFunctionReturn(PETSC_SUCCESS);
     148             : }
     149             : 
     150         331 : static PetscErrorCode BVScale_Svec(BV bv,PetscInt j,PetscScalar alpha)
     151             : {
     152         331 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     153         331 :   PetscScalar    *array;
     154             : 
     155         331 :   PetscFunctionBegin;
     156         331 :   if (!bv->n) PetscFunctionReturn(PETSC_SUCCESS);
     157         331 :   PetscCall(VecGetArray(ctx->v,&array));
     158         331 :   if (PetscUnlikely(j<0)) PetscCall(BVScale_BLAS_Private(bv,(bv->k-bv->l)*bv->ld,array+(bv->nc+bv->l)*bv->ld,alpha));
     159         327 :   else PetscCall(BVScale_BLAS_Private(bv,bv->n,array+(bv->nc+j)*bv->ld,alpha));
     160         331 :   PetscCall(VecRestoreArray(ctx->v,&array));
     161         331 :   PetscFunctionReturn(PETSC_SUCCESS);
     162             : }
     163             : 
     164         109 : static PetscErrorCode BVNorm_Svec(BV bv,PetscInt j,NormType type,PetscReal *val)
     165             : {
     166         109 :   BV_SVEC           *ctx = (BV_SVEC*)bv->data;
     167         109 :   const PetscScalar *array;
     168             : 
     169         109 :   PetscFunctionBegin;
     170         109 :   PetscCall(VecGetArrayRead(ctx->v,&array));
     171         109 :   if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->ld,bv->ld,type,val,ctx->mpi));
     172          41 :   else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->ld,bv->ld,type,val,ctx->mpi));
     173         109 :   PetscCall(VecRestoreArrayRead(ctx->v,&array));
     174         109 :   PetscFunctionReturn(PETSC_SUCCESS);
     175             : }
     176             : 
     177           2 : static PetscErrorCode BVNorm_Local_Svec(BV bv,PetscInt j,NormType type,PetscReal *val)
     178             : {
     179           2 :   BV_SVEC           *ctx = (BV_SVEC*)bv->data;
     180           2 :   const PetscScalar *array;
     181             : 
     182           2 :   PetscFunctionBegin;
     183           2 :   PetscCall(VecGetArrayRead(ctx->v,&array));
     184           2 :   if (PetscUnlikely(j<0)) PetscCall(BVNorm_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->ld,bv->ld,type,val,PETSC_FALSE));
     185           2 :   else PetscCall(BVNorm_LAPACK_Private(bv,bv->n,1,array+(bv->nc+j)*bv->ld,bv->ld,type,val,PETSC_FALSE));
     186           2 :   PetscCall(VecRestoreArrayRead(ctx->v,&array));
     187           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     188             : }
     189             : 
     190           3 : static PetscErrorCode BVNormalize_Svec(BV bv,PetscScalar *eigi)
     191             : {
     192           3 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     193           3 :   PetscScalar    *array,*wi=NULL;
     194             : 
     195           3 :   PetscFunctionBegin;
     196           3 :   PetscCall(VecGetArray(ctx->v,&array));
     197           3 :   if (eigi) wi = eigi+bv->l;
     198           3 :   PetscCall(BVNormalize_LAPACK_Private(bv,bv->n,bv->k-bv->l,array+(bv->nc+bv->l)*bv->ld,bv->ld,wi,ctx->mpi));
     199           3 :   PetscCall(VecRestoreArray(ctx->v,&array));
     200           3 :   PetscFunctionReturn(PETSC_SUCCESS);
     201             : }
     202             : 
     203         103 : static PetscErrorCode BVMatMult_Svec(BV V,Mat A,BV W)
     204             : {
     205         103 :   PetscInt       j;
     206         103 :   Mat            Vmat,Wmat;
     207         103 :   Vec            vv,ww;
     208             : 
     209         103 :   PetscFunctionBegin;
     210         103 :   if (V->vmm) {
     211          96 :     PetscCall(BVGetMat(V,&Vmat));
     212          96 :     PetscCall(BVGetMat(W,&Wmat));
     213          96 :     PetscCall(MatProductCreateWithMat(A,Vmat,NULL,Wmat));
     214          96 :     PetscCall(MatProductSetType(Wmat,MATPRODUCT_AB));
     215          96 :     PetscCall(MatProductSetFromOptions(Wmat));
     216          96 :     PetscCall(MatProductSymbolic(Wmat));
     217          96 :     PetscCall(MatProductNumeric(Wmat));
     218          96 :     PetscCall(MatProductClear(Wmat));
     219          96 :     PetscCall(BVRestoreMat(V,&Vmat));
     220          96 :     PetscCall(BVRestoreMat(W,&Wmat));
     221             :   } else {
     222          53 :     for (j=0;j<V->k-V->l;j++) {
     223          46 :       PetscCall(BVGetColumn(V,V->l+j,&vv));
     224          46 :       PetscCall(BVGetColumn(W,W->l+j,&ww));
     225          46 :       PetscCall(MatMult(A,vv,ww));
     226          46 :       PetscCall(BVRestoreColumn(V,V->l+j,&vv));
     227          46 :       PetscCall(BVRestoreColumn(W,W->l+j,&ww));
     228             :     }
     229             :   }
     230         103 :   PetscFunctionReturn(PETSC_SUCCESS);
     231             : }
     232             : 
     233          97 : static PetscErrorCode BVCopy_Svec(BV V,BV W)
     234             : {
     235          97 :   BV_SVEC           *v = (BV_SVEC*)V->data,*w = (BV_SVEC*)W->data;
     236          97 :   const PetscScalar *pv;
     237          97 :   PetscScalar       *pw;
     238          97 :   PetscInt          j;
     239             : 
     240          97 :   PetscFunctionBegin;
     241          97 :   PetscCall(VecGetArrayRead(v->v,&pv));
     242          97 :   PetscCall(VecGetArray(w->v,&pw));
     243         726 :   for (j=0;j<V->k-V->l;j++) PetscCall(PetscArraycpy(pw+(W->nc+W->l+j)*W->ld,pv+(V->nc+V->l+j)*V->ld,V->n));
     244          97 :   PetscCall(VecRestoreArrayRead(v->v,&pv));
     245          97 :   PetscCall(VecRestoreArray(w->v,&pw));
     246          97 :   PetscFunctionReturn(PETSC_SUCCESS);
     247             : }
     248             : 
     249           2 : static PetscErrorCode BVCopyColumn_Svec(BV V,PetscInt j,PetscInt i)
     250             : {
     251           2 :   BV_SVEC        *v = (BV_SVEC*)V->data;
     252           2 :   PetscScalar    *pv;
     253             : 
     254           2 :   PetscFunctionBegin;
     255           2 :   PetscCall(VecGetArray(v->v,&pv));
     256           2 :   PetscCall(PetscArraycpy(pv+(V->nc+i)*V->ld,pv+(V->nc+j)*V->ld,V->n));
     257           2 :   PetscCall(VecRestoreArray(v->v,&pv));
     258           2 :   PetscFunctionReturn(PETSC_SUCCESS);
     259             : }
     260             : 
     261          11 : static PetscErrorCode BVResize_Svec(BV bv,PetscInt m,PetscBool copy)
     262             : {
     263          11 :   BV_SVEC           *ctx = (BV_SVEC*)bv->data;
     264          11 :   PetscScalar       *pnew;
     265          11 :   const PetscScalar *pv;
     266          11 :   PetscInt          bs;
     267          11 :   Vec               vnew;
     268          11 :   char              str[50];
     269             : 
     270          11 :   PetscFunctionBegin;
     271          11 :   PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
     272          11 :   PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),&vnew));
     273          11 :   PetscCall(VecSetType(vnew,bv->vtype));
     274          11 :   PetscCall(VecSetSizes(vnew,m*bv->ld,PETSC_DECIDE));
     275          11 :   PetscCall(VecSetBlockSize(vnew,bs));
     276          11 :   if (((PetscObject)bv)->name) {
     277          11 :     PetscCall(PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name));
     278          11 :     PetscCall(PetscObjectSetName((PetscObject)vnew,str));
     279             :   }
     280          11 :   if (copy) {
     281           2 :     PetscCall(VecGetArrayRead(ctx->v,&pv));
     282           2 :     PetscCall(VecGetArray(vnew,&pnew));
     283           2 :     PetscCall(PetscArraycpy(pnew,pv,PetscMin(m,bv->m)*bv->ld));
     284           2 :     PetscCall(VecRestoreArrayRead(ctx->v,&pv));
     285           2 :     PetscCall(VecRestoreArray(vnew,&pnew));
     286             :   }
     287          11 :   PetscCall(VecDestroy(&ctx->v));
     288          11 :   ctx->v = vnew;
     289          11 :   PetscFunctionReturn(PETSC_SUCCESS);
     290             : }
     291             : 
     292        2698 : static PetscErrorCode BVGetColumn_Svec(BV bv,PetscInt j,Vec *v)
     293             : {
     294        2698 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     295        2698 :   PetscScalar    *pv;
     296        2698 :   PetscInt       l;
     297             : 
     298        2698 :   PetscFunctionBegin;
     299        2698 :   l = BVAvailableVec;
     300        2698 :   PetscCall(VecGetArray(ctx->v,&pv));
     301        2698 :   PetscCall(VecPlaceArray(bv->cv[l],pv+(bv->nc+j)*bv->ld));
     302        2698 :   PetscFunctionReturn(PETSC_SUCCESS);
     303             : }
     304             : 
     305        2698 : static PetscErrorCode BVRestoreColumn_Svec(BV bv,PetscInt j,Vec *v)
     306             : {
     307        2698 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     308        2698 :   PetscInt       l;
     309             : 
     310        2698 :   PetscFunctionBegin;
     311        2698 :   l = (j==bv->ci[0])? 0: 1;
     312        2698 :   PetscCall(VecResetArray(bv->cv[l]));
     313        2698 :   PetscCall(VecRestoreArray(ctx->v,NULL));
     314        2698 :   PetscFunctionReturn(PETSC_SUCCESS);
     315             : }
     316             : 
     317         239 : static PetscErrorCode BVGetArray_Svec(BV bv,PetscScalar **a)
     318             : {
     319         239 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     320             : 
     321         239 :   PetscFunctionBegin;
     322         239 :   PetscCall(VecGetArray(ctx->v,a));
     323         239 :   PetscFunctionReturn(PETSC_SUCCESS);
     324             : }
     325             : 
     326         239 : static PetscErrorCode BVRestoreArray_Svec(BV bv,PetscScalar **a)
     327             : {
     328         239 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     329             : 
     330         239 :   PetscFunctionBegin;
     331         239 :   PetscCall(VecRestoreArray(ctx->v,a));
     332         239 :   PetscFunctionReturn(PETSC_SUCCESS);
     333             : }
     334             : 
     335           8 : static PetscErrorCode BVGetArrayRead_Svec(BV bv,const PetscScalar **a)
     336             : {
     337           8 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     338             : 
     339           8 :   PetscFunctionBegin;
     340           8 :   PetscCall(VecGetArrayRead(ctx->v,a));
     341           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     342             : }
     343             : 
     344           8 : static PetscErrorCode BVRestoreArrayRead_Svec(BV bv,const PetscScalar **a)
     345             : {
     346           8 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     347             : 
     348           8 :   PetscFunctionBegin;
     349           8 :   PetscCall(VecRestoreArrayRead(ctx->v,a));
     350           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     351             : }
     352             : 
     353          12 : static PetscErrorCode BVView_Svec(BV bv,PetscViewer viewer)
     354             : {
     355          12 :   PetscInt          j;
     356          12 :   Vec               v;
     357          12 :   PetscViewerFormat format;
     358          12 :   PetscBool         isascii,ismatlab=PETSC_FALSE;
     359          12 :   const char        *bvname,*name;
     360             : 
     361          12 :   PetscFunctionBegin;
     362          12 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     363          12 :   if (isascii) {
     364          12 :     PetscCall(PetscViewerGetFormat(viewer,&format));
     365          12 :     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
     366          10 :     if (format == PETSC_VIEWER_ASCII_MATLAB) ismatlab = PETSC_TRUE;
     367             :   }
     368           0 :   if (ismatlab) {
     369           0 :     PetscCall(PetscObjectGetName((PetscObject)bv,&bvname));
     370           0 :     PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[];\n",bvname));
     371             :   }
     372          48 :   for (j=0;j<bv->m;j++) {
     373          38 :     PetscCall(BVGetColumn(bv,j,&v));
     374          38 :     PetscCall(VecView(v,viewer));
     375          38 :     if (ismatlab) {
     376           0 :       PetscCall(PetscObjectGetName((PetscObject)v,&name));
     377           0 :       PetscCall(PetscViewerASCIIPrintf(viewer,"%s=[%s,%s];clear %s\n",bvname,bvname,name,name));
     378             :     }
     379          38 :     PetscCall(BVRestoreColumn(bv,j,&v));
     380             :   }
     381          10 :   PetscFunctionReturn(PETSC_SUCCESS);
     382             : }
     383             : 
     384         350 : static PetscErrorCode BVDestroy_Svec(BV bv)
     385             : {
     386         350 :   BV_SVEC        *ctx = (BV_SVEC*)bv->data;
     387             : 
     388         350 :   PetscFunctionBegin;
     389         350 :   PetscCall(VecDestroy(&ctx->v));
     390         350 :   PetscCall(VecDestroy(&bv->cv[0]));
     391         350 :   PetscCall(VecDestroy(&bv->cv[1]));
     392         350 :   PetscCall(PetscFree(bv->data));
     393         350 :   bv->cuda = PETSC_FALSE;
     394         350 :   PetscFunctionReturn(PETSC_SUCCESS);
     395             : }
     396             : 
     397         350 : SLEPC_EXTERN PetscErrorCode BVCreate_Svec(BV bv)
     398             : {
     399         350 :   BV_SVEC           *ctx;
     400         350 :   PetscInt          nloc,N,bs,tglobal=0,tlocal,lsplit,j,lda;
     401         350 :   PetscBool         seq,isdense;
     402         350 :   PetscScalar       *vv;
     403         350 :   const PetscScalar *aa,*array,*ptr;
     404         350 :   char              str[50];
     405         350 :   BV                parent;
     406         350 :   Vec               vpar;
     407         350 :   MatType           mtype;
     408             : 
     409         350 :   PetscFunctionBegin;
     410         350 :   PetscCall(PetscNew(&ctx));
     411         350 :   bv->data = (void*)ctx;
     412             : 
     413         350 :   PetscCall(PetscStrcmpAny(bv->vtype,&bv->cuda,VECSEQCUDA,VECMPICUDA,""));
     414         350 :   PetscCall(PetscStrcmpAny(bv->vtype,&bv->hip,VECSEQHIP,VECMPIHIP,""));
     415         350 :   PetscCall(PetscStrcmpAny(bv->vtype,&ctx->mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
     416             : 
     417         350 :   PetscCall(PetscStrcmp(bv->vtype,VECSEQ,&seq));
     418         350 :   PetscCheck(seq || ctx->mpi || bv->cuda || bv->hip,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"BVSVEC does not support the requested vector type: %s",bv->vtype);
     419             : 
     420         350 :   PetscCall(PetscLayoutGetLocalSize(bv->map,&nloc));
     421         350 :   PetscCall(PetscLayoutGetSize(bv->map,&N));
     422         350 :   PetscCall(PetscLayoutGetBlockSize(bv->map,&bs));
     423         350 :   PetscCall(BV_SetDefaultLD(bv,nloc));
     424         350 :   tlocal = bv->m*bv->ld;
     425         350 :   PetscCall(PetscIntMultError(bv->m,N,&tglobal));  /* just to check integer overflow */
     426             : 
     427         350 :   if (PetscUnlikely(bv->issplit)) {
     428             :     /* split BV: create Vec sharing the memory of the parent BV */
     429          92 :     parent = bv->splitparent;
     430          92 :     lsplit = parent->lsplit;
     431          92 :     vpar = ((BV_SVEC*)parent->data)->v;
     432          92 :     if (bv->cuda) {
     433             : #if defined(PETSC_HAVE_CUDA)
     434             :       PetscCall(VecCUDAGetArrayRead(vpar,&array));
     435             :       if (bv->issplit>0) ptr = (bv->issplit==1)? array: array+lsplit*bv->ld;
     436             :       else ptr = (bv->issplit==1)? array: array-lsplit;
     437             :       PetscCall(VecCUDARestoreArrayRead(vpar,&array));
     438             :       if (ctx->mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,PETSC_DECIDE,NULL,&ctx->v));
     439             :       else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,NULL,&ctx->v));
     440             :       PetscCall(VecCUDAPlaceArray(ctx->v,ptr));
     441             : #endif
     442          92 :     } else if (bv->hip) {
     443             : #if defined(PETSC_HAVE_HIP)
     444             :       PetscCall(VecHIPGetArrayRead(vpar,&array));
     445             :       if (bv->issplit>0) ptr = (bv->issplit==1)? array: array+lsplit*bv->ld;
     446             :       else ptr = (bv->issplit==1)? array: array-lsplit;
     447             :       PetscCall(VecHIPRestoreArrayRead(vpar,&array));
     448             :       if (ctx->mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,PETSC_DECIDE,NULL,&ctx->v));
     449             :       else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,NULL,&ctx->v));
     450             :       PetscCall(VecHIPPlaceArray(ctx->v,ptr));
     451             : #endif
     452             :     } else {
     453          92 :       PetscCall(VecGetArrayRead(vpar,&array));
     454          92 :       if (bv->issplit>0) ptr = (bv->issplit==1)? array: array+lsplit*bv->ld;
     455          12 :       else ptr = (bv->issplit==1)? array: array-lsplit;
     456          92 :       PetscCall(VecRestoreArrayRead(vpar,&array));
     457          92 :       if (ctx->mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,PETSC_DECIDE,NULL,&ctx->v));
     458           4 :       else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)bv),bs,tlocal,NULL,&ctx->v));
     459          92 :       PetscCall(VecPlaceArray(ctx->v,ptr));
     460             :     }
     461             :   } else {
     462             :     /* regular BV: create Vec to store the BV entries */
     463         258 :     PetscCall(VecCreate(PetscObjectComm((PetscObject)bv),&ctx->v));
     464         258 :     PetscCall(VecSetType(ctx->v,bv->vtype));
     465         258 :     PetscCall(VecSetSizes(ctx->v,tlocal,PETSC_DECIDE));
     466         258 :     PetscCall(VecSetBlockSize(ctx->v,bs));
     467             :   }
     468         350 :   if (((PetscObject)bv)->name) {
     469         118 :     PetscCall(PetscSNPrintf(str,sizeof(str),"%s_0",((PetscObject)bv)->name));
     470         118 :     PetscCall(PetscObjectSetName((PetscObject)ctx->v,str));
     471             :   }
     472             : 
     473         350 :   if (PetscUnlikely(bv->Acreate)) {
     474           2 :     PetscCall(MatGetType(bv->Acreate,&mtype));
     475           2 :     PetscCall(PetscStrcmpAny(mtype,&isdense,MATSEQDENSE,MATMPIDENSE,""));
     476           2 :     PetscCheck(isdense,PetscObjectComm((PetscObject)bv->Acreate),PETSC_ERR_SUP,"BVSVEC requires a dense matrix in BVCreateFromMat()");
     477           2 :     PetscCall(MatDenseGetArrayRead(bv->Acreate,&aa));
     478           2 :     PetscCall(MatDenseGetLDA(bv->Acreate,&lda));
     479           2 :     PetscCall(VecGetArray(ctx->v,&vv));
     480          18 :     for (j=0;j<bv->m;j++) PetscCall(PetscArraycpy(vv+j*bv->ld,aa+j*lda,bv->n));
     481           2 :     PetscCall(VecRestoreArray(ctx->v,&vv));
     482           2 :     PetscCall(MatDenseRestoreArrayRead(bv->Acreate,&aa));
     483           2 :     PetscCall(MatDestroy(&bv->Acreate));
     484             :   }
     485             : 
     486         350 :   PetscCall(BVCreateVecEmpty(bv,&bv->cv[0]));
     487         350 :   PetscCall(BVCreateVecEmpty(bv,&bv->cv[1]));
     488             : 
     489         350 :   if (bv->cuda) {
     490             : #if defined(PETSC_HAVE_CUDA)
     491             :     bv->ops->mult             = BVMult_Svec_CUDA;
     492             :     bv->ops->multvec          = BVMultVec_Svec_CUDA;
     493             :     bv->ops->multinplace      = BVMultInPlace_Svec_CUDA;
     494             :     bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Svec_CUDA;
     495             :     bv->ops->dot              = BVDot_Svec_CUDA;
     496             :     bv->ops->dotvec           = BVDotVec_Svec_CUDA;
     497             :     bv->ops->dotvec_local     = BVDotVec_Local_Svec_CUDA;
     498             :     bv->ops->scale            = BVScale_Svec_CUDA;
     499             :     bv->ops->norm             = BVNorm_Svec_CUDA;
     500             :     bv->ops->norm_local       = BVNorm_Local_Svec_CUDA;
     501             :     bv->ops->normalize        = BVNormalize_Svec_CUDA;
     502             :     bv->ops->matmult          = BVMatMult_Svec_CUDA;
     503             :     bv->ops->copy             = BVCopy_Svec_CUDA;
     504             :     bv->ops->copycolumn       = BVCopyColumn_Svec_CUDA;
     505             :     bv->ops->resize           = BVResize_Svec_CUDA;
     506             :     bv->ops->getcolumn        = BVGetColumn_Svec_CUDA;
     507             :     bv->ops->restorecolumn    = BVRestoreColumn_Svec_CUDA;
     508             :     bv->ops->restoresplit     = BVRestoreSplit_Svec_CUDA;
     509             :     bv->ops->restoresplitrows = BVRestoreSplitRows_Svec_CUDA;
     510             :     bv->ops->getmat           = BVGetMat_Svec_CUDA;
     511             :     bv->ops->restoremat       = BVRestoreMat_Svec_CUDA;
     512             : #endif
     513         350 :   } else if (bv->hip) {
     514             : #if defined(PETSC_HAVE_HIP)
     515             :     bv->ops->mult             = BVMult_Svec_HIP;
     516             :     bv->ops->multvec          = BVMultVec_Svec_HIP;
     517             :     bv->ops->multinplace      = BVMultInPlace_Svec_HIP;
     518             :     bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Svec_HIP;
     519             :     bv->ops->dot              = BVDot_Svec_HIP;
     520             :     bv->ops->dotvec           = BVDotVec_Svec_HIP;
     521             :     bv->ops->dotvec_local     = BVDotVec_Local_Svec_HIP;
     522             :     bv->ops->scale            = BVScale_Svec_HIP;
     523             :     bv->ops->norm             = BVNorm_Svec_HIP;
     524             :     bv->ops->norm_local       = BVNorm_Local_Svec_HIP;
     525             :     bv->ops->normalize        = BVNormalize_Svec_HIP;
     526             :     bv->ops->matmult          = BVMatMult_Svec_HIP;
     527             :     bv->ops->copy             = BVCopy_Svec_HIP;
     528             :     bv->ops->copycolumn       = BVCopyColumn_Svec_HIP;
     529             :     bv->ops->resize           = BVResize_Svec_HIP;
     530             :     bv->ops->getcolumn        = BVGetColumn_Svec_HIP;
     531             :     bv->ops->restorecolumn    = BVRestoreColumn_Svec_HIP;
     532             :     bv->ops->restoresplit     = BVRestoreSplit_Svec_HIP;
     533             :     bv->ops->restoresplitrows = BVRestoreSplitRows_Svec_HIP;
     534             :     bv->ops->getmat           = BVGetMat_Svec_HIP;
     535             :     bv->ops->restoremat       = BVRestoreMat_Svec_HIP;
     536             : #endif
     537             :   } else {
     538         350 :     bv->ops->mult             = BVMult_Svec;
     539         350 :     bv->ops->multvec          = BVMultVec_Svec;
     540         350 :     bv->ops->multinplace      = BVMultInPlace_Svec;
     541         350 :     bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Svec;
     542         350 :     bv->ops->dot              = BVDot_Svec;
     543         350 :     bv->ops->dotvec           = BVDotVec_Svec;
     544         350 :     bv->ops->dotvec_local     = BVDotVec_Local_Svec;
     545         350 :     bv->ops->scale            = BVScale_Svec;
     546         350 :     bv->ops->norm             = BVNorm_Svec;
     547         350 :     bv->ops->norm_local       = BVNorm_Local_Svec;
     548         350 :     bv->ops->normalize        = BVNormalize_Svec;
     549         350 :     bv->ops->matmult          = BVMatMult_Svec;
     550         350 :     bv->ops->copy             = BVCopy_Svec;
     551         350 :     bv->ops->copycolumn       = BVCopyColumn_Svec;
     552         350 :     bv->ops->resize           = BVResize_Svec;
     553         350 :     bv->ops->getcolumn        = BVGetColumn_Svec;
     554         350 :     bv->ops->restorecolumn    = BVRestoreColumn_Svec;
     555         350 :     bv->ops->getmat           = BVGetMat_Default;
     556         350 :     bv->ops->restoremat       = BVRestoreMat_Default;
     557             :   }
     558         350 :   bv->ops->getarray         = BVGetArray_Svec;
     559         350 :   bv->ops->restorearray     = BVRestoreArray_Svec;
     560         350 :   bv->ops->getarrayread     = BVGetArrayRead_Svec;
     561         350 :   bv->ops->restorearrayread = BVRestoreArrayRead_Svec;
     562         350 :   bv->ops->destroy          = BVDestroy_Svec;
     563         350 :   if (!ctx->mpi) bv->ops->view = BVView_Svec;
     564         350 :   PetscFunctionReturn(PETSC_SUCCESS);
     565             : }