Actual source code: ex15f.F

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: !
  2: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3: !  SLEPc - Scalable Library for Eigenvalue Problem Computations
  4: !  Copyright (c) 2002-2021, Universitat Politecnica de Valencia, Spain
  5: !
  6: !  This file is part of SLEPc.
  7: !  SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: !
 10: !  Program usage: mpiexec -n <np> ./ex15f [-help] [-n <n>] [-mu <mu>] [all SLEPc options]
 11: !
 12: !  Description: Singular value decomposition of the Lauchli matrix.
 13: !
 14: !  The command line options are:
 15: !    -n <n>, where <n> = matrix dimension.
 16: !    -mu <mu>, where <mu> = subdiagonal value.
 17: !
 18: ! ----------------------------------------------------------------------
 19: !
 20:       program main
 21: #include <slepc/finclude/slepcsvd.h>
 22:       use slepcsvd
 23:       implicit none

 25: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 26: !     Declarations
 27: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 28: !
 29: !  Variables:
 30: !     A     operator matrix
 31: !     svd   singular value solver context

 33:       Mat            A
 34:       SVD            svd
 35:       SVDType        tname
 36:       PetscReal      tol, error, sigma, mu
 37:       PetscInt       n, i, j, Istart, Iend
 38:       PetscInt       nsv, maxit, its, nconv
 39:       PetscMPIInt    rank
 40:       PetscErrorCode ierr
 41:       PetscBool      flg
 42:       PetscScalar    one, alpha

 44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 45: !     Beginning of program
 46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 48:       call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
 49:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 50:       n = 100
 51:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,  &
 52:      &                        '-n',n,flg,ierr)
 53:       mu = PETSC_SQRT_MACHINE_EPSILON
 54:       call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
 55:      &                        '-mu',mu,flg,ierr)

 57:       if (rank .eq. 0) then
 58:         write(*,100) n, mu
 59:       endif
 60:  100  format (/'Lauchli SVD, n =',I3,', mu=',E12.4,' (Fortran)')

 62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 63: !     Build the Lauchli matrix
 64: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 66:       call MatCreate(PETSC_COMM_WORLD,A,ierr)
 67:       call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n+1,n,ierr)
 68:       call MatSetFromOptions(A,ierr)
 69:       call MatSetUp(A,ierr)

 71:       call MatGetOwnershipRange(A,Istart,Iend,ierr)
 72:       one = 1.0
 73:       do i=Istart,Iend-1
 74:         if (i .eq. 0) then
 75:           do j=0,n-1
 76:             call MatSetValue(A,i,j,one,INSERT_VALUES,ierr)
 77:           end do
 78:         else
 79:           alpha = mu
 80:           call MatSetValue(A,i,i-1,alpha,INSERT_VALUES,ierr)
 81:         end if
 82:       enddo

 84:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 85:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

 87: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 88: !     Create the singular value solver and display info
 89: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 91: !     ** Create singular value solver context
 92:       call SVDCreate(PETSC_COMM_WORLD,svd,ierr)

 94: !     ** Set operators and problem type
 95:       call SVDSetOperators(svd,A,PETSC_NULL_MAT,ierr)
 96:       call SVDSetProblemType(svd,SVD_STANDARD,ierr)

 98: !     ** Use thick-restart Lanczos as default solver
 99:       call SVDSetType(svd,SVDTRLANCZOS,ierr)

101: !     ** Set solver parameters at runtime
102:       call SVDSetFromOptions(svd,ierr)

104: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: !     Solve the singular value system
106: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

108:       call SVDSolve(svd,ierr)
109:       call SVDGetIterationNumber(svd,its,ierr)
110:       if (rank .eq. 0) then
111:         write(*,110) its
112:       endif
113:  110  format (/' Number of iterations of the method:',I4)

115: !     ** Optional: Get some information from the solver and display it
116:       call SVDGetType(svd,tname,ierr)
117:       if (rank .eq. 0) then
118:         write(*,120) tname
119:       endif
120:  120  format (' Solution method: ',A)
121:       call SVDGetDimensions(svd,nsv,PETSC_NULL_INTEGER,                 &
122:      &                      PETSC_NULL_INTEGER,ierr)
123:       if (rank .eq. 0) then
124:         write(*,130) nsv
125:       endif
126:  130  format (' Number of requested singular values:',I2)
127:       call SVDGetTolerances(svd,tol,maxit,ierr)
128:       if (rank .eq. 0) then
129:         write(*,140) tol, maxit
130:       endif
131:  140  format (' Stopping condition: tol=',1P,E11.4,', maxit=',I4)

133: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134: !     Display solution and clean up
135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

137: !     ** Get number of converged singular triplets
138:       call SVDGetConverged(svd,nconv,ierr)
139:       if (rank .eq. 0) then
140:         write(*,150) nconv
141:       endif
142:  150  format (' Number of converged approximate singular triplets:',I2/)

144: !     ** Display singular values and relative errors
145:       if (nconv.gt.0) then
146:         if (rank .eq. 0) then
147:           write(*,*) '       sigma          relative error'
148:           write(*,*) ' ----------------- ------------------'
149:         endif
150:         do i=0,nconv-1
151: !         ** Get i-th singular value
152:           call SVDGetSingularTriplet(svd,i,sigma,PETSC_NULL_VEC,        &
153:      &         PETSC_NULL_VEC,ierr)

155: !         ** Compute the relative error for each singular triplet
156:           call SVDComputeError(svd,i,SVD_ERROR_RELATIVE,error,ierr)
157:           if (rank .eq. 0) then
158:             write(*,160) sigma, error
159:           endif
160:  160      format (1P,'   ',E12.4,'       ',E12.4)

162:         enddo
163:         if (rank .eq. 0) then
164:           write(*,*)
165:         endif
166:       endif

168: !     ** Free work space
169:       call SVDDestroy(svd,ierr)
170:       call MatDestroy(A,ierr)

172:       call SlepcFinalize(ierr)
173:       end

175: !/*TEST
176: !
177: !   test:
178: !      suffix: 1
179: !      filter: sed -e "s/[0-9]\.[0-9]*E[+-]\([0-9]*\)/removed/g"
180: !
181: !TEST*/