Actual source code: rsvd.c
slepc-3.21.1 2024-04-26
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: SLEPc singular value solver: "randomized"
13: Method: RSVD
15: Algorithm:
17: Randomized singular value decomposition.
19: References:
21: [1] N. Halko, P.-G. Martinsson, and J. A. Tropp, "Finding
22: structure with randomness: Probabilistic algorithms for
23: constructing approximate matrix decompositions", SIAM Rev.,
24: 53(2):217-288, 2011.
25: */
27: #include <slepc/private/svdimpl.h>
29: static PetscErrorCode SVDSetUp_Randomized(SVD svd)
30: {
31: PetscInt N;
33: PetscFunctionBegin;
34: SVDCheckStandard(svd);
35: SVDCheckDefinite(svd);
36: PetscCheck(svd->which==SVD_LARGEST,PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"This solver supports only largest singular values");
37: PetscCall(MatGetSize(svd->A,NULL,&N));
38: PetscCall(SVDSetDimensions_Default(svd));
39: PetscCheck(svd->ncv>=svd->nsv,PetscObjectComm((PetscObject)svd),PETSC_ERR_USER_INPUT,"The value of ncv must not be smaller than nsv");
40: if (svd->max_it==PETSC_DEFAULT) svd->max_it = PetscMax(N/svd->ncv,100);
41: svd->leftbasis = PETSC_TRUE;
42: svd->mpd = svd->ncv;
43: PetscCall(SVDAllocateSolution(svd,0));
44: PetscCall(DSSetType(svd->ds,DSSVD));
45: PetscCall(DSAllocate(svd->ds,svd->ncv));
46: PetscCall(SVDSetWorkVecs(svd,1,1));
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode SVDRandomizedResidualNorm(SVD svd,PetscInt i,PetscScalar sigma,PetscReal *res)
51: {
52: PetscReal norm1,norm2;
53: Vec u,v,wu,wv;
55: PetscFunctionBegin;
56: *res = 1.0;
57: if (svd->conv!=SVD_CONV_MAXIT) {
58: wu = svd->swapped? svd->workr[0]: svd->workl[0];
59: wv = svd->swapped? svd->workl[0]: svd->workr[0];
60: PetscCall(BVGetColumn(svd->V,i,&v));
61: PetscCall(BVGetColumn(svd->U,i,&u));
62: /* norm1 = ||A*v-sigma*u||_2 */
63: PetscCall(MatMult(svd->A,v,wu));
64: PetscCall(VecAXPY(wu,-sigma,u));
65: PetscCall(VecNorm(wu,NORM_2,&norm1));
66: /* norm2 = ||A^T*u-sigma*v||_2 */
67: PetscCall(MatMult(svd->AT,u,wv));
68: PetscCall(VecAXPY(wv,-sigma,v));
69: PetscCall(VecNorm(wv,NORM_2,&norm2));
70: PetscCall(BVRestoreColumn(svd->V,i,&v));
71: PetscCall(BVRestoreColumn(svd->U,i,&u));
72: *res = SlepcAbs(norm1,norm2);
73: }
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: /* If A is a virtual Hermitian transpose, then BVMatMult will fail if PRODUCT_AhB is not implemented */
78: static PetscErrorCode BlockMatMult(BV V,Mat A,BV Y,Mat AT)
79: {
80: PetscFunctionBegin;
81: if (!PetscDefined(USE_COMPLEX)) PetscCall(BVMatMult(V,A,Y));
82: else {
83: PetscBool flg=PETSC_FALSE;
84: PetscCall(PetscObjectTypeCompare((PetscObject)A,MATHERMITIANTRANSPOSEVIRTUAL,&flg));
85: if (flg) PetscCall(BVMatMultHermitianTranspose(V,AT,Y));
86: else PetscCall(BVMatMult(V,A,Y));
87: }
88: PetscFunctionReturn(PETSC_SUCCESS);
89: }
91: static PetscErrorCode SVDSolve_Randomized(SVD svd)
92: {
93: PetscScalar *w;
94: PetscReal res=1.0;
95: PetscInt i,k=0;
96: Mat A,U,V;
98: PetscFunctionBegin;
99: /* Form random matrix, G. Complete the initial basis with random vectors */
100: PetscCall(BVSetActiveColumns(svd->V,svd->nini,svd->ncv));
101: PetscCall(BVSetRandomNormal(svd->V));
102: PetscCall(PetscCalloc1(svd->ncv,&w));
104: /* Subspace Iteration */
105: do {
106: svd->its++;
107: PetscCall(BVSetActiveColumns(svd->V,svd->nconv,svd->ncv));
108: PetscCall(BVSetActiveColumns(svd->U,svd->nconv,svd->ncv));
109: /* Form AG */
110: PetscCall(BlockMatMult(svd->V,svd->A,svd->U,svd->AT));
111: /* Orthogonalization Q=qr(AG)*/
112: PetscCall(BVOrthogonalize(svd->U,NULL));
113: /* Form B^*= AQ */
114: PetscCall(BlockMatMult(svd->U,svd->AT,svd->V,svd->A));
116: PetscCall(DSSetDimensions(svd->ds,svd->ncv,svd->nconv,svd->ncv));
117: PetscCall(DSSVDSetDimensions(svd->ds,svd->ncv));
118: PetscCall(DSGetMat(svd->ds,DS_MAT_A,&A));
119: PetscCall(MatZeroEntries(A));
120: PetscCall(BVOrthogonalize(svd->V,A));
121: PetscCall(DSRestoreMat(svd->ds,DS_MAT_A,&A));
122: PetscCall(DSSetState(svd->ds,DS_STATE_RAW));
123: PetscCall(DSSolve(svd->ds,w,NULL));
124: PetscCall(DSSort(svd->ds,w,NULL,NULL,NULL,NULL));
125: PetscCall(DSSynchronize(svd->ds,w,NULL));
126: PetscCall(DSGetMat(svd->ds,DS_MAT_U,&U));
127: PetscCall(DSGetMat(svd->ds,DS_MAT_V,&V));
128: PetscCall(BVMultInPlace(svd->U,V,svd->nconv,svd->ncv));
129: PetscCall(BVMultInPlace(svd->V,U,svd->nconv,svd->ncv));
130: PetscCall(DSRestoreMat(svd->ds,DS_MAT_U,&U));
131: PetscCall(DSRestoreMat(svd->ds,DS_MAT_V,&V));
132: /* Check convergence */
133: k = 0;
134: for (i=svd->nconv;i<svd->ncv;i++) {
135: PetscCall(SVDRandomizedResidualNorm(svd,i,w[i],&res));
136: svd->sigma[i] = PetscRealPart(w[i]);
137: PetscCall((*svd->converged)(svd,svd->sigma[i],res,&svd->errest[i],svd->convergedctx));
138: if (svd->errest[i] < svd->tol) k++;
139: else break;
140: }
141: if (svd->conv == SVD_CONV_MAXIT && svd->its >= svd->max_it) {
142: k = svd->nsv;
143: for (i=0;i<svd->ncv;i++) svd->sigma[i] = PetscRealPart(w[i]);
144: }
145: PetscCall((*svd->stopping)(svd,svd->its,svd->max_it,svd->nconv+k,svd->nsv,&svd->reason,svd->stoppingctx));
146: svd->nconv += k;
147: PetscCall(SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,svd->ncv));
148: } while (svd->reason == SVD_CONVERGED_ITERATING);
149: PetscCall(PetscFree(w));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: SLEPC_EXTERN PetscErrorCode SVDCreate_Randomized(SVD svd)
154: {
155: PetscFunctionBegin;
156: svd->ops->setup = SVDSetUp_Randomized;
157: svd->ops->solve = SVDSolve_Randomized;
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }