Actual source code: svdsolve.c

slepc-3.18.1 2022-11-02
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,omega2,u,v,w;
 23:   PetscInt           i,n,N,oldsize;
 24:   MatType            Atype;
 25:   VecType            vtype;
 26:   Mat                Omega;
 27:   const PetscScalar* varray;

 29:   if (!svd->leftbasis) {
 30:     /* generate left singular vectors on U */
 31:     if (!svd->U) SVDGetBV(svd,NULL,&svd->U);
 32:     BVGetSizes(svd->U,NULL,NULL,&oldsize);
 33:     if (!oldsize) {
 34:       if (!((PetscObject)(svd->U))->type_name) BVSetType(svd->U,((PetscObject)(svd->V))->type_name);
 35:       MatCreateVecsEmpty(svd->A,NULL,&tl);
 36:       BVSetSizesFromVec(svd->U,tl,svd->ncv);
 37:       VecDestroy(&tl);
 38:     }
 39:     BVSetActiveColumns(svd->V,0,svd->nconv);
 40:     BVSetActiveColumns(svd->U,0,svd->nconv);
 41:     if (!svd->ishyperbolic) BVMatMult(svd->V,svd->A,svd->U);
 42:     else if (svd->swapped) {  /* compute right singular vectors as V=A'*Omega*U */
 43:       MatCreateVecs(svd->A,&w,NULL);
 44:       for (i=0;i<svd->nconv;i++) {
 45:         BVGetColumn(svd->V,i,&v);
 46:         BVGetColumn(svd->U,i,&u);
 47:         VecPointwiseMult(w,v,svd->omega);
 48:         MatMult(svd->A,w,u);
 49:         BVRestoreColumn(svd->V,i,&v);
 50:         BVRestoreColumn(svd->U,i,&u);
 51:       }
 52:       VecDestroy(&w);
 53:     } else {  /* compute left singular vectors as usual U=A*V, and set-up Omega-orthogonalization of U */
 54:       MatGetType(svd->A,&Atype);
 55:       BVGetSizes(svd->U,&n,&N,NULL);
 56:       MatCreate(PetscObjectComm((PetscObject)svd),&Omega);
 57:       MatSetSizes(Omega,n,n,N,N);
 58:       MatSetType(Omega,Atype);
 59:       MatSetUp(Omega);
 60:       MatDiagonalSet(Omega,svd->omega,INSERT_VALUES);
 61:       BVMatMult(svd->V,svd->A,svd->U);
 62:       BVSetMatrix(svd->U,Omega,PETSC_TRUE);
 63:       MatDestroy(&Omega);
 64:     }
 65:     BVOrthogonalize(svd->U,NULL);
 66:     if (svd->ishyperbolic && !svd->swapped) {  /* store signature after Omega-orthogonalization */
 67:       MatGetVecType(svd->A,&vtype);
 68:       VecCreate(PETSC_COMM_SELF,&omega2);
 69:       VecSetSizes(omega2,svd->nconv,svd->nconv);
 70:       VecSetType(omega2,vtype);
 71:       BVGetSignature(svd->U,omega2);
 72:       VecGetArrayRead(omega2,&varray);
 73:       for (i=0;i<svd->nconv;i++) {
 74:         svd->sign[i] = PetscRealPart(varray[i]);
 75:         if (PetscRealPart(varray[i])<0.0) BVScaleColumn(svd->U,i,-1.0);
 76:       }
 77:       VecRestoreArrayRead(omega2,&varray);
 78:       VecDestroy(&omega2);
 79:     }
 80:   }
 81:   return 0;
 82: }

 84: PetscErrorCode SVDComputeVectors(SVD svd)
 85: {
 86:   SVDCheckSolved(svd,1);
 87:   if (svd->state==SVD_STATE_SOLVED) PetscTryTypeMethod(svd,computevectors);
 88:   svd->state = SVD_STATE_VECTORS;
 89:   return 0;
 90: }

 92: /*@
 93:    SVDSolve - Solves the singular value problem.

 95:    Collective on svd

 97:    Input Parameter:
 98: .  svd - singular value solver context obtained from SVDCreate()

100:    Options Database Keys:
101: +  -svd_view - print information about the solver used
102: .  -svd_view_mat0 - view the first matrix (A)
103: .  -svd_view_mat1 - view the second matrix (B)
104: .  -svd_view_signature - view the signature matrix (omega)
105: .  -svd_view_vectors - view the computed singular vectors
106: .  -svd_view_values - view the computed singular values
107: .  -svd_converged_reason - print reason for convergence, and number of iterations
108: .  -svd_error_absolute - print absolute errors of each singular triplet
109: .  -svd_error_relative - print relative errors of each singular triplet
110: -  -svd_error_norm     - print errors relative to the matrix norms of each singular triplet

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

119:    Level: beginner

121: .seealso: SVDCreate(), SVDSetUp(), SVDDestroy()
122: @*/
123: PetscErrorCode SVDSolve(SVD svd)
124: {
125:   PetscInt       i,*workperm;

128:   if (svd->state>=SVD_STATE_SOLVED) return 0;
129:   PetscLogEventBegin(SVD_Solve,svd,0,0,0);

131:   /* call setup */
132:   SVDSetUp(svd);
133:   svd->its = 0;
134:   svd->nconv = 0;
135:   for (i=0;i<svd->ncv;i++) {
136:     svd->sigma[i]  = 0.0;
137:     svd->errest[i] = 0.0;
138:     svd->perm[i]   = i;
139:   }
140:   SVDViewFromOptions(svd,NULL,"-svd_view_pre");

142:   switch (svd->problem_type) {
143:     case SVD_STANDARD:
144:       PetscUseTypeMethod(svd,solve);
145:       break;
146:     case SVD_GENERALIZED:
147:       PetscUseTypeMethod(svd,solveg);
148:       break;
149:     case SVD_HYPERBOLIC:
150:       PetscUseTypeMethod(svd,solveh);
151:       break;
152:   }
153:   svd->state = SVD_STATE_SOLVED;

155:   /* sort singular triplets */
156:   if (svd->which == SVD_SMALLEST) PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm);
157:   else {
158:     PetscMalloc1(svd->nconv,&workperm);
159:     for (i=0;i<svd->nconv;i++) workperm[i] = i;
160:     PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm);
161:     for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
162:     PetscFree(workperm);
163:   }
164:   PetscLogEventEnd(SVD_Solve,svd,0,0,0);

