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

 10:   PetscErrorCode ierr

 12:   PetscCallA(PetscInitialize(ierr))

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

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

 20:   PetscCallA(PetscFinalize(ierr))
 21: end

 23: ! -----------------------------------------------------------------
 24: !
 25: !  Demo1 -  This subroutine demonstrates the use of PETSc-allocated dense
 26: !  matrix storage.  Here MatDenseGetArray() is used for direct access to the
 27: !  array that stores the dense matrix.
 28: !
 29: !  Note the use of PETSC_NULL_SCALAR_ARRAY in MatCreateSeqDense() to indicate that no
 30: !  storage is being provided by the user. (Do NOT pass a zero in that
 31: !  location.)
 32: !
 33: subroutine Demo1()
 34: #include <petsc/finclude/petscmat.h>
 35:   use petscmat
 36:   implicit none

 38:   Mat A
 39:   PetscInt n, m
 40:   PetscErrorCode ierr
 41:   PetscScalar, pointer :: aa(:, :)

 43:   n = 4
 44:   m = 5

 46: !  Create matrix

 48:   PetscCall(MatCreate(PETSC_COMM_SELF, A, ierr))
 49:   PetscCall(MatSetSizes(A, m, n, m, n, ierr))
 50:   PetscCall(MatSetType(A, MATSEQDENSE, ierr))
 51:   PetscCall(MatSetUp(A, ierr))

 53: !  Access array storage
 54:   PetscCall(MatDenseGetArray(A, aa, ierr))

 56: !  Set matrix values directly
 57:   PetscCall(FillUpMatrix(m, n, aa))

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

 61: !  Finalize matrix assembly
 62:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 63:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

 65: !  View matrix
 66:   PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))

 68: !  Clean up
 69:   PetscCall(MatDestroy(A, ierr))
 70: end

 72: ! -----------------------------------------------------------------
 73: !
 74: !  Demo2 -  This subroutine demonstrates the use of user-allocated dense
 75: !  matrix storage.
 76: !
 77: subroutine Demo2()
 78: #include <petsc/finclude/petscmat.h>
 79:   use petscmat
 80:   implicit none

 82:   PetscInt n, m
 83:   PetscErrorCode ierr
 84:   parameter(m=5, n=4)
 85:   Mat A
 86:   PetscScalar aa(m, n)

 88: !  Create matrix
 89:   PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, n, aa, A, ierr))

 91: !  Set matrix values directly
 92:   PetscCall(FillUpMatrix(m, n, aa))

 94: !  Finalize matrix assembly
 95:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 96:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

 98: !  View matrix
 99:   PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr))

101: !  Clean up
102:   PetscCall(MatDestroy(A, ierr))
103: end

105: ! -----------------------------------------------------------------

107: subroutine FillUpMatrix(m, n, X)
108:   PetscInt m, n, i, j
109:   PetscScalar X(m, n)

111:   do 10, j = 1, n
112:     do 20, i = 1, m
113:       X(i, j) = 1.0/real(i + j - 1)
114: 20    continue
115: 10    continue
116:     end