Actual source code: svdsolve.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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: SVD routines related to the solution process
12: */
14: #include <slepc/private/svdimpl.h>
15: #include <slepc/private/bvimpl.h>
17: /*
18: SVDComputeVectors_Left - Compute left singular vectors as U=A*V.
19: Only done if the leftbasis flag is false. Assumes V is available.
20: */
21: PetscErrorCode SVDComputeVectors_Left(SVD svd)
22: {
23: Vec tl,omega2,u,v,w;
24: PetscInt i,oldsize;
25: VecType vtype;
26: const PetscScalar* varray;
28: PetscFunctionBegin;
29: if (!svd->leftbasis) {
30: /* generate left singular vectors on U */
31: if (!svd->U) PetscCall(SVDGetBV(svd,NULL,&svd->U));
32: PetscCall(BVGetSizes(svd->U,NULL,NULL,&oldsize));
33: if (!oldsize) {
34: if (!((PetscObject)svd->U)->type_name) PetscCall(BVSetType(svd->U,((PetscObject)svd->V)->type_name));
35: PetscCall(MatCreateVecsEmpty(svd->A,NULL,&tl));
36: PetscCall(BVSetSizesFromVec(svd->U,tl,svd->ncv));
37: PetscCall(VecDestroy(&tl));
38: }
39: PetscCall(BVSetActiveColumns(svd->V,0,svd->nconv));
40: PetscCall(BVSetActiveColumns(svd->U,0,svd->nconv));
41: if (!svd->ishyperbolic) PetscCall(BVMatMult(svd->V,svd->A,svd->U));
42: else if (svd->swapped) { /* compute right singular vectors as V=A'*Omega*U */
43: PetscCall(MatCreateVecs(svd->A,&w,NULL));
44: for (i=0;i<svd->nconv;i++) {
45: PetscCall(BVGetColumn(svd->V,i,&v));
46: PetscCall(BVGetColumn(svd->U,i,&u));
47: PetscCall(VecPointwiseMult(w,v,svd->omega));
48: PetscCall(MatMult(svd->A,w,u));
49: PetscCall(BVRestoreColumn(svd->V,i,&v));
50: PetscCall(BVRestoreColumn(svd->U,i,&u));
51: }
52: PetscCall(VecDestroy(&w));
53: } else { /* compute left singular vectors as usual U=A*V, and set-up Omega-orthogonalization of U */
54: PetscCall(BV_SetMatrixDiagonal(svd->U,svd->omega,svd->A));
55: PetscCall(BVMatMult(svd->V,svd->A,svd->U));
56: }
57: PetscCall(BVOrthogonalize(svd->U,NULL));
58: if (svd->ishyperbolic && !svd->swapped) { /* store signature after Omega-orthogonalization */
59: PetscCall(MatGetVecType(svd->A,&vtype));
60: PetscCall(VecCreate(PETSC_COMM_SELF,&omega2));
61: PetscCall(VecSetSizes(omega2,svd->nconv,svd->nconv));
62: PetscCall(VecSetType(omega2,vtype));
63: PetscCall(BVGetSignature(svd->U,omega2));
64: PetscCall(VecGetArrayRead(omega2,&varray));
65: for (i=0;i<svd->nconv;i++) {
66: svd->sign[i] = PetscRealPart(varray[i]);
67: if (PetscRealPart(varray[i])<0.0) PetscCall(BVScaleColumn(svd->U,i,-1.0));
68: }
69: PetscCall(VecRestoreArrayRead(omega2,&varray));
70: PetscCall(VecDestroy(&omega2));
71: }
72: }
73: PetscFunctionReturn(PETSC_SUCCESS);
74: }
76: PetscErrorCode SVDComputeVectors(SVD svd)
77: {
78: PetscFunctionBegin;
79: SVDCheckSolved(svd,1);
80: if (svd->state==SVD_STATE_SOLVED) PetscTryTypeMethod(svd,computevectors);
81: svd->state = SVD_STATE_VECTORS;
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: /*@
86: SVDSolve - Solves the singular value problem.
88: Collective
90: Input Parameter:
91: . svd - singular value solver context obtained from SVDCreate()
93: Options Database Keys:
94: + -svd_view - print information about the solver used
95: . -svd_view_mat0 - view the first matrix (A)
96: . -svd_view_mat1 - view the second matrix (B)
97: . -svd_view_signature - view the signature matrix (omega)
98: . -svd_view_vectors - view the computed singular vectors
99: . -svd_view_values - view the computed singular values
100: . -svd_converged_reason - print reason for convergence, and number of iterations
101: . -svd_error_absolute - print absolute errors of each singular triplet
102: . -svd_error_relative - print relative errors of each singular triplet
103: - -svd_error_norm - print errors relative to the matrix norms of each singular triplet
105: Notes:
106: All the command-line options listed above admit an optional argument specifying
107: the viewer type and options. For instance, use '-svd_view_mat0 binary:amatrix.bin'
108: to save the A matrix to a binary file, '-svd_view_values draw' to draw the computed
109: singular values graphically, or '-svd_error_relative :myerr.m:ascii_matlab' to save
110: the errors in a file that can be executed in Matlab.
112: Level: beginner
114: .seealso: `SVDCreate()`, `SVDSetUp()`, `SVDDestroy()`
115: @*/
116: PetscErrorCode SVDSolve(SVD svd)
117: {
118: PetscInt i,m,n,*workperm;
120: PetscFunctionBegin;
122: if (svd->state>=SVD_STATE_SOLVED) PetscFunctionReturn(PETSC_SUCCESS);
123: PetscCall(PetscLogEventBegin(SVD_Solve,svd,0,0,0));
125: /* call setup */
126: PetscCall(SVDSetUp(svd));
128: /* safeguard for matrices with zero rows or columns */
129: PetscCall(MatGetSize(svd->OP,&m,&n));
130: if (m == 0 || n == 0) {
131: svd->nconv = 0;
132: svd->reason = SVD_CONVERGED_TOL;
133: svd->state = SVD_STATE_SOLVED;
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: svd->its = 0;
138: svd->nconv = 0;
139: for (i=0;i<svd->ncv;i++) {
140: svd->sigma[i] = 0.0;
141: svd->errest[i] = 0.0;
142: svd->perm[i] = i;
143: }
144: PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view_pre"));
146: switch (svd->problem_type) {
147: case SVD_STANDARD:
148: PetscUseTypeMethod(svd,solve);
149: break;
150: case SVD_GENERALIZED:
151: PetscUseTypeMethod(svd,solveg);
152: break;
153: case SVD_HYPERBOLIC:
154: PetscUseTypeMethod(svd,solveh);
155: break;
156: }
157: svd->state = SVD_STATE_SOLVED;
159: /* sort singular triplets */
160: if (svd->which == SVD_SMALLEST) PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,svd->perm));
161: else {
162: PetscCall(PetscMalloc1(svd->nconv,&workperm));
163: for (i=0;i<svd->nconv;i++) workperm[i] = i;
164: PetscCall(PetscSortRealWithPermutation(svd->nconv,svd->sigma,workperm));
165: for (i=0;i<svd->nconv;i++) svd->perm[i] = workperm[svd->nconv-i-1];
166: PetscCall(PetscFree(workperm));
167: }
168: PetscCall(PetscLogEventEnd(SVD_Solve,svd,0,0,0));
170: /* various viewers */
171: PetscCall(SVDViewFromOptions(svd,NULL,"-svd_view"));
172: PetscCall(SVDConvergedReasonViewFromOptions(svd));
173: PetscCall(SVDErrorViewFromOptions(svd));
174: PetscCall(SVDValuesViewFromOptions(svd));
175: PetscCall(SVDVectorsViewFromOptions(svd));
176: PetscCall(MatViewFromOptions(svd->OP,(PetscObject)svd,"-svd_view_mat0"));
177: if (svd->isgeneralized) PetscCall(MatViewFromOptions(svd->OPb,(PetscObject)svd,"-svd_view_mat1"));
178: if (svd->ishyperbolic) PetscCall(VecViewFromOptions(svd->omega,(PetscObject)svd,"-svd_view_signature"));
180: /* Remove the initial subspaces */
181: svd->nini = 0;
182: svd->ninil = 0;
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: /*@
187: SVDGetIterationNumber - Gets the current iteration number. If the
188: call to SVDSolve() is complete, then it returns the number of iterations
189: carried out by the solution method.
191: Not Collective
193: Input Parameter:
194: . svd - the singular value solver context
196: Output Parameter:
197: . its - number of iterations
199: Note:
200: During the i-th iteration this call returns i-1. If SVDSolve() is
201: complete, then parameter "its" contains either the iteration number at
202: which convergence was successfully reached, or failure was detected.
203: Call SVDGetConvergedReason() to determine if the solver converged or
204: failed and why.
206: Level: intermediate
208: .seealso: `SVDGetConvergedReason()`, `SVDSetTolerances()`
209: @*/
210: PetscErrorCode SVDGetIterationNumber(SVD svd,PetscInt *its)
211: {
212: PetscFunctionBegin;
214: PetscAssertPointer(its,2);
215: *its = svd->its;
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: /*@
220: SVDGetConvergedReason - Gets the reason why the SVDSolve() iteration was
221: stopped.
223: Not Collective
225: Input Parameter:
226: . svd - the singular value solver context
228: Output Parameter:
229: . reason - negative value indicates diverged, positive value converged
230: (see SVDConvergedReason)
232: Options Database Key:
233: . -svd_converged_reason - print the reason to a viewer
235: Notes:
236: Possible values for reason are
237: + SVD_CONVERGED_TOL - converged up to tolerance
238: . SVD_CONVERGED_USER - converged due to a user-defined condition
239: . SVD_CONVERGED_MAXIT - reached the maximum number of iterations with SVD_CONV_MAXIT criterion
240: . SVD_DIVERGED_ITS - required more than max_it iterations to reach convergence
241: . SVD_DIVERGED_BREAKDOWN - generic breakdown in method
242: - SVD_DIVERGED_SYMMETRY_LOST - underlying indefinite eigensolver was not able to keep symmetry
244: Can only be called after the call to SVDSolve() is complete.
246: Level: intermediate
248: .seealso: `SVDSetTolerances()`, `SVDSolve()`, `SVDConvergedReason`
249: @*/
250: PetscErrorCode SVDGetConvergedReason(SVD svd,SVDConvergedReason *reason)
251: {
252: PetscFunctionBegin;
254: PetscAssertPointer(reason,2);
255: SVDCheckSolved(svd,1);
256: *reason = svd->reason;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: /*@
261: SVDGetConverged - Gets the number of converged singular values.
263: Not Collective
265: Input Parameter:
266: . svd - the singular value solver context
268: Output Parameter:
269: . nconv - number of converged singular values
271: Note:
272: This function should be called after SVDSolve() has finished.
274: Level: beginner
276: .seealso: `SVDSetDimensions()`, `SVDSolve()`, `SVDGetSingularTriplet()`
277: @*/
278: PetscErrorCode SVDGetConverged(SVD svd,PetscInt *nconv)
279: {
280: PetscFunctionBegin;
282: PetscAssertPointer(nconv,2);
283: SVDCheckSolved(svd,1);
284: *nconv = svd->nconv;
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*@
289: SVDGetSingularTriplet - Gets the i-th triplet of the singular value decomposition
290: as computed by SVDSolve(). The solution consists in the singular value and its left
291: and right singular vectors.
293: Collective
295: Input Parameters:
296: + svd - singular value solver context
297: - i - index of the solution
299: Output Parameters:
300: + sigma - singular value
301: . u - left singular vector
302: - v - right singular vector
304: Note:
305: Both u or v can be NULL if singular vectors are not required.
306: Otherwise, the caller must provide valid Vec objects, i.e.,
307: they must be created by the calling program with e.g. MatCreateVecs().
309: The index i should be a value between 0 and nconv-1 (see SVDGetConverged()).
310: Singular triplets are indexed according to the ordering criterion established
311: with SVDSetWhichSingularTriplets().
313: In the case of GSVD, the solution consists in three vectors u,v,x that are
314: returned as follows. Vector x is returned in the right singular vector
315: (argument v) and has length equal to the number of columns of A and B.
316: The other two vectors are returned stacked on top of each other [u;v] in
317: the left singular vector argument, with length equal to m+n (number of rows
318: of A plus number of rows of B).
320: Level: beginner
322: .seealso: `SVDSolve()`, `SVDGetConverged()`, `SVDSetWhichSingularTriplets()`
323: @*/
324: PetscErrorCode SVDGetSingularTriplet(SVD svd,PetscInt i,PetscReal *sigma,Vec u,Vec v)
325: {
326: PetscInt M,N;
327: Vec w;
329: PetscFunctionBegin;
332: SVDCheckSolved(svd,1);
335: PetscCheck(i>=0,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index cannot be negative");
336: PetscCheck(i<svd->nconv,PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"The index can be nconv-1 at most, see SVDGetConverged()");
337: if (sigma) *sigma = svd->sigma[svd->perm[i]];
338: if (u || v) {
339: if (!svd->isgeneralized) {
340: PetscCall(MatGetSize(svd->OP,&M,&N));
341: if (M<N) { w = u; u = v; v = w; }
342: }
343: PetscCall(SVDComputeVectors(svd));
344: if (u) PetscCall(BVCopyVec(svd->U,svd->perm[i],u));
345: if (v) PetscCall(BVCopyVec(svd->V,svd->perm[i],v));
346: }
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /*
351: SVDComputeResidualNorms_Standard - Computes the norms of the left and
352: right residuals associated with the i-th computed singular triplet.
354: Input Parameters:
355: sigma - singular value
356: u,v - singular vectors
357: x,y - two work vectors with the same dimensions as u,v
358: */
359: static PetscErrorCode SVDComputeResidualNorms_Standard(SVD svd,PetscReal sigma,Vec u,Vec v,Vec x,Vec y,PetscReal *norm1,PetscReal *norm2)
360: {
361: PetscInt M,N;
363: PetscFunctionBegin;
364: /* norm1 = ||A*v-sigma*u||_2 */
365: if (norm1) {
366: PetscCall(MatMult(svd->OP,v,x));
367: PetscCall(VecAXPY(x,-sigma,u));
368: PetscCall(VecNorm(x,NORM_2,norm1));
369: }
370: /* norm2 = ||A^T*u-sigma*v||_2 */
371: if (norm2) {
372: PetscCall(MatGetSize(svd->OP,&M,&N));
373: if (M<N) PetscCall(MatMult(svd->A,u,y));
374: else PetscCall(MatMult(svd->AT,u,y));
375: PetscCall(VecAXPY(y,-sigma,v));
376: PetscCall(VecNorm(y,NORM_2,norm2));
377: }
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }
381: /*
382: SVDComputeResidualNorms_Generalized - In GSVD, compute the residual norms
383: norm1 = ||s^2*A'*u-c*B'*B*x||_2 and norm2 = ||c^2*B'*v-s*A'*A*x||_2.
385: Input Parameters:
386: sigma - singular value
387: uv - left singular vectors [u;v]
388: x - right singular vector
389: y,z - two work vectors with the same dimension as x
390: */
391: static PetscErrorCode SVDComputeResidualNorms_Generalized(SVD svd,PetscReal sigma,Vec uv,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
392: {
393: Vec u,v,ax,bx,nest,aux[2];
394: PetscReal c,s;
396: PetscFunctionBegin;
397: PetscCall(MatCreateVecs(svd->OP,NULL,&u));
398: PetscCall(MatCreateVecs(svd->OPb,NULL,&v));
399: aux[0] = u;
400: aux[1] = v;
401: PetscCall(VecCreateNest(PetscObjectComm((PetscObject)svd),2,NULL,aux,&nest));
402: PetscCall(VecCopy(uv,nest));
404: s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
405: c = sigma*s;
407: /* norm1 = ||s^2*A'*u-c*B'*B*x||_2 */
408: if (norm1) {
409: PetscCall(VecDuplicate(v,&bx));
410: PetscCall(MatMultHermitianTranspose(svd->OP,u,z));
411: PetscCall(MatMult(svd->OPb,x,bx));
412: PetscCall(MatMultHermitianTranspose(svd->OPb,bx,y));
413: PetscCall(VecAXPBY(y,s*s,-c,z));
414: PetscCall(VecNorm(y,NORM_2,norm1));
415: PetscCall(VecDestroy(&bx));
416: }
417: /* norm2 = ||c^2*B'*v-s*A'*A*x||_2 */
418: if (norm2) {
419: PetscCall(VecDuplicate(u,&ax));
420: PetscCall(MatMultHermitianTranspose(svd->OPb,v,z));
421: PetscCall(MatMult(svd->OP,x,ax));
422: PetscCall(MatMultHermitianTranspose(svd->OP,ax,y));
423: PetscCall(VecAXPBY(y,c*c,-s,z));
424: PetscCall(VecNorm(y,NORM_2,norm2));
425: PetscCall(VecDestroy(&ax));
426: }
428: PetscCall(VecDestroy(&v));
429: PetscCall(VecDestroy(&u));
430: PetscCall(VecDestroy(&nest));
431: PetscFunctionReturn(PETSC_SUCCESS);
432: }
434: /*
435: SVDComputeResidualNorms_Hyperbolic - Computes the norms of the left and
436: right residuals associated with the i-th computed singular triplet.
438: Input Parameters:
439: sigma - singular value
440: sign - corresponding element of the signature Omega2
441: u,v - singular vectors
442: x,y,z - three work vectors with the same dimensions as u,v,u
443: */
444: static PetscErrorCode SVDComputeResidualNorms_Hyperbolic(SVD svd,PetscReal sigma,PetscReal sign,Vec u,Vec v,Vec x,Vec y,Vec z,PetscReal *norm1,PetscReal *norm2)
445: {
446: PetscInt M,N;
448: PetscFunctionBegin;
449: /* norm1 = ||A*v-sigma*u||_2 */
450: if (norm1) {
451: PetscCall(MatMult(svd->OP,v,x));
452: PetscCall(VecAXPY(x,-sigma,u));
453: PetscCall(VecNorm(x,NORM_2,norm1));
454: }
455: /* norm2 = ||A^T*Omega*u-sigma*sign*v||_2 */
456: if (norm2) {
457: PetscCall(MatGetSize(svd->OP,&M,&N));
458: PetscCall(VecPointwiseMult(z,u,svd->omega));
459: if (M<N) PetscCall(MatMult(svd->A,z,y));
460: else PetscCall(MatMult(svd->AT,z,y));
461: PetscCall(VecAXPY(y,-sigma*sign,v));
462: PetscCall(VecNorm(y,NORM_2,norm2));
463: }
464: PetscFunctionReturn(PETSC_SUCCESS);
465: }
467: /*@
468: SVDComputeError - Computes the error (based on the residual norm) associated
469: with the i-th singular triplet.
471: Collective
473: Input Parameters:
474: + svd - the singular value solver context
475: . i - the solution index
476: - type - the type of error to compute
478: Output Parameter:
479: . error - the error
481: Notes:
482: The error can be computed in various ways, all of them based on the residual
483: norm obtained as sqrt(n1^2+n2^2) with n1 = ||A*v-sigma*u||_2 and
484: n2 = ||A^T*u-sigma*v||_2, where sigma is the singular value, u is the left
485: singular vector and v is the right singular vector.
487: In the case of the GSVD, the two components of the residual norm are
488: n1 = ||s^2*A'*u-c*B'*B*x||_2 and n2 = ||c^2*B'*v-s*A'*A*x||_2, where [u;v]
489: are the left singular vectors and x is the right singular vector, with
490: sigma=c/s.
492: Level: beginner
494: .seealso: `SVDErrorType`, `SVDSolve()`
495: @*/
496: PetscErrorCode SVDComputeError(SVD svd,PetscInt i,SVDErrorType type,PetscReal *error)
497: {
498: PetscReal sigma,norm1,norm2,c,s;
499: Vec u=NULL,v=NULL,x=NULL,y=NULL,z=NULL;
500: PetscReal vecnorm=1.0;
502: PetscFunctionBegin;
506: PetscAssertPointer(error,4);
507: SVDCheckSolved(svd,1);
509: /* allocate work vectors */
510: switch (svd->problem_type) {
511: case SVD_STANDARD:
512: PetscCall(SVDSetWorkVecs(svd,2,2));
513: u = svd->workl[0];
514: v = svd->workr[0];
515: x = svd->workl[1];
516: y = svd->workr[1];
517: break;
518: case SVD_GENERALIZED:
519: PetscCall(SVDSetWorkVecs(svd,1,3));
520: u = svd->workl[0];
521: v = svd->workr[0];
522: x = svd->workr[1];
523: y = svd->workr[2];
524: break;
525: case SVD_HYPERBOLIC:
526: PetscCall(SVDSetWorkVecs(svd,3,2));
527: u = svd->workl[0];
528: v = svd->workr[0];
529: x = svd->workl[1];
530: y = svd->workr[1];
531: z = svd->workl[2];
532: break;
533: }
535: /* compute residual norm */
536: PetscCall(SVDGetSingularTriplet(svd,i,&sigma,u,v));
537: switch (svd->problem_type) {
538: case SVD_STANDARD:
539: PetscCall(SVDComputeResidualNorms_Standard(svd,sigma,u,v,x,y,&norm1,&norm2));
540: break;
541: case SVD_GENERALIZED:
542: PetscCall(SVDComputeResidualNorms_Generalized(svd,sigma,u,v,x,y,&norm1,&norm2));
543: break;
544: case SVD_HYPERBOLIC:
545: PetscCall(SVDComputeResidualNorms_Hyperbolic(svd,sigma,svd->sign[svd->perm[i]],u,v,x,y,z,&norm1,&norm2));
546: break;
547: }
548: *error = SlepcAbs(norm1,norm2);
550: /* compute 2-norm of eigenvector of the cyclic form */
551: if (type!=SVD_ERROR_ABSOLUTE) {
552: switch (svd->problem_type) {
553: case SVD_STANDARD:
554: vecnorm = PETSC_SQRT2;
555: break;
556: case SVD_GENERALIZED:
557: PetscCall(VecNorm(v,NORM_2,&vecnorm));
558: vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
559: break;
560: case SVD_HYPERBOLIC:
561: PetscCall(VecNorm(u,NORM_2,&vecnorm));
562: vecnorm = PetscSqrtReal(1.0+vecnorm*vecnorm);
563: break;
564: }
565: }
567: /* compute error */
568: switch (type) {
569: case SVD_ERROR_ABSOLUTE:
570: break;
571: case SVD_ERROR_RELATIVE:
572: if (svd->isgeneralized) {
573: s = 1.0/PetscSqrtReal(1.0+sigma*sigma);
574: c = sigma*s;
575: norm1 /= c*vecnorm;
576: norm2 /= s*vecnorm;
577: *error = PetscMax(norm1,norm2);
578: } else *error /= sigma*vecnorm;
579: break;
580: case SVD_ERROR_NORM:
581: if (!svd->nrma) PetscCall(MatNorm(svd->OP,NORM_INFINITY,&svd->nrma));
582: if (svd->isgeneralized && !svd->nrmb) PetscCall(MatNorm(svd->OPb,NORM_INFINITY,&svd->nrmb));
583: *error /= PetscMax(svd->nrma,svd->nrmb)*vecnorm;
584: break;
585: default:
586: SETERRQ(PetscObjectComm((PetscObject)svd),PETSC_ERR_ARG_OUTOFRANGE,"Invalid error type");
587: }
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }