Actual source code: bvtensor.c
slepc-3.11.2 2019-07-30
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 BVMult_Tensor(BV Y,PetscScalar alpha,PetscScalar beta,BV X,Mat Q)
29: {
31: SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
32: }
34: PetscErrorCode BVMultVec_Tensor(BV X,PetscScalar alpha,PetscScalar beta,Vec y,PetscScalar *q)
35: {
37: SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
38: }
40: PetscErrorCode BVMultInPlace_Tensor(BV V,Mat Q,PetscInt s,PetscInt e)
41: {
43: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
44: PetscScalar *pS,*q;
45: PetscInt ldq,lds = ctx->ld*ctx->d;
48: MatGetSize(Q,&ldq,NULL);
49: MatDenseGetArray(ctx->S,&pS);
50: MatDenseGetArray(Q,&q);
51: 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);
52: MatDenseRestoreArray(Q,&q);
53: MatDenseRestoreArray(ctx->S,&pS);
54: return(0);
55: }
57: PetscErrorCode BVMultInPlaceTranspose_Tensor(BV V,Mat Q,PetscInt s,PetscInt e)
58: {
60: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
61: PetscScalar *pS,*q;
62: PetscInt ldq,lds = ctx->ld*ctx->d;
65: MatGetSize(Q,&ldq,NULL);
66: MatDenseGetArray(ctx->S,&pS);
67: MatDenseGetArray(Q,&q);
68: 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);
69: MatDenseRestoreArray(Q,&q);
70: MatDenseRestoreArray(ctx->S,&pS);
71: return(0);
72: }
74: PetscErrorCode BVDot_Tensor(BV X,BV Y,Mat M)
75: {
77: BV_TENSOR *x = (BV_TENSOR*)X->data,*y = (BV_TENSOR*)Y->data;
78: PetscScalar *m,*px,*py;
79: PetscInt ldm,lds = x->ld*x->d;
82: if (x->U!=y->U) SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"BVDot() in BVTENSOR requires that both operands have the same U factor");
83: 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);
84: MatGetSize(M,&ldm,NULL);
85: MatDenseGetArray(x->S,&px);
86: MatDenseGetArray(y->S,&py);
87: MatDenseGetArray(M,&m);
88: 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);
89: MatDenseRestoreArray(M,&m);
90: MatDenseRestoreArray(x->S,&px);
91: MatDenseRestoreArray(y->S,&py);
92: return(0);
93: }
95: PetscErrorCode BVDotVec_Tensor(BV X,Vec y,PetscScalar *q)
96: {
98: SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
99: }
101: PetscErrorCode BVDotVec_Local_Tensor(BV X,Vec y,PetscScalar *m)
102: {
104: SETERRQ(PetscObjectComm((PetscObject)X),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
105: }
107: PetscErrorCode BVScale_Tensor(BV bv,PetscInt j,PetscScalar alpha)
108: {
110: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
111: PetscScalar *pS;
112: PetscInt lds = ctx->ld*ctx->d;
115: MatDenseGetArray(ctx->S,&pS);
116: if (j<0) {
117: BVScale_BLAS_Private(bv,(bv->k-bv->l)*lds,pS+(bv->nc+bv->l)*lds,alpha);
118: } else {
119: BVScale_BLAS_Private(bv,lds,pS+(bv->nc+j)*lds,alpha);
120: }
121: MatDenseRestoreArray(ctx->S,&pS);
122: return(0);
123: }
125: PetscErrorCode BVNorm_Tensor(BV bv,PetscInt j,NormType type,PetscReal *val)
126: {
128: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
129: PetscScalar *pS;
130: PetscInt lds = ctx->ld*ctx->d;
133: MatDenseGetArray(ctx->S,&pS);
134: if (j<0) {
135: BVNorm_LAPACK_Private(bv,lds,bv->k-bv->l,pS+(bv->nc+bv->l)*lds,type,val,PETSC_FALSE);
136: } else {
137: BVNorm_LAPACK_Private(bv,lds,1,pS+(bv->nc+j)*lds,type,val,PETSC_FALSE);
138: }
139: MatDenseRestoreArray(ctx->S,&pS);
140: return(0);
141: }
143: PetscErrorCode BVNorm_Local_Tensor(BV bv,PetscInt j,NormType type,PetscReal *val)
144: {
146: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
147: }
149: PetscErrorCode BVMatMult_Tensor(BV V,Mat A,BV W)
150: {
152: SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
153: }
155: PetscErrorCode BVCopy_Tensor(BV V,BV W)
156: {
158: SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
159: }
161: PetscErrorCode BVCopyColumn_Tensor(BV V,PetscInt j,PetscInt i)
162: {
164: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
165: PetscScalar *pS;
166: PetscInt lds = ctx->ld*ctx->d;
169: MatDenseGetArray(ctx->S,&pS);
170: PetscMemcpy(pS+(V->nc+i)*lds,pS+(V->nc+j)*lds,lds*sizeof(PetscScalar));
171: MatDenseRestoreArray(ctx->S,&pS);
172: return(0);
173: }
175: PetscErrorCode BVResize_Tensor(BV bv,PetscInt m,PetscBool copy)
176: {
178: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
179: }
181: PetscErrorCode BVGetColumn_Tensor(BV bv,PetscInt j,Vec *v)
182: {
184: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
185: }
187: PetscErrorCode BVRestoreColumn_Tensor(BV bv,PetscInt j,Vec *v)
188: {
190: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
191: }
193: PetscErrorCode BVGetArray_Tensor(BV bv,PetscScalar **a)
194: {
196: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
197: }
199: PetscErrorCode BVRestoreArray_Tensor(BV bv,PetscScalar **a)
200: {
202: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
203: }
205: PetscErrorCode BVGetArrayRead_Tensor(BV bv,const PetscScalar **a)
206: {
208: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
209: }
211: PetscErrorCode BVRestoreArrayRead_Tensor(BV bv,const PetscScalar **a)
212: {
214: SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Operation not implemented in BVTENSOR");
215: }
217: static PetscErrorCode BVTensorNormColumn(BV bv,PetscInt j,PetscReal *norm)
218: {
220: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
221: PetscBLASInt one=1,lds_;
222: PetscScalar sone=1.0,szero=0.0,*S,*x,dot;
223: PetscReal alpha=1.0,scale=0.0,aval;
224: PetscInt i,lds,ld=ctx->ld;
227: lds = ld*ctx->d;
228: MatDenseGetArray(ctx->S,&S);
229: PetscBLASIntCast(lds,&lds_);
230: if (ctx->qB) {
231: x = ctx->sw;
232: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,S+j*lds,&one,&szero,x,&one));
233: dot = BLASdot_(&lds_,S+j*lds,&one,x,&one);
234: BV_SafeSqrt(bv,dot,norm);
235: } else {
236: /* Compute *norm = BLASnrm2_(&lds_,S+j*lds,&one); */
237: if (lds==1) *norm = PetscAbsScalar(S[j*lds]);
238: else {
239: for (i=0;i<lds;i++) {
240: aval = PetscAbsScalar(S[i+j*lds]);
241: if (aval!=0.0) {
242: if (scale<aval) {
243: alpha = 1.0 + alpha*PetscSqr(scale/aval);
244: scale = aval;
245: } else alpha += PetscSqr(aval/scale);
246: }
247: }
248: *norm = scale*PetscSqrtReal(alpha);
249: }
250: }
251: return(0);
252: }
254: PetscErrorCode BVOrthogonalizeGS1_Tensor(BV bv,PetscInt k,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onorm,PetscReal *norm)
255: {
256: PetscErrorCode ierr;
257: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
258: PetscScalar *pS,*cc,*x,dot,sonem=-1.0,sone=1.0,szero=0.0;
259: PetscInt i,lds = ctx->ld*ctx->d;
260: PetscBLASInt lds_,k_,one=1;
261: const PetscScalar *omega;
264: if (v) SETERRQ(PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Orthogonalization against an external vector is not allowed in BVTENSOR");
265: MatDenseGetArray(ctx->S,&pS);
266: if (!c) {
267: VecGetArray(bv->buffer,&cc);
268: } else cc = c;
269: PetscBLASIntCast(lds,&lds_);
270: PetscBLASIntCast(k,&k_);
272: if (onorm) { BVTensorNormColumn(bv,k,onorm); }
274: if (ctx->qB) x = ctx->sw;
275: else x = pS+k*lds;
277: if (bv->orthog_type==BV_ORTHOG_MGS) { /* modified Gram-Schmidt */
279: if (bv->indef) { /* signature */
280: VecGetArrayRead(bv->omega,&omega);
281: }
282: for (i=-bv->nc;i<k;i++) {
283: if (which && i>=0 && !which[i]) continue;
284: if (ctx->qB) PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,pS+k*lds,&one,&szero,x,&one));
285: /* c_i = ( s_k, s_i ) */
286: dot = BLASdot_(&lds_,pS+i*lds,&one,x,&one);
287: if (bv->indef) dot /= PetscRealPart(omega[i]);
288: BV_SetValue(bv,i,0,cc,dot);
289: /* s_k = s_k - c_i s_i */
290: dot = -dot;
291: PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&lds_,&dot,pS+i*lds,&one,pS+k*lds,&one));
292: }
293: if (bv->indef) {
294: VecRestoreArrayRead(bv->omega,&omega);
295: }
297: } else { /* classical Gram-Schmidt */
298: if (ctx->qB) PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,pS+k*lds,&one,&szero,x,&one));
300: /* cc = S_{0:k-1}^* s_k */
301: PetscStackCallBLAS("BLASgemv",BLASgemv_("C",&lds_,&k_,&sone,pS,&lds_,x,&one,&szero,cc,&one));
303: /* s_k = s_k - S_{0:k-1} cc */
304: if (bv->indef) { BV_ApplySignature(bv,k,cc,PETSC_TRUE); }
305: PetscStackCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&k_,&sonem,pS,&lds_,cc,&one,&sone,pS+k*lds,&one));
306: if (bv->indef) { BV_ApplySignature(bv,k,cc,PETSC_FALSE); }
307: }
309: if (norm) { BVTensorNormColumn(bv,k,norm); }
310: BV_AddCoefficients(bv,k,h,cc);
311: MatDenseRestoreArray(ctx->S,&pS);
312: VecRestoreArray(bv->buffer,&cc);
313: return(0);
314: }
316: PetscErrorCode BVView_Tensor(BV bv,PetscViewer viewer)
317: {
318: PetscErrorCode ierr;
319: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
320: PetscViewerFormat format;
321: PetscBool isascii;
322: const char *bvname,*uname,*sname;
325: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
326: if (isascii) {
327: PetscViewerGetFormat(viewer,&format);
328: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
329: PetscViewerASCIIPrintf(viewer,"number of tensor blocks (degree): %D\n",ctx->d);
330: PetscViewerASCIIPrintf(viewer,"number of columns of U factor: %D\n",ctx->ld);
331: return(0);
332: }
333: BVView(ctx->U,viewer);
334: MatView(ctx->S,viewer);
335: if (format == PETSC_VIEWER_ASCII_MATLAB) {
336: PetscObjectGetName((PetscObject)bv,&bvname);
337: PetscObjectGetName((PetscObject)ctx->U,&uname);
338: PetscObjectGetName((PetscObject)ctx->S,&sname);
339: PetscViewerASCIIPrintf(viewer,"%s=kron(eye(%D),%s)*%s(:,1:%D);\n",bvname,ctx->d,uname,sname,bv->k);
340: }
341: } else {
342: BVView(ctx->U,viewer);
343: MatView(ctx->S,viewer);
344: }
345: return(0);
346: }
348: static PetscErrorCode BVTensorUpdateMatrix(BV V,PetscInt ini,PetscInt end)
349: {
351: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
352: PetscInt i,j,r,c,l,k,ld=ctx->ld,lds=ctx->d*ctx->ld;
353: PetscScalar *qB,*sqB;
354: Vec u;
355: Mat A;
358: if (!V->matrix) return(0);
359: l = ctx->U->l; k = ctx->U->k;
360: /* update inner product matrix */
361: if (!ctx->qB) {
362: PetscCalloc2(lds*lds,&ctx->qB,lds,&ctx->sw);
363: VecDuplicate(ctx->U->t,&ctx->u);
364: }
365: ctx->U->l = 0;
366: for (r=0;r<ctx->d;r++) {
367: for (c=0;c<=r;c++) {
368: MatNestGetSubMat(V->matrix,r,c,&A);
369: if (A) {
370: qB = ctx->qB+c*ld*lds+r*ld;
371: for (i=ini;i<end;i++) {
372: BVGetColumn(ctx->U,i,&u);
373: MatMult(A,u,ctx->u);
374: ctx->U->k = i+1;
375: BVDotVec(ctx->U,ctx->u,qB+i*lds);
376: BVRestoreColumn(ctx->U,i,&u);
377: for (j=0;j<i;j++) qB[i+j*lds] = PetscConj(qB[j+i*lds]);
378: }
379: if (c!=r) {
380: sqB = ctx->qB+r*ld*lds+c*ld;
381: for (i=ini;i<end;i++) for (j=0;j<=i;j++) {
382: sqB[i+j*lds] = PetscConj(qB[j+i*lds]);
383: sqB[j+i*lds] = qB[j+i*lds];
384: }
385: }
386: }
387: }
388: }
389: ctx->U->l = l; ctx->U->k = k;
390: return(0);
391: }
393: static PetscErrorCode BVTensorBuildFirstColumn_Tensor(BV V,PetscInt k)
394: {
396: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
397: PetscInt i,nq=0;
398: PetscScalar *pS,*omega;
399: PetscReal norm;
400: PetscBool breakdown=PETSC_FALSE;
403: MatDenseGetArray(ctx->S,&pS);
404: for (i=0;i<ctx->d;i++) {
405: if (i>=k) {
406: BVSetRandomColumn(ctx->U,nq);
407: } else {
408: BVCopyColumn(ctx->U,i,nq);
409: }
410: BVOrthogonalizeColumn(ctx->U,nq,pS+i*ctx->ld,&norm,&breakdown);
411: if (!breakdown) {
412: BVScaleColumn(ctx->U,nq,1.0/norm);
413: pS[nq+i*ctx->ld] = norm;
414: nq++;
415: }
416: }
417: MatDenseRestoreArray(ctx->S,&pS);
418: if (!nq) SETERRQ1(PetscObjectComm((PetscObject)V),1,"Cannot build first column of tensor BV; U should contain k=%D nonzero columns",k);
419: BVTensorUpdateMatrix(V,0,nq);
420: BVTensorNormColumn(V,0,&norm);
421: BVScale_Tensor(V,0,1.0/norm);
422: if (V->indef) {
423: BV_AllocateSignature(V);
424: VecGetArray(V->omega,&omega);
425: omega[0] = (norm<0.0)? -1.0: 1.0;
426: VecRestoreArray(V->omega,&omega);
427: }
428: /* set active columns */
429: ctx->U->l = 0;
430: ctx->U->k = nq;
431: return(0);
432: }
434: /*@
435: BVTensorBuildFirstColumn - Builds the first column of the tensor basis vectors
436: V from the data contained in the first k columns of U.
438: Collective on BV
440: Input Parameters:
441: + V - the basis vectors context
442: - k - the number of columns of U with relevant information
444: Notes:
445: At most d columns are considered, where d is the degree of the tensor BV.
446: Given V = (I otimes U) S, this function computes the first column of V, that
447: is, it computes the coefficients of the first column of S by orthogonalizing
448: the first d columns of U. If k is less than d (or linearly dependent columns
449: are found) then additional random columns are used.
451: The computed column has unit norm.
453: Level: advanced
455: .seealso: BVCreateTensor()
456: @*/
457: PetscErrorCode BVTensorBuildFirstColumn(BV V,PetscInt k)
458: {
464: PetscUseMethod(V,"BVTensorBuildFirstColumn_C",(BV,PetscInt),(V,k));
465: return(0);
466: }
468: static PetscErrorCode BVTensorCompress_Tensor(BV V,PetscInt newc)
469: {
470: #if defined(PETSC_MISSING_LAPACK_GESVD) || defined(PETSC_MISSING_LAPACK_GEQRF) || defined(PETSC_MISSING_LAPACK_ORGQR)
472: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"GESVD/GEQRF/ORGQR - Lapack routine is unavailable");
473: #else
475: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
476: PetscInt nwu=0,nnc,nrow,lwa,r,c;
477: PetscInt i,j,k,n,lds=ctx->ld*ctx->d,deg=ctx->d,lock,cs1=V->k,rs1=ctx->U->k,rk=0,offu;
478: PetscScalar *S,*M,*Z,*pQ,*SS,*SS2,t,sone=1.0,zero=0.0,mone=-1.0,*p,*tau,*work,*qB,*sqB;
479: PetscReal *sg,tol,*rwork;
480: PetscBLASInt ld_,cs1_,rs1_,cs1tdeg,n_,info,lw_,newc_,newctdeg,nnc_,nrow_,nnctdeg,lds_,rk_;
481: Mat Q,A;
484: if (!cs1) return(0);
485: lwa = 6*ctx->ld*lds+2*cs1;
486: n = PetscMin(rs1,deg*cs1);
487: lock = ctx->U->l;
488: nnc = cs1-lock-newc;
489: nrow = rs1-lock;
490: PetscCalloc6(deg*newc*nnc,&SS,newc*nnc,&SS2,(rs1+lock+newc)*n,&pQ,deg*rs1,&tau,lwa,&work,6*n,&rwork);
491: offu = lock*(rs1+1);
492: M = work+nwu;
493: nwu += rs1*cs1*deg;
494: sg = rwork;
495: Z = work+nwu;
496: nwu += deg*cs1*n;
497: PetscBLASIntCast(n,&n_);
498: PetscBLASIntCast(nnc,&nnc_);
499: PetscBLASIntCast(cs1,&cs1_);
500: PetscBLASIntCast(rs1,&rs1_);
501: PetscBLASIntCast(newc,&newc_);
502: PetscBLASIntCast(newc*deg,&newctdeg);
503: PetscBLASIntCast(nnc*deg,&nnctdeg);
504: PetscBLASIntCast(cs1*deg,&cs1tdeg);
505: PetscBLASIntCast(lwa-nwu,&lw_);
506: PetscBLASIntCast(nrow,&nrow_);
507: PetscBLASIntCast(lds,&lds_);
508: MatDenseGetArray(ctx->S,&S);
510: if (newc>0) {
511: /* truncate columns associated with new converged eigenpairs */
512: for (j=0;j<deg;j++) {
513: for (i=lock;i<lock+newc;i++) {
514: PetscMemcpy(M+(i-lock+j*newc)*nrow,S+i*lds+j*ctx->ld+lock,nrow*sizeof(PetscScalar));
515: }
516: }
517: #if !defined (PETSC_USE_COMPLEX)
518: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&newctdeg,M,&nrow_,sg,pQ+offu,&rs1_,Z,&n_,work+nwu,&lw_,&info));
519: #else
520: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&newctdeg,M,&nrow_,sg,pQ+offu,&rs1_,Z,&n_,work+nwu,&lw_,rwork+n,&info));
521: #endif
522: SlepcCheckLapackInfo("gesvd",info);
523: /* SVD has rank min(newc,nrow) */
524: rk = PetscMin(newc,nrow);
525: for (i=0;i<rk;i++) {
526: t = sg[i];
527: PetscStackCallBLAS("BLASscal",BLASscal_(&newctdeg,&t,Z+i,&n_));
528: }
529: for (i=0;i<deg;i++) {
530: for (j=lock;j<lock+newc;j++) {
531: PetscMemcpy(S+j*lds+i*ctx->ld+lock,Z+(newc*i+j-lock)*n,rk*sizeof(PetscScalar));
532: PetscMemzero(S+j*lds+i*ctx->ld+lock+rk,(ctx->ld-lock-rk)*sizeof(PetscScalar));
533: }
534: }
535: /*
536: update columns associated with non-converged vectors, orthogonalize
537: against pQ so that next M has rank nnc+d-1 insted of nrow+d-1
538: */
539: for (i=0;i<deg;i++) {
540: 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_));
541: 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_));
542: /* repeat orthogonalization step */
543: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&newc_,&nnc_,&nrow_,&sone,pQ+offu,&rs1_,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_,&zero,SS2,&newc_));
544: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&nrow_,&nnc_,&newc_,&mone,pQ+offu,&rs1_,SS2,&newc_,&sone,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_));
545: for (j=0;j<newc*nnc;j++) *(SS+i*newc*nnc+j) += SS2[j];
546: }
547: }
549: /* truncate columns associated with non-converged eigenpairs */
550: for (j=0;j<deg;j++) {
551: for (i=lock+newc;i<cs1;i++) {
552: PetscMemcpy(M+(i-lock-newc+j*nnc)*nrow,S+i*lds+j*ctx->ld+lock,nrow*sizeof(PetscScalar));
553: }
554: }
555: #if !defined (PETSC_USE_COMPLEX)
556: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&nnctdeg,M,&nrow_,sg,pQ+offu+newc*rs1,&rs1_,Z,&n_,work+nwu,&lw_,&info));
557: #else
558: PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&nnctdeg,M,&nrow_,sg,pQ+offu+newc*rs1,&rs1_,Z,&n_,work+nwu,&lw_,rwork+n,&info));
559: #endif
560: SlepcCheckLapackInfo("gesvd",info);
561: tol = PetscMax(rs1,deg*cs1)*PETSC_MACHINE_EPSILON*sg[0];
562: rk = 0;
563: for (i=0;i<PetscMin(nrow,nnctdeg);i++) if (sg[i]>tol) rk++;
564: rk = PetscMin(nnc+deg-1,rk);
565: /* the SVD has rank (at most) nnc+deg-1 */
566: for (i=0;i<rk;i++) {
567: t = sg[i];
568: PetscStackCallBLAS("BLASscal",BLASscal_(&nnctdeg,&t,Z+i,&n_));
569: }
570: /* update S */
571: PetscMemzero(S+cs1*lds,(V->m-cs1)*lds*sizeof(PetscScalar));
572: k = ctx->ld-lock-newc-rk;
573: for (i=0;i<deg;i++) {
574: for (j=lock+newc;j<cs1;j++) {
575: PetscMemcpy(S+j*lds+i*ctx->ld+lock+newc,Z+(nnc*i+j-lock-newc)*n,rk*sizeof(PetscScalar));
576: PetscMemzero(S+j*lds+i*ctx->ld+lock+newc+rk,k*sizeof(PetscScalar));
577: }
578: }
579: if (newc>0) {
580: for (i=0;i<deg;i++) {
581: p = SS+nnc*newc*i;
582: for (j=lock+newc;j<cs1;j++) {
583: for (k=0;k<newc;k++) S[j*lds+i*ctx->ld+lock+k] = *(p++);
584: }
585: }
586: }
588: /* orthogonalize pQ */
589: rk = rk+newc;
590: PetscBLASIntCast(rk,&rk_);
591: PetscBLASIntCast(cs1-lock,&nnc_);
592: PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&nrow_,&rk_,pQ+offu,&rs1_,tau,work+nwu,&lw_,&info));
593: SlepcCheckLapackInfo("geqrf",info);
594: for (i=0;i<deg;i++) {
595: PetscStackCallBLAS("BLAStrmm",BLAStrmm_("L","U","N","N",&rk_,&nnc_,&sone,pQ+offu,&rs1_,S+lock*lds+lock+i*ctx->ld,&lds_));
596: }
597: PetscStackCallBLAS("LAPACKorgqr",LAPACKorgqr_(&nrow_,&rk_,&rk_,pQ+offu,&rs1_,tau,work+nwu,&lw_,&info));
598: SlepcCheckLapackInfo("orgqr",info);
600: /* update vectors U(:,idx) = U*Q(:,idx) */
601: rk = rk+lock;
602: for (i=0;i<lock;i++) pQ[i*(1+rs1)] = 1.0;
603: MatCreateSeqDense(PETSC_COMM_SELF,rs1,rk,pQ,&Q);
604: ctx->U->k = rs1;
605: BVMultInPlace(ctx->U,Q,lock,rk);
606: MatDestroy(&Q);
608: if (ctx->qB) {
609: /* update matrix qB */
610: PetscBLASIntCast(ctx->ld,&ld_);
611: PetscBLASIntCast(rk,&rk_);
612: for (r=0;r<ctx->d;r++) {
613: for (c=0;c<=r;c++) {
614: MatNestGetSubMat(V->matrix,r,c,&A);
615: if (A) {
616: qB = ctx->qB+r*ctx->ld+c*ctx->ld*lds;
617: PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&rs1_,&rk_,&rs1_,&sone,qB,&lds_,pQ,&rs1_,&zero,work+nwu,&rs1_));
618: PetscStackCallBLAS("BLASgemm",BLASgemm_("C","N",&rk_,&rk_,&rs1_,&sone,pQ,&rs1_,work+nwu,&rs1_,&zero,qB,&lds_));
619: for (i=0;i<rk;i++) for (j=0;j<i;j++) qB[i+j*lds] = PetscConj(qB[j+i*lds]);
620: for (i=rk;i<ctx->ld;i++) {
621: PetscMemzero(qB+i*lds,ctx->ld*sizeof(PetscScalar));
622: }
623: for (i=0;i<rk;i++) {
624: PetscMemzero(qB+i*lds+rk,(ctx->ld-rk)*sizeof(PetscScalar));
625: }
626: if (c!=r) {
627: sqB = ctx->qB+r*ctx->ld*lds+c*ctx->ld;
628: for (i=0;i<ctx->ld;i++) for (j=0;j<ctx->ld;j++) sqB[i+j*lds] = PetscConj(qB[j+i*lds]);
629: }
630: }
631: }
632: }
633: }
635: /* free work space */
636: PetscFree6(SS,SS2,pQ,tau,work,rwork);
637: MatDenseRestoreArray(ctx->S,&S);
639: /* set active columns */
640: if (newc) ctx->U->l += newc;
641: ctx->U->k = rk;
642: return(0);
643: #endif
644: }
646: /*@
647: BVTensorCompress - Updates the U and S factors of the tensor basis vectors
648: object V by means of an SVD, removing redundant information.
650: Collective on BV
652: Input Parameters:
653: + V - the tensor basis vectors context
654: - newc - additional columns to be locked
656: Notes:
657: This function is typically used when restarting Krylov solvers. Truncating a
658: tensor BV V = (I otimes U) S to its leading columns amounts to keeping the
659: leading columns of S. However, to effectively reduce the size of the
660: decomposition, it is necessary to compress it in a way that fewer columns of
661: U are employed. This can be achieved by means of an update that involves the
662: SVD of the low-rank matrix [S_0 S_1 ... S_{d-1}], where S_i are the pieces of S.
664: If newc is nonzero, then newc columns are added to the leading columns of V.
665: This means that the corresponding columns of the U and S factors will remain
666: invariant in subsequent operations.
668: Level: advanced
670: .seealso: BVCreateTensor(), BVSetActiveColumns()
671: @*/
672: PetscErrorCode BVTensorCompress(BV V,PetscInt newc)
673: {
679: PetscUseMethod(V,"BVTensorCompress_C",(BV,PetscInt),(V,newc));
680: return(0);
681: }
683: static PetscErrorCode BVTensorGetDegree_Tensor(BV bv,PetscInt *d)
684: {
685: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
688: *d = ctx->d;
689: return(0);
690: }
692: /*@
693: BVTensorGetDegree - Returns the number of blocks (degree) of the tensor BV.
695: Not collective
697: Input Parameter:
698: . bv - the basis vectors context
700: Output Parameter:
701: . d - the degree
703: Level: advanced
705: .seealso: BVCreateTensor()
706: @*/
707: PetscErrorCode BVTensorGetDegree(BV bv,PetscInt *d)
708: {
714: PetscUseMethod(bv,"BVTensorGetDegree_C",(BV,PetscInt*),(bv,d));
715: return(0);
716: }
718: static PetscErrorCode BVTensorGetFactors_Tensor(BV V,BV *U,Mat *S)
719: {
720: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
723: ctx->puk = ctx->U->k;
724: if (U) *U = ctx->U;
725: if (S) *S = ctx->S;
726: return(0);
727: }
729: /*@C
730: BVTensorGetFactors - Returns the two factors involved in the definition of the
731: tensor basis vectors object, V = (I otimes U) S.
733: Logically Collective on BV
735: Input Parameter:
736: . V - the basis vectors context
738: Output Parameters:
739: + U - the BV factor
740: - S - the Mat factor
742: Notes:
743: The returned factors are references (not copies) of the internal factors,
744: so modifying them will change the tensor BV as well. Some operations of the
745: tensor BV assume that U has orthonormal columns, so if the user modifies U
746: this restriction must be taken into account.
748: The returned factors must not be destroyed. BVTensorRestoreFactors() must
749: be called when they are no longer needed.
751: Pass a NULL vector for any of the arguments that is not needed.
753: Level: advanced
755: .seealso: BVTensorRestoreFactors()
756: @*/
757: PetscErrorCode BVTensorGetFactors(BV V,BV *U,Mat *S)
758: {
760: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
764: if (ctx->puk>-1) SETERRQ(PetscObjectComm((PetscObject)V),PETSC_ERR_ORDER,"Previous call to BVTensonGetFactors without a BVTensorRestoreFactors call");
765: PetscUseMethod(V,"BVTensorGetFactors_C",(BV,BV*,Mat*),(V,U,S));
766: return(0);
767: }
769: static PetscErrorCode BVTensorRestoreFactors_Tensor(BV V,BV *U,Mat *S)
770: {
772: BV_TENSOR *ctx = (BV_TENSOR*)V->data;
775: PetscObjectStateIncrease((PetscObject)V);
776: if (U) *U = NULL;
777: if (S) *S = NULL;
778: BVTensorUpdateMatrix(V,ctx->puk,ctx->U->k);
779: ctx->puk = -1;
780: return(0);
781: }
783: /*@C
784: BVTensorRestoreFactors - Restore the two factors that were obtained with
785: BVTensorGetFactors().
787: Logically Collective on BV
789: Input Parameters:
790: + V - the basis vectors context
791: . U - the BV factor (or NULL)
792: - S - the Mat factor (or NULL)
794: Nots:
795: The arguments must match the corresponding call to BVTensorGetFactors().
797: Level: advanced
799: .seealso: BVTensorGetFactors()
800: @*/
801: PetscErrorCode BVTensorRestoreFactors(BV V,BV *U,Mat *S)
802: {
809: PetscUseMethod(V,"BVTensorRestoreFactors_C",(BV,BV*,Mat*),(V,U,S));
810: return(0);
811: }
813: PetscErrorCode BVDestroy_Tensor(BV bv)
814: {
816: BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
819: BVDestroy(&ctx->U);
820: MatDestroy(&ctx->S);
821: if (ctx->u) {
822: PetscFree2(ctx->qB,ctx->sw);
823: VecDestroy(&ctx->u);
824: }
825: VecDestroy(&bv->cv[0]);
826: VecDestroy(&bv->cv[1]);
827: PetscFree(bv->data);
828: PetscObjectComposeFunction((PetscObject)bv,"BVTensorBuildFirstColumn_C",NULL);
829: PetscObjectComposeFunction((PetscObject)bv,"BVTensorCompress_C",NULL);
830: PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetDegree_C",NULL);
831: PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetFactors_C",NULL);
832: PetscObjectComposeFunction((PetscObject)bv,"BVTensorRestoreFactors_C",NULL);
833: return(0);
834: }
836: SLEPC_EXTERN PetscErrorCode BVCreate_Tensor(BV bv)
837: {
839: BV_TENSOR *ctx;
842: PetscNewLog(bv,&ctx);
843: bv->data = (void*)ctx;
844: ctx->puk = -1;
846: VecDuplicateEmpty(bv->t,&bv->cv[0]);
847: VecDuplicateEmpty(bv->t,&bv->cv[1]);
849: bv->ops->mult = BVMult_Tensor;
850: bv->ops->multvec = BVMultVec_Tensor;
851: bv->ops->multinplace = BVMultInPlace_Tensor;
852: bv->ops->multinplacetrans = BVMultInPlaceTranspose_Tensor;
853: bv->ops->dot = BVDot_Tensor;
854: bv->ops->dotvec = BVDotVec_Tensor;
855: bv->ops->dotvec_local = BVDotVec_Local_Tensor;
856: bv->ops->scale = BVScale_Tensor;
857: bv->ops->norm = BVNorm_Tensor;
858: bv->ops->norm_local = BVNorm_Local_Tensor;
859: bv->ops->matmult = BVMatMult_Tensor;
860: bv->ops->copy = BVCopy_Tensor;
861: bv->ops->copycolumn = BVCopyColumn_Tensor;
862: bv->ops->resize = BVResize_Tensor;
863: bv->ops->getcolumn = BVGetColumn_Tensor;
864: bv->ops->restorecolumn = BVRestoreColumn_Tensor;
865: bv->ops->getarray = BVGetArray_Tensor;
866: bv->ops->restorearray = BVRestoreArray_Tensor;
867: bv->ops->getarrayread = BVGetArrayRead_Tensor;
868: bv->ops->restorearrayread = BVRestoreArrayRead_Tensor;
869: bv->ops->gramschmidt = BVOrthogonalizeGS1_Tensor;
870: bv->ops->destroy = BVDestroy_Tensor;
871: bv->ops->view = BVView_Tensor;
873: PetscObjectComposeFunction((PetscObject)bv,"BVTensorBuildFirstColumn_C",BVTensorBuildFirstColumn_Tensor);
874: PetscObjectComposeFunction((PetscObject)bv,"BVTensorCompress_C",BVTensorCompress_Tensor);
875: PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetDegree_C",BVTensorGetDegree_Tensor);
876: PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetFactors_C",BVTensorGetFactors_Tensor);
877: PetscObjectComposeFunction((PetscObject)bv,"BVTensorRestoreFactors_C",BVTensorRestoreFactors_Tensor);
878: return(0);
879: }
881: /*@
882: BVCreateTensor - Creates a tensor BV that is represented in compact form
883: as V = (I otimes U) S, where U has orthonormal columns.
885: Collective on BV
887: Input Parameters:
888: + U - a basis vectors object
889: - d - the number of blocks (degree) of the tensor BV
891: Output Parameter:
892: . V - the new basis vectors context
894: Notes:
895: The new basis vectors object is V = (I otimes U) S, where otimes denotes
896: the Kronecker product, I is the identity matrix of order d, and S is a
897: sequential matrix allocated internally. This compact representation is
898: used e.g. to represent the Krylov basis generated with the linearization
899: of a matrix polynomial of degree d.
901: The size of V (number of rows) is equal to d times n, where n is the size
902: of U. The dimensions of S are d times m rows and m-d+1 columns, where m is
903: the number of columns of U, so m should be at least d.
905: The communicator of V will be the same as U.
907: On input, the content of U is irrelevant. Alternatively, it may contain
908: some nonzero columns that will be used by BVTensorBuildFirstColumn().
910: Level: advanced
912: .seealso: BVTensorGetDegree(), BVTensorGetFactors(), BVTensorBuildFirstColumn()
913: @*/
914: PetscErrorCode BVCreateTensor(BV U,PetscInt d,BV *V)
915: {
917: PetscBool match;
918: PetscInt n,N,m;
919: BV_TENSOR *ctx;
924: PetscObjectTypeCompare((PetscObject)U,BVTENSOR,&match);
925: if (match) SETERRQ(PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"U cannot be of type tensor");
927: BVCreate(PetscObjectComm((PetscObject)U),V);
928: BVGetSizes(U,&n,&N,&m);
929: if (m<d) SETERRQ2(PetscObjectComm((PetscObject)U),PETSC_ERR_ARG_SIZ,"U has %D columns, it should have at least d=%D",m,d);
930: BVSetSizes(*V,d*n,d*N,m-d+1);
931: PetscObjectChangeTypeName((PetscObject)*V,BVTENSOR);
932: PetscLogEventBegin(BV_Create,*V,0,0,0);
933: BVCreate_Tensor(*V);
934: PetscLogEventEnd(BV_Create,*V,0,0,0);
936: ctx = (BV_TENSOR*)(*V)->data;
937: ctx->U = U;
938: ctx->d = d;
939: ctx->ld = m;
940: PetscObjectReference((PetscObject)U);
941: PetscLogObjectParent((PetscObject)*V,(PetscObject)U);
942: MatCreateSeqDense(PETSC_COMM_SELF,d*m,m-d+1,NULL,&ctx->S);
943: PetscLogObjectParent((PetscObject)*V,(PetscObject)ctx->S);
944: PetscObjectSetName((PetscObject)ctx->S,"S");
945: return(0);
946: }