Actual source code: ex1f90.F90

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

      6: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      7: !                   Variable declarations
      8: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: !
     10: !  Variables:
     11: !     x, y, w - vectors
     12: !     z       - array of vectors
     13: !
     14:   Vec x, y, w
     15:   Vec, pointer :: z(:)
     16:   PetscInt, pointer :: ranges(:)
     17:   PetscReal norm, v, v1, v2
     18:   PetscInt n, ithree
     19:   PetscErrorCode ierr
     20:   PetscMPIInt rank
     21:   PetscBool flg
     22:   PetscScalar one, two, three
     23:   PetscScalar dots(3), dot
     24:   PetscReal nfloat

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

     30:   PetscCallA(PetscInitialize(ierr))
     31:   one = 1.0
     32:   two = 2.0
     33:   three = 3.0
     34:   n = 20
     35:   ithree = 3

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

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

     47:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
     48:   PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr))
     49:   PetscCallA(VecSetFromOptions(x, ierr))

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

     54:   PetscCallA(VecDuplicate(x, y, ierr))
     55:   PetscCallA(VecDuplicate(x, w, ierr))

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

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

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

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

     72: !  Demonstrate various basic vector routines.

     74:   PetscCallA(VecDot(x, x, dot, ierr))
     75:   PetscCallA(VecMDot(x, ithree, z, dots, ierr))

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

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

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

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

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

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

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

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

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

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

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

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

    174:   PetscCallA(VecGetOwnershipRanges(x, ranges, ierr))
    175: !      PetscCallA(VecRestoreOwnershipRanges(x, ranges, ierr))

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

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

    187: !
    188: !/*TEST
    189: !
    190: !     test:
    191: !       nsize: 2
    192: !
    193: !TEST*/