166:   /* various viewers */
167:   SVDViewFromOptions(svd,NULL,"-svd_view");
168:   SVDConvergedReasonViewFromOptions(svd);
169:   SVDErrorViewFromOptions(svd);
170:   SVDValuesViewFromOptions(svd);
171:   SVDVectorsViewFromOptions(svd);
172:   MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0");
173:   if (svd->isgeneralized) MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1");
174:   if (svd->ishyperbolic) VecViewFromOptions(svd->omega,(PetscObject)svd,"-svd_view_signature");

176:   /* Remove the initial subspaces */
177:   svd->nini = 0;
178:   svd->ninil = 0;
179:   return 0;
180: }

182: /*@
183:    SVDGetIterationNumber - Gets the current iteration number. If the
184:    call to SVDSolve() is complete, then it returns the number of iterations
185:    carried out by the solution method.

187:    Not Collective

189:    Input Parameter:
190: .  svd - the singular value solver context

192:    Output Parameter:
193: .  its - number of iterations

195:    Note:
196:    During the i-th iteration this call returns i-1. If SVDSolve() is
197:    complete, then parameter "its" contains either the iteration number at
198:    which convergence was successfully reached, or failure was detected.
199:    Call SVDGetConvergedReason() to determine if the solver converged or
200:    failed and why.

202:    Level: intermediate

204: .seealso: SVDGetConvergedReason(), SVDSetTolerances()
205: @*/
206: PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)
207: {
210:   *its = svd->its;
211:   return 0;
212: }

214: /*@
215:    SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
216:    stopped.

218:    Not Collective

220:    Input Parameter:
221: .  svd - the singular value solver context

223:    Output Parameter:
224: .  reason - negative value indicates diverged, positive value converged
225:    (see SVDConvergedReason)

227:    Options Database Key:
228: .  -svd_converged_reason - print the reason to a viewer

230:    Notes:
231:    Possible values for reason are
232: +  SVD_CONVERGED_TOL - converged up to tolerance
233: .  SVD_CONVERGED_USER - converged due to a user-defined condition
234: .  SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
235: .  SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
236: -  SVD_DIVERGED_BREAKDOWN - generic breakdown in method

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

240:    Level: intermediate

242: .seealso: SVDSetTolerances(), SVDSolve(), SVDConvergedReason
243: @*/
244: PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
245: {
248:   SVDCheckSolved(svd,1);
249:   *reason = svd->reason;
250:   return 0;
251: }

