/*
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
SLEPc - Scalable Library for Eigenvalue Problem Computations
Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain
This file is part of SLEPc.
SLEPc is free software: you can redistribute it and/or modify it under the
terms of version 3 of the GNU Lesser General Public License as published by
the Free Software Foundation.
SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
more details.
You should have received a copy of the GNU Lesser General Public License
along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
*/
static char help[] = "Solves a singular value problem with the matrix loaded from a file.\n"
"This example works for both real and complex numbers.\n\n"
"The command line options are:\n"
" -file <filename>, where <filename> = matrix file in PETSc binary form.\n\n";
#include "slepcsvd.h"
#undef __FUNCT__
#define __FUNCT__ "main"
int main( int argc, char **argv )
{
Mat A; /* operator matrix */
SVD svd; /* singular value problem solver context */
const SVDType type;
PetscReal error, tol, sigma;
PetscErrorCode ierr;
PetscInt nsv, maxit, i, its, nconv;
char filename[256];
PetscViewer viewer;
PetscTruth flg;
SlepcInitialize(&argc,&argv,(char*)0,help);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Load the operator matrix that defines the singular value problem
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = PetscPrintf(PETSC_COMM_WORLD,"\nSingular value problem stored in file.\n\n");CHKERRQ(ierr);
ierr = PetscOptionsGetString(PETSC_NULL,"-file",filename,256,&flg);CHKERRQ(ierr);
if (!flg) {
SETERRQ(1,"Must indicate a file name with the -file option.");
}
#if defined(PETSC_USE_COMPLEX)
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");CHKERRQ(ierr);
#else
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");CHKERRQ(ierr);
#endif
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
ierr = MatLoad(viewer,MATAIJ,&A);CHKERRQ(ierr);
ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the singular value solver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create singular value solver context
*/
ierr = SVDCreate(PETSC_COMM_WORLD,&svd);CHKERRQ(ierr);
/*
Set operator
*/
ierr = SVDSetOperator(svd,A);CHKERRQ(ierr);
/*
Set solver parameters at runtime
*/
ierr = SVDSetFromOptions(svd);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve the singular value system
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = SVDSolve(svd);CHKERRQ(ierr);
ierr = SVDGetIterationNumber(svd, &its);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);CHKERRQ(ierr);
/*
Optional: Get some information from the solver and display it
*/
ierr = SVDGetType(svd,&type);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);CHKERRQ(ierr);
ierr = SVDGetDimensions(svd,&nsv,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of requested singular values: %d\n",nsv);CHKERRQ(ierr);
ierr = SVDGetTolerances(svd,&tol,&maxit);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Display solution and clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Get number of converged singular triplets
*/
ierr = SVDGetConverged(svd,&nconv);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate singular triplets: %d\n\n",nconv);CHKERRQ(ierr);
if (nconv>0) {
/*
Display singular values and relative errors
*/
ierr = PetscPrintf(PETSC_COMM_WORLD,
" sigma residual norm\n"
" --------------------- ------------------\n" );CHKERRQ(ierr);
for( i=0; i<nconv; i++ ) {
/*
Get converged singular triplets: i-th singular value is stored in sigma
*/
ierr = SVDGetSingularTriplet(svd,i,&sigma,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
/*
Compute the error associated to each singular triplet
*/
ierr = SVDComputeRelativeError(svd,i,&error);CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," % 6f ",sigma); CHKERRQ(ierr);
ierr = PetscPrintf(PETSC_COMM_WORLD," % 12g\n",error);CHKERRQ(ierr);
}
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n" );CHKERRQ(ierr);
}
/*
Free work space
*/
ierr = SVDDestroy(svd);CHKERRQ(ierr);
ierr = MatDestroy(A);CHKERRQ(ierr);
ierr = SlepcFinalize();CHKERRQ(ierr);
return 0;
}