Actual source code: ex15.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[] = "Singular value decomposition of the Lauchli matrix.\n"
 23:   "The command line options are:\n"
 24:   "  -n <n>, where <n> = matrix dimension.\n"
 25:   "  -mu <mu>, where <mu> = subdiagonal value.\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, mu=PETSC_SQRT_MACHINE_EPSILON;
 38:   PetscInt       n=100, i, j, Istart, Iend, nsv, maxit, its, nconv;

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

 42:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 43:   PetscOptionsGetReal(PETSC_NULL,"-mu",&mu,PETSC_NULL);
 44:   PetscPrintf(PETSC_COMM_WORLD,"\nLauchli singular value decomposition, (%d x %d) mu=%g\n\n",n+1,n,mu);

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 47:                           Build the Lauchli matrix
 48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 50:   MatCreate(PETSC_COMM_WORLD,&A);
 51:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n+1,n);
 52:   MatSetFromOptions(A);

 54:   MatGetOwnershipRange(A,&Istart,&Iend);
 55:   for (i=Istart;i<Iend;i++) {
 56:     if (i == 0) {
 57:       for (j=0;j<n;j++) {
 58:         MatSetValue(A,0,j,1.0,INSERT_VALUES);
 59:       }
 60:     } else {
 61:       MatSetValue(A,i,i-1,mu,INSERT_VALUES);
 62:     }
 63:   }
 64: 
 65:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 66:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 68:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 69:           Create the singular value solver and set various options
 70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 72:   /* 
 73:      Create singular value solver context
 74:   */
 75:   SVDCreate(PETSC_COMM_WORLD,&svd);

 77:   /* 
 78:      Set operator
 79:   */
 80:   SVDSetOperator(svd,A);
 81: 
 82:   /*
 83:      Use thick-restart Lanczos as default solver
 84:   */
 85:   SVDSetType(svd,SVDTRLANCZOS);

 87:   /*
 88:      Set solver parameters at runtime
 89:   */
 90:   SVDSetFromOptions(svd);

 92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 93:                       Solve the singular value system
 94:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 96:   SVDSolve(svd);
 97:   SVDGetIterationNumber(svd, &its);
 98:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

100:   /*
101:      Optional: Get some information from the solver and display it
102:   */
103:   SVDGetType(svd,&type);
104:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
105:   SVDGetDimensions(svd,&nsv,PETSC_NULL,PETSC_NULL);
106:   PetscPrintf(PETSC_COMM_WORLD," Number of requested singular values: %d\n",nsv);
107:   SVDGetTolerances(svd,&tol,&maxit);
108:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
111:                     Display solution and clean up
112:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

114:   /* 
115:      Get number of converged singular triplets
116:   */
117:   SVDGetConverged(svd,&nconv);
118:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate singular triplets: %d\n\n",nconv);

120:   if (nconv>0) {
121:     /*
122:        Display singular values and relative errors
123:     */
124:     PetscPrintf(PETSC_COMM_WORLD,
125:          "          sigma           residual norm\n"
126:          "  --------------------- ------------------\n" );
127:     for( i=0; i<nconv; i++ ) {
128:       /* 
129:          Get converged singular triplets: i-th singular value is stored in sigma
130:       */
131:       SVDGetSingularTriplet(svd,i,&sigma,PETSC_NULL,PETSC_NULL);

133:       /*
134:          Compute the error associated to each singular triplet 
135:       */
136:       SVDComputeRelativeError(svd,i,&error);

138:       PetscPrintf(PETSC_COMM_WORLD,"       % 6f      ",sigma);
139:       PetscPrintf(PETSC_COMM_WORLD," % 12g\n",error);
140:     }
141:     PetscPrintf(PETSC_COMM_WORLD,"\n" );
142:   }
143: 
144:   /* 
145:      Free work space
146:   */
147:   SVDDestroy(svd);
148:   MatDestroy(A);
149:   SlepcFinalize();
150:   return 0;
151: }