253: /*@
254:    SVDGetConverged - Gets the number of converged singular values.

256:    Not Collective

258:    Input Parameter:
259: .  svd - the singular value solver context

261:    Output Parameter:
262: .  nconv - number of converged singular values

264:    Note:
265:    This function should be called after SVDSolve() has finished.

267:    Level: beginner

269: .seealso: SVDSetDimensions(), SVDSolve(), SVDGetSingularTriplet()
270: @*/
271: PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)
272: {
275:   SVDCheckSolved(svd,1);
276:   *nconv = svd->nconv;
277:   return 0;
278: }

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

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

287:    Input Parameters:
288: +  svd - singular value solver context
289: -  i   - index of the solution

291:    Output Parameters:
292: +  sigma - singular value
293: .  u     - left singular vector
294: -  v     - right singular vector

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

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

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

312:    Level: beginner

314: .seealso: SVDSolve(), SVDGetConverged(), SVDSetWhichSingularTriplets()
315: @*/
316: PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)
317: {
318:   PetscInt       M,N;
319:   Vec            w;

323:   SVDCheckSolved(svd,1);
328:   if (sigma) *sigma = svd->sigma[svd->perm[i]];
329:   if (u || v) {
330:     if (!svd->isgeneralized) {
331:       MatGetSize(svd->OP,&M,&N);
332:       if (M<N) { w = u; u = v; v = w; }
333:     }
334:     SVDComputeVectors(svd);
335:     if (u) BVCopyVec(svd->U,svd->perm[i],u);
336:     if (v) BVCopyVec(svd->V,svd->perm[i],v);
337:   }
338:   return 0;
339: }

341: /*
342:    SVDComputeResidualNorms_Standard - Computes the norms of the left and
343:    right residuals associated with the i-th computed singular triplet.

345:    Input Parameters:
346:      sigma - singular value
347:      u,v   - singular vectors
348:      x,y   - two work vectors with the same dimensions as u,v
349: */
350: static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)
351: {
352:   PetscInt       M,N;

354:   /* norm1 = ||A*v-sigma*u||_2 */
355:   if (norm1) {
356:     MatMult(svd->OP,v,x);
357:     VecAXPY(x,-sigma,u);
358:     VecNorm(x,NORM_2,norm1);
359:   }
360:   /* norm2 = ||A^T*u-sigma*v||_2 */
361:   if (norm2) {
362:     MatGetSize(svd->OP,&M,&N);
363:     if (M<N) MatMult(svd->A,u,y);
364:     else MatMult(svd->AT,u,y);
365:     VecAXPY(y,-sigma,v);
366:     VecNorm(y,NORM_2,norm2);
367:   }
368:   return 0;
369: }

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.

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: static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
382: {
383:   Vec            u,v,ax,bx,nest,aux[2];
384:   PetscReal      c,s;

386:   MatCreateVecs(svd->OP,NULL,&u);
387:   MatCreateVecs(svd->OPb,NULL,&v);
388:   aux[0] = u;
389:   aux[1] = v;
390:   VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest);
391:   VecCopy(uv,nest);

393:   s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
394:   c = sigma*s;

396:   /* norm1 = ||s^2*A'*u-c*B'*B*x||_2 */
397:   if (norm1) {
398:     VecDuplicate(v,&bx);
399:     MatMultHermitianTranspose(svd->OP,u,z);
400:     MatMult(svd->OPb,x,bx);
401:     MatMultHermitianTranspose(svd->OPb,bx,y);
402:     VecAXPBY(y,s*s,-c,z);
403:     VecNorm(y,NORM_2,norm1);
404:     VecDestroy(&bx);
405:   }
406:   /* norm2 = ||c^2*B'*v-s*A'*A*x||_2 */
407:   if (norm2) {
408:     VecDuplicate(u,&ax);
409:     MatMultHermitianTranspose(svd->OPb,v,z);
410:     MatMult(svd->OP,x,ax);
411:     MatMultHermitianTranspose(svd->OP,ax,y);
412:     VecAXPBY(y,c*c,-s,z);
413:     VecNorm(y,NORM_2,norm2);
414:     VecDestroy(&ax);
415:   }

