Actual source code: ex14.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2009, Universidad Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:       
  8:    SLEPc is free software: you can redistribute it and/or modify it under  the
  9:    terms of version 3 of the GNU Lesser General Public License as published by
 10:    the Free Software Foundation.

 12:    SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY 
 13:    WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS 
 14:    FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for 
 15:    more details.

 17:    You  should have received a copy of the GNU Lesser General  Public  License
 18:    along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
 19:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 20: */

 22: static char help[] = "Solves a singular value problem with the matrix loaded from a file.\n"
 23:   "This example works for both real and complex numbers.\n\n"
 24:   "The command line options are:\n"
 25:   "  -file <filename>, where <filename> = matrix file in PETSc binary form.\n\n";

 27:  #include slepcsvd.h

 31: int main( int argc, char **argv )
 32: {
 33:   Mat                  A;                  /* operator matrix */
 34:   SVD                  svd;                  /* singular value problem solver context */
 35:   const SVDType  type;
 36:   PetscReal            error, tol, sigma;
 38:   PetscInt             nsv, maxit, i, its, nconv;
 39:   char                 filename[256];
 40:   PetscViewer          viewer;
 41:   PetscTruth           flg;


 44:   SlepcInitialize(&argc,&argv,(char*)0,help);

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 47:         Load the operator matrix that defines the singular value problem
 48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 50:   PetscPrintf(PETSC_COMM_WORLD,"\nSingular value problem stored in file.\n\n");
 51:   PetscOptionsGetString(PETSC_NULL,"-file",filename,256,&flg);
 52:   if (!flg) {
 53:     SETERRQ(1,"Must indicate a file name with the -file option.");
 54:   }

 56: #if defined(PETSC_USE_COMPLEX)
 57:   PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");
 58: #else
 59:   PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");
 60: #endif
 61:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 62:   MatLoad(viewer,MATAIJ,&A);
 63:   PetscViewerDestroy(viewer);

 65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 66:                 Create the singular value solver and set various options
 67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 69:   /* 
 70:      Create singular value solver context
 71:   */
 72:   SVDCreate(PETSC_COMM_WORLD,&svd);

 74:   /* 
 75:      Set operator
 76:   */
 77:   SVDSetOperator(svd,A);

 79:   /*
 80:      Set solver parameters at runtime
 81:   */
 82:   SVDSetFromOptions(svd);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 85:                       Solve the singular value system
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 88:   SVDSolve(svd);
 89:   SVDGetIterationNumber(svd, &its);
 90:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

 92:   /*
 93:      Optional: Get some information from the solver and display it
 94:   */
 95:   SVDGetType(svd,&type);
 96:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
 97:   SVDGetDimensions(svd,&nsv,PETSC_NULL,PETSC_NULL);
 98:   PetscPrintf(PETSC_COMM_WORLD," Number of requested singular values: %d\n",nsv);
 99:   SVDGetTolerances(svd,&tol,&maxit);
100:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
103:                     Display solution and clean up
104:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

106:   /* 
107:      Get number of converged singular triplets
108:   */
109:   SVDGetConverged(svd,&nconv);
110:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate singular triplets: %d\n\n",nconv);

112:   if (nconv>0) {
113:     /*
114:        Display singular values and relative errors
115:     */
116:     PetscPrintf(PETSC_COMM_WORLD,
117:          "          sigma           residual norm\n"
118:          "  --------------------- ------------------\n" );
119:     for( i=0; i<nconv; i++ ) {
120:       /* 
121:          Get converged singular triplets: i-th singular value is stored in sigma
122:       */
123:       SVDGetSingularTriplet(svd,i,&sigma,PETSC_NULL,PETSC_NULL);

125:       /*
126:          Compute the error associated to each singular triplet 
127:       */
128:       SVDComputeRelativeError(svd,i,&error);
129:       PetscPrintf(PETSC_COMM_WORLD,"       % 6f      ",sigma);
130:       PetscPrintf(PETSC_COMM_WORLD," % 12g\n",error);
131:     }
132:     PetscPrintf(PETSC_COMM_WORLD,"\n" );
133:   }
134: 
135:   /* 
136:      Free work space
137:   */
138:   SVDDestroy(svd);
139:   MatDestroy(A);
140:   SlepcFinalize();
141:   return 0;
142: }