Actual source code: ex15f.F90
slepc-3.21.1 2024-04-26
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-, 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: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER,ierr))
49: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr))
50: n = 100
51: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))
52: mu = PETSC_SQRT_MACHINE_EPSILON
53: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-mu',mu,flg,ierr))
55: if (rank .eq. 0) then
56: write(*,100) n, mu
57: endif
58: 100 format (/'Lauchli SVD, n =',I3,', mu=',E12.4,' (Fortran)')
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: ! Build the Lauchli matrix
62: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
64: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
65: PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,n+1,n,ierr))
66: PetscCallA(MatSetFromOptions(A,ierr))
68: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
69: one = 1.0
70: do i=Istart,Iend-1
71: if (i .eq. 0) then
72: do j=0,n-1
73: PetscCallA(MatSetValue(A,i,j,one,INSERT_VALUES,ierr))
74: end do
75: else
76: alpha = mu
77: PetscCallA(MatSetValue(A,i,i-1,alpha,INSERT_VALUES,ierr))
78: end if
79: enddo
81: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
82: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
84: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: ! Create the singular value solver and display info
86: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
88: ! ** Create singular value solver context
89: PetscCallA(SVDCreate(PETSC_COMM_WORLD,svd,ierr))
91: ! ** Set operators and problem type
92: PetscCallA(SVDSetOperators(svd,A,PETSC_NULL_MAT,ierr))
93: PetscCallA(SVDSetProblemType(svd,SVD_STANDARD,ierr))
95: ! ** Use thick-restart Lanczos as default solver
96: PetscCallA(SVDSetType(svd,SVDTRLANCZOS,ierr))
98: ! ** Set solver parameters at runtime
99: PetscCallA(SVDSetFromOptions(svd,ierr))
101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: ! Solve the singular value system
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: PetscCallA(SVDSolve(svd,ierr))
106: PetscCallA(SVDGetIterationNumber(svd,its,ierr))
107: if (rank .eq. 0) then
108: write(*,110) its
109: endif
110: 110 format (/' Number of iterations of the method:',I4)
112: ! ** Optional: Get some information from the solver and display it
113: PetscCallA(SVDGetType(svd,tname,ierr))
114: if (rank .eq. 0) then
115: write(*,120) tname
116: endif
117: 120 format (' Solution method: ',A)
118: PetscCallA(SVDGetDimensions(svd,nsv,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr))
119: if (rank .eq. 0) then
120: write(*,130) nsv
121: endif
122: 130 format (' Number of requested singular values:',I2)
123: PetscCallA(SVDGetTolerances(svd,tol,maxit,ierr))
124: if (rank .eq. 0) then
125: write(*,140) tol, maxit
126: endif
127: 140 format (' Stopping condition: tol=',1P,E11.4,', maxit=',I4)
129: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: ! Display solution and clean up
131: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
133: ! ** Get number of converged singular triplets
134: PetscCallA(SVDGetConverged(svd,nconv,ierr))
135: if (rank .eq. 0) then
136: write(*,150) nconv
137: endif
138: 150 format (' Number of converged approximate singular triplets:',I2/)
140: ! ** Display singular values and relative errors
141: if (nconv.gt.0) then
142: if (rank .eq. 0) then
143: write(*,*) ' sigma relative error'
144: write(*,*) ' ----------------- ------------------'
145: endif
146: do i=0,nconv-1
147: ! ** Get i-th singular value
148: PetscCallA(SVDGetSingularTriplet(svd,i,sigma,PETSC_NULL_VEC,PETSC_NULL_VEC,ierr))
150: ! ** Compute the relative error for each singular triplet
151: PetscCallA(SVDComputeError(svd,i,SVD_ERROR_RELATIVE,error,ierr))
152: if (rank .eq. 0) then
153: write(*,160) sigma, error
154: endif
155: 160 format (1P,' ',E12.4,' ',E12.4)
157: enddo
158: if (rank .eq. 0) then
159: write(*,*)
160: endif
161: endif
163: ! ** Free work space
164: PetscCallA(SVDDestroy(svd,ierr))
165: PetscCallA(MatDestroy(A,ierr))
167: PetscCallA(SlepcFinalize(ierr))
168: end
170: !/*TEST
171: !
172: ! test:
173: ! suffix: 1
174: ! filter: sed -e "s/[0-9]\.[0-9]*E[+-]\([0-9]*\)/removed/g"
175: !
176: !TEST*/