Actual source code: ex36f.F90

  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.  The user declares an array (aa(1))
 28: !  and index variable (ia), which are then used together to manipulate
 29: !  the array contents.
 30: !
 31: !  Note the use of PETSC_NULL_SCALAR in MatCreateSeqDense() to indicate that no
 32: !  storage is being provided by the user. (Do NOT pass a zero in that
 33: !  location.)
 34: !
 35:       subroutine Demo1()
 36: #include <petsc/finclude/petscmat.h>
 37:       use petscmat
 38:       implicit none

 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:       PetscCall(MatCreate(PETSC_COMM_SELF,A,ierr))
 52:       PetscCall(MatSetSizes(A,m,n,m,n,ierr))
 53:       PetscCall(MatSetType(A,MATSEQDENSE,ierr))
 54:       PetscCall(MatSetUp(A,ierr))

 56: !  Access array storage
 57:       PetscCall(MatDenseGetArray(A,aa,ia,ierr))

 59: !  Set matrix values directly
 60:       PetscCall(FillUpMatrix(m,n,aa(ia+1))

 62:       PetscCall(MatDenseRestoreArray(A,aa,ia,ierr))

 64: !  Finalize matrix assembly
 65:       PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
 66:       PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))

 68: !  View matrix
 69:       PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr))

 71: !  Clean up
 72:       PetscCall(MatDestroy(A,ierr))
 73:       return
 74:       end

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

 86:       PetscInt   n,m
 87:       PetscErrorCode ierr
 88:       parameter (m=5,n=4)
 89:       Mat       A
 90:       PetscScalar    aa(m,n)

 92: !  Create matrix
 93:       PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,m,n,aa,A,ierr))

 95: !  Set matrix values directly
 96:       PetscCall(FillUpMatrix(m,n,aa)

 98: !  Finalize matrix assembly
 99:       PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
100:       PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))

102: !  View matrix
103:       PetscCall(MatView(A,PETSC_VIEWER_STDOUT_SELF,ierr))

105: !  Clean up
106:       PetscCall(MatDestroy(A,ierr))
107:       return
108:       end

110: ! -----------------------------------------------------------------

112:       subroutine FillUpMatrix(m,n,X)
113:       PetscInt          m,n,i,j
114:       PetscScalar      X(m,n)

116:       do 10, j=1,n
117:         do 20, i=1,m
118:           X(i,j) = 1.0/real(i+j-1)
119:  20     continue
120:  10   continue
121:       return
122:       end