Subversion Repositories slepc-dev

Compare Revisions

Regard whitespace Rev 1253 → Rev 1252

/trunk/include/slepcsvd.h
4,7 → 4,6
#if !defined(__SLEPCSVD_H)
#define __SLEPCSVD_H
#include "slepc.h"
#include "slepceps.h"
PETSC_EXTERN_CXX_BEGIN
 
extern PetscCookie SVD_COOKIE;
38,12 → 37,4
EXTERN PetscErrorCode SVDDestroy(SVD);
EXTERN PetscErrorCode SVDInitializePackage(char*);
 
typedef enum { SVDEIGENSOLVER_DIRECT, SVDEIGENSOLVER_TRANSPOSE,
SVDEIGENSOLVER_CYCLIC } SVDEigensolverMode;
 
EXTERN PetscErrorCode SVDEigensolverSetMode(SVD,SVDEigensolverMode);
EXTERN PetscErrorCode SVDEigensolverGetMode(SVD,SVDEigensolverMode*);
EXTERN PetscErrorCode SVDEigensolverSetEPS(SVD,EPS);
EXTERN PetscErrorCode SVDEigensolverGetEPS(SVD,EPS*);
 
#endif
/trunk/src/svd/interface/svdsetup.c
88,7 → 88,6
 
PetscFunctionBegin;
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
svd->setupcalled = 0;
ierr = PetscOptionsBegin(svd->comm,svd->prefix,"Singular Value Solver (SVD) Options","SVD");CHKERRQ(ierr);
 
ierr = PetscOptionsList("-svd_type","Singular Value Solver method","SVDSetType",SVDList,(char*)(svd->type_name?svd->type_name:SVDEIGENSOLVER),type,256,&flg);CHKERRQ(ierr);
/trunk/src/svd/impls/eigensolver/eigensolver.c
1,6 → 1,6
/*
 
SLEPc singular value solver: "eigensolver"
SLEPc singular value solver: "eigen"
 
Method: Uses an Hermitian eigensolver for A^T*A, A*A^T or H(A)
 
11,7 → 11,6
#include "slepceps.h"
 
typedef struct {
SVDEigensolverMode mode;
EPS eps;
Mat mat;
Vec w;
28,20 → 27,8
PetscFunctionBegin;
ierr = MatShellGetContext(B,(void**)&svd);CHKERRQ(ierr);
eigen = (SVD_EIGENSOLVER *)svd->data;
switch (eigen->mode) {
case SVDEIGENSOLVER_DIRECT:
ierr = MatMult(svd->A,x,eigen->w);CHKERRQ(ierr);
ierr = MatMultTranspose(svd->A,eigen->w,y);CHKERRQ(ierr);
break;
case SVDEIGENSOLVER_TRANSPOSE:
ierr = MatMultTranspose(svd->A,x,eigen->w);CHKERRQ(ierr);
ierr = MatMult(svd->A,eigen->w,y);CHKERRQ(ierr);
break;
case SVDEIGENSOLVER_CYCLIC:
SETERRQ(1,"Not implemented :-)");
default:
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid SVD type");
}
PetscFunctionReturn(0);
}
 
53,34 → 40,24
PetscErrorCode ierr;
SVD_EIGENSOLVER *eigen = (SVD_EIGENSOLVER *)svd->data;
PetscInt m,n,M,N;
Vec x;
 
PetscFunctionBegin;
if (eigen->w) { ierr = VecDestroy(eigen->w);CHKERRQ(ierr); }
if (eigen->mat) { ierr = MatDestroy(eigen->mat);CHKERRQ(ierr); }
ierr = MatGetVecs(svd->A,&x,&eigen->w);CHKERRQ(ierr);
ierr = VecGetSize(x,&N);CHKERRQ(ierr);
ierr = VecGetLocalSize(x,&n);CHKERRQ(ierr);
ierr = VecDestroy(x);CHKERRQ(ierr);
 
ierr = MatGetSize(svd->A,&M,&N);CHKERRQ(ierr);
ierr = MatGetLocalSize(svd->A,&m,&n);CHKERRQ(ierr);
switch (eigen->mode) {
case SVDEIGENSOLVER_DIRECT:
ierr = MatGetVecs(svd->A,PETSC_NULL,&eigen->w);CHKERRQ(ierr);
if (eigen->mat) { ierr = MatDestroy(eigen->mat);CHKERRQ(ierr); }
ierr = MatCreateShell(svd->comm,n,n,N,N,svd,&eigen->mat);CHKERRQ(ierr);
break;
case SVDEIGENSOLVER_TRANSPOSE:
ierr = MatGetVecs(svd->A,&eigen->w,PETSC_NULL);CHKERRQ(ierr);
ierr = MatCreateShell(svd->comm,m,m,M,M,svd,&eigen->mat);CHKERRQ(ierr);
break;
case SVDEIGENSOLVER_CYCLIC:
SETERRQ(1,"Not implemented :-)");
default:
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid SVD type");
}
ierr = MatShellSetOperation(eigen->mat,MATOP_MULT,(void(*)(void))ShellMatMult_EIGENSOLVER);CHKERRQ(ierr);
 
ierr = EPSSetOperators(eigen->eps,eigen->mat,PETSC_NULL);CHKERRQ(ierr);
/* EPSSetProblemType deberia estar en SVDSetFromOptions
PENDIENTE: ARREGLAR EL PROBLEMA EN EPS */
ierr = EPSSetProblemType(eigen->eps,EPS_HEP);CHKERRQ(ierr);
ierr = EPSSetProblemType(eigen->eps,EPS_NHEP);CHKERRQ(ierr);
ierr = EPSSetUp(eigen->eps);CHKERRQ(ierr);
 
PetscFunctionReturn(0);
118,25 → 95,11
ierr = PetscMalloc(svd->nconv*sizeof(Vec),&svd->V);CHKERRQ(ierr);
for (i=0;i<svd->nconv;i++) {
ierr = MatGetVecs(svd->A,svd->V+i,svd->U+i);CHKERRQ(ierr);
switch (eigen->mode) {
case SVDEIGENSOLVER_DIRECT:
ierr = EPSGetEigenpair(eigen->eps,i,&sigma,PETSC_NULL,svd->V[i],PETSC_NULL);CHKERRQ(ierr);
svd->sigma[i] = sqrt(PetscRealPart(sigma));
ierr = MatMult(svd->A,svd->V[i],svd->U[i]);CHKERRQ(ierr);
ierr = VecScale(svd->U[i],1.0/svd->sigma[i]);CHKERRQ(ierr);
break;
case SVDEIGENSOLVER_TRANSPOSE:
ierr = EPSGetEigenpair(eigen->eps,i,&sigma,PETSC_NULL,svd->U[i],PETSC_NULL);CHKERRQ(ierr);
svd->sigma[i] = sqrt(PetscRealPart(sigma));
ierr = MatMultTranspose(svd->A,svd->U[i],svd->V[i]);CHKERRQ(ierr);
ierr = VecScale(svd->V[i],1.0/svd->sigma[i]);CHKERRQ(ierr);
break;
case SVDEIGENSOLVER_CYCLIC:
SETERRQ(1,"Not implemented :-)");
default:
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid SVD type");
}
}
PetscFunctionReturn(0);
}
 
146,14 → 109,9
{
PetscErrorCode ierr;
SVD_EIGENSOLVER *eigen = (SVD_EIGENSOLVER *)svd->data;
PetscTruth flg;
const char *mode_list[3] = { "direct" , "transpose", "cyclic" };
PetscInt mode;
 
PetscFunctionBegin;
ierr = PetscOptionsHead("EIGEN options");CHKERRQ(ierr);
ierr = PetscOptionsEList("-svd_eigensolver_mode","Eigensolver SVD mode","SVDEigenSolverSetMode",mode_list,3,mode_list[eigen->mode],&mode,&flg);CHKERRQ(ierr);
if (flg) { eigen->mode = (SVDEigensolverMode)mode; }
ierr = EPSSetFromOptions(eigen->eps);
ierr = PetscOptionsTail();CHKERRQ(ierr);
PetscFunctionReturn(0);
161,112 → 119,9
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "SVDEigensolverSetMode_EIGENSOLVER"
PetscErrorCode SVDEigensolverSetMode_EIGENSOLVER(SVD svd,SVDEigensolverMode mode)
#define __FUNCT__ "SVDEigenSetEPS_EIGENSOLVER"
PetscErrorCode SVDEigenSetEPS_EIGENSOLVER(SVD svd,EPS eps)
{
SVD_EIGENSOLVER *eigen = (SVD_EIGENSOLVER *)svd->data;
 
PetscFunctionBegin;
switch (eigen->mode) {
case SVDEIGENSOLVER_DIRECT:
case SVDEIGENSOLVER_TRANSPOSE:
case SVDEIGENSOLVER_CYCLIC:
eigen->mode = mode;
svd->setupcalled = 0;
break;
default:
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Invalid SVD type");
}
PetscFunctionReturn(0);
}
EXTERN_C_BEGIN
 
#undef __FUNCT__
#define __FUNCT__ "SVDEigensolverSetMode"
/*@
SVDEigensolverSetMode - Sets the transformation used in the eigensolver.
 
Collective on SVD
 
Input Parameters:
+ svd - singular value solver context obtained from SVDCreate()
- mode - the mode flag, one of SVDEIGENSOLVER_DIRECT,
SVDEIGENSOLVER_TRANSPOSE or SVDEIGENSOLVER_CYCLIC
 
Options Database Key:
. -svd_eigensolver_mode <mode> - Indicates the mode flag, where <mode>
is one of 'direct', 'transpose' or 'cyclic' (see explanation below).
 
Notes:
This parameter selects the eigensystem used to compute the SVD:
A^T*A (SVDEIGENSOLVER_DIRECT), A*A^T (SVDEIGENSOLVER_TRANSPOSE)
or H(A) = [ 0 A ; A^T 0 ] (SVDEIGENSOLVER_CYCLIC).
 
Level: beginner
 
.seealso: SVDEigensolverGetMode()
@*/
PetscErrorCode SVDEigensolverSetMode(SVD svd,SVDEigensolverMode mode)
{
PetscErrorCode ierr, (*f)(SVD,SVDEigensolverMode);
 
PetscFunctionBegin;
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDEigensolverSetMode_C",(void (**)())&f);CHKERRQ(ierr);
if (f) {
ierr = (*f)(svd,mode);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "SVDEigensolverGetMode_EIGENSOLVER"
PetscErrorCode SVDEigensolverGetMode_EIGENSOLVER(SVD svd,SVDEigensolverMode *mode)
{
SVD_EIGENSOLVER *eigen = (SVD_EIGENSOLVER *)svd->data;
 
PetscFunctionBegin;
PetscValidPointer(mode,2);
*mode = eigen->mode;
PetscFunctionReturn(0);
}
EXTERN_C_BEGIN
 
#undef __FUNCT__
#define __FUNCT__ "SVDEigensolverGetMode"
/*@C
SVDEigensolverGetMode - Gets the transformation used by the eigensolver.
 
Not collective
 
Input Parameters:
+ svd - singular value solver context obtained from SVDCreate()
Output Parameters:
- mode - the mode flag
 
Level: beginner
 
.seealso: SVDEigensolverSetMode()
@*/
PetscErrorCode SVDEigensolverGetMode(SVD svd,SVDEigensolverMode *mode)
{
PetscErrorCode ierr, (*f)(SVD,SVDEigensolverMode*);
 
PetscFunctionBegin;
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDEigensolverGetMode_C",(void (**)())&f);CHKERRQ(ierr);
if (f) {
ierr = (*f)(svd,mode);CHKERRQ(ierr);
}
PetscFunctionReturn(0);
}
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "SVDEigensolverSetEPS_EIGENSOLVER"
PetscErrorCode SVDEigensolverSetEPS_EIGENSOLVER(SVD svd,EPS eps)
{
PetscErrorCode ierr;
SVD_EIGENSOLVER *eigen = (SVD_EIGENSOLVER *)svd->data;
 
282,9 → 137,9
EXTERN_C_END
 
#undef __FUNCT__
#define __FUNCT__ "SVDEigensolverSetEPS"
#define __FUNCT__ "SVDEigenSetEPS"
/*@
SVDEigensolverSetEPS - Associates a eigensolver object to the
SVDEigenSetEPS - Associates a eigensolver object to the
singular value solver.
 
Collective on SVD
295,15 → 150,15
 
Level: advanced
 
.seealso: SVDEigensolverGetEPS()
.seealso: SVDEigenGetEPS()
@*/
PetscErrorCode SVDEigensolverSetEPS(SVD svd,EPS eps)
PetscErrorCode SVDEigenSetEPS(SVD svd,EPS eps)
{
PetscErrorCode ierr, (*f)(SVD,EPS eps);
 
PetscFunctionBegin;
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDEigensolverSetEPS_C",(void (**)())&f);CHKERRQ(ierr);
ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDEigenSetEPS_C",(void (**)())&f);CHKERRQ(ierr);
if (f) {
ierr = (*f)(svd,eps);CHKERRQ(ierr);
}
312,8 → 167,8
 
EXTERN_C_BEGIN
#undef __FUNCT__
#define __FUNCT__ "SVDEigensolverGetEPS_EIGENSOLVER"
PetscErrorCode SVDEigensolverGetEPS_EIGENSOLVER(SVD svd,EPS *eps)
#define __FUNCT__ "SVDEigenGetEPS_EIGENSOLVER"
PetscErrorCode SVDEigenGetEPS_EIGENSOLVER(SVD svd,EPS *eps)
{
SVD_EIGENSOLVER *eigen = (SVD_EIGENSOLVER *)svd->data;
 
325,9 → 180,9
EXTERN_C_END
 
#undef __FUNCT__
#define __FUNCT__ "SVDEigensolverGetEPS"
/*@C
SVDEigensolverGetEPS - Obtain the eigensolver (EPS) object associated
#define __FUNCT__ "SVDEigenGetEPS"
/*@
SVDEigenGetEPS - Obtain the eigensolver (EPS) object associated
to the singular value solver object.
 
Not Collective
338,17 → 193,15
Output Parameter:
. eps - the eigensolver object
 
Level: advanced
 
.seealso: SVDEigensolverSetEPS()
Level: beginner
@*/
PetscErrorCode SVDEigensolverGetEPS(SVD svd,EPS *eps)
PetscErrorCode SVDEigenGetEPS(SVD svd,EPS *eps)
{
PetscErrorCode ierr, (*f)(SVD,EPS *eps);
 
PetscFunctionBegin;
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDEigensolverGetEPS_C",(void (**)())&f);CHKERRQ(ierr);
ierr = PetscObjectQueryFunction((PetscObject)svd,"SVDEigenGetEPS_C",(void (**)())&f);CHKERRQ(ierr);
if (f) {
ierr = (*f)(svd,eps);CHKERRQ(ierr);
}
398,16 → 251,13
svd->ops->setfromoptions = SVDSetFromOptions_EIGENSOLVER;
svd->ops->destroy = SVDDestroy_EIGENSOLVER;
svd->ops->view = SVDView_EIGENSOLVER;
ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDEigensolverSetEPS_C","SVDEigensolverSetEPS_EIGENSOLVER",SVDEigensolverSetEPS_EIGENSOLVER);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDEigensolverGetEPS_C","SVDEigensolverGetEPS_EIGENSOLVER",SVDEigensolverGetEPS_EIGENSOLVER);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDEigensolverSetMode_C","SVDEigensolverSetMode_EIGENSOLVER",SVDEigensolverSetMode_EIGENSOLVER);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDEigensolverGetMode_C","SVDEigensolverGetMode_EIGENSOLVER",SVDEigensolverGetMode_EIGENSOLVER);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDEigenSetEPS_C","SVDEigenSetEPS_EIGENSOLVER",SVDEigenSetEPS_EIGENSOLVER);CHKERRQ(ierr);
ierr = PetscObjectComposeFunctionDynamic((PetscObject)svd,"SVDEigenGetEPS_C","SVDEigenGetEPS_EIGENSOLVER",SVDEigenGetEPS_EIGENSOLVER);CHKERRQ(ierr);
 
ierr = EPSCreate(svd->comm,&eigen->eps);CHKERRQ(ierr);
ierr = EPSSetOptionsPrefix(eigen->eps,svd->prefix);CHKERRQ(ierr);
ierr = EPSAppendOptionsPrefix(eigen->eps,"svd_");CHKERRQ(ierr);
PetscLogObjectParent(svd,eigen->eps);
eigen->mode = SVDEIGENSOLVER_DIRECT;
eigen->mat = PETSC_NULL;
eigen->w = PETSC_NULL;
PetscFunctionReturn(0);
/trunk/src/examples/ex8.c
4,13 → 4,14
"The command line options are:\n"
" -n <n>, where <n> = matrix dimension.\n\n";
 
#include "slepcsvd.h"
#include "slepceps.h"
 
/*
This example computes the singular values of an nxn Grcar matrix,
which is a nonsymmetric Toeplitz matrix:
This example computes the singular values of A by computing the eigenvalues
of A^T*A, where A^T denotes the transpose of A.
 
An nxn Grcar matrix is a nonsymmetric Toeplitz matrix:
 
| 1 1 1 1 |
| -1 1 1 1 1 |
| -1 1 1 1 1 |
21,19 → 22,49
| -1 1 1 |
| -1 1 |
 
Note that working with A^T*A can lead to poor accuracy of the computed
singular values when A is ill-conditioned (which is not the case here).
Another alternative would be to compute the eigenvalues of
 
| 0 A |
| A^T 0 |
 
but this significantly increases the cost of the solution process.
 
*/
 
 
/*
Matrix multiply routine
*/
#undef __FUNCT__
#define __FUNCT__ "MatSVD_Mult"
PetscErrorCode MatSVD_Mult(Mat H,Vec x,Vec y)
{
Mat A;
Vec w;
PetscErrorCode ierr;
 
PetscFunctionBegin;
ierr = MatShellGetContext(H,(void**)&A);CHKERRQ(ierr);
ierr = MatGetVecs(A,PETSC_NULL,&w);CHKERRQ(ierr);
ierr = MatMult(A,x,w);CHKERRQ(ierr);
ierr = MatMultTranspose(A,w,y);CHKERRQ(ierr);
ierr = VecDestroy(w);CHKERRQ(ierr);
PetscFunctionReturn(0);
}
 
#undef __FUNCT__
#define __FUNCT__ "main"
int main( int argc, char **argv )
{
PetscErrorCode ierr;
Mat A; /* Grcar matrix */
SVD svd; /* singular value solver context */
Mat H; /* eigenvalue problem matrix, H=A^T*A */
EPS eps; /* eigenproblem solver context */
PetscInt N=30, Istart, Iend, i, col[5];
PetscInt N=30, n, Istart, Iend, i, col[5];
int nconv1, nconv2;
PetscScalar value[] = { -1, 1, 1, 1, 1 };
PetscScalar kl, ks, value[] = { -1, 1, 1, 1, 1 };
PetscReal sigma_1, sigma_n;
 
SlepcInitialize(&argc,&argv,(char*)0,help);
63,27 → 94,35
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
 
/*
Now create a symmetric shell matrix H=A^T*A
*/
ierr = MatGetLocalSize(A,PETSC_NULL,&n);CHKERRQ(ierr);
ierr = MatCreateShell(PETSC_COMM_WORLD,n,n,PETSC_DETERMINE,PETSC_DETERMINE,(void*)A,&H);CHKERRQ(ierr);
ierr = MatShellSetOperation(H,MATOP_MULT,(void(*)())MatSVD_Mult);CHKERRQ(ierr);
ierr = MatShellSetOperation(H,MATOP_MULT_TRANSPOSE,(void(*)())MatSVD_Mult);CHKERRQ(ierr);
 
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the singular value solver and set the solution method
Create the eigensolver and set the solution method
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 
/*
Create singular value context
Create eigensolver context
*/
ierr = SVDCreate(PETSC_COMM_WORLD,&svd);CHKERRQ(ierr);
ierr = EPSCreate(PETSC_COMM_WORLD,&eps);CHKERRQ(ierr);
 
/*
Set operator
Set operators. In this case, it is a standard symmetric eigenvalue problem
*/
ierr = SVDSetOperator(svd,A);CHKERRQ(ierr);
ierr = EPSSetOperators(eps,H,PETSC_NULL);CHKERRQ(ierr);
ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);
 
/*
Set solver parameters at runtime
*/
ierr = SVDSetFromOptions(svd);CHKERRQ(ierr);
/* PENDIENTE: utilizar funciones SVD en lugar de EPS*/
ierr = SVDEigensolverGetEPS(svd,&eps);CHKERRQ(ierr);
ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
ierr = EPSSetDimensions(eps,1,PETSC_DEFAULT);CHKERRQ(ierr);
ierr = EPSSetTolerances(eps,PETSC_DEFAULT,1000);CHKERRQ(ierr);
 
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the eigensystem
93,40 → 132,43
First request an eigenvalue from one end of the spectrum
*/
ierr = EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);CHKERRQ(ierr);
ierr = SVDSolve(svd);CHKERRQ(ierr);
ierr = EPSSolve(eps);CHKERRQ(ierr);
/*
Get number of converged singular values
Get number of converged eigenpairs
*/
ierr = SVDGetConverged(svd,&nconv1);CHKERRQ(ierr);
ierr = EPSGetConverged(eps,&nconv1);CHKERRQ(ierr);
/*
Get converged singular values: largest singular value is stored in sigma_1.
In this example, we are not interested in the singular vectors
Get converged eigenpairs: largest eigenvalue is stored in kl. In this
example, we are not interested in the eigenvectors
*/
if (nconv1 > 0) {
ierr = SVDGetSingularTriplet(svd,0,&sigma_1,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
ierr = EPSGetEigenpair(eps,0,&kl,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
}
 
/*
Request an eigenvalue from the other end of the spectrum
*/
ierr = EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);CHKERRQ(ierr);
ierr = SVDSolve(svd);CHKERRQ(ierr);
ierr = EPSSolve(eps);CHKERRQ(ierr);
/*
Get number of converged eigenpairs
*/
ierr = SVDGetConverged(svd,&nconv2);CHKERRQ(ierr);
ierr = EPSGetConverged(eps,&nconv2);CHKERRQ(ierr);
/*
Get converged singular values: smallest singular value is stored in sigma_n.
As before, we are not interested in the singular vectors
Get converged eigenpairs: smallest eigenvalue is stored in ks. In this
example, we are not interested in the eigenvectors
*/
if (nconv2 > 0) {
ierr = SVDGetSingularTriplet(svd,0,&sigma_n,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
ierr = EPSGetEigenpair(eps,0,&ks,PETSC_NULL,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
}
 
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Display solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
if (nconv1 > 0 && nconv2 > 0) {
sigma_1 = sqrt(PetscRealPart(kl));
sigma_n = sqrt(PetscRealPart(ks));
 
ierr = PetscPrintf(PETSC_COMM_WORLD," Computed singular values: sigma_1=%6f, sigma_n=%6f\n",sigma_1,sigma_n);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Estimated condition number: sigma_1/sigma_n=%6f\n\n",sigma_1/sigma_n);CHKERRQ(ierr);
} else {
136,8 → 178,9
/*
Free work space
*/
ierr = SVDDestroy(svd);CHKERRQ(ierr);
ierr = EPSDestroy(eps);CHKERRQ(ierr);
ierr = MatDestroy(A);CHKERRQ(ierr);
ierr = MatDestroy(H);CHKERRQ(ierr);
ierr = SlepcFinalize();CHKERRQ(ierr);
return 0;
}