Actual source code: svd.c
petsc-3.3-p7 2013-05-11
2: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
3: #include <petscblaslapack.h>
5: /*
6: Private context (data structure) for the SVD preconditioner.
7: */
8: typedef struct {
9: Vec diag,work;
10: Mat A,U,V;
11: PetscInt nzero;
12: PetscReal zerosing; /* measure of smallest singular value treated as nonzero */
13: PetscInt essrank; /* essential rank of operator */
14: PetscViewer monitor;
15: } PC_SVD;
18: /* -------------------------------------------------------------------------- */
19: /*
20: PCSetUp_SVD - Prepares for the use of the SVD preconditioner
21: by setting data structures and options.
23: Input Parameter:
24: . pc - the preconditioner context
26: Application Interface Routine: PCSetUp()
28: Notes:
29: The interface routine PCSetUp() is not usually called directly by
30: the user, but instead is called by PCApply() if necessary.
31: */
34: static PetscErrorCode PCSetUp_SVD(PC pc)
35: {
36: #if defined(PETSC_MISSING_LAPACK_GESVD)
37: SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_SUP,"GESVD - Lapack routine is unavailable\nNot able to provide singular value estimates.");
38: #else
39: PC_SVD *jac = (PC_SVD*)pc->data;
41: PetscScalar *a,*u,*v,*d,*work;
42: PetscBLASInt nb,lwork;
43: PetscInt i,n;
46: if (!jac->diag) {
47: /* assume square matrices */
48: MatGetVecs(pc->pmat,&jac->diag,&jac->work);
49: }
50: MatDestroy(&jac->A);
51: MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);
52: if (!jac->U) {
53: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);
54: MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->V);
55: }
56: MatGetSize(pc->pmat,&n,PETSC_NULL);
57: nb = PetscBLASIntCast(n);
58: lwork = 5*nb;
59: PetscMalloc(lwork*sizeof(PetscScalar),&work);
60: MatGetArray(jac->A,&a);
61: MatGetArray(jac->U,&u);
62: MatGetArray(jac->V,&v);
63: VecGetArray(jac->diag,&d);
64: #if !defined(PETSC_USE_COMPLEX)
65: {
66: PetscBLASInt lierr;
67: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
68: LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr);
69: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
70: PetscFPTrapPop();
71: }
72: #else
73: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not coded for complex");
74: #endif
75: MatRestoreArray(jac->A,&a);
76: MatRestoreArray(jac->U,&u);
77: MatRestoreArray(jac->V,&v);
78: for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
79: jac->nzero = n-1-i;
80: if (jac->monitor) {
81: PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);
82: PetscViewerASCIIPrintf(jac->monitor," SVD: condition number %14.12e, %D of %D singular values are (nearly) zero\n",(double)PetscRealPart(d[0]/d[n-1]),jac->nzero,n);
83: if (n >= 10) { /* print 5 smallest and 5 largest */
84: PetscViewerASCIIPrintf(jac->monitor," SVD: smallest singular values: %14.12e %14.12e %14.12e %14.12e %14.12e\n",(double)PetscRealPart(d[n-1]),(double)PetscRealPart(d[n-2]),(double)PetscRealPart(d[n-3]),(double)PetscRealPart(d[n-4]),(double)PetscRealPart(d[n-5]));
85: PetscViewerASCIIPrintf(jac->monitor," SVD: largest singular values : %14.12e %14.12e %14.12e %14.12e %14.12e\n",(double)PetscRealPart(d[4]),(double)PetscRealPart(d[3]),(double)PetscRealPart(d[2]),(double)PetscRealPart(d[1]),(double)PetscRealPart(d[0]));
86: } else { /* print all singular values */
87: char buf[256],*p;
88: size_t left = sizeof buf,used;
89: PetscInt thisline;
90: for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
91: PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));
92: left -= used;
93: p += used;
94: if (thisline > 4 || i==0) {
95: PetscViewerASCIIPrintf(jac->monitor," SVD: singular values:%s\n",buf);
96: p = buf;
97: thisline = 0;
98: }
99: }
100: }
101: PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);
102: }
103: PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));
104: for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
105: for (; i<n; i++) d[i] = 0.0;
106: if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
107: PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);
108: VecRestoreArray(jac->diag,&d);
109: #if defined(foo)
110: {
111: PetscViewer viewer;
112: PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);
113: MatView(jac->A,viewer);
114: MatView(jac->U,viewer);
115: MatView(jac->V,viewer);
116: VecView(jac->diag,viewer);
117: PetscViewerDestroy(viewer);
118: }
119: #endif
120: PetscFree(work);
121: return(0);
122: #endif
123: }
125: /* -------------------------------------------------------------------------- */
126: /*
127: PCApply_SVD - Applies the SVD preconditioner to a vector.
129: Input Parameters:
130: . pc - the preconditioner context
131: . x - input vector
133: Output Parameter:
134: . y - output vector
136: Application Interface Routine: PCApply()
137: */
140: static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
141: {
142: PC_SVD *jac = (PC_SVD*)pc->data;
143: Vec work = jac->work;
147: MatMultTranspose(jac->U,x,work);
148: VecPointwiseMult(work,work,jac->diag);
149: MatMultTranspose(jac->V,work,y);
150: return(0);
151: }
155: static PetscErrorCode PCReset_SVD(PC pc)
156: {
157: PC_SVD *jac = (PC_SVD*)pc->data;
161: MatDestroy(&jac->A);
162: MatDestroy(&jac->U);
163: MatDestroy(&jac->V);
164: VecDestroy(&jac->diag);
165: VecDestroy(&jac->work);
166: return(0);
167: }
169: /* -------------------------------------------------------------------------- */
170: /*
171: PCDestroy_SVD - Destroys the private context for the SVD preconditioner
172: that was created with PCCreate_SVD().
174: Input Parameter:
175: . pc - the preconditioner context
177: Application Interface Routine: PCDestroy()
178: */
181: static PetscErrorCode PCDestroy_SVD(PC pc)
182: {
183: PC_SVD *jac = (PC_SVD*)pc->data;
187: PCReset_SVD(pc);
188: PetscViewerDestroy(&jac->monitor);
189: PetscFree(pc->data);
190: return(0);
191: }
195: static PetscErrorCode PCSetFromOptions_SVD(PC pc)
196: {
198: PC_SVD *jac = (PC_SVD*)pc->data;
199: PetscBool flg,set;
202: PetscOptionsHead("SVD options");
203: PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,PETSC_NULL);
204: PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,PETSC_NULL);
205: PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor?PETSC_TRUE:PETSC_FALSE,&flg,&set);
206: if (set) { /* Should make PCSVDSetMonitor() */
207: if (flg && !jac->monitor) {
208: PetscViewerASCIIOpen(((PetscObject)pc)->comm,"stdout",&jac->monitor);
209: } else if (!flg) {
210: PetscViewerDestroy(&jac->monitor);
211: }
212: }
213: PetscOptionsTail();
214: return(0);
215: }
217: /* -------------------------------------------------------------------------- */
218: /*
219: PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
220: and sets this as the private data within the generic preconditioning
221: context, PC, that was created within PCCreate().
223: Input Parameter:
224: . pc - the preconditioner context
226: Application Interface Routine: PCCreate()
227: */
229: /*MC
230: PCSVD - Use pseudo inverse defined by SVD of operator
232: Level: advanced
234: Concepts: SVD
236: Options Database:
237: - -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
238: + -pc_svd_monitor Print information on the extreme singular values of the operator
240: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
241: M*/
243: EXTERN_C_BEGIN
246: PetscErrorCode PCCreate_SVD(PC pc)
247: {
248: PC_SVD *jac;
252: /*
253: Creates the private data structure for this preconditioner and
254: attach it to the PC object.
255: */
256: PetscNewLog(pc,PC_SVD,&jac);
257: jac->zerosing = 1.e-12;
258: pc->data = (void*)jac;
260: /*
261: Set the pointers for the functions that are provided above.
262: Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
263: are called, they will automatically call these functions. Note we
264: choose not to provide a couple of these functions since they are
265: not needed.
266: */
267: pc->ops->apply = PCApply_SVD;
268: pc->ops->applytranspose = PCApply_SVD;
269: pc->ops->setup = PCSetUp_SVD;
270: pc->ops->reset = PCReset_SVD;
271: pc->ops->destroy = PCDestroy_SVD;
272: pc->ops->setfromoptions = PCSetFromOptions_SVD;
273: pc->ops->view = 0;
274: pc->ops->applyrichardson = 0;
275: return(0);
276: }
277: EXTERN_C_END