1: !
  2: !
  3: !   This program demonstrates use of PETSc dense matrices.
  4: !
  5: #include <petsc/finclude/petscmat.h>
  6: module ex36fmodule
  7:   use petscmat
  8:   implicit none
  9: ! -----------------------------------------------------------------
 10: !
 11: !  Demo1 -  This subroutine demonstrates the use of PETSc-allocated dense
 12: !  matrix storage.  Here MatDenseGetArray() is used for direct access to the
 13: !  array that stores the dense matrix.
 14: !
 15: !  Note the use of PETSC_NULL_SCALAR_ARRAY in MatCreateSeqDense() to indicate that no
 16: !  storage is being provided by the user. (Do NOT pass a zero in that
 17: !  location.)
 18: !
 19: contains
 20:   subroutine Demo1()

 22:     Mat A
 23:     PetscInt n, m
 24:     PetscErrorCode ierr
 25:     PetscScalar, pointer :: aa(:, :)

 27:     n = 4
 28:     m = 5

 30: !  Create matrix

 32:     PetscCall(MatCreate(PETSC_COMM_SELF, A, ierr))
 33:     PetscCall(MatSetSizes(A, m, n, m, n, ierr))
 34:     PetscCall(MatSetType(A, MATSEQDENSE, ierr))
 35:     PetscCall(MatSetUp(A, ierr))

 37: !  Access array storage
 38:     PetscCall(MatDenseGetArray(A, aa, ierr))

 40: !  Set matrix values directly
 41:     PetscCall(FillUpMatrix(m, n, aa))

 43:     PetscCall(MatDenseRestoreArray(A, aa, ierr))

 45: !  Finalize matrix assembly
 46:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 47:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

 49: !  View matrix
 50:     PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))

 52: !  Clean up
 53:     PetscCall(MatDestroy(A, ierr))
 54:   end

 56: ! -----------------------------------------------------------------
 57: !
 58: !  Demo2 -  This subroutine demonstrates the use of user-allocated dense
 59: !  matrix storage.
 60: !
 61:   subroutine Demo2()

 63:     PetscInt n, m
 64:     PetscErrorCode ierr
 65:     parameter(m=5, n=4)
 66:     Mat A
 67:     PetscScalar aa(m, n)

 69: !  Create matrix
 70:     PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, n, aa, A, ierr))

 72: !  Set matrix values directly
 73:     PetscCall(FillUpMatrix(m, n, aa))

 75: !  Finalize matrix assembly
 76:     PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 77:     PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

 79: !  View matrix
 80:     PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))

 82: !  Clean up
 83:     PetscCall(MatDestroy(A, ierr))
 84:   end

 86: ! -----------------------------------------------------------------

 88:   subroutine FillUpMatrix(m, n, X)
 89:     PetscInt m, n, i, j
 90:     PetscScalar X(m, n)

 92:     do j = 1, n
 93:       do i = 1, m
 94:         X(i, j) = 1.0/real(i + j - 1)
 95:       end do
 96:     end do
 97:     end module ex36fmodule

 99:     end program main
100:     use ex36fmodule
101:     implicit none

103:     PetscErrorCode ierr

105:     PetscCallA(PetscInitialize(ierr))

107: !  Demo of PETSc-allocated dense matrix storage
108:     call Demo1()

110: !  Demo of user-allocated dense matrix storage
111:     call Demo2()

113:     PetscCallA(PetscFinalize(ierr))
114:   end