417:   VecDestroy(&v);
418:   VecDestroy(&u);
419:   VecDestroy(&nest);
420:   return 0;
421: }

423: /*
424:    SVDComputeResidualNorms_Hyperbolic - Computes the norms of the left and
425:    right residuals associated with the i-th computed singular triplet.

427:    Input Parameters:
428:      sigma - singular value
429:      sign  - corresponding element of the signature Omega2
430:      u,v   - singular vectors
431:      x,y,z - three work vectors with the same dimensions as u,v,u
432: */
433: static PetscErrorCode SVDComputeResidualNorms_Hyperbolic(SVD svd,PetscReal sigma,PetscReal sign,Vec u,Vec v,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
434: {
435:   PetscInt       M,N;

437:   /* norm1 = ||A*v-sigma*u||_2 */
438:   if (norm1) {
439:     MatMult(svd->OP,v,x);
440:     VecAXPY(x,-sigma,u);
441:     VecNorm(x,NORM_2,norm1);
442:   }
443:   /* norm2 = ||A^T*Omega*u-sigma*sign*v||_2 */
444:   if (norm2) {
445:     MatGetSize(svd->OP,&M,&N);
446:     VecPointwiseMult(z,u,svd->omega);
447:     if (M<N) MatMult(svd->A,z,y);
448:     else MatMult(svd->AT,z,y);
449:     VecAXPY(y,-sigma*sign,v);
450:     VecNorm(y,NORM_2,norm2);
451:   }
452:   return 0;
453: }

455: /*@
456:    SVDComputeError - Computes the error (based on the residual norm) associated
457:    with the i-th singular triplet.

459:    Collective on svd

461:    Input Parameters:
462: +  svd  - the singular value solver context
463: .  i    - the solution index
464: -  type - the type of error to compute

466:    Output Parameter:
467: .  error - the error

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

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

480:    Level: beginner

482: .seealso: SVDErrorType, SVDSolve()
483: @*/
484: PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)
485: {
486:   PetscReal      sigma,norm1,norm2;
487:   Vec            u=NULL,v=NULL,x=NULL,y=NULL,z=NULL;

493:   SVDCheckSolved(svd,1);

495:   /* allocate work vectors */
496:   switch (svd->problem_type) {
497:     case SVD_STANDARD:
498:       SVDSetWorkVecs(svd,2,2);
499:       u = svd->workl[0];
500:       v = svd->workr[0];
501:       x = svd->workl[1];
502:       y = svd->workr[1];
503:       break;
504:     case SVD_GENERALIZED:
506:       SVDSetWorkVecs(svd,1,3);
507:       u = svd->workl[0];
508:       v = svd->workr[0];
509:       x = svd->workr[1];
510:       y = svd->workr[2];
511:       break;
512:     case SVD_HYPERBOLIC:
513:       SVDSetWorkVecs(svd,3,2);
514:       u = svd->workl[0];
515:       v = svd->workr[0];
516:       x = svd->workl[1];
517:       y = svd->workr[1];
518:       z = svd->workl[2];
519:       break;
520:   }

522:   /* compute residual norm and error */
523:   SVDGetSingularTriplet(svd,i,&sigma,u,v);
524:   switch (svd->problem_type) {
525:     case SVD_STANDARD:
526:       SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2);
527:       break;
528:     case SVD_GENERALIZED:
529:       SVDComputeResidualNorms_Generalized(svd,sigma,u,v,x,y,&norm1,&norm2);
530:       break;
531:     case SVD_HYPERBOLIC:
532:       SVDComputeResidualNorms_Hyperbolic(svd,sigma,svd->sign[svd->perm[i]],u,v,x,y,z,&norm1,&norm2);
533:       break;
534:   }
535:   *error = SlepcAbs(norm1,norm2);
536:   switch (type) {
537:     case SVD_ERROR_ABSOLUTE:
538:       break;
539:     case SVD_ERROR_RELATIVE:
540:       *error /= sigma;
541:       break;
542:     case SVD_ERROR_NORM:
543:       if (!svd->nrma) MatNorm(svd->OP,NORM_INFINITY,&svd->nrma);
544:       if (svd->isgeneralized && !svd->nrmb) MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb);
545:       *error /= PetscMax(svd->nrma,svd->nrmb);
546:       break;
547:     default:
548:       SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
549:   }
550:   return 0;
551: }