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