1: #include <petsc/finclude/petscvec.h>
  2: program main
  3:   use petscvec
  4:   implicit none

  6: !
  7: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  8: !                   Variable declarations
  9: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 10: !
 11: !  Variables:
 12: !     x, y, w - vectors
 13: !     z       - array of vectors
 14: !
 15:   type(tVec) x, y, w
 16:   type(tVec), pointer :: z(:)

 18:   PetscReal norm, v, v1, v2, tol
 19:   PetscInt n, ithree
 20:   PetscErrorCode ierr
 21:   PetscMPIInt rank
 22:   PetscBool flg
 23:   PetscScalar one, two, three
 24:   PetscScalar dots(3), dot
 25:   PetscReal nfloat

 27: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 28: !                 Beginning of program
 29: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 31:   PetscCallA(PetscInitialize(ierr))
 32:   tol = 1.e-10_PETSC_REAL_KIND
 33:   one = 1.0
 34:   two = 2.0
 35:   three = 3.0
 36:   n = 20
 37:   ithree = 3

 39:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
 40:   nfloat = n
 41:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))

 43: !  Create a vector, specifying only its global dimension.
 44: !  When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
 45: !  the vector format (currently parallel
 46: !  or sequential) is determined at runtime.  Also, the parallel
 47: !  partitioning of the vector is determined by PETSc at runtime.
 48: !
 49:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
 50:   PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr))
 51:   PetscCallA(VecSetFromOptions(x, ierr))

 53: !  Duplicate some work vectors (of the same format and
 54: !  partitioning as the initial vector).

 56:   PetscCallA(VecDuplicate(x, y, ierr))
 57:   PetscCallA(VecDuplicate(x, w, ierr))

 59: !  Duplicate more work vectors (of the same format and
 60: !  partitioning as the initial vector).  Here we duplicate
 61: !  an array of vectors, which is often more convenient than
 62: !  duplicating individual ones.

 64:   PetscCallA(VecDuplicateVecs(x, ithree, z, ierr))

 66: !  Set the vectors to entries to a constant value.

 68:   PetscCallA(VecSet(x, one, ierr))
 69:   PetscCallA(VecSet(y, two, ierr))
 70:   PetscCallA(VecSet(z(1), one, ierr))
 71:   PetscCallA(VecSet(z(2), two, ierr))
 72:   PetscCallA(VecSet(z(3), three, ierr))

 74: !  Demonstrate various basic vector routines.

 76:   PetscCallA(VecDot(x, x, dot, ierr))
 77:   PetscCallA(VecMDot(x, ithree, z, dots, ierr))

 79: !  Note: If using a complex numbers version of PETSc, then
 80: !  PETSC_USE_COMPLEX is defined in the makefiles; otherwise,
 81: !  (when using real numbers) it is undefined.

 83:   if (rank == 0) then
 84: #if defined(PETSC_USE_COMPLEX)
 85:     write (6, 100) int(PetscRealPart(dot))
 86:     write (6, 110) int(PetscRealPart(dots(1))), int(PetscRealPart(dots(2))), int(PetscRealPart(dots(3)))
 87: #else
 88:     write (6, 100) int(dot)
 89:     write (6, 110) int(dots(1)), int(dots(2)), int(dots(3))
 90: #endif
 91:     write (6, 120)
 92:   end if
 93: 100 format('Vector length ', i6)
 94: 110 format('Vector length ', 3(i6))
 95: 120 format('All other values should be near zero')

 97:   PetscCallA(VecScale(x, two, ierr))
 98:   PetscCallA(VecNorm(x, NORM_2, norm, ierr))
 99:   v = abs(norm - 2.0*sqrt(nfloat))
100:   if (v > -tol .and. v < tol) v = 0.0
101:   if (rank == 0) write (6, 130) v
102: 130 format('VecScale ', 1pe9.2)

104:   PetscCallA(VecCopy(x, w, ierr))
105:   PetscCallA(VecNorm(w, NORM_2, norm, ierr))
106:   v = abs(norm - 2.0*sqrt(nfloat))
107:   if (v > -tol .and. v < tol) v = 0.0
108:   if (rank == 0) write (6, 140) v
109: 140 format('VecCopy ', 1pe9.2)

