Actual source code: ex1f.F
petsc-3.8.4 2018-03-24
1: !
2: !
3: !/*T
4: ! Concepts: vectors^basic routines
5: ! Processors: n
6: !T*/
7: !
8: ! -----------------------------------------------------------------------
10: program main
11: #include <petsc/finclude/petscvec.h>
12: use petscvec
13: implicit none
15: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
16: ! Variable declarations
17: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
18: !
19: ! Variables:
20: ! x, y, w - vectors
21: ! z - array of vectors
23: Vec x,y,w,z(5)
24: PetscReal norm,v,v1,v2
25: PetscInt n,ithree
26: PetscBool flg
27: PetscErrorCode ierr
28: PetscMPIInt rank
29: PetscScalar one,two,three
30: PetscScalar dots(3),dot
31: character*(40) name
32: PetscReal nfloat
34: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
35: ! Beginning of program
36: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
39: if (ierr .ne. 0) then
40: print*,'Unable to initialize PETSc'
41: stop
42: endif
43: one = 1.0
44: two = 2.0
45: three = 3.0
46: n = 20
47: ithree = 3
48: call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, &
49: & '-n',n,flg,ierr)
50: nfloat = n
51: call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
53: ! Create a vector, specifying only its global dimension.
54: ! When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
55: ! the vector format (currently parallel
56: ! or sequential) is determined at runtime. Also, the parallel
57: ! partitioning of the vector is determined by PETSc at runtime.
58: !
59: ! Routines for creating particular vector types directly are:
60: ! VecCreateSeq() - uniprocessor vector
61: ! VecCreateMPI() - distributed vector, where the user can
62: ! determine the parallel partitioning
63: ! VecCreateShared() - parallel vector that uses shared memory
64: ! (available only on the SGI); otherwise,
65: ! is the same as VecCreateMPI()
66: !
67: ! VecCreate(), VecSetSizes() and VecSetFromOptions() allows one
68: ! to determine at runtime which version to use
69: ! with the options -vec_type mpi or -vec_type shared
70: !
71: call VecCreate(PETSC_COMM_WORLD,x,ierr)
72: call VecSetSizes(x,PETSC_DECIDE,n,ierr)
73: call VecSetFromOptions(x,ierr)
74: call VecGetType(x,name,ierr)
76: ! Duplicate some work vectors (of the same format and
77: ! partitioning as the initial vector).
79: call VecDuplicate(x,y,ierr)
80: call VecDuplicate(x,w,ierr)
82: ! Duplicate more work vectors (of the same format and
83: ! partitioning as the initial vector). Here we duplicate
84: ! an array of vectors, which is often more convenient than
85: ! duplicating individual ones.
87: call VecDuplicateVecs(x,ithree,z,ierr)
89: ! Set the vectors to entries to a constant value.
91: call VecSet(x,one,ierr)
93: call VecSet(y,two,ierr)
94: call VecSet(z(1),one,ierr)
95: call VecSet(z(2),two,ierr)
96: call VecSet(z(3),three,ierr)
98: ! Demonstrate various basic vector routines.
100: call VecDot(x,x,dot,ierr)
101: call VecMDot(x,ithree,z,dots,ierr)
103: ! Note: If using a complex numbers version of PETSc, then
104: ! PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
105: ! (when using real numbers) it is undefined.
107: if (rank .eq. 0) then
108: #if defined(PETSC_USE_COMPLEX)
109: write(6,100) int(PetscRealPart(dot))
110: write(6,110) int(PetscRealPart(dots(1))), &
111: & int(PetscRealPart(dots(2))), &
112: & int(PetscRealPart(dots(3)))
113: #else
114: write(6,100) int(dot)
115: write(6,110) int(dots(1)),int(dots(2)),int(dots(3))
116: #endif
117: write(6,120)
118: endif
119: 100 format ('Vector length ',i6)
120: 110 format ('Vector length ',3(i6))
121: 120 format ('All other values should be near zero')
123: call VecScale(x,two,ierr)
124: call VecNorm(x,NORM_2,norm,ierr)
125: v = abs(norm-2.0*sqrt(nfloat))
126: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
127: if (rank .eq. 0) write(6,130) v
128: 130 format ('VecScale ',1pe8.2)
130: call VecCopy(x,w,ierr)
131: call VecNorm(w,NORM_2,norm,ierr)
132: v = abs(norm-2.0*sqrt(nfloat))
133: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
134: if (rank .eq. 0) write(6,140) v
135: 140 format ('VecCopy ',1pe8.2)
137: call VecAXPY(y,three,x,ierr)
138: call VecNorm(y,NORM_2,norm,ierr)
139: v = abs(norm-8.0*sqrt(nfloat))
140: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
141: if (rank .eq. 0) write(6,150) v
142: 150 format ('VecAXPY ',1pe8.2)
144: call VecAYPX(y,two,x,ierr)
145: call VecNorm(y,NORM_2,norm,ierr)
146: v = abs(norm-18.0*sqrt(nfloat))
147: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
148: if (rank .eq. 0) write(6,160) v
149: 160 format ('VecAYXP ',1pe8.2)
151: call VecSwap(x,y,ierr)
152: call VecNorm(y,NORM_2,norm,ierr)
153: v = abs(norm-2.0*sqrt(nfloat))
154: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
155: if (rank .eq. 0) write(6,170) v
156: 170 format ('VecSwap ',1pe8.2)
158: call VecNorm(x,NORM_2,norm,ierr)
159: v = abs(norm-18.0*sqrt(nfloat))
160: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
161: if (rank .eq. 0) write(6,180) v
162: 180 format ('VecSwap ',1pe8.2)
164: call VecWAXPY(w,two,x,y,ierr)
165: call VecNorm(w,NORM_2,norm,ierr)
166: v = abs(norm-38.0*sqrt(nfloat))
167: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
168: if (rank .eq. 0) write(6,190) v
169: 190 format ('VecWAXPY ',1pe8.2)
171: call VecPointwiseMult(w,y,x,ierr)
172: call VecNorm(w,NORM_2,norm,ierr)
173: v = abs(norm-36.0*sqrt(nfloat))
174: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
175: if (rank .eq. 0) write(6,200) v
176: 200 format ('VecPointwiseMult ',1pe8.2)
178: call VecPointwiseDivide(w,x,y,ierr)
179: call VecNorm(w,NORM_2,norm,ierr)
180: v = abs(norm-9.0*sqrt(nfloat))
181: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
182: if (rank .eq. 0) write(6,210) v
183: 210 format ('VecPointwiseDivide ',1pe8.2)
186: dots(1) = one
187: dots(2) = three
188: dots(3) = two
189: call VecSet(x,one,ierr)
190: call VecMAXPY(x,ithree,dots,z,ierr)
191: call VecNorm(z(1),NORM_2,norm,ierr)
192: v = abs(norm-sqrt(nfloat))
193: if (v .gt. -PETSC_SMALL .and. v .lt. PETSC_SMALL) v = 0.0
194: call VecNorm(z(2),NORM_2,norm,ierr)
195: v1 = abs(norm-2.0*sqrt(nfloat))
196: if (v1 .gt. -PETSC_SMALL .and. v1 .lt. PETSC_SMALL) v1 = 0.0
197: call VecNorm(z(3),NORM_2,norm,ierr)
198: v2 = abs(norm-3.0*sqrt(nfloat))
199: if (v2 .gt. -PETSC_SMALL .and. v2 .lt. PETSC_SMALL) v2 = 0.0
200: if (rank .eq. 0) write(6,220) v,v1,v2
201: 220 format ('VecMAXPY ',3(2x,1pe8.2))
203: ! Free work space. All PETSc objects should be destroyed when they
204: ! are no longer needed.
206: call VecDestroy(x,ierr)
207: call VecDestroy(y,ierr)
208: call VecDestroy(w,ierr)
209: call VecDestroyVecs(ithree,z,ierr)
210: call PetscFinalize(ierr)
212: end