Actual source code: ex15f.F90
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: #include <slepc/finclude/slepcsvd.h>
21: program ex15f
22: use slepcsvd
23: implicit none
25: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
26: ! Declarations
27: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
29: Mat :: A ! operator matrix
30: SVD :: svd ! singular value solver context
31: SVDType :: tname
32: PetscReal :: tol, error, sigma, mu
33: PetscInt :: n, i, j, Istart, Iend
34: PetscInt :: nsv, maxit, its, nconv
35: PetscMPIInt :: rank
36: PetscErrorCode :: ierr
37: PetscBool :: flg
38: PetscScalar :: alpha
39: PetscScalar, parameter :: one = 1.0
41: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
42: ! Beginning of program
43: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
46: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
47: n = 100
48: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
49: mu = PETSC_SQRT_MACHINE_EPSILON
50: PetscCallA(PetscOptionsGetReal(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-mu', mu, flg, ierr))
52: if (rank == 0) then
53: write (*, '(/a,i3,a,e12.4,a)') 'Lauchli SVD, n =', n, ', mu=', mu, ' (Fortran)'
54: end if
56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: ! Build the Lauchli matrix
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
61: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n + 1, n, ierr))
62: PetscCallA(MatSetFromOptions(A, ierr))
64: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
65: do i = Istart, Iend - 1
66: if (i == 0) then
67: do j = 0, n - 1
68: PetscCallA(MatSetValue(A, i, j, one, INSERT_VALUES, ierr))
69: end do
70: else
71: alpha = mu
72: PetscCallA(MatSetValue(A, i, i - 1, alpha, INSERT_VALUES, ierr))
73: end if
74: end do
76: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
77: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
79: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
80: ! Create the singular value solver and display info
81: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: ! ** Create singular value solver context
84: PetscCallA(SVDCreate(PETSC_COMM_WORLD, svd, ierr))
86: ! ** Set operators and problem type
87: PetscCallA(SVDSetOperators(svd, A, PETSC_NULL_MAT, ierr))
88: PetscCallA(SVDSetProblemType(svd, SVD_STANDARD, ierr))
90: ! ** Use thick-restart Lanczos as default solver
91: PetscCallA(SVDSetType(svd, SVDTRLANCZOS, ierr))
93: ! ** Set solver parameters at runtime
94: PetscCallA(SVDSetFromOptions(svd, ierr))
96: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
97: ! Solve the singular value system
98: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: PetscCallA(SVDSolve(svd, ierr))
101: PetscCallA(SVDGetIterationNumber(svd, its, ierr))
102: if (rank == 0) then
103: write (*, '(/a,i4)') ' Number of iterations of the method:', its
104: end if
106: ! ** Optional: Get some information from the solver and display it
107: PetscCallA(SVDGetType(svd, tname, ierr))
108: if (rank == 0) then
109: write (*, '(a,a)') ' Solution method: ', tname
110: end if
111: PetscCallA(SVDGetDimensions(svd, nsv, PETSC_NULL_INTEGER, PETSC_NULL_INTEGER, ierr))
112: if (rank == 0) then
113: write (*, '(a,i2)') ' Number of requested singular values:', nsv
114: end if
115: PetscCallA(SVDGetTolerances(svd, tol, maxit, ierr))
116: if (rank == 0) then
117: write (*, '(a,1pe11.4,a,i4)') ' Stopping condition: tol=', tol, ', maxit=', maxit
118: end if
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121: ! Display solution and clean up
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: ! ** Get number of converged singular triplets
125: PetscCallA(SVDGetConverged(svd, nconv, ierr))
126: if (rank == 0) then
127: write (*, '(a,i2/)') ' Number of converged approximate singular triplets:', nconv
128: end if
130: ! ** Display singular values and relative errors
131: if (nconv > 0) then
132: if (rank == 0) then
133: write (*, *) ' sigma relative error'
134: write (*, *) ' ----------------- ------------------'
135: end if
136: do i = 0, nconv - 1
137: ! ** Get i-th singular value
138: PetscCallA(SVDGetSingularTriplet(svd, i, sigma, PETSC_NULL_VEC, PETSC_NULL_VEC, ierr))
140: ! ** Compute the relative error for each singular triplet
141: PetscCallA(SVDComputeError(svd, i, SVD_ERROR_RELATIVE, error, ierr))
142: if (rank == 0) then
143: write (*, '(1p,a,e12.4,a,e12.4)') ' ', sigma, ' ', error
144: end if
146: end do
147: if (rank == 0) then
148: write (*, *)
149: end if
150: end if
152: ! ** Free work space
153: PetscCallA(SVDDestroy(svd, ierr))
154: PetscCallA(MatDestroy(A, ierr))
156: PetscCallA(SlepcFinalize(ierr))
157: end program ex15f
159: !/*TEST
160: !
161: ! test:
162: ! suffix: 1
163: ! filter: sed -e "s/[0-9]\.[0-9]*E[+-]\([0-9]*\)/removed/g"
164: !
165: !TEST*/