| Line 120... |
Line 120... |
}
|
}
|
PetscFunctionReturn(0);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
#undef __FUNCT__
|
#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.
|
the i-th computed singular triplet.
|
|
|
Collective on EPS
|
Collective on EPS
|
|
|
Input Parameter:
|
Input Parameter:
|
. svd - the eigensolver context
|
+ svd - the eigensolver context
|
. i - the solution index
|
- i - the solution index
|
|
|
Output Parameter:
|
Output Parameter:
|
. norm - the residual norm, computed as ||A*v-sigma*u||_2 where sigma is the
|
+ norm1 - the norm ||A*v-sigma*u||_2 where sigma is the
|
singular value, u and v are the left and right singular vectors.
|
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:
|
Notes:
|
The index i should be a value between 0 and nconv (see SVDGetConverged()).
|
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
|
Level: beginner
|
|
|
.seealso: SVDSolve(), SVDGetConverged()
|
.seealso: SVDSolve(), SVDGetConverged()
|
@*/
|
@*/
|
PetscErrorCode SVDComputeResidualNorm(SVD svd, int i, PetscReal *norm)
|
PetscErrorCode SVDComputeResidualNorms(SVD svd, int i, PetscReal *norm1, PetscReal *norm2)
|
{
|
{
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
Vec u,v,x;
|
Vec u,v,x = PETSC_NULL,y = PETSC_NULL;
|
PetscReal sigma;
|
PetscReal sigma;
|
|
|
PetscFunctionBegin;
|
PetscFunctionBegin;
|
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
|
PetscValidHeaderSpecific(svd,SVD_COOKIE,1);
|
PetscValidPointer(sigma,3);
|
|
if (svd->nconv < 0) {
|
if (svd->nconv < 0) {
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "SVDSolve must be called first");
|
SETERRQ(PETSC_ERR_ARG_WRONGSTATE, "SVDSolve must be called first");
|
}
|
}
|
if (i<0 || i>=svd->nconv) {
|
if (i<0 || i>=svd->nconv) {
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
SETERRQ(PETSC_ERR_ARG_OUTOFRANGE, "Argument 2 out of range");
|
}
|
}
|
|
|
ierr = MatGetVecs(svd->A,&v,&u);CHKERRQ(ierr);
|
ierr = MatGetVecs(svd->A,&v,&u);CHKERRQ(ierr);
|
ierr = VecDuplicate(u,&x);CHKERRQ(ierr);
|
|
ierr = SVDGetSingularTriplet(svd,i,&sigma,u,v);
|
ierr = SVDGetSingularTriplet(svd,i,&sigma,u,v);
|
ierr = MatMult(svd->A,v,x);CHKERRQ(ierr);
|
if (norm1) {
|
ierr = VecAXPY(x,-sigma,u);CHKERRQ(ierr);
|
ierr = VecDuplicate(u,&x);CHKERRQ(ierr);
|
ierr = VecNorm(x,NORM_2,norm);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(v);CHKERRQ(ierr);
|
ierr = VecDestroy(u);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);
|
PetscFunctionReturn(0);
|
}
|
}
|
|
|
|
|