| 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__
|
| 1281 |
slepc |
34 |
#define __FUNCT__ "cgs"
|
| 1298 |
slepc |
35 |
PetscErrorCode cgs(Vec v, int n, Vec *V, PetscReal *norm)
|
| 1278 |
slepc |
36 |
{
|
|
|
37 |
PetscErrorCode ierr;
|
|
|
38 |
Vec w;
|
| 1298 |
slepc |
39 |
PetscScalar *h,alpha;
|
|
|
40 |
PetscReal sum;
|
|
|
41 |
int j;
|
| 1278 |
slepc |
42 |
|
|
|
43 |
PetscFunctionBegin;
|
|
|
44 |
if (n>0) {
|
|
|
45 |
ierr = VecDuplicate(v,&w);CHKERRQ(ierr);
|
|
|
46 |
ierr = PetscMalloc(sizeof(PetscScalar)*n,&h);CHKERRQ(ierr);
|
|
|
47 |
|
| 1298 |
slepc |
48 |
if (norm) {
|
|
|
49 |
ierr = VecDotBegin(v,v,&alpha);CHKERRQ(ierr);
|
|
|
50 |
ierr = VecMDotBegin(v,n,V,h);CHKERRQ(ierr);
|
|
|
51 |
ierr = VecDotEnd(v,v,&alpha);CHKERRQ(ierr);
|
|
|
52 |
ierr = VecMDotEnd(v,n,V,h);CHKERRQ(ierr);
|
|
|
53 |
} else {
|
|
|
54 |
ierr = VecMDot(v,n,V,h);CHKERRQ(ierr);
|
|
|
55 |
}
|
| 1278 |
slepc |
56 |
ierr = VecSet(w,0.0);CHKERRQ(ierr);
|
|
|
57 |
ierr = VecMAXPY(w,n,h,V);CHKERRQ(ierr);
|
|
|
58 |
ierr = VecAXPY(v,-1.0,w);CHKERRQ(ierr);
|
|
|
59 |
|
| 1298 |
slepc |
60 |
if (norm) {
|
|
|
61 |
sum = 0.0;
|
|
|
62 |
for (j=0; j<n; j++)
|
|
|
63 |
sum += PetscRealPart(h[j] * PetscConj(h[j]));
|
|
|
64 |
*norm = PetscRealPart(alpha)-sum;
|
|
|
65 |
if (*norm < 0.0) {
|
|
|
66 |
ierr = VecNorm(v,NORM_2,norm);CHKERRQ(ierr);
|
|
|
67 |
} else *norm = sqrt(*norm);
|
|
|
68 |
}
|
|
|
69 |
|
| 1278 |
slepc |
70 |
ierr = VecDestroy(w);CHKERRQ(ierr);
|
|
|
71 |
ierr = PetscFree(h);CHKERRQ(ierr);
|
| 1298 |
slepc |
72 |
} else if (norm) {
|
|
|
73 |
ierr = VecNorm(v,NORM_2,norm);CHKERRQ(ierr);
|
| 1278 |
slepc |
74 |
}
|
|
|
75 |
PetscFunctionReturn(0);
|
|
|
76 |
}
|
|
|
77 |
|
|
|
78 |
#undef __FUNCT__
|
|
|
79 |
#define __FUNCT__ "SVDSolve_LANCZOS"
|
|
|
80 |
PetscErrorCode SVDSolve_LANCZOS(SVD svd)
|
|
|
81 |
{
|
|
|
82 |
PetscErrorCode ierr;
|
| 1298 |
slepc |
83 |
SVD_LANCZOS *lanczos = (SVD_LANCZOS *)svd->data;
|
| 1293 |
slepc |
84 |
PetscReal *alpha,*beta,norm,*work;
|
| 1278 |
slepc |
85 |
PetscScalar *Q,*PT;
|
| 1293 |
slepc |
86 |
PetscInt *perm;
|
| 1285 |
slepc |
87 |
int i,j,k,l,n,zero=0,info;
|
| 1278 |
slepc |
88 |
Vec *V,*U;
|
| 1293 |
slepc |
89 |
PetscTruth conv;
|
| 1278 |
slepc |
90 |
|
|
|
91 |
PetscFunctionBegin;
|
| 1293 |
slepc |
92 |
/* allocate working space */
|
| 1278 |
slepc |
93 |
ierr = PetscMalloc(sizeof(PetscReal)*svd->n,&alpha);CHKERRQ(ierr);
|
|
|
94 |
ierr = PetscMalloc(sizeof(PetscReal)*svd->n,&beta);CHKERRQ(ierr);
|
|
|
95 |
ierr = PetscMalloc(sizeof(PetscScalar)*svd->n*svd->n,&Q);CHKERRQ(ierr);
|
|
|
96 |
ierr = PetscMalloc(sizeof(PetscScalar)*svd->n*svd->n,&PT);CHKERRQ(ierr);
|
|
|
97 |
ierr = PetscMalloc(sizeof(PetscReal)*4*svd->n,&work);CHKERRQ(ierr);
|
|
|
98 |
ierr = VecDuplicateVecs(svd->V[0],svd->n+1,&V);CHKERRQ(ierr);
|
|
|
99 |
ierr = VecDuplicateVecs(svd->U[0],svd->n,&U);CHKERRQ(ierr);
|
|
|
100 |
|
| 1293 |
slepc |
101 |
/* normalize start vector */
|
| 1278 |
slepc |
102 |
ierr = VecCopy(svd->vec_initial,V[0]);CHKERRQ(ierr);
|
| 1293 |
slepc |
103 |
ierr = VecNormalize(V[0],&norm);CHKERRQ(ierr);
|
| 1278 |
slepc |
104 |
|
| 1283 |
slepc |
105 |
while (svd->reason == SVD_CONVERGED_ITERATING) {
|
|
|
106 |
svd->its++;
|
|
|
107 |
|
| 1293 |
slepc |
108 |
/* inner loop */
|
| 1280 |
slepc |
109 |
for (i=svd->nconv;i<svd->n;i++) {
|
| 1278 |
slepc |
110 |
ierr = MatMult(svd->A,V[i],U[i]);CHKERRQ(ierr);
|
| 1298 |
slepc |
111 |
if (lanczos->oneside) {
|
|
|
112 |
if (i>svd->nconv) { ierr = VecAXPY(U[i],-beta[i-svd->nconv-1],U[i-1]);CHKERRQ(ierr); }
|
|
|
113 |
} else {
|
| 1302 |
slepc |
114 |
ierr = IPOrthogonalize(svd->ip,i,PETSC_NULL,U,U[i],PT,alpha+i-svd->nconv,PETSC_NULL);CHKERRQ(ierr);
|
| 1298 |
slepc |
115 |
ierr = VecScale(U[i],1.0/alpha[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
116 |
}
|
| 1278 |
slepc |
117 |
|
|
|
118 |
if (svd->AT) {
|
|
|
119 |
ierr = MatMult(svd->AT,U[i],V[i+1]);CHKERRQ(ierr);
|
|
|
120 |
} else {
|
|
|
121 |
ierr = MatMultTranspose(svd->A,U[i],V[i+1]);CHKERRQ(ierr);
|
|
|
122 |
}
|
| 1298 |
slepc |
123 |
if (lanczos->oneside) {
|
|
|
124 |
ierr = VecNormBegin(U[i],NORM_2,alpha+i-svd->nconv);CHKERRQ(ierr);
|
|
|
125 |
ierr = VecMDotBegin(V[i+1],i+1,V,PT);CHKERRQ(ierr);
|
|
|
126 |
ierr = VecNormEnd(U[i],NORM_2,alpha+i-svd->nconv);CHKERRQ(ierr);
|
|
|
127 |
ierr = VecMDotEnd(V[i+1],i+1,V,PT);CHKERRQ(ierr);
|
|
|
128 |
|
|
|
129 |
ierr = VecScale(U[i],1.0/alpha[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
130 |
ierr = VecScale(V[i+1],1.0/alpha[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
131 |
for (j=0;j<=i;j++) PT[j] = - PT[j] / alpha[i-svd->nconv];
|
|
|
132 |
ierr = VecMAXPY(V[i+1],i+1,PT,V);CHKERRQ(ierr);
|
|
|
133 |
|
| 1302 |
slepc |
134 |
ierr = IPOrthogonalizeGS(svd->ip,i+1,PETSC_NULL,V,V[i+1],PT,PETSC_NULL,beta+i-svd->nconv,svd->V[0]);CHKERRQ(ierr);
|
| 1298 |
slepc |
135 |
ierr = VecScale(V[i+1],1.0/beta[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
136 |
} else {
|
| 1302 |
slepc |
137 |
ierr = IPOrthogonalize(svd->ip,i+1,PETSC_NULL,V,V[i+1],PT,beta+i-svd->nconv,PETSC_NULL);CHKERRQ(ierr);
|
| 1298 |
slepc |
138 |
ierr = VecScale(V[i+1],1.0/beta[i-svd->nconv]);CHKERRQ(ierr);
|
|
|
139 |
}
|
| 1278 |
slepc |
140 |
}
|
|
|
141 |
|
| 1293 |
slepc |
142 |
/* compute SVD of bidiagonal matrix */
|
| 1281 |
slepc |
143 |
n = svd->n - svd->nconv;
|
| 1278 |
slepc |
144 |
ierr = PetscMemzero(PT,sizeof(PetscScalar)*n*n);CHKERRQ(ierr);
|
|
|
145 |
ierr = PetscMemzero(Q,sizeof(PetscScalar)*n*n);CHKERRQ(ierr);
|
|
|
146 |
for (i=0;i<n;i++)
|
|
|
147 |
PT[i*n+i] = Q[i*n+i] = 1.0;
|
| 1285 |
slepc |
148 |
LAPACKbdsqr_("U",&n,&n,&n,&zero,alpha,beta,PT,&n,Q,&n,PETSC_NULL,&n,work,&info,1);
|
| 1278 |
slepc |
149 |
|
| 1293 |
slepc |
150 |
/* compute error estimates and converged singular vectors */
|
| 1278 |
slepc |
151 |
k = svd->nconv;
|
| 1293 |
slepc |
152 |
conv = PETSC_TRUE;
|
| 1280 |
slepc |
153 |
for (i=svd->nconv;i<svd->n;i++) {
|
| 1285 |
slepc |
154 |
if (svd->which == SVD_SMALLEST) j = n-i+svd->nconv-1;
|
|
|
155 |
else j = i-svd->nconv;
|
|
|
156 |
svd->sigma[i] = alpha[j];
|
| 1293 |
slepc |
157 |
svd->errest[i] = PetscAbsReal(Q[j*n+n-1])*beta[n-1];
|
|
|
158 |
if (conv) {
|
|
|
159 |
if (svd->errest[i] < svd->tol) {
|
|
|
160 |
ierr = VecSet(svd->V[i],0.0);CHKERRQ(ierr);
|
|
|
161 |
for (l=0;l<n;l++) {
|
|
|
162 |
ierr = VecAXPY(svd->V[i],PT[l*n+j],V[l+svd->nconv]);CHKERRQ(ierr);
|
|
|
163 |
}
|
|
|
164 |
ierr = VecSet(svd->U[i],0.0);CHKERRQ(ierr);
|
|
|
165 |
ierr = VecMAXPY(svd->U[i],n,Q+j*n,U+svd->nconv);CHKERRQ(ierr);
|
|
|
166 |
k++;
|
|
|
167 |
} else conv = PETSC_FALSE;
|
| 1278 |
slepc |
168 |
}
|
|
|
169 |
}
|
| 1293 |
slepc |
170 |
|
|
|
171 |
if (svd->its > svd->max_it) svd->reason = SVD_DIVERGED_ITS;
|
| 1295 |
slepc |
172 |
if (k >= svd->nsv) svd->reason = SVD_CONVERGED_TOL;
|
| 1293 |
slepc |
173 |
if (svd->reason == SVD_CONVERGED_ITERATING) {
|
|
|
174 |
/* compute restart vector */
|
|
|
175 |
if (svd->which == SVD_SMALLEST) j = n-k+svd->nconv-1;
|
|
|
176 |
else j = k-svd->nconv;
|
|
|
177 |
ierr = VecSet(svd->V[k],0.0);CHKERRQ(ierr);
|
|
|
178 |
for (l=0;l<n;l++) {
|
|
|
179 |
ierr = VecAXPY(svd->V[k],PT[l*n+j],V[l+svd->nconv]);CHKERRQ(ierr);
|
|
|
180 |
}
|
|
|
181 |
ierr = VecCopy(svd->V[k],V[k]);CHKERRQ(ierr);
|
|
|
182 |
}
|
|
|
183 |
|
|
|
184 |
/* copy converged singular vectors from temporary space */
|
| 1281 |
slepc |
185 |
for (i=svd->nconv;i<k;i++) {
|
|
|
186 |
ierr = VecCopy(svd->V[i],V[i]);CHKERRQ(ierr);
|
|
|
187 |
ierr = VecCopy(svd->U[i],U[i]);CHKERRQ(ierr);
|
|
|
188 |
}
|
| 1278 |
slepc |
189 |
svd->nconv = k;
|
| 1293 |
slepc |
190 |
|
|
|
191 |
SVDMonitor(svd,svd->its,svd->nconv,svd->sigma,svd->errest,svd->n);
|
| 1278 |
slepc |
192 |
}
|
|
|
193 |
|
| 1293 |
slepc |
194 |
/* sort singular triplets */
|
|
|
195 |
ierr = PetscMalloc(sizeof(PetscInt)*svd->nconv,&perm);CHKERRQ(ierr);
|
|
|
196 |
for (i=0;i<svd->nconv;i++) {
|
|
|
197 |
alpha[i] = svd->sigma[i];
|
|
|
198 |
beta[i] = svd->errest[i];
|
|
|
199 |
perm[i] = i;
|
|
|
200 |
}
|
|
|
201 |
ierr = PetscSortRealWithPermutation(svd->nconv,svd->sigma,perm);CHKERRQ(ierr);
|
|
|
202 |
for (i=0;i<svd->nconv;i++) {
|
|
|
203 |
if (svd->which == SVD_SMALLEST) j = perm[i];
|
|
|
204 |
else j = perm[svd->nconv-i-1];
|
|
|
205 |
svd->sigma[i] = alpha[j];
|
|
|
206 |
svd->errest[i] = beta[j];
|
|
|
207 |
ierr = VecCopy(V[j],svd->V[i]);CHKERRQ(ierr);
|
|
|
208 |
ierr = VecCopy(U[j],svd->U[i]);CHKERRQ(ierr);
|
|
|
209 |
}
|
|
|
210 |
|
|
|
211 |
/* free working space */
|
| 1302 |
slepc |
212 |
ierr = VecDestroyVecs(V,svd->n+1);CHKERRQ(ierr);
|
|
|
213 |
ierr = VecDestroyVecs(U,svd->n);CHKERRQ(ierr);
|
| 1278 |
slepc |
214 |
|
|
|
215 |
ierr = PetscFree(alpha);CHKERRQ(ierr);
|
|
|
216 |
ierr = PetscFree(beta);CHKERRQ(ierr);
|
|
|
217 |
ierr = PetscFree(Q);CHKERRQ(ierr);
|
|
|
218 |
ierr = PetscFree(PT);CHKERRQ(ierr);
|
|
|
219 |
ierr = PetscFree(work);CHKERRQ(ierr);
|
| 1293 |
slepc |
220 |
ierr = PetscFree(perm);CHKERRQ(ierr);
|
| 1278 |
slepc |
221 |
PetscFunctionReturn(0);
|
|
|
222 |
}
|
|
|
223 |
|
| 1298 |
slepc |
224 |
#undef __FUNCT__
|
|
|
225 |
#define __FUNCT__ "SVDSetFromOptions_LANCZOS"
|
|
|
226 |
PetscErrorCode SVDSetFromOptions_LANCZOS(SVD svd)
|
|
|
227 |
{
|
|
|
228 |
PetscErrorCode ierr;
|
|
|
229 |
SVD_LANCZOS *lanczos = (SVD_LANCZOS *)svd->data;
|
|
|
230 |
|
|
|
231 |
PetscFunctionBegin;
|
|
|
232 |
ierr = PetscOptionsBegin(svd->comm,svd->prefix,"LANCZOS Singular Value Solver Options","SVD");CHKERRQ(ierr);
|
|
|
233 |
ierr = PetscOptionsName("-svd_lanczos_oneside","Lanczos one-side reorthogonalization","SVDLanczosSetOneSideReorthogonalization",&lanczos->oneside);CHKERRQ(ierr);
|
|
|
234 |
ierr = PetscOptionsEnd();CHKERRQ(ierr);
|
|
|
235 |
PetscFunctionReturn(0);
|
|
|
236 |
}
|
| 1278 |
slepc |
237 |
EXTERN_C_BEGIN
|
| 1298 |
slepc |
238 |
|
| 1278 |
slepc |
239 |
#undef __FUNCT__
|
| 1298 |
slepc |
240 |
#define __FUNCT__ "SVDLanczosSetOneSideReorthogonalization_LANCZOS"
|
|
|
241 |
PetscErrorCode SVDLanczosSetOneSideReorthogonalization_LANCZOS(SVD svd,PetscTruth oneside)
|
|
|
242 |
{
|
|
|
243 |
SVD_LANCZOS *lanczos = (SVD_LANCZOS *)svd->data;
|
|
|
244 |
|
|
|
245 |
PetscFunctionBegin;
|
|
|
246 |
lanczos->oneside = oneside;
|
|
|
247 |
PetscFunctionReturn(0);
|
|
|
248 |
}
|
|
|
249 |
EXTERN_C_BEGIN
|
|
|
250 |
|
|
|
251 |
#undef __FUNCT__
|
|
|
252 |
#define __FUNCT__ "SVDLanczosSetOneSideReorthogonalization"
|
|
|
253 |
PetscErrorCode SVDLanczosSetOneSideReorthogonalization(SVD svd,PetscTruth oneside)
|
|
|
254 |
{
|
|
|
255 |
PetscErrorCode ierr, (*f)(SVD,PetscTruth);
|
|
|
256 |
|
|
|
257 |
PetscFunctionBegin;
|
|
|
258 |
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
|
|
|
259 |
ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDLanczosSetOneSideReorthogonalization_C",(void (**)())&f);CHKERRQ(ierr);
|
|
|
260 |
if (f) {
|
|
|
261 |
ierr = (*f)(svd,oneside);CHKERRQ(ierr);
|
|
|
262 |
}
|
|
|
263 |
PetscFunctionReturn(0);
|
|
|
264 |
}
|
|
|
265 |
|
|
|
266 |
#undef __FUNCT__
|
|
|
267 |
#define __FUNCT__ "SVDView_LANCZOS"
|
|
|
268 |
PetscErrorCode SVDView_LANCZOS(SVD svd,PetscViewer viewer)
|
|
|
269 |
{
|
|
|
270 |
PetscErrorCode ierr;
|
|
|
271 |
SVD_LANCZOS *lanczos = (SVD_LANCZOS *)svd->data;
|
|
|
272 |
|
|
|
273 |
PetscFunctionBegin;
|
|
|
274 |
ierr = PetscViewerASCIIPrintf(viewer,"Lanczos reorthogonalization: %s\n",lanczos->oneside ? "one-side" : "two-side");CHKERRQ(ierr);
|
|
|
275 |
PetscFunctionReturn(0);
|
|
|
276 |
}
|
|
|
277 |
|
|
|
278 |
EXTERN_C_BEGIN
|
|
|
279 |
#undef __FUNCT__
|
| 1278 |
slepc |
280 |
#define __FUNCT__ "SVDCreate_LANCZOS"
|
|
|
281 |
PetscErrorCode SVDCreate_LANCZOS(SVD svd)
|
|
|
282 |
{
|
| 1298 |
slepc |
283 |
PetscErrorCode ierr;
|
|
|
284 |
SVD_LANCZOS *lanczos;
|
|
|
285 |
|
| 1278 |
slepc |
286 |
PetscFunctionBegin;
|
| 1298 |
slepc |
287 |
ierr = PetscNew(SVD_LANCZOS,&lanczos);CHKERRQ(ierr);
|
|
|
288 |
PetscLogObjectMemory(svd,sizeof(SVD_LANCZOS));
|
|
|
289 |
svd->data = (void *)lanczos;
|
|
|
290 |
svd->ops->setup = SVDSetUp_LANCZOS;
|
|
|
291 |
svd->ops->solve = SVDSolve_LANCZOS;
|
|
|
292 |
svd->ops->setfromoptions = SVDSetFromOptions_LANCZOS;
|
|
|
293 |
svd->ops->view = SVDView_LANCZOS;
|
|
|
294 |
lanczos->oneside = PETSC_FALSE;
|
|
|
295 |
ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDLanczosSetOneSideReorthogonalization_C","SVDLanczosSetOneSideReorthogonalization_LANCZOS",SVDLanczosSetOneSideReorthogonalization_LANCZOS);CHKERRQ(ierr);
|
| 1278 |
slepc |
296 |
PetscFunctionReturn(0);
|
|
|
297 |
}
|
|
|
298 |
EXTERN_C_END
|