Actual source code: svdsolve.c

slepc-3.17.2 2022-08-09
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

  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: */

 14: #include <slepc/private/svdimpl.h>

 16: /*
 17:   SVDComputeVectors_Left - Compute left singular vectors as U=A*V.
 18:   Only done if the leftbasis flag is false. Assumes V is available.
 19:  */
 20: PetscErrorCode SVDComputeVectors_Left(SVD svd)
 21: {
 22:   Vec            tl;
 23:   PetscInt       oldsize;

 25:   if (!svd->leftbasis) {
 26:     /* generate left singular vectors on U */
 27:     if (!svd->U) SVDGetBV(svd,NULL,&svd->U);
 28:     BVGetSizes(svd->U,NULL,NULL,&oldsize);
 29:     if (!oldsize) {
 30:       if (!((PetscObject)(svd->U))->type_name) BVSetType(svd->U,((PetscObject)(svd->V))->type_name);
 31:       MatCreateVecsEmpty(svd->A,NULL,&tl);
 32:       BVSetSizesFromVec(svd->U,tl,svd->ncv);
 33:       VecDestroy(&tl);
 34:     }
 35:     BVSetActiveColumns(svd->V,0,svd->nconv);
 36:     BVSetActiveColumns(svd->U,0,svd->nconv);
 37:     BVMatMult(svd->V,svd->A,svd->U);
 38:     BVOrthogonalize(svd->U,NULL);
 39:   }
 40:   PetscFunctionReturn(0);
 41: }

 43: PetscErrorCode SVDComputeVectors(SVD svd)
 44: {
 45:   SVDCheckSolved(svd,1);
 46:   if (svd->state==SVD_STATE_SOLVED && svd->ops->computevectors) (*svd->ops->computevectors)(svd);
 47:   svd->state = SVD_STATE_VECTORS;
 48:   PetscFunctionReturn(0);
 49: }

 51: /*@
 52:    SVDSolve - Solves the singular value problem.

 54:    Collective on svd

 56:    Input Parameter:
 57: .  svd - singular value solver context obtained from SVDCreate()

 59:    Options Database Keys:
 60: +  -svd_view - print information about the solver used
 61: .  -svd_view_mat0 - view the first matrix (A)
 62: .  -svd_view_mat1 - view the second matrix (B)
 63: .  -svd_view_vectors - view the computed singular vectors
 64: .  -svd_view_values - view the computed singular values
 65: .  -svd_converged_reason - print reason for convergence, and number of iterations
 66: .  -svd_error_absolute - print absolute errors of each singular triplet
 67: .  -svd_error_relative - print relative errors of each singular triplet
 68: -  -svd_error_norm     - print errors relative to the matrix norms of each singular triplet

 70:    Notes:
 71:    All the command-line options listed above admit an optional argument specifying
 72:    the viewer type and options. For instance, use '-svd_view_mat0 binary:amatrix.bin'
 73:    to save the A matrix to a binary file, '-svd_view_values draw' to draw the computed
 74:    singular values graphically, or '-svd_error_relative :myerr.m:ascii_matlab' to save
 75:    the errors in a file that can be executed in Matlab.

 77:    Level: beginner

 79: .seealso: SVDCreate(), SVDSetUp(), SVDDestroy()
 80: @*/
 81: PetscErrorCode SVDSolve(SVD svd)
 82: {
 83:   PetscInt       i,*workperm;

 86:   if (svd->state>=SVD_STATE_SOLVED) PetscFunctionReturn(0);
 87:   PetscLogEventBegin(SVD_Solve,svd,0,0,0);

 89:   /* call setup */
 90:   SVDSetUp(svd);
 91:   svd->its = 0;
 92:   svd->nconv = 0;
 93:   for (i=0;i<svd->ncv;i++) {
 94:     svd->sigma[i]  = 0.0;
 95:     svd->errest[i] = 0.0;
 96:     svd->perm[i]   = i;
 97:   }
 98:   SVDViewFromOptions(svd,NULL,"-svd_view_pre");

100:   switch (svd->problem_type) {
101:     case SVD_STANDARD:
102:       (*svd->ops->solve)(svd);
103:       break;
104:     case SVD_GENERALIZED:
105:       (*svd->ops->solveg)(svd);
106:       break;
107:   }
108:   svd->state = SVD_STATE_SOLVED;

110:   /* sort singular triplets */
111:   if (svd->which == SVD_SMALLEST) PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm);
112:   else {
113:     PetscMalloc1(svd->nconv,&workperm);
114:     for (i=0;i<svd->nconv;i++) workperm[i] = i;
115:     PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm);
116:     for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
117:     PetscFree(workperm);
118:   }
119:   PetscLogEventEnd(SVD_Solve,svd,0,0,0);

121:   /* various viewers */
122:   SVDViewFromOptions(svd,NULL,"-svd_view");
123:   SVDConvergedReasonViewFromOptions(svd);
124:   SVDErrorViewFromOptions(svd);
125:   SVDValuesViewFromOptions(svd);
126:   SVDVectorsViewFromOptions(svd);
127:   MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0");
128:   if (svd->isgeneralized) MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1");

