Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
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: "cross"
12 :
13 : Method: Uses a Hermitian eigensolver for A^T*A
14 : */
15 :
16 : #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
17 :
18 : typedef struct {
19 : PetscBool explicitmatrix;
20 : EPS eps;
21 : PetscBool usereps;
22 : Mat C,D;
23 : } SVD_CROSS;
24 :
25 : typedef struct {
26 : Mat A,AT;
27 : Vec w,diag,omega;
28 : PetscBool swapped;
29 : } SVD_CROSS_SHELL;
30 :
31 9217 : static PetscErrorCode MatMult_Cross(Mat B,Vec x,Vec y)
32 : {
33 9217 : SVD_CROSS_SHELL *ctx;
34 :
35 9217 : PetscFunctionBegin;
36 9217 : PetscCall(MatShellGetContext(B,&ctx));
37 9217 : PetscCall(MatMult(ctx->A,x,ctx->w));
38 9217 : if (ctx->omega && !ctx->swapped) PetscCall(VecPointwiseMult(ctx->w,ctx->w,ctx->omega));
39 9217 : PetscCall(MatMult(ctx->AT,ctx->w,y));
40 9217 : PetscFunctionReturn(PETSC_SUCCESS);
41 : }
42 :
43 8 : static PetscErrorCode MatGetDiagonal_Cross(Mat B,Vec d)
44 : {
45 8 : SVD_CROSS_SHELL *ctx;
46 8 : PetscMPIInt len;
47 8 : PetscInt N,n,i,j,start,end,ncols;
48 8 : PetscScalar *work1,*work2,*diag;
49 8 : const PetscInt *cols;
50 8 : const PetscScalar *vals;
51 :
52 8 : PetscFunctionBegin;
53 8 : PetscCall(MatShellGetContext(B,&ctx));
54 8 : if (!ctx->diag) {
55 : /* compute diagonal from rows and store in ctx->diag */
56 8 : PetscCall(VecDuplicate(d,&ctx->diag));
57 8 : PetscCall(MatGetSize(ctx->A,NULL,&N));
58 8 : PetscCall(MatGetLocalSize(ctx->A,NULL,&n));
59 8 : PetscCall(PetscCalloc2(N,&work1,N,&work2));
60 8 : if (ctx->swapped) {
61 2 : PetscCall(MatGetOwnershipRange(ctx->AT,&start,&end));
62 42 : for (i=start;i<end;i++) {
63 40 : PetscCall(MatGetRow(ctx->AT,i,&ncols,NULL,&vals));
64 120 : for (j=0;j<ncols;j++) work1[i] += vals[j]*vals[j];
65 40 : PetscCall(MatRestoreRow(ctx->AT,i,&ncols,NULL,&vals));
66 : }
67 : } else {
68 6 : PetscCall(MatGetOwnershipRange(ctx->A,&start,&end));
69 137 : for (i=start;i<end;i++) {
70 131 : PetscCall(MatGetRow(ctx->A,i,&ncols,&cols,&vals));
71 894 : for (j=0;j<ncols;j++) work1[cols[j]] += vals[j]*vals[j];
72 131 : PetscCall(MatRestoreRow(ctx->A,i,&ncols,&cols,&vals));
73 : }
74 : }
75 8 : PetscCall(PetscMPIIntCast(N,&len));
76 32 : PetscCall(MPIU_Allreduce(work1,work2,len,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)B)));
77 8 : PetscCall(VecGetOwnershipRange(ctx->diag,&start,&end));
78 8 : PetscCall(VecGetArrayWrite(ctx->diag,&diag));
79 163 : for (i=start;i<end;i++) diag[i-start] = work2[i];
80 8 : PetscCall(VecRestoreArrayWrite(ctx->diag,&diag));
81 8 : PetscCall(PetscFree2(work1,work2));
82 : }
83 8 : PetscCall(VecCopy(ctx->diag,d));
84 8 : PetscFunctionReturn(PETSC_SUCCESS);
85 : }
86 :
87 42 : static PetscErrorCode MatDestroy_Cross(Mat B)
88 : {
89 42 : SVD_CROSS_SHELL *ctx;
90 :
91 42 : PetscFunctionBegin;
92 42 : PetscCall(MatShellGetContext(B,&ctx));
93 42 : PetscCall(VecDestroy(&ctx->w));
94 42 : PetscCall(VecDestroy(&ctx->diag));
95 42 : PetscCall(PetscFree(ctx));
96 42 : PetscFunctionReturn(PETSC_SUCCESS);
97 : }
98 :
99 62 : static PetscErrorCode SVDCrossGetProductMat(SVD svd,Mat A,Mat AT,Mat *C)
100 : {
101 62 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
102 62 : SVD_CROSS_SHELL *ctx;
103 62 : PetscInt n;
104 62 : VecType vtype;
105 62 : Mat B;
106 :
107 62 : PetscFunctionBegin;
108 62 : if (cross->explicitmatrix) {
109 39 : if (!svd->ishyperbolic || svd->swapped) B = (!svd->expltrans && svd->swapped)? AT: A;
110 : else { /* duplicate A and scale by signature */
111 1 : PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B));
112 1 : PetscCall(MatDiagonalScale(B,svd->omega,NULL));
113 : }
114 20 : if (svd->expltrans) { /* explicit transpose */
115 20 : PetscCall(MatProductCreate(AT,B,NULL,C));
116 20 : PetscCall(MatProductSetType(*C,MATPRODUCT_AB));
117 : } else { /* implicit transpose */
118 : #if defined(PETSC_USE_COMPLEX)
119 0 : SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Must use explicit transpose with complex scalars");
120 : #else
121 : if (!svd->swapped) {
122 : PetscCall(MatProductCreate(A,B,NULL,C));
123 : PetscCall(MatProductSetType(*C,MATPRODUCT_AtB));
124 : } else {
125 : PetscCall(MatProductCreate(B,AT,NULL,C));
126 : PetscCall(MatProductSetType(*C,MATPRODUCT_ABt));
127 : }
128 : #endif
129 : }
130 20 : PetscCall(MatProductSetFromOptions(*C));
131 20 : PetscCall(MatProductSymbolic(*C));
132 20 : PetscCall(MatProductNumeric(*C));
133 20 : if (svd->ishyperbolic && !svd->swapped) PetscCall(MatDestroy(&B));
134 : } else {
135 42 : PetscCall(PetscNew(&ctx));
136 42 : ctx->A = A;
137 42 : ctx->AT = AT;
138 42 : ctx->omega = svd->omega;
139 42 : ctx->swapped = svd->swapped;
140 42 : PetscCall(MatCreateVecs(A,NULL,&ctx->w));
141 42 : PetscCall(MatGetLocalSize(A,NULL,&n));
142 42 : PetscCall(MatCreateShell(PetscObjectComm((PetscObject)svd),n,n,PETSC_DETERMINE,PETSC_DETERMINE,(void*)ctx,C));
143 42 : PetscCall(MatShellSetOperation(*C,MATOP_MULT,(void(*)(void))MatMult_Cross));
144 42 : if (!svd->ishyperbolic || svd->swapped) PetscCall(MatShellSetOperation(*C,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cross));
145 42 : PetscCall(MatShellSetOperation(*C,MATOP_DESTROY,(void(*)(void))MatDestroy_Cross));
146 42 : PetscCall(MatGetVecType(A,&vtype));
147 42 : PetscCall(MatSetVecType(*C,vtype));
148 : }
149 62 : PetscFunctionReturn(PETSC_SUCCESS);
150 : }
151 :
152 : /* Convergence test relative to the norm of R (used in GSVD only) */
153 122 : static PetscErrorCode EPSConv_Cross(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
154 : {
155 122 : SVD svd = (SVD)ctx;
156 :
157 122 : PetscFunctionBegin;
158 122 : *errest = res/PetscMax(svd->nrma,svd->nrmb);
159 122 : PetscFunctionReturn(PETSC_SUCCESS);
160 : }
161 :
162 50 : static PetscErrorCode SVDSetUp_Cross(SVD svd)
163 : {
164 50 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
165 50 : ST st;
166 50 : PetscBool trackall,issinv,isks;
167 50 : EPSProblemType ptype;
168 50 : EPSWhich which;
169 50 : Mat Omega;
170 50 : MatType Atype;
171 50 : PetscInt n,N;
172 :
173 50 : PetscFunctionBegin;
174 50 : if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
175 50 : PetscCall(MatDestroy(&cross->C));
176 50 : PetscCall(MatDestroy(&cross->D));
177 50 : PetscCall(SVDCrossGetProductMat(svd,svd->A,svd->AT,&cross->C));
178 50 : if (svd->isgeneralized) {
179 12 : PetscCall(SVDCrossGetProductMat(svd,svd->B,svd->BT,&cross->D));
180 12 : PetscCall(EPSSetOperators(cross->eps,cross->C,cross->D));
181 12 : PetscCall(EPSGetProblemType(cross->eps,&ptype));
182 12 : if (!ptype) PetscCall(EPSSetProblemType(cross->eps,EPS_GHEP));
183 38 : } else if (svd->ishyperbolic && svd->swapped) {
184 2 : PetscCall(MatGetType(svd->OP,&Atype));
185 2 : PetscCall(MatGetSize(svd->A,NULL,&N));
186 2 : PetscCall(MatGetLocalSize(svd->A,NULL,&n));
187 2 : PetscCall(MatCreate(PetscObjectComm((PetscObject)svd),&Omega));
188 2 : PetscCall(MatSetSizes(Omega,n,n,N,N));
189 2 : PetscCall(MatSetType(Omega,Atype));
190 2 : PetscCall(MatDiagonalSet(Omega,svd->omega,INSERT_VALUES));
191 2 : PetscCall(EPSSetOperators(cross->eps,cross->C,Omega));
192 2 : PetscCall(EPSSetProblemType(cross->eps,EPS_GHIEP));
193 2 : PetscCall(MatDestroy(&Omega));
194 : } else {
195 36 : PetscCall(EPSSetOperators(cross->eps,cross->C,NULL));
196 36 : PetscCall(EPSSetProblemType(cross->eps,EPS_HEP));
197 : }
198 50 : if (!cross->usereps) {
199 48 : PetscCall(EPSGetST(cross->eps,&st));
200 48 : PetscCall(PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv));
201 48 : PetscCall(PetscObjectTypeCompare((PetscObject)cross->eps,EPSKRYLOVSCHUR,&isks));
202 48 : if (svd->isgeneralized && svd->which==SVD_SMALLEST) {
203 2 : if (cross->explicitmatrix && isks && !issinv) { /* default to shift-and-invert */
204 1 : PetscCall(STSetType(st,STSINVERT));
205 1 : PetscCall(EPSSetTarget(cross->eps,0.0));
206 : which = EPS_TARGET_REAL;
207 1 : } else which = issinv?EPS_TARGET_REAL:EPS_SMALLEST_REAL;
208 : } else {
209 46 : if (issinv) which = EPS_TARGET_MAGNITUDE;
210 46 : else if (svd->ishyperbolic) which = svd->which==SVD_LARGEST?EPS_LARGEST_MAGNITUDE:EPS_SMALLEST_MAGNITUDE;
211 39 : else which = svd->which==SVD_LARGEST?EPS_LARGEST_REAL:EPS_SMALLEST_REAL;
212 : }
213 48 : PetscCall(EPSSetWhichEigenpairs(cross->eps,which));
214 48 : PetscCall(EPSSetDimensions(cross->eps,svd->nsv,svd->ncv,svd->mpd));
215 65 : PetscCall(EPSSetTolerances(cross->eps,svd->tol==(PetscReal)PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it));
216 48 : switch (svd->conv) {
217 2 : case SVD_CONV_ABS:
218 2 : PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_ABS));break;
219 34 : case SVD_CONV_REL:
220 34 : PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_REL));break;
221 12 : case SVD_CONV_NORM:
222 12 : if (svd->isgeneralized) {
223 12 : if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
224 12 : if (!svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
225 12 : PetscCall(EPSSetConvergenceTestFunction(cross->eps,EPSConv_Cross,svd,NULL));
226 : } else {
227 0 : PetscCall(EPSSetConvergenceTest(cross->eps,EPS_CONV_NORM));break;
228 : }
229 : break;
230 0 : case SVD_CONV_MAXIT:
231 0 : SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Maxit convergence test not supported in this solver");
232 0 : case SVD_CONV_USER:
233 0 : SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
234 : }
235 : }
236 50 : SVDCheckUnsupported(svd,SVD_FEATURE_STOPPING);
237 : /* Transfer the trackall option from svd to eps */
238 50 : PetscCall(SVDGetTrackAll(svd,&trackall));
239 50 : PetscCall(EPSSetTrackAll(cross->eps,trackall));
240 : /* Transfer the initial space from svd to eps */
241 50 : if (svd->nini<0) {
242 5 : PetscCall(EPSSetInitialSpace(cross->eps,-svd->nini,svd->IS));
243 5 : PetscCall(SlepcBasisDestroy_Private(&svd->nini,&svd->IS));
244 : }
245 50 : PetscCall(EPSSetUp(cross->eps));
246 50 : PetscCall(EPSGetDimensions(cross->eps,NULL,&svd->ncv,&svd->mpd));
247 50 : PetscCall(EPSGetTolerances(cross->eps,NULL,&svd->max_it));
248 50 : if (svd->tol==(PetscReal)PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
249 :
250 50 : svd->leftbasis = PETSC_FALSE;
251 50 : PetscCall(SVDAllocateSolution(svd,0));
252 50 : PetscFunctionReturn(PETSC_SUCCESS);
253 : }
254 :
255 50 : static PetscErrorCode SVDSolve_Cross(SVD svd)
256 : {
257 50 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
258 50 : PetscInt i;
259 50 : PetscScalar lambda;
260 50 : PetscReal sigma;
261 :
262 50 : PetscFunctionBegin;
263 50 : PetscCall(EPSSolve(cross->eps));
264 50 : PetscCall(EPSGetConverged(cross->eps,&svd->nconv));
265 50 : PetscCall(EPSGetIterationNumber(cross->eps,&svd->its));
266 50 : PetscCall(EPSGetConvergedReason(cross->eps,(EPSConvergedReason*)&svd->reason));
267 620 : for (i=0;i<svd->nconv;i++) {
268 570 : PetscCall(EPSGetEigenvalue(cross->eps,i,&lambda,NULL));
269 570 : sigma = PetscRealPart(lambda);
270 570 : if (svd->ishyperbolic) svd->sigma[i] = PetscSqrtReal(PetscAbsReal(sigma));
271 : else {
272 332 : PetscCheck(sigma>-10*PETSC_MACHINE_EPSILON,PetscObjectComm((PetscObject)svd),PETSC_ERR_FP,"Negative eigenvalue computed by EPS: %g",(double)sigma);
273 332 : if (sigma<0.0) {
274 2 : PetscCall(PetscInfo(svd,"Negative eigenvalue computed by EPS: %g, resetting to 0\n",(double)sigma));
275 : sigma = 0.0;
276 : }
277 332 : svd->sigma[i] = PetscSqrtReal(sigma);
278 : }
279 : }
280 50 : PetscFunctionReturn(PETSC_SUCCESS);
281 : }
282 :
283 44 : static PetscErrorCode SVDComputeVectors_Cross(SVD svd)
284 : {
285 44 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
286 44 : PetscInt i,mloc,ploc;
287 44 : Vec u,v,x,uv,w,omega2=NULL;
288 44 : Mat Omega;
289 44 : PetscScalar *dst,alpha,lambda,*varray;
290 44 : const PetscScalar *src;
291 44 : PetscReal nrm;
292 :
293 44 : PetscFunctionBegin;
294 44 : if (svd->isgeneralized) {
295 12 : PetscCall(MatCreateVecs(svd->A,NULL,&u));
296 12 : PetscCall(VecGetLocalSize(u,&mloc));
297 12 : PetscCall(MatCreateVecs(svd->B,NULL,&v));
298 12 : PetscCall(VecGetLocalSize(v,&ploc));
299 114 : for (i=0;i<svd->nconv;i++) {
300 102 : PetscCall(BVGetColumn(svd->V,i,&x));
301 102 : PetscCall(EPSGetEigenpair(cross->eps,i,&lambda,NULL,x,NULL));
302 102 : PetscCall(MatMult(svd->A,x,u)); /* u_i*c_i/alpha = A*x_i */
303 102 : PetscCall(VecNormalize(u,NULL));
304 102 : PetscCall(MatMult(svd->B,x,v)); /* v_i*s_i/alpha = B*x_i */
305 102 : PetscCall(VecNormalize(v,&nrm)); /* ||v||_2 = s_i/alpha */
306 102 : alpha = 1.0/(PetscSqrtReal(1.0+PetscRealPart(lambda))*nrm); /* alpha=s_i/||v||_2 */
307 102 : PetscCall(VecScale(x,alpha));
308 102 : PetscCall(BVRestoreColumn(svd->V,i,&x));
309 : /* copy [u;v] to U[i] */
310 102 : PetscCall(BVGetColumn(svd->U,i,&uv));
311 102 : PetscCall(VecGetArrayWrite(uv,&dst));
312 102 : PetscCall(VecGetArrayRead(u,&src));
313 102 : PetscCall(PetscArraycpy(dst,src,mloc));
314 102 : PetscCall(VecRestoreArrayRead(u,&src));
315 102 : PetscCall(VecGetArrayRead(v,&src));
316 102 : PetscCall(PetscArraycpy(dst+mloc,src,ploc));
317 102 : PetscCall(VecRestoreArrayRead(v,&src));
318 102 : PetscCall(VecRestoreArrayWrite(uv,&dst));
319 102 : PetscCall(BVRestoreColumn(svd->U,i,&uv));
320 : }
321 12 : PetscCall(VecDestroy(&v));
322 12 : PetscCall(VecDestroy(&u));
323 32 : } else if (svd->ishyperbolic && svd->swapped) { /* was solved as GHIEP, set u=Omega*u and normalize */
324 2 : PetscCall(EPSGetOperators(cross->eps,NULL,&Omega));
325 2 : PetscCall(MatCreateVecs(Omega,&w,NULL));
326 2 : PetscCall(VecCreateSeq(PETSC_COMM_SELF,svd->ncv,&omega2));
327 2 : PetscCall(VecGetArrayWrite(omega2,&varray));
328 58 : for (i=0;i<svd->nconv;i++) {
329 56 : PetscCall(BVGetColumn(svd->V,i,&v));
330 56 : PetscCall(EPSGetEigenvector(cross->eps,i,v,NULL));
331 56 : PetscCall(MatMult(Omega,v,w));
332 56 : PetscCall(VecDot(v,w,&alpha));
333 56 : svd->sign[i] = PetscSign(PetscRealPart(alpha));
334 56 : varray[i] = svd->sign[i];
335 56 : alpha = 1.0/PetscSqrtScalar(PetscAbsScalar(alpha));
336 56 : PetscCall(VecScale(w,alpha));
337 56 : PetscCall(VecCopy(w,v));
338 56 : PetscCall(BVRestoreColumn(svd->V,i,&v));
339 : }
340 2 : PetscCall(BVSetSignature(svd->V,omega2));
341 2 : PetscCall(VecRestoreArrayWrite(omega2,&varray));
342 2 : PetscCall(VecDestroy(&omega2));
343 2 : PetscCall(VecDestroy(&w));
344 2 : PetscCall(SVDComputeVectors_Left(svd));
345 : } else {
346 426 : for (i=0;i<svd->nconv;i++) {
347 396 : PetscCall(BVGetColumn(svd->V,i,&v));
348 396 : PetscCall(EPSGetEigenvector(cross->eps,i,v,NULL));
349 396 : PetscCall(BVRestoreColumn(svd->V,i,&v));
350 : }
351 30 : PetscCall(SVDComputeVectors_Left(svd));
352 : }
353 44 : PetscFunctionReturn(PETSC_SUCCESS);
354 : }
355 :
356 198 : static PetscErrorCode EPSMonitor_Cross(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
357 : {
358 198 : PetscInt i;
359 198 : SVD svd = (SVD)ctx;
360 198 : PetscScalar er,ei;
361 198 : ST st;
362 :
363 198 : PetscFunctionBegin;
364 198 : PetscCall(EPSGetST(eps,&st));
365 3052 : for (i=0;i<PetscMin(nest,svd->ncv);i++) {
366 2854 : er = eigr[i]; ei = eigi[i];
367 2854 : PetscCall(STBackTransform(st,1,&er,&ei));
368 2854 : svd->sigma[i] = PetscSqrtReal(PetscAbsReal(PetscRealPart(er)));
369 2854 : svd->errest[i] = errest[i];
370 : }
371 198 : PetscCall(SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest));
372 198 : PetscFunctionReturn(PETSC_SUCCESS);
373 : }
374 :
375 37 : static PetscErrorCode SVDSetFromOptions_Cross(SVD svd,PetscOptionItems *PetscOptionsObject)
376 : {
377 37 : PetscBool set,val;
378 37 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
379 37 : ST st;
380 :
381 37 : PetscFunctionBegin;
382 37 : PetscOptionsHeadBegin(PetscOptionsObject,"SVD Cross Options");
383 :
384 37 : PetscCall(PetscOptionsBool("-svd_cross_explicitmatrix","Use cross explicit matrix","SVDCrossSetExplicitMatrix",cross->explicitmatrix,&val,&set));
385 37 : if (set) PetscCall(SVDCrossSetExplicitMatrix(svd,val));
386 :
387 37 : PetscOptionsHeadEnd();
388 :
389 37 : if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
390 37 : if (!cross->explicitmatrix && !cross->usereps) {
391 : /* use as default an ST with shell matrix and Jacobi */
392 24 : PetscCall(EPSGetST(cross->eps,&st));
393 24 : PetscCall(STSetMatMode(st,ST_MATMODE_SHELL));
394 : }
395 37 : PetscCall(EPSSetFromOptions(cross->eps));
396 37 : PetscFunctionReturn(PETSC_SUCCESS);
397 : }
398 :
399 18 : static PetscErrorCode SVDCrossSetExplicitMatrix_Cross(SVD svd,PetscBool explicitmatrix)
400 : {
401 18 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
402 :
403 18 : PetscFunctionBegin;
404 18 : if (cross->explicitmatrix != explicitmatrix) {
405 11 : cross->explicitmatrix = explicitmatrix;
406 11 : svd->state = SVD_STATE_INITIAL;
407 : }
408 18 : PetscFunctionReturn(PETSC_SUCCESS);
409 : }
410 :
411 : /*@
412 : SVDCrossSetExplicitMatrix - Indicate if the eigensolver operator A^T*A must
413 : be computed explicitly.
414 :
415 : Logically Collective
416 :
417 : Input Parameters:
418 : + svd - singular value solver
419 : - explicitmat - boolean flag indicating if A^T*A is built explicitly
420 :
421 : Options Database Key:
422 : . -svd_cross_explicitmatrix <boolean> - Indicates the boolean flag
423 :
424 : Level: advanced
425 :
426 : .seealso: SVDCrossGetExplicitMatrix()
427 : @*/
428 18 : PetscErrorCode SVDCrossSetExplicitMatrix(SVD svd,PetscBool explicitmat)
429 : {
430 18 : PetscFunctionBegin;
431 18 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
432 72 : PetscValidLogicalCollectiveBool(svd,explicitmat,2);
433 18 : PetscTryMethod(svd,"SVDCrossSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmat));
434 18 : PetscFunctionReturn(PETSC_SUCCESS);
435 : }
436 :
437 2 : static PetscErrorCode SVDCrossGetExplicitMatrix_Cross(SVD svd,PetscBool *explicitmat)
438 : {
439 2 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
440 :
441 2 : PetscFunctionBegin;
442 2 : *explicitmat = cross->explicitmatrix;
443 2 : PetscFunctionReturn(PETSC_SUCCESS);
444 : }
445 :
446 : /*@
447 : SVDCrossGetExplicitMatrix - Returns the flag indicating if A^T*A is built explicitly.
448 :
449 : Not Collective
450 :
451 : Input Parameter:
452 : . svd - singular value solver
453 :
454 : Output Parameter:
455 : . explicitmat - the mode flag
456 :
457 : Level: advanced
458 :
459 : .seealso: SVDCrossSetExplicitMatrix()
460 : @*/
461 2 : PetscErrorCode SVDCrossGetExplicitMatrix(SVD svd,PetscBool *explicitmat)
462 : {
463 2 : PetscFunctionBegin;
464 2 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
465 2 : PetscAssertPointer(explicitmat,2);
466 2 : PetscUseMethod(svd,"SVDCrossGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmat));
467 2 : PetscFunctionReturn(PETSC_SUCCESS);
468 : }
469 :
470 2 : static PetscErrorCode SVDCrossSetEPS_Cross(SVD svd,EPS eps)
471 : {
472 2 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
473 :
474 2 : PetscFunctionBegin;
475 2 : PetscCall(PetscObjectReference((PetscObject)eps));
476 2 : PetscCall(EPSDestroy(&cross->eps));
477 2 : cross->eps = eps;
478 2 : cross->usereps = PETSC_TRUE;
479 2 : svd->state = SVD_STATE_INITIAL;
480 2 : PetscFunctionReturn(PETSC_SUCCESS);
481 : }
482 :
483 : /*@
484 : SVDCrossSetEPS - Associate an eigensolver object (EPS) to the
485 : singular value solver.
486 :
487 : Collective
488 :
489 : Input Parameters:
490 : + svd - singular value solver
491 : - eps - the eigensolver object
492 :
493 : Level: advanced
494 :
495 : .seealso: SVDCrossGetEPS()
496 : @*/
497 2 : PetscErrorCode SVDCrossSetEPS(SVD svd,EPS eps)
498 : {
499 2 : PetscFunctionBegin;
500 2 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
501 2 : PetscValidHeaderSpecific(eps,EPS_CLASSID,2);
502 2 : PetscCheckSameComm(svd,1,eps,2);
503 2 : PetscTryMethod(svd,"SVDCrossSetEPS_C",(SVD,EPS),(svd,eps));
504 2 : PetscFunctionReturn(PETSC_SUCCESS);
505 : }
506 :
507 37 : static PetscErrorCode SVDCrossGetEPS_Cross(SVD svd,EPS *eps)
508 : {
509 37 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
510 :
511 37 : PetscFunctionBegin;
512 37 : if (!cross->eps) {
513 37 : PetscCall(EPSCreate(PetscObjectComm((PetscObject)svd),&cross->eps));
514 37 : PetscCall(PetscObjectIncrementTabLevel((PetscObject)cross->eps,(PetscObject)svd,1));
515 37 : PetscCall(EPSSetOptionsPrefix(cross->eps,((PetscObject)svd)->prefix));
516 37 : PetscCall(EPSAppendOptionsPrefix(cross->eps,"svd_cross_"));
517 37 : PetscCall(PetscObjectSetOptions((PetscObject)cross->eps,((PetscObject)svd)->options));
518 37 : PetscCall(EPSSetWhichEigenpairs(cross->eps,EPS_LARGEST_REAL));
519 37 : PetscCall(EPSMonitorSet(cross->eps,EPSMonitor_Cross,svd,NULL));
520 : }
521 37 : *eps = cross->eps;
522 37 : PetscFunctionReturn(PETSC_SUCCESS);
523 : }
524 :
525 : /*@
526 : SVDCrossGetEPS - Retrieve the eigensolver object (EPS) associated
527 : to the singular value solver.
528 :
529 : Collective
530 :
531 : Input Parameter:
532 : . svd - singular value solver
533 :
534 : Output Parameter:
535 : . eps - the eigensolver object
536 :
537 : Level: advanced
538 :
539 : .seealso: SVDCrossSetEPS()
540 : @*/
541 37 : PetscErrorCode SVDCrossGetEPS(SVD svd,EPS *eps)
542 : {
543 37 : PetscFunctionBegin;
544 37 : PetscValidHeaderSpecific(svd,SVD_CLASSID,1);
545 37 : PetscAssertPointer(eps,2);
546 37 : PetscUseMethod(svd,"SVDCrossGetEPS_C",(SVD,EPS*),(svd,eps));
547 37 : PetscFunctionReturn(PETSC_SUCCESS);
548 : }
549 :
550 1 : static PetscErrorCode SVDView_Cross(SVD svd,PetscViewer viewer)
551 : {
552 1 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
553 1 : PetscBool isascii;
554 :
555 1 : PetscFunctionBegin;
556 1 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
557 1 : if (isascii) {
558 1 : if (!cross->eps) PetscCall(SVDCrossGetEPS(svd,&cross->eps));
559 2 : PetscCall(PetscViewerASCIIPrintf(viewer," %s matrix\n",cross->explicitmatrix?"explicit":"implicit"));
560 1 : PetscCall(PetscViewerASCIIPushTab(viewer));
561 1 : PetscCall(EPSView(cross->eps,viewer));
562 1 : PetscCall(PetscViewerASCIIPopTab(viewer));
563 : }
564 1 : PetscFunctionReturn(PETSC_SUCCESS);
565 : }
566 :
567 40 : static PetscErrorCode SVDReset_Cross(SVD svd)
568 : {
569 40 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
570 :
571 40 : PetscFunctionBegin;
572 40 : PetscCall(EPSReset(cross->eps));
573 40 : PetscCall(MatDestroy(&cross->C));
574 40 : PetscCall(MatDestroy(&cross->D));
575 40 : PetscFunctionReturn(PETSC_SUCCESS);
576 : }
577 :
578 39 : static PetscErrorCode SVDDestroy_Cross(SVD svd)
579 : {
580 39 : SVD_CROSS *cross = (SVD_CROSS*)svd->data;
581 :
582 39 : PetscFunctionBegin;
583 39 : PetscCall(EPSDestroy(&cross->eps));
584 39 : PetscCall(PetscFree(svd->data));
585 39 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",NULL));
586 39 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",NULL));
587 39 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",NULL));
588 39 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",NULL));
589 39 : PetscFunctionReturn(PETSC_SUCCESS);
590 : }
591 :
592 39 : SLEPC_EXTERN PetscErrorCode SVDCreate_Cross(SVD svd)
593 : {
594 39 : SVD_CROSS *cross;
595 :
596 39 : PetscFunctionBegin;
597 39 : PetscCall(PetscNew(&cross));
598 39 : svd->data = (void*)cross;
599 :
600 39 : svd->ops->solve = SVDSolve_Cross;
601 39 : svd->ops->solveg = SVDSolve_Cross;
602 39 : svd->ops->solveh = SVDSolve_Cross;
603 39 : svd->ops->setup = SVDSetUp_Cross;
604 39 : svd->ops->setfromoptions = SVDSetFromOptions_Cross;
605 39 : svd->ops->destroy = SVDDestroy_Cross;
606 39 : svd->ops->reset = SVDReset_Cross;
607 39 : svd->ops->view = SVDView_Cross;
608 39 : svd->ops->computevectors = SVDComputeVectors_Cross;
609 39 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetEPS_C",SVDCrossSetEPS_Cross));
610 39 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetEPS_C",SVDCrossGetEPS_Cross));
611 39 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossSetExplicitMatrix_C",SVDCrossSetExplicitMatrix_Cross));
612 39 : PetscCall(PetscObjectComposeFunction((PetscObject)svd,"SVDCrossGetExplicitMatrix_C",SVDCrossGetExplicitMatrix_Cross));
613 39 : PetscFunctionReturn(PETSC_SUCCESS);
614 : }
|