Actual source code: test3f.F
slepc-3.11.2 2019-07-30
1: !
2: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: ! SLEPc - Scalable Library for Eigenvalue Problem Computations
4: ! Copyright (c) 2002-2019, 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: call SlepcInitialize(PETSC_NULL_CHARACTER,ierr)
47: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
48: n = 4
49: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
50: & '-n',n,flg,ierr)
51: Nt = n*n
53: if (rank .eq. 0) then
54: write(*,100) n
55: endif
56: 100 format (/'nSquare root of Laplacian, n=',I3,' (Fortran)')
58: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: ! Compute the discrete 2-D Laplacian
60: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: call MatCreate(PETSC_COMM_WORLD,A,ierr)
63: call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,Nt,Nt,ierr)
64: call MatSetFromOptions(A,ierr)
65: call MatSetUp(A,ierr)
67: call MatGetOwnershipRange(A,Istart,Iend,ierr)
68: do II=Istart,Iend-1
69: i = II/n
70: j = II-i*n
71: val = -1.0
72: if (i .gt. 0) then
73: col = II-n
74: call MatSetValue(A,II,col,val,INSERT_VALUES,ierr)
75: end if
76: if (i .lt. n-1) then
77: col = II+n
78: call MatSetValue(A,II,col,val,INSERT_VALUES,ierr)
79: end if
80: if (j .gt. 0) then
81: col = II-1
82: call MatSetValue(A,II,col,val,INSERT_VALUES,ierr)
83: end if
84: if (j .lt. n-1) then
85: col = II+1
86: call MatSetValue(A,II,col,val,INSERT_VALUES,ierr)
87: end if
88: val = 4.0
89: call MatSetValue(A,II,II,val,INSERT_VALUES,ierr)
90: enddo
92: call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
93: call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)
95: call MatCreateVecs(A,PETSC_NULL_VEC,v,ierr)
96: i = 0
97: val = 1.0
98: call VecSetValue(v,i,val,INSERT_VALUES,ierr)
99: call VecAssemblyBegin(v,ierr)
100: call VecAssemblyEnd(v,ierr)
101: call VecDuplicate(v,y,ierr)
103: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: ! Compute singular values
105: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: call MFNCreate(PETSC_COMM_WORLD,mfn,ierr)
108: call MFNSetOperator(mfn,A,ierr)
109: call MFNGetFN(mfn,f,ierr)
110: call FNSetType(f,FNSQRT,ierr)
112: ! ** test some interface functions
113: call MFNGetOperator(mfn,B,ierr)
114: call MatView(B,PETSC_VIEWER_STDOUT_WORLD,ierr)
115: call MFNSetOptionsPrefix(mfn,'myprefix_',ierr)
116: tol = 1e-4
117: maxit = 500
118: call MFNSetTolerances(mfn,tol,maxit,ierr)
119: ncv = 6
120: call MFNSetDimensions(mfn,ncv,ierr)
121: call MFNSetErrorIfNotConverged(mfn,PETSC_TRUE,ierr)
122: call MFNSetFromOptions(mfn,ierr)
124: ! ** query properties and print them
125: call MFNGetTolerances(mfn,tol,maxit,ierr)
126: if (rank .eq. 0) then
127: write(*,110) tol,maxit
128: endif
129: 110 format (/' Tolerance: ',F7.4,', maxit: ',I4)
130: call MFNGetDimensions(mfn,ncv,ierr)
131: if (rank .eq. 0) then
132: write(*,120) ncv
133: endif
134: 120 format (' Subspace dimension: ',I3)
135: call MFNGetErrorIfNotConverged(mfn,flg,ierr)
136: if (rank .eq. 0 .and. flg) then
137: write(*,*) 'Erroring out if convergence fails'
138: endif
140: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: ! Call the solver
142: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144: call MFNSolve(mfn,v,y,ierr)
145: call MFNGetConvergedReason(mfn,reason,ierr)
146: if (rank .eq. 0) then
147: write(*,130) reason
148: endif
149: 130 format (' Converged reason:',I2)
150: call MFNGetIterationNumber(mfn,its,ierr)
151: ! if (rank .eq. 0) then
152: ! write(*,140) its
153: ! endif
154: !140 format (' Number of iterations of the method:',I4)
156: call VecNorm(y,NORM_2,norm,ierr)
157: if (rank .eq. 0) then
158: write(*,150) norm
159: endif
160: 150 format (' sqrt(A)*v has norm ',F7.4)
162: call MFNDestroy(mfn,ierr)
163: call MatDestroy(A,ierr)
164: call VecDestroy(v,ierr)
165: call VecDestroy(y,ierr)
167: call SlepcFinalize(ierr)
168: end
170: !/*TEST
171: !
172: ! test:
173: ! suffix: 1
174: ! args: -info_exclude mfn -log_exclude mfn
175: !
176: !TEST*/