Actual source code: ex36f.F

petsc-3.12.2 2019-11-22
Report Typos and Errors
  1: !
  2: !
  3: !   This program demonstrates use of PETSc dense matrices.
  4: !
  5:       program main
  6:  #include <petsc/finclude/petscsys.h>
  7:       use petscsys
  8:       implicit none
  9: 
 10:       PetscErrorCode ierr

 12:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 13:       if (ierr .ne. 0) then
 14:         print*,'Unable to initialize PETSc'
 15:         stop
 16:       endif

 18: !  Demo of PETSc-allocated dense matrix storage
 19:       call Demo1()

 21: !  Demo of user-allocated dense matrix storage
 22:       call Demo2()

 24:       call PetscFinalize(ierr)
 25:       end

 27: ! -----------------------------------------------------------------
 28: !
 29: !  Demo1 -  This subroutine demonstrates the use of PETSc-allocated dense
 30: !  matrix storage.  Here MatDenseGetArray() is used for direct access to the
 31: !  array that stores the dense matrix.  The user declares an array (aa(1))
 32: !  and index variable (ia), which are then used together to manipulate
 33: !  the array contents.
 34: !
 35: !  Note the use of PETSC_NULL_SCALAR in MatCreateSeqDense() to indicate that no
 36: !  storage is being provided by the user. (Do NOT pass a zero in that
 37: !  location.)
 38: !
 39:       subroutine Demo1()
 40:  #include <petsc/finclude/petscmat.h>
 41:       use petscmat
 42:       implicit none

 44:       Mat         A
 45:       PetscInt   n,m
 46:       PetscErrorCode ierr
 47:       PetscScalar aa(1)
 48:       PetscOffset ia

 50:       n = 4
 51:       m = 5

 53: !  Create matrix

 55: !      call MatCreateSeqDense(PETSC_COMM_SELF,m,n,PETSC_NULL_SCALAR,     &
 56: !     &     A,ierr)

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

 61:       call MatCreate(PETSC_COMM_SELF,A,ierr)
 62:       call MatSetSizes(A,m,n,m,n,ierr)
 63:       call MatSetType(A,MATSEQDENSE,ierr)
 64:       call MatSetUp(A,ierr)

 66: !  Access array storage
 67:       call MatDenseGetArray(A,aa,ia,ierr)

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

 72:       call MatDenseRestoreArray(A,aa,ia,ierr)

 74: !  Finalize matrix assembly
 75:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
 76:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

 78: !  View matrix
 79:       call MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)

 81: !  Clean up
 82:       call MatDestroy(A,ierr)
 83:       return
 84:       end

 86: ! -----------------------------------------------------------------
 87: !
 88: !  Demo2 -  This subroutine demonstrates the use of user-allocated dense
 89: !  matrix storage.
 90: !
 91:       subroutine Demo2()
 92:  #include <petsc/finclude/petscmat.h>
 93:       use petscmat
 94:       implicit none

 96:       PetscInt   n,m
 97:       PetscErrorCode ierr
 98:       parameter (m=5,n=4)
 99:       Mat       A
100:       PetscScalar    aa(m,n)

102: !  Create matrix
103:       call MatCreateSeqDense(PETSC_COMM_SELF,m,n,aa,A,ierr)

105: !  Set matrix values directly
106:       call FillUpMatrix(m,n,aa)

108: !  Finalize matrix assembly
109:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)
110:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)

112: !  View matrix
113:       call MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr)

115: !  Clean up
116:       call MatDestroy(A,ierr)
117:       return
118:       end

120: ! -----------------------------------------------------------------

122:       subroutine FillUpMatrix(m,n,X)
123:       PetscInt          m,n,i,j
124:       PetscScalar      X(m,n)

126:       do 10, j=1,n
127:         do 20, i=1,m
128:           X(i,j) = 1.0/real(i+j-1)
129:  20     continue
130:  10   continue
131:       return
132:       end