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