Actual source code: ex36f.F

  1: !
  2: !
  3: !   This program demonstrates use of PETSc dense matrices.
  4: !
  5:       program main

 7:  #include finclude/petscsys.h

  9:       PetscErrorCode ierr

 11:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)

 13: !  Demo of PETSc-allocated dense matrix storage
 14:       call Demo1()

 16: !  Demo of user-allocated dense matrix storage
 17:       call Demo2()

 19:       call PetscFinalize(ierr)
 20:       end

 22: ! -----------------------------------------------------------------
 23: !
 24: !  Demo1 -  This subroutine demonstrates the use of PETSc-allocated dense
 25: !  matrix storage.  Here MatGetArray() is used for direct access to the
 26: !  array that stores the dense matrix.  The user declares an array (aa(1))
 27: !  and index variable (ia), which are then used together to manipulate
 28: !  the array contents.
 29: !
 30: !  Note the use of PETSC_NULL_SCALAR in MatCreateSeqDense() to indicate that no
 31: !  storage is being provided by the user. (Do NOT pass a zero in that
 32: !  location.)
 33: !
 34:       subroutine Demo1()

 36:  #include finclude/petscsys.h
 37:  #include finclude/petscmat.h
 38:  #include finclude/petscviewer.h

 40:       Mat         A
 41:       PetscInt   n,m
 42:       PetscErrorCode ierr
 43:       PetscScalar aa(1)
 44:       PetscOffset ia

 46:       n = 4
 47:       m = 5

 49: !  Create matrix

 51: !      call MatCreateSeqDense(PETSC_COMM_SELF,m,n,PETSC_NULL_SCALAR,     &
 52: !     &     A,ierr)

 54: !  Using MatCreate() instead of  MatCreateSeqDense() as above to avoid Nag F90 errors
 55: !  However both cases are equivalent

 57:       call MatCreate(PETSC_COMM_SELF,A,ierr)
 58:       call MatSetSizes(A,m,n,m,n,ierr)
 59:       call MatSetType(A,MATSEQDENSE,ierr)

 61: !  Access array storage
 62:       call MatGetArray(A,aa,ia,ierr)

 64: !  Set matrix values directly
 65:       call FillUpMatrix(m,n,aa(ia+1))

 67:       call MatRestoreArray(A,aa,ia,ierr)

 69: !  Finalize matrix assembly
 70:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 71:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

 73: !  View matrix
 74:       call MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)

 76: !  Clean up
 77:       call MatDestroy(A,ierr)
 78:       return
 79:       end

 81: ! -----------------------------------------------------------------
 82: !
 83: !  Demo2 -  This subroutine demonstrates the use of user-allocated dense
 84: !  matrix storage.
 85: !
 86:       subroutine Demo2()

 88:  #include finclude/petscsys.h
 89:  #include finclude/petscmat.h
 90:  #include finclude/petscviewer.h

 92:       PetscInt   n,m
 93:       PetscErrorCode ierr
 94:       parameter (m=5,n=4)
 95:       Mat       A
 96:       PetscScalar    aa(m,n)

 98: !  Create matrix
 99:       call MatCreateSeqDense(PETSC_COMM_SELF,m,n,aa,A,ierr)

101: !  Set matrix values directly
102:       call FillUpMatrix(m,n,aa)

104: !  Finalize matrix assembly
105:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
106:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

108: !  View matrix
109:       call MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)

111: !  Clean up
112:       call MatDestroy(A,ierr)
113:       return
114:       end

116: ! -----------------------------------------------------------------

118:       subroutine FillUpMatrix(m,n,X)
119:       PetscInt          m,n,i,j
120:       PetscScalar      X(m,n)

122:       do 10, j=1,n
123:         do 20, i=1,m
124:           X(i,j) = 1.d0/(i+j-1)
125:  20     continue
126:  10   continue
127:       return
128:       end