130:   /* Remove the initial subspaces */
131:   svd->nini = 0;
132:   svd->ninil = 0;
133:   PetscFunctionReturn(0);
134: }

136: /*@
137:    SVDGetIterationNumber - Gets the current iteration number. If the
138:    call to SVDSolve() is complete, then it returns the number of iterations
139:    carried out by the solution method.

141:    Not Collective

143:    Input Parameter:
144: .  svd - the singular value solver context

146:    Output Parameter:
147: .  its - number of iterations

149:    Note:
150:    During the i-th iteration this call returns i-1. If SVDSolve() is
151:    complete, then parameter "its" contains either the iteration number at
152:    which convergence was successfully reached, or failure was detected.
153:    Call SVDGetConvergedReason() to determine if the solver converged or
154:    failed and why.

156:    Level: intermediate

158: .seealso: SVDGetConvergedReason(), SVDSetTolerances()
159: @*/
160: PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)
161: {
164:   *its = svd->its;
165:   PetscFunctionReturn(0);
166: }

168: /*@
169:    SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
170:    stopped.

172:    Not Collective

174:    Input Parameter:
175: .  svd - the singular value solver context

177:    Output Parameter:
178: .  reason - negative value indicates diverged, positive value converged
179:    (see SVDConvergedReason)

181:    Options Database Key:
182: .  -svd_converged_reason - print the reason to a viewer

184:    Notes:
185:    Possible values for reason are
186: +  SVD_CONVERGED_TOL - converged up to tolerance
187: .  SVD_CONVERGED_USER - converged due to a user-defined condition
188: .  SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
189: .  SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
190: -  SVD_DIVERGED_BREAKDOWN - generic breakdown in method

192:    Can only be called after the call to SVDSolve() is complete.

194:    Level: intermediate

196: .seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason
197: @*/
198: PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
199: {
202:   SVDCheckSolved(svd,1);
203:   *reason = svd->reason;
204:   PetscFunctionReturn(0);
205: }

207: /*@
208:    SVDGetConverged - Gets the number of converged singular values.

210:    Not Collective

212:    Input Parameter:
213: .  svd - the singular value solver context

215:    Output Parameter:
216: .  nconv - number of converged singular values

218:    Note:
219:    This function should be called after SVDSolve() has finished.

221:    Level: beginner

223: .seealso: SVDSetDimensions(), SVDSolve(), SVDGetSingularTriplet()
224: @*/
225: PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)
226: {
229:   SVDCheckSolved(svd,1);
230:   *nconv = svd->nconv;
231:   PetscFunctionReturn(0);
232: }

234: /*@C
235:    SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
236:    as computed by SVDSolve(). The solution consists in the singular value and its left
237:    and right singular vectors.

239:    Not Collective, but vectors are shared by all processors that share the SVD

241:    Input Parameters:
242: +  svd - singular value solver context
243: -  i   - index of the solution

245:    Output Parameters:
246: +  sigma - singular value
247: .  u     - left singular vector
248: -  v     - right singular vector

250:    Note:
251:    Both u or v can be NULL if singular vectors are not required.
252:    Otherwise, the caller must provide valid Vec objects, i.e.,
253:    they must be created by the calling program with e.g. MatCreateVecs().

255:    The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
256:    Singular triplets are indexed according to the ordering criterion established
257:    with SVDSetWhichSingularTriplets().

259:    In the case of GSVD, the solution consists in three vectors u,v,x that are
260:    returned as follows. Vector x is returned in the right singular vector
261:    (argument v) and has length equal to the number of columns of A and B.
262:    The other two vectors are returned stacked on top of each other [u;v] in
263:    the left singular vector argument, with length equal to m+n (number of rows
264:    of A plus number of rows of B).

266:    Level: beginner

268: .seealso: SVDSolve(), SVDGetConverged(), SVDSetWhichSingularTriplets()
269: @*/
270: PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)
271: {
272:   PetscInt       M,N;
273:   Vec            w;

277:   SVDCheckSolved(svd,1);
282:   if (sigma) *sigma = svd->sigma[svd->perm[i]];
283:   if (u || v) {
284:     if (!svd->isgeneralized) {
285:       MatGetSize(svd->OP,&M,&N);
286:       if (M<N) { w = u; u = v; v = w; }
287:     }
288:     SVDComputeVectors(svd);
289:     if (u) BVCopyVec(svd->U,svd->perm[i],u);
290:     if (v) BVCopyVec(svd->V,svd->perm[i],v);
291:   }
292:   PetscFunctionReturn(0);
293: }

295: /*
296:    SVDComputeResidualNorms_Standard - Computes the norms of the left and
297:    right residuals associated with the i-th computed singular triplet.

299:    Input Parameters:
300:      sigma - singular value
301:      u,v   - singular vectors
302:      x,y   - two work vectors with the same dimensions as u,v
303: */
304: static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)
305: {
306:   PetscInt       M,N;

308:   /* norm1 = ||A*v-sigma*u||_2 */
309:   if (norm1) {
310:     MatMult(svd->OP,v,x);
311:     VecAXPY(x,-sigma,u);
312:     VecNorm(x,NORM_2,norm1);
313:   }
314:   /* norm2 = ||A^T*u-sigma*v||_2 */
315:   if (norm2) {
316:     MatGetSize(svd->OP,&M,&N);
317:     if (M<N) MatMult(svd->A,u,y);
318:     else MatMult(svd->AT,u,y);
319:     VecAXPY(y,-sigma,v);
320:     VecNorm(y,NORM_2,norm2);
321:   }
322:   PetscFunctionReturn(0);
323: }

325: /*
326:    SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
327:    norm1 = ||s^2*A'*u-c*B'*B*x||_2 and norm2 = ||c^2*B'*v-s*A'*A*x||_2.

329:    Input Parameters:
330:      sigma - singular value
331:      uv    - left singular vectors [u;v]
332:      x     - right singular vector
333:      y,z   - two work vectors with the same dimension as x
334: */
335: static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
336: {
337:   Vec            u,v,ax,bx,nest,aux[2];
338:   PetscReal      c,s;

340:   MatCreateVecs(svd->OP,NULL,&u);
341:   MatCreateVecs(svd->OPb,NULL,&v);
342:   aux[0] = u;
343:   aux[1] = v;
344:   VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest);
345:   VecCopy(uv,nest);

347:   s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
348:   c = sigma*s;

350:   /* norm1 = ||s^2*A'*u-c*B'*B*x||_2 */
351:   if (norm1) {
352:     VecDuplicate(v,&bx);
353:     MatMultHermitianTranspose(svd->OP,u,z);
354:     MatMult(svd->OPb,x,bx);
355:     MatMultHermitianTranspose(svd->OPb,bx,y);
356:     VecAXPBY(y,s*s,-c,z);
357:     VecNorm(y,NORM_2,norm1);
358:     VecDestroy(&bx);
359:   }
360:   /* norm2 = ||c^2*B'*v-s*A'*A*x||_2 */
361:   if (norm2) {
362:     VecDuplicate(u,&ax);
363:     MatMultHermitianTranspose(svd->OPb,v,z);
364:     MatMult(svd->OP,x,ax);
365:     MatMultHermitianTranspose(svd->OP,ax,y);
366:     VecAXPBY(y,c*c,-s,z);
367:     VecNorm(y,NORM_2,norm2);
368:     VecDestroy(&ax);
369:   }

371:   VecDestroy(&v);
372:   VecDestroy(&u);
373:   VecDestroy(&nest);
374:   PetscFunctionReturn(0);
375: }

377: /*@
378:    SVDComputeError - Computes the error (based on the residual norm) associated
379:    with the i-th singular triplet.

381:    Collective on svd

383:    Input Parameters:
384: +  svd  - the singular value solver context
385: .  i    - the solution index
386: -  type - the type of error to compute

388:    Output Parameter:
389: .  error - the error

391:    Notes:
392:    The error can be computed in various ways, all of them based on the residual
393:    norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
394:    n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
395:    singular vector and v is the right singular vector.

397:    In the case of the GSVD, the two components of the residual norm are
398:    n1 = ||s^2*A'*u-c*B'*B*x||_2 and n2 = ||c^2*B'*v-s*A'*A*x||_2, where [u;v]
399:    are the left singular vectors and x is the right singular vector, with
400:    sigma=c/s.

402:    Level: beginner

404: .seealso: SVDErrorType, SVDSolve()
405: @*/
406: PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)
407: {
408:   PetscReal      sigma,norm1,norm2;
409:   Vec            u=NULL,v=NULL,x=NULL,y=NULL;

415:   SVDCheckSolved(svd,1);

417:   /* allocate work vectors */
418:   switch (svd->problem_type) {
419:     case SVD_STANDARD:
420:       SVDSetWorkVecs(svd,2,2);
421:       u = svd->workl[0];
422:       v = svd->workr[0];
423:       x = svd->workl[1];
424:       y = svd->workr[1];
425:       break;
426:     case SVD_GENERALIZED:
428:       SVDSetWorkVecs(svd,1,3);
429:       u = svd->workl[0];
430:       v = svd->workr[0];
431:       x = svd->workr[1];
432:       y = svd->workr[2];
433:       break;
434:   }

436:   /* compute residual norm and error */
437:   SVDGetSingularTriplet(svd,i,&sigma,u,v);
438:   switch (svd->problem_type) {
439:     case SVD_STANDARD:
440:       SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2);
441:       break;
442:     case SVD_GENERALIZED:
443:       SVDComputeResidualNorms_Generalized(svd,sigma,u,v,x,y,&norm1,&norm2);
444:       break;
445:   }
446:   *error = PetscSqrtReal(norm1*norm1+norm2*norm2);
447:   switch (type) {
448:     case SVD_ERROR_ABSOLUTE:
449:       break;
450:     case SVD_ERROR_RELATIVE:
451:       *error /= sigma;
452:       break;
453:     case SVD_ERROR_NORM:
454:       if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
455:       if (svd->isgeneralized && !svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
456:       *error /= PetscMax(svd->nrma,svd->nrmb);
457:       break;
458:     default:
459:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
460:   }
461:   PetscFunctionReturn(0);
462: }