Actual source code: cyclic.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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: "cyclic"
13: Method: Uses a Hermitian eigensolver for H(A) = [ 0 A ; A^T 0 ]
14: */
16: #include <slepc/private/svdimpl.h> /*I "slepcsvd.h" I*/
17: #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
18: #include "./cyclicimpl.h"
20: static PetscErrorCode MatMult_Cyclic(Mat B,Vec x,Vec y)
21: {
22: PetscErrorCode ierr;
23: SVD svd;
24: SVD_CYCLIC *cyclic;
25: const PetscScalar *px;
26: PetscScalar *py;
27: PetscInt m;
30: MatShellGetContext(B,(void**)&svd);
31: cyclic = (SVD_CYCLIC*)svd->data;
32: SVDMatGetLocalSize(svd,&m,NULL);
33: VecGetArrayRead(x,&px);
34: VecGetArray(y,&py);
35: VecPlaceArray(cyclic->x1,px);
36: VecPlaceArray(cyclic->x2,px+m);
37: VecPlaceArray(cyclic->y1,py);
38: VecPlaceArray(cyclic->y2,py+m);
39: SVDMatMult(svd,PETSC_FALSE,cyclic->x2,cyclic->y1);
40: SVDMatMult(svd,PETSC_TRUE,cyclic->x1,cyclic->y2);
41: VecResetArray(cyclic->x1);
42: VecResetArray(cyclic->x2);
43: VecResetArray(cyclic->y1);
44: VecResetArray(cyclic->y2);
45: VecRestoreArrayRead(x,&px);
46: VecRestoreArray(y,&py);
47: return(0);
48: }
50: static PetscErrorCode MatGetDiagonal_Cyclic(Mat B,Vec diag)
51: {
55: VecSet(diag,0.0);
56: return(0);
57: }
59: PetscErrorCode SVDSetUp_Cyclic(SVD svd)
60: {
61: PetscErrorCode ierr;
62: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
63: PetscInt M,N,m,n,i,isl,Istart,Iend;
64: const PetscScalar *isa;
65: PetscScalar *va;
66: PetscBool trackall,issinv;
67: Vec v;
68: Mat Zm,Zn;
69: ST st;
70: #if defined(PETSC_HAVE_CUDA)
71: PetscBool cuda;
72: #endif
75: SVDMatGetSize(svd,&M,&N);
76: SVDMatGetLocalSize(svd,&m,&n);
77: if (!cyclic->mat) {
78: if (cyclic->explicitmatrix) {
79: if (!svd->A || !svd->AT) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Cannot use explicit cyclic matrix with implicit transpose");
80: MatCreate(PetscObjectComm((PetscObject)svd),&Zm);
81: MatSetSizes(Zm,m,m,M,M);
82: MatSetFromOptions(Zm);
83: MatSetUp(Zm);
84: MatGetOwnershipRange(Zm,&Istart,&Iend);
85: for (i=Istart;i<Iend;i++) {
86: MatSetValue(Zm,i,i,0.0,INSERT_VALUES);
87: }
88: MatAssemblyBegin(Zm,MAT_FINAL_ASSEMBLY);
89: MatAssemblyEnd(Zm,MAT_FINAL_ASSEMBLY);
90: MatCreate(PetscObjectComm((PetscObject)svd),&Zn);
91: MatSetSizes(Zn,n,n,N,N);
92: MatSetFromOptions(Zn);
93: MatSetUp(Zn);
94: MatGetOwnershipRange(Zn,&Istart,&Iend);
95: for (i=Istart;i<Iend;i++) {
96: MatSetValue(Zn,i,i,0.0,INSERT_VALUES);
97: }
98: MatAssemblyBegin(Zn,MAT_FINAL_ASSEMBLY);
99: MatAssemblyEnd(Zn,MAT_FINAL_ASSEMBLY);
100: MatCreateTile(1.0,Zm,1.0,svd->A,1.0,svd->AT,1.0,Zn,&cyclic->mat);
101: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->mat);
102: MatDestroy(&Zm);
103: MatDestroy(&Zn);
104: } else {
105: SVDMatCreateVecsEmpty(svd,&cyclic->x2,&cyclic->x1);
106: SVDMatCreateVecsEmpty(svd,&cyclic->y2,&cyclic->y1);
107: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->x1);
108: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->x2);
109: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->y1);
110: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->y2);
111: MatCreateShell(PetscObjectComm((PetscObject)svd),m+n,m+n,M+N,M+N,svd,&cyclic->mat);
112: MatShellSetOperation(cyclic->mat,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Cyclic);
113: #if defined(PETSC_HAVE_CUDA)
114: PetscObjectTypeCompareAny((PetscObject)(svd->A?svd->A:svd->AT),&cuda,MATSEQAIJCUSPARSE,MATMPIAIJCUSPARSE,"");
115: if (cuda) {
116: MatShellSetOperation(cyclic->mat,MATOP_CREATE_VECS,(void(*)(void))MatCreateVecs_Cyclic_CUDA);
117: MatShellSetOperation(cyclic->mat,MATOP_MULT,(void(*)(void))MatMult_Cyclic_CUDA);
118: } else
119: #endif
120: {
121: MatShellSetOperation(cyclic->mat,MATOP_MULT,(void(*)(void))MatMult_Cyclic);
122: }
123: }
124: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->mat);
125: }
127: if (!cyclic->eps) { SVDCyclicGetEPS(svd,&cyclic->eps); }
128: EPSSetOperators(cyclic->eps,cyclic->mat,NULL);
129: EPSSetProblemType(cyclic->eps,EPS_HEP);
130: if (!cyclic->usereps) {
131: if (svd->which == SVD_LARGEST) {
132: EPSGetST(cyclic->eps,&st);
133: PetscObjectTypeCompare((PetscObject)st,STSINVERT,&issinv);
134: if (issinv) {
135: EPSSetWhichEigenpairs(cyclic->eps,EPS_TARGET_MAGNITUDE);
136: } else {
137: EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
138: }
139: } else {
140: EPSSetEigenvalueComparison(cyclic->eps,SlepcCompareSmallestPosReal,NULL);
141: EPSSetTarget(cyclic->eps,0.0);
142: }
143: EPSSetDimensions(cyclic->eps,svd->nsv,svd->ncv?svd->ncv:PETSC_DEFAULT,svd->mpd?svd->mpd:PETSC_DEFAULT);
144: EPSSetTolerances(cyclic->eps,svd->tol==PETSC_DEFAULT?SLEPC_DEFAULT_TOL/10.0:svd->tol,svd->max_it?svd->max_it:PETSC_DEFAULT);
145: switch (svd->conv) {
146: case SVD_CONV_ABS:
147: EPSSetConvergenceTest(cyclic->eps,EPS_CONV_ABS);break;
148: case SVD_CONV_REL:
149: EPSSetConvergenceTest(cyclic->eps,EPS_CONV_REL);break;
150: case SVD_CONV_USER:
151: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined convergence test not supported in this solver");
152: }
153: }
154: if (svd->stop!=SVD_STOP_BASIC) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"User-defined stopping test not supported in this solver");
155: /* Transfer the trackall option from svd to eps */
156: SVDGetTrackAll(svd,&trackall);
157: EPSSetTrackAll(cyclic->eps,trackall);
158: /* Transfer the initial subspace from svd to eps */
159: if (svd->nini<0 || svd->ninil<0) {
160: for (i=0;i<-PetscMin(svd->nini,svd->ninil);i++) {
161: MatCreateVecs(cyclic->mat,&v,NULL);
162: VecGetArray(v,&va);
163: if (i<-svd->ninil) {
164: VecGetSize(svd->ISL[i],&isl);
165: if (isl!=m) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for left initial vector");
166: VecGetArrayRead(svd->ISL[i],&isa);
167: PetscMemcpy(va,isa,sizeof(PetscScalar)*m);
168: VecRestoreArrayRead(svd->IS[i],&isa);
169: } else {
170: PetscMemzero(&va,sizeof(PetscScalar)*m);
171: }
172: if (i<-svd->nini) {
173: VecGetSize(svd->IS[i],&isl);
174: if (isl!=n) SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_SUP,"Size mismatch for right initial vector");
175: VecGetArrayRead(svd->IS[i],&isa);
176: PetscMemcpy(va+m,isa,sizeof(PetscScalar)*n);
177: VecRestoreArrayRead(svd->IS[i],&isa);
178: } else {
179: PetscMemzero(va+m,sizeof(PetscScalar)*n);
180: }
181: VecRestoreArray(v,&va);
182: VecDestroy(&svd->IS[i]);
183: svd->IS[i] = v;
184: }
185: svd->nini = PetscMin(svd->nini,svd->ninil);
186: EPSSetInitialSpace(cyclic->eps,-svd->nini,svd->IS);
187: SlepcBasisDestroy_Private(&svd->nini,&svd->IS);
188: SlepcBasisDestroy_Private(&svd->ninil,&svd->ISL);
189: }
190: EPSSetUp(cyclic->eps);
191: EPSGetDimensions(cyclic->eps,NULL,&svd->ncv,&svd->mpd);
192: svd->ncv = PetscMin(svd->ncv,PetscMin(M,N));
193: EPSGetTolerances(cyclic->eps,NULL,&svd->max_it);
194: if (svd->tol==PETSC_DEFAULT) svd->tol = SLEPC_DEFAULT_TOL;
196: svd->leftbasis = PETSC_TRUE;
197: SVDAllocateSolution(svd,0);
198: return(0);
199: }
201: PetscErrorCode SVDSolve_Cyclic(SVD svd)
202: {
203: PetscErrorCode ierr;
204: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
205: PetscInt i,j,M,N,m,n;
206: PetscScalar sigma;
207: const PetscScalar *px;
208: Vec x,x1,x2;
211: EPSSolve(cyclic->eps);
212: EPSGetConverged(cyclic->eps,&svd->nconv);
213: EPSGetIterationNumber(cyclic->eps,&svd->its);
214: EPSGetConvergedReason(cyclic->eps,(EPSConvergedReason*)&svd->reason);
216: MatCreateVecs(cyclic->mat,&x,NULL);
217: SVDMatGetSize(svd,&M,&N);
218: SVDMatGetLocalSize(svd,&m,&n);
219: SVDMatCreateVecsEmpty(svd,&x2,&x1);
220: for (i=0,j=0;i<svd->nconv;i++) {
221: EPSGetEigenpair(cyclic->eps,i,&sigma,NULL,x,NULL);
222: if (PetscRealPart(sigma) > 0.0) {
223: svd->sigma[j] = PetscRealPart(sigma);
224: VecGetArrayRead(x,&px);
225: VecPlaceArray(x1,px);
226: VecPlaceArray(x2,px+m);
227: BVInsertVec(svd->U,j,x1);
228: BVScaleColumn(svd->U,j,1.0/PetscSqrtReal(2.0));
229: BVInsertVec(svd->V,j,x2);
230: BVScaleColumn(svd->V,j,1.0/PetscSqrtReal(2.0));
231: VecResetArray(x1);
232: VecResetArray(x2);
233: VecRestoreArrayRead(x,&px);
234: j++;
235: }
236: }
237: svd->nconv = j;
239: VecDestroy(&x);
240: VecDestroy(&x1);
241: VecDestroy(&x2);
242: return(0);
243: }
245: static PetscErrorCode EPSMonitor_Cyclic(EPS eps,PetscInt its,PetscInt nconv,PetscScalar *eigr,PetscScalar *eigi,PetscReal *errest,PetscInt nest,void *ctx)
246: {
247: PetscInt i,j;
248: SVD svd = (SVD)ctx;
249: PetscScalar er,ei;
253: nconv = 0;
254: for (i=0,j=0;i<PetscMin(nest,svd->ncv);i++) {
255: er = eigr[i]; ei = eigi[i];
256: STBackTransform(eps->st,1,&er,&ei);
257: if (PetscRealPart(er) > 0.0) {
258: svd->sigma[j] = PetscRealPart(er);
259: svd->errest[j] = errest[i];
260: if (errest[i] && errest[i] < svd->tol) nconv++;
261: j++;
262: }
263: }
264: nest = j;
265: SVDMonitor(svd,its,nconv,svd->sigma,svd->errest,nest);
266: return(0);
267: }
269: PetscErrorCode SVDSetFromOptions_Cyclic(PetscOptionItems *PetscOptionsObject,SVD svd)
270: {
272: PetscBool set,val;
273: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
274: ST st;
277: PetscOptionsHead(PetscOptionsObject,"SVD Cyclic Options");
279: PetscOptionsBool("-svd_cyclic_explicitmatrix","Use cyclic explicit matrix","SVDCyclicSetExplicitMatrix",cyclic->explicitmatrix,&val,&set);
280: if (set) { SVDCyclicSetExplicitMatrix(svd,val); }
282: PetscOptionsTail();
284: if (!cyclic->eps) { SVDCyclicGetEPS(svd,&cyclic->eps); }
285: if (!cyclic->explicitmatrix && !cyclic->usereps) {
286: /* use as default an ST with shell matrix and Jacobi */
287: EPSGetST(cyclic->eps,&st);
288: STSetMatMode(st,ST_MATMODE_SHELL);
289: }
290: EPSSetFromOptions(cyclic->eps);
291: return(0);
292: }
294: static PetscErrorCode SVDCyclicSetExplicitMatrix_Cyclic(SVD svd,PetscBool explicitmatrix)
295: {
296: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
299: if (cyclic->explicitmatrix != explicitmatrix) {
300: cyclic->explicitmatrix = explicitmatrix;
301: svd->state = SVD_STATE_INITIAL;
302: }
303: return(0);
304: }
306: /*@
307: SVDCyclicSetExplicitMatrix - Indicate if the eigensolver operator
308: H(A) = [ 0 A ; A^T 0 ] must be computed explicitly.
310: Logically Collective on SVD
312: Input Parameters:
313: + svd - singular value solver
314: - explicit - boolean flag indicating if H(A) is built explicitly
316: Options Database Key:
317: . -svd_cyclic_explicitmatrix <boolean> - Indicates the boolean flag
319: Level: advanced
321: .seealso: SVDCyclicGetExplicitMatrix()
322: @*/
323: PetscErrorCode SVDCyclicSetExplicitMatrix(SVD svd,PetscBool explicitmatrix)
324: {
330: PetscTryMethod(svd,"SVDCyclicSetExplicitMatrix_C",(SVD,PetscBool),(svd,explicitmatrix));
331: return(0);
332: }
334: static PetscErrorCode SVDCyclicGetExplicitMatrix_Cyclic(SVD svd,PetscBool *explicitmatrix)
335: {
336: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
339: *explicitmatrix = cyclic->explicitmatrix;
340: return(0);
341: }
343: /*@
344: SVDCyclicGetExplicitMatrix - Returns the flag indicating if H(A) is built explicitly.
346: Not Collective
348: Input Parameter:
349: . svd - singular value solver
351: Output Parameter:
352: . explicit - the mode flag
354: Level: advanced
356: .seealso: SVDCyclicSetExplicitMatrix()
357: @*/
358: PetscErrorCode SVDCyclicGetExplicitMatrix(SVD svd,PetscBool *explicitmatrix)
359: {
365: PetscUseMethod(svd,"SVDCyclicGetExplicitMatrix_C",(SVD,PetscBool*),(svd,explicitmatrix));
366: return(0);
367: }
369: static PetscErrorCode SVDCyclicSetEPS_Cyclic(SVD svd,EPS eps)
370: {
371: PetscErrorCode ierr;
372: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
375: PetscObjectReference((PetscObject)eps);
376: EPSDestroy(&cyclic->eps);
377: cyclic->eps = eps;
378: cyclic->usereps = PETSC_TRUE;
379: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->eps);
380: svd->state = SVD_STATE_INITIAL;
381: return(0);
382: }
384: /*@
385: SVDCyclicSetEPS - Associate an eigensolver object (EPS) to the
386: singular value solver.
388: Collective on SVD
390: Input Parameters:
391: + svd - singular value solver
392: - eps - the eigensolver object
394: Level: advanced
396: .seealso: SVDCyclicGetEPS()
397: @*/
398: PetscErrorCode SVDCyclicSetEPS(SVD svd,EPS eps)
399: {
406: PetscTryMethod(svd,"SVDCyclicSetEPS_C",(SVD,EPS),(svd,eps));
407: return(0);
408: }
410: static PetscErrorCode SVDCyclicGetEPS_Cyclic(SVD svd,EPS *eps)
411: {
413: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
416: if (!cyclic->eps) {
417: EPSCreate(PetscObjectComm((PetscObject)svd),&cyclic->eps);
418: PetscObjectIncrementTabLevel((PetscObject)cyclic->eps,(PetscObject)svd,1);
419: EPSSetOptionsPrefix(cyclic->eps,((PetscObject)svd)->prefix);
420: EPSAppendOptionsPrefix(cyclic->eps,"svd_cyclic_");
421: PetscLogObjectParent((PetscObject)svd,(PetscObject)cyclic->eps);
422: PetscObjectSetOptions((PetscObject)cyclic->eps,((PetscObject)svd)->options);
423: EPSSetWhichEigenpairs(cyclic->eps,EPS_LARGEST_REAL);
424: EPSMonitorSet(cyclic->eps,EPSMonitor_Cyclic,svd,NULL);
425: }
426: *eps = cyclic->eps;
427: return(0);
428: }
430: /*@
431: SVDCyclicGetEPS - Retrieve the eigensolver object (EPS) associated
432: to the singular value solver.
434: Not Collective
436: Input Parameter:
437: . svd - singular value solver
439: Output Parameter:
440: . eps - the eigensolver object
442: Level: advanced
444: .seealso: SVDCyclicSetEPS()
445: @*/
446: PetscErrorCode SVDCyclicGetEPS(SVD svd,EPS *eps)
447: {
453: PetscUseMethod(svd,"SVDCyclicGetEPS_C",(SVD,EPS*),(svd,eps));
454: return(0);
455: }
457: PetscErrorCode SVDView_Cyclic(SVD svd,PetscViewer viewer)
458: {
460: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
461: PetscBool isascii;
464: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
465: if (isascii) {
466: if (!cyclic->eps) { SVDCyclicGetEPS(svd,&cyclic->eps); }
467: PetscViewerASCIIPrintf(viewer," %s matrix\n",cyclic->explicitmatrix?"explicit":"implicit");
468: PetscViewerASCIIPushTab(viewer);
469: EPSView(cyclic->eps,viewer);
470: PetscViewerASCIIPopTab(viewer);
471: }
472: return(0);
473: }
475: PetscErrorCode SVDReset_Cyclic(SVD svd)
476: {
478: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
481: EPSReset(cyclic->eps);
482: MatDestroy(&cyclic->mat);
483: VecDestroy(&cyclic->x1);
484: VecDestroy(&cyclic->x2);
485: VecDestroy(&cyclic->y1);
486: VecDestroy(&cyclic->y2);
487: return(0);
488: }
490: PetscErrorCode SVDDestroy_Cyclic(SVD svd)
491: {
493: SVD_CYCLIC *cyclic = (SVD_CYCLIC*)svd->data;
496: EPSDestroy(&cyclic->eps);
497: PetscFree(svd->data);
498: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",NULL);
499: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",NULL);
500: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",NULL);
501: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",NULL);
502: return(0);
503: }
505: SLEPC_EXTERN PetscErrorCode SVDCreate_Cyclic(SVD svd)
506: {
508: SVD_CYCLIC *cyclic;
511: PetscNewLog(svd,&cyclic);
512: svd->data = (void*)cyclic;
513: svd->ops->solve = SVDSolve_Cyclic;
514: svd->ops->setup = SVDSetUp_Cyclic;
515: svd->ops->setfromoptions = SVDSetFromOptions_Cyclic;
516: svd->ops->destroy = SVDDestroy_Cyclic;
517: svd->ops->reset = SVDReset_Cyclic;
518: svd->ops->view = SVDView_Cyclic;
519: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetEPS_C",SVDCyclicSetEPS_Cyclic);
520: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetEPS_C",SVDCyclicGetEPS_Cyclic);
521: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicSetExplicitMatrix_C",SVDCyclicSetExplicitMatrix_Cyclic);
522: PetscObjectComposeFunction((PetscObject)svd,"SVDCyclicGetExplicitMatrix_C",SVDCyclicGetExplicitMatrix_Cyclic);
523: return(0);
524: }