Subversion Repositories slepc-dev

Compare Revisions

Ignore whitespace Rev 1256 → Rev 1257

/trunk/src/svd/interface/svdbasic.c
92,6 → 92,7
PetscErrorCode ierr;
const char *type;
PetscTruth isascii;
const char *mode_list[3] = { "default" , "explicit", "user" };
 
PetscFunctionBegin;
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
104,10 → 105,11
ierr = PetscViewerASCIIPrintf(viewer,"SVD Object:\n");CHKERRQ(ierr);
ierr = SVDGetType(svd,&type);CHKERRQ(ierr);
if (type) {
ierr = PetscViewerASCIIPrintf(viewer," method: %s",type);CHKERRQ(ierr);
ierr = PetscViewerASCIIPrintf(viewer," method: %s\n",type);CHKERRQ(ierr);
} else {
ierr = PetscViewerASCIIPrintf(viewer," method: not yet set\n");CHKERRQ(ierr);
}
ierr = PetscViewerASCIIPrintf(viewer," transpose mode: %s\n",mode_list[svd->transmode]);CHKERRQ(ierr);
if (svd->ops->view) {
ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
ierr = (*svd->ops->view)(svd,viewer);CHKERRQ(ierr);
153,12 → 155,9
{
PetscErrorCode ierr;
SVD svd;
PetscMPIInt size;
 
 
PetscFunctionBegin;
PetscValidPointer(outsvd,2);
ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
 
PetscHeaderCreate(svd,_p_SVD,struct _SVDOps,SVD_COOKIE,-1,"SVD",comm,SVDDestroy,SVDView);
PetscLogObjectCreate(svd);
170,7 → 169,7
svd->type_name = PETSC_NULL;
svd->A = PETSC_NULL;
svd->AT = PETSC_NULL;
svd->transmode = size == 1 ? SVD_TRANSPOSE_DEFAULT : SVD_TRANSPOSE_EXPLICIT;
svd->transmode = SVD_TRANSPOSE_EXPLICIT;
svd->sigma = PETSC_NULL;
svd->U = PETSC_NULL;
svd->V = PETSC_NULL;
/trunk/src/svd/interface/svdsolve.c
122,37 → 122,38
}
 
#undef __FUNCT__
#define __FUNCT__ "SVDComputeResidualNorm"
#define __FUNCT__ "SVDComputeResidualNorms"
/*@
SVDComputeResidualNorm - Computes the norm of the residual vector associated with
SVDComputeResidualNorm - Computes the norms of the residuals vectors associated with
the i-th computed singular triplet.
 
Collective on EPS
 
Input Parameter:
. svd - the eigensolver context
. i - the solution index
+ svd - the eigensolver context
- i - the solution index
 
Output Parameter:
. norm - the residual norm, computed as ||A*v-sigma*u||_2 where sigma is the
singular value, u and v are the left and right singular vectors.
+ norm1 - the norm ||A*v-sigma*u||_2 where sigma is the
singular value, u and v are the left and right singular vectors.
- norm2 - the norm ||A*u-sigma*v||_2 with the same sigma, u and v
 
Notes:
The index i should be a value between 0 and nconv (see SVDGetConverged()).
Both output parameters can be PETSC_NULL on input if not needed.
 
Level: beginner
 
.seealso: SVDSolve(), SVDGetConverged()
@*/
PetscErrorCode SVDComputeResidualNorm(SVD svd, int i, PetscReal *norm)
PetscErrorCode SVDComputeResidualNorms(SVD svd, int i, PetscReal *norm1, PetscReal *norm2)
{
PetscErrorCode ierr;
Vec u,v,x;
Vec u,v,x = PETSC_NULL,y = PETSC_NULL;
PetscReal sigma;
 
PetscFunctionBegin;
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
PetscValidPointer(sigma,3);
if (svd->nconv < 0) {
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "SVDSolve must be called first");
}
161,14 → 162,27
}
ierr = MatGetVecs(svd->A,&v,&u);CHKERRQ(ierr);
ierr = VecDuplicate(u,&x);CHKERRQ(ierr);
ierr = SVDGetSingularTriplet(svd,i,&sigma,u,v);
ierr = MatMult(svd->A,v,x);CHKERRQ(ierr);
ierr = VecAXPY(x,-sigma,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,norm);CHKERRQ(ierr);
if (norm1) {
ierr = VecDuplicate(u,&x);CHKERRQ(ierr);
ierr = MatMult(svd->A,v,x);CHKERRQ(ierr);
ierr = VecAXPY(x,-sigma,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,norm1);CHKERRQ(ierr);
}
if (norm2) {
ierr = VecDuplicate(v,&y);CHKERRQ(ierr);
if (svd->AT) {
ierr = MatMult(svd->AT,u,y);CHKERRQ(ierr);
} else {
ierr = MatMultTranspose(svd->A,u,y);CHKERRQ(ierr);
}
ierr = VecAXPY(y,-sigma,v);CHKERRQ(ierr);
ierr = VecNorm(y,NORM_2,norm2);CHKERRQ(ierr);
}
 
ierr = VecDestroy(v);CHKERRQ(ierr);
ierr = VecDestroy(u);CHKERRQ(ierr);
ierr = VecDestroy(x);CHKERRQ(ierr);
if (x) { ierr = VecDestroy(x);CHKERRQ(ierr); }
if (y) { ierr = VecDestroy(y);CHKERRQ(ierr); }
PetscFunctionReturn(0);
}