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 : SLEPc singular value solver: "trlanczos"
12 :
13 : Method: Thick-restart Lanczos
14 :
15 : Algorithm:
16 :
17 : Golub-Kahan-Lanczos bidiagonalization with thick-restart.
18 :
19 : References:
20 :
21 : [1] G.H. Golub and W. Kahan, "Calculating the singular values
22 : and pseudo-inverse of a matrix", SIAM J. Numer. Anal. Ser.
23 : B 2:205-224, 1965.
24 :
25 : [2] V. Hernandez, J.E. Roman, and A. Tomas, "A robust and
26 : efficient parallel SVD solver based on restarted Lanczos
27 : bidiagonalization", Elec. Trans. Numer. Anal. 31:68-85,
28 : 2008.
29 : */
30 :
31 : #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
32 : #include <slepc/private/bvimpl.h>
33 :
34 : static PetscBool cited = PETSC_FALSE,citedg = PETSC_FALSE;
35 : static const char citation[] =
36 : "@Article{slepc-svd,\n"
37 : " author = \"V. Hern{\\'a}ndez and J. E. Rom{\\'a}n and A. Tom{\\'a}s\",\n"
38 : " title = \"A robust and efficient parallel {SVD} solver based on restarted {Lanczos} bidiagonalization\",\n"
39 : " journal = \"Electron. Trans. Numer. Anal.\",\n"
40 : " volume = \"31\",\n"
41 : " pages = \"68--85\",\n"
42 : " year = \"2008\"\n"
43 : "}\n";
44 : static const char citationg[] =
45 : "@Article{slepc-gsvd,\n"
46 : " author = \"F. Alvarruiz and C. Campos and J. E. Roman\",\n"
47 : " title = \"Thick-restarted {Lanczos} bidiagonalization methods for the {GSVD}\",\n"
48 : " journal = \"J. Comput. Appl. Math.\",\n"
49 : " volume = \"440\",\n"
50 : " pages = \"115506\",\n"
51 : " year = \"2024\"\n"
52 : "}\n";
53 :
54 : typedef struct {
55 : /* user parameters */
56 : PetscBool oneside; /* one-sided variant */
57 : PetscReal keep; /* restart parameter */
58 : PetscBool lock; /* locking/non-locking variant */
59 : KSP ksp; /* solver for least-squares problem in GSVD */
60 : SVDTRLanczosGBidiag bidiag; /* bidiagonalization variant for GSVD */
61 : PetscReal scalef; /* scale factor for matrix B */
62 : PetscReal scaleth; /* scale threshold for automatic scaling */
63 : PetscBool explicitmatrix;
64 : /* auxiliary variables */
65 : Mat Z; /* aux matrix for GSVD, Z=[A;B] */
66 : } SVD_TRLANCZOS;
67 :
68 : /* Context for shell matrix [A; B] */
69 : typedef struct {
70 : Mat A,B,AT,BT;
71 : Vec y1,y2,y;
72 : PetscInt m;
73 : PetscReal scalef;
74 : } MatZData;
75 :
76 55 : static PetscErrorCode MatZCreateContext(SVD svd,MatZData **zdata)
77 : {
78 55 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
79 :
80 55 : PetscFunctionBegin;
81 55 : PetscCall(PetscNew(zdata));
82 55 : (*zdata)->A = svd->A;
83 55 : (*zdata)->B = svd->B;
84 55 : (*zdata)->AT = svd->AT;
85 55 : (*zdata)->BT = svd->BT;
86 55 : (*zdata)->scalef = lanczos->scalef;
87 55 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&(*zdata)->y1));
88 55 : PetscCall(MatCreateVecsEmpty(svd->B,NULL,&(*zdata)->y2));
89 55 : PetscCall(VecGetLocalSize((*zdata)->y1,&(*zdata)->m));
90 55 : PetscCall(BVCreateVec(svd->U,&(*zdata)->y));
91 55 : PetscFunctionReturn(PETSC_SUCCESS);
92 : }
93 :
94 : /* Update scale factor for B in Z=[A;B]
95 : If matrices are swapped, the scale factor is inverted.*/
96 78 : static PetscErrorCode MatZUpdateScale(SVD svd)
97 : {
98 78 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
99 78 : MatZData *zdata;
100 78 : Mat mats[2],normal;
101 78 : MatType Atype;
102 78 : PetscBool sametype;
103 78 : PetscReal scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
104 :
105 78 : PetscFunctionBegin;
106 78 : if (lanczos->explicitmatrix) {
107 : /* Destroy the matrix Z and create it again */
108 13 : PetscCall(MatDestroy(&lanczos->Z));
109 13 : mats[0] = svd->A;
110 13 : if (scalef == 1.0) {
111 8 : mats[1] = svd->B;
112 : } else {
113 5 : PetscCall(MatDuplicate(svd->B,MAT_COPY_VALUES,&mats[1]));
114 5 : PetscCall(MatScale(mats[1],scalef));
115 : }
116 13 : PetscCall(MatCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,1,NULL,mats,&lanczos->Z));
117 13 : PetscCall(MatGetType(svd->A,&Atype));
118 13 : PetscCall(PetscObjectTypeCompare((PetscObject)svd->B,Atype,&sametype));
119 13 : if (!sametype) Atype = MATAIJ;
120 13 : PetscCall(MatConvert(lanczos->Z,Atype,MAT_INPLACE_MATRIX,&lanczos->Z));
121 13 : if (scalef != 1.0) PetscCall(MatDestroy(&mats[1]));
122 : } else {
123 65 : PetscCall(MatShellGetContext(lanczos->Z,&zdata));
124 65 : zdata->scalef = scalef;
125 : }
126 :
127 : /* create normal equations matrix, to build the preconditioner in LSQR */
128 78 : PetscCall(MatCreateNormalHermitian(lanczos->Z,&normal));
129 :
130 78 : if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
131 78 : PetscCall(SVD_KSPSetOperators(lanczos->ksp,lanczos->Z,normal));
132 78 : PetscCall(KSPSetUp(lanczos->ksp));
133 78 : PetscCall(MatDestroy(&normal));
134 78 : PetscFunctionReturn(PETSC_SUCCESS);
135 : }
136 :
137 55 : static PetscErrorCode MatDestroy_Z(Mat Z)
138 : {
139 55 : MatZData *zdata;
140 :
141 55 : PetscFunctionBegin;
142 55 : PetscCall(MatShellGetContext(Z,&zdata));
143 55 : PetscCall(VecDestroy(&zdata->y1));
144 55 : PetscCall(VecDestroy(&zdata->y2));
145 55 : PetscCall(VecDestroy(&zdata->y));
146 55 : PetscCall(PetscFree(zdata));
147 55 : PetscFunctionReturn(PETSC_SUCCESS);
148 : }
149 :
150 132465 : static PetscErrorCode MatMult_Z(Mat Z,Vec x,Vec y)
151 : {
152 132465 : MatZData *zdata;
153 132465 : PetscScalar *y_elems;
154 :
155 132465 : PetscFunctionBegin;
156 132465 : PetscCall(MatShellGetContext(Z,&zdata));
157 132465 : PetscCall(VecGetArray(y,&y_elems));
158 132465 : PetscCall(VecPlaceArray(zdata->y1,y_elems));
159 132465 : PetscCall(VecPlaceArray(zdata->y2,y_elems+zdata->m));
160 :
161 132465 : PetscCall(MatMult(zdata->A,x,zdata->y1));
162 132465 : PetscCall(MatMult(zdata->B,x,zdata->y2));
163 132465 : PetscCall(VecScale(zdata->y2,zdata->scalef));
164 :
165 132465 : PetscCall(VecResetArray(zdata->y1));
166 132465 : PetscCall(VecResetArray(zdata->y2));
167 132465 : PetscCall(VecRestoreArray(y,&y_elems));
168 132465 : PetscFunctionReturn(PETSC_SUCCESS);
169 : }
170 :
171 132679 : static PetscErrorCode MatMultTranspose_Z(Mat Z,Vec y,Vec x)
172 : {
173 132679 : MatZData *zdata;
174 132679 : const PetscScalar *y_elems;
175 :
176 132679 : PetscFunctionBegin;
177 132679 : PetscCall(MatShellGetContext(Z,&zdata));
178 132679 : PetscCall(VecGetArrayRead(y,&y_elems));
179 132679 : PetscCall(VecPlaceArray(zdata->y1,y_elems));
180 132679 : PetscCall(VecPlaceArray(zdata->y2,y_elems+zdata->m));
181 :
182 132679 : PetscCall(MatMult(zdata->BT,zdata->y2,x));
183 132679 : PetscCall(VecScale(x,zdata->scalef));
184 132679 : PetscCall(MatMultAdd(zdata->AT,zdata->y1,x,x));
185 :
186 132679 : PetscCall(VecResetArray(zdata->y1));
187 132679 : PetscCall(VecResetArray(zdata->y2));
188 132679 : PetscCall(VecRestoreArrayRead(y,&y_elems));
189 132679 : PetscFunctionReturn(PETSC_SUCCESS);
190 : }
191 :
192 163 : static PetscErrorCode MatCreateVecs_Z(Mat Z,Vec *right,Vec *left)
193 : {
194 163 : MatZData *zdata;
195 :
196 163 : PetscFunctionBegin;
197 163 : PetscCall(MatShellGetContext(Z,&zdata));
198 163 : if (right) PetscCall(MatCreateVecs(zdata->A,right,NULL));
199 163 : if (left) PetscCall(VecDuplicate(zdata->y,left));
200 163 : PetscFunctionReturn(PETSC_SUCCESS);
201 : }
202 :
203 : #define SWAP(a,b,t) do {t=a;a=b;b=t;} while (0)
204 :
205 108 : static PetscErrorCode SVDSetUp_TRLanczos(SVD svd)
206 : {
207 108 : PetscInt M,N,P,m,n,p;
208 108 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
209 108 : MatZData *zdata;
210 108 : Mat aux;
211 :
212 108 : PetscFunctionBegin;
213 108 : PetscCall(MatGetSize(svd->A,&M,&N));
214 108 : PetscCall(SVDSetDimensions_Default(svd));
215 108 : PetscCheck(svd->ncv<=svd->nsv+svd->mpd,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must not be larger than nsv+mpd");
216 108 : PetscCheck(lanczos->lock || svd->mpd>=svd->ncv,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Should not use mpd parameter in non-locking variant");
217 108 : if (svd->max_it==PETSC_DEFAULT) svd->max_it = PetscMax(N/svd->ncv,100);
218 108 : if (!lanczos->keep) lanczos->keep = 0.5;
219 108 : svd->leftbasis = PETSC_TRUE;
220 108 : PetscCall(SVDAllocateSolution(svd,1));
221 108 : if (svd->isgeneralized) {
222 64 : PetscCall(MatGetSize(svd->B,&P,NULL));
223 64 : if (lanczos->bidiag == SVD_TRLANCZOS_GBIDIAG_LOWER && ((svd->which==SVD_LARGEST && P<=N) || (svd->which==SVD_SMALLEST && M>N && P<=N))) {
224 5 : SWAP(svd->A,svd->B,aux);
225 5 : SWAP(svd->AT,svd->BT,aux);
226 5 : svd->swapped = PETSC_TRUE;
227 59 : } else svd->swapped = PETSC_FALSE;
228 :
229 64 : PetscCall(SVDSetWorkVecs(svd,1,1));
230 :
231 64 : if (svd->conv==SVD_CONV_ABS) { /* Residual norms are multiplied by matrix norms */
232 0 : if (!svd->nrma) PetscCall(MatNorm(svd->A,NORM_INFINITY,&svd->nrma));
233 0 : if (!svd->nrmb) PetscCall(MatNorm(svd->B,NORM_INFINITY,&svd->nrmb));
234 : }
235 :
236 : /* Create the matrix Z=[A;B] */
237 64 : PetscCall(MatGetLocalSize(svd->A,&m,&n));
238 64 : PetscCall(MatGetLocalSize(svd->B,&p,NULL));
239 64 : if (!lanczos->explicitmatrix) {
240 55 : PetscCall(MatDestroy(&lanczos->Z));
241 55 : PetscCall(MatZCreateContext(svd,&zdata));
242 55 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),m+p,n,PETSC_DECIDE,PETSC_DECIDE,zdata,&lanczos->Z));
243 55 : PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT,(void(*)(void))MatMult_Z));
244 : #if defined(PETSC_USE_COMPLEX)
245 : PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT_HERMITIAN_TRANSPOSE,(void(*)(void))MatMultTranspose_Z));
246 : #else
247 55 : PetscCall(MatShellSetOperation(lanczos->Z,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Z));
248 : #endif
249 55 : PetscCall(MatShellSetOperation(lanczos->Z,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_Z));
250 55 : PetscCall(MatShellSetOperation(lanczos->Z,MATOP_DESTROY,(void(*)(void))MatDestroy_Z));
251 : }
252 : /* Explicit matrix is created here, when updating the scale */
253 64 : PetscCall(MatZUpdateScale(svd));
254 :
255 44 : } else if (svd->ishyperbolic) {
256 14 : PetscCall(BV_SetMatrixDiagonal(svd->swapped?svd->V:svd->U,svd->omega,svd->OP));
257 14 : PetscCall(SVDSetWorkVecs(svd,1,0));
258 : }
259 108 : PetscCall(DSSetCompact(svd->ds,PETSC_TRUE));
260 108 : PetscCall(DSSetExtraRow(svd->ds,PETSC_TRUE));
261 108 : PetscCall(DSAllocate(svd->ds,svd->ncv+1));
262 108 : PetscFunctionReturn(PETSC_SUCCESS);
263 : }
264 :
265 20 : static PetscErrorCode SVDOneSideTRLanczosMGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
266 : {
267 20 : PetscReal a,b;
268 20 : PetscInt i,k=nconv+l;
269 20 : Vec ui,ui1,vi;
270 :
271 20 : PetscFunctionBegin;
272 20 : PetscCall(BVGetColumn(V,k,&vi));
273 20 : PetscCall(BVGetColumn(U,k,&ui));
274 20 : PetscCall(MatMult(svd->A,vi,ui));
275 20 : PetscCall(BVRestoreColumn(V,k,&vi));
276 20 : PetscCall(BVRestoreColumn(U,k,&ui));
277 20 : if (l>0) {
278 18 : PetscCall(BVSetActiveColumns(U,nconv,n));
279 102 : for (i=0;i<l;i++) work[i]=beta[i+nconv];
280 18 : PetscCall(BVMultColumn(U,-1.0,1.0,k,work));
281 : }
282 20 : PetscCall(BVNormColumn(U,k,NORM_2,&a));
283 20 : PetscCall(BVScaleColumn(U,k,1.0/a));
284 20 : alpha[k] = a;
285 :
286 104 : for (i=k+1;i<n;i++) {
287 84 : PetscCall(BVGetColumn(V,i,&vi));
288 84 : PetscCall(BVGetColumn(U,i-1,&ui1));
289 84 : PetscCall(MatMult(svd->AT,ui1,vi));
290 84 : PetscCall(BVRestoreColumn(V,i,&vi));
291 84 : PetscCall(BVRestoreColumn(U,i-1,&ui1));
292 84 : PetscCall(BVOrthonormalizeColumn(V,i,PETSC_FALSE,&b,NULL));
293 84 : beta[i-1] = b;
294 :
295 84 : PetscCall(BVGetColumn(V,i,&vi));
296 84 : PetscCall(BVGetColumn(U,i,&ui));
297 84 : PetscCall(MatMult(svd->A,vi,ui));
298 84 : PetscCall(BVRestoreColumn(V,i,&vi));
299 84 : PetscCall(BVGetColumn(U,i-1,&ui1));
300 84 : PetscCall(VecAXPY(ui,-b,ui1));
301 84 : PetscCall(BVRestoreColumn(U,i-1,&ui1));
302 84 : PetscCall(BVRestoreColumn(U,i,&ui));
303 84 : PetscCall(BVNormColumn(U,i,NORM_2,&a));
304 84 : PetscCall(BVScaleColumn(U,i,1.0/a));
305 84 : alpha[i] = a;
306 : }
307 :
308 20 : PetscCall(BVGetColumn(V,n,&vi));
309 20 : PetscCall(BVGetColumn(U,n-1,&ui1));
310 20 : PetscCall(MatMult(svd->AT,ui1,vi));
311 20 : PetscCall(BVRestoreColumn(V,n,&vi));
312 20 : PetscCall(BVRestoreColumn(U,n-1,&ui1));
313 20 : PetscCall(BVOrthogonalizeColumn(V,n,NULL,&b,NULL));
314 20 : beta[n-1] = b;
315 20 : PetscFunctionReturn(PETSC_SUCCESS);
316 : }
317 :
318 : /*
319 : Custom CGS orthogonalization, preprocess after first orthogonalization
320 : */
321 208 : static PetscErrorCode SVDOrthogonalizeCGS(BV V,PetscInt i,PetscScalar* h,PetscReal a,BVOrthogRefineType refine,PetscReal eta,PetscReal *norm)
322 : {
323 208 : PetscReal sum,onorm;
324 208 : PetscScalar dot;
325 208 : PetscInt j;
326 :
327 208 : PetscFunctionBegin;
328 208 : switch (refine) {
329 0 : case BV_ORTHOG_REFINE_NEVER:
330 0 : PetscCall(BVNormColumn(V,i,NORM_2,norm));
331 : break;
332 104 : case BV_ORTHOG_REFINE_ALWAYS:
333 104 : PetscCall(BVSetActiveColumns(V,0,i));
334 104 : PetscCall(BVDotColumn(V,i,h));
335 104 : PetscCall(BVMultColumn(V,-1.0,1.0,i,h));
336 104 : PetscCall(BVNormColumn(V,i,NORM_2,norm));
337 : break;
338 104 : case BV_ORTHOG_REFINE_IFNEEDED:
339 104 : dot = h[i];
340 104 : onorm = PetscSqrtReal(PetscRealPart(dot)) / a;
341 104 : sum = 0.0;
342 898 : for (j=0;j<i;j++) {
343 794 : sum += PetscRealPart(h[j] * PetscConj(h[j]));
344 : }
345 104 : *norm = PetscRealPart(dot)/(a*a) - sum;
346 104 : if (*norm>0.0) *norm = PetscSqrtReal(*norm);
347 0 : else PetscCall(BVNormColumn(V,i,NORM_2,norm));
348 104 : if (*norm < eta*onorm) {
349 104 : PetscCall(BVSetActiveColumns(V,0,i));
350 104 : PetscCall(BVDotColumn(V,i,h));
351 104 : PetscCall(BVMultColumn(V,-1.0,1.0,i,h));
352 104 : PetscCall(BVNormColumn(V,i,NORM_2,norm));
353 : }
354 : break;
355 : }
356 208 : PetscFunctionReturn(PETSC_SUCCESS);
357 : }
358 :
359 40 : static PetscErrorCode SVDOneSideTRLanczosCGS(SVD svd,PetscReal *alpha,PetscReal *beta,BV V,BV U,PetscInt nconv,PetscInt l,PetscInt n,PetscScalar* work)
360 : {
361 40 : PetscReal a,b,eta;
362 40 : PetscInt i,j,k=nconv+l;
363 40 : Vec ui,ui1,vi;
364 40 : BVOrthogRefineType refine;
365 :
366 40 : PetscFunctionBegin;
367 40 : PetscCall(BVGetColumn(V,k,&vi));
368 40 : PetscCall(BVGetColumn(U,k,&ui));
369 40 : PetscCall(MatMult(svd->A,vi,ui));
370 40 : PetscCall(BVRestoreColumn(V,k,&vi));
371 40 : PetscCall(BVRestoreColumn(U,k,&ui));
372 40 : if (l>0) {
373 36 : PetscCall(BVSetActiveColumns(U,nconv,n));
374 204 : for (i=0;i<l;i++) work[i]=beta[i+nconv];
375 36 : PetscCall(BVMultColumn(U,-1.0,1.0,k,work));
376 : }
377 40 : PetscCall(BVGetOrthogonalization(V,NULL,&refine,&eta,NULL));
378 :
379 208 : for (i=k+1;i<n;i++) {
380 168 : PetscCall(BVGetColumn(V,i,&vi));
381 168 : PetscCall(BVGetColumn(U,i-1,&ui1));
382 168 : PetscCall(MatMult(svd->AT,ui1,vi));
383 168 : PetscCall(BVRestoreColumn(V,i,&vi));
384 168 : PetscCall(BVRestoreColumn(U,i-1,&ui1));
385 168 : PetscCall(BVNormColumnBegin(U,i-1,NORM_2,&a));
386 168 : if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
387 84 : PetscCall(BVSetActiveColumns(V,0,i+1));
388 84 : PetscCall(BVGetColumn(V,i,&vi));
389 84 : PetscCall(BVDotVecBegin(V,vi,work));
390 : } else {
391 84 : PetscCall(BVSetActiveColumns(V,0,i));
392 84 : PetscCall(BVDotColumnBegin(V,i,work));
393 : }
394 168 : PetscCall(BVNormColumnEnd(U,i-1,NORM_2,&a));
395 168 : if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
396 84 : PetscCall(BVDotVecEnd(V,vi,work));
397 84 : PetscCall(BVRestoreColumn(V,i,&vi));
398 84 : PetscCall(BVSetActiveColumns(V,0,i));
399 84 : } else PetscCall(BVDotColumnEnd(V,i,work));
400 :
401 168 : PetscCall(BVScaleColumn(U,i-1,1.0/a));
402 1356 : for (j=0;j<i;j++) work[j] = work[j] / a;
403 168 : PetscCall(BVMultColumn(V,-1.0,1.0/a,i,work));
404 168 : PetscCall(SVDOrthogonalizeCGS(V,i,work,a,refine,eta,&b));
405 168 : PetscCall(BVScaleColumn(V,i,1.0/b));
406 168 : PetscCheck(PetscAbsReal(b)>10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)svd),PETSC_ERR_PLIB,"Recurrence generated a zero vector; use a two-sided variant");
407 :
408 168 : PetscCall(BVGetColumn(V,i,&vi));
409 168 : PetscCall(BVGetColumn(U,i,&ui));
410 168 : PetscCall(BVGetColumn(U,i-1,&ui1));
411 168 : PetscCall(MatMult(svd->A,vi,ui));
412 168 : PetscCall(VecAXPY(ui,-b,ui1));
413 168 : PetscCall(BVRestoreColumn(V,i,&vi));
414 168 : PetscCall(BVRestoreColumn(U,i,&ui));
415 168 : PetscCall(BVRestoreColumn(U,i-1,&ui1));
416 :
417 168 : alpha[i-1] = a;
418 168 : beta[i-1] = b;
419 : }
420 :
421 40 : PetscCall(BVGetColumn(V,n,&vi));
422 40 : PetscCall(BVGetColumn(U,n-1,&ui1));
423 40 : PetscCall(MatMult(svd->AT,ui1,vi));
424 40 : PetscCall(BVRestoreColumn(V,n,&vi));
425 40 : PetscCall(BVRestoreColumn(U,n-1,&ui1));
426 :
427 40 : PetscCall(BVNormColumnBegin(svd->U,n-1,NORM_2,&a));
428 40 : if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
429 20 : PetscCall(BVSetActiveColumns(V,0,n+1));
430 20 : PetscCall(BVGetColumn(V,n,&vi));
431 20 : PetscCall(BVDotVecBegin(V,vi,work));
432 : } else {
433 20 : PetscCall(BVSetActiveColumns(V,0,n));
434 20 : PetscCall(BVDotColumnBegin(V,n,work));
435 : }
436 40 : PetscCall(BVNormColumnEnd(svd->U,n-1,NORM_2,&a));
437 40 : if (refine == BV_ORTHOG_REFINE_IFNEEDED) {
438 20 : PetscCall(BVDotVecEnd(V,vi,work));
439 20 : PetscCall(BVRestoreColumn(V,n,&vi));
440 20 : } else PetscCall(BVDotColumnEnd(V,n,work));
441 :
442 40 : PetscCall(BVScaleColumn(U,n-1,1.0/a));
443 440 : for (j=0;j<n;j++) work[j] = work[j] / a;
444 40 : PetscCall(BVMultColumn(V,-1.0,1.0/a,n,work));
445 40 : PetscCall(SVDOrthogonalizeCGS(V,n,work,a,refine,eta,&b));
446 40 : PetscCall(BVSetActiveColumns(V,nconv,n));
447 40 : alpha[n-1] = a;
448 40 : beta[n-1] = b;
449 40 : PetscFunctionReturn(PETSC_SUCCESS);
450 : }
451 :
452 30 : static PetscErrorCode SVDSolve_TRLanczos(SVD svd)
453 : {
454 30 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
455 30 : PetscReal *alpha,*beta;
456 30 : PetscScalar *swork=NULL,*w;
457 30 : PetscInt i,k,l,nv,ld;
458 30 : Mat U,V;
459 30 : PetscBool breakdown=PETSC_FALSE;
460 30 : BVOrthogType orthog;
461 :
462 30 : PetscFunctionBegin;
463 30 : PetscCall(PetscCitationsRegister(citation,&cited));
464 : /* allocate working space */
465 30 : PetscCall(DSGetLeadingDimension(svd->ds,&ld));
466 30 : PetscCall(BVGetOrthogonalization(svd->V,&orthog,NULL,NULL,NULL));
467 30 : PetscCall(PetscMalloc1(ld,&w));
468 30 : if (lanczos->oneside) PetscCall(PetscMalloc1(svd->ncv+1,&swork));
469 :
470 : /* normalize start vector */
471 30 : if (!svd->nini) {
472 21 : PetscCall(BVSetRandomColumn(svd->V,0));
473 21 : PetscCall(BVOrthonormalizeColumn(svd->V,0,PETSC_TRUE,NULL,NULL));
474 : }
475 :
476 : l = 0;
477 518 : while (svd->reason == SVD_CONVERGED_ITERATING) {
478 488 : svd->its++;
479 :
480 : /* inner loop */
481 488 : nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
482 488 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
483 488 : beta = alpha + ld;
484 488 : if (lanczos->oneside) {
485 60 : if (orthog == BV_ORTHOG_MGS) PetscCall(SVDOneSideTRLanczosMGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork));
486 40 : else PetscCall(SVDOneSideTRLanczosCGS(svd,alpha,beta,svd->V,svd->U,svd->nconv,l,nv,swork));
487 428 : } else PetscCall(SVDTwoSideLanczos(svd,alpha,beta,svd->V,svd->U,svd->nconv+l,&nv,&breakdown));
488 488 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
489 488 : PetscCall(BVScaleColumn(svd->V,nv,1.0/beta[nv-1]));
490 488 : PetscCall(BVSetActiveColumns(svd->V,svd->nconv,nv));
491 488 : PetscCall(BVSetActiveColumns(svd->U,svd->nconv,nv));
492 :
493 : /* solve projected problem */
494 488 : PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
495 488 : PetscCall(DSSVDSetDimensions(svd->ds,nv));
496 488 : PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
497 488 : PetscCall(DSSolve(svd->ds,w,NULL));
498 488 : PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
499 488 : PetscCall(DSUpdateExtraRow(svd->ds));
500 488 : PetscCall(DSSynchronize(svd->ds,w,NULL));
501 5368 : for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
502 :
503 : /* check convergence */
504 488 : PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k));
505 488 : PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
506 :
507 : /* update l */
508 488 : if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
509 458 : else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
510 488 : if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
511 488 : if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
512 :
513 488 : if (svd->reason == SVD_CONVERGED_ITERATING) {
514 458 : if (PetscUnlikely(breakdown || k==nv)) {
515 : /* Start a new bidiagonalization */
516 0 : PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
517 0 : if (k<svd->nsv) {
518 0 : PetscCall(BVSetRandomColumn(svd->V,k));
519 0 : PetscCall(BVOrthonormalizeColumn(svd->V,k,PETSC_FALSE,NULL,&breakdown));
520 0 : if (breakdown) {
521 0 : svd->reason = SVD_DIVERGED_BREAKDOWN;
522 0 : PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
523 : }
524 : }
525 458 : } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
526 : }
527 :
528 : /* compute converged singular vectors and restart vectors */
529 488 : PetscCall(DSGetMat(svd->ds,DS_MAT_V,&V));
530 488 : PetscCall(BVMultInPlace(svd->V,V,svd->nconv,k+l));
531 488 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&V));
532 488 : PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
533 488 : PetscCall(BVMultInPlace(svd->U,U,svd->nconv,k+l));
534 488 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
535 :
536 : /* copy the last vector to be the next initial vector */
537 488 : if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(svd->V,nv,k+l));
538 :
539 488 : svd->nconv = k;
540 518 : PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
541 : }
542 :
543 : /* orthonormalize U columns in one side method */
544 30 : if (lanczos->oneside) {
545 33 : for (i=0;i<svd->nconv;i++) PetscCall(BVOrthonormalizeColumn(svd->U,i,PETSC_FALSE,NULL,NULL));
546 : }
547 :
548 : /* free working space */
549 30 : PetscCall(PetscFree(w));
550 30 : if (swork) PetscCall(PetscFree(swork));
551 30 : PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
552 30 : PetscFunctionReturn(PETSC_SUCCESS);
553 : }
554 :
555 376 : static PetscErrorCode SVDLanczosHSVD(SVD svd,PetscReal *alpha,PetscReal *beta,PetscReal *omega,Mat A,Mat AT,BV V,BV U,PetscInt k,PetscInt *n,PetscBool *breakdown)
556 : {
557 376 : PetscInt i;
558 376 : Vec u,v,ou=svd->workl[0];
559 376 : PetscBool lindep=PETSC_FALSE;
560 376 : PetscReal norm;
561 :
562 376 : PetscFunctionBegin;
563 5252 : for (i=k;i<*n;i++) {
564 4876 : PetscCall(BVGetColumn(V,i,&v));
565 4876 : PetscCall(BVGetColumn(U,i,&u));
566 4876 : PetscCall(MatMult(A,v,u));
567 4876 : PetscCall(BVRestoreColumn(V,i,&v));
568 4876 : PetscCall(BVRestoreColumn(U,i,&u));
569 4876 : PetscCall(BVOrthonormalizeColumn(U,i,PETSC_FALSE,alpha+i,&lindep));
570 4876 : omega[i] = PetscSign(alpha[i]);
571 4876 : if (PetscUnlikely(lindep)) {
572 0 : *n = i;
573 0 : break;
574 : }
575 :
576 4876 : PetscCall(BVGetColumn(V,i+1,&v));
577 4876 : PetscCall(BVGetColumn(U,i,&u));
578 4876 : PetscCall(VecPointwiseMult(ou,svd->omega,u));
579 4876 : PetscCall(MatMult(AT,ou,v));
580 4876 : PetscCall(BVRestoreColumn(V,i+1,&v));
581 4876 : PetscCall(BVRestoreColumn(U,i,&u));
582 4876 : PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,&norm,&lindep));
583 4876 : beta[i] = omega[i]*norm;
584 4876 : if (PetscUnlikely(lindep)) {
585 0 : *n = i+1;
586 0 : break;
587 : }
588 : }
589 :
590 376 : if (breakdown) *breakdown = lindep;
591 376 : PetscFunctionReturn(PETSC_SUCCESS);
592 : }
593 :
594 14 : static PetscErrorCode SVDSolve_TRLanczos_HSVD(SVD svd)
595 : {
596 14 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
597 14 : PetscReal *alpha,*beta,*omega;
598 14 : PetscScalar *w;
599 14 : PetscInt i,k,l,nv,ld,nini;
600 14 : Mat UU,VV,D,A,AT;
601 14 : BV U,V;
602 14 : PetscBool breakdown=PETSC_FALSE;
603 14 : BVOrthogType orthog;
604 14 : Vec vomega;
605 :
606 14 : PetscFunctionBegin;
607 : /* undo the effect of swapping in this function */
608 14 : if (svd->swapped) {
609 4 : A = svd->AT;
610 4 : AT = svd->A;
611 4 : U = svd->V;
612 4 : V = svd->U;
613 4 : nini = svd->ninil;
614 : } else {
615 10 : A = svd->A;
616 10 : AT = svd->AT;
617 10 : U = svd->U;
618 10 : V = svd->V;
619 10 : nini = svd->nini;
620 : }
621 : /* allocate working space */
622 14 : PetscCall(DSGetLeadingDimension(svd->ds,&ld));
623 14 : PetscCall(BVGetOrthogonalization(V,&orthog,NULL,NULL,NULL));
624 14 : PetscCall(PetscMalloc1(ld,&w));
625 14 : PetscCheck(!lanczos->oneside,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Oneside orthogonalization not supported for HSVD");
626 :
627 : /* normalize start vector */
628 14 : if (!nini) {
629 14 : PetscCall(BVSetRandomColumn(V,0));
630 14 : PetscCall(BVOrthonormalizeColumn(V,0,PETSC_TRUE,NULL,NULL));
631 : }
632 :
633 : l = 0;
634 390 : while (svd->reason == SVD_CONVERGED_ITERATING) {
635 376 : svd->its++;
636 :
637 : /* inner loop */
638 376 : nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
639 376 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
640 376 : beta = alpha + ld;
641 376 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&omega));
642 376 : PetscCall(SVDLanczosHSVD(svd,alpha,beta,omega,A,AT,V,U,svd->nconv+l,&nv,&breakdown));
643 376 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
644 376 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&omega));
645 376 : PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
646 376 : PetscCall(BVSetActiveColumns(U,svd->nconv,nv));
647 :
648 : /* solve projected problem */
649 376 : PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
650 376 : PetscCall(DSHSVDSetDimensions(svd->ds,nv));
651 376 : PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
652 376 : PetscCall(DSSolve(svd->ds,w,NULL));
653 376 : PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
654 376 : PetscCall(DSUpdateExtraRow(svd->ds));
655 376 : PetscCall(DSSynchronize(svd->ds,w,NULL));
656 376 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&omega));
657 9704 : for (i=svd->nconv;i<nv;i++) {
658 9328 : svd->sigma[i] = PetscRealPart(w[i]);
659 9328 : svd->sign[i] = omega[i];
660 : }
661 376 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&omega));
662 :
663 : /* check convergence */
664 376 : PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,1.0,&k));
665 376 : PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
666 :
667 : /* update l */
668 376 : if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
669 362 : else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
670 376 : if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
671 376 : if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
672 :
673 376 : if (svd->reason == SVD_CONVERGED_ITERATING) {
674 362 : if (PetscUnlikely(breakdown || k==nv)) {
675 : /* Start a new bidiagonalization */
676 0 : PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
677 0 : if (k<svd->nsv) {
678 0 : PetscCall(BVSetRandomColumn(V,k));
679 0 : PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,&breakdown));
680 0 : if (breakdown) {
681 0 : svd->reason = SVD_DIVERGED_BREAKDOWN;
682 0 : PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
683 : }
684 : }
685 362 : } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
686 : }
687 :
688 : /* compute converged singular vectors and restart vectors */
689 376 : PetscCall(DSGetMat(svd->ds,DS_MAT_V,&VV));
690 376 : PetscCall(BVMultInPlace(V,VV,svd->nconv,k+l));
691 376 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&VV));
692 376 : PetscCall(DSGetMat(svd->ds,DS_MAT_U,&UU));
693 376 : PetscCall(BVMultInPlace(U,UU,svd->nconv,k+l));
694 376 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&UU));
695 :
696 : /* copy the last vector of V to be the next initial vector
697 : and change signature matrix of U */
698 376 : if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) {
699 362 : PetscCall(BVCopyColumn(V,nv,k+l));
700 362 : PetscCall(BVSetActiveColumns(U,0,k+l));
701 362 : PetscCall(DSGetMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
702 362 : PetscCall(BVSetSignature(U,vomega));
703 362 : PetscCall(DSRestoreMatAndColumn(svd->ds,DS_MAT_D,0,&D,&vomega));
704 : }
705 :
706 376 : svd->nconv = k;
707 390 : PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,nv));
708 : }
709 :
710 : /* free working space */
711 14 : PetscCall(PetscFree(w));
712 14 : PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
713 14 : PetscFunctionReturn(PETSC_SUCCESS);
714 : }
715 :
716 : /* Given n computed generalized singular values in sigmain, backtransform them
717 : in sigmaout by undoing scaling and reciprocating if swapped=true. Also updates vectors V
718 : if given. If sigmaout=NULL then the result overwrites sigmain. */
719 879 : static PetscErrorCode SVDLanczosBackTransform(SVD svd,PetscInt n,PetscReal *sigmain,PetscReal *sigmaout,BV V)
720 : {
721 879 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
722 879 : PetscInt i;
723 879 : PetscReal c,s,r,f,scalef;
724 :
725 879 : PetscFunctionBegin;
726 879 : scalef = svd->swapped? 1.0/lanczos->scalef: lanczos->scalef;
727 8902 : for (i=0;i<n;i++) {
728 8023 : if (V && scalef != 1.0) {
729 60 : s = 1.0/PetscSqrtReal(1.0+sigmain[i]*sigmain[i]);
730 60 : c = sigmain[i]*s;
731 60 : r = s/scalef;
732 60 : f = 1.0/PetscSqrtReal(c*c+r*r);
733 60 : PetscCall(BVScaleColumn(V,i,f));
734 : }
735 8023 : if (sigmaout) {
736 7740 : if (svd->swapped) sigmaout[i] = 1.0/(sigmain[i]*scalef);
737 7100 : else sigmaout[i] = sigmain[i]*scalef;
738 : } else {
739 283 : sigmain[i] *= scalef;
740 283 : if (svd->swapped) sigmain[i] = 1.0/sigmain[i];
741 : }
742 : }
743 879 : PetscFunctionReturn(PETSC_SUCCESS);
744 : }
745 :
746 195 : static PetscErrorCode SVDLanczosGSingle(SVD svd,PetscReal *alpha,PetscReal *beta,Mat Z,BV V,BV U,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
747 : {
748 195 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
749 195 : PetscInt i,j,m;
750 195 : const PetscScalar *carr;
751 195 : PetscScalar *arr;
752 195 : Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
753 195 : PetscBool lindep=PETSC_FALSE;
754 :
755 195 : PetscFunctionBegin;
756 195 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
757 195 : PetscCall(BVGetColumn(V,k,&v));
758 195 : PetscCall(BVGetColumn(U,k,&u));
759 :
760 : /* Form ut=[u;0] */
761 195 : PetscCall(VecZeroEntries(ut));
762 195 : PetscCall(VecGetLocalSize(u,&m));
763 195 : PetscCall(VecGetArrayRead(u,&carr));
764 195 : PetscCall(VecGetArray(ut,&arr));
765 3459 : for (j=0; j<m; j++) arr[j] = carr[j];
766 195 : PetscCall(VecRestoreArrayRead(u,&carr));
767 195 : PetscCall(VecRestoreArray(ut,&arr));
768 :
769 : /* Solve least squares problem */
770 195 : PetscCall(KSPSolve(ksp,ut,x));
771 :
772 195 : PetscCall(MatMult(Z,x,v));
773 :
774 195 : PetscCall(BVRestoreColumn(U,k,&u));
775 195 : PetscCall(BVRestoreColumn(V,k,&v));
776 195 : PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,&lindep));
777 195 : if (PetscUnlikely(lindep)) {
778 0 : *n = k;
779 0 : if (breakdown) *breakdown = lindep;
780 0 : PetscFunctionReturn(PETSC_SUCCESS);
781 : }
782 :
783 1069 : for (i=k+1; i<*n; i++) {
784 :
785 : /* Compute vector i of BV U */
786 874 : PetscCall(BVGetColumn(V,i-1,&v));
787 874 : PetscCall(VecGetArray(v,&arr));
788 874 : PetscCall(VecPlaceArray(v1,arr));
789 874 : PetscCall(VecRestoreArray(v,&arr));
790 874 : PetscCall(BVRestoreColumn(V,i-1,&v));
791 874 : PetscCall(BVInsertVec(U,i,v1));
792 874 : PetscCall(VecResetArray(v1));
793 874 : PetscCall(BVOrthonormalizeColumn(U,i,PETSC_FALSE,beta+i-1,&lindep));
794 874 : if (PetscUnlikely(lindep)) {
795 0 : *n = i;
796 0 : break;
797 : }
798 :
799 : /* Compute vector i of BV V */
800 :
801 874 : PetscCall(BVGetColumn(V,i,&v));
802 874 : PetscCall(BVGetColumn(U,i,&u));
803 :
804 : /* Form ut=[u;0] */
805 874 : PetscCall(VecGetArrayRead(u,&carr));
806 874 : PetscCall(VecGetArray(ut,&arr));
807 15635 : for (j=0; j<m; j++) arr[j] = carr[j];
808 874 : PetscCall(VecRestoreArrayRead(u,&carr));
809 874 : PetscCall(VecRestoreArray(ut,&arr));
810 :
811 : /* Solve least squares problem */
812 874 : PetscCall(KSPSolve(ksp,ut,x));
813 :
814 874 : PetscCall(MatMult(Z,x,v));
815 :
816 874 : PetscCall(BVRestoreColumn(U,i,&u));
817 874 : PetscCall(BVRestoreColumn(V,i,&v));
818 874 : if (!lanczos->oneside || i==k+1) PetscCall(BVOrthonormalizeColumn(V,i,PETSC_FALSE,alpha+i,&lindep));
819 : else { /* cheap computation of V[i], if restart (i==k+1) do a full reorthogonalization */
820 267 : PetscCall(BVGetColumn(V,i,&u2));
821 267 : PetscCall(BVGetColumn(V,i-1,&u1));
822 267 : PetscCall(VecAXPY(u2,-beta[i-1],u1));
823 267 : PetscCall(BVRestoreColumn(V,i-1,&u1));
824 267 : PetscCall(VecNorm(u2,NORM_2,&alpha[i]));
825 267 : if (alpha[i]==0.0) lindep = PETSC_TRUE;
826 267 : else PetscCall(VecScale(u2,1.0/alpha[i]));
827 267 : PetscCall(BVRestoreColumn(V,i,&u2));
828 : }
829 874 : if (PetscUnlikely(lindep)) {
830 0 : *n = i;
831 0 : break;
832 : }
833 : }
834 :
835 : /* Compute vector n of BV U */
836 195 : if (!lindep) {
837 195 : PetscCall(BVGetColumn(V,*n-1,&v));
838 195 : PetscCall(VecGetArray(v,&arr));
839 195 : PetscCall(VecPlaceArray(v1,arr));
840 195 : PetscCall(VecRestoreArray(v,&arr));
841 195 : PetscCall(BVRestoreColumn(V,*n-1,&v));
842 195 : PetscCall(BVInsertVec(U,*n,v1));
843 195 : PetscCall(VecResetArray(v1));
844 195 : PetscCall(BVOrthonormalizeColumn(U,*n,PETSC_FALSE,beta+*n-1,&lindep));
845 : }
846 195 : if (breakdown) *breakdown = lindep;
847 195 : PetscCall(VecDestroy(&v1));
848 195 : PetscFunctionReturn(PETSC_SUCCESS);
849 : }
850 :
851 : /* solve generalized problem with single bidiagonalization of Q_A */
852 10 : static PetscErrorCode SVDSolve_TRLanczosGSingle(SVD svd,BV U1,BV V)
853 : {
854 10 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
855 10 : PetscReal *alpha,*beta,normr,scaleth,sigma0,*sigma;
856 10 : PetscScalar *w;
857 10 : PetscInt i,k,l,nv,ld;
858 10 : Mat U,VV;
859 10 : PetscBool breakdown=PETSC_FALSE;
860 :
861 10 : PetscFunctionBegin;
862 10 : PetscCall(DSGetLeadingDimension(svd->ds,&ld));
863 10 : PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
864 10 : normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*lanczos->scalef): 1.0;
865 : /* Convert scale threshold th=c/s to the corresponding c */
866 10 : scaleth = (lanczos->scaleth!=0)? lanczos->scaleth/PetscSqrtReal(lanczos->scaleth*lanczos->scaleth+1): 0.0;
867 :
868 : /* normalize start vector */
869 10 : if (!svd->ninil) {
870 9 : PetscCall(BVSetRandomColumn(U1,0));
871 9 : PetscCall(BVOrthonormalizeColumn(U1,0,PETSC_TRUE,NULL,NULL));
872 : }
873 :
874 : l = 0;
875 205 : while (svd->reason == SVD_CONVERGED_ITERATING) {
876 195 : svd->its++;
877 :
878 : /* inner loop */
879 195 : nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
880 195 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
881 195 : beta = alpha + ld;
882 195 : PetscCall(SVDLanczosGSingle(svd,alpha,beta,lanczos->Z,V,U1,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
883 195 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
884 195 : PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
885 195 : PetscCall(BVSetActiveColumns(U1,svd->nconv,nv));
886 :
887 : /* solve projected problem */
888 195 : PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
889 195 : PetscCall(DSSVDSetDimensions(svd->ds,nv));
890 195 : PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
891 195 : PetscCall(DSSolve(svd->ds,w,NULL));
892 195 : PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
893 195 : PetscCall(DSUpdateExtraRow(svd->ds));
894 195 : PetscCall(DSSynchronize(svd->ds,w,NULL));
895 1955 : for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
896 :
897 : /* check convergence */
898 195 : PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
899 195 : PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
900 :
901 195 : sigma0 = svd->which==SVD_LARGEST? svd->sigma[0] : 1.0/svd->sigma[0];
902 195 : if (scaleth!=0 && k==0 && sigma0>scaleth) {
903 :
904 : /* Scale and start from scratch */
905 0 : lanczos->scalef *= svd->sigma[0]/PetscSqrtReal(1-svd->sigma[0]*svd->sigma[0]);
906 0 : PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
907 0 : PetscCall(MatZUpdateScale(svd));
908 0 : if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*lanczos->scalef);
909 : l = 0;
910 :
911 : } else {
912 :
913 : /* update l */
914 195 : if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
915 185 : else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
916 195 : if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
917 195 : if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
918 :
919 195 : if (svd->reason == SVD_CONVERGED_ITERATING) {
920 185 : if (PetscUnlikely(breakdown || k==nv)) {
921 : /* Start a new bidiagonalization */
922 0 : PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
923 0 : if (k<svd->nsv) {
924 0 : PetscCall(BVSetRandomColumn(U1,k));
925 0 : PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,&breakdown));
926 0 : if (breakdown) {
927 0 : svd->reason = SVD_DIVERGED_BREAKDOWN;
928 0 : PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
929 : }
930 : }
931 185 : } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
932 : }
933 :
934 : /* compute converged singular vectors and restart vectors */
935 195 : PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
936 195 : PetscCall(BVMultInPlace(V,U,svd->nconv,k+l));
937 195 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
938 195 : PetscCall(DSGetMat(svd->ds,DS_MAT_V,&VV));
939 195 : PetscCall(BVMultInPlace(U1,VV,svd->nconv,k+l));
940 195 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&VV));
941 :
942 : /* copy the last vector to be the next initial vector */
943 195 : if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(U1,nv,k+l));
944 : }
945 :
946 195 : svd->nconv = k;
947 195 : PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
948 205 : PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
949 : }
950 :
951 10 : PetscCall(PetscFree2(w,sigma));
952 10 : PetscFunctionReturn(PETSC_SUCCESS);
953 : }
954 :
955 : /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (single variant) */
956 10 : static inline PetscErrorCode SVDLeftSingularVectors_Single(SVD svd,BV U1,BV U2)
957 : {
958 10 : PetscInt i,k,m,p;
959 10 : Vec u,u1,u2;
960 10 : PetscScalar *ua,*u2a;
961 10 : const PetscScalar *u1a;
962 10 : PetscReal s;
963 :
964 10 : PetscFunctionBegin;
965 10 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
966 10 : PetscCall(MatGetLocalSize(svd->B,&p,NULL));
967 51 : for (i=0;i<svd->nconv;i++) {
968 41 : PetscCall(BVGetColumn(U1,i,&u1));
969 41 : PetscCall(BVGetColumn(U2,i,&u2));
970 41 : PetscCall(BVGetColumn(svd->U,i,&u));
971 41 : PetscCall(VecGetArrayRead(u1,&u1a));
972 41 : PetscCall(VecGetArray(u,&ua));
973 41 : PetscCall(VecGetArray(u2,&u2a));
974 : /* Copy column from U1 to upper part of u */
975 717 : for (k=0;k<m;k++) ua[k] = u1a[k];
976 : /* Copy column from lower part of U to U2. Orthogonalize column in U2 and copy back to U */
977 905 : for (k=0;k<p;k++) u2a[k] = ua[m+k];
978 41 : PetscCall(VecRestoreArray(u2,&u2a));
979 41 : PetscCall(BVRestoreColumn(U2,i,&u2));
980 41 : PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,&s,NULL));
981 41 : PetscCall(BVGetColumn(U2,i,&u2));
982 41 : PetscCall(VecGetArray(u2,&u2a));
983 905 : for (k=0;k<p;k++) ua[m+k] = u2a[k];
984 : /* Update singular value */
985 41 : svd->sigma[i] /= s;
986 41 : PetscCall(VecRestoreArrayRead(u1,&u1a));
987 41 : PetscCall(VecRestoreArray(u,&ua));
988 41 : PetscCall(VecRestoreArray(u2,&u2a));
989 41 : PetscCall(BVRestoreColumn(U1,i,&u1));
990 41 : PetscCall(BVRestoreColumn(U2,i,&u2));
991 41 : PetscCall(BVRestoreColumn(svd->U,i,&u));
992 : }
993 10 : PetscFunctionReturn(PETSC_SUCCESS);
994 : }
995 :
996 247 : static PetscErrorCode SVDLanczosGUpper(SVD svd,PetscReal *alpha,PetscReal *beta,PetscReal *alphah,PetscReal *betah,Mat Z,BV U1,BV U2,BV V,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
997 : {
998 247 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
999 247 : PetscInt i,j,m,p;
1000 247 : const PetscScalar *carr;
1001 247 : PetscScalar *arr,*u2arr;
1002 247 : Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
1003 247 : PetscBool lindep=PETSC_FALSE,lindep1=PETSC_FALSE,lindep2=PETSC_FALSE;
1004 :
1005 247 : PetscFunctionBegin;
1006 247 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
1007 247 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1008 247 : PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1009 :
1010 1572 : for (i=k; i<*n; i++) {
1011 : /* Compute vector i of BV U1 */
1012 1325 : PetscCall(BVGetColumn(V,i,&v));
1013 1325 : PetscCall(VecGetArrayRead(v,&carr));
1014 1325 : PetscCall(VecPlaceArray(v1,carr));
1015 1325 : PetscCall(BVInsertVec(U1,i,v1));
1016 1325 : PetscCall(VecResetArray(v1));
1017 1325 : if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U1,i,PETSC_FALSE,alpha+i,&lindep1));
1018 : else { /* cheap computation of U1[i], if restart (i==k) do a full reorthogonalization */
1019 484 : PetscCall(BVGetColumn(U1,i,&u2));
1020 484 : if (i>0) {
1021 484 : PetscCall(BVGetColumn(U1,i-1,&u1));
1022 484 : PetscCall(VecAXPY(u2,-beta[i-1],u1));
1023 484 : PetscCall(BVRestoreColumn(U1,i-1,&u1));
1024 : }
1025 484 : PetscCall(VecNorm(u2,NORM_2,&alpha[i]));
1026 484 : if (alpha[i]==0.0) lindep = PETSC_TRUE;
1027 484 : else PetscCall(VecScale(u2,1.0/alpha[i]));
1028 484 : PetscCall(BVRestoreColumn(U1,i,&u2));
1029 : }
1030 :
1031 : /* Compute vector i of BV U2 */
1032 1325 : PetscCall(BVGetColumn(U2,i,&u2));
1033 1325 : PetscCall(VecGetArray(u2,&u2arr));
1034 1325 : if (i%2) {
1035 23996 : for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
1036 : } else {
1037 25607 : for (j=0; j<p; j++) u2arr[j] = carr[m+j];
1038 : }
1039 1325 : PetscCall(VecRestoreArray(u2,&u2arr));
1040 1325 : PetscCall(VecRestoreArrayRead(v,&carr));
1041 1325 : PetscCall(BVRestoreColumn(V,i,&v));
1042 1325 : if (lanczos->oneside && i>k) { /* cheap computation of U2[i], if restart (i==k) do a full reorthogonalization */
1043 484 : if (i>0) {
1044 484 : PetscCall(BVGetColumn(U2,i-1,&u1));
1045 484 : PetscCall(VecAXPY(u2,(i%2)?betah[i-1]:-betah[i-1],u1));
1046 484 : PetscCall(BVRestoreColumn(U2,i-1,&u1));
1047 : }
1048 484 : PetscCall(VecNorm(u2,NORM_2,&alphah[i]));
1049 484 : if (alphah[i]==0.0) lindep = PETSC_TRUE;
1050 484 : else PetscCall(VecScale(u2,1.0/alphah[i]));
1051 : }
1052 1325 : PetscCall(BVRestoreColumn(U2,i,&u2));
1053 1325 : if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep2));
1054 1325 : if (i%2) alphah[i] = -alphah[i];
1055 1325 : if (PetscUnlikely(lindep1 || lindep2)) {
1056 0 : lindep = PETSC_TRUE;
1057 0 : *n = i;
1058 0 : break;
1059 : }
1060 :
1061 : /* Compute vector i+1 of BV V */
1062 1325 : PetscCall(BVGetColumn(V,i+1,&v));
1063 : /* Form ut=[u;0] */
1064 1325 : PetscCall(BVGetColumn(U1,i,&u));
1065 1325 : PetscCall(VecZeroEntries(ut));
1066 1325 : PetscCall(VecGetArrayRead(u,&carr));
1067 1325 : PetscCall(VecGetArray(ut,&arr));
1068 37307 : for (j=0; j<m; j++) arr[j] = carr[j];
1069 1325 : PetscCall(VecRestoreArrayRead(u,&carr));
1070 1325 : PetscCall(VecRestoreArray(ut,&arr));
1071 : /* Solve least squares problem */
1072 1325 : PetscCall(KSPSolve(ksp,ut,x));
1073 1325 : PetscCall(MatMult(Z,x,v));
1074 1325 : PetscCall(BVRestoreColumn(U1,i,&u));
1075 1325 : PetscCall(BVRestoreColumn(V,i+1,&v));
1076 1325 : PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,beta+i,&lindep));
1077 1325 : betah[i] = -alpha[i]*beta[i]/alphah[i];
1078 1325 : if (PetscUnlikely(lindep)) {
1079 0 : *n = i;
1080 0 : break;
1081 : }
1082 : }
1083 247 : if (breakdown) *breakdown = lindep;
1084 247 : PetscCall(VecDestroy(&v1));
1085 247 : PetscFunctionReturn(PETSC_SUCCESS);
1086 : }
1087 :
1088 : /* generate random initial vector in column k for joint upper-upper bidiagonalization */
1089 21 : static inline PetscErrorCode SVDInitialVectorGUpper(SVD svd,BV V,BV U1,PetscInt k,PetscBool *breakdown)
1090 : {
1091 21 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1092 21 : Vec u,v,ut=svd->workl[0],x=svd->workr[0];
1093 21 : PetscInt m,j;
1094 21 : PetscScalar *arr;
1095 21 : const PetscScalar *carr;
1096 :
1097 21 : PetscFunctionBegin;
1098 : /* Form ut=[u;0] where u is the k-th column of U1 */
1099 21 : PetscCall(VecZeroEntries(ut));
1100 21 : PetscCall(BVGetColumn(U1,k,&u));
1101 21 : PetscCall(VecGetLocalSize(u,&m));
1102 21 : PetscCall(VecGetArrayRead(u,&carr));
1103 21 : PetscCall(VecGetArray(ut,&arr));
1104 1111 : for (j=0; j<m; j++) arr[j] = carr[j];
1105 21 : PetscCall(VecRestoreArrayRead(u,&carr));
1106 21 : PetscCall(VecRestoreArray(ut,&arr));
1107 21 : PetscCall(BVRestoreColumn(U1,k,&u));
1108 : /* Solve least squares problem Z*x=ut for x. Then set v=Z*x */
1109 21 : PetscCall(KSPSolve(lanczos->ksp,ut,x));
1110 21 : PetscCall(BVGetColumn(V,k,&v));
1111 21 : PetscCall(MatMult(lanczos->Z,x,v));
1112 21 : PetscCall(BVRestoreColumn(V,k,&v));
1113 21 : if (breakdown) PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,NULL,breakdown));
1114 21 : else PetscCall(BVOrthonormalizeColumn(V,k,PETSC_TRUE,NULL,NULL));
1115 21 : PetscFunctionReturn(PETSC_SUCCESS);
1116 : }
1117 :
1118 : /* solve generalized problem with joint upper-upper bidiagonalization */
1119 16 : static PetscErrorCode SVDSolve_TRLanczosGUpper(SVD svd,BV U1,BV U2,BV V)
1120 : {
1121 16 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1122 16 : PetscReal *alpha,*beta,*alphah,*betah,normr,sigma0,*sigma;
1123 16 : PetscScalar *w;
1124 16 : PetscInt i,k,l,nv,ld;
1125 16 : Mat U,Vmat,X;
1126 16 : PetscBool breakdown=PETSC_FALSE;
1127 :
1128 16 : PetscFunctionBegin;
1129 16 : PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1130 16 : PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
1131 16 : normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*lanczos->scalef): 1.0;
1132 :
1133 : /* normalize start vector */
1134 16 : if (!svd->ninil) PetscCall(BVSetRandomColumn(U1,0));
1135 16 : PetscCall(SVDInitialVectorGUpper(svd,V,U1,0,NULL));
1136 :
1137 : l = 0;
1138 263 : while (svd->reason == SVD_CONVERGED_ITERATING) {
1139 247 : svd->its++;
1140 :
1141 : /* inner loop */
1142 247 : nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1143 247 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1144 247 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&alphah));
1145 247 : beta = alpha + ld;
1146 247 : betah = alpha + 2*ld;
1147 247 : PetscCall(SVDLanczosGUpper(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
1148 247 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1149 247 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah));
1150 247 : PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
1151 247 : PetscCall(BVSetActiveColumns(U1,svd->nconv,nv));
1152 247 : PetscCall(BVSetActiveColumns(U2,svd->nconv,nv));
1153 :
1154 : /* solve projected problem */
1155 247 : PetscCall(DSSetDimensions(svd->ds,nv,svd->nconv,svd->nconv+l));
1156 247 : PetscCall(DSGSVDSetDimensions(svd->ds,nv,nv));
1157 247 : PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
1158 247 : PetscCall(DSSolve(svd->ds,w,NULL));
1159 247 : PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
1160 247 : PetscCall(DSUpdateExtraRow(svd->ds));
1161 247 : PetscCall(DSSynchronize(svd->ds,w,NULL));
1162 2390 : for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
1163 :
1164 : /* check convergence */
1165 247 : PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
1166 247 : PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
1167 :
1168 247 : sigma0 = svd->which==SVD_LARGEST? svd->sigma[0] : 1.0/svd->sigma[0];
1169 247 : if (lanczos->scaleth!=0 && k==0 && sigma0>lanczos->scaleth) {
1170 :
1171 : /* Scale and start from scratch */
1172 5 : lanczos->scalef *= svd->sigma[0];
1173 5 : PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
1174 5 : PetscCall(MatZUpdateScale(svd));
1175 5 : if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*lanczos->scalef);
1176 5 : l = 0;
1177 5 : if (!svd->ninil) PetscCall(BVSetRandomColumn(U1,0));
1178 5 : PetscCall(SVDInitialVectorGUpper(svd,V,U1,0,NULL));
1179 :
1180 : } else {
1181 :
1182 : /* update l */
1183 242 : if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1184 226 : else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1185 242 : if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1186 242 : if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1187 :
1188 242 : if (svd->reason == SVD_CONVERGED_ITERATING) {
1189 226 : if (PetscUnlikely(breakdown || k==nv)) {
1190 : /* Start a new bidiagonalization */
1191 0 : PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
1192 0 : if (k<svd->nsv) {
1193 0 : PetscCall(BVSetRandomColumn(U1,k));
1194 0 : PetscCall(SVDInitialVectorGUpper(svd,V,U1,k,&breakdown));
1195 0 : if (breakdown) {
1196 0 : svd->reason = SVD_DIVERGED_BREAKDOWN;
1197 0 : PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
1198 : }
1199 : }
1200 226 : } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
1201 : }
1202 : /* compute converged singular vectors and restart vectors */
1203 242 : PetscCall(DSGetMat(svd->ds,DS_MAT_X,&X));
1204 242 : PetscCall(BVMultInPlace(V,X,svd->nconv,k+l));
1205 242 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_X,&X));
1206 242 : PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
1207 242 : PetscCall(BVMultInPlace(U1,U,svd->nconv,k+l));
1208 242 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
1209 242 : PetscCall(DSGetMat(svd->ds,DS_MAT_V,&Vmat));
1210 242 : PetscCall(BVMultInPlace(U2,Vmat,svd->nconv,k+l));
1211 242 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&Vmat));
1212 :
1213 : /* copy the last vector to be the next initial vector */
1214 242 : if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(V,nv,k+l));
1215 : }
1216 :
1217 247 : svd->nconv = k;
1218 247 : PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
1219 263 : PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
1220 : }
1221 :
1222 16 : PetscCall(PetscFree2(w,sigma));
1223 16 : PetscFunctionReturn(PETSC_SUCCESS);
1224 : }
1225 :
1226 : /* Move generalized left singular vectors (0..nconv) from U1 and U2 to its final destination svd->U (upper and lower variants) */
1227 54 : static inline PetscErrorCode SVDLeftSingularVectors(SVD svd,BV U1,BV U2)
1228 : {
1229 54 : PetscInt i,k,m,p;
1230 54 : Vec u,u1,u2;
1231 54 : PetscScalar *ua;
1232 54 : const PetscScalar *u1a,*u2a;
1233 :
1234 54 : PetscFunctionBegin;
1235 54 : PetscCall(BVGetSizes(U1,&m,NULL,NULL));
1236 54 : PetscCall(BVGetSizes(U2,&p,NULL,NULL));
1237 296 : for (i=0;i<svd->nconv;i++) {
1238 242 : PetscCall(BVGetColumn(U1,i,&u1));
1239 242 : PetscCall(BVGetColumn(U2,i,&u2));
1240 242 : PetscCall(BVGetColumn(svd->U,i,&u));
1241 242 : PetscCall(VecGetArrayRead(u1,&u1a));
1242 242 : PetscCall(VecGetArrayRead(u2,&u2a));
1243 242 : PetscCall(VecGetArray(u,&ua));
1244 : /* Copy column from u1 to upper part of u */
1245 12744 : for (k=0;k<m;k++) ua[k] = u1a[k];
1246 : /* Copy column from u2 to lower part of u */
1247 17996 : for (k=0;k<p;k++) ua[m+k] = u2a[k];
1248 242 : PetscCall(VecRestoreArrayRead(u1,&u1a));
1249 242 : PetscCall(VecRestoreArrayRead(u2,&u2a));
1250 242 : PetscCall(VecRestoreArray(u,&ua));
1251 242 : PetscCall(BVRestoreColumn(U1,i,&u1));
1252 242 : PetscCall(BVRestoreColumn(U2,i,&u2));
1253 242 : PetscCall(BVRestoreColumn(svd->U,i,&u));
1254 : }
1255 54 : PetscFunctionReturn(PETSC_SUCCESS);
1256 : }
1257 :
1258 373 : static PetscErrorCode SVDLanczosGLower(SVD svd,PetscReal *alpha,PetscReal *beta,PetscReal *alphah,PetscReal *betah,Mat Z,BV U1,BV U2,BV V,KSP ksp,PetscInt k,PetscInt *n,PetscBool *breakdown)
1259 : {
1260 373 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1261 373 : PetscInt i,j,m,p;
1262 373 : const PetscScalar *carr;
1263 373 : PetscScalar *arr,*u2arr;
1264 373 : Vec u,v,ut=svd->workl[0],x=svd->workr[0],v1,u1,u2;
1265 373 : PetscBool lindep=PETSC_FALSE;
1266 :
1267 373 : PetscFunctionBegin;
1268 373 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&v1));
1269 373 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1270 373 : PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1271 :
1272 2585 : for (i=k; i<*n; i++) {
1273 : /* Compute vector i of BV U2 */
1274 2212 : PetscCall(BVGetColumn(V,i,&v));
1275 2212 : PetscCall(VecGetArrayRead(v,&carr));
1276 2212 : PetscCall(BVGetColumn(U2,i,&u2));
1277 2212 : PetscCall(VecGetArray(u2,&u2arr));
1278 2212 : if (i%2) {
1279 54451 : for (j=0; j<p; j++) u2arr[j] = -carr[m+j];
1280 : } else {
1281 52515 : for (j=0; j<p; j++) u2arr[j] = carr[m+j];
1282 : }
1283 2212 : PetscCall(VecRestoreArray(u2,&u2arr));
1284 2212 : if (lanczos->oneside && i>k) { /* cheap computation of U2[i], if restart (i==k) do a full reorthogonalization */
1285 450 : if (i>0) {
1286 450 : PetscCall(BVGetColumn(U2,i-1,&u1));
1287 450 : PetscCall(VecAXPY(u2,(i%2)?betah[i-1]:-betah[i-1],u1));
1288 450 : PetscCall(BVRestoreColumn(U2,i-1,&u1));
1289 : }
1290 450 : PetscCall(VecNorm(u2,NORM_2,&alphah[i]));
1291 450 : if (alphah[i]==0.0) lindep = PETSC_TRUE;
1292 450 : else PetscCall(VecScale(u2,1.0/alphah[i]));
1293 : }
1294 2212 : PetscCall(BVRestoreColumn(U2,i,&u2));
1295 2212 : if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(U2,i,PETSC_FALSE,alphah+i,&lindep));
1296 2212 : if (i%2) alphah[i] = -alphah[i];
1297 2212 : if (PetscUnlikely(lindep)) {
1298 0 : PetscCall(BVRestoreColumn(V,i,&v));
1299 0 : *n = i;
1300 0 : break;
1301 : }
1302 :
1303 : /* Compute vector i+1 of BV U1 */
1304 2212 : PetscCall(VecPlaceArray(v1,carr));
1305 2212 : PetscCall(BVInsertVec(U1,i+1,v1));
1306 2212 : PetscCall(VecResetArray(v1));
1307 2212 : PetscCall(BVOrthonormalizeColumn(U1,i+1,PETSC_FALSE,beta+i,&lindep));
1308 2212 : PetscCall(VecRestoreArrayRead(v,&carr));
1309 2212 : PetscCall(BVRestoreColumn(V,i,&v));
1310 2212 : if (PetscUnlikely(lindep)) {
1311 0 : *n = i+1;
1312 0 : break;
1313 : }
1314 :
1315 : /* Compute vector i+1 of BV V */
1316 2212 : PetscCall(BVGetColumn(V,i+1,&v));
1317 : /* Form ut=[u;0] where u is column i+1 of BV U1 */
1318 2212 : PetscCall(BVGetColumn(U1,i+1,&u));
1319 2212 : PetscCall(VecZeroEntries(ut));
1320 2212 : PetscCall(VecGetArrayRead(u,&carr));
1321 2212 : PetscCall(VecGetArray(ut,&arr));
1322 84420 : for (j=0; j<m; j++) arr[j] = carr[j];
1323 2212 : PetscCall(VecRestoreArrayRead(u,&carr));
1324 2212 : PetscCall(VecRestoreArray(ut,&arr));
1325 : /* Solve least squares problem */
1326 2212 : PetscCall(KSPSolve(ksp,ut,x));
1327 2212 : PetscCall(MatMult(Z,x,v));
1328 2212 : PetscCall(BVRestoreColumn(U1,i+1,&u));
1329 2212 : PetscCall(BVRestoreColumn(V,i+1,&v));
1330 2212 : if (!lanczos->oneside || i==k) PetscCall(BVOrthonormalizeColumn(V,i+1,PETSC_FALSE,alpha+i+1,&lindep));
1331 : else { /* cheap computation of V[i+1], if restart (i==k) do a full reorthogonalization */
1332 450 : PetscCall(BVGetColumn(V,i+1,&u2));
1333 450 : PetscCall(BVGetColumn(V,i,&u1));
1334 450 : PetscCall(VecAXPY(u2,-beta[i],u1));
1335 450 : PetscCall(BVRestoreColumn(V,i,&u1));
1336 450 : PetscCall(VecNorm(u2,NORM_2,&alpha[i+1]));
1337 450 : if (alpha[i+1]==0.0) lindep = PETSC_TRUE;
1338 450 : else PetscCall(VecScale(u2,1.0/alpha[i+1]));
1339 450 : PetscCall(BVRestoreColumn(V,i+1,&u2));
1340 : }
1341 2212 : betah[i] = -alpha[i+1]*beta[i]/alphah[i];
1342 2212 : if (PetscUnlikely(lindep)) {
1343 0 : *n = i+1;
1344 0 : break;
1345 : }
1346 : }
1347 373 : if (breakdown) *breakdown = lindep;
1348 373 : PetscCall(VecDestroy(&v1));
1349 373 : PetscFunctionReturn(PETSC_SUCCESS);
1350 : }
1351 :
1352 : /* generate random initial vector in column k for joint lower-upper bidiagonalization */
1353 47 : static inline PetscErrorCode SVDInitialVectorGLower(SVD svd,BV V,BV U1,BV U2,PetscInt k,PetscBool *breakdown)
1354 : {
1355 47 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1356 47 : const PetscScalar *carr;
1357 47 : PetscScalar *arr;
1358 47 : PetscReal *alpha;
1359 47 : PetscInt j,m,p;
1360 47 : Vec u,uh,v,ut=svd->workl[0],x=svd->workr[0];
1361 :
1362 47 : PetscFunctionBegin;
1363 47 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1364 47 : PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1365 : /* Form ut=[0;uh], where uh is the k-th column of U2 */
1366 47 : PetscCall(BVGetColumn(U2,k,&uh));
1367 47 : PetscCall(VecZeroEntries(ut));
1368 47 : PetscCall(VecGetArrayRead(uh,&carr));
1369 47 : PetscCall(VecGetArray(ut,&arr));
1370 3498 : for (j=0; j<p; j++) arr[m+j] = carr[j];
1371 47 : PetscCall(VecRestoreArrayRead(uh,&carr));
1372 47 : PetscCall(VecRestoreArray(ut,&arr));
1373 47 : PetscCall(BVRestoreColumn(U2,k,&uh));
1374 : /* Solve least squares problem Z*x=ut for x. Then set ut=Z*x */
1375 47 : PetscCall(KSPSolve(lanczos->ksp,ut,x));
1376 47 : PetscCall(MatMult(lanczos->Z,x,ut));
1377 : /* Form u, column k of BV U1, as the upper part of ut and orthonormalize */
1378 47 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&u));
1379 47 : PetscCall(VecGetArrayRead(ut,&carr));
1380 47 : PetscCall(VecPlaceArray(u,carr));
1381 47 : PetscCall(BVInsertVec(U1,k,u));
1382 47 : PetscCall(VecResetArray(u));
1383 47 : PetscCall(VecRestoreArrayRead(ut,&carr));
1384 47 : PetscCall(VecDestroy(&u));
1385 47 : if (breakdown) PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_FALSE,NULL,breakdown));
1386 47 : else PetscCall(BVOrthonormalizeColumn(U1,k,PETSC_TRUE,NULL,NULL));
1387 :
1388 47 : if (!breakdown || !*breakdown) {
1389 47 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1390 : /* Compute k-th vector of BV V */
1391 47 : PetscCall(BVGetColumn(V,k,&v));
1392 : /* Form ut=[u;0] where u is the 1st column of U1 */
1393 47 : PetscCall(BVGetColumn(U1,k,&u));
1394 47 : PetscCall(VecZeroEntries(ut));
1395 47 : PetscCall(VecGetArrayRead(u,&carr));
1396 47 : PetscCall(VecGetArray(ut,&arr));
1397 2408 : for (j=0; j<m; j++) arr[j] = carr[j];
1398 47 : PetscCall(VecRestoreArrayRead(u,&carr));
1399 47 : PetscCall(VecRestoreArray(ut,&arr));
1400 : /* Solve least squares problem */
1401 47 : PetscCall(KSPSolve(lanczos->ksp,ut,x));
1402 47 : PetscCall(MatMult(lanczos->Z,x,v));
1403 47 : PetscCall(BVRestoreColumn(U1,k,&u));
1404 47 : PetscCall(BVRestoreColumn(V,k,&v));
1405 47 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1406 47 : if (breakdown) PetscCall(BVOrthonormalizeColumn(V,k,PETSC_FALSE,alpha+k,breakdown));
1407 47 : else PetscCall(BVOrthonormalizeColumn(V,k,PETSC_TRUE,alpha+k,NULL));
1408 47 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1409 : }
1410 47 : PetscFunctionReturn(PETSC_SUCCESS);
1411 : }
1412 :
1413 : /* solve generalized problem with joint lower-upper bidiagonalization */
1414 38 : static PetscErrorCode SVDSolve_TRLanczosGLower(SVD svd,BV U1,BV U2,BV V)
1415 : {
1416 38 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1417 38 : PetscReal *alpha,*beta,*alphah,*betah,normr,scalef,*sigma,sigma0;
1418 38 : PetscScalar *w;
1419 38 : PetscInt i,k,l,nv,ld;
1420 38 : Mat U,Vmat,X;
1421 38 : PetscBool breakdown=PETSC_FALSE,inverted;
1422 :
1423 38 : PetscFunctionBegin;
1424 38 : PetscCall(DSGetLeadingDimension(svd->ds,&ld));
1425 38 : PetscCall(PetscMalloc2(ld,&w,ld,&sigma));
1426 38 : inverted = ((svd->which==SVD_LARGEST && svd->swapped) || (svd->which==SVD_SMALLEST && !svd->swapped))? PETSC_TRUE: PETSC_FALSE;
1427 38 : scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
1428 38 : normr = (svd->conv==SVD_CONV_ABS)? PetscMax(svd->nrma,svd->nrmb*scalef): 1.0;
1429 :
1430 : /* normalize start vector */
1431 38 : if (!svd->ninil) PetscCall(BVSetRandomColumn(U2,0));
1432 38 : PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,0,NULL));
1433 :
1434 : l = 0;
1435 411 : while (svd->reason == SVD_CONVERGED_ITERATING) {
1436 373 : svd->its++;
1437 :
1438 : /* inner loop */
1439 373 : nv = PetscMin(svd->nconv+svd->mpd,svd->ncv);
1440 373 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_T,&alpha));
1441 373 : PetscCall(DSGetArrayReal(svd->ds,DS_MAT_D,&alphah));
1442 373 : beta = alpha + ld;
1443 373 : betah = alpha + 2*ld;
1444 373 : PetscCall(SVDLanczosGLower(svd,alpha,beta,alphah,betah,lanczos->Z,U1,U2,V,lanczos->ksp,svd->nconv+l,&nv,&breakdown));
1445 373 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_T,&alpha));
1446 373 : PetscCall(DSRestoreArrayReal(svd->ds,DS_MAT_D,&alphah));
1447 373 : PetscCall(BVSetActiveColumns(V,svd->nconv,nv));
1448 373 : PetscCall(BVSetActiveColumns(U1,svd->nconv,nv+1));
1449 373 : PetscCall(BVSetActiveColumns(U2,svd->nconv,nv));
1450 :
1451 : /* solve projected problem */
1452 373 : PetscCall(DSSetDimensions(svd->ds,nv+1,svd->nconv,svd->nconv+l));
1453 373 : PetscCall(DSGSVDSetDimensions(svd->ds,nv,nv));
1454 373 : PetscCall(DSSetState(svd->ds,l?DS_STATE_RAW:DS_STATE_INTERMEDIATE));
1455 373 : PetscCall(DSSolve(svd->ds,w,NULL));
1456 373 : PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
1457 373 : PetscCall(DSUpdateExtraRow(svd->ds));
1458 373 : PetscCall(DSSynchronize(svd->ds,w,NULL));
1459 3941 : for (i=svd->nconv;i<nv;i++) svd->sigma[i] = PetscRealPart(w[i]);
1460 :
1461 : /* check convergence */
1462 373 : PetscCall(SVDKrylovConvergence(svd,PETSC_FALSE,svd->nconv,nv-svd->nconv,normr,&k));
1463 373 : PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,k,svd->nsv,&svd->reason,svd->stoppingctx));
1464 :
1465 373 : sigma0 = inverted? 1.0/svd->sigma[0] : svd->sigma[0];
1466 373 : if (lanczos->scaleth!=0 && k==0 && sigma0>lanczos->scaleth) {
1467 :
1468 : /* Scale and start from scratch */
1469 9 : lanczos->scalef *= svd->swapped? 1.0/svd->sigma[0] : svd->sigma[0];
1470 9 : PetscCall(PetscInfo(svd,"Scaling by factor %g and starting from scratch\n",(double)lanczos->scalef));
1471 9 : PetscCall(MatZUpdateScale(svd));
1472 9 : scalef = svd->swapped? 1.0/lanczos->scalef : lanczos->scalef;
1473 9 : if (svd->conv==SVD_CONV_ABS) normr = PetscMax(svd->nrma,svd->nrmb*scalef);
1474 9 : l = 0;
1475 9 : if (!svd->ninil) PetscCall(BVSetRandomColumn(U2,0));
1476 9 : PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,0,NULL));
1477 :
1478 : } else {
1479 :
1480 : /* update l */
1481 364 : if (svd->reason != SVD_CONVERGED_ITERATING || breakdown || k==nv) l = 0;
1482 326 : else l = PetscMax(1,(PetscInt)((nv-k)*lanczos->keep));
1483 364 : if (!lanczos->lock && l>0) { l += k; k = 0; } /* non-locking variant: reset no. of converged triplets */
1484 364 : if (l) PetscCall(PetscInfo(svd,"Preparing to restart keeping l=%" PetscInt_FMT " vectors\n",l));
1485 :
1486 364 : if (svd->reason == SVD_CONVERGED_ITERATING) {
1487 326 : if (PetscUnlikely(breakdown || k==nv)) {
1488 : /* Start a new bidiagonalization */
1489 0 : PetscCall(PetscInfo(svd,"Breakdown in bidiagonalization (it=%" PetscInt_FMT ")\n",svd->its));
1490 0 : if (k<svd->nsv) {
1491 0 : PetscCall(BVSetRandomColumn(U2,k));
1492 0 : PetscCall(SVDInitialVectorGLower(svd,V,U1,U2,k,&breakdown));
1493 0 : if (breakdown) {
1494 0 : svd->reason = SVD_DIVERGED_BREAKDOWN;
1495 0 : PetscCall(PetscInfo(svd,"Unable to generate more start vectors\n"));
1496 : }
1497 : }
1498 326 : } else PetscCall(DSTruncate(svd->ds,k+l,PETSC_FALSE));
1499 : }
1500 :
1501 : /* compute converged singular vectors and restart vectors */
1502 364 : PetscCall(DSGetMat(svd->ds,DS_MAT_X,&X));
1503 364 : PetscCall(BVMultInPlace(V,X,svd->nconv,k+l));
1504 364 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_X,&X));
1505 364 : PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
1506 364 : PetscCall(BVMultInPlace(U1,U,svd->nconv,k+l+1));
1507 364 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
1508 364 : PetscCall(DSGetMat(svd->ds,DS_MAT_V,&Vmat));
1509 364 : PetscCall(BVMultInPlace(U2,Vmat,svd->nconv,k+l));
1510 364 : PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&Vmat));
1511 :
1512 : /* copy the last vector to be the next initial vector */
1513 364 : if (svd->reason == SVD_CONVERGED_ITERATING && !breakdown) PetscCall(BVCopyColumn(V,nv,k+l));
1514 : }
1515 :
1516 373 : svd->nconv = k;
1517 373 : PetscCall(SVDLanczosBackTransform(svd,nv,svd->sigma,sigma,NULL));
1518 411 : PetscCall(SVDMonitor(svd,svd->its,svd->nconv,sigma,svd->errest,nv));
1519 : }
1520 :
1521 38 : PetscCall(PetscFree2(w,sigma));
1522 38 : PetscFunctionReturn(PETSC_SUCCESS);
1523 : }
1524 :
1525 64 : static PetscErrorCode SVDSolve_TRLanczos_GSVD(SVD svd)
1526 : {
1527 64 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1528 64 : PetscInt k,m,p;
1529 64 : PetscBool convchg=PETSC_FALSE;
1530 64 : BV U1,U2,UU;
1531 64 : BVType type;
1532 64 : VecType vtype;
1533 64 : Mat U,V;
1534 64 : SlepcSC sc;
1535 :
1536 64 : PetscFunctionBegin;
1537 64 : PetscCall(PetscCitationsRegister(citationg,&citedg));
1538 :
1539 64 : if (svd->swapped) {
1540 5 : PetscCall(DSGetSlepcSC(svd->ds,&sc));
1541 5 : if (svd->which==SVD_LARGEST) sc->comparison = SlepcCompareSmallestReal;
1542 0 : else sc->comparison = SlepcCompareLargestReal;
1543 : }
1544 64 : if (svd->converged==SVDConvergedNorm) { /* override temporarily since computed residual is already relative to the norms */
1545 64 : svd->converged = SVDConvergedAbsolute;
1546 64 : convchg = PETSC_TRUE;
1547 : }
1548 64 : PetscCall(MatGetLocalSize(svd->A,&m,NULL));
1549 64 : PetscCall(MatGetLocalSize(svd->B,&p,NULL));
1550 :
1551 : /* Create BV for U1 */
1552 64 : PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&U1));
1553 64 : PetscCall(BVGetType(svd->U,&type));
1554 64 : PetscCall(BVSetType(U1,type));
1555 64 : PetscCall(BVGetSizes(svd->U,NULL,NULL,&k));
1556 64 : PetscCall(BVSetSizes(U1,m,PETSC_DECIDE,k));
1557 64 : PetscCall(BVGetVecType(svd->U,&vtype));
1558 64 : PetscCall(BVSetVecType(U1,vtype));
1559 :
1560 : /* Create BV for U2 */
1561 64 : PetscCall(BVCreate(PetscObjectComm((PetscObject)svd),&U2));
1562 64 : PetscCall(BVSetType(U2,type));
1563 64 : PetscCall(BVSetSizes(U2,p,PETSC_DECIDE,k));
1564 64 : PetscCall(BVSetVecType(U2,vtype));
1565 :
1566 : /* Copy initial vectors from svd->U to U1 and U2 */
1567 64 : if (svd->ninil) {
1568 4 : Vec u, uh, nest, aux[2];
1569 4 : PetscCall(BVGetColumn(U1,0,&u));
1570 4 : PetscCall(BVGetColumn(U2,0,&uh));
1571 4 : aux[0] = u;
1572 4 : aux[1] = uh;
1573 4 : PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
1574 4 : PetscCall(BVCopyVec(svd->U,0,nest));
1575 4 : PetscCall(BVRestoreColumn(U1,0,&u));
1576 4 : PetscCall(BVRestoreColumn(U2,0,&uh));
1577 4 : PetscCall(VecDestroy(&nest));
1578 : }
1579 :
1580 64 : switch (lanczos->bidiag) {
1581 10 : case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1582 10 : PetscCall(SVDSolve_TRLanczosGSingle(svd,U1,svd->U));
1583 : break;
1584 16 : case SVD_TRLANCZOS_GBIDIAG_UPPER:
1585 16 : PetscCall(SVDSolve_TRLanczosGUpper(svd,U1,U2,svd->U));
1586 : break;
1587 38 : case SVD_TRLANCZOS_GBIDIAG_LOWER:
1588 38 : PetscCall(SVDSolve_TRLanczosGLower(svd,U1,U2,svd->U));
1589 : break;
1590 : }
1591 :
1592 : /* Compute converged right singular vectors */
1593 64 : PetscCall(BVSetActiveColumns(svd->U,0,svd->nconv));
1594 64 : PetscCall(BVSetActiveColumns(svd->V,0,svd->nconv));
1595 64 : PetscCall(BVGetMat(svd->U,&U));
1596 64 : PetscCall(BVGetMat(svd->V,&V));
1597 64 : PetscCall(KSPMatSolve(lanczos->ksp,U,V));
1598 64 : PetscCall(BVRestoreMat(svd->U,&U));
1599 64 : PetscCall(BVRestoreMat(svd->V,&V));
1600 :
1601 : /* Finish computing left singular vectors and move them to its place */
1602 64 : if (svd->swapped) SWAP(U1,U2,UU);
1603 64 : switch (lanczos->bidiag) {
1604 10 : case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1605 10 : PetscCall(SVDLeftSingularVectors_Single(svd,U1,U2));
1606 : break;
1607 54 : case SVD_TRLANCZOS_GBIDIAG_UPPER:
1608 : case SVD_TRLANCZOS_GBIDIAG_LOWER:
1609 54 : PetscCall(SVDLeftSingularVectors(svd,U1,U2));
1610 : break;
1611 : }
1612 :
1613 : /* undo scaling and compute the reciprocals of sigma if matrices were swapped */
1614 64 : PetscCall(SVDLanczosBackTransform(svd,svd->nconv,svd->sigma,NULL,svd->V));
1615 :
1616 64 : PetscCall(BVDestroy(&U1));
1617 64 : PetscCall(BVDestroy(&U2));
1618 64 : PetscCall(DSTruncate(svd->ds,svd->nconv,PETSC_TRUE));
1619 64 : if (convchg) svd->converged = SVDConvergedNorm;
1620 64 : PetscFunctionReturn(PETSC_SUCCESS);
1621 : }
1622 :
1623 95 : static PetscErrorCode SVDSetFromOptions_TRLanczos(SVD svd,PetscOptionItems *PetscOptionsObject)
1624 : {
1625 95 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1626 95 : PetscBool flg,val,lock;
1627 95 : PetscReal keep,scale;
1628 95 : SVDTRLanczosGBidiag bidiag;
1629 :
1630 95 : PetscFunctionBegin;
1631 95 : PetscOptionsHeadBegin(PetscOptionsObject,"SVD TRLanczos Options");
1632 :
1633 95 : PetscCall(PetscOptionsBool("-svd_trlanczos_oneside","Use one-side reorthogonalization","SVDTRLanczosSetOneSide",lanczos->oneside,&val,&flg));
1634 95 : if (flg) PetscCall(SVDTRLanczosSetOneSide(svd,val));
1635 :
1636 95 : PetscCall(PetscOptionsReal("-svd_trlanczos_restart","Proportion of vectors kept after restart","SVDTRLanczosSetRestart",0.5,&keep,&flg));
1637 95 : if (flg) PetscCall(SVDTRLanczosSetRestart(svd,keep));
1638 :
1639 95 : PetscCall(PetscOptionsBool("-svd_trlanczos_locking","Choose between locking and non-locking variants","SVDTRLanczosSetLocking",PETSC_TRUE,&lock,&flg));
1640 95 : if (flg) PetscCall(SVDTRLanczosSetLocking(svd,lock));
1641 :
1642 95 : PetscCall(PetscOptionsEnum("-svd_trlanczos_gbidiag","Bidiagonalization choice for Generalized Problem","SVDTRLanczosSetGBidiag",SVDTRLanczosGBidiags,(PetscEnum)lanczos->bidiag,(PetscEnum*)&bidiag,&flg));
1643 95 : if (flg) PetscCall(SVDTRLanczosSetGBidiag(svd,bidiag));
1644 :
1645 95 : PetscCall(PetscOptionsBool("-svd_trlanczos_explicitmatrix","Build explicit matrix for KSP solver","SVDTRLanczosSetExplicitMatrix",lanczos->explicitmatrix,&val,&flg));
1646 95 : if (flg) PetscCall(SVDTRLanczosSetExplicitMatrix(svd,val));
1647 :
1648 95 : PetscCall(SVDTRLanczosGetScale(svd,&scale));
1649 95 : PetscCall(PetscOptionsReal("-svd_trlanczos_scale","Scale parameter for matrix B","SVDTRLanczosSetScale",scale,&scale,&flg));
1650 95 : if (flg) PetscCall(SVDTRLanczosSetScale(svd,scale));
1651 :
1652 95 : PetscOptionsHeadEnd();
1653 :
1654 95 : if (svd->OPb) {
1655 58 : if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
1656 58 : PetscCall(KSPSetFromOptions(lanczos->ksp));
1657 : }
1658 95 : PetscFunctionReturn(PETSC_SUCCESS);
1659 : }
1660 :
1661 34 : static PetscErrorCode SVDTRLanczosSetOneSide_TRLanczos(SVD svd,PetscBool oneside)
1662 : {
1663 34 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1664 :
1665 34 : PetscFunctionBegin;
1666 34 : if (lanczos->oneside != oneside) {
1667 21 : lanczos->oneside = oneside;
1668 21 : svd->state = SVD_STATE_INITIAL;
1669 : }
1670 34 : PetscFunctionReturn(PETSC_SUCCESS);
1671 : }
1672 :
1673 : /*@
1674 : SVDTRLanczosSetOneSide - Indicate if the variant of the Lanczos method
1675 : to be used is one-sided or two-sided.
1676 :
1677 : Logically Collective
1678 :
1679 : Input Parameters:
1680 : + svd - singular value solver
1681 : - oneside - boolean flag indicating if the method is one-sided or not
1682 :
1683 : Options Database Key:
1684 : . -svd_trlanczos_oneside <boolean> - Indicates the boolean flag
1685 :
1686 : Notes:
1687 : By default, a two-sided variant is selected, which is sometimes slightly
1688 : more robust. However, the one-sided variant is faster because it avoids
1689 : the orthogonalization associated to left singular vectors.
1690 :
1691 : One-sided orthogonalization is also available for the GSVD, in which case
1692 : two orthogonalizations out of three are avoided.
1693 :
1694 : Level: advanced
1695 :
1696 : .seealso: SVDLanczosSetOneSide()
1697 : @*/
1698 34 : PetscErrorCode SVDTRLanczosSetOneSide(SVD svd,PetscBool oneside)
1699 : {
1700 34 : PetscFunctionBegin;
1701 34 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1702 136 : PetscValidLogicalCollectiveBool(svd,oneside,2);
1703 34 : PetscTryMethod(svd,"SVDTRLanczosSetOneSide_C",(SVD,PetscBool),(svd,oneside));
1704 34 : PetscFunctionReturn(PETSC_SUCCESS);
1705 : }
1706 :
1707 7 : static PetscErrorCode SVDTRLanczosGetOneSide_TRLanczos(SVD svd,PetscBool *oneside)
1708 : {
1709 7 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1710 :
1711 7 : PetscFunctionBegin;
1712 7 : *oneside = lanczos->oneside;
1713 7 : PetscFunctionReturn(PETSC_SUCCESS);
1714 : }
1715 :
1716 : /*@
1717 : SVDTRLanczosGetOneSide - Gets if the variant of the Lanczos method
1718 : to be used is one-sided or two-sided.
1719 :
1720 : Not Collective
1721 :
1722 : Input Parameters:
1723 : . svd - singular value solver
1724 :
1725 : Output Parameters:
1726 : . oneside - boolean flag indicating if the method is one-sided or not
1727 :
1728 : Level: advanced
1729 :
1730 : .seealso: SVDTRLanczosSetOneSide()
1731 : @*/
1732 7 : PetscErrorCode SVDTRLanczosGetOneSide(SVD svd,PetscBool *oneside)
1733 : {
1734 7 : PetscFunctionBegin;
1735 7 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1736 7 : PetscAssertPointer(oneside,2);
1737 7 : PetscUseMethod(svd,"SVDTRLanczosGetOneSide_C",(SVD,PetscBool*),(svd,oneside));
1738 7 : PetscFunctionReturn(PETSC_SUCCESS);
1739 : }
1740 :
1741 48 : static PetscErrorCode SVDTRLanczosSetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag bidiag)
1742 : {
1743 48 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1744 :
1745 48 : PetscFunctionBegin;
1746 48 : switch (bidiag) {
1747 48 : case SVD_TRLANCZOS_GBIDIAG_SINGLE:
1748 : case SVD_TRLANCZOS_GBIDIAG_UPPER:
1749 : case SVD_TRLANCZOS_GBIDIAG_LOWER:
1750 48 : if (lanczos->bidiag != bidiag) {
1751 32 : lanczos->bidiag = bidiag;
1752 32 : svd->state = SVD_STATE_INITIAL;
1753 : }
1754 48 : break;
1755 0 : default:
1756 0 : SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid bidiagonalization choice");
1757 : }
1758 48 : PetscFunctionReturn(PETSC_SUCCESS);
1759 : }
1760 :
1761 : /*@
1762 : SVDTRLanczosSetGBidiag - Sets the bidiagonalization choice to use in
1763 : the GSVD TRLanczos solver.
1764 :
1765 : Logically Collective
1766 :
1767 : Input Parameters:
1768 : + svd - the singular value solver
1769 : - bidiag - the bidiagonalization choice
1770 :
1771 : Options Database Key:
1772 : . -svd_trlanczos_gbidiag - Sets the bidiagonalization choice (either 's' or 'juu'
1773 : or 'jlu')
1774 :
1775 : Level: advanced
1776 :
1777 : .seealso: SVDTRLanczosGetGBidiag(), SVDTRLanczosGBidiag
1778 : @*/
1779 48 : PetscErrorCode SVDTRLanczosSetGBidiag(SVD svd,SVDTRLanczosGBidiag bidiag)
1780 : {
1781 48 : PetscFunctionBegin;
1782 48 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1783 192 : PetscValidLogicalCollectiveEnum(svd,bidiag,2);
1784 48 : PetscTryMethod(svd,"SVDTRLanczosSetGBidiag_C",(SVD,SVDTRLanczosGBidiag),(svd,bidiag));
1785 48 : PetscFunctionReturn(PETSC_SUCCESS);
1786 : }
1787 :
1788 6 : static PetscErrorCode SVDTRLanczosGetGBidiag_TRLanczos(SVD svd,SVDTRLanczosGBidiag *bidiag)
1789 : {
1790 6 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
1791 :
1792 6 : PetscFunctionBegin;
1793 6 : *bidiag = lanczos->bidiag;
1794 6 : PetscFunctionReturn(PETSC_SUCCESS);
1795 : }
1796 :
1797 : /*@
1798 : SVDTRLanczosGetGBidiag - Gets the bidiagonalization choice used in the GSVD
1799 : TRLanczos solver.
1800 :
1801 : Not Collective
1802 :
1803 : Input Parameter:
1804 : . svd - the singular value solver
1805 :
1806 : Output Parameter:
1807 : . bidiag - the bidiagonalization choice
1808 :
1809 : Level: advanced
1810 :
1811 : .seealso: SVDTRLanczosSetGBidiag(), SVDTRLanczosGBidiag
1812 : @*/
1813 6 : PetscErrorCode SVDTRLanczosGetGBidiag(SVD svd,SVDTRLanczosGBidiag *bidiag)
1814 : {
1815 6 : PetscFunctionBegin;
1816 6 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1817 6 : PetscAssertPointer(bidiag,2);
1818 6 : PetscUseMethod(svd,"SVDTRLanczosGetGBidiag_C",(SVD,SVDTRLanczosGBidiag*),(svd,bidiag));
1819 6 : PetscFunctionReturn(PETSC_SUCCESS);
1820 : }
1821 :
1822 6 : static PetscErrorCode SVDTRLanczosSetKSP_TRLanczos(SVD svd,KSP ksp)
1823 : {
1824 6 : SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1825 :
1826 6 : PetscFunctionBegin;
1827 6 : PetscCall(PetscObjectReference((PetscObject)ksp));
1828 6 : PetscCall(KSPDestroy(&ctx->ksp));
1829 6 : ctx->ksp = ksp;
1830 6 : svd->state = SVD_STATE_INITIAL;
1831 6 : PetscFunctionReturn(PETSC_SUCCESS);
1832 : }
1833 :
1834 : /*@
1835 : SVDTRLanczosSetKSP - Associate a linear solver object (KSP) to the SVD solver.
1836 :
1837 : Collective
1838 :
1839 : Input Parameters:
1840 : + svd - SVD solver
1841 : - ksp - the linear solver object
1842 :
1843 : Note:
1844 : Only used for the GSVD problem.
1845 :
1846 : Level: advanced
1847 :
1848 : .seealso: SVDTRLanczosGetKSP()
1849 : @*/
1850 6 : PetscErrorCode SVDTRLanczosSetKSP(SVD svd,KSP ksp)
1851 : {
1852 6 : PetscFunctionBegin;
1853 6 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1854 6 : PetscValidHeaderSpecific(ksp,KSP_CLASSID,2);
1855 6 : PetscCheckSameComm(svd,1,ksp,2);
1856 6 : PetscTryMethod(svd,"SVDTRLanczosSetKSP_C",(SVD,KSP),(svd,ksp));
1857 6 : PetscFunctionReturn(PETSC_SUCCESS);
1858 : }
1859 :
1860 52 : static PetscErrorCode SVDTRLanczosGetKSP_TRLanczos(SVD svd,KSP *ksp)
1861 : {
1862 52 : SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1863 52 : PC pc;
1864 :
1865 52 : PetscFunctionBegin;
1866 52 : if (!ctx->ksp) {
1867 : /* Create linear solver */
1868 52 : PetscCall(KSPCreate(PetscObjectComm((PetscObject)svd),&ctx->ksp));
1869 52 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)ctx->ksp,(PetscObject)svd,1));
1870 52 : PetscCall(KSPSetOptionsPrefix(ctx->ksp,((PetscObject)svd)->prefix));
1871 52 : PetscCall(KSPAppendOptionsPrefix(ctx->ksp,"svd_trlanczos_"));
1872 52 : PetscCall(PetscObjectSetOptions((PetscObject)ctx->ksp,((PetscObject)svd)->options));
1873 52 : PetscCall(KSPSetType(ctx->ksp,KSPLSQR));
1874 52 : PetscCall(KSPGetPC(ctx->ksp,&pc));
1875 52 : PetscCall(PCSetType(pc,PCNONE));
1876 52 : PetscCall(KSPSetErrorIfNotConverged(ctx->ksp,PETSC_TRUE));
1877 104 : PetscCall(KSPSetTolerances(ctx->ksp,SlepcDefaultTol(svd->tol)/10.0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT));
1878 : }
1879 52 : *ksp = ctx->ksp;
1880 52 : PetscFunctionReturn(PETSC_SUCCESS);
1881 : }
1882 :
1883 : /*@
1884 : SVDTRLanczosGetKSP - Retrieve the linear solver object (KSP) associated with
1885 : the SVD solver.
1886 :
1887 : Collective
1888 :
1889 : Input Parameter:
1890 : . svd - SVD solver
1891 :
1892 : Output Parameter:
1893 : . ksp - the linear solver object
1894 :
1895 : Level: advanced
1896 :
1897 : .seealso: SVDTRLanczosSetKSP()
1898 : @*/
1899 52 : PetscErrorCode SVDTRLanczosGetKSP(SVD svd,KSP *ksp)
1900 : {
1901 52 : PetscFunctionBegin;
1902 52 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1903 52 : PetscAssertPointer(ksp,2);
1904 52 : PetscUseMethod(svd,"SVDTRLanczosGetKSP_C",(SVD,KSP*),(svd,ksp));
1905 52 : PetscFunctionReturn(PETSC_SUCCESS);
1906 : }
1907 :
1908 14 : static PetscErrorCode SVDTRLanczosSetRestart_TRLanczos(SVD svd,PetscReal keep)
1909 : {
1910 14 : SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1911 :
1912 14 : PetscFunctionBegin;
1913 14 : if (keep==(PetscReal)PETSC_DEFAULT) ctx->keep = 0.5;
1914 : else {
1915 14 : PetscCheck(keep>=0.1 && keep<=0.9,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The keep argument %g must be in the range [0.1,0.9]",(double)keep);
1916 14 : ctx->keep = keep;
1917 : }
1918 14 : PetscFunctionReturn(PETSC_SUCCESS);
1919 : }
1920 :
1921 : /*@
1922 : SVDTRLanczosSetRestart - Sets the restart parameter for the thick-restart
1923 : Lanczos method, in particular the proportion of basis vectors that must be
1924 : kept after restart.
1925 :
1926 : Logically Collective
1927 :
1928 : Input Parameters:
1929 : + svd - the singular value solver
1930 : - keep - the number of vectors to be kept at restart
1931 :
1932 : Options Database Key:
1933 : . -svd_trlanczos_restart - Sets the restart parameter
1934 :
1935 : Notes:
1936 : Allowed values are in the range [0.1,0.9]. The default is 0.5.
1937 :
1938 : Level: advanced
1939 :
1940 : .seealso: SVDTRLanczosGetRestart()
1941 : @*/
1942 14 : PetscErrorCode SVDTRLanczosSetRestart(SVD svd,PetscReal keep)
1943 : {
1944 14 : PetscFunctionBegin;
1945 14 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1946 56 : PetscValidLogicalCollectiveReal(svd,keep,2);
1947 14 : PetscTryMethod(svd,"SVDTRLanczosSetRestart_C",(SVD,PetscReal),(svd,keep));
1948 14 : PetscFunctionReturn(PETSC_SUCCESS);
1949 : }
1950 :
1951 6 : static PetscErrorCode SVDTRLanczosGetRestart_TRLanczos(SVD svd,PetscReal *keep)
1952 : {
1953 6 : SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1954 :
1955 6 : PetscFunctionBegin;
1956 6 : *keep = ctx->keep;
1957 6 : PetscFunctionReturn(PETSC_SUCCESS);
1958 : }
1959 :
1960 : /*@
1961 : SVDTRLanczosGetRestart - Gets the restart parameter used in the thick-restart
1962 : Lanczos method.
1963 :
1964 : Not Collective
1965 :
1966 : Input Parameter:
1967 : . svd - the singular value solver
1968 :
1969 : Output Parameter:
1970 : . keep - the restart parameter
1971 :
1972 : Level: advanced
1973 :
1974 : .seealso: SVDTRLanczosSetRestart()
1975 : @*/
1976 6 : PetscErrorCode SVDTRLanczosGetRestart(SVD svd,PetscReal *keep)
1977 : {
1978 6 : PetscFunctionBegin;
1979 6 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
1980 6 : PetscAssertPointer(keep,2);
1981 6 : PetscUseMethod(svd,"SVDTRLanczosGetRestart_C",(SVD,PetscReal*),(svd,keep));
1982 6 : PetscFunctionReturn(PETSC_SUCCESS);
1983 : }
1984 :
1985 20 : static PetscErrorCode SVDTRLanczosSetLocking_TRLanczos(SVD svd,PetscBool lock)
1986 : {
1987 20 : SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
1988 :
1989 20 : PetscFunctionBegin;
1990 20 : ctx->lock = lock;
1991 20 : PetscFunctionReturn(PETSC_SUCCESS);
1992 : }
1993 :
1994 : /*@
1995 : SVDTRLanczosSetLocking - Choose between locking and non-locking variants of
1996 : the thick-restart Lanczos method.
1997 :
1998 : Logically Collective
1999 :
2000 : Input Parameters:
2001 : + svd - the singular value solver
2002 : - lock - true if the locking variant must be selected
2003 :
2004 : Options Database Key:
2005 : . -svd_trlanczos_locking - Sets the locking flag
2006 :
2007 : Notes:
2008 : The default is to lock converged singular triplets when the method restarts.
2009 : This behaviour can be changed so that all directions are kept in the
2010 : working subspace even if already converged to working accuracy (the
2011 : non-locking variant).
2012 :
2013 : Level: advanced
2014 :
2015 : .seealso: SVDTRLanczosGetLocking()
2016 : @*/
2017 20 : PetscErrorCode SVDTRLanczosSetLocking(SVD svd,PetscBool lock)
2018 : {
2019 20 : PetscFunctionBegin;
2020 20 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
2021 80 : PetscValidLogicalCollectiveBool(svd,lock,2);
2022 20 : PetscTryMethod(svd,"SVDTRLanczosSetLocking_C",(SVD,PetscBool),(svd,lock));
2023 20 : PetscFunctionReturn(PETSC_SUCCESS);
2024 : }
2025 :
2026 6 : static PetscErrorCode SVDTRLanczosGetLocking_TRLanczos(SVD svd,PetscBool *lock)
2027 : {
2028 6 : SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2029 :
2030 6 : PetscFunctionBegin;
2031 6 : *lock = ctx->lock;
2032 6 : PetscFunctionReturn(PETSC_SUCCESS);
2033 : }
2034 :
2035 : /*@
2036 : SVDTRLanczosGetLocking - Gets the locking flag used in the thick-restart
2037 : Lanczos method.
2038 :
2039 : Not Collective
2040 :
2041 : Input Parameter:
2042 : . svd - the singular value solver
2043 :
2044 : Output Parameter:
2045 : . lock - the locking flag
2046 :
2047 : Level: advanced
2048 :
2049 : .seealso: SVDTRLanczosSetLocking()
2050 : @*/
2051 6 : PetscErrorCode SVDTRLanczosGetLocking(SVD svd,PetscBool *lock)
2052 : {
2053 6 : PetscFunctionBegin;
2054 6 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
2055 6 : PetscAssertPointer(lock,2);
2056 6 : PetscUseMethod(svd,"SVDTRLanczosGetLocking_C",(SVD,PetscBool*),(svd,lock));
2057 6 : PetscFunctionReturn(PETSC_SUCCESS);
2058 : }
2059 :
2060 16 : static PetscErrorCode SVDTRLanczosSetExplicitMatrix_TRLanczos(SVD svd,PetscBool explicitmat)
2061 : {
2062 16 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2063 :
2064 16 : PetscFunctionBegin;
2065 16 : if (lanczos->explicitmatrix != explicitmat) {
2066 12 : lanczos->explicitmatrix = explicitmat;
2067 12 : svd->state = SVD_STATE_INITIAL;
2068 : }
2069 16 : PetscFunctionReturn(PETSC_SUCCESS);
2070 : }
2071 :
2072 : /*@
2073 : SVDTRLanczosSetExplicitMatrix - Indicate if the matrix Z=[A;B] must
2074 : be built explicitly.
2075 :
2076 : Logically Collective
2077 :
2078 : Input Parameters:
2079 : + svd - singular value solver
2080 : - explicitmat - Boolean flag indicating if Z=[A;B] is built explicitly
2081 :
2082 : Options Database Key:
2083 : . -svd_trlanczos_explicitmatrix <boolean> - Indicates the boolean flag
2084 :
2085 : Notes:
2086 : This option is relevant for the GSVD case only.
2087 : Z is the coefficient matrix of the KSP solver used internally.
2088 :
2089 : Level: advanced
2090 :
2091 : .seealso: SVDTRLanczosGetExplicitMatrix()
2092 : @*/
2093 16 : PetscErrorCode SVDTRLanczosSetExplicitMatrix(SVD svd,PetscBool explicitmat)
2094 : {
2095 16 : PetscFunctionBegin;
2096 16 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
2097 64 : PetscValidLogicalCollectiveBool(svd,explicitmat,2);
2098 16 : PetscTryMethod(svd,"SVDTRLanczosSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
2099 16 : PetscFunctionReturn(PETSC_SUCCESS);
2100 : }
2101 :
2102 0 : static PetscErrorCode SVDTRLanczosGetExplicitMatrix_TRLanczos(SVD svd,PetscBool *explicitmat)
2103 : {
2104 0 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2105 :
2106 0 : PetscFunctionBegin;
2107 0 : *explicitmat = lanczos->explicitmatrix;
2108 0 : PetscFunctionReturn(PETSC_SUCCESS);
2109 : }
2110 :
2111 : /*@
2112 : SVDTRLanczosGetExplicitMatrix - Returns the flag indicating if Z=[A;B] is built explicitly.
2113 :
2114 : Not Collective
2115 :
2116 : Input Parameter:
2117 : . svd - singular value solver
2118 :
2119 : Output Parameter:
2120 : . explicitmat - the mode flag
2121 :
2122 : Level: advanced
2123 :
2124 : .seealso: SVDTRLanczosSetExplicitMatrix()
2125 : @*/
2126 0 : PetscErrorCode SVDTRLanczosGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
2127 : {
2128 0 : PetscFunctionBegin;
2129 0 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
2130 0 : PetscAssertPointer(explicitmat,2);
2131 0 : PetscUseMethod(svd,"SVDTRLanczosGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
2132 0 : PetscFunctionReturn(PETSC_SUCCESS);
2133 : }
2134 :
2135 20 : static PetscErrorCode SVDTRLanczosSetScale_TRLanczos(SVD svd,PetscReal scale)
2136 : {
2137 20 : SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2138 :
2139 20 : PetscFunctionBegin;
2140 20 : if (scale<0) {
2141 11 : ctx->scalef = 1.0;
2142 11 : ctx->scaleth = -scale;
2143 : } else {
2144 9 : ctx->scalef = scale;
2145 9 : ctx->scaleth = 0.0;
2146 : }
2147 20 : PetscFunctionReturn(PETSC_SUCCESS);
2148 : }
2149 :
2150 : /*@
2151 : SVDTRLanczosSetScale - Sets the scale parameter for the GSVD.
2152 :
2153 : Logically Collective
2154 :
2155 : Input Parameters:
2156 : + svd - singular value solver
2157 : - scale - scale parameter
2158 :
2159 : Options Database Key:
2160 : . -svd_trlanczos_scale <real> - scale factor/threshold
2161 :
2162 : Notes:
2163 : This parameter is relevant for the GSVD case only. If the parameter is
2164 : positive, it indicates the scale factor for B in matrix Z=[A;B]. If
2165 : negative, its absolute value is the threshold for automatic scaling.
2166 : In automatic scaling, whenever the largest approximate generalized singular
2167 : value (or the inverse of the smallest value, if SVD_SMALLEST is used)
2168 : exceeds the threshold, the computation is restarted with matrix B
2169 : scaled by that value.
2170 :
2171 : Level: advanced
2172 :
2173 : .seealso: SVDTRLanczosGetScale()
2174 : @*/
2175 20 : PetscErrorCode SVDTRLanczosSetScale(SVD svd,PetscReal scale)
2176 : {
2177 20 : PetscFunctionBegin;
2178 20 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
2179 80 : PetscValidLogicalCollectiveReal(svd,scale,2);
2180 20 : PetscTryMethod(svd,"SVDTRLanczosSetScale_C",(SVD,PetscReal),(svd,scale));
2181 20 : PetscFunctionReturn(PETSC_SUCCESS);
2182 : }
2183 :
2184 101 : static PetscErrorCode SVDTRLanczosGetScale_TRLanczos(SVD svd,PetscReal *scale)
2185 : {
2186 101 : SVD_TRLANCZOS *ctx = (SVD_TRLANCZOS*)svd->data;
2187 :
2188 101 : PetscFunctionBegin;
2189 101 : if (ctx->scaleth==0) *scale = ctx->scalef;
2190 6 : else *scale = -ctx->scaleth;
2191 101 : PetscFunctionReturn(PETSC_SUCCESS);
2192 : }
2193 :
2194 : /*@
2195 : SVDTRLanczosGetScale - Gets the scale parameter for the GSVD.
2196 :
2197 : Not Collective
2198 :
2199 : Input Parameter:
2200 : . svd - the singular value solver
2201 :
2202 : Output Parameter:
2203 : . scale - the scale parameter
2204 :
2205 : Notes:
2206 : This parameter is relevant for the GSVD case only. If the parameter is
2207 : positive, it indicates the scale factor for B in matrix Z=[A;B]. If
2208 : negative, its absolute value is the threshold for automatic scaling.
2209 :
2210 : Level: advanced
2211 :
2212 : .seealso: SVDTRLanczosSetScale()
2213 : @*/
2214 101 : PetscErrorCode SVDTRLanczosGetScale(SVD svd,PetscReal *scale)
2215 : {
2216 101 : PetscFunctionBegin;
2217 101 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
2218 101 : PetscAssertPointer(scale,2);
2219 101 : PetscUseMethod(svd,"SVDTRLanczosGetScale_C",(SVD,PetscReal*),(svd,scale));
2220 101 : PetscFunctionReturn(PETSC_SUCCESS);
2221 : }
2222 :
2223 97 : static PetscErrorCode SVDReset_TRLanczos(SVD svd)
2224 : {
2225 97 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2226 :
2227 97 : PetscFunctionBegin;
2228 97 : if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) {
2229 58 : PetscCall(KSPReset(lanczos->ksp));
2230 58 : PetscCall(MatDestroy(&lanczos->Z));
2231 : }
2232 97 : PetscFunctionReturn(PETSC_SUCCESS);
2233 : }
2234 :
2235 96 : static PetscErrorCode SVDDestroy_TRLanczos(SVD svd)
2236 : {
2237 96 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2238 :
2239 96 : PetscFunctionBegin;
2240 96 : if (svd->isgeneralized || (!svd->problem_type && svd->OPb)) PetscCall(KSPDestroy(&lanczos->ksp));
2241 96 : PetscCall(PetscFree(svd->data));
2242 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",NULL));
2243 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",NULL));
2244 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",NULL));
2245 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",NULL));
2246 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",NULL));
2247 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",NULL));
2248 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",NULL));
2249 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",NULL));
2250 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",NULL));
2251 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",NULL));
2252 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",NULL));
2253 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",NULL));
2254 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetScale_C",NULL));
2255 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetScale_C",NULL));
2256 96 : PetscFunctionReturn(PETSC_SUCCESS);
2257 : }
2258 :
2259 1 : static PetscErrorCode SVDView_TRLanczos(SVD svd,PetscViewer viewer)
2260 : {
2261 1 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2262 1 : PetscBool isascii;
2263 :
2264 1 : PetscFunctionBegin;
2265 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
2266 1 : if (isascii) {
2267 1 : PetscCall(PetscViewerASCIIPrintf(viewer," %d%% of basis vectors kept after restart\n",(int)(100*lanczos->keep)));
2268 1 : PetscCall(PetscViewerASCIIPrintf(viewer," using the %slocking variant\n",lanczos->lock?"":"non-"));
2269 1 : if (svd->isgeneralized) {
2270 0 : const char *bidiag="";
2271 :
2272 0 : switch (lanczos->bidiag) {
2273 0 : case SVD_TRLANCZOS_GBIDIAG_SINGLE: bidiag = "single"; break;
2274 0 : case SVD_TRLANCZOS_GBIDIAG_UPPER: bidiag = "joint upper-upper"; break;
2275 0 : case SVD_TRLANCZOS_GBIDIAG_LOWER: bidiag = "joint lower-upper"; break;
2276 : }
2277 0 : PetscCall(PetscViewerASCIIPrintf(viewer," bidiagonalization choice: %s\n",bidiag));
2278 0 : PetscCall(PetscViewerASCIIPrintf(viewer," %s matrix\n",lanczos->explicitmatrix?"explicit":"implicit"));
2279 0 : if (lanczos->scaleth==0) PetscCall(PetscViewerASCIIPrintf(viewer," scale factor for matrix B: %g\n",(double)lanczos->scalef));
2280 0 : else PetscCall(PetscViewerASCIIPrintf(viewer," automatic scaling for matrix B with threshold: %g\n",(double)lanczos->scaleth));
2281 0 : if (!lanczos->ksp) PetscCall(SVDTRLanczosGetKSP(svd,&lanczos->ksp));
2282 0 : PetscCall(PetscViewerASCIIPushTab(viewer));
2283 0 : PetscCall(KSPView(lanczos->ksp,viewer));
2284 0 : PetscCall(PetscViewerASCIIPopTab(viewer));
2285 2 : } else PetscCall(PetscViewerASCIIPrintf(viewer," %s-sided reorthogonalization\n",lanczos->oneside? "one": "two"));
2286 : }
2287 1 : PetscFunctionReturn(PETSC_SUCCESS);
2288 : }
2289 :
2290 203 : static PetscErrorCode SVDSetDSType_TRLanczos(SVD svd)
2291 : {
2292 203 : SVD_TRLANCZOS *lanczos = (SVD_TRLANCZOS*)svd->data;
2293 203 : DSType dstype;
2294 :
2295 203 : PetscFunctionBegin;
2296 203 : dstype = svd->ishyperbolic? DSHSVD: DSSVD;
2297 203 : if (svd->OPb && (lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_UPPER || lanczos->bidiag==SVD_TRLANCZOS_GBIDIAG_LOWER)) dstype = DSGSVD;
2298 203 : PetscCall(DSSetType(svd->ds,dstype));
2299 203 : PetscFunctionReturn(PETSC_SUCCESS);
2300 : }
2301 :
2302 96 : SLEPC_EXTERN PetscErrorCode SVDCreate_TRLanczos(SVD svd)
2303 : {
2304 96 : SVD_TRLANCZOS *ctx;
2305 :
2306 96 : PetscFunctionBegin;
2307 96 : PetscCall(PetscNew(&ctx));
2308 96 : svd->data = (void*)ctx;
2309 :
2310 96 : ctx->lock = PETSC_TRUE;
2311 96 : ctx->bidiag = SVD_TRLANCZOS_GBIDIAG_LOWER;
2312 96 : ctx->scalef = 1.0;
2313 96 : ctx->scaleth = 0.0;
2314 :
2315 96 : svd->ops->setup = SVDSetUp_TRLanczos;
2316 96 : svd->ops->solve = SVDSolve_TRLanczos;
2317 96 : svd->ops->solveg = SVDSolve_TRLanczos_GSVD;
2318 96 : svd->ops->solveh = SVDSolve_TRLanczos_HSVD;
2319 96 : svd->ops->destroy = SVDDestroy_TRLanczos;
2320 96 : svd->ops->reset = SVDReset_TRLanczos;
2321 96 : svd->ops->setfromoptions = SVDSetFromOptions_TRLanczos;
2322 96 : svd->ops->view = SVDView_TRLanczos;
2323 96 : svd->ops->setdstype = SVDSetDSType_TRLanczos;
2324 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetOneSide_C",SVDTRLanczosSetOneSide_TRLanczos));
2325 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetOneSide_C",SVDTRLanczosGetOneSide_TRLanczos));
2326 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetGBidiag_C",SVDTRLanczosSetGBidiag_TRLanczos));
2327 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetGBidiag_C",SVDTRLanczosGetGBidiag_TRLanczos));
2328 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetKSP_C",SVDTRLanczosSetKSP_TRLanczos));
2329 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetKSP_C",SVDTRLanczosGetKSP_TRLanczos));
2330 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetRestart_C",SVDTRLanczosSetRestart_TRLanczos));
2331 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetRestart_C",SVDTRLanczosGetRestart_TRLanczos));
2332 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetLocking_C",SVDTRLanczosSetLocking_TRLanczos));
2333 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetLocking_C",SVDTRLanczosGetLocking_TRLanczos));
2334 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetExplicitMatrix_C",SVDTRLanczosSetExplicitMatrix_TRLanczos));
2335 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetExplicitMatrix_C",SVDTRLanczosGetExplicitMatrix_TRLanczos));
2336 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosSetScale_C",SVDTRLanczosSetScale_TRLanczos));
2337 96 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDTRLanczosGetScale_C",SVDTRLanczosGetScale_TRLanczos));
2338 96 : PetscFunctionReturn(PETSC_SUCCESS);
2339 : }
|
