LCOV - code coverage report
Current view: top level - sys/classes/bv/impls/tensor - bvtensor.c (source / functions) Hit Total Coverage
Test: SLEPc Lines: 436 462 94.4 %
Date: 2024-03-28 00:52:16 Functions: 22 23 95.7 %
Legend: Lines: hit not hit

          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             :    Tensor BV that is represented in compact form as V = (I otimes U) S
      12             : */
      13             : 
      14             : #include <slepc/private/bvimpl.h>      /*I "slepcbv.h" I*/
      15             : #include <slepcblaslapack.h>
      16             : 
      17             : typedef struct {
      18             :   BV          U;        /* first factor */
      19             :   Mat         S;        /* second factor */
      20             :   PetscScalar *qB;      /* auxiliary matrix used in non-standard inner products */
      21             :   PetscScalar *sw;      /* work space */
      22             :   PetscInt    d;        /* degree of the tensor BV */
      23             :   PetscInt    ld;       /* leading dimension of a single block in S */
      24             :   PetscInt    puk;      /* copy of the k value */
      25             :   Vec         u;        /* auxiliary work vector */
      26             : } BV_TENSOR;
      27             : 
      28        1187 : static PetscErrorCode BVMultInPlace_Tensor(BV V,Mat Q,PetscInt s,PetscInt e)
      29             : {
      30        1187 :   BV_TENSOR         *ctx = (BV_TENSOR*)V->data;
      31        1187 :   PetscScalar       *pS;
      32        1187 :   const PetscScalar *q;
      33        1187 :   PetscInt          ldq,lds = ctx->ld*ctx->d;
      34             : 
      35        1187 :   PetscFunctionBegin;
      36        1187 :   PetscCall(MatDenseGetLDA(Q,&ldq));
      37        1187 :   PetscCall(MatDenseGetArray(ctx->S,&pS));
      38        1187 :   PetscCall(MatDenseGetArrayRead(Q,&q));
      39        1187 :   PetscCall(BVMultInPlace_BLAS_Private(V,lds,V->k-V->l,s-V->l,e-V->l,pS+(V->nc+V->l)*lds,lds,q+V->l*ldq+V->l,ldq,PETSC_FALSE));
      40        1187 :   PetscCall(MatDenseRestoreArrayRead(Q,&q));
      41        1187 :   PetscCall(MatDenseRestoreArray(ctx->S,&pS));
      42        1187 :   PetscFunctionReturn(PETSC_SUCCESS);
      43             : }
      44             : 
      45           0 : static PetscErrorCode BVMultInPlaceHermitianTranspose_Tensor(BV V,Mat Q,PetscInt s,PetscInt e)
      46             : {
      47           0 :   BV_TENSOR         *ctx = (BV_TENSOR*)V->data;
      48           0 :   PetscScalar       *pS;
      49           0 :   const PetscScalar *q;
      50           0 :   PetscInt          ldq,lds = ctx->ld*ctx->d;
      51             : 
      52           0 :   PetscFunctionBegin;
      53           0 :   PetscCall(MatDenseGetLDA(Q,&ldq));
      54           0 :   PetscCall(MatDenseGetArray(ctx->S,&pS));
      55           0 :   PetscCall(MatDenseGetArrayRead(Q,&q));
      56           0 :   PetscCall(BVMultInPlace_BLAS_Private(V,lds,V->k-V->l,s-V->l,e-V->l,pS+(V->nc+V->l)*lds,lds,q+V->l*ldq+V->l,ldq,PETSC_TRUE));
      57           0 :   PetscCall(MatDenseRestoreArrayRead(Q,&q));
      58           0 :   PetscCall(MatDenseRestoreArray(ctx->S,&pS));
      59           0 :   PetscFunctionReturn(PETSC_SUCCESS);
      60             : }
      61             : 
      62           8 : static PetscErrorCode BVDot_Tensor(BV X,BV Y,Mat M)
      63             : {
      64           8 :   BV_TENSOR         *x = (BV_TENSOR*)X->data,*y = (BV_TENSOR*)Y->data;
      65           8 :   PetscScalar       *m;
      66           8 :   const PetscScalar *px,*py;
      67           8 :   PetscInt          ldm,lds = x->ld*x->d;
      68             : 
      69           8 :   PetscFunctionBegin;
      70           8 :   PetscCheck(x->U==y->U,PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"BVDot() in BVTENSOR requires that both operands have the same U factor");
      71           8 :   PetscCheck(lds==y->ld*y->d,PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mismatching dimensions ld*d %" PetscInt_FMT " %" PetscInt_FMT,lds,y->ld*y->d);
      72           8 :   PetscCall(MatDenseGetLDA(M,&ldm));
      73           8 :   PetscCall(MatDenseGetArrayRead(x->S,&px));
      74           8 :   PetscCall(MatDenseGetArrayRead(y->S,&py));
      75           8 :   PetscCall(MatDenseGetArray(M,&m));
      76           8 :   PetscCall(BVDot_BLAS_Private(X,Y->k-Y->l,X->k-X->l,lds,py+(Y->nc+Y->l)*lds,lds,px+(X->nc+X->l)*lds,lds,m+X->l*ldm+Y->l,ldm,PETSC_FALSE));
      77           8 :   PetscCall(MatDenseRestoreArray(M,&m));
      78           8 :   PetscCall(MatDenseRestoreArrayRead(x->S,&px));
      79           8 :   PetscCall(MatDenseRestoreArrayRead(y->S,&py));
      80           8 :   PetscFunctionReturn(PETSC_SUCCESS);
      81             : }
      82             : 
      83       13955 : static PetscErrorCode BVScale_Tensor(BV bv,PetscInt j,PetscScalar alpha)
      84             : {
      85       13955 :   BV_TENSOR      *ctx = (BV_TENSOR*)bv->data;
      86       13955 :   PetscScalar    *pS;
      87       13955 :   PetscInt       lds = ctx->ld*ctx->d;
      88             : 
      89       13955 :   PetscFunctionBegin;
      90       13955 :   PetscCall(MatDenseGetArray(ctx->S,&pS));
      91       13955 :   if (PetscUnlikely(j<0)) PetscCall(BVScale_BLAS_Private(bv,(bv->k-bv->l)*lds,pS+(bv->nc+bv->l)*lds,alpha));
      92       13955 :   else PetscCall(BVScale_BLAS_Private(bv,lds,pS+(bv->nc+j)*lds,alpha));
      93       13955 :   PetscCall(MatDenseRestoreArray(ctx->S,&pS));
      94       13955 :   PetscFunctionReturn(PETSC_SUCCESS);
      95             : }
      96             : 
      97           8 : static PetscErrorCode BVNorm_Tensor(BV bv,PetscInt j,NormType type,PetscReal *val)
      98             : {
      99           8 :   BV_TENSOR         *ctx = (BV_TENSOR*)bv->data;
     100           8 :   const PetscScalar *pS;
     101           8 :   PetscInt          lds = ctx->ld*ctx->d;
     102             : 
     103           8 :   PetscFunctionBegin;
     104           8 :   PetscCall(MatDenseGetArrayRead(ctx->S,&pS));
     105           8 :   if (j<0) PetscCall(BVNorm_LAPACK_Private(bv,lds,bv->k-bv->l,pS+(bv->nc+bv->l)*lds,lds,type,val,PETSC_FALSE));
     106           0 :   else PetscCall(BVNorm_LAPACK_Private(bv,lds,1,pS+(bv->nc+j)*lds,lds,type,val,PETSC_FALSE));
     107           8 :   PetscCall(MatDenseRestoreArrayRead(ctx->S,&pS));
     108           8 :   PetscFunctionReturn(PETSC_SUCCESS);
     109             : }
     110             : 
     111        1160 : static PetscErrorCode BVCopyColumn_Tensor(BV V,PetscInt j,PetscInt i)
     112             : {
     113        1160 :   BV_TENSOR      *ctx = (BV_TENSOR*)V->data;
     114        1160 :   PetscScalar    *pS;
     115        1160 :   PetscInt       lds = ctx->ld*ctx->d;
     116             : 
     117        1160 :   PetscFunctionBegin;
     118        1160 :   PetscCall(MatDenseGetArray(ctx->S,&pS));
     119        1160 :   PetscCall(PetscArraycpy(pS+(V->nc+i)*lds,pS+(V->nc+j)*lds,lds));
     120        1160 :   PetscCall(MatDenseRestoreArray(ctx->S,&pS));
     121        1160 :   PetscFunctionReturn(PETSC_SUCCESS);
     122             : }
     123             : 
     124       40451 : static PetscErrorCode BVTensorNormColumn(BV bv,PetscInt j,PetscReal *norm)
     125             : {
     126       40451 :   BV_TENSOR         *ctx = (BV_TENSOR*)bv->data;
     127       40451 :   PetscBLASInt      one=1,lds_;
     128       40451 :   PetscScalar       sone=1.0,szero=0.0,*x,dot;
     129       40451 :   const PetscScalar *S;
     130       40451 :   PetscReal         alpha=1.0,scale=0.0,aval;
     131       40451 :   PetscInt          i,lds,ld=ctx->ld;
     132             : 
     133       40451 :   PetscFunctionBegin;
     134       40451 :   lds = ld*ctx->d;
     135       40451 :   PetscCall(MatDenseGetArrayRead(ctx->S,&S));
     136       40451 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     137       40451 :   if (PetscUnlikely(ctx->qB)) {
     138        4475 :     x = ctx->sw;
     139        4475 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,S+j*lds,&one,&szero,x,&one));
     140        4475 :     dot = PetscRealPart(BLASdot_(&lds_,S+j*lds,&one,x,&one));
     141        4475 :     PetscCall(BV_SafeSqrt(bv,dot,norm));
     142             :   } else {
     143             :     /* Compute *norm = BLASnrm2_(&lds_,S+j*lds,&one); */
     144       35976 :     if (lds==1) *norm = PetscAbsScalar(S[j*lds]);
     145             :     else {
     146     6220196 :       for (i=0;i<lds;i++) {
     147     6184220 :         aval = PetscAbsScalar(S[i+j*lds]);
     148     6184220 :         if (aval!=0.0) {
     149     4536691 :           if (PetscUnlikely(scale<aval)) {
     150      208520 :             alpha = 1.0 + alpha*PetscSqr(scale/aval);
     151      208520 :             scale = aval;
     152     4328171 :           } else alpha += PetscSqr(aval/scale);
     153             :         }
     154             :       }
     155       35976 :       *norm = scale*PetscSqrtReal(alpha);
     156             :     }
     157             :   }
     158       40451 :   PetscCall(MatDenseRestoreArrayRead(ctx->S,&S));
     159       40451 :   PetscFunctionReturn(PETSC_SUCCESS);
     160             : }
     161             : 
     162       26837 : static PetscErrorCode BVOrthogonalizeGS1_Tensor(BV bv,PetscInt k,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onorm,PetscReal *norm)
     163             : {
     164       26837 :   BV_TENSOR         *ctx = (BV_TENSOR*)bv->data;
     165       26837 :   PetscScalar       *pS,*cc,*x,dot,sonem=-1.0,sone=1.0,szero=0.0;
     166       26837 :   PetscInt          i,lds = ctx->ld*ctx->d;
     167       26837 :   PetscBLASInt      lds_,k_,one=1;
     168       26837 :   const PetscScalar *omega;
     169             : 
     170       26837 :   PetscFunctionBegin;
     171       26837 :   PetscCheck(!v,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Orthogonalization against an external vector is not allowed in BVTENSOR");
     172       26837 :   PetscCall(MatDenseGetArray(ctx->S,&pS));
     173       26837 :   if (!c) PetscCall(VecGetArray(bv->buffer,&cc));
     174           0 :   else cc = c;
     175       26837 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     176       26837 :   PetscCall(PetscBLASIntCast(k,&k_));
     177             : 
     178       26837 :   if (onorm) PetscCall(BVTensorNormColumn(bv,k,onorm));
     179             : 
     180       26837 :   if (ctx->qB) x = ctx->sw;
     181       17955 :   else x = pS+k*lds;
     182             : 
     183       26837 :   if (PetscUnlikely(bv->orthog_type==BV_ORTHOG_MGS)) {  /* modified Gram-Schmidt */
     184             : 
     185         116 :     if (PetscUnlikely(bv->indef)) { /* signature */
     186           0 :       PetscCall(VecGetArrayRead(bv->omega,&omega));
     187             :     }
     188        1645 :     for (i=-bv->nc;i<k;i++) {
     189        1529 :       if (which && i>=0 && !which[i]) continue;
     190        1529 :       if (ctx->qB) PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,pS+k*lds,&one,&szero,x,&one));
     191             :       /* c_i = (s_k, s_i) */
     192        1529 :       dot = PetscRealPart(BLASdot_(&lds_,pS+i*lds,&one,x,&one));
     193        1529 :       if (bv->indef) dot /= PetscRealPart(omega[i]);
     194        1529 :       PetscCall(BV_SetValue(bv,i,0,cc,dot));
     195             :       /* s_k = s_k - c_i s_i */
     196        1529 :       dot = -dot;
     197        1529 :       PetscCallBLAS("BLASaxpy",BLASaxpy_(&lds_,&dot,pS+i*lds,&one,pS+k*lds,&one));
     198             :     }
     199         116 :     if (PetscUnlikely(bv->indef)) PetscCall(VecRestoreArrayRead(bv->omega,&omega));
     200             : 
     201             :   } else {  /* classical Gram-Schmidt */
     202       26721 :     if (ctx->qB) PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,pS+k*lds,&one,&szero,x,&one));
     203             : 
     204             :     /* cc = S_{0:k-1}^* s_k */
     205       26721 :     PetscCallBLAS("BLASgemv",BLASgemv_("C",&lds_,&k_,&sone,pS,&lds_,x,&one,&szero,cc,&one));
     206             : 
     207             :     /* s_k = s_k - S_{0:k-1} cc */
     208       26721 :     if (PetscUnlikely(bv->indef)) PetscCall(BV_ApplySignature(bv,k,cc,PETSC_TRUE));
     209       26721 :     PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&k_,&sonem,pS,&lds_,cc,&one,&sone,pS+k*lds,&one));
     210       26721 :     if (PetscUnlikely(bv->indef)) PetscCall(BV_ApplySignature(bv,k,cc,PETSC_FALSE));
     211             :   }
     212             : 
     213       26837 :   if (norm) PetscCall(BVTensorNormColumn(bv,k,norm));
     214       26837 :   PetscCall(BV_AddCoefficients(bv,k,h,cc));
     215       26837 :   PetscCall(MatDenseRestoreArray(ctx->S,&pS));
     216       26837 :   PetscCall(VecRestoreArray(bv->buffer,&cc));
     217       26837 :   PetscFunctionReturn(PETSC_SUCCESS);
     218             : }
     219             : 
     220           8 : static PetscErrorCode BVView_Tensor(BV bv,PetscViewer viewer)
     221             : {
     222           8 :   BV_TENSOR         *ctx = (BV_TENSOR*)bv->data;
     223           8 :   PetscViewerFormat format;
     224           8 :   PetscBool         isascii;
     225           8 :   const char        *bvname,*uname,*sname;
     226             : 
     227           8 :   PetscFunctionBegin;
     228           8 :   PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
     229           8 :   if (isascii) {
     230           8 :     PetscCall(PetscViewerGetFormat(viewer,&format));
     231           8 :     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
     232           8 :       PetscCall(PetscViewerASCIIPrintf(viewer,"number of tensor blocks (degree): %" PetscInt_FMT "\n",ctx->d));
     233           8 :       PetscCall(PetscViewerASCIIPrintf(viewer,"number of columns of U factor: %" PetscInt_FMT "\n",ctx->ld));
     234           8 :       PetscFunctionReturn(PETSC_SUCCESS);
     235             :     }
     236           0 :     PetscCall(BVView(ctx->U,viewer));
     237           0 :     PetscCall(MatView(ctx->S,viewer));
     238           0 :     if (format == PETSC_VIEWER_ASCII_MATLAB) {
     239           0 :       PetscCall(PetscObjectGetName((PetscObject)bv,&bvname));
     240           0 :       PetscCall(PetscObjectGetName((PetscObject)ctx->U,&uname));
     241           0 :       PetscCall(PetscObjectGetName((PetscObject)ctx->S,&sname));
     242           0 :       PetscCall(PetscViewerASCIIPrintf(viewer,"%s=kron(eye(%" PetscInt_FMT "),%s)*%s(:,1:%" PetscInt_FMT ");\n",bvname,ctx->d,uname,sname,bv->k));
     243             :     }
     244             :   } else {
     245           0 :     PetscCall(BVView(ctx->U,viewer));
     246           0 :     PetscCall(MatView(ctx->S,viewer));
     247             :   }
     248           0 :   PetscFunctionReturn(PETSC_SUCCESS);
     249             : }
     250             : 
     251        5572 : static PetscErrorCode BVTensorUpdateMatrix(BV V,PetscInt ini,PetscInt end)
     252             : {
     253        5572 :   BV_TENSOR      *ctx = (BV_TENSOR*)V->data;
     254        5572 :   PetscInt       i,j,r,c,l,k,ld=ctx->ld,lds=ctx->d*ctx->ld;
     255        5572 :   PetscScalar    *qB,*sqB;
     256        5572 :   Vec            u;
     257        5572 :   Mat            A;
     258             : 
     259        5572 :   PetscFunctionBegin;
     260        5572 :   if (!V->matrix) PetscFunctionReturn(PETSC_SUCCESS);
     261        4458 :   l = ctx->U->l; k = ctx->U->k;
     262             :   /* update inner product matrix */
     263        4458 :   if (!ctx->qB) {
     264          34 :     PetscCall(PetscCalloc2(lds*lds,&ctx->qB,lds,&ctx->sw));
     265          34 :     PetscCall(BVCreateVec(ctx->U,&ctx->u));
     266             :   }
     267        4458 :   ctx->U->l = 0;
     268       13374 :   for (r=0;r<ctx->d;r++) {
     269       22290 :     for (c=0;c<=r;c++) {
     270       13374 :       PetscCall(MatNestGetSubMat(V->matrix,r,c,&A));
     271       13374 :       if (A) {
     272        8944 :         qB = ctx->qB+c*ld*lds+r*ld;
     273       18008 :         for (i=ini;i<end;i++) {
     274        9064 :           PetscCall(BVGetColumn(ctx->U,i,&u));
     275        9064 :           PetscCall(MatMult(A,u,ctx->u));
     276        9064 :           ctx->U->k = i+1;
     277        9064 :           PetscCall(BVDotVec(ctx->U,ctx->u,qB+i*lds));
     278        9064 :           PetscCall(BVRestoreColumn(ctx->U,i,&u));
     279      210704 :           for (j=0;j<i;j++) qB[i+j*lds] = PetscConj(qB[j+i*lds]);
     280             :           qB[i*lds+i] = PetscRealPart(qB[i+i*lds]);
     281             :         }
     282        8944 :         if (PetscUnlikely(c!=r)) {
     283          56 :           sqB = ctx->qB+r*ld*lds+c*ld;
     284         796 :           for (i=ini;i<end;i++) for (j=0;j<=i;j++) {
     285         684 :             sqB[i+j*lds] = PetscConj(qB[j+i*lds]);
     286         684 :             sqB[j+i*lds] = qB[j+i*lds];
     287             :           }
     288             :         }
     289             :       }
     290             :     }
     291             :   }
     292        4458 :   ctx->U->l = l; ctx->U->k = k;
     293        4458 :   PetscFunctionReturn(PETSC_SUCCESS);
     294             : }
     295             : 
     296         157 : static PetscErrorCode BVTensorBuildFirstColumn_Tensor(BV V,PetscInt k)
     297             : {
     298         157 :   BV_TENSOR      *ctx = (BV_TENSOR*)V->data;
     299         157 :   PetscInt       i,nq=0;
     300         157 :   PetscScalar    *pS,*omega;
     301         157 :   PetscReal      norm;
     302         157 :   PetscBool      breakdown=PETSC_FALSE;
     303             : 
     304         157 :   PetscFunctionBegin;
     305         157 :   PetscCall(MatDenseGetArray(ctx->S,&pS));
     306         847 :   for (i=0;i<ctx->d;i++) {
     307         690 :     if (i>=k) PetscCall(BVSetRandomColumn(ctx->U,nq));
     308          29 :     else PetscCall(BVCopyColumn(ctx->U,i,nq));
     309         690 :     PetscCall(BVOrthogonalizeColumn(ctx->U,nq,pS+i*ctx->ld,&norm,&breakdown));
     310         690 :     if (!breakdown) {
     311         690 :       PetscCall(BVScaleColumn(ctx->U,nq,1.0/norm));
     312         690 :       pS[nq+i*ctx->ld] = norm;
     313         690 :       nq++;
     314             :     }
     315             :   }
     316         157 :   PetscCall(MatDenseRestoreArray(ctx->S,&pS));
     317         157 :   PetscCheck(nq,PetscObjectComm((PetscObject)V),PETSC_ERR_ARG_SIZ,"Cannot build first column of tensor BV; U should contain k=%" PetscInt_FMT " nonzero columns",k);
     318         157 :   PetscCall(BVTensorUpdateMatrix(V,0,nq));
     319         157 :   PetscCall(BVTensorNormColumn(V,0,&norm));
     320         157 :   PetscCall(BVScale_Tensor(V,0,1.0/norm));
     321         157 :   if (V->indef) {
     322          34 :     PetscCall(BV_AllocateSignature(V));
     323          34 :     PetscCall(VecGetArray(V->omega,&omega));
     324          34 :     omega[0] = (norm<0.0)? -1.0: 1.0;
     325          34 :     PetscCall(VecRestoreArray(V->omega,&omega));
     326             :   }
     327             :   /* set active columns */
     328         157 :   ctx->U->l = 0;
     329         157 :   ctx->U->k = nq;
     330         157 :   PetscFunctionReturn(PETSC_SUCCESS);
     331             : }
     332             : 
     333             : /*@
     334             :    BVTensorBuildFirstColumn - Builds the first column of the tensor basis vectors
     335             :    V from the data contained in the first k columns of U.
     336             : 
     337             :    Collective
     338             : 
     339             :    Input Parameters:
     340             : +  V - the basis vectors context
     341             : -  k - the number of columns of U with relevant information
     342             : 
     343             :    Notes:
     344             :    At most d columns are considered, where d is the degree of the tensor BV.
     345             :    Given V = (I otimes U) S, this function computes the first column of V, that
     346             :    is, it computes the coefficients of the first column of S by orthogonalizing
     347             :    the first d columns of U. If k is less than d (or linearly dependent columns
     348             :    are found) then additional random columns are used.
     349             : 
     350             :    The computed column has unit norm.
     351             : 
     352             :    Level: advanced
     353             : 
     354             : .seealso: BVCreateTensor()
     355             : @*/
     356         157 : PetscErrorCode BVTensorBuildFirstColumn(BV V,PetscInt k)
     357             : {
     358         157 :   PetscFunctionBegin;
     359         157 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     360         628 :   PetscValidLogicalCollectiveInt(V,k,2);
     361         157 :   PetscUseMethod(V,"BVTensorBuildFirstColumn_C",(BV,PetscInt),(V,k));
     362         157 :   PetscFunctionReturn(PETSC_SUCCESS);
     363             : }
     364             : 
     365        1210 : static PetscErrorCode BVTensorCompress_Tensor(BV V,PetscInt newc)
     366             : {
     367        1210 :   BV_TENSOR      *ctx = (BV_TENSOR*)V->data;
     368        1210 :   PetscInt       nwu=0,nnc,nrow,lwa,r,c;
     369        1210 :   PetscInt       i,j,k,n,lds=ctx->ld*ctx->d,deg=ctx->d,lock,cs1=V->k,rs1=ctx->U->k,rk=0,offu;
     370        1210 :   PetscScalar    *S,*M,*Z,*pQ,*SS,*SS2,t,sone=1.0,zero=0.0,mone=-1.0,*p,*tau,*work,*qB,*sqB;
     371        1210 :   PetscReal      *sg,tol,*rwork;
     372        1210 :   PetscBLASInt   ld_,cs1_,rs1_,cs1tdeg,n_,info,lw_,newc_,newctdeg,nnc_,nrow_,nnctdeg,lds_,rk_;
     373        1210 :   Mat            Q,A;
     374             : 
     375        1210 :   PetscFunctionBegin;
     376        1210 :   if (!cs1) PetscFunctionReturn(PETSC_SUCCESS);
     377        1206 :   lwa = 6*ctx->ld*lds+2*cs1;
     378        1206 :   n = PetscMin(rs1,deg*cs1);
     379        1206 :   lock = ctx->U->l;
     380        1206 :   nnc = cs1-lock-newc;
     381        1206 :   nrow = rs1-lock;
     382        1206 :   PetscCall(PetscCalloc6(deg*newc*nnc,&SS,newc*nnc,&SS2,(rs1+lock+newc)*n,&pQ,deg*rs1,&tau,lwa,&work,6*n,&rwork));
     383        1206 :   offu = lock*(rs1+1);
     384        1206 :   M = work+nwu;
     385        1206 :   nwu += rs1*cs1*deg;
     386        1206 :   sg = rwork;
     387        1206 :   Z = work+nwu;
     388        1206 :   nwu += deg*cs1*n;
     389        1206 :   PetscCall(PetscBLASIntCast(n,&n_));
     390        1206 :   PetscCall(PetscBLASIntCast(nnc,&nnc_));
     391        1206 :   PetscCall(PetscBLASIntCast(cs1,&cs1_));
     392        1206 :   PetscCall(PetscBLASIntCast(rs1,&rs1_));
     393        1206 :   PetscCall(PetscBLASIntCast(newc,&newc_));
     394        1206 :   PetscCall(PetscBLASIntCast(newc*deg,&newctdeg));
     395        1206 :   PetscCall(PetscBLASIntCast(nnc*deg,&nnctdeg));
     396        1206 :   PetscCall(PetscBLASIntCast(cs1*deg,&cs1tdeg));
     397        1206 :   PetscCall(PetscBLASIntCast(lwa-nwu,&lw_));
     398        1206 :   PetscCall(PetscBLASIntCast(nrow,&nrow_));
     399        1206 :   PetscCall(PetscBLASIntCast(lds,&lds_));
     400        1206 :   PetscCall(MatDenseGetArray(ctx->S,&S));
     401             : 
     402        1206 :   if (newc>0) {
     403             :     /* truncate columns associated with new converged eigenpairs */
     404        1856 :     for (j=0;j<deg;j++) {
     405        6447 :       for (i=lock;i<lock+newc;i++) PetscCall(PetscArraycpy(M+(i-lock+j*newc)*nrow,S+i*lds+j*ctx->ld+lock,nrow));
     406             :     }
     407         289 :     PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     408             : #if !defined (PETSC_USE_COMPLEX)
     409         289 :     PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&newctdeg,M,&nrow_,sg,pQ+offu,&rs1_,Z,&n_,work+nwu,&lw_,&info));
     410             : #else
     411             :     PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&newctdeg,M,&nrow_,sg,pQ+offu,&rs1_,Z,&n_,work+nwu,&lw_,rwork+n,&info));
     412             : #endif
     413         289 :     SlepcCheckLapackInfo("gesvd",info);
     414         289 :     PetscCall(PetscFPTrapPop());
     415             :     /* SVD has rank min(newc,nrow) */
     416         289 :     rk = PetscMin(newc,nrow);
     417        1790 :     for (i=0;i<rk;i++) {
     418        1501 :       t = sg[i];
     419        1501 :       PetscCallBLAS("BLASscal",BLASscal_(&newctdeg,&t,Z+i,&n_));
     420             :     }
     421        1856 :     for (i=0;i<deg;i++) {
     422        6447 :       for (j=lock;j<lock+newc;j++) {
     423        4880 :         PetscCall(PetscArraycpy(S+j*lds+i*ctx->ld+lock,Z+(newc*i+j-lock)*n,rk));
     424        4880 :         PetscCall(PetscArrayzero(S+j*lds+i*ctx->ld+lock+rk,(ctx->ld-lock-rk)));
     425             :       }
     426             :     }
     427             :     /*
     428             :       update columns associated with non-converged vectors, orthogonalize
     429             :       against pQ so that next M has rank nnc+d-1 instead of nrow+d-1
     430             :     */
     431        1856 :     for (i=0;i<deg;i++) {
     432        1567 :       PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&newc_,&nnc_,&nrow_,&sone,pQ+offu,&rs1_,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_,&zero,SS+i*newc*nnc,&newc_));
     433        1567 :       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&nrow_,&nnc_,&newc_,&mone,pQ+offu,&rs1_,SS+i*newc*nnc,&newc_,&sone,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_));
     434             :       /* repeat orthogonalization step */
     435        1567 :       PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&newc_,&nnc_,&nrow_,&sone,pQ+offu,&rs1_,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_,&zero,SS2,&newc_));
     436        1567 :       PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&nrow_,&nnc_,&newc_,&mone,pQ+offu,&rs1_,SS2,&newc_,&sone,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_));
     437       14755 :       for (j=0;j<newc*nnc;j++) *(SS+i*newc*nnc+j) += SS2[j];
     438             :     }
     439             :   }
     440             : 
     441             :   /* truncate columns associated with non-converged eigenpairs */
     442        6274 :   for (j=0;j<deg;j++) {
     443       51392 :     for (i=lock+newc;i<cs1;i++) PetscCall(PetscArraycpy(M+(i-lock-newc+j*nnc)*nrow,S+i*lds+j*ctx->ld+lock,nrow));
     444             :   }
     445        1206 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     446             : #if !defined (PETSC_USE_COMPLEX)
     447        1206 :   PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&nnctdeg,M,&nrow_,sg,pQ+offu+newc*rs1,&rs1_,Z,&n_,work+nwu,&lw_,&info));
     448             : #else
     449             :   PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&nnctdeg,M,&nrow_,sg,pQ+offu+newc*rs1,&rs1_,Z,&n_,work+nwu,&lw_,rwork+n,&info));
     450             : #endif
     451        1206 :   SlepcCheckLapackInfo("gesvd",info);
     452        1206 :   PetscCall(PetscFPTrapPop());
     453        1206 :   tol = PetscMax(rs1,deg*cs1)*PETSC_MACHINE_EPSILON*sg[0];
     454        1206 :   rk = 0;
     455       26921 :   for (i=0;i<PetscMin(nrow,nnctdeg);i++) if (sg[i]>tol) rk++;
     456        1206 :   rk = PetscMin(nnc+deg-1,rk);
     457             :   /* the SVD has rank (at most) nnc+deg-1 */
     458       16698 :   for (i=0;i<rk;i++) {
     459       15492 :     t = sg[i];
     460       15492 :     PetscCallBLAS("BLASscal",BLASscal_(&nnctdeg,&t,Z+i,&n_));
     461             :   }
     462             :   /* update S */
     463        1206 :   PetscCall(PetscArrayzero(S+cs1*lds,(V->m-cs1)*lds));
     464        1206 :   k = ctx->ld-lock-newc-rk;
     465        6274 :   for (i=0;i<deg;i++) {
     466       51392 :     for (j=lock+newc;j<cs1;j++) {
     467       46324 :       PetscCall(PetscArraycpy(S+j*lds+i*ctx->ld+lock+newc,Z+(nnc*i+j-lock-newc)*n,rk));
     468       46324 :       PetscCall(PetscArrayzero(S+j*lds+i*ctx->ld+lock+newc+rk,k));
     469             :     }
     470             :   }
     471        1206 :   if (newc>0) {
     472        1856 :     for (i=0;i<deg;i++) {
     473        1567 :       p = SS+nnc*newc*i;
     474       10037 :       for (j=lock+newc;j<cs1;j++) {
     475       21658 :         for (k=0;k<newc;k++) S[j*lds+i*ctx->ld+lock+k] = *(p++);
     476             :       }
     477             :     }
     478             :   }
     479             : 
     480             :   /* orthogonalize pQ */
     481        1206 :   rk = rk+newc;
     482        1206 :   PetscCall(PetscBLASIntCast(rk,&rk_));
     483        1206 :   PetscCall(PetscBLASIntCast(cs1-lock,&nnc_));
     484        1206 :   PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
     485        1206 :   PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&nrow_,&rk_,pQ+offu,&rs1_,tau,work+nwu,&lw_,&info));
     486        1206 :   SlepcCheckLapackInfo("geqrf",info);
     487        6274 :   for (i=0;i<deg;i++) {
     488        5068 :     PetscCallBLAS("BLAStrmm",BLAStrmm_("L","U","N","N",&rk_,&nnc_,&sone,pQ+offu,&rs1_,S+lock*lds+lock+i*ctx->ld,&lds_));
     489             :   }
     490        1206 :   PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&nrow_,&rk_,&rk_,pQ+offu,&rs1_,tau,work+nwu,&lw_,&info));
     491        1206 :   SlepcCheckLapackInfo("orgqr",info);
     492        1206 :   PetscCall(PetscFPTrapPop());
     493             : 
     494             :   /* update vectors U(:,idx) = U*Q(:,idx) */
     495        1206 :   rk = rk+lock;
     496        2534 :   for (i=0;i<lock;i++) pQ[i*(1+rs1)] = 1.0;
     497        1206 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,rs1,rk,pQ,&Q));
     498        1206 :   ctx->U->k = rs1;
     499        1206 :   PetscCall(BVMultInPlace(ctx->U,Q,lock,rk));
     500        1206 :   PetscCall(MatDestroy(&Q));
     501             : 
     502        1206 :   if (ctx->qB) {
     503             :    /* update matrix qB */
     504         368 :     PetscCall(PetscBLASIntCast(ctx->ld,&ld_));
     505         368 :     PetscCall(PetscBLASIntCast(rk,&rk_));
     506        1104 :     for (r=0;r<ctx->d;r++) {
     507        1840 :       for (c=0;c<=r;c++) {
     508        1104 :         PetscCall(MatNestGetSubMat(V->matrix,r,c,&A));
     509        1104 :         if (A) {
     510         739 :           qB = ctx->qB+r*ctx->ld+c*ctx->ld*lds;
     511         739 :           PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&rs1_,&rk_,&rs1_,&sone,qB,&lds_,pQ,&rs1_,&zero,work+nwu,&rs1_));
     512         739 :           PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&rk_,&rk_,&rs1_,&sone,pQ,&rs1_,work+nwu,&rs1_,&zero,qB,&lds_));
     513       12536 :           for (i=0;i<rk;i++) {
     514      118997 :             for (j=0;j<i;j++) qB[i+j*lds] = PetscConj(qB[j+i*lds]);
     515       11797 :             qB[i+i*lds] = PetscRealPart(qB[i+i*lds]);
     516             :           }
     517       12054 :           for (i=rk;i<ctx->ld;i++) PetscCall(PetscArrayzero(qB+i*lds,ctx->ld));
     518       12536 :           for (i=0;i<rk;i++) PetscCall(PetscArrayzero(qB+i*lds+rk,(ctx->ld-rk)));
     519         739 :           if (c!=r) {
     520           6 :             sqB = ctx->qB+r*ctx->ld*lds+c*ctx->ld;
     521        2526 :             for (i=0;i<ctx->ld;i++) for (j=0;j<ctx->ld;j++) sqB[i+j*lds] = PetscConj(qB[j+i*lds]);
     522             :           }
     523             :         }
     524             :       }
     525             :     }
     526             :   }
     527             : 
     528             :   /* free work space */
     529        1206 :   PetscCall(PetscFree6(SS,SS2,pQ,tau,work,rwork));
     530        1206 :   PetscCall(MatDenseRestoreArray(ctx->S,&S));
     531             : 
     532             :   /* set active columns */
     533        1206 :   if (newc) ctx->U->l += newc;
     534        1206 :   ctx->U->k = rk;
     535        1206 :   PetscFunctionReturn(PETSC_SUCCESS);
     536             : }
     537             : 
     538             : /*@
     539             :    BVTensorCompress - Updates the U and S factors of the tensor basis vectors
     540             :    object V by means of an SVD, removing redundant information.
     541             : 
     542             :    Collective
     543             : 
     544             :    Input Parameters:
     545             : +  V - the tensor basis vectors context
     546             : -  newc - additional columns to be locked
     547             : 
     548             :    Notes:
     549             :    This function is typically used when restarting Krylov solvers. Truncating a
     550             :    tensor BV V = (I otimes U) S to its leading columns amounts to keeping the
     551             :    leading columns of S. However, to effectively reduce the size of the
     552             :    decomposition, it is necessary to compress it in a way that fewer columns of
     553             :    U are employed. This can be achieved by means of an update that involves the
     554             :    SVD of the low-rank matrix [S_0 S_1 ... S_{d-1}], where S_i are the pieces of S.
     555             : 
     556             :    If newc is nonzero, then newc columns are added to the leading columns of V.
     557             :    This means that the corresponding columns of the U and S factors will remain
     558             :    invariant in subsequent operations.
     559             : 
     560             :    Level: advanced
     561             : 
     562             : .seealso: BVCreateTensor(), BVSetActiveColumns()
     563             : @*/
     564        1210 : PetscErrorCode BVTensorCompress(BV V,PetscInt newc)
     565             : {
     566        1210 :   PetscFunctionBegin;
     567        1210 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     568        4840 :   PetscValidLogicalCollectiveInt(V,newc,2);
     569        1210 :   PetscUseMethod(V,"BVTensorCompress_C",(BV,PetscInt),(V,newc));
     570        1210 :   PetscFunctionReturn(PETSC_SUCCESS);
     571             : }
     572             : 
     573         315 : static PetscErrorCode BVTensorGetDegree_Tensor(BV bv,PetscInt *d)
     574             : {
     575         315 :   BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
     576             : 
     577         315 :   PetscFunctionBegin;
     578         315 :   *d = ctx->d;
     579         315 :   PetscFunctionReturn(PETSC_SUCCESS);
     580             : }
     581             : 
     582             : /*@
     583             :    BVTensorGetDegree - Returns the number of blocks (degree) of the tensor BV.
     584             : 
     585             :    Not Collective
     586             : 
     587             :    Input Parameter:
     588             : .  bv - the basis vectors context
     589             : 
     590             :    Output Parameter:
     591             : .  d - the degree
     592             : 
     593             :    Level: advanced
     594             : 
     595             : .seealso: BVCreateTensor()
     596             : @*/
     597         315 : PetscErrorCode BVTensorGetDegree(BV bv,PetscInt *d)
     598             : {
     599         315 :   PetscFunctionBegin;
     600         315 :   PetscValidHeaderSpecific(bv,BV_CLASSID,1);
     601         315 :   PetscAssertPointer(d,2);
     602         315 :   PetscUseMethod(bv,"BVTensorGetDegree_C",(BV,PetscInt*),(bv,d));
     603         315 :   PetscFunctionReturn(PETSC_SUCCESS);
     604             : }
     605             : 
     606        5415 : static PetscErrorCode BVTensorGetFactors_Tensor(BV V,BV *U,Mat *S)
     607             : {
     608        5415 :   BV_TENSOR *ctx = (BV_TENSOR*)V->data;
     609             : 
     610        5415 :   PetscFunctionBegin;
     611        5415 :   PetscCheck(ctx->puk==-1,PetscObjectComm((PetscObject)V),PETSC_ERR_ORDER,"Previous call to BVTensonGetFactors without a BVTensorRestoreFactors call");
     612        5415 :   ctx->puk = ctx->U->k;
     613        5415 :   if (U) *U = ctx->U;
     614        5415 :   if (S) *S = ctx->S;
     615        5415 :   PetscFunctionReturn(PETSC_SUCCESS);
     616             : }
     617             : 
     618             : /*@C
     619             :    BVTensorGetFactors - Returns the two factors involved in the definition of the
     620             :    tensor basis vectors object, V = (I otimes U) S.
     621             : 
     622             :    Logically Collective
     623             : 
     624             :    Input Parameter:
     625             : .  V - the basis vectors context
     626             : 
     627             :    Output Parameters:
     628             : +  U - the BV factor
     629             : -  S - the Mat factor
     630             : 
     631             :    Notes:
     632             :    The returned factors are references (not copies) of the internal factors,
     633             :    so modifying them will change the tensor BV as well. Some operations of the
     634             :    tensor BV assume that U has orthonormal columns, so if the user modifies U
     635             :    this restriction must be taken into account.
     636             : 
     637             :    The returned factors must not be destroyed. BVTensorRestoreFactors() must
     638             :    be called when they are no longer needed.
     639             : 
     640             :    Pass a NULL vector for any of the arguments that is not needed.
     641             : 
     642             :    Level: advanced
     643             : 
     644             : .seealso: BVTensorRestoreFactors()
     645             : @*/
     646        5415 : PetscErrorCode BVTensorGetFactors(BV V,BV *U,Mat *S)
     647             : {
     648        5415 :   PetscFunctionBegin;
     649        5415 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     650        5415 :   PetscUseMethod(V,"BVTensorGetFactors_C",(BV,BV*,Mat*),(V,U,S));
     651        5415 :   PetscFunctionReturn(PETSC_SUCCESS);
     652             : }
     653             : 
     654        5415 : static PetscErrorCode BVTensorRestoreFactors_Tensor(BV V,BV *U,Mat *S)
     655             : {
     656        5415 :   BV_TENSOR      *ctx = (BV_TENSOR*)V->data;
     657             : 
     658        5415 :   PetscFunctionBegin;
     659        5415 :   PetscCall(PetscObjectStateIncrease((PetscObject)V));
     660        5415 :   if (U) *U = NULL;
     661        5415 :   if (S) *S = NULL;
     662        5415 :   PetscCall(BVTensorUpdateMatrix(V,ctx->puk,ctx->U->k));
     663        5415 :   ctx->puk = -1;
     664        5415 :   PetscFunctionReturn(PETSC_SUCCESS);
     665             : }
     666             : 
     667             : /*@C
     668             :    BVTensorRestoreFactors - Restore the two factors that were obtained with
     669             :    BVTensorGetFactors().
     670             : 
     671             :    Logically Collective
     672             : 
     673             :    Input Parameters:
     674             : +  V - the basis vectors context
     675             : .  U - the BV factor (or NULL)
     676             : -  S - the Mat factor (or NULL)
     677             : 
     678             :    Notes:
     679             :    The arguments must match the corresponding call to BVTensorGetFactors().
     680             : 
     681             :    Level: advanced
     682             : 
     683             : .seealso: BVTensorGetFactors()
     684             : @*/
     685        5415 : PetscErrorCode BVTensorRestoreFactors(BV V,BV *U,Mat *S)
     686             : {
     687        5415 :   PetscFunctionBegin;
     688        5415 :   PetscValidHeaderSpecific(V,BV_CLASSID,1);
     689        5415 :   if (U) PetscValidHeaderSpecific(*U,BV_CLASSID,2);
     690        5415 :   if (S) PetscValidHeaderSpecific(*S,MAT_CLASSID,3);
     691        5415 :   PetscUseMethod(V,"BVTensorRestoreFactors_C",(BV,BV*,Mat*),(V,U,S));
     692        5415 :   PetscFunctionReturn(PETSC_SUCCESS);
     693             : }
     694             : 
     695         165 : static PetscErrorCode BVDestroy_Tensor(BV bv)
     696             : {
     697         165 :   BV_TENSOR      *ctx = (BV_TENSOR*)bv->data;
     698             : 
     699         165 :   PetscFunctionBegin;
     700         165 :   PetscCall(BVDestroy(&ctx->U));
     701         165 :   PetscCall(MatDestroy(&ctx->S));
     702         165 :   if (ctx->u) {
     703          34 :     PetscCall(PetscFree2(ctx->qB,ctx->sw));
     704          34 :     PetscCall(VecDestroy(&ctx->u));
     705             :   }
     706         165 :   PetscCall(PetscFree(bv->data));
     707         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorBuildFirstColumn_C",NULL));
     708         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorCompress_C",NULL));
     709         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetDegree_C",NULL));
     710         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetFactors_C",NULL));
     711         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorRestoreFactors_C",NULL));
     712         165 :   PetscFunctionReturn(PETSC_SUCCESS);
     713             : }
     714             : 
     715         165 : SLEPC_EXTERN PetscErrorCode BVCreate_Tensor(BV bv)
     716             : {
     717         165 :   BV_TENSOR      *ctx;
     718             : 
     719         165 :   PetscFunctionBegin;
     720         165 :   PetscCall(PetscNew(&ctx));
     721         165 :   bv->data = (void*)ctx;
     722         165 :   ctx->puk = -1;
     723             : 
     724         165 :   bv->ops->multinplace      = BVMultInPlace_Tensor;
     725         165 :   bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Tensor;
     726         165 :   bv->ops->dot              = BVDot_Tensor;
     727         165 :   bv->ops->scale            = BVScale_Tensor;
     728         165 :   bv->ops->norm             = BVNorm_Tensor;
     729         165 :   bv->ops->copycolumn       = BVCopyColumn_Tensor;
     730         165 :   bv->ops->gramschmidt      = BVOrthogonalizeGS1_Tensor;
     731         165 :   bv->ops->destroy          = BVDestroy_Tensor;
     732         165 :   bv->ops->view             = BVView_Tensor;
     733             : 
     734         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorBuildFirstColumn_C",BVTensorBuildFirstColumn_Tensor));
     735         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorCompress_C",BVTensorCompress_Tensor));
     736         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetDegree_C",BVTensorGetDegree_Tensor));
     737         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetFactors_C",BVTensorGetFactors_Tensor));
     738         165 :   PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorRestoreFactors_C",BVTensorRestoreFactors_Tensor));
     739         165 :   PetscFunctionReturn(PETSC_SUCCESS);
     740             : }
     741             : 
     742             : /*@
     743             :    BVCreateTensor - Creates a tensor BV that is represented in compact form
     744             :    as V = (I otimes U) S, where U has orthonormal columns.
     745             : 
     746             :    Collective
     747             : 
     748             :    Input Parameters:
     749             : +  U - a basis vectors object
     750             : -  d - the number of blocks (degree) of the tensor BV
     751             : 
     752             :    Output Parameter:
     753             : .  V - the new basis vectors context
     754             : 
     755             :    Notes:
     756             :    The new basis vectors object is V = (I otimes U) S, where otimes denotes
     757             :    the Kronecker product, I is the identity matrix of order d, and S is a
     758             :    sequential matrix allocated internally. This compact representation is
     759             :    used e.g. to represent the Krylov basis generated with the linearization
     760             :    of a matrix polynomial of degree d.
     761             : 
     762             :    The size of V (number of rows) is equal to d times n, where n is the size
     763             :    of U. The dimensions of S are d times m rows and m-d+1 columns, where m is
     764             :    the number of columns of U, so m should be at least d.
     765             : 
     766             :    The communicator of V will be the same as U.
     767             : 
     768             :    On input, the content of U is irrelevant. Alternatively, it may contain
     769             :    some nonzero columns that will be used by BVTensorBuildFirstColumn().
     770             : 
     771             :    Level: advanced
     772             : 
     773             : .seealso: BVTensorGetDegree(), BVTensorGetFactors(), BVTensorBuildFirstColumn()
     774             : @*/
     775         165 : PetscErrorCode BVCreateTensor(BV U,PetscInt d,BV *V)
     776             : {
     777         165 :   PetscBool      match;
     778         165 :   PetscInt       n,N,m;
     779         165 :   VecType        vtype;
     780         165 :   BV_TENSOR      *ctx;
     781             : 
     782         165 :   PetscFunctionBegin;
     783         165 :   PetscValidHeaderSpecific(U,BV_CLASSID,1);
     784         660 :   PetscValidLogicalCollectiveInt(U,d,2);
     785         165 :   PetscCall(PetscObjectTypeCompare((PetscObject)U,BVTENSOR,&match));
     786         165 :   PetscCheck(!match,PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"U cannot be of type tensor");
     787             : 
     788         165 :   PetscCall(BVCreate(PetscObjectComm((PetscObject)U),V));
     789         165 :   PetscCall(BVGetSizes(U,&n,&N,&m));
     790         165 :   PetscCheck(m>=d,PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_SIZ,"U has %" PetscInt_FMT " columns, it should have at least d=%" PetscInt_FMT,m,d);
     791         165 :   PetscCall(BVSetSizes(*V,d*n,d*N,m-d+1));
     792         165 :   PetscCall(BVGetVecType(U,&vtype));
     793         165 :   PetscCall(BVSetVecType(*V,vtype));
     794         165 :   PetscCall(PetscObjectChangeTypeName((PetscObject)*V,BVTENSOR));
     795         165 :   PetscCall(PetscLogEventBegin(BV_Create,*V,0,0,0));
     796         165 :   PetscCall(BVCreate_Tensor(*V));
     797         165 :   PetscCall(PetscLogEventEnd(BV_Create,*V,0,0,0));
     798             : 
     799         165 :   ctx = (BV_TENSOR*)(*V)->data;
     800         165 :   ctx->U  = U;
     801         165 :   ctx->d  = d;
     802         165 :   ctx->ld = m;
     803         165 :   PetscCall(PetscObjectReference((PetscObject)U));
     804         165 :   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,d*m,m-d+1,NULL,&ctx->S));
     805         165 :   PetscCall(PetscObjectSetName((PetscObject)ctx->S,"S"));
     806             : 
     807             :   /* Copy user-provided attributes of U */
     808         165 :   (*V)->orthog_type  = U->orthog_type;
     809         165 :   (*V)->orthog_ref   = U->orthog_ref;
     810         165 :   (*V)->orthog_eta   = U->orthog_eta;
     811         165 :   (*V)->orthog_block = U->orthog_block;
     812         165 :   (*V)->vmm          = U->vmm;
     813         165 :   (*V)->rrandom      = U->rrandom;
     814         165 :   PetscFunctionReturn(PETSC_SUCCESS);
     815             : }

Generated by: LCOV version 1.14