111:   PetscCallA(VecAXPY(y, three, x, ierr))
112:   PetscCallA(VecNorm(y, NORM_2, norm, ierr))
113:   v = abs(norm - 8.0*sqrt(nfloat))
114:   if (v > -tol .and. v < tol) v = 0.0
115:   if (rank == 0) write (6, 150) v
116: 150 format('VecAXPY ', 1pe9.2)

118:   PetscCallA(VecAYPX(y, two, x, ierr))
119:   PetscCallA(VecNorm(y, NORM_2, norm, ierr))
120:   v = abs(norm - 18.0*sqrt(nfloat))
121:   if (v > -tol .and. v < tol) v = 0.0
122:   if (rank == 0) write (6, 160) v
123: 160 format('VecAYXP ', 1pe9.2)

125:   PetscCallA(VecSwap(x, y, ierr))
126:   PetscCallA(VecNorm(y, NORM_2, norm, ierr))
127:   v = abs(norm - 2.0*sqrt(nfloat))
128:   if (v > -tol .and. v < tol) v = 0.0
129:   if (rank == 0) write (6, 170) v
130: 170 format('VecSwap ', 1pe9.2)

132:   PetscCallA(VecNorm(x, NORM_2, norm, ierr))
133:   v = abs(norm - 18.0*sqrt(nfloat))
134:   if (v > -tol .and. v < tol) v = 0.0
135:   if (rank == 0) write (6, 180) v
136: 180 format('VecSwap ', 1pe9.2)

138:   PetscCallA(VecWAXPY(w, two, x, y, ierr))
139:   PetscCallA(VecNorm(w, NORM_2, norm, ierr))
140:   v = abs(norm - 38.0*sqrt(nfloat))
141:   if (v > -tol .and. v < tol) v = 0.0
142:   if (rank == 0) write (6, 190) v
143: 190 format('VecWAXPY ', 1pe9.2)

145:   PetscCallA(VecPointwiseMult(w, y, x, ierr))
146:   PetscCallA(VecNorm(w, NORM_2, norm, ierr))
147:   v = abs(norm - 36.0*sqrt(nfloat))
148:   if (v > -tol .and. v < tol) v = 0.0
149:   if (rank == 0) write (6, 200) v
150: 200 format('VecPointwiseMult ', 1pe9.2)

152:   PetscCallA(VecPointwiseDivide(w, x, y, ierr))
153:   PetscCallA(VecNorm(w, NORM_2, norm, ierr))
154:   v = abs(norm - 9.0*sqrt(nfloat))
155:   if (v > -tol .and. v < tol) v = 0.0
156:   if (rank == 0) write (6, 210) v
157: 210 format('VecPointwiseDivide ', 1pe9.2)

159:   dots(1) = one
160:   dots(2) = three
161:   dots(3) = two
162:   PetscCallA(VecSet(x, one, ierr))
163:   PetscCallA(VecMAXPY(x, ithree, dots, z, ierr))
164:   PetscCallA(VecNorm(z(1), NORM_2, norm, ierr))
165:   v = abs(norm - sqrt(nfloat))
166:   if (v > -tol .and. v < tol) v = 0.0
167:   PetscCallA(VecNorm(z(2), NORM_2, norm, ierr))
168:   v1 = abs(norm - 2.0*sqrt(nfloat))
169:   if (v1 > -tol .and. v1 < tol) v1 = 0.0
170:   PetscCallA(VecNorm(z(3), NORM_2, norm, ierr))
171:   v2 = abs(norm - 3.0*sqrt(nfloat))
172:   if (v2 > -tol .and. v2 < tol) v2 = 0.0
173:   if (rank == 0) write (6, 220) v, v1, v2
174: 220 format('VecMAXPY ', 3(1pe9.2))

176: !  Free work space.  All PETSc objects should be destroyed when they
177: !  are no longer needed.

179:   PetscCallA(VecDestroy(x, ierr))
180:   PetscCallA(VecDestroy(y, ierr))
181:   PetscCallA(VecDestroy(w, ierr))
182:   PetscCallA(VecDestroyVecs(ithree, z, ierr))
183:   PetscCallA(PetscFinalize(ierr))

185: end

187: !/*TEST
188: !
189: !     test:
190: !
191: !TEST*/