1! 2! 3! This program demonstrates use of PETSc dense matrices. 4! 5program main 6#include <petsc/finclude/petscsys.h> 7 use petscsys 8 implicit none 9 10 PetscErrorCode ierr 11 12 PetscCallA(PetscInitialize(ierr)) 13 14! Demo of PETSc-allocated dense matrix storage 15 call Demo1() 16 17! Demo of user-allocated dense matrix storage 18 call Demo2() 19 20 PetscCallA(PetscFinalize(ierr)) 21end 22 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! 33subroutine Demo1() 34#include <petsc/finclude/petscmat.h> 35 use petscmat 36 implicit none 37 38 Mat A 39 PetscInt n, m 40 PetscErrorCode ierr 41 PetscScalar, pointer :: aa(:, :) 42 43 n = 4 44 m = 5 45 46! Create matrix 47 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)) 52 53! Access array storage 54 PetscCall(MatDenseGetArray(A, aa, ierr)) 55 56! Set matrix values directly 57 PetscCall(FillUpMatrix(m, n, aa)) 58 59 PetscCall(MatDenseRestoreArray(A, aa, ierr)) 60 61! Finalize matrix assembly 62 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) 63 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) 64 65! View matrix 66 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr)) 67 68! Clean up 69 PetscCall(MatDestroy(A, ierr)) 70end 71 72! ----------------------------------------------------------------- 73! 74! Demo2 - This subroutine demonstrates the use of user-allocated dense 75! matrix storage. 76! 77subroutine Demo2() 78#include <petsc/finclude/petscmat.h> 79 use petscmat 80 implicit none 81 82 PetscInt n, m 83 PetscErrorCode ierr 84 parameter(m=5, n=4) 85 Mat A 86 PetscScalar aa(m, n) 87 88! Create matrix 89 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, n, aa, A, ierr)) 90 91! Set matrix values directly 92 PetscCall(FillUpMatrix(m, n, aa)) 93 94! Finalize matrix assembly 95 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) 96 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) 97 98! View matrix 99 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_SELF, ierr)) 100 101! Clean up 102 PetscCall(MatDestroy(A, ierr)) 103end 104 105! ----------------------------------------------------------------- 106 107subroutine FillUpMatrix(m, n, X) 108 PetscInt m, n, i, j 109 PetscScalar X(m, n) 110 111 do 10, j = 1, n 112 do 20, i = 1, m 113 X(i, j) = 1.0/real(i + j - 1) 11420 continue 11510 continue 116 end 117