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 : SVD routines related to the solution process
12 : */
13 :
14 : #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
15 : #include <slepc/private/bvimpl.h>
16 :
17 : /*
18 : SVDComputeVectors_Left - Compute left singular vectors as U=A*V.
19 : Only done if the leftbasis flag is false. Assumes V is available.
20 : */
21 60 : PetscErrorCode SVDComputeVectors_Left(SVD svd)
22 : {
23 60 : Vec tl,omega2,u,v,w;
24 60 : PetscInt i,oldsize;
25 60 : VecType vtype;
26 60 : const PetscScalar* varray;
27 :
28 60 : PetscFunctionBegin;
29 60 : if (!svd->leftbasis) {
30 : /* generate left singular vectors on U */
31 43 : if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
32 43 : PetscCall(BVGetSizes(svd->U,NULL,NULL,&oldsize));
33 43 : if (!oldsize) {
34 38 : if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
35 38 : PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
36 38 : PetscCall(BVSetSizesFromVec(svd->U,tl,svd->ncv));
37 38 : PetscCall(VecDestroy(&tl));
38 : }
39 43 : PetscCall(BVSetActiveColumns(svd->V,0,svd->nconv));
40 43 : PetscCall(BVSetActiveColumns(svd->U,0,svd->nconv));
41 43 : if (!svd->ishyperbolic) PetscCall(BVMatMult(svd->V,svd->A,svd->U));
42 13 : else if (svd->swapped) { /* compute right singular vectors as V=A'*Omega*U */
43 4 : PetscCall(MatCreateVecs(svd->A,&w,NULL));
44 120 : for (i=0;i<svd->nconv;i++) {
45 116 : PetscCall(BVGetColumn(svd->V,i,&v));
46 116 : PetscCall(BVGetColumn(svd->U,i,&u));
47 116 : PetscCall(VecPointwiseMult(w,v,svd->omega));
48 116 : PetscCall(MatMult(svd->A,w,u));
49 116 : PetscCall(BVRestoreColumn(svd->V,i,&v));
50 116 : PetscCall(BVRestoreColumn(svd->U,i,&u));
51 : }
52 4 : PetscCall(VecDestroy(&w));
53 : } else { /* compute left singular vectors as usual U=A*V, and set-up Omega-orthogonalization of U */
54 9 : PetscCall(BV_SetMatrixDiagonal(svd->U,svd->omega,svd->A));
55 9 : PetscCall(BVMatMult(svd->V,svd->A,svd->U));
56 : }
57 43 : PetscCall(BVOrthogonalize(svd->U,NULL));
58 43 : if (svd->ishyperbolic && !svd->swapped) { /* store signature after Omega-orthogonalization */
59 9 : PetscCall(MatGetVecType(svd->A,&vtype));
60 9 : PetscCall(VecCreate(PETSC_COMM_SELF,&omega2));
61 9 : PetscCall(VecSetSizes(omega2,svd->nconv,svd->nconv));
62 9 : PetscCall(VecSetType(omega2,vtype));
63 9 : PetscCall(BVGetSignature(svd->U,omega2));
64 9 : PetscCall(VecGetArrayRead(omega2,&varray));
65 374 : for (i=0;i<svd->nconv;i++) {
66 365 : svd->sign[i] = PetscRealPart(varray[i]);
67 365 : if (PetscRealPart(varray[i])<0.0) PetscCall(BVScaleColumn(svd->U,i,-1.0));
68 : }
69 9 : PetscCall(VecRestoreArrayRead(omega2,&varray));
70 9 : PetscCall(VecDestroy(&omega2));
71 : }
72 : }
73 60 : PetscFunctionReturn(PETSC_SUCCESS);
74 : }
75 :
76 2804 : PetscErrorCode SVDComputeVectors(SVD svd)
77 : {
78 2804 : PetscFunctionBegin;
79 2804 : SVDCheckSolved(svd,1);
80 2804 : if (svd->state==SVD_STATE_SOLVED) PetscTryTypeMethod(svd,computevectors);
81 2804 : svd->state = SVD_STATE_VECTORS;
82 2804 : PetscFunctionReturn(PETSC_SUCCESS);
83 : }
84 :
85 : /*@
86 : SVDSolve - Solves the singular value problem.
87 :
88 : Collective
89 :
90 : Input Parameter:
91 : . svd - singular value solver context obtained from SVDCreate()
92 :
93 : Options Database Keys:
94 : + -svd_view - print information about the solver used
95 : . -svd_view_mat0 - view the first matrix (A)
96 : . -svd_view_mat1 - view the second matrix (B)
97 : . -svd_view_signature - view the signature matrix (omega)
98 : . -svd_view_vectors - view the computed singular vectors
99 : . -svd_view_values - view the computed singular values
100 : . -svd_converged_reason - print reason for convergence, and number of iterations
101 : . -svd_error_absolute - print absolute errors of each singular triplet
102 : . -svd_error_relative - print relative errors of each singular triplet
103 : - -svd_error_norm - print errors relative to the matrix norms of each singular triplet
104 :
105 : Notes:
106 : All the command-line options listed above admit an optional argument specifying
107 : the viewer type and options. For instance, use '-svd_view_mat0 binary:amatrix.bin'
108 : to save the A matrix to a binary file, '-svd_view_values draw' to draw the computed
109 : singular values graphically, or '-svd_error_relative :myerr.m:ascii_matlab' to save
110 : the errors in a file that can be executed in Matlab.
111 :
112 : Level: beginner
113 :
114 : .seealso: SVDCreate(), SVDSetUp(), SVDDestroy()
115 : @*/
116 286 : PetscErrorCode SVDSolve(SVD svd)
117 : {
118 286 : PetscInt i,*workperm;
119 :
120 286 : PetscFunctionBegin;
121 286 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
122 286 : if (svd->state>=SVD_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
123 286 : PetscCall(PetscLogEventBegin(SVD_Solve,svd,0,0,0));
124 :
125 : /* call setup */
126 286 : PetscCall(SVDSetUp(svd));
127 286 : svd->its = 0;
128 286 : svd->nconv = 0;
129 6417 : for (i=0;i<svd->ncv;i++) {
130 6131 : svd->sigma[i] = 0.0;
131 6131 : svd->errest[i] = 0.0;
132 6131 : svd->perm[i] = i;
133 : }
134 286 : PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view_pre"));
135 :
136 286 : switch (svd->problem_type) {
137 146 : case SVD_STANDARD:
138 146 : PetscUseTypeMethod(svd,solve);
139 : break;
140 102 : case SVD_GENERALIZED:
141 102 : PetscUseTypeMethod(svd,solveg);
142 : break;
143 38 : case SVD_HYPERBOLIC:
144 38 : PetscUseTypeMethod(svd,solveh);
145 : break;
146 : }
147 286 : svd->state = SVD_STATE_SOLVED;
148 :
149 : /* sort singular triplets */
150 286 : if (svd->which == SVD_SMALLEST) PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm));
151 : else {
152 247 : PetscCall(PetscMalloc1(svd->nconv,&workperm));
153 3094 : for (i=0;i<svd->nconv;i++) workperm[i] = i;
154 247 : PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm));
155 3094 : for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
156 247 : PetscCall(PetscFree(workperm));
157 : }
158 286 : PetscCall(PetscLogEventEnd(SVD_Solve,svd,0,0,0));
159 :
160 : /* various viewers */
161 286 : PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view"));
162 286 : PetscCall(SVDConvergedReasonViewFromOptions(svd));
163 286 : PetscCall(SVDErrorViewFromOptions(svd));
164 286 : PetscCall(SVDValuesViewFromOptions(svd));
165 286 : PetscCall(SVDVectorsViewFromOptions(svd));
166 286 : PetscCall(MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0"));
167 286 : if (svd->isgeneralized) PetscCall(MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1"));
168 286 : if (svd->ishyperbolic) PetscCall(VecViewFromOptions(svd->omega,(PetscObject)svd,"-svd_view_signature"));
169 :
170 : /* Remove the initial subspaces */
171 286 : svd->nini = 0;
172 286 : svd->ninil = 0;
173 286 : PetscFunctionReturn(PETSC_SUCCESS);
174 : }
175 :
176 : /*@
177 : SVDGetIterationNumber - Gets the current iteration number. If the
178 : call to SVDSolve() is complete, then it returns the number of iterations
179 : carried out by the solution method.
180 :
181 : Not Collective
182 :
183 : Input Parameter:
184 : . svd - the singular value solver context
185 :
186 : Output Parameter:
187 : . its - number of iterations
188 :
189 : Note:
190 : During the i-th iteration this call returns i-1. If SVDSolve() is
191 : complete, then parameter "its" contains either the iteration number at
192 : which convergence was successfully reached, or failure was detected.
193 : Call SVDGetConvergedReason() to determine if the solver converged or
194 : failed and why.
195 :
196 : Level: intermediate
197 :
198 : .seealso: SVDGetConvergedReason(), SVDSetTolerances()
199 : @*/
200 68 : PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)
201 : {
202 68 : PetscFunctionBegin;
203 68 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
204 68 : PetscAssertPointer(its,2);
205 68 : *its = svd->its;
206 68 : PetscFunctionReturn(PETSC_SUCCESS);
207 : }
208 :
209 : /*@
210 : SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
211 : stopped.
212 :
213 : Not Collective
214 :
215 : Input Parameter:
216 : . svd - the singular value solver context
217 :
218 : Output Parameter:
219 : . reason - negative value indicates diverged, positive value converged
220 : (see SVDConvergedReason)
221 :
222 : Options Database Key:
223 : . -svd_converged_reason - print the reason to a viewer
224 :
225 : Notes:
226 : Possible values for reason are
227 : + SVD_CONVERGED_TOL - converged up to tolerance
228 : . SVD_CONVERGED_USER - converged due to a user-defined condition
229 : . SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
230 : . SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
231 : . SVD_DIVERGED_BREAKDOWN - generic breakdown in method
232 : - SVD_DIVERGED_SYMMETRY_LOST - underlying indefinite eigensolver was not able to keep symmetry
233 :
234 : Can only be called after the call to SVDSolve() is complete.
235 :
236 : Level: intermediate
237 :
238 : .seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason
239 : @*/
240 13 : PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
241 : {
242 13 : PetscFunctionBegin;
243 13 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
244 13 : PetscAssertPointer(reason,2);
245 13 : SVDCheckSolved(svd,1);
246 13 : *reason = svd->reason;
247 13 : PetscFunctionReturn(PETSC_SUCCESS);
248 : }
249 :
250 : /*@
251 : SVDGetConverged - Gets the number of converged singular values.
252 :
253 : Not Collective
254 :
255 : Input Parameter:
256 : . svd - the singular value solver context
257 :
258 : Output Parameter:
259 : . nconv - number of converged singular values
260 :
261 : Note:
262 : This function should be called after SVDSolve() has finished.
263 :
264 : Level: beginner
265 :
266 : .seealso: SVDSetDimensions(), SVDSolve(), SVDGetSingularTriplet()
267 : @*/
268 129 : PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)
269 : {
270 129 : PetscFunctionBegin;
271 129 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
272 129 : PetscAssertPointer(nconv,2);
273 129 : SVDCheckSolved(svd,1);
274 129 : *nconv = svd->nconv;
275 129 : PetscFunctionReturn(PETSC_SUCCESS);
276 : }
277 :
278 : /*@
279 : SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
280 : as computed by SVDSolve(). The solution consists in the singular value and its left
281 : and right singular vectors.
282 :
283 : Collective
284 :
285 : Input Parameters:
286 : + svd - singular value solver context
287 : - i - index of the solution
288 :
289 : Output Parameters:
290 : + sigma - singular value
291 : . u - left singular vector
292 : - v - right singular vector
293 :
294 : Note:
295 : Both u or v can be NULL if singular vectors are not required.
296 : Otherwise, the caller must provide valid Vec objects, i.e.,
297 : they must be created by the calling program with e.g. MatCreateVecs().
298 :
299 : The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
300 : Singular triplets are indexed according to the ordering criterion established
301 : with SVDSetWhichSingularTriplets().
302 :
303 : In the case of GSVD, the solution consists in three vectors u,v,x that are
304 : returned as follows. Vector x is returned in the right singular vector
305 : (argument v) and has length equal to the number of columns of A and B.
306 : The other two vectors are returned stacked on top of each other [u;v] in
307 : the left singular vector argument, with length equal to m+n (number of rows
308 : of A plus number of rows of B).
309 :
310 : Level: beginner
311 :
312 : .seealso: SVDSolve(), SVDGetConverged(), SVDSetWhichSingularTriplets()
313 : @*/
314 4246 : PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)
315 : {
316 4246 : PetscInt M,N;
317 4246 : Vec w;
318 :
319 4246 : PetscFunctionBegin;
320 4246 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
321 12738 : PetscValidLogicalCollectiveInt(svd,i,2);
322 4246 : SVDCheckSolved(svd,1);
323 4246 : if (u) { PetscValidHeaderSpecific(u,VEC_CLASSID,4); PetscCheckSameComm(svd,1,u,4); }
324 4246 : if (v) { PetscValidHeaderSpecific(v,VEC_CLASSID,5); PetscCheckSameComm(svd,1,v,5); }
325 4246 : PetscCheck(i>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
326 4246 : PetscCheck(i<svd->nconv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see SVDGetConverged()");
327 4246 : if (sigma) *sigma = svd->sigma[svd->perm[i]];
328 4246 : if (u || v) {
329 2803 : if (!svd->isgeneralized) {
330 2171 : PetscCall(MatGetSize(svd->OP,&M,&N));
331 2171 : if (M<N) { w = u; u = v; v = w; }
332 : }
333 2803 : PetscCall(SVDComputeVectors(svd));
334 2803 : if (u) PetscCall(BVCopyVec(svd->U,svd->perm[i],u));
335 2803 : if (v) PetscCall(BVCopyVec(svd->V,svd->perm[i],v));
336 : }
337 4246 : PetscFunctionReturn(PETSC_SUCCESS);
338 : }
339 :
340 : /*
341 : SVDComputeResidualNorms_Standard - Computes the norms of the left and
342 : right residuals associated with the i-th computed singular triplet.
343 :
344 : Input Parameters:
345 : sigma - singular value
346 : u,v - singular vectors
347 : x,y - two work vectors with the same dimensions as u,v
348 : */
349 479 : static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)
350 : {
351 479 : PetscInt M,N;
352 :
353 479 : PetscFunctionBegin;
354 : /* norm1 = ||A*v-sigma*u||_2 */
355 479 : if (norm1) {
356 479 : PetscCall(MatMult(svd->OP,v,x));
357 479 : PetscCall(VecAXPY(x,-sigma,u));
358 479 : PetscCall(VecNorm(x,NORM_2,norm1));
359 : }
360 : /* norm2 = ||A^T*u-sigma*v||_2 */
361 479 : if (norm2) {
362 479 : PetscCall(MatGetSize(svd->OP,&M,&N));
363 479 : if (M<N) PetscCall(MatMult(svd->A,u,y));
364 408 : else PetscCall(MatMult(svd->AT,u,y));
365 479 : PetscCall(VecAXPY(y,-sigma,v));
366 479 : PetscCall(VecNorm(y,NORM_2,norm2));
367 : }
368 479 : PetscFunctionReturn(PETSC_SUCCESS);
369 : }
370 :
371 : /*
372 : SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
373 : norm1 = ||s^2*A'*u-c*B'*B*x||_2 and norm2 = ||c^2*B'*v-s*A'*A*x||_2.
374 :
375 : Input Parameters:
376 : sigma - singular value
377 : uv - left singular vectors [u;v]
378 : x - right singular vector
379 : y,z - two work vectors with the same dimension as x
380 : */
381 431 : static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
382 : {
383 431 : Vec u,v,ax,bx,nest,aux[2];
384 431 : PetscReal c,s;
385 :
386 431 : PetscFunctionBegin;
387 431 : PetscCall(MatCreateVecs(svd->OP,NULL,&u));
388 431 : PetscCall(MatCreateVecs(svd->OPb,NULL,&v));
389 431 : aux[0] = u;
390 431 : aux[1] = v;
391 431 : PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
392 431 : PetscCall(VecCopy(uv,nest));
393 :
394 431 : s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
395 431 : c = sigma*s;
396 :
397 : /* norm1 = ||s^2*A'*u-c*B'*B*x||_2 */
398 431 : if (norm1) {
399 431 : PetscCall(VecDuplicate(v,&bx));
400 431 : PetscCall(MatMultHermitianTranspose(svd->OP,u,z));
401 431 : PetscCall(MatMult(svd->OPb,x,bx));
402 431 : PetscCall(MatMultHermitianTranspose(svd->OPb,bx,y));
403 431 : PetscCall(VecAXPBY(y,s*s,-c,z));
404 431 : PetscCall(VecNorm(y,NORM_2,norm1));
405 431 : PetscCall(VecDestroy(&bx));
406 : }
407 : /* norm2 = ||c^2*B'*v-s*A'*A*x||_2 */
408 431 : if (norm2) {
409 431 : PetscCall(VecDuplicate(u,&ax));
410 431 : PetscCall(MatMultHermitianTranspose(svd->OPb,v,z));
411 431 : PetscCall(MatMult(svd->OP,x,ax));
412 431 : PetscCall(MatMultHermitianTranspose(svd->OP,ax,y));
413 431 : PetscCall(VecAXPBY(y,c*c,-s,z));
414 431 : PetscCall(VecNorm(y,NORM_2,norm2));
415 431 : PetscCall(VecDestroy(&ax));
416 : }
417 :
418 431 : PetscCall(VecDestroy(&v));
419 431 : PetscCall(VecDestroy(&u));
420 431 : PetscCall(VecDestroy(&nest));
421 431 : PetscFunctionReturn(PETSC_SUCCESS);
422 : }
423 :
424 : /*
425 : SVDComputeResidualNorms_Hyperbolic - Computes the norms of the left and
426 : right residuals associated with the i-th computed singular triplet.
427 :
428 : Input Parameters:
429 : sigma - singular value
430 : sign - corresponding element of the signature Omega2
431 : u,v - singular vectors
432 : x,y,z - three work vectors with the same dimensions as u,v,u
433 : */
434 730 : static PetscErrorCode SVDComputeResidualNorms_Hyperbolic(SVD svd,PetscReal sigma,PetscReal sign,Vec u,Vec v,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
435 : {
436 730 : PetscInt M,N;
437 :
438 730 : PetscFunctionBegin;
439 : /* norm1 = ||A*v-sigma*u||_2 */
440 730 : if (norm1) {
441 730 : PetscCall(MatMult(svd->OP,v,x));
442 730 : PetscCall(VecAXPY(x,-sigma,u));
443 730 : PetscCall(VecNorm(x,NORM_2,norm1));
444 : }
445 : /* norm2 = ||A^T*Omega*u-sigma*sign*v||_2 */
446 730 : if (norm2) {
447 730 : PetscCall(MatGetSize(svd->OP,&M,&N));
448 730 : PetscCall(VecPointwiseMult(z,u,svd->omega));
449 730 : if (M<N) PetscCall(MatMult(svd->A,z,y));
450 670 : else PetscCall(MatMult(svd->AT,z,y));
451 730 : PetscCall(VecAXPY(y,-sigma*sign,v));
452 730 : PetscCall(VecNorm(y,NORM_2,norm2));
453 : }
454 730 : PetscFunctionReturn(PETSC_SUCCESS);
455 : }
456 :
457 : /*@
458 : SVDComputeError - Computes the error (based on the residual norm) associated
459 : with the i-th singular triplet.
460 :
461 : Collective
462 :
463 : Input Parameters:
464 : + svd - the singular value solver context
465 : . i - the solution index
466 : - type - the type of error to compute
467 :
468 : Output Parameter:
469 : . error - the error
470 :
471 : Notes:
472 : The error can be computed in various ways, all of them based on the residual
473 : norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
474 : n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
475 : singular vector and v is the right singular vector.
476 :
477 : In the case of the GSVD, the two components of the residual norm are
478 : n1 = ||s^2*A'*u-c*B'*B*x||_2 and n2 = ||c^2*B'*v-s*A'*A*x||_2, where [u;v]
479 : are the left singular vectors and x is the right singular vector, with
480 : sigma=c/s.
481 :
482 : Level: beginner
483 :
484 : .seealso: SVDErrorType, SVDSolve()
485 : @*/
486 1640 : PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)
487 : {
488 1640 : PetscReal sigma,norm1,norm2,c,s;
489 1640 : Vec u=NULL,v=NULL,x=NULL,y=NULL,z=NULL;
490 1640 : PetscReal vecnorm=1.0;
491 :
492 1640 : PetscFunctionBegin;
493 1640 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
494 4920 : PetscValidLogicalCollectiveInt(svd,i,2);
495 4920 : PetscValidLogicalCollectiveEnum(svd,type,3);
496 1640 : PetscAssertPointer(error,4);
497 1640 : SVDCheckSolved(svd,1);
498 :
499 : /* allocate work vectors */
500 1640 : switch (svd->problem_type) {
501 479 : case SVD_STANDARD:
502 479 : PetscCall(SVDSetWorkVecs(svd,2,2));
503 479 : u = svd->workl[0];
504 479 : v = svd->workr[0];
505 479 : x = svd->workl[1];
506 479 : y = svd->workr[1];
507 479 : break;
508 431 : case SVD_GENERALIZED:
509 431 : PetscCall(SVDSetWorkVecs(svd,1,3));
510 431 : u = svd->workl[0];
511 431 : v = svd->workr[0];
512 431 : x = svd->workr[1];
513 431 : y = svd->workr[2];
514 431 : break;
515 730 : case SVD_HYPERBOLIC:
516 730 : PetscCall(SVDSetWorkVecs(svd,3,2));
517 730 : u = svd->workl[0];
518 730 : v = svd->workr[0];
519 730 : x = svd->workl[1];
520 730 : y = svd->workr[1];
521 730 : z = svd->workl[2];
522 730 : break;
523 : }
524 :
525 : /* compute residual norm */
526 1640 : PetscCall(SVDGetSingularTriplet(svd,i,&sigma,u,v));
527 1640 : switch (svd->problem_type) {
528 479 : case SVD_STANDARD:
529 479 : PetscCall(SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2));
530 : break;
531 431 : case SVD_GENERALIZED:
532 431 : PetscCall(SVDComputeResidualNorms_Generalized(svd,sigma,u,v,x,y,&norm1,&norm2));
533 : break;
534 730 : case SVD_HYPERBOLIC:
535 730 : PetscCall(SVDComputeResidualNorms_Hyperbolic(svd,sigma,svd->sign[svd->perm[i]],u,v,x,y,z,&norm1,&norm2));
536 : break;
537 : }
538 1640 : *error = SlepcAbs(norm1,norm2);
539 :
540 : /* compute 2-norm of eigenvector of the cyclic form */
541 1640 : if (type!=SVD_ERROR_ABSOLUTE) {
542 1627 : switch (svd->problem_type) {
543 466 : case SVD_STANDARD:
544 466 : vecnorm = PETSC_SQRT2;
545 466 : break;
546 431 : case SVD_GENERALIZED:
547 431 : PetscCall(VecNorm(v,NORM_2,&vecnorm));
548 431 : vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
549 431 : break;
550 730 : case SVD_HYPERBOLIC:
551 730 : PetscCall(VecNorm(u,NORM_2,&vecnorm));
552 730 : vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
553 730 : break;
554 : }
555 13 : }
556 :
557 : /* compute error */
558 1640 : switch (type) {
559 : case SVD_ERROR_ABSOLUTE:
560 : break;
561 459 : case SVD_ERROR_RELATIVE:
562 459 : if (svd->isgeneralized) {
563 0 : s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
564 0 : c = sigma*s;
565 0 : norm1 /= c*vecnorm;
566 0 : norm2 /= s*vecnorm;
567 0 : *error = PetscMax(norm1,norm2);
568 459 : } else *error /= sigma*vecnorm;
569 : break;
570 1168 : case SVD_ERROR_NORM:
571 1168 : if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
572 1168 : if (svd->isgeneralized && !svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
573 1168 : *error /= PetscMax(svd->nrma,svd->nrmb)*vecnorm;
574 1168 : break;
575 0 : default:
576 0 : SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
577 : }
578 1640 : PetscFunctionReturn(PETSC_SUCCESS);
579 : }
|
