Actual source code: ex20f90.F90

      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, nfloat
     19:   PetscInt n
     20:   PetscErrorCode ierr
     21:   PetscMPIInt rank
     22:   PetscBool flg
     23:   PetscScalar, parameter :: one = 1.0, two = 2.0, three = 3.0
     24:   PetscScalar dots(3), dot
     25:   PetscReal, parameter :: tol = 1.e-10_PETSC_REAL_KIND

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

     31:   PetscCallA(PetscInitialize(ierr))

     33:   n = 20
     34:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))
     35:   nfloat = real(n, PETSC_REAL_KIND)
     36:   PetscCallMPIA(MPI_Comm_rank(PETSC_COMM_WORLD, rank, ierr))

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

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

     51:   PetscCallA(VecDuplicate(x, y, ierr))
     52:   PetscCallA(VecDuplicate(x, w, ierr))

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

     59:   PetscCallA(VecDuplicateVecs(x, 3_PETSC_INT_KIND, z, ierr))

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

     63:   PetscCallA(VecSet(x, one, ierr))
     64:   PetscCallA(VecSet(y, two, ierr))
     65:   PetscCallA(VecSet(z(1), one, ierr))
     66:   PetscCallA(VecSet(z(2), two, ierr))
     67:   PetscCallA(VecSet(z(3), three, ierr))

     69: !  Demonstrate various basic vector routines.

     71:   PetscCallA(VecDot(x, x, dot, ierr))
     72:   PetscCallA(VecMDot(x, 3_PETSC_INT_KIND, z, dots, ierr))

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

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

     92:   PetscCallA(VecScale(x, two, ierr))
     93:   PetscCallA(VecNorm(x, NORM_2, norm, ierr))
     94:   v = abs(norm - 2.0*sqrt(nfloat))
     95:   if (v > -tol .and. v < tol) v = 0.0
     96:   if (rank == 0) write (6, 130) v
     97: 130 format('VecScale ', 1pe9.2)

     99:   PetscCallA(VecCopy(x, w, ierr))
    100:   PetscCallA(VecNorm(w, NORM_2, norm, ierr))
    101:   v = abs(norm - 2.0*sqrt(nfloat))
    102:   if (v > -tol .and. v < tol) v = 0.0
    103:   if (rank == 0) write (6, 140) v
    104: 140 format('VecCopy ', 1pe9.2)

    106:   PetscCallA(VecAXPY(y, three, x, ierr))
    107:   PetscCallA(VecNorm(y, NORM_2, norm, ierr))
    108:   v = abs(norm - 8.0*sqrt(nfloat))
    109:   if (v > -tol .and. v < tol) v = 0.0
    110:   if (rank == 0) write (6, 150) v
    111: 150 format('VecAXPY ', 1pe9.2)

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

    120:   PetscCallA(VecSwap(x, y, ierr))
    121:   PetscCallA(VecNorm(y, NORM_2, norm, ierr))
    122:   v = abs(norm - 2.0*sqrt(nfloat))
    123:   if (v > -tol .and. v < tol) v = 0.0
    124:   if (rank == 0) write (6, 170) v
    125: 170 format('VecSwap ', 1pe9.2)

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

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

    140:   PetscCallA(VecPointwiseMult(w, y, x, ierr))
    141:   PetscCallA(VecNorm(w, NORM_2, norm, ierr))
    142:   v = abs(norm - 36.0*sqrt(nfloat))
    143:   if (v > -tol .and. v < tol) v = 0.0
    144:   if (rank == 0) write (6, 200) v
    145: 200 format('VecPointwiseMult ', 1pe9.2)

    147:   PetscCallA(VecPointwiseDivide(w, x, y, ierr))
    148:   PetscCallA(VecNorm(w, NORM_2, norm, ierr))
    149:   v = abs(norm - 9.0*sqrt(nfloat))
    150:   if (v > -tol .and. v < tol) v = 0.0
    151:   if (rank == 0) write (6, 210) v
    152: 210 format('VecPointwiseDivide ', 1pe9.2)

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

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

    174:   PetscCallA(VecDestroy(x, ierr))
    175:   PetscCallA(VecDestroy(y, ierr))
    176:   PetscCallA(VecDestroy(w, ierr))
    177:   PetscCallA(VecDestroyVecs(3_PETSC_INT_KIND, z, ierr))
    178:   PetscCallA(PetscFinalize(ierr))

    180: end

    182: !/*TEST
    183: !
    184: !     test:
    185: !
    186: !TEST*/