1
2 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
3 #include <petscblaslapack.h>
4
5 /*
6 Private context (data structure) for the SVD preconditioner.
7 */
8 typedef struct {
9 Vec diag,work;
10 Mat A,U,Vt;
11 PetscInt nzero;
12 PetscReal zerosing; /* measure of smallest singular value treated as nonzero */
13 PetscInt essrank; /* essential rank of operator */
14 VecScatter left2red,right2red;
15 Vec leftred,rightred;
16 PetscViewer monitor;
17 } PC_SVD;
18
19 typedef enum {READ=1, WRITE=2, READ_WRITE=3} AccessMode;
20
21 /* -------------------------------------------------------------------------- */
22 /*
23 PCSetUp_SVD - Prepares for the use of the SVD preconditioner
24 by setting data structures and options.
25
26 Input Parameter:
27 . pc - the preconditioner context
28
29 Application Interface Routine: PCSetUp()
30
31 Notes:
32 The interface routine PCSetUp() is not usually called directly by
33 the user, but instead is called by PCApply() if necessary.
34 */
PCSetUp_SVD(PC pc)35 static PetscErrorCode PCSetUp_SVD(PC pc)
36 {
37 PC_SVD *jac = (PC_SVD*)pc->data;
38 PetscErrorCode ierr;
39 PetscScalar *a,*u,*v,*d,*work;
40 PetscBLASInt nb,lwork;
41 PetscInt i,n;
42 PetscMPIInt size;
43
44 PetscFunctionBegin;
45 ierr = MatDestroy(&jac->A);CHKERRQ(ierr);
46 ierr = MPI_Comm_size(((PetscObject)pc->pmat)->comm,&size);CHKERRQ(ierr);
47 if (size > 1) {
48 Mat redmat;
49
50 ierr = MatCreateRedundantMatrix(pc->pmat,size,PETSC_COMM_SELF,MAT_INITIAL_MATRIX,&redmat);CHKERRQ(ierr);
51 ierr = MatConvert(redmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);CHKERRQ(ierr);
52 ierr = MatDestroy(&redmat);CHKERRQ(ierr);
53 } else {
54 ierr = MatConvert(pc->pmat,MATSEQDENSE,MAT_INITIAL_MATRIX,&jac->A);CHKERRQ(ierr);
55 }
56 if (!jac->diag) { /* assume square matrices */
57 ierr = MatCreateVecs(jac->A,&jac->diag,&jac->work);CHKERRQ(ierr);
58 }
59 if (!jac->U) {
60 ierr = MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->U);CHKERRQ(ierr);
61 ierr = MatDuplicate(jac->A,MAT_DO_NOT_COPY_VALUES,&jac->Vt);CHKERRQ(ierr);
62 }
63 ierr = MatGetSize(jac->A,&n,NULL);CHKERRQ(ierr);
64 if (!n) {
65 ierr = PetscInfo(pc,"Matrix has zero rows, skipping svd\n");CHKERRQ(ierr);
66 PetscFunctionReturn(0);
67 }
68 ierr = PetscBLASIntCast(n,&nb);CHKERRQ(ierr);
69 lwork = 5*nb;
70 ierr = PetscMalloc1(lwork,&work);CHKERRQ(ierr);
71 ierr = MatDenseGetArray(jac->A,&a);CHKERRQ(ierr);
72 ierr = MatDenseGetArray(jac->U,&u);CHKERRQ(ierr);
73 ierr = MatDenseGetArray(jac->Vt,&v);CHKERRQ(ierr);
74 ierr = VecGetArray(jac->diag,&d);CHKERRQ(ierr);
75 #if !defined(PETSC_USE_COMPLEX)
76 {
77 PetscBLASInt lierr;
78 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
79 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,d,u,&nb,v,&nb,work,&lwork,&lierr));
80 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
81 ierr = PetscFPTrapPop();CHKERRQ(ierr);
82 }
83 #else
84 {
85 PetscBLASInt lierr;
86 PetscReal *rwork,*dd;
87 ierr = PetscMalloc1(5*nb,&rwork);CHKERRQ(ierr);
88 ierr = PetscMalloc1(nb,&dd);CHKERRQ(ierr);
89 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr);
90 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("A","A",&nb,&nb,a,&nb,dd,u,&nb,v,&nb,work,&lwork,rwork,&lierr));
91 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"gesv() error %d",lierr);
92 ierr = PetscFree(rwork);CHKERRQ(ierr);
93 for (i=0; i<n; i++) d[i] = dd[i];
94 ierr = PetscFree(dd);CHKERRQ(ierr);
95 ierr = PetscFPTrapPop();CHKERRQ(ierr);
96 }
97 #endif
98 ierr = MatDenseRestoreArray(jac->A,&a);CHKERRQ(ierr);
99 ierr = MatDenseRestoreArray(jac->U,&u);CHKERRQ(ierr);
100 ierr = MatDenseRestoreArray(jac->Vt,&v);CHKERRQ(ierr);
101 for (i=n-1; i>=0; i--) if (PetscRealPart(d[i]) > jac->zerosing) break;
102 jac->nzero = n-1-i;
103 if (jac->monitor) {
104 ierr = PetscViewerASCIIAddTab(jac->monitor,((PetscObject)pc)->tablevel);CHKERRQ(ierr);
105 ierr = 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);CHKERRQ(ierr);
106 if (n >= 10) { /* print 5 smallest and 5 largest */
107 ierr = 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]));CHKERRQ(ierr);
108 ierr = 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]));CHKERRQ(ierr);
109 } else { /* print all singular values */
110 char buf[256],*p;
111 size_t left = sizeof(buf),used;
112 PetscInt thisline;
113 for (p=buf,i=n-1,thisline=1; i>=0; i--,thisline++) {
114 ierr = PetscSNPrintfCount(p,left," %14.12e",&used,(double)PetscRealPart(d[i]));CHKERRQ(ierr);
115 left -= used;
116 p += used;
117 if (thisline > 4 || i==0) {
118 ierr = PetscViewerASCIIPrintf(jac->monitor," SVD: singular values:%s\n",buf);CHKERRQ(ierr);
119 p = buf;
120 thisline = 0;
121 }
122 }
123 }
124 ierr = PetscViewerASCIISubtractTab(jac->monitor,((PetscObject)pc)->tablevel);CHKERRQ(ierr);
125 }
126 ierr = PetscInfo2(pc,"Largest and smallest singular values %14.12e %14.12e\n",(double)PetscRealPart(d[0]),(double)PetscRealPart(d[n-1]));CHKERRQ(ierr);
127 for (i=0; i<n-jac->nzero; i++) d[i] = 1.0/d[i];
128 for (; i<n; i++) d[i] = 0.0;
129 if (jac->essrank > 0) for (i=0; i<n-jac->nzero-jac->essrank; i++) d[i] = 0.0; /* Skip all but essrank eigenvalues */
130 ierr = PetscInfo1(pc,"Number of zero or nearly singular values %D\n",jac->nzero);CHKERRQ(ierr);
131 ierr = VecRestoreArray(jac->diag,&d);CHKERRQ(ierr);
132 #if defined(foo)
133 {
134 PetscViewer viewer;
135 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"joe",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
136 ierr = MatView(jac->A,viewer);CHKERRQ(ierr);
137 ierr = MatView(jac->U,viewer);CHKERRQ(ierr);
138 ierr = MatView(jac->Vt,viewer);CHKERRQ(ierr);
139 ierr = VecView(jac->diag,viewer);CHKERRQ(ierr);
140 ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
141 }
142 #endif
143 ierr = PetscFree(work);CHKERRQ(ierr);
144 PetscFunctionReturn(0);
145 }
146
PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec * xred)147 static PetscErrorCode PCSVDGetVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
148 {
149 PC_SVD *jac = (PC_SVD*)pc->data;
150 PetscErrorCode ierr;
151 PetscMPIInt size;
152
153 PetscFunctionBegin;
154 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
155 *xred = NULL;
156 switch (side) {
157 case PC_LEFT:
158 if (size == 1) *xred = x;
159 else {
160 if (!jac->left2red) {ierr = VecScatterCreateToAll(x,&jac->left2red,&jac->leftred);CHKERRQ(ierr);}
161 if (amode & READ) {
162 ierr = VecScatterBegin(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
163 ierr = VecScatterEnd(jac->left2red,x,jac->leftred,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
164 }
165 *xred = jac->leftred;
166 }
167 break;
168 case PC_RIGHT:
169 if (size == 1) *xred = x;
170 else {
171 if (!jac->right2red) {ierr = VecScatterCreateToAll(x,&jac->right2red,&jac->rightred);CHKERRQ(ierr);}
172 if (amode & READ) {
173 ierr = VecScatterBegin(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
174 ierr = VecScatterEnd(jac->right2red,x,jac->rightred,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
175 }
176 *xred = jac->rightred;
177 }
178 break;
179 default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
180 }
181 PetscFunctionReturn(0);
182 }
183
PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec * xred)184 static PetscErrorCode PCSVDRestoreVec(PC pc,PCSide side,AccessMode amode,Vec x,Vec *xred)
185 {
186 PC_SVD *jac = (PC_SVD*)pc->data;
187 PetscErrorCode ierr;
188 PetscMPIInt size;
189
190 PetscFunctionBegin;
191 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);CHKERRQ(ierr);
192 switch (side) {
193 case PC_LEFT:
194 if (size != 1 && amode & WRITE) {
195 ierr = VecScatterBegin(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
196 ierr = VecScatterEnd(jac->left2red,jac->leftred,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
197 }
198 break;
199 case PC_RIGHT:
200 if (size != 1 && amode & WRITE) {
201 ierr = VecScatterBegin(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
202 ierr = VecScatterEnd(jac->right2red,jac->rightred,x,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr);
203 }
204 break;
205 default: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"Side must be LEFT or RIGHT");
206 }
207 *xred = NULL;
208 PetscFunctionReturn(0);
209 }
210
211 /* -------------------------------------------------------------------------- */
212 /*
213 PCApply_SVD - Applies the SVD preconditioner to a vector.
214
215 Input Parameters:
216 . pc - the preconditioner context
217 . x - input vector
218
219 Output Parameter:
220 . y - output vector
221
222 Application Interface Routine: PCApply()
223 */
PCApply_SVD(PC pc,Vec x,Vec y)224 static PetscErrorCode PCApply_SVD(PC pc,Vec x,Vec y)
225 {
226 PC_SVD *jac = (PC_SVD*)pc->data;
227 Vec work = jac->work,xred,yred;
228 PetscErrorCode ierr;
229
230 PetscFunctionBegin;
231 ierr = PCSVDGetVec(pc,PC_RIGHT,READ,x,&xred);CHKERRQ(ierr);
232 ierr = PCSVDGetVec(pc,PC_LEFT,WRITE,y,&yred);CHKERRQ(ierr);
233 #if !defined(PETSC_USE_COMPLEX)
234 ierr = MatMultTranspose(jac->U,xred,work);CHKERRQ(ierr);
235 #else
236 ierr = MatMultHermitianTranspose(jac->U,xred,work);CHKERRQ(ierr);
237 #endif
238 ierr = VecPointwiseMult(work,work,jac->diag);CHKERRQ(ierr);
239 #if !defined(PETSC_USE_COMPLEX)
240 ierr = MatMultTranspose(jac->Vt,work,yred);CHKERRQ(ierr);
241 #else
242 ierr = MatMultHermitianTranspose(jac->Vt,work,yred);CHKERRQ(ierr);
243 #endif
244 ierr = PCSVDRestoreVec(pc,PC_RIGHT,READ,x,&xred);CHKERRQ(ierr);
245 ierr = PCSVDRestoreVec(pc,PC_LEFT,WRITE,y,&yred);CHKERRQ(ierr);
246 PetscFunctionReturn(0);
247 }
248
PCApplyTranspose_SVD(PC pc,Vec x,Vec y)249 static PetscErrorCode PCApplyTranspose_SVD(PC pc,Vec x,Vec y)
250 {
251 PC_SVD *jac = (PC_SVD*)pc->data;
252 Vec work = jac->work,xred,yred;
253 PetscErrorCode ierr;
254
255 PetscFunctionBegin;
256 ierr = PCSVDGetVec(pc,PC_LEFT,READ,x,&xred);CHKERRQ(ierr);
257 ierr = PCSVDGetVec(pc,PC_RIGHT,WRITE,y,&yred);CHKERRQ(ierr);
258 ierr = MatMult(jac->Vt,xred,work);CHKERRQ(ierr);
259 ierr = VecPointwiseMult(work,work,jac->diag);CHKERRQ(ierr);
260 ierr = MatMult(jac->U,work,yred);CHKERRQ(ierr);
261 ierr = PCSVDRestoreVec(pc,PC_LEFT,READ,x,&xred);CHKERRQ(ierr);
262 ierr = PCSVDRestoreVec(pc,PC_RIGHT,WRITE,y,&yred);CHKERRQ(ierr);
263 PetscFunctionReturn(0);
264 }
265
PCReset_SVD(PC pc)266 static PetscErrorCode PCReset_SVD(PC pc)
267 {
268 PC_SVD *jac = (PC_SVD*)pc->data;
269 PetscErrorCode ierr;
270
271 PetscFunctionBegin;
272 ierr = MatDestroy(&jac->A);CHKERRQ(ierr);
273 ierr = MatDestroy(&jac->U);CHKERRQ(ierr);
274 ierr = MatDestroy(&jac->Vt);CHKERRQ(ierr);
275 ierr = VecDestroy(&jac->diag);CHKERRQ(ierr);
276 ierr = VecDestroy(&jac->work);CHKERRQ(ierr);
277 ierr = VecScatterDestroy(&jac->right2red);CHKERRQ(ierr);
278 ierr = VecScatterDestroy(&jac->left2red);CHKERRQ(ierr);
279 ierr = VecDestroy(&jac->rightred);CHKERRQ(ierr);
280 ierr = VecDestroy(&jac->leftred);CHKERRQ(ierr);
281 PetscFunctionReturn(0);
282 }
283
284 /* -------------------------------------------------------------------------- */
285 /*
286 PCDestroy_SVD - Destroys the private context for the SVD preconditioner
287 that was created with PCCreate_SVD().
288
289 Input Parameter:
290 . pc - the preconditioner context
291
292 Application Interface Routine: PCDestroy()
293 */
PCDestroy_SVD(PC pc)294 static PetscErrorCode PCDestroy_SVD(PC pc)
295 {
296 PC_SVD *jac = (PC_SVD*)pc->data;
297 PetscErrorCode ierr;
298
299 PetscFunctionBegin;
300 ierr = PCReset_SVD(pc);CHKERRQ(ierr);
301 ierr = PetscViewerDestroy(&jac->monitor);CHKERRQ(ierr);
302 ierr = PetscFree(pc->data);CHKERRQ(ierr);
303 PetscFunctionReturn(0);
304 }
305
PCSetFromOptions_SVD(PetscOptionItems * PetscOptionsObject,PC pc)306 static PetscErrorCode PCSetFromOptions_SVD(PetscOptionItems *PetscOptionsObject,PC pc)
307 {
308 PetscErrorCode ierr;
309 PC_SVD *jac = (PC_SVD*)pc->data;
310 PetscBool flg,set;
311
312 PetscFunctionBegin;
313 ierr = PetscOptionsHead(PetscOptionsObject,"SVD options");CHKERRQ(ierr);
314 ierr = PetscOptionsReal("-pc_svd_zero_sing","Singular values smaller than this treated as zero","None",jac->zerosing,&jac->zerosing,NULL);CHKERRQ(ierr);
315 ierr = PetscOptionsInt("-pc_svd_ess_rank","Essential rank of operator (0 to use entire operator)","None",jac->essrank,&jac->essrank,NULL);CHKERRQ(ierr);
316 ierr = PetscOptionsBool("-pc_svd_monitor","Monitor the conditioning, and extremal singular values","None",jac->monitor ? PETSC_TRUE : PETSC_FALSE,&flg,&set);CHKERRQ(ierr);
317 if (set) { /* Should make PCSVDSetMonitor() */
318 if (flg && !jac->monitor) {
319 ierr = PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),"stdout",&jac->monitor);CHKERRQ(ierr);
320 } else if (!flg) {
321 ierr = PetscViewerDestroy(&jac->monitor);CHKERRQ(ierr);
322 }
323 }
324 ierr = PetscOptionsTail();CHKERRQ(ierr);
325 PetscFunctionReturn(0);
326 }
327
PCView_SVD(PC pc,PetscViewer viewer)328 static PetscErrorCode PCView_SVD(PC pc,PetscViewer viewer)
329 {
330 PC_SVD *svd = (PC_SVD*)pc->data;
331 PetscErrorCode ierr;
332 PetscBool iascii;
333
334 PetscFunctionBegin;
335 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
336 if (iascii) {
337 ierr = PetscViewerASCIIPrintf(viewer," All singular values smaller than %g treated as zero\n",(double)svd->zerosing);CHKERRQ(ierr);
338 ierr = PetscViewerASCIIPrintf(viewer," Provided essential rank of the matrix %D (all other eigenvalues are zeroed)\n",svd->essrank);CHKERRQ(ierr);
339 }
340 PetscFunctionReturn(0);
341 }
342 /* -------------------------------------------------------------------------- */
343 /*
344 PCCreate_SVD - Creates a SVD preconditioner context, PC_SVD,
345 and sets this as the private data within the generic preconditioning
346 context, PC, that was created within PCCreate().
347
348 Input Parameter:
349 . pc - the preconditioner context
350
351 Application Interface Routine: PCCreate()
352 */
353
354 /*MC
355 PCSVD - Use pseudo inverse defined by SVD of operator
356
357 Level: advanced
358
359 Options Database:
360 + -pc_svd_zero_sing <rtol> Singular values smaller than this are treated as zero
361 - -pc_svd_monitor Print information on the extreme singular values of the operator
362
363 Developer Note:
364 This implementation automatically creates a redundant copy of the
365 matrix on each process and uses a sequential SVD solve. Why does it do this instead
366 of using the composable PCREDUNDANT object?
367
368 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
369 M*/
370
PCCreate_SVD(PC pc)371 PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC pc)
372 {
373 PC_SVD *jac;
374 PetscErrorCode ierr;
375
376 PetscFunctionBegin;
377 /*
378 Creates the private data structure for this preconditioner and
379 attach it to the PC object.
380 */
381 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr);
382 jac->zerosing = 1.e-12;
383 pc->data = (void*)jac;
384
385 /*
386 Set the pointers for the functions that are provided above.
387 Now when the user-level routines (such as PCApply(), PCDestroy(), etc.)
388 are called, they will automatically call these functions. Note we
389 choose not to provide a couple of these functions since they are
390 not needed.
391 */
392 pc->ops->apply = PCApply_SVD;
393 pc->ops->applytranspose = PCApplyTranspose_SVD;
394 pc->ops->setup = PCSetUp_SVD;
395 pc->ops->reset = PCReset_SVD;
396 pc->ops->destroy = PCDestroy_SVD;
397 pc->ops->setfromoptions = PCSetFromOptions_SVD;
398 pc->ops->view = PCView_SVD;
399 pc->ops->applyrichardson = NULL;
400 PetscFunctionReturn(0);
401 }
402
403