Actual source code: bvtensor.c
slepc-3.12.2 2020-01-13
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: Tensor BV that is represented in compact form as V = (I otimes U) S
12: */
14: #include <slepc/private/bvimpl.h> /*I "slepcbv.h" I*/
15: #include <slepcblaslapack.h>
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;
28: PetscErrorCode BVMultInPlace_Tensor(BV V,Mat Q,PetscInt s,PetscInt e)
29: {
31: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
32: PetscScalar *pS,*q;
33: PetscInt ldq,lds = ctx->ld*ctx->d;
36: MatGetSize(Q,&ldq,NULL);
37: MatDenseGetArray(ctx->S,&pS);
38: MatDenseGetArray(Q,&q);
39: BVMultInPlace_BLAS_Private(V,lds,V->k-V->l,ldq,s-V->l,e-V->l,pS+(V->nc+V->l)*lds,q+V->l*ldq+V->l,PETSC_FALSE);
40: MatDenseRestoreArray(Q,&q);
41: MatDenseRestoreArray(ctx->S,&pS);
42: return(0);
43: }
45: PetscErrorCode BVMultInPlaceTranspose_Tensor(BV V,Mat Q,PetscInt s,PetscInt e)
46: {
48: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
49: PetscScalar *pS,*q;
50: PetscInt ldq,lds = ctx->ld*ctx->d;
53: MatGetSize(Q,&ldq,NULL);
54: MatDenseGetArray(ctx->S,&pS);
55: MatDenseGetArray(Q,&q);
56: BVMultInPlace_BLAS_Private(V,lds,V->k-V->l,ldq,s-V->l,e-V->l,pS+(V->nc+V->l)*lds,q+V->l*ldq+V->l,PETSC_TRUE);
57: MatDenseRestoreArray(Q,&q);
58: MatDenseRestoreArray(ctx->S,&pS);
59: return(0);
60: }
62: PetscErrorCode BVDot_Tensor(BV X,BV Y,Mat M)
63: {
65: BV_TENSOR *x = (BV_TENSOR*)X->data,*y = (BV_TENSOR*)Y->data;
66: PetscScalar *m,*px,*py;
67: PetscInt ldm,lds = x->ld*x->d;
70: if (x->U!=y->U) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"BVDot() in BVTENSOR requires that both operands have the same U factor");
71: if (lds!=y->ld*y->d) SETERRQ2(PetscObjectComm((PetscObject)X),PETSC_ERR_ARG_SIZ,"Mismatching dimensions ld*d %D %D",lds,y->ld*y->d);
72: MatGetSize(M,&ldm,NULL);
73: MatDenseGetArray(x->S,&px);
74: MatDenseGetArray(y->S,&py);
75: MatDenseGetArray(M,&m);
76: BVDot_BLAS_Private(X,Y->k-Y->l,X->k-X->l,lds,ldm,py+(Y->nc+Y->l)*lds,px+(X->nc+X->l)*lds,m+X->l*ldm+Y->l,PETSC_FALSE);
77: MatDenseRestoreArray(M,&m);
78: MatDenseRestoreArray(x->S,&px);
79: MatDenseRestoreArray(y->S,&py);
80: return(0);
81: }
83: PetscErrorCode BVScale_Tensor(BV bv,PetscInt j,PetscScalar alpha)
84: {
86: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
87: PetscScalar *pS;
88: PetscInt lds = ctx->ld*ctx->d;
91: MatDenseGetArray(ctx->S,&pS);
92: if (j<0) {
93: BVScale_BLAS_Private(bv,(bv->k-bv->l)*lds,pS+(bv->nc+bv->l)*lds,alpha);
94: } else {
95: BVScale_BLAS_Private(bv,lds,pS+(bv->nc+j)*lds,alpha);
96: }
97: MatDenseRestoreArray(ctx->S,&pS);
98: return(0);
99: }
101: PetscErrorCode BVNorm_Tensor(BV bv,PetscInt j,NormType type,PetscReal *val)
102: {
104: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
105: PetscScalar *pS;
106: PetscInt lds = ctx->ld*ctx->d;
109: MatDenseGetArray(ctx->S,&pS);
110: if (j<0) {
111: BVNorm_LAPACK_Private(bv,lds,bv->k-bv->l,pS+(bv->nc+bv->l)*lds,type,val,PETSC_FALSE);
112: } else {
113: BVNorm_LAPACK_Private(bv,lds,1,pS+(bv->nc+j)*lds,type,val,PETSC_FALSE);
114: }
115: MatDenseRestoreArray(ctx->S,&pS);
116: return(0);
117: }
119: PetscErrorCode BVCopyColumn_Tensor(BV V,PetscInt j,PetscInt i)
120: {
122: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
123: PetscScalar *pS;
124: PetscInt lds = ctx->ld*ctx->d;
127: MatDenseGetArray(ctx->S,&pS);
128: PetscArraycpy(pS+(V->nc+i)*lds,pS+(V->nc+j)*lds,lds);
129: MatDenseRestoreArray(ctx->S,&pS);
130: return(0);
131: }
133: static PetscErrorCode BVTensorNormColumn(BV bv,PetscInt j,PetscReal *norm)
134: {
136: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
137: PetscBLASInt one=1,lds_;
138: PetscScalar sone=1.0,szero=0.0,*S,*x,dot;
139: PetscReal alpha=1.0,scale=0.0,aval;
140: PetscInt i,lds,ld=ctx->ld;
143: lds = ld*ctx->d;
144: MatDenseGetArray(ctx->S,&S);
145: PetscBLASIntCast(lds,&lds_);
146: if (ctx->qB) {
147: x = ctx->sw;
148: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,S+j*lds,&one,&szero,x,&one));
149: dot = PetscRealPart(BLASdot_(&lds_,S+j*lds,&one,x,&one));
150: BV_SafeSqrt(bv,dot,norm);
151: } else {
152: /* Compute *norm = BLASnrm2_(&lds_,S+j*lds,&one); */
153: if (lds==1) *norm = PetscAbsScalar(S[j*lds]);
154: else {
155: for (i=0;i<lds;i++) {
156: aval = PetscAbsScalar(S[i+j*lds]);
157: if (aval!=0.0) {
158: if (scale<aval) {
159: alpha = 1.0 + alpha*PetscSqr(scale/aval);
160: scale = aval;
161: } else alpha += PetscSqr(aval/scale);
162: }
163: }
164: *norm = scale*PetscSqrtReal(alpha);
165: }
166: }
167: return(0);
168: }
170: PetscErrorCode BVOrthogonalizeGS1_Tensor(BV bv,PetscInt k,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onorm,PetscReal *norm)
171: {
172: PetscErrorCode ierr;
173: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
174: PetscScalar *pS,*cc,*x,dot,sonem=-1.0,sone=1.0,szero=0.0;
175: PetscInt i,lds = ctx->ld*ctx->d;
176: PetscBLASInt lds_,k_,one=1;
177: const PetscScalar *omega;
180: if (v) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Orthogonalization against an external vector is not allowed in BVTENSOR");
181: MatDenseGetArray(ctx->S,&pS);
182: if (!c) {
183: VecGetArray(bv->buffer,&cc);
184: } else cc = c;
185: PetscBLASIntCast(lds,&lds_);
186: PetscBLASIntCast(k,&k_);
188: if (onorm) { BVTensorNormColumn(bv,k,onorm); }
190: if (ctx->qB) x = ctx->sw;
191: else x = pS+k*lds;
193: if (bv->orthog_type==BV_ORTHOG_MGS) { /* modified Gram-Schmidt */
195: if (bv->indef) { /* signature */
196: VecGetArrayRead(bv->omega,&omega);
197: }
198: for (i=-bv->nc;i<k;i++) {
199: if (which && i>=0 && !which[i]) continue;
200: if (ctx->qB) PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,pS+k*lds,&one,&szero,x,&one));
201: /* c_i = (s_k, s_i) */
202: dot = PetscRealPart(BLASdot_(&lds_,pS+i*lds,&one,x,&one));
203: if (bv->indef) dot /= PetscRealPart(omega[i]);
204: BV_SetValue(bv,i,0,cc,dot);
205: /* s_k = s_k - c_i s_i */
206: dot = -dot;
207: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&lds_,&dot,pS+i*lds,&one,pS+k*lds,&one));
208: }
209: if (bv->indef) {
210: VecRestoreArrayRead(bv->omega,&omega);
211: }
213: } else { /* classical Gram-Schmidt */
214: if (ctx->qB) PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,pS+k*lds,&one,&szero,x,&one));
216: /* cc = S_{0:k-1}^* s_k */
217: PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&lds_,&k_,&sone,pS,&lds_,x,&one,&szero,cc,&one));
219: /* s_k = s_k - S_{0:k-1} cc */
220: if (bv->indef) { BV_ApplySignature(bv,k,cc,PETSC_TRUE); }
221: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&k_,&sonem,pS,&lds_,cc,&one,&sone,pS+k*lds,&one));
222: if (bv->indef) { BV_ApplySignature(bv,k,cc,PETSC_FALSE); }
223: }
225: if (norm) { BVTensorNormColumn(bv,k,norm); }
226: BV_AddCoefficients(bv,k,h,cc);
227: MatDenseRestoreArray(ctx->S,&pS);
228: VecRestoreArray(bv->buffer,&cc);
229: return(0);
230: }
232: PetscErrorCode BVView_Tensor(BV bv,PetscViewer viewer)
233: {
234: PetscErrorCode ierr;
235: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
236: PetscViewerFormat format;
237: PetscBool isascii;
238: const char *bvname,*uname,*sname;
241: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
242: if (isascii) {
243: PetscViewerGetFormat(viewer,&format);
244: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
245: PetscViewerASCIIPrintf(viewer,"number of tensor blocks (degree): %D\n",ctx->d);
246: PetscViewerASCIIPrintf(viewer,"number of columns of U factor: %D\n",ctx->ld);
247: return(0);
248: }
249: BVView(ctx->U,viewer);
250: MatView(ctx->S,viewer);
251: if (format == PETSC_VIEWER_ASCII_MATLAB) {
252: PetscObjectGetName((PetscObject)bv,&bvname);
253: PetscObjectGetName((PetscObject)ctx->U,&uname);
254: PetscObjectGetName((PetscObject)ctx->S,&sname);
255: PetscViewerASCIIPrintf(viewer,"%s=kron(eye(%D),%s)*%s(:,1:%D);\n",bvname,ctx->d,uname,sname,bv->k);
256: }
257: } else {
258: BVView(ctx->U,viewer);
259: MatView(ctx->S,viewer);
260: }
261: return(0);
262: }
264: static PetscErrorCode BVTensorUpdateMatrix(BV V,PetscInt ini,PetscInt end)
265: {
267: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
268: PetscInt i,j,r,c,l,k,ld=ctx->ld,lds=ctx->d*ctx->ld;
269: PetscScalar *qB,*sqB;
270: Vec u;
271: Mat A;
274: if (!V->matrix) return(0);
275: l = ctx->U->l; k = ctx->U->k;
276: /* update inner product matrix */
277: if (!ctx->qB) {
278: PetscCalloc2(lds*lds,&ctx->qB,lds,&ctx->sw);
279: VecDuplicate(ctx->U->t,&ctx->u);
280: }
281: ctx->U->l = 0;
282: for (r=0;r<ctx->d;r++) {
283: for (c=0;c<=r;c++) {
284: MatNestGetSubMat(V->matrix,r,c,&A);
285: if (A) {
286: qB = ctx->qB+c*ld*lds+r*ld;
287: for (i=ini;i<end;i++) {
288: BVGetColumn(ctx->U,i,&u);
289: MatMult(A,u,ctx->u);
290: ctx->U->k = i+1;
291: BVDotVec(ctx->U,ctx->u,qB+i*lds);
292: BVRestoreColumn(ctx->U,i,&u);
293: for (j=0;j<i;j++) qB[i+j*lds] = PetscConj(qB[j+i*lds]);
294: qB[i*lds+i] = PetscRealPart(qB[i+i*lds]);
295: }
296: if (c!=r) {
297: sqB = ctx->qB+r*ld*lds+c*ld;
298: for (i=ini;i<end;i++) for (j=0;j<=i;j++) {
299: sqB[i+j*lds] = PetscConj(qB[j+i*lds]);
300: sqB[j+i*lds] = qB[j+i*lds];
301: }
302: }
303: }
304: }
305: }
306: ctx->U->l = l; ctx->U->k = k;
307: return(0);
308: }
310: static PetscErrorCode BVTensorBuildFirstColumn_Tensor(BV V,PetscInt k)
311: {
313: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
314: PetscInt i,nq=0;
315: PetscScalar *pS,*omega;
316: PetscReal norm;
317: PetscBool breakdown=PETSC_FALSE;
320: MatDenseGetArray(ctx->S,&pS);
321: for (i=0;i<ctx->d;i++) {
322: if (i>=k) {
323: BVSetRandomColumn(ctx->U,nq);
324: } else {
325: BVCopyColumn(ctx->U,i,nq);
326: }
327: BVOrthogonalizeColumn(ctx->U,nq,pS+i*ctx->ld,&norm,&breakdown);
328: if (!breakdown) {
329: BVScaleColumn(ctx->U,nq,1.0/norm);
330: pS[nq+i*ctx->ld] = norm;
331: nq++;
332: }
333: }
334: MatDenseRestoreArray(ctx->S,&pS);
335: if (!nq) SETERRQ1(PetscObjectComm((PetscObject)V),1,"Cannot build first column of tensor BV; U should contain k=%D nonzero columns",k);
336: BVTensorUpdateMatrix(V,0,nq);
337: BVTensorNormColumn(V,0,&norm);
338: BVScale_Tensor(V,0,1.0/norm);
339: if (V->indef) {
340: BV_AllocateSignature(V);
341: VecGetArray(V->omega,&omega);
342: omega[0] = (norm<0.0)? -1.0: 1.0;
343: VecRestoreArray(V->omega,&omega);
344: }
345: /* set active columns */
346: ctx->U->l = 0;
347: ctx->U->k = nq;
348: return(0);
349: }
351: /*@
352: BVTensorBuildFirstColumn - Builds the first column of the tensor basis vectors
353: V from the data contained in the first k columns of U.
355: Collective on V
357: Input Parameters:
358: + V - the basis vectors context
359: - k - the number of columns of U with relevant information
361: Notes:
362: At most d columns are considered, where d is the degree of the tensor BV.
363: Given V = (I otimes U) S, this function computes the first column of V, that
364: is, it computes the coefficients of the first column of S by orthogonalizing
365: the first d columns of U. If k is less than d (or linearly dependent columns
366: are found) then additional random columns are used.
368: The computed column has unit norm.
370: Level: advanced
372: .seealso: BVCreateTensor()
373: @*/
374: PetscErrorCode BVTensorBuildFirstColumn(BV V,PetscInt k)
375: {
381: PetscUseMethod(V,"BVTensorBuildFirstColumn_C",(BV,PetscInt),(V,k));
382: return(0);
383: }
385: static PetscErrorCode BVTensorCompress_Tensor(BV V,PetscInt newc)
386: {
387: #if defined(PETSC_MISSING_LAPACK_GESVD) || defined(PETSC_MISSING_LAPACK_GEQRF) || defined(PETSC_MISSING_LAPACK_ORGQR)
389: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESVD/GEQRF/ORGQR - Lapack routine is unavailable");
390: #else
392: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
393: PetscInt nwu=0,nnc,nrow,lwa,r,c;
394: PetscInt i,j,k,n,lds=ctx->ld*ctx->d,deg=ctx->d,lock,cs1=V->k,rs1=ctx->U->k,rk=0,offu;
395: PetscScalar *S,*M,*Z,*pQ,*SS,*SS2,t,sone=1.0,zero=0.0,mone=-1.0,*p,*tau,*work,*qB,*sqB;
396: PetscReal *sg,tol,*rwork;
397: PetscBLASInt ld_,cs1_,rs1_,cs1tdeg,n_,info,lw_,newc_,newctdeg,nnc_,nrow_,nnctdeg,lds_,rk_;
398: Mat Q,A;
401: if (!cs1) return(0);
402: lwa = 6*ctx->ld*lds+2*cs1;
403: n = PetscMin(rs1,deg*cs1);
404: lock = ctx->U->l;
405: nnc = cs1-lock-newc;
406: nrow = rs1-lock;
407: PetscCalloc6(deg*newc*nnc,&SS,newc*nnc,&SS2,(rs1+lock+newc)*n,&pQ,deg*rs1,&tau,lwa,&work,6*n,&rwork);
408: offu = lock*(rs1+1);
409: M = work+nwu;
410: nwu += rs1*cs1*deg;
411: sg = rwork;
412: Z = work+nwu;
413: nwu += deg*cs1*n;
414: PetscBLASIntCast(n,&n_);
415: PetscBLASIntCast(nnc,&nnc_);
416: PetscBLASIntCast(cs1,&cs1_);
417: PetscBLASIntCast(rs1,&rs1_);
418: PetscBLASIntCast(newc,&newc_);
419: PetscBLASIntCast(newc*deg,&newctdeg);
420: PetscBLASIntCast(nnc*deg,&nnctdeg);
421: PetscBLASIntCast(cs1*deg,&cs1tdeg);
422: PetscBLASIntCast(lwa-nwu,&lw_);
423: PetscBLASIntCast(nrow,&nrow_);
424: PetscBLASIntCast(lds,&lds_);
425: MatDenseGetArray(ctx->S,&S);
427: if (newc>0) {
428: /* truncate columns associated with new converged eigenpairs */
429: for (j=0;j<deg;j++) {
430: for (i=lock;i<lock+newc;i++) {
431: PetscArraycpy(M+(i-lock+j*newc)*nrow,S+i*lds+j*ctx->ld+lock,nrow);
432: }
433: }
434: #if !defined (PETSC_USE_COMPLEX)
435: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&newctdeg,M,&nrow_,sg,pQ+offu,&rs1_,Z,&n_,work+nwu,&lw_,&info));
436: #else
437: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&newctdeg,M,&nrow_,sg,pQ+offu,&rs1_,Z,&n_,work+nwu,&lw_,rwork+n,&info));
438: #endif
439: SlepcCheckLapackInfo("gesvd",info);
440: /* SVD has rank min(newc,nrow) */
441: rk = PetscMin(newc,nrow);
442: for (i=0;i<rk;i++) {
443: t = sg[i];
444: PetscStackCallBLAS("BLASscal",BLASscal_(&newctdeg,&t,Z+i,&n_));
445: }
446: for (i=0;i<deg;i++) {
447: for (j=lock;j<lock+newc;j++) {
448: PetscArraycpy(S+j*lds+i*ctx->ld+lock,Z+(newc*i+j-lock)*n,rk);
449: PetscArrayzero(S+j*lds+i*ctx->ld+lock+rk,(ctx->ld-lock-rk));
450: }
451: }
452: /*
453: update columns associated with non-converged vectors, orthogonalize
454: against pQ so that next M has rank nnc+d-1 insted of nrow+d-1
455: */
456: for (i=0;i<deg;i++) {
457: PetscStackCallBLAS("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_));
458: PetscStackCallBLAS("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_));
459: /* repeat orthogonalization step */
460: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&newc_,&nnc_,&nrow_,&sone,pQ+offu,&rs1_,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_,&zero,SS2,&newc_));
461: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&nrow_,&nnc_,&newc_,&mone,pQ+offu,&rs1_,SS2,&newc_,&sone,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_));
462: for (j=0;j<newc*nnc;j++) *(SS+i*newc*nnc+j) += SS2[j];
463: }
464: }
466: /* truncate columns associated with non-converged eigenpairs */
467: for (j=0;j<deg;j++) {
468: for (i=lock+newc;i<cs1;i++) {
469: PetscArraycpy(M+(i-lock-newc+j*nnc)*nrow,S+i*lds+j*ctx->ld+lock,nrow);
470: }
471: }
472: #if !defined (PETSC_USE_COMPLEX)
473: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&nnctdeg,M,&nrow_,sg,pQ+offu+newc*rs1,&rs1_,Z,&n_,work+nwu,&lw_,&info));
474: #else
475: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&nnctdeg,M,&nrow_,sg,pQ+offu+newc*rs1,&rs1_,Z,&n_,work+nwu,&lw_,rwork+n,&info));
476: #endif
477: SlepcCheckLapackInfo("gesvd",info);
478: tol = PetscMax(rs1,deg*cs1)*PETSC_MACHINE_EPSILON*sg[0];
479: rk = 0;
480: for (i=0;i<PetscMin(nrow,nnctdeg);i++) if (sg[i]>tol) rk++;
481: rk = PetscMin(nnc+deg-1,rk);
482: /* the SVD has rank (at most) nnc+deg-1 */
483: for (i=0;i<rk;i++) {
484: t = sg[i];
485: PetscStackCallBLAS("BLASscal",BLASscal_(&nnctdeg,&t,Z+i,&n_));
486: }
487: /* update S */
488: PetscArrayzero(S+cs1*lds,(V->m-cs1)*lds);
489: k = ctx->ld-lock-newc-rk;
490: for (i=0;i<deg;i++) {
491: for (j=lock+newc;j<cs1;j++) {
492: PetscArraycpy(S+j*lds+i*ctx->ld+lock+newc,Z+(nnc*i+j-lock-newc)*n,rk);
493: PetscArrayzero(S+j*lds+i*ctx->ld+lock+newc+rk,k);
494: }
495: }
496: if (newc>0) {
497: for (i=0;i<deg;i++) {
498: p = SS+nnc*newc*i;
499: for (j=lock+newc;j<cs1;j++) {
500: for (k=0;k<newc;k++) S[j*lds+i*ctx->ld+lock+k] = *(p++);
501: }
502: }
503: }
505: /* orthogonalize pQ */
506: rk = rk+newc;
507: PetscBLASIntCast(rk,&rk_);
508: PetscBLASIntCast(cs1-lock,&nnc_);
509: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&nrow_,&rk_,pQ+offu,&rs1_,tau,work+nwu,&lw_,&info));
510: SlepcCheckLapackInfo("geqrf",info);
511: for (i=0;i<deg;i++) {
512: PetscStackCallBLAS("BLAStrmm",BLAStrmm_("L","U","N","N",&rk_,&nnc_,&sone,pQ+offu,&rs1_,S+lock*lds+lock+i*ctx->ld,&lds_));
513: }
514: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&nrow_,&rk_,&rk_,pQ+offu,&rs1_,tau,work+nwu,&lw_,&info));
515: SlepcCheckLapackInfo("orgqr",info);
517: /* update vectors U(:,idx) = U*Q(:,idx) */
518: rk = rk+lock;
519: for (i=0;i<lock;i++) pQ[i*(1+rs1)] = 1.0;
520: MatCreateSeqDense(PETSC_COMM_SELF,rs1,rk,pQ,&Q);
521: ctx->U->k = rs1;
522: BVMultInPlace(ctx->U,Q,lock,rk);
523: MatDestroy(&Q);
525: if (ctx->qB) {
526: /* update matrix qB */
527: PetscBLASIntCast(ctx->ld,&ld_);
528: PetscBLASIntCast(rk,&rk_);
529: for (r=0;r<ctx->d;r++) {
530: for (c=0;c<=r;c++) {
531: MatNestGetSubMat(V->matrix,r,c,&A);
532: if (A) {
533: qB = ctx->qB+r*ctx->ld+c*ctx->ld*lds;
534: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&rs1_,&rk_,&rs1_,&sone,qB,&lds_,pQ,&rs1_,&zero,work+nwu,&rs1_));
535: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&rk_,&rk_,&rs1_,&sone,pQ,&rs1_,work+nwu,&rs1_,&zero,qB,&lds_));
536: for (i=0;i<rk;i++) {
537: for (j=0;j<i;j++) qB[i+j*lds] = PetscConj(qB[j+i*lds]);
538: qB[i+i*lds] = PetscRealPart(qB[i+i*lds]);
539: }
540: for (i=rk;i<ctx->ld;i++) {
541: PetscArrayzero(qB+i*lds,ctx->ld);
542: }
543: for (i=0;i<rk;i++) {
544: PetscArrayzero(qB+i*lds+rk,(ctx->ld-rk));
545: }
546: if (c!=r) {
547: sqB = ctx->qB+r*ctx->ld*lds+c*ctx->ld;
548: for (i=0;i<ctx->ld;i++) for (j=0;j<ctx->ld;j++) sqB[i+j*lds] = PetscConj(qB[j+i*lds]);
549: }
550: }
551: }
552: }
553: }
555: /* free work space */
556: PetscFree6(SS,SS2,pQ,tau,work,rwork);
557: MatDenseRestoreArray(ctx->S,&S);
559: /* set active columns */
560: if (newc) ctx->U->l += newc;
561: ctx->U->k = rk;
562: return(0);
563: #endif
564: }
566: /*@
567: BVTensorCompress - Updates the U and S factors of the tensor basis vectors
568: object V by means of an SVD, removing redundant information.
570: Collective on V
572: Input Parameters:
573: + V - the tensor basis vectors context
574: - newc - additional columns to be locked
576: Notes:
577: This function is typically used when restarting Krylov solvers. Truncating a
578: tensor BV V = (I otimes U) S to its leading columns amounts to keeping the
579: leading columns of S. However, to effectively reduce the size of the
580: decomposition, it is necessary to compress it in a way that fewer columns of
581: U are employed. This can be achieved by means of an update that involves the
582: SVD of the low-rank matrix [S_0 S_1 ... S_{d-1}], where S_i are the pieces of S.
584: If newc is nonzero, then newc columns are added to the leading columns of V.
585: This means that the corresponding columns of the U and S factors will remain
586: invariant in subsequent operations.
588: Level: advanced
590: .seealso: BVCreateTensor(), BVSetActiveColumns()
591: @*/
592: PetscErrorCode BVTensorCompress(BV V,PetscInt newc)
593: {
599: PetscUseMethod(V,"BVTensorCompress_C",(BV,PetscInt),(V,newc));
600: return(0);
601: }
603: static PetscErrorCode BVTensorGetDegree_Tensor(BV bv,PetscInt *d)
604: {
605: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
608: *d = ctx->d;
609: return(0);
610: }
612: /*@
613: BVTensorGetDegree - Returns the number of blocks (degree) of the tensor BV.
615: Not collective
617: Input Parameter:
618: . bv - the basis vectors context
620: Output Parameter:
621: . d - the degree
623: Level: advanced
625: .seealso: BVCreateTensor()
626: @*/
627: PetscErrorCode BVTensorGetDegree(BV bv,PetscInt *d)
628: {
634: PetscUseMethod(bv,"BVTensorGetDegree_C",(BV,PetscInt*),(bv,d));
635: return(0);
636: }
638: static PetscErrorCode BVTensorGetFactors_Tensor(BV V,BV *U,Mat *S)
639: {
640: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
643: if (ctx->puk>-1) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ORDER,"Previous call to BVTensonGetFactors without a BVTensorRestoreFactors call");
644: ctx->puk = ctx->U->k;
645: if (U) *U = ctx->U;
646: if (S) *S = ctx->S;
647: return(0);
648: }
650: /*@C
651: BVTensorGetFactors - Returns the two factors involved in the definition of the
652: tensor basis vectors object, V = (I otimes U) S.
654: Logically Collective on V
656: Input Parameter:
657: . V - the basis vectors context
659: Output Parameters:
660: + U - the BV factor
661: - S - the Mat factor
663: Notes:
664: The returned factors are references (not copies) of the internal factors,
665: so modifying them will change the tensor BV as well. Some operations of the
666: tensor BV assume that U has orthonormal columns, so if the user modifies U
667: this restriction must be taken into account.
669: The returned factors must not be destroyed. BVTensorRestoreFactors() must
670: be called when they are no longer needed.
672: Pass a NULL vector for any of the arguments that is not needed.
674: Level: advanced
676: .seealso: BVTensorRestoreFactors()
677: @*/
678: PetscErrorCode BVTensorGetFactors(BV V,BV *U,Mat *S)
679: {
684: PetscUseMethod(V,"BVTensorGetFactors_C",(BV,BV*,Mat*),(V,U,S));
685: return(0);
686: }
688: static PetscErrorCode BVTensorRestoreFactors_Tensor(BV V,BV *U,Mat *S)
689: {
691: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
694: PetscObjectStateIncrease((PetscObject)V);
695: if (U) *U = NULL;
696: if (S) *S = NULL;
697: BVTensorUpdateMatrix(V,ctx->puk,ctx->U->k);
698: ctx->puk = -1;
699: return(0);
700: }
702: /*@C
703: BVTensorRestoreFactors - Restore the two factors that were obtained with
704: BVTensorGetFactors().
706: Logically Collective on V
708: Input Parameters:
709: + V - the basis vectors context
710: . U - the BV factor (or NULL)
711: - S - the Mat factor (or NULL)
713: Nots:
714: The arguments must match the corresponding call to BVTensorGetFactors().
716: Level: advanced
718: .seealso: BVTensorGetFactors()
719: @*/
720: PetscErrorCode BVTensorRestoreFactors(BV V,BV *U,Mat *S)
721: {
728: PetscUseMethod(V,"BVTensorRestoreFactors_C",(BV,BV*,Mat*),(V,U,S));
729: return(0);
730: }
732: PetscErrorCode BVDestroy_Tensor(BV bv)
733: {
735: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
738: BVDestroy(&ctx->U);
739: MatDestroy(&ctx->S);
740: if (ctx->u) {
741: PetscFree2(ctx->qB,ctx->sw);
742: VecDestroy(&ctx->u);
743: }
744: PetscFree(bv->data);
745: PetscObjectComposeFunction((PetscObject)bv,"BVTensorBuildFirstColumn_C",NULL);
746: PetscObjectComposeFunction((PetscObject)bv,"BVTensorCompress_C",NULL);
747: PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetDegree_C",NULL);
748: PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetFactors_C",NULL);
749: PetscObjectComposeFunction((PetscObject)bv,"BVTensorRestoreFactors_C",NULL);
750: return(0);
751: }
753: SLEPC_EXTERN PetscErrorCode BVCreate_Tensor(BV bv)
754: {
756: BV_TENSOR *ctx;
759: PetscNewLog(bv,&ctx);
760: bv->data = (void*)ctx;
761: ctx->puk = -1;
763: bv->ops->multinplace = BVMultInPlace_Tensor;
764: bv->ops->multinplacetrans = BVMultInPlaceTranspose_Tensor;
765: bv->ops->dot = BVDot_Tensor;
766: bv->ops->scale = BVScale_Tensor;
767: bv->ops->norm = BVNorm_Tensor;
768: bv->ops->copycolumn = BVCopyColumn_Tensor;
769: bv->ops->gramschmidt = BVOrthogonalizeGS1_Tensor;
770: bv->ops->destroy = BVDestroy_Tensor;
771: bv->ops->view = BVView_Tensor;
773: PetscObjectComposeFunction((PetscObject)bv,"BVTensorBuildFirstColumn_C",BVTensorBuildFirstColumn_Tensor);
774: PetscObjectComposeFunction((PetscObject)bv,"BVTensorCompress_C",BVTensorCompress_Tensor);
775: PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetDegree_C",BVTensorGetDegree_Tensor);
776: PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetFactors_C",BVTensorGetFactors_Tensor);
777: PetscObjectComposeFunction((PetscObject)bv,"BVTensorRestoreFactors_C",BVTensorRestoreFactors_Tensor);
778: return(0);
779: }
781: /*@
782: BVCreateTensor - Creates a tensor BV that is represented in compact form
783: as V = (I otimes U) S, where U has orthonormal columns.
785: Collective on U
787: Input Parameters:
788: + U - a basis vectors object
789: - d - the number of blocks (degree) of the tensor BV
791: Output Parameter:
792: . V - the new basis vectors context
794: Notes:
795: The new basis vectors object is V = (I otimes U) S, where otimes denotes
796: the Kronecker product, I is the identity matrix of order d, and S is a
797: sequential matrix allocated internally. This compact representation is
798: used e.g. to represent the Krylov basis generated with the linearization
799: of a matrix polynomial of degree d.
801: The size of V (number of rows) is equal to d times n, where n is the size
802: of U. The dimensions of S are d times m rows and m-d+1 columns, where m is
803: the number of columns of U, so m should be at least d.
805: The communicator of V will be the same as U.
807: On input, the content of U is irrelevant. Alternatively, it may contain
808: some nonzero columns that will be used by BVTensorBuildFirstColumn().
810: Level: advanced
812: .seealso: BVTensorGetDegree(), BVTensorGetFactors(), BVTensorBuildFirstColumn()
813: @*/
814: PetscErrorCode BVCreateTensor(BV U,PetscInt d,BV *V)
815: {
817: PetscBool match;
818: PetscInt n,N,m;
819: BV_TENSOR *ctx;
824: PetscObjectTypeCompare((PetscObject)U,BVTENSOR,&match);
825: if (match) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"U cannot be of type tensor");
827: BVCreate(PetscObjectComm((PetscObject)U),V);
828: BVGetSizes(U,&n,&N,&m);
829: if (m<d) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_SIZ,"U has %D columns, it should have at least d=%D",m,d);
830: BVSetSizes(*V,d*n,d*N,m-d+1);
831: PetscObjectChangeTypeName((PetscObject)*V,BVTENSOR);
832: PetscLogEventBegin(BV_Create,*V,0,0,0);
833: BVCreate_Tensor(*V);
834: PetscLogEventEnd(BV_Create,*V,0,0,0);
836: ctx = (BV_TENSOR*)(*V)->data;
837: ctx->U = U;
838: ctx->d = d;
839: ctx->ld = m;
840: PetscObjectReference((PetscObject)U);
841: PetscLogObjectParent((PetscObject)*V,(PetscObject)U);
842: MatCreateSeqDense(PETSC_COMM_SELF,d*m,m-d+1,NULL,&ctx->S);
843: PetscLogObjectParent((PetscObject)*V,(PetscObject)ctx->S);
844: PetscObjectSetName((PetscObject)ctx->S,"S");
846: /* Copy user-provided attributes of U */
847: (*V)->orthog_type = U->orthog_type;
848: (*V)->orthog_ref = U->orthog_ref;
849: (*V)->orthog_eta = U->orthog_eta;
850: (*V)->orthog_block = U->orthog_block;
851: (*V)->vmm = U->vmm;
852: (*V)->rrandom = U->rrandom;
853: return(0);
854: }