Actual source code: test3f.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> ./test3f [-help] [-n <n>] [all SLEPc options]
11: !
12: ! Description: square root of the 2-D Laplacian.
13: !
14: ! The command line options are:
15: ! -n <n>, where <n> = matrix rows and columns
16: !
17: ! ----------------------------------------------------------------------
18: !
19: #include <slepc/finclude/slepcmfn.h>
20: program test3f
21: use slepcmfn
22: implicit none
24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
25: ! Declarations
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
28: Mat :: A, B
29: MFN :: mfn
30: FN :: f
31: MFNConvergedReason :: reason
32: Vec :: v, y
33: PetscInt :: Nt, n, i, j, II
34: PetscInt :: Istart, maxit, ncv
35: PetscInt :: col, its, Iend
36: PetscScalar :: val
37: PetscReal :: tol, norm
38: PetscMPIInt :: rank
39: PetscErrorCode :: ierr
40: PetscBool :: flg
42: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
43: ! Beginning of program
44: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
46: PetscCallA(SlepcInitialize(PETSC_NULL_CHARACTER, ierr))
47: PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))
48: n = 4
49: PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
50: Nt = n*n
52: if (rank == 0) then
53: write (*, '(/a,i3,a)') 'Square root of Laplacian, n=', n, ' (Fortran)'
54: end if
56: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
57: ! Compute the discrete 2-D Laplacian
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
61: PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, Nt, Nt, ierr))
62: PetscCallA(MatSetFromOptions(A, ierr))
64: PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))
65: do II = Istart, Iend - 1
66: i = II/n
67: j = II - i*n
68: val = -1.0
69: if (i > 0) then
70: col = II - n
71: PetscCallA(MatSetValue(A, II, col, val, INSERT_VALUES, ierr))
72: end if
73: if (i < n - 1) then
74: col = II + n
75: PetscCallA(MatSetValue(A, II, col, val, INSERT_VALUES, ierr))
76: end if
77: if (j > 0) then
78: col = II - 1
79: PetscCallA(MatSetValue(A, II, col, val, INSERT_VALUES, ierr))
80: end if
81: if (j < n - 1) then
82: col = II + 1
83: PetscCallA(MatSetValue(A, II, col, val, INSERT_VALUES, ierr))
84: end if
85: val = 4.0
86: PetscCallA(MatSetValue(A, II, II, val, INSERT_VALUES, ierr))
87: end do
89: PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
90: PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))
92: PetscCallA(MatCreateVecs(A, PETSC_NULL_VEC, v, ierr))
93: i = 0
94: val = 1.0
95: PetscCallA(VecSetValue(v, i, val, INSERT_VALUES, ierr))
96: PetscCallA(VecAssemblyBegin(v, ierr))
97: PetscCallA(VecAssemblyEnd(v, ierr))
98: PetscCallA(VecDuplicate(v, y, ierr))
100: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: ! Compute singular values
102: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: PetscCallA(MFNCreate(PETSC_COMM_WORLD, mfn, ierr))
105: PetscCallA(MFNSetOperator(mfn, A, ierr))
106: PetscCallA(MFNGetFN(mfn, f, ierr))
107: PetscCallA(FNSetType(f, FNSQRT, ierr))
109: ! ** test some interface functions
110: PetscCallA(MFNGetOperator(mfn, B, ierr))
111: PetscCallA(MatView(B, PETSC_VIEWER_STDOUT_WORLD, ierr))
112: PetscCallA(MFNSetOptionsPrefix(mfn, 'myprefix_', ierr))
113: tol = 1e-4
114: maxit = 500
115: PetscCallA(MFNSetTolerances(mfn, tol, maxit, ierr))
116: ncv = 6
117: PetscCallA(MFNSetDimensions(mfn, ncv, ierr))
118: PetscCallA(MFNSetErrorIfNotConverged(mfn, PETSC_TRUE, ierr))
119: PetscCallA(MFNSetFromOptions(mfn, ierr))
121: ! ** query properties and print them
122: PetscCallA(MFNGetTolerances(mfn, tol, maxit, ierr))
123: if (rank == 0) then
124: write (*, '(/a,f7.4,a,i4)') ' Tolerance: ', tol, ', maxit: ', maxit
125: end if
126: PetscCallA(MFNGetDimensions(mfn, ncv, ierr))
127: if (rank == 0) then
128: write (*, '(a,i3)') ' Subspace dimension: ', ncv
129: end if
130: PetscCallA(MFNGetErrorIfNotConverged(mfn, flg, ierr))
131: if (rank == 0 .and. flg) then
132: write (*, *) 'Erroring out if convergence fails'
133: end if
135: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: ! Call the solver
137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: PetscCallA(MFNSolve(mfn, v, y, ierr))
140: PetscCallA(MFNGetConvergedReason(mfn, reason, ierr))
141: if (rank == 0) then
142: write (*, '(a,i2)') ' Converged reason:', reason
143: end if
144: PetscCallA(MFNGetIterationNumber(mfn, its, ierr))
145: ! if (rank==0) then
146: ! write(*, '(a,i4)') ' Number of iterations of the method:', its
147: ! end if
149: PetscCallA(VecNorm(y, NORM_2, norm, ierr))
150: if (rank == 0) then
151: write (*, '(a,f7.4)') ' sqrt(A)*v has norm ', norm
152: end if
154: PetscCallA(MFNDestroy(mfn, ierr))
155: PetscCallA(MatDestroy(A, ierr))
156: PetscCallA(VecDestroy(v, ierr))
157: PetscCallA(VecDestroy(y, ierr))
159: PetscCallA(SlepcFinalize(ierr))
160: end program test3f
162: !/*TEST
163: !
164: ! test:
165: ! suffix: 1
166: ! args: -log_exclude mfn
167: !
168: !TEST*/