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