/*
|
/*
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
SLEPc - Scalable Library for Eigenvalue Problem Computations
|
Copyright (c) 2002-2007, Universidad Politecnica de Valencia, Spain
|
Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain
|
|
|
This file is part of SLEPc. See the README file for conditions of use
|
This file is part of SLEPc.
|
and additional information.
|
|
|
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"
|
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"
|
"This example works for both real and complex numbers.\n\n"
|
"The command line options are:\n"
|
"The command line options are:\n"
|
" -file <filename>, where <filename> = matrix file in PETSc binary form.\n\n";
|
" -file <filename>, where <filename> = matrix file in PETSc binary form.\n\n";
|
|
|
#include "slepcsvd.h"
|
#include "slepcsvd.h"
|
|
|
#undef __FUNCT__
|
#undef __FUNCT__
|
#define __FUNCT__ "main"
|
#define __FUNCT__ "main"
|
int main( int argc, char **argv )
|
int main( int argc, char **argv )
|
{
|
{
|
Mat A; /* operator matrix */
|
Mat A; /* operator matrix */
|
SVD svd; /* singular value problem solver context */
|
SVD svd; /* singular value problem solver context */
|
const SVDType type;
|
const SVDType type;
|
PetscReal error, tol, sigma;
|
PetscReal error, tol, sigma;
|
PetscErrorCode ierr;
|
PetscErrorCode ierr;
|
PetscInt nsv, maxit, i, its, nconv;
|
PetscInt nsv, maxit, i, its, nconv;
|
char filename[256];
|
char filename[256];
|
PetscViewer viewer;
|
PetscViewer viewer;
|
PetscTruth flg;
|
PetscTruth flg;
|
|
|
|
|
SlepcInitialize(&argc,&argv,(char*)0,help);
|
SlepcInitialize(&argc,&argv,(char*)0,help);
|
|
|
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
Load the operator matrix that defines the singular value problem
|
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 = PetscPrintf(PETSC_COMM_WORLD,"\nSingular value problem stored in file.\n\n");CHKERRQ(ierr);
|
ierr = PetscOptionsGetString(PETSC_NULL,"-file",filename,256,&flg);CHKERRQ(ierr);
|
ierr = PetscOptionsGetString(PETSC_NULL,"-file",filename,256,&flg);CHKERRQ(ierr);
|
if (!flg) {
|
if (!flg) {
|
SETERRQ(1,"Must indicate a file name with the -file option.");
|
SETERRQ(1,"Must indicate a file name with the -file option.");
|
}
|
}
|
|
|
#if defined(PETSC_USE_COMPLEX)
|
#if defined(PETSC_USE_COMPLEX)
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");CHKERRQ(ierr);
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");CHKERRQ(ierr);
|
#else
|
#else
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");CHKERRQ(ierr);
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");CHKERRQ(ierr);
|
#endif
|
#endif
|
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
|
ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
|
ierr = MatLoad(viewer,MATAIJ,&A);CHKERRQ(ierr);
|
ierr = MatLoad(viewer,MATAIJ,&A);CHKERRQ(ierr);
|
ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
|
ierr = PetscViewerDestroy(viewer);CHKERRQ(ierr);
|
|
|
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
Create the singular value solver and set various options
|
Create the singular value solver and set various options
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
/*
|
/*
|
Create singular value solver context
|
Create singular value solver context
|
*/
|
*/
|
ierr = SVDCreate(PETSC_COMM_WORLD,&svd);CHKERRQ(ierr);
|
ierr = SVDCreate(PETSC_COMM_WORLD,&svd);CHKERRQ(ierr);
|
|
|
/*
|
/*
|
Set operator
|
Set operator
|
*/
|
*/
|
ierr = SVDSetOperator(svd,A);CHKERRQ(ierr);
|
ierr = SVDSetOperator(svd,A);CHKERRQ(ierr);
|
|
|
/*
|
/*
|
Set solver parameters at runtime
|
Set solver parameters at runtime
|
*/
|
*/
|
ierr = SVDSetFromOptions(svd);CHKERRQ(ierr);
|
ierr = SVDSetFromOptions(svd);CHKERRQ(ierr);
|
|
|
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
Solve the singular value system
|
Solve the singular value system
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
ierr = SVDSolve(svd);CHKERRQ(ierr);
|
ierr = SVDSolve(svd);CHKERRQ(ierr);
|
ierr = SVDGetIterationNumber(svd, &its);CHKERRQ(ierr);
|
ierr = SVDGetIterationNumber(svd, &its);CHKERRQ(ierr);
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",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
|
Optional: Get some information from the solver and display it
|
*/
|
*/
|
ierr = SVDGetType(svd,&type);CHKERRQ(ierr);
|
ierr = SVDGetType(svd,&type);CHKERRQ(ierr);
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",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 = 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 = PetscPrintf(PETSC_COMM_WORLD," Number of requested singular values: %d\n",nsv);CHKERRQ(ierr);
|
ierr = SVDGetTolerances(svd,&tol,&maxit);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);
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);CHKERRQ(ierr);
|
|
|
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
|
Display solution and clean up
|
Display solution and clean up
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
|
|
|
/*
|
/*
|
Get number of converged singular triplets
|
Get number of converged singular triplets
|
*/
|
*/
|
ierr = SVDGetConverged(svd,&nconv);CHKERRQ(ierr);
|
ierr = SVDGetConverged(svd,&nconv);CHKERRQ(ierr);
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate singular triplets: %d\n\n",nconv);CHKERRQ(ierr);
|
ierr = PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate singular triplets: %d\n\n",nconv);CHKERRQ(ierr);
|
|
|
if (nconv>0) {
|
if (nconv>0) {
|
/*
|
/*
|
Display singular values and relative errors
|
Display singular values and relative errors
|
*/
|
*/
|
ierr = PetscPrintf(PETSC_COMM_WORLD,
|
ierr = PetscPrintf(PETSC_COMM_WORLD,
|
" sigma residual norm\n"
|
" sigma residual norm\n"
|
" --------------------- ------------------\n" );CHKERRQ(ierr);
|
" --------------------- ------------------\n" );CHKERRQ(ierr);
|
for( i=0; i<nconv; i++ ) {
|
for( i=0; i<nconv; i++ ) {
|
/*
|
/*
|
Get converged singular triplets: i-th singular value is stored in sigma
|
Get converged singular triplets: i-th singular value is stored in sigma
|
*/
|
*/
|
ierr = SVDGetSingularTriplet(svd,i,&sigma,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
ierr = SVDGetSingularTriplet(svd,i,&sigma,PETSC_NULL,PETSC_NULL);CHKERRQ(ierr);
|
|
|
/*
|
/*
|
Compute the error associated to each singular triplet
|
Compute the error associated to each singular triplet
|
*/
|
*/
|
ierr = SVDComputeRelativeError(svd,i,&error);CHKERRQ(ierr);
|
ierr = SVDComputeRelativeError(svd,i,&error);CHKERRQ(ierr);
|
ierr = PetscPrintf(PETSC_COMM_WORLD," % 6f ",sigma); 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," % 12g\n",error);CHKERRQ(ierr);
|
}
|
}
|
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n" );CHKERRQ(ierr);
|
ierr = PetscPrintf(PETSC_COMM_WORLD,"\n" );CHKERRQ(ierr);
|
}
|
}
|
|
|
/*
|
/*
|
Free work space
|
Free work space
|
*/
|
*/
|
ierr = SVDDestroy(svd);CHKERRQ(ierr);
|
ierr = SVDDestroy(svd);CHKERRQ(ierr);
|
ierr = MatDestroy(A);CHKERRQ(ierr);
|
ierr = MatDestroy(A);CHKERRQ(ierr);
|
ierr = SlepcFinalize();CHKERRQ(ierr);
|
ierr = SlepcFinalize();CHKERRQ(ierr);
|
return 0;
|
return 0;
|
}
|
}
|
|
|
|
|