Actual source code: test3f.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> ./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: program main
20: #include <slepc/finclude/slepcmfn.h>
21: use slepcmfn
22: implicit none
24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
25: ! Declarations
26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
27: !
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 .eq. 0) then
53: write(*,100) n
54: endif
55: 100 format (/'nSquare root of Laplacian, n=',I3,' (Fortran)')
57: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58: ! Compute the discrete 2-D Laplacian
59: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
62: PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,Nt,Nt,ierr))
63: PetscCallA(MatSetFromOptions(A,ierr))
65: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
66: do II=Istart,Iend-1
67: i = II/n
68: j = II-i*n
69: val = -1.0
70: if (i .gt. 0) then
71: col = II-n
72: PetscCallA(MatSetValue(A,II,col,val,INSERT_VALUES,ierr))
73: end if
74: if (i .lt. n-1) then
75: col = II+n
76: PetscCallA(MatSetValue(A,II,col,val,INSERT_VALUES,ierr))
77: end if
78: if (j .gt. 0) then
79: col = II-1
80: PetscCallA(MatSetValue(A,II,col,val,INSERT_VALUES,ierr))
81: end if
82: if (j .lt. n-1) then
83: col = II+1
84: PetscCallA(MatSetValue(A,II,col,val,INSERT_VALUES,ierr))
85: end if
86: val = 4.0
87: PetscCallA(MatSetValue(A,II,II,val,INSERT_VALUES,ierr))
88: enddo
90: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
91: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
93: PetscCallA(MatCreateVecs(A,PETSC_NULL_VEC,v,ierr))
94: i = 0
95: val = 1.0
96: PetscCallA(VecSetValue(v,i,val,INSERT_VALUES,ierr))
97: PetscCallA(VecAssemblyBegin(v,ierr))
98: PetscCallA(VecAssemblyEnd(v,ierr))
99: PetscCallA(VecDuplicate(v,y,ierr))
101: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: ! Compute singular values
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: PetscCallA(MFNCreate(PETSC_COMM_WORLD,mfn,ierr))
106: PetscCallA(MFNSetOperator(mfn,A,ierr))
107: PetscCallA(MFNGetFN(mfn,f,ierr))
108: PetscCallA(FNSetType(f,FNSQRT,ierr))
110: ! ** test some interface functions
111: PetscCallA(MFNGetOperator(mfn,B,ierr))
112: PetscCallA(MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr))
113: PetscCallA(MFNSetOptionsPrefix(mfn,'myprefix_',ierr))
114: tol = 1e-4
115: maxit = 500
116: PetscCallA(MFNSetTolerances(mfn,tol,maxit,ierr))
117: ncv = 6
118: PetscCallA(MFNSetDimensions(mfn,ncv,ierr))
119: PetscCallA(MFNSetErrorIfNotConverged(mfn,PETSC_TRUE,ierr))
120: PetscCallA(MFNSetFromOptions(mfn,ierr))
122: ! ** query properties and print them
123: PetscCallA(MFNGetTolerances(mfn,tol,maxit,ierr))
124: if (rank .eq. 0) then
125: write(*,110) tol,maxit
126: endif
127: 110 format (/' Tolerance: ',F7.4,', maxit: ',I4)
128: PetscCallA(MFNGetDimensions(mfn,ncv,ierr))
129: if (rank .eq. 0) then
130: write(*,120) ncv
131: endif
132: 120 format (' Subspace dimension: ',I3)
133: PetscCallA(MFNGetErrorIfNotConverged(mfn,flg,ierr))
134: if (rank .eq. 0 .and. flg) then
135: write(*,*) 'Erroring out if convergence fails'
136: endif
138: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: ! Call the solver
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142: PetscCallA(MFNSolve(mfn,v,y,ierr))
143: PetscCallA(MFNGetConvergedReason(mfn,reason,ierr))
144: if (rank .eq. 0) then
145: write(*,130) reason
146: endif
147: 130 format (' Converged reason:',I2)
148: PetscCallA(MFNGetIterationNumber(mfn,its,ierr))
149: ! if (rank .eq. 0) then
150: ! write(*,140) its
151: ! endif
152: !140 format (' Number of iterations of the method:',I4)
154: PetscCallA(VecNorm(y,NORM_2,norm,ierr))
155: if (rank .eq. 0) then
156: write(*,150) norm
157: endif
158: 150 format (' sqrt(A)*v has norm ',F7.4)
160: PetscCallA(MFNDestroy(mfn,ierr))
161: PetscCallA(MatDestroy(A,ierr))
162: PetscCallA(VecDestroy(v,ierr))
163: PetscCallA(VecDestroy(y,ierr))
165: PetscCallA(SlepcFinalize(ierr))
166: end
168: !/*TEST
169: !
170: ! test:
171: ! suffix: 1
172: ! args: -log_exclude mfn
173: !
174: !TEST*/