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 1272 : static PetscErrorCode BVMultInPlace_Tensor(BV V,Mat Q,PetscInt s,PetscInt e)
29 : {
30 1272 : BV_TENSOR *ctx = (BV_TENSOR*)V->data;
31 1272 : PetscScalar *pS;
32 1272 : const PetscScalar *q;
33 1272 : PetscInt ldq,lds = ctx->ld*ctx->d;
34 :
35 1272 : PetscFunctionBegin;
36 1272 : PetscCall(MatDenseGetLDA(Q,&ldq));
37 1272 : PetscCall(MatDenseGetArray(ctx->S,&pS));
38 1272 : PetscCall(MatDenseGetArrayRead(Q,&q));
39 1272 : 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 1272 : PetscCall(MatDenseRestoreArrayRead(Q,&q));
41 1272 : PetscCall(MatDenseRestoreArray(ctx->S,&pS));
42 1272 : 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 15072 : static PetscErrorCode BVScale_Tensor(BV bv,PetscInt j,PetscScalar alpha)
84 : {
85 15072 : BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
86 15072 : PetscScalar *pS;
87 15072 : PetscInt lds = ctx->ld*ctx->d;
88 :
89 15072 : PetscFunctionBegin;
90 15072 : PetscCall(MatDenseGetArray(ctx->S,&pS));
91 15072 : if (PetscUnlikely(j<0)) PetscCall(BVScale_BLAS_Private(bv,(bv->k-bv->l)*lds,pS+(bv->nc+bv->l)*lds,alpha));
92 15072 : else PetscCall(BVScale_BLAS_Private(bv,lds,pS+(bv->nc+j)*lds,alpha));
93 15072 : PetscCall(MatDenseRestoreArray(ctx->S,&pS));
94 15072 : 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 1245 : static PetscErrorCode BVCopyColumn_Tensor(BV V,PetscInt j,PetscInt i)
112 : {
113 1245 : BV_TENSOR *ctx = (BV_TENSOR*)V->data;
114 1245 : PetscScalar *pS;
115 1245 : PetscInt lds = ctx->ld*ctx->d;
116 :
117 1245 : PetscFunctionBegin;
118 1245 : PetscCall(MatDenseGetArray(ctx->S,&pS));
119 1245 : PetscCall(PetscArraycpy(pS+(V->nc+i)*lds,pS+(V->nc+j)*lds,lds));
120 1245 : PetscCall(MatDenseRestoreArray(ctx->S,&pS));
121 1245 : PetscFunctionReturn(PETSC_SUCCESS);
122 : }
123 :
124 44817 : static PetscErrorCode BVTensorNormColumn(BV bv,PetscInt j,PetscReal *norm)
125 : {
126 44817 : BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
127 44817 : PetscBLASInt one=1,lds_;
128 44817 : PetscScalar sone=1.0,szero=0.0,*x,dot;
129 44817 : const PetscScalar *S;
130 44817 : PetscReal alpha=1.0,scale=0.0,aval;
131 44817 : PetscInt i,lds,ld=ctx->ld;
132 :
133 44817 : PetscFunctionBegin;
134 44817 : lds = ld*ctx->d;
135 44817 : PetscCall(MatDenseGetArrayRead(ctx->S,&S));
136 44817 : PetscCall(PetscBLASIntCast(lds,&lds_));
137 44817 : if (PetscUnlikely(ctx->qB)) {
138 4451 : x = ctx->sw;
139 4451 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&lds_,&sone,ctx->qB,&lds_,S+j*lds,&one,&szero,x,&one));
140 4451 : dot = PetscRealPart(BLASdot_(&lds_,S+j*lds,&one,x,&one));
141 4451 : PetscCall(BV_SafeSqrt(bv,dot,norm));
142 : } else {
143 : /* Compute *norm = BLASnrm2_(&lds_,S+j*lds,&one); */
144 40366 : if (lds==1) *norm = PetscAbsScalar(S[j*lds]);
145 : else {
146 8633370 : for (i=0;i<lds;i++) {
147 8593004 : aval = PetscAbsScalar(S[i+j*lds]);
148 8593004 : if (aval!=0.0) {
149 5928067 : if (PetscUnlikely(scale<aval)) {
150 257996 : alpha = 1.0 + alpha*PetscSqr(scale/aval);
151 257996 : scale = aval;
152 5670071 : } else alpha += PetscSqr(aval/scale);
153 : }
154 : }
155 40366 : *norm = scale*PetscSqrtReal(alpha);
156 : }
157 : }
158 44817 : PetscCall(MatDenseRestoreArrayRead(ctx->S,&S));
159 44817 : PetscFunctionReturn(PETSC_SUCCESS);
160 : }
161 :
162 28982 : static PetscErrorCode BVOrthogonalizeGS1_Tensor(BV bv,PetscInt k,Vec v,PetscBool *which,PetscScalar *h,PetscScalar *c,PetscReal *onorm,PetscReal *norm)
163 : {
164 28982 : BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
165 28982 : PetscScalar *pS,*cc,*x,dot,sonem=-1.0,sone=1.0,szero=0.0;
166 28982 : PetscInt i,lds = ctx->ld*ctx->d;
167 28982 : PetscBLASInt lds_,k_,one=1;
168 28982 : const PetscScalar *omega;
169 :
170 28982 : PetscFunctionBegin;
171 28982 : PetscCheck(!v,PetscObjectComm((PetscObject)bv),PETSC_ERR_SUP,"Orthogonalization against an external vector is not allowed in BVTENSOR");
172 28982 : PetscCall(MatDenseGetArray(ctx->S,&pS));
173 28982 : if (!c) PetscCall(VecGetArray(bv->buffer,&cc));
174 0 : else cc = c;
175 28982 : PetscCall(PetscBLASIntCast(lds,&lds_));
176 28982 : PetscCall(PetscBLASIntCast(k,&k_));
177 :
178 28982 : if (onorm) PetscCall(BVTensorNormColumn(bv,k,onorm));
179 :
180 28982 : if (ctx->qB) x = ctx->sw;
181 20150 : else x = pS+k*lds;
182 :
183 28982 : 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 28866 : 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 28866 : 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 28866 : if (PetscUnlikely(bv->indef)) PetscCall(BV_ApplySignature(bv,k,cc,PETSC_TRUE));
209 28866 : PetscCallBLAS("BLASgemv",BLASgemv_("N",&lds_,&k_,&sonem,pS,&lds_,cc,&one,&sone,pS+k*lds,&one));
210 28866 : if (PetscUnlikely(bv->indef)) PetscCall(BV_ApplySignature(bv,k,cc,PETSC_FALSE));
211 : }
212 :
213 28982 : if (norm) PetscCall(BVTensorNormColumn(bv,k,norm));
214 28982 : PetscCall(BV_AddCoefficients(bv,k,h,cc));
215 28982 : PetscCall(MatDenseRestoreArray(ctx->S,&pS));
216 28982 : PetscCall(VecRestoreArray(bv->buffer,&cc));
217 28982 : 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 5637 : static PetscErrorCode BVTensorUpdateMatrix(BV V,PetscInt ini,PetscInt end)
252 : {
253 5637 : BV_TENSOR *ctx = (BV_TENSOR*)V->data;
254 5637 : PetscInt i,j,r,c,l,k,ld=ctx->ld,lds=ctx->d*ctx->ld;
255 5637 : PetscScalar *qB,*sqB;
256 5637 : Vec u;
257 5637 : Mat A;
258 :
259 5637 : PetscFunctionBegin;
260 5637 : if (!V->matrix) PetscFunctionReturn(PETSC_SUCCESS);
261 4435 : l = ctx->U->l; k = ctx->U->k;
262 : /* update inner product matrix */
263 4435 : if (!ctx->qB) {
264 35 : PetscCall(PetscCalloc2(lds*lds,&ctx->qB,lds,&ctx->sw));
265 35 : PetscCall(BVCreateVec(ctx->U,&ctx->u));
266 : }
267 4435 : ctx->U->l = 0;
268 13305 : for (r=0;r<ctx->d;r++) {
269 22175 : for (c=0;c<=r;c++) {
270 13305 : PetscCall(MatNestGetSubMat(V->matrix,r,c,&A));
271 13305 : if (A) {
272 8898 : qB = ctx->qB+c*ld*lds+r*ld;
273 17918 : for (i=ini;i<end;i++) {
274 9020 : PetscCall(BVGetColumn(ctx->U,i,&u));
275 9020 : PetscCall(MatMult(A,u,ctx->u));
276 9020 : ctx->U->k = i+1;
277 9020 : PetscCall(BVDotVec(ctx->U,ctx->u,qB+i*lds));
278 9020 : PetscCall(BVRestoreColumn(ctx->U,i,&u));
279 208494 : for (j=0;j<i;j++) qB[i+j*lds] = PetscConj(qB[j+i*lds]);
280 9020 : qB[i*lds+i] = PetscRealPart(qB[i+i*lds]);
281 : }
282 8898 : 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 4435 : ctx->U->l = l; ctx->U->k = k;
293 4435 : PetscFunctionReturn(PETSC_SUCCESS);
294 : }
295 :
296 158 : static PetscErrorCode BVTensorBuildFirstColumn_Tensor(BV V,PetscInt k)
297 : {
298 158 : BV_TENSOR *ctx = (BV_TENSOR*)V->data;
299 158 : PetscInt i,nq=0;
300 158 : PetscScalar *pS,*omega;
301 158 : PetscReal norm;
302 158 : PetscBool breakdown=PETSC_FALSE;
303 :
304 158 : PetscFunctionBegin;
305 158 : PetscCall(MatDenseGetArray(ctx->S,&pS));
306 859 : for (i=0;i<ctx->d;i++) {
307 701 : if (i>=k) PetscCall(BVSetRandomColumn(ctx->U,nq));
308 29 : else PetscCall(BVCopyColumn(ctx->U,i,nq));
309 701 : PetscCall(BVOrthogonalizeColumn(ctx->U,nq,pS+i*ctx->ld,&norm,&breakdown));
310 701 : if (!breakdown) {
311 700 : PetscCall(BVScaleColumn(ctx->U,nq,1.0/norm));
312 700 : pS[nq+i*ctx->ld] = norm;
313 700 : nq++;
314 : }
315 : }
316 158 : PetscCall(MatDenseRestoreArray(ctx->S,&pS));
317 158 : 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 158 : PetscCall(BVTensorUpdateMatrix(V,0,nq));
319 158 : PetscCall(BVTensorNormColumn(V,0,&norm));
320 158 : PetscCall(BVScale_Tensor(V,0,1.0/norm));
321 158 : if (V->indef) {
322 35 : PetscCall(BV_AllocateSignature(V));
323 35 : PetscCall(VecGetArray(V->omega,&omega));
324 35 : omega[0] = (norm<0.0)? -1.0: 1.0;
325 35 : PetscCall(VecRestoreArray(V->omega,&omega));
326 : }
327 : /* set active columns */
328 158 : ctx->U->l = 0;
329 158 : ctx->U->k = nq;
330 158 : 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 158 : PetscErrorCode BVTensorBuildFirstColumn(BV V,PetscInt k)
357 : {
358 158 : PetscFunctionBegin;
359 158 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
360 474 : PetscValidLogicalCollectiveInt(V,k,2);
361 158 : PetscUseMethod(V,"BVTensorBuildFirstColumn_C",(BV,PetscInt),(V,k));
362 158 : PetscFunctionReturn(PETSC_SUCCESS);
363 : }
364 :
365 1230 : static PetscErrorCode BVTensorCompress_Tensor(BV V,PetscInt newc)
366 : {
367 1230 : BV_TENSOR *ctx = (BV_TENSOR*)V->data;
368 1230 : PetscInt nwu=0,nnc,nrow,lwa,r,c;
369 1230 : 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 1230 : PetscScalar *S,*M,*Z,*pQ,*SS,*SS2,t,sone=1.0,zero=0.0,mone=-1.0,*p,*tau,*work,*qB,*sqB;
371 1230 : PetscReal *sg,tol,*rwork;
372 1230 : PetscBLASInt ld_,cs1_,rs1_,cs1tdeg,n_,info,lw_,newc_,newctdeg,nnc_,nrow_,nnctdeg,lds_,rk_;
373 1230 : Mat Q,A;
374 :
375 1230 : PetscFunctionBegin;
376 1230 : if (!cs1) PetscFunctionReturn(PETSC_SUCCESS);
377 1226 : lwa = 6*ctx->ld*lds+2*cs1;
378 1226 : n = PetscMin(rs1,deg*cs1);
379 1226 : lock = ctx->U->l;
380 1226 : nnc = cs1-lock-newc;
381 1226 : nrow = rs1-lock;
382 1226 : PetscCall(PetscCalloc6(deg*newc*nnc,&SS,newc*nnc,&SS2,(rs1+lock+newc)*n,&pQ,deg*rs1,&tau,lwa,&work,6*n,&rwork));
383 1226 : offu = lock*(rs1+1);
384 1226 : M = work+nwu;
385 1226 : nwu += rs1*cs1*deg;
386 1226 : sg = rwork;
387 1226 : Z = work+nwu;
388 1226 : nwu += deg*cs1*n;
389 1226 : PetscCall(PetscBLASIntCast(n,&n_));
390 1226 : PetscCall(PetscBLASIntCast(nnc,&nnc_));
391 1226 : PetscCall(PetscBLASIntCast(cs1,&cs1_));
392 1226 : PetscCall(PetscBLASIntCast(rs1,&rs1_));
393 1226 : PetscCall(PetscBLASIntCast(newc,&newc_));
394 1226 : PetscCall(PetscBLASIntCast(newc*deg,&newctdeg));
395 1226 : PetscCall(PetscBLASIntCast(nnc*deg,&nnctdeg));
396 1226 : PetscCall(PetscBLASIntCast(cs1*deg,&cs1tdeg));
397 1226 : PetscCall(PetscBLASIntCast(lwa-nwu,&lw_));
398 1226 : PetscCall(PetscBLASIntCast(nrow,&nrow_));
399 1226 : PetscCall(PetscBLASIntCast(lds,&lds_));
400 1226 : PetscCall(MatDenseGetArray(ctx->S,&S));
401 :
402 1226 : if (newc>0) {
403 : /* truncate columns associated with new converged eigenpairs */
404 1993 : for (j=0;j<deg;j++) {
405 6867 : 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 299 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
408 : #if !defined (PETSC_USE_COMPLEX)
409 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&newctdeg,M,&nrow_,sg,pQ+offu,&rs1_,Z,&n_,work+nwu,&lw_,&info));
410 : #else
411 299 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&newctdeg,M,&nrow_,sg,pQ+offu,&rs1_,Z,&n_,work+nwu,&lw_,rwork+n,&info));
412 : #endif
413 299 : SlepcCheckLapackInfo("gesvd",info);
414 299 : PetscCall(PetscFPTrapPop());
415 : /* SVD has rank min(newc,nrow) */
416 299 : rk = PetscMin(newc,nrow);
417 1802 : for (i=0;i<rk;i++) {
418 1503 : t = sg[i];
419 1503 : PetscCallBLAS("BLASscal",BLASscal_(&newctdeg,&t,Z+i,&n_));
420 : }
421 1993 : for (i=0;i<deg;i++) {
422 6867 : for (j=lock;j<lock+newc;j++) {
423 5173 : PetscCall(PetscArraycpy(S+j*lds+i*ctx->ld+lock,Z+(newc*i+j-lock)*n,rk));
424 5173 : 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 1993 : for (i=0;i<deg;i++) {
432 1694 : PetscCallBLAS("BLASgemm",BLASgemm_("C","N",&newc_,&nnc_,&nrow_,&sone,pQ+offu,&rs1_,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_,&zero,PetscSafePointerPlusOffset(SS,i*nnc*newc),&newc_));
433 1694 : PetscCallBLAS("BLASgemm",BLASgemm_("N","N",&nrow_,&nnc_,&newc_,&mone,pQ+offu,&rs1_,PetscSafePointerPlusOffset(SS,i*nnc*newc),&newc_,&sone,S+(lock+newc)*lds+i*ctx->ld+lock,&lds_));
434 : /* repeat orthogonalization step */
435 1694 : 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 1694 : 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 19949 : 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 6516 : for (j=0;j<deg;j++) {
443 52905 : 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 1226 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
446 : #if !defined (PETSC_USE_COMPLEX)
447 : PetscCallBLAS("LAPACKgesvd",LAPACKgesvd_("S","S",&nrow_,&nnctdeg,M,&nrow_,sg,pQ+offu+newc*rs1,&rs1_,Z,&n_,work+nwu,&lw_,&info));
448 : #else
449 1226 : 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 1226 : SlepcCheckLapackInfo("gesvd",info);
452 1226 : PetscCall(PetscFPTrapPop());
453 1226 : tol = PetscMax(rs1,deg*cs1)*PETSC_MACHINE_EPSILON*sg[0];
454 1226 : rk = 0;
455 27718 : for (i=0;i<PetscMin(nrow,nnctdeg);i++) if (sg[i]>tol) rk++;
456 1226 : rk = PetscMin(nnc+deg-1,rk);
457 : /* the SVD has rank (at most) nnc+deg-1 */
458 17017 : for (i=0;i<rk;i++) {
459 15791 : t = sg[i];
460 15791 : PetscCallBLAS("BLASscal",BLASscal_(&nnctdeg,&t,Z+i,&n_));
461 : }
462 : /* update S */
463 1226 : PetscCall(PetscArrayzero(S+cs1*lds,(V->m-cs1)*lds));
464 1226 : k = ctx->ld-lock-newc-rk;
465 6516 : for (i=0;i<deg;i++) {
466 52905 : for (j=lock+newc;j<cs1;j++) {
467 47615 : PetscCall(PetscArraycpy(S+j*lds+i*ctx->ld+lock+newc,Z+(nnc*i+j-lock-newc)*n,rk));
468 47615 : PetscCall(PetscArrayzero(S+j*lds+i*ctx->ld+lock+newc+rk,k));
469 : }
470 : }
471 1226 : if (newc>0) {
472 1993 : for (i=0;i<deg;i++) {
473 1694 : p = PetscSafePointerPlusOffset(SS,i*nnc*newc);
474 11425 : for (j=lock+newc;j<cs1;j++) {
475 27986 : for (k=0;k<newc;k++) S[j*lds+i*ctx->ld+lock+k] = *(p++);
476 : }
477 : }
478 : }
479 :
480 : /* orthogonalize pQ */
481 1226 : rk = rk+newc;
482 1226 : PetscCall(PetscBLASIntCast(rk,&rk_));
483 1226 : PetscCall(PetscBLASIntCast(cs1-lock,&nnc_));
484 1226 : PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
485 1226 : PetscCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&nrow_,&rk_,pQ+offu,&rs1_,tau,work+nwu,&lw_,&info));
486 1226 : SlepcCheckLapackInfo("geqrf",info);
487 6516 : for (i=0;i<deg;i++) {
488 5290 : PetscCallBLAS("BLAStrmm",BLAStrmm_("L","U","N","N",&rk_,&nnc_,&sone,pQ+offu,&rs1_,S+lock*lds+lock+i*ctx->ld,&lds_));
489 : }
490 1226 : PetscCallBLAS("LAPACKorgqr",LAPACKorgqr_(&nrow_,&rk_,&rk_,pQ+offu,&rs1_,tau,work+nwu,&lw_,&info));
491 1226 : SlepcCheckLapackInfo("orgqr",info);
492 1226 : PetscCall(PetscFPTrapPop());
493 :
494 : /* update vectors U(:,idx) = U*Q(:,idx) */
495 1226 : rk = rk+lock;
496 2722 : for (i=0;i<lock;i++) pQ[i*(1+rs1)] = 1.0;
497 1226 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,rs1,rk,pQ,&Q));
498 1226 : ctx->U->k = rs1;
499 1226 : PetscCall(BVMultInPlace(ctx->U,Q,lock,rk));
500 1226 : PetscCall(MatDestroy(&Q));
501 :
502 1226 : 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 12470 : for (i=0;i<rk;i++) {
514 119305 : for (j=0;j<i;j++) qB[i+j*lds] = PetscConj(qB[j+i*lds]);
515 11731 : qB[i+i*lds] = PetscRealPart(qB[i+i*lds]);
516 : }
517 12170 : for (i=rk;i<ctx->ld;i++) PetscCall(PetscArrayzero(qB+i*lds,ctx->ld));
518 12470 : 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 1226 : PetscCall(PetscFree6(SS,SS2,pQ,tau,work,rwork));
530 1226 : PetscCall(MatDenseRestoreArray(ctx->S,&S));
531 :
532 : /* set active columns */
533 1226 : if (newc) ctx->U->l += newc;
534 1226 : ctx->U->k = rk;
535 1226 : 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 1230 : PetscErrorCode BVTensorCompress(BV V,PetscInt newc)
565 : {
566 1230 : PetscFunctionBegin;
567 1230 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
568 3690 : PetscValidLogicalCollectiveInt(V,newc,2);
569 1230 : PetscUseMethod(V,"BVTensorCompress_C",(BV,PetscInt),(V,newc));
570 1230 : PetscFunctionReturn(PETSC_SUCCESS);
571 : }
572 :
573 313 : static PetscErrorCode BVTensorGetDegree_Tensor(BV bv,PetscInt *d)
574 : {
575 313 : BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
576 :
577 313 : PetscFunctionBegin;
578 313 : *d = ctx->d;
579 313 : 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 313 : PetscErrorCode BVTensorGetDegree(BV bv,PetscInt *d)
598 : {
599 313 : PetscFunctionBegin;
600 313 : PetscValidHeaderSpecific(bv,BV_CLASSID,1);
601 313 : PetscAssertPointer(d,2);
602 313 : PetscUseMethod(bv,"BVTensorGetDegree_C",(BV,PetscInt*),(bv,d));
603 313 : PetscFunctionReturn(PETSC_SUCCESS);
604 : }
605 :
606 5479 : static PetscErrorCode BVTensorGetFactors_Tensor(BV V,BV *U,Mat *S)
607 : {
608 5479 : BV_TENSOR *ctx = (BV_TENSOR*)V->data;
609 :
610 5479 : PetscFunctionBegin;
611 5479 : PetscCheck(ctx->puk==-1,PetscObjectComm((PetscObject)V),PETSC_ERR_ORDER,"Previous call to BVTensonGetFactors without a BVTensorRestoreFactors call");
612 5479 : ctx->puk = ctx->U->k;
613 5479 : if (U) *U = ctx->U;
614 5479 : if (S) *S = ctx->S;
615 5479 : PetscFunctionReturn(PETSC_SUCCESS);
616 : }
617 :
618 : /*@
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 5479 : PetscErrorCode BVTensorGetFactors(BV V,BV *U,Mat *S)
647 : {
648 5479 : PetscFunctionBegin;
649 5479 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
650 5479 : PetscUseMethod(V,"BVTensorGetFactors_C",(BV,BV*,Mat*),(V,U,S));
651 5479 : PetscFunctionReturn(PETSC_SUCCESS);
652 : }
653 :
654 5479 : static PetscErrorCode BVTensorRestoreFactors_Tensor(BV V,BV *U,Mat *S)
655 : {
656 5479 : BV_TENSOR *ctx = (BV_TENSOR*)V->data;
657 :
658 5479 : PetscFunctionBegin;
659 5479 : PetscCall(PetscObjectStateIncrease((PetscObject)V));
660 5479 : if (U) *U = NULL;
661 5479 : if (S) *S = NULL;
662 5479 : PetscCall(BVTensorUpdateMatrix(V,ctx->puk,ctx->U->k));
663 5479 : ctx->puk = -1;
664 5479 : PetscFunctionReturn(PETSC_SUCCESS);
665 : }
666 :
667 : /*@
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 5479 : PetscErrorCode BVTensorRestoreFactors(BV V,BV *U,Mat *S)
686 : {
687 5479 : PetscFunctionBegin;
688 5479 : PetscValidHeaderSpecific(V,BV_CLASSID,1);
689 5479 : if (U) PetscValidHeaderSpecific(*U,BV_CLASSID,2);
690 5479 : if (S) PetscValidHeaderSpecific(*S,MAT_CLASSID,3);
691 5479 : PetscUseMethod(V,"BVTensorRestoreFactors_C",(BV,BV*,Mat*),(V,U,S));
692 5479 : PetscFunctionReturn(PETSC_SUCCESS);
693 : }
694 :
695 166 : static PetscErrorCode BVDestroy_Tensor(BV bv)
696 : {
697 166 : BV_TENSOR *ctx = (BV_TENSOR*)bv->data;
698 :
699 166 : PetscFunctionBegin;
700 166 : PetscCall(BVDestroy(&ctx->U));
701 166 : PetscCall(MatDestroy(&ctx->S));
702 166 : if (ctx->u) {
703 35 : PetscCall(PetscFree2(ctx->qB,ctx->sw));
704 35 : PetscCall(VecDestroy(&ctx->u));
705 : }
706 166 : PetscCall(PetscFree(bv->data));
707 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorBuildFirstColumn_C",NULL));
708 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorCompress_C",NULL));
709 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetDegree_C",NULL));
710 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetFactors_C",NULL));
711 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorRestoreFactors_C",NULL));
712 166 : PetscFunctionReturn(PETSC_SUCCESS);
713 : }
714 :
715 166 : SLEPC_EXTERN PetscErrorCode BVCreate_Tensor(BV bv)
716 : {
717 166 : BV_TENSOR *ctx;
718 :
719 166 : PetscFunctionBegin;
720 166 : PetscCall(PetscNew(&ctx));
721 166 : bv->data = (void*)ctx;
722 166 : ctx->puk = -1;
723 :
724 166 : bv->ops->multinplace = BVMultInPlace_Tensor;
725 166 : bv->ops->multinplacetrans = BVMultInPlaceHermitianTranspose_Tensor;
726 166 : bv->ops->dot = BVDot_Tensor;
727 166 : bv->ops->scale = BVScale_Tensor;
728 166 : bv->ops->norm = BVNorm_Tensor;
729 166 : bv->ops->copycolumn = BVCopyColumn_Tensor;
730 166 : bv->ops->gramschmidt = BVOrthogonalizeGS1_Tensor;
731 166 : bv->ops->destroy = BVDestroy_Tensor;
732 166 : bv->ops->view = BVView_Tensor;
733 :
734 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorBuildFirstColumn_C",BVTensorBuildFirstColumn_Tensor));
735 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorCompress_C",BVTensorCompress_Tensor));
736 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetDegree_C",BVTensorGetDegree_Tensor));
737 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorGetFactors_C",BVTensorGetFactors_Tensor));
738 166 : PetscCall(PetscObjectComposeFunction((PetscObject)bv,"BVTensorRestoreFactors_C",BVTensorRestoreFactors_Tensor));
739 166 : 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 166 : PetscErrorCode BVCreateTensor(BV U,PetscInt d,BV *V)
776 : {
777 166 : PetscBool match;
778 166 : PetscInt n,N,m;
779 166 : VecType vtype;
780 166 : BV_TENSOR *ctx;
781 :
782 166 : PetscFunctionBegin;
783 166 : PetscValidHeaderSpecific(U,BV_CLASSID,1);
784 498 : PetscValidLogicalCollectiveInt(U,d,2);
785 166 : PetscCall(PetscObjectTypeCompare((PetscObject)U,BVTENSOR,&match));
786 166 : PetscCheck(!match,PetscObjectComm((PetscObject)U),PETSC_ERR_SUP,"U cannot be of type tensor");
787 :
788 166 : PetscCall(BVCreate(PetscObjectComm((PetscObject)U),V));
789 166 : PetscCall(BVGetSizes(U,&n,&N,&m));
790 166 : 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 166 : PetscCall(BVSetSizes(*V,d*n,d*N,m-d+1));
792 166 : PetscCall(BVGetVecType(U,&vtype));
793 166 : PetscCall(BVSetVecType(*V,vtype));
794 166 : PetscCall(PetscObjectChangeTypeName((PetscObject)*V,BVTENSOR));
795 166 : PetscCall(PetscLogEventBegin(BV_Create,*V,0,0,0));
796 166 : PetscCall(BVCreate_Tensor(*V));
797 166 : PetscCall(PetscLogEventEnd(BV_Create,*V,0,0,0));
798 :
799 166 : ctx = (BV_TENSOR*)(*V)->data;
800 166 : ctx->U = U;
801 166 : ctx->d = d;
802 166 : ctx->ld = m;
803 166 : PetscCall(PetscObjectReference((PetscObject)U));
804 166 : PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,d*m,m-d+1,NULL,&ctx->S));
805 166 : PetscCall(PetscObjectSetName((PetscObject)ctx->S,"S"));
806 :
807 : /* Copy user-provided attributes of U */
808 166 : (*V)->orthog_type = U->orthog_type;
809 166 : (*V)->orthog_ref = U->orthog_ref;
810 166 : (*V)->orthog_eta = U->orthog_eta;
811 166 : (*V)->orthog_block = U->orthog_block;
812 166 : (*V)->vmm = U->vmm;
813 166 : (*V)->rrandom = U->rrandom;
814 166 : PetscFunctionReturn(PETSC_SUCCESS);
815 : }
|
