| 1278 |
slepc |
1 |
/*
|
|
|
2 |
|
|
|
3 |
SLEPc singular value solver: "lanczos"
|
|
|
4 |
|
| 1281 |
slepc |
5 |
Method: Golub-Kahan-Lanczos bidiagonalization
|
| 1278 |
slepc |
6 |
|
|
|
7 |
Last update: Nov 2006
|
|
|
8 |
|
|
|
9 |
*/
|
|
|
10 |
#include "src/svd/svdimpl.h" /*I "slepcsvd.h" I*/
|
| 1283 |
slepc |
11 |
#include "slepcblaslapack.h"
|
| 1278 |
slepc |
12 |
|
| 1298 |
slepc |
13 |
typedef struct {
|
|
|
14 |
PetscTruth oneside;
|
|
|
15 |
} SVD_LANCZOS;
|
|
|
16 |
|
| 1278 |
slepc |
17 |
#undef __FUNCT__
|
|
|
18 |
#define __FUNCT__ "SVDSetUp_LANCZOS"
|
|
|
19 |
PetscErrorCode SVDSetUp_LANCZOS(SVD svd)
|
|
|
20 |
{
|
|
|
21 |
PetscErrorCode ierr;
|
|
|
22 |
PetscInt M,N;
|
|
|
23 |
|
|
|
24 |
PetscFunctionBegin;
|
|
|
25 |
ierr = MatGetSize(svd->A,&M,&N);CHKERRQ(ierr);
|
|
|
26 |
if (svd->ncv == PETSC_DECIDE)
|
| 1283 |
slepc |
27 |
svd->ncv = PetscMin(PetscMin(M,N),PetscMax(2*svd->nsv,10));
|
|
|
28 |
if (svd->max_it == PETSC_DECIDE)
|
|
|
29 |
svd->max_it = PetscMax(PetscMin(M,N)/svd->ncv,100);
|
| 1278 |
slepc |
30 |
PetscFunctionReturn(0);
|
|
|
31 |
}
|
|
|
32 |
|
|
|
33 |
#undef __FUNCT__
|
|
|
34 |
#define __FUNCT__ "SVDSolve_LANCZOS"
|
|
|
35 |
PetscErrorCode SVDSolve_LANCZOS(SVD svd)
|
|
|
36 |
{
|
|
|
37 |
PetscErrorCode ierr;
|
| 1298 |
slepc |
38 |
SVD_LANCZOS *lanczos = (SVD_LANCZOS *)svd->data;
|
| 1293 |
slepc |
39 |
PetscReal *alpha,*beta,norm,*work;
|
| 1278 |
slepc |
40 |
PetscScalar *Q,*PT;
|
| 1293 |
slepc |
41 |
PetscInt *perm;
|
| 1285 |
slepc |
42 |
int i,j,k,l,n,zero=0,info;
|
| 1278 |
slepc |
43 |
Vec *V,*U;
|
| 1293 |
slepc |
44 |
PetscTruth conv;
|
| 1278 |
slepc |
45 |
|
|
|
46 |
PetscFunctionBegin;
|
| 1293 |
slepc |
47 |
/* allocate working space */
|
| 1278 |
slepc |
48 |
ierr = PetscMalloc(sizeof(PetscReal)*svd->n,&alpha);CHKERRQ(ierr);
|
|
|
49 |
ierr = PetscMalloc(sizeof(PetscReal)*svd->n,&beta);CHKERRQ(ierr);
|
|
|
50 |
ierr = PetscMalloc(sizeof(PetscScalar)*svd->n*svd->n,&Q);CHKERRQ(ierr);
|
|
|
51 |
ierr = PetscMalloc(sizeof(PetscScalar)*svd->n*svd->n,&PT);CHKERRQ(ierr);
|
|
|
52 |
ierr = PetscMalloc(sizeof(PetscReal)*4*svd->n,&work);CHKERRQ(ierr);
|
|
|
53 |
ierr = VecDuplicateVecs(svd->V[0],svd->n+1,&V);CHKERRQ(ierr);
|
|
|
54 |
ierr = VecDuplicateVecs(svd->U[0],svd->n,&U);CHKERRQ(ierr);
|
|
|
55 |
|
| 1293 |
slepc |
56 |
/* normalize start vector */
|
| 1278 |
slepc |
57 |
ierr = VecCopy(svd->vec_initial,V[0]);CHKERRQ(ierr);
|
| 1293 |
slepc |
58 |
ierr = VecNormalize(V[0],&norm);CHKERRQ(ierr);
|
| 1278 |
slepc |
59 |
|
| 1283 |
slepc |
60 |
while (svd->reason == SVD_CONVERGED_ITERATING) {
|
|
|
61 |
svd->its++;
|
|
|
62 |
|
| 1293 |
slepc |
63 |
/* inner loop */
|
| 1280 |
slepc |
64 |
for (i=svd->nconv;i<svd->n;i++) {
|
| 1307 |
slepc |
65 |
svd->matvecs++;
|
| 1278 |
slepc |
66 |
ierr = MatMult(svd->A,V[i],U[i]);CHKERRQ(ierr);
|
| 1298 |
slepc |
67 |
if (lanczos->oneside) {
|
|
|
68 |
if (i>svd->nconv) { ierr = VecAXPY(U[i],-beta[i-svd->nconv-1],U[i-1]);CHKERRQ(ierr); }
|
|
|
69 |
} else {
|
| 1307 |
slepc |
70 |
svd->dots += i;
|
| 1302 |
slepc |
71 |
ierr = IPOrthogonalize(svd->ip,i,PETSC_NULL,U,U[i],PT,alpha+i-svd->nconv,PETSC_NULL);CHKERRQ(ierr);
|
| 1298 |
slepc |
72 |
ierr = VecScale(U[i],1.0/alpha[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
73 |
}
|
| 1278 |
slepc |
74 |
|
| 1307 |
slepc |
75 |
svd->matvecs++;
|
| 1278 |
slepc |
76 |
if (svd->AT) {
|
|
|
77 |
ierr = MatMult(svd->AT,U[i],V[i+1]);CHKERRQ(ierr);
|
|
|
78 |
} else {
|
|
|
79 |
ierr = MatMultTranspose(svd->A,U[i],V[i+1]);CHKERRQ(ierr);
|
|
|
80 |
}
|
| 1298 |
slepc |
81 |
if (lanczos->oneside) {
|
| 1307 |
slepc |
82 |
svd->dots += i+1;
|
| 1298 |
slepc |
83 |
ierr = VecNormBegin(U[i],NORM_2,alpha+i-svd->nconv);CHKERRQ(ierr);
|
|
|
84 |
ierr = VecMDotBegin(V[i+1],i+1,V,PT);CHKERRQ(ierr);
|
|
|
85 |
ierr = VecNormEnd(U[i],NORM_2,alpha+i-svd->nconv);CHKERRQ(ierr);
|
|
|
86 |
ierr = VecMDotEnd(V[i+1],i+1,V,PT);CHKERRQ(ierr);
|
|
|
87 |
|
|
|
88 |
ierr = VecScale(U[i],1.0/alpha[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
89 |
ierr = VecScale(V[i+1],1.0/alpha[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
90 |
for (j=0;j<=i;j++) PT[j] = - PT[j] / alpha[i-svd->nconv];
|
|
|
91 |
ierr = VecMAXPY(V[i+1],i+1,PT,V);CHKERRQ(ierr);
|
|
|
92 |
|
| 1307 |
slepc |
93 |
ierr = IPOrthogonalizeGS(svd->ip,i+1,PETSC_NULL,V,V[i+1],PT,PETSC_NULL,beta+i-svd->nconv);CHKERRQ(ierr);
|
| 1298 |
slepc |
94 |
ierr = VecScale(V[i+1],1.0/beta[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
95 |
} else {
|
| 1307 |
slepc |
96 |
svd->dots += i+1;
|
| 1302 |
slepc |
97 |
ierr = IPOrthogonalize(svd->ip,i+1,PETSC_NULL,V,V[i+1],PT,beta+i-svd->nconv,PETSC_NULL);CHKERRQ(ierr);
|
| 1298 |
slepc |
98 |
ierr = VecScale(V[i+1],1.0/beta[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
99 |
}
|
| 1278 |
slepc |
100 |
}
|
|
|
101 |
|
| 1293 |
slepc |
102 |
/* compute SVD of bidiagonal matrix */
|
| 1281 |
slepc |
103 |
n = svd->n - svd->nconv;
|
| 1278 |
slepc |
104 |
ierr = PetscMemzero(PT,sizeof(PetscScalar)*n*n);CHKERRQ(ierr);
|
|
|
105 |
ierr = PetscMemzero(Q,sizeof(PetscScalar)*n*n);CHKERRQ(ierr);
|
|
|
106 |
for (i=0;i<n;i++)
|
|
|
107 |
PT[i*n+i] = Q[i*n+i] = 1.0;
|
| 1285 |
slepc |
108 |
LAPACKbdsqr_("U",&n,&n,&n,&zero,alpha,beta,PT,&n,Q,&n,PETSC_NULL,&n,work,&info,1);
|
| 1278 |
slepc |
109 |
|
| 1293 |
slepc |
110 |
/* compute error estimates and converged singular vectors */
|
| 1278 |
slepc |
111 |
k = svd->nconv;
|
| 1293 |
slepc |
112 |
conv = PETSC_TRUE;
|
| 1280 |
slepc |
113 |
for (i=svd->nconv;i<svd->n;i++) {
|
| 1285 |
slepc |
114 |
if (svd->which == SVD_SMALLEST) j = n-i+svd->nconv-1;
|
|
|
115 |
else j = i-svd->nconv;
|
|
|
116 |
svd->sigma[i] = alpha[j];
|
| 1312 |
slepc |
117 |
svd->errest[i] = PetscAbsReal(Q[j*n+n-1])*beta[n-1] / alpha[j];
|
| 1293 |
slepc |
118 |
if (conv) {
|
|
|
119 |
if (svd->errest[i] < svd->tol) {
|
|
|
120 |
ierr = VecSet(svd->V[i],0.0);CHKERRQ(ierr);
|
|
|
121 |
for (l=0;l<n;l++) {
|
|
|
122 |
ierr = VecAXPY(svd->V[i],PT[l*n+j],V[l+svd->nconv]);CHKERRQ(ierr);
|
|
|
123 |
}
|
|
|
124 |
ierr = VecSet(svd->U[i],0.0);CHKERRQ(ierr);
|
|
|
125 |
ierr = VecMAXPY(svd->U[i],n,Q+j*n,U+svd->nconv);CHKERRQ(ierr);
|
|
|
126 |
k++;
|
|
|
127 |
} else conv = PETSC_FALSE;
|
| 1278 |
slepc |
128 |
}
|
|
|
129 |
}
|
| 1293 |
slepc |
130 |
|
|
|
131 |
if (svd->its > svd->max_it) svd->reason = SVD_DIVERGED_ITS;
|
| 1295 |
slepc |
132 |
if (k >= svd->nsv) svd->reason = SVD_CONVERGED_TOL;
|
| 1293 |
slepc |
133 |
if (svd->reason == SVD_CONVERGED_ITERATING) {
|
|
|
134 |
/* compute restart vector */
|
|
|
135 |
if (svd->which == SVD_SMALLEST) j = n-k+svd->nconv-1;
|
|
|
136 |
else j = k-svd->nconv;
|
|
|
137 |
ierr = VecSet(svd->V[k],0.0);CHKERRQ(ierr);
|
|
|
138 |
for (l=0;l<n;l++) {
|
|
|
139 |
ierr = VecAXPY(svd->V[k],PT[l*n+j],V[l+svd->nconv]);CHKERRQ(ierr);
|
|
|
140 |
}
|
|
|
141 |
ierr = VecCopy(svd->V[k],V[k]);CHKERRQ(ierr);
|
|
|
142 |
}
|
|
|
143 |
|
|
|
144 |
/* copy converged singular vectors from temporary space */
|
| 1281 |
slepc |
145 |
for (i=svd->nconv;i<k;i++) {
|
|
|
146 |
ierr = VecCopy(svd->V[i],V[i]);CHKERRQ(ierr);
|
|
|
147 |
ierr = VecCopy(svd->U[i],U[i]);CHKERRQ(ierr);
|
|
|
148 |
}
|
| 1278 |
slepc |
149 |
svd->nconv = k;
|
| 1293 |
slepc |
150 |
|
|
|
151 |
SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,svd->n);
|
| 1278 |
slepc |
152 |
}
|
|
|
153 |
|
| 1293 |
slepc |
154 |
/* sort singular triplets */
|
|
|
155 |
ierr = PetscMalloc(sizeof(PetscInt)*svd->nconv,&perm);CHKERRQ(ierr);
|
|
|
156 |
for (i=0;i<svd->nconv;i++) {
|
|
|
157 |
alpha[i] = svd->sigma[i];
|
|
|
158 |
beta[i] = svd->errest[i];
|
|
|
159 |
perm[i] = i;
|
|
|
160 |
}
|
|
|
161 |
ierr = PetscSortRealWithPermutation(svd->nconv,svd->sigma,perm);CHKERRQ(ierr);
|
|
|
162 |
for (i=0;i<svd->nconv;i++) {
|
|
|
163 |
if (svd->which == SVD_SMALLEST) j = perm[i];
|
|
|
164 |
else j = perm[svd->nconv-i-1];
|
|
|
165 |
svd->sigma[i] = alpha[j];
|
|
|
166 |
svd->errest[i] = beta[j];
|
|
|
167 |
ierr = VecCopy(V[j],svd->V[i]);CHKERRQ(ierr);
|
|
|
168 |
ierr = VecCopy(U[j],svd->U[i]);CHKERRQ(ierr);
|
|
|
169 |
}
|
|
|
170 |
|
|
|
171 |
/* free working space */
|
| 1302 |
slepc |
172 |
ierr = VecDestroyVecs(V,svd->n+1);CHKERRQ(ierr);
|
|
|
173 |
ierr = VecDestroyVecs(U,svd->n);CHKERRQ(ierr);
|
| 1278 |
slepc |
174 |
|
|
|
175 |
ierr = PetscFree(alpha);CHKERRQ(ierr);
|
|
|
176 |
ierr = PetscFree(beta);CHKERRQ(ierr);
|
|
|
177 |
ierr = PetscFree(Q);CHKERRQ(ierr);
|
|
|
178 |
ierr = PetscFree(PT);CHKERRQ(ierr);
|
|
|
179 |
ierr = PetscFree(work);CHKERRQ(ierr);
|
| 1293 |
slepc |
180 |
ierr = PetscFree(perm);CHKERRQ(ierr);
|
| 1278 |
slepc |
181 |
PetscFunctionReturn(0);
|
|
|
182 |
}
|
|
|
183 |
|
| 1298 |
slepc |
184 |
#undef __FUNCT__
|
|
|
185 |
#define __FUNCT__ "SVDSetFromOptions_LANCZOS"
|
|
|
186 |
PetscErrorCode SVDSetFromOptions_LANCZOS(SVD svd)
|
|
|
187 |
{
|
|
|
188 |
PetscErrorCode ierr;
|
|
|
189 |
SVD_LANCZOS *lanczos = (SVD_LANCZOS *)svd->data;
|
|
|
190 |
|
|
|
191 |
PetscFunctionBegin;
|
|
|
192 |
ierr = PetscOptionsBegin(svd->comm,svd->prefix,"LANCZOS Singular Value Solver Options","SVD");CHKERRQ(ierr);
|
| 1312 |
slepc |
193 |
ierr = PetscOptionsTruth("-svd_lanczos_oneside","Lanczos one-side reorthogonalization","SVDLanczosSetOneSideReorthogonalization",PETSC_FALSE,&lanczos->oneside,PETSC_NULL);CHKERRQ(ierr);
|
| 1298 |
slepc |
194 |
ierr = PetscOptionsEnd();CHKERRQ(ierr);
|
|
|
195 |
PetscFunctionReturn(0);
|
|
|
196 |
}
|
| 1278 |
slepc |
197 |
EXTERN_C_BEGIN
|
| 1298 |
slepc |
198 |
|
| 1278 |
slepc |
199 |
#undef __FUNCT__
|
| 1298 |
slepc |
200 |
#define __FUNCT__ "SVDLanczosSetOneSideReorthogonalization_LANCZOS"
|
|
|
201 |
PetscErrorCode SVDLanczosSetOneSideReorthogonalization_LANCZOS(SVD svd,PetscTruth oneside)
|
|
|
202 |
{
|
|
|
203 |
SVD_LANCZOS *lanczos = (SVD_LANCZOS *)svd->data;
|
|
|
204 |
|
|
|
205 |
PetscFunctionBegin;
|
|
|
206 |
lanczos->oneside = oneside;
|
|
|
207 |
PetscFunctionReturn(0);
|
|
|
208 |
}
|
|
|
209 |
EXTERN_C_BEGIN
|
|
|
210 |
|
|
|
211 |
#undef __FUNCT__
|
|
|
212 |
#define __FUNCT__ "SVDLanczosSetOneSideReorthogonalization"
|
|
|
213 |
PetscErrorCode SVDLanczosSetOneSideReorthogonalization(SVD svd,PetscTruth oneside)
|
|
|
214 |
{
|
|
|
215 |
PetscErrorCode ierr, (*f)(SVD,PetscTruth);
|
|
|
216 |
|
|
|
217 |
PetscFunctionBegin;
|
|
|
218 |
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
|
|
|
219 |
ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDLanczosSetOneSideReorthogonalization_C",(void (**)())&f);CHKERRQ(ierr);
|
|
|
220 |
if (f) {
|
|
|
221 |
ierr = (*f)(svd,oneside);CHKERRQ(ierr);
|
|
|
222 |
}
|
|
|
223 |
PetscFunctionReturn(0);
|
|
|
224 |
}
|
|
|
225 |
|
|
|
226 |
#undef __FUNCT__
|
|
|
227 |
#define __FUNCT__ "SVDView_LANCZOS"
|
|
|
228 |
PetscErrorCode SVDView_LANCZOS(SVD svd,PetscViewer viewer)
|
|
|
229 |
{
|
|
|
230 |
PetscErrorCode ierr;
|
|
|
231 |
SVD_LANCZOS *lanczos = (SVD_LANCZOS *)svd->data;
|
|
|
232 |
|
|
|
233 |
PetscFunctionBegin;
|
|
|
234 |
ierr = PetscViewerASCIIPrintf(viewer,"Lanczos reorthogonalization: %s\n",lanczos->oneside ? "one-side" : "two-side");CHKERRQ(ierr);
|
|
|
235 |
PetscFunctionReturn(0);
|
|
|
236 |
}
|
|
|
237 |
|
|
|
238 |
EXTERN_C_BEGIN
|
|
|
239 |
#undef __FUNCT__
|
| 1278 |
slepc |
240 |
#define __FUNCT__ "SVDCreate_LANCZOS"
|
|
|
241 |
PetscErrorCode SVDCreate_LANCZOS(SVD svd)
|
|
|
242 |
{
|
| 1298 |
slepc |
243 |
PetscErrorCode ierr;
|
|
|
244 |
SVD_LANCZOS *lanczos;
|
|
|
245 |
|
| 1278 |
slepc |
246 |
PetscFunctionBegin;
|
| 1298 |
slepc |
247 |
ierr = PetscNew(SVD_LANCZOS,&lanczos);CHKERRQ(ierr);
|
|
|
248 |
PetscLogObjectMemory(svd,sizeof(SVD_LANCZOS));
|
|
|
249 |
svd->data = (void *)lanczos;
|
|
|
250 |
svd->ops->setup = SVDSetUp_LANCZOS;
|
|
|
251 |
svd->ops->solve = SVDSolve_LANCZOS;
|
|
|
252 |
svd->ops->setfromoptions = SVDSetFromOptions_LANCZOS;
|
|
|
253 |
svd->ops->view = SVDView_LANCZOS;
|
|
|
254 |
lanczos->oneside = PETSC_FALSE;
|
|
|
255 |
ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDLanczosSetOneSideReorthogonalization_C","SVDLanczosSetOneSideReorthogonalization_LANCZOS",SVDLanczosSetOneSideReorthogonalization_LANCZOS);CHKERRQ(ierr);
|
| 1278 |
slepc |
256 |
PetscFunctionReturn(0);
|
|
|
257 |
}
|
|
|
258 |
EXTERN_